(*  Ray Trace Through the IONOSPHERE  Mathematica Notebook  Post Processsing  *)
(*  Copyright (C)2018 Sasan Ardalan                                           *)
(*  This program is free software:you can redistribute it and/or modify       *)
(*  it under the terms of the GNU General Public License as published by      *)
(*  the Free Software Foundation,either version 3 of the License,or           *)
(*  at your option) any later version.This program is distributed             *)
(*  in the hope that it will be useful,but WITHOUT ANY WARRANTY               *)
(*  without even the implied warranty of                                      *)
(*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.See the               *)
(*  GNU General Public License for more details.You should have               *)
(*  received a copy of the GNU General Public License                         *)
(*  along with this program.If not,see<http://www.gnu.org/licenses/>.         *)
(*  Author:Sasan Ardalan                                                      *)
(*  Date:December 28,2018                                                     *)
(*  http://www.radiocalc.com                                                  *)
(*                                                                            *)
(* The Mathematica Code processes the Ray Matrix computed using               *)
(* The Mathematica Notebook for Ray Tracing which is  based on                *)
(* the FORTRAN Code by :                                                      *)
(* R. Michael Jones, Judith J. Stephenson,   A Versatile Three-Dimensional    *)
(* Ray Tracing Compute Program for Radio Waves in the Ionosphere,             *) 
(* OT Report 75-75, US Department of Commerce,                                *)
(* October 1975                                                                *)
(* Visit http://www.radiocalc.com for details                                  *)

(*  The program reads in the Ray Matrix rayMatrix.wdx                         *)
(* It also reads in and uses the parameters stored in W.wdx                   *)
(*  Plots the Propagation Vector as Arrows Showing Direction in 3D            *)

(*                                                                             *)
(***************************************************************************)
(***********************   RAY_PATH Ray Vector 3D Plot      *****************)
(***************************************************************************)

(* *********************************************************************)
(* *************** PolCarAnalysis   *************************************)
(***   CONVERTS SPHERICAL COORDINATES TO CARTESIAN                     ***)
(* Based on Jones and Stephenson "POL CAR" Routine with Modifications  ***)


(*
 * Global
* S
* XO, RO, XXX
* R
*)
(* ****************************************************************************)

PolCarAnalysis[  IENTR_]:=
 Module [{ SINA, COSA,SINP,COSP,TEMP, Done},
Done=False;

(* Print ["Ray in Analysis=",Ray]; *)
If [IENTR==0,
(*  CONVERTS SPHERICAL COORDINATES TO CARTESIAN   *)
If [!(Re[Ray[[5]]]==0.&& Re[Ray[[6]]]==0.),  
VERT=0.;
                 SINA=Sin[Re[Ray[[2]]]];
                COSA=Cos[ Re[Ray[[2]]]];
                SINP=Sin[Re[Ray[[3]]]];
                COSP=Cos[ Re[Ray[[3]]]];

 
                 XO[[1]]=Re[Ray[[1]]]*SINA*COSP;
                 XO[[2]]=Re[Ray[[1]]]*SINA*SINP;
                 XO[[3]]=Re[Ray[[1]]]*COSA;

 

    XXX[[4]]=Re[Ray[[4]]]*SINA*COSP+Re[Ray[[5]]]*COSA*COSP-Re[Ray[[6]]]*SINP;
    XXX[[5]]=Re[Ray[[4]]]*SINA*SINP+Re[Ray[[5]]]*COSA*SINP+Re[Ray[[6]]]*COSP;
                 XXX[[6]]=Re[Ray[[4]]]*COSA-Re[Ray[[5]]]*SINA;

 
     (*Print["PolCar 0 000000 R(5)=", R[[5]], " R(6)=",R[[6]]];*)
 Return[];
 
];
 
(*  VERTICAL INCIDENCE   *)
VERT=1.;
             RO[[1]]=Ray[[1]];
             RO[[2]]=Ray[[2]];
             RO[[3]]=Ray[[3]];
              RO[[4]]=Sign [Ray[[4]]];

Print["PolCar 0 VERTICAL   1111 R(5)=", R[[5]], " R(6)=",R[[6]]];
Return[];

];
If[ VERT ==0.,
    XXX[[1]]=XO [[1]]+S*XXX[[4]];

     XXX[[2]]=XO[[2]]+S*XXX[[5]];
              XXX[[3]]=XO[[3]]+S*XXX[[6]];
              TEMP=Sqrt[XXX[[1]]^2+XXX[[2]]^2];
               R[[1]]=Sqrt[XXX[[1]]^2+XXX[[2]]^2+XXX[[3]]^2];
      R[[2]]=ArcTan[Re[TEMP],Re[XXX[[3]]]];
                 R[[3]]=ArcTan[Re[XXX[[2]]],Re[XXX[[1]]]];
                 R[[4]]=(XXX[[1]]*XXX[[4]]+XXX[[2]]*XXX[[5]]+XXX[[3]]*XXX[[6]])/R[[1]];
                R[[5]]=(XXX[[3]]*(XXX[[1]]*XXX[[4]]+XXX[[2]]*XXX[[5]])-(XXX[[1]]^2 +XXX[[2]]^2)*XXX[[6]])/(R[[1]]*TEMP);
                  R[[6]]=(XXX[[1]]*XXX[[5]]-XXX[[2]]*XXX[[4]])/TEMP;
        For [ii=1, ii≤3,ii++,
              XO[[ii]]=XXX[[ii]];
      ];

       Print["PolCar 8888 R(1)=", R[[1]], "  S=",S];

                  Return[];
];

(*  VERTICAL INCIDENCE   *)
R[[1]]=RO[[1]]+RO[[4]]*S;
        R[[2]]=RO[[2]];
        R[[3]]=RO[[3]];
        R[[4]]=RO[[4]];
        R[[5]]=0.;
        R[[6]]=0.;
  Print["PolCar 99999 R(1)=", R[[1]]," RO(1) =",RO[[1]], " RO(4) = ",RO[[4]], "  S= ",S];

 For [ii=1, ii≤4,ii++,
              RO[[ii]]=R[[ii]];
      ];
Return[];


];

ClearAll;
XO =Table[ 0,20];
XXX=Table[ 0,20];


(* Single Ray *)
rayNumber=120;
Print["Ray Number=",rayNumber];
 RayMatrix=Import["rayMatrix.wdx"];
lengthMatrix=Length[RayMatrix]
Print ["Length Matrix=",lengthMatrix];
Ray= RayMatrix[[rayNumber] ];
 (* Print["RAY=",Ray];  *)



PolCarAnalysis[0];

(* Print["XXX=",XXX];
Print["XO=",XO];
Print["RO=",RO]; *)


kx= Re[XXX[[4]]] ;
ky=Re[XXX[[5]]];
kz=Re[XXX[[6]]];
a={0,0,0};
b={kx,ky,kz};
Print["Mag b=",Sqrt[kx^2+ky^2+kz^2]];
Print["b=",b]
Print["Ray Number=",rayNumber];
Graphics3D[{Red,Arrow[{a,b}]},Axes->True,Boxed->True,AxesLabel->{x,y,z}]

(* Test *)
a={0,0,0};
b={kx,ky,kz};
c={kx,ky,-kz};
Print["Mag[b]=",Sqrt[kx^2+ky^2+kz^2]];
Print["b=",b]
Print["Ray Number=",rayNumber];
Graphics3D[{{Red,Arrow[{a,b} ]},{Red,Arrow[{a,c} ]}},Axes->True,Boxed->True,AxesLabel->{x,y,z}]

(* Multiple rays   *)

lenghRayMatrix=Length[RayMatrix];
arrowArray={};

(* For[ii=2, ii< 600, ii++, *)
For[ii=20, ii< 300, ii++, 
Ray=Re[RayMatrix[[ii]] ];
(*Print["RAY=",Ray]; *)
XO =Table[ 0,20];
XXX=Table[ 0,20];
PolCarAnalysis[0];
kx= Re[XXX[[4]]] ;
ky=Re[XXX[[5]]];
kz=Re[XXX[[6]]];
a={0,0,0};
  b={kx,ky,kz};
AppendTo[arrowArray,{Red,Arrow[{a,b} ]}];
AppendTo[arrowArray,{Red,Text[i,b] }];

magb= Sqrt[b[[1]]^2+b[[2]]^2+b[[3]]^2];
Print ["ii=", ii, "  magb =",magb];

];

Graphics3D[arrowArray,Axes->True,Boxed->True,AxesLabel->{x,y,z}]

801