(*  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                                                  *)
(*        WORK IN PROGRESS , INCOMPLETE                                       *)

                                                                    
(* 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                   *)
(*  Translation and Modification of:  *) 
(*           PRINTR FORTRAN Routine                                           *)
(*     Purpose is to generate the ray path ( ray trace) Height versus Range   *)
(*            and other data                                                  *)
(* The Ray path is then plotted ( range versus height)                        *)
(*  This Mathematica Code is semi working but range is off as well as height  *)
(*  due to differences found n ArcTan versus Atan2 (FORTRAN)                  *)
(* and other descrepencies.                                                   *)
(* Ultamately this code needs to be debugged.                                 *)
(* We are providing it as a start. IT IS NOT WORKING PROPERLY                *)
(*                                                                            *)
(*                                                                            *)
(* **************************************************************************)
(* *************** RayPath   ************************************************)
(*   Adopted from PRINTR FORTRAN ROUTINE                                     *)
(*                                                                           *)
(*         Hard Coded for 45 Degree Elevation and 45 Degree Azimuth          *)
(* **************************************************************************)

RayPath[  ]:=
 Module [{R0, SPL, CPL,SP,CP,SL, CL,XR,YR,ZR,CTHR,STHR,PHIR,ALPH,G,G1,NR,NP,PDEV,ABSORB,DOPP,TLOND,TLATD,AZ,EL,NHOP,V,H,STH,CTH,
XP,YP,ZP,EPS,ETA,ZETA,RCE2,RCE, RANGE ,ELL,SR,ANGA,AZDEV,ANA,  SANA, SPHI,CPHI, AZA,
 F,THETA,PHI,AZ1,RCVRH,HOP,BETA,
PLAT,PLON,TLAT,TLON,EARTHR,XMTRH},

(* Print ["RayPath Ray=",Ray]; *)

THETA=Re[Ray[[2]]];
PHI=Re[Ray[[3]]];
F=WPost[[7]];
AZ1=WPost[[10]];

AZ1= 45/DEGS; (* for now*)
BETA=WPost[[14]];

BETA= 45/DEGS; (* for now *)

RCVRH=WPost[[20]];
HOP=WPost[[22]];

PLAT=WPost[[24]];
PLON=WPost[[25]];
TLAT=WPost[[4]];
 TLON=WPost[[5]];

 (* Print["Ray Path Degrees PLAT=",PLAT, " PLON=",PLON];
Print["Ray Path Degrees TLAT=",TLAT, " TLON=",TLON]; *)
PLAT=PLAT/DEGS;
PLON=PLON/DEGS;
TLAT=TLAT/DEGS;
TLON=TLON/DEGS;


 (* Print["Ray Path PLAT=",PLAT, " PLON=",PLON, " AZ1= ", AZ1, "BETA=", BETA]; *)
(* Print["Ray Path THETA=",THETA, " PHI=",PHI, " F= ", F, "AZ1=", AZ1];  *)

 EARTHR=WPost[[2]];
 XMTRH=WPost[[3]];

  (* Print["Ray Path ", "EARTHR= ", EARTHR," XMTRH=", XMTRH];   *)


NWHY="TYPE SPEC";

SPL=Sin[PLON-TLON];
CPL=Cos[ (PLON-TLON)];
SP=Sin[PLAT];
CP=Cos[PLAT];
SL=Sin[ TLAT];
CL=Cos[TLAT];
(**********MATRIX TO ROTATE COORDINATES  *)
G=Table[0,3,3];
G1=Table[0,3,3];
 


G[[1,1]]=CPL*SP*CL-CP*SL;
G[[1,2]]=SPL*SP;
G [[1,3]]=-SL*SP*CPL-CL*CP;
G[[2,1]]=-SPL*CL;
G[[2,2]]=CPL;
G[[2,3]]=SL*SPL;
G [[3,1]]=CL*CP*CPL+SP*SL;
G[[3,2]]=CP*SPL;
G[[3,3]]=-SL*CP*CPL+SP*CL;

  (* Print["G matrix=",G];  *) 


R0=Re[EARTHR+XMTRH];
(*********CARTESIAN COORDINATES OF TRANSMITTER *)
XR=R0*G[[1,1]];
YR=R0*G[[2,1]];
ZR=R0*G[[3,1]];
CTHR=G[[3,1]];
STHR=Sin [ArcCos [CTHR]];
PHIR=ArcTan[YR,XR];
ALPH=ArcTan [G[[3,2]],G[[3,3]]];
G1=Inverse[G];


 (* Print["G1 matrix=",G1];   *)

NR=6;
NP=0;
PDEV=0;
ABSORB=0;
DOPP=0.;
TLOND=TLON*DEGS;
If[TLOND <0., TLOND=TLOND+360.];
TLATD=TLAT*DEGS;
If[TLATD <0., TLATD=TLATD+360.];
AZ=AZ1*DEGS;
EL=BETA*DEGS;
NHOP=HOP;
V=0;

H=Re[Ray[[1]]]-EARTHR;
STH=Sin [THETA];
CTH=Cos [ THETA];

(* PHI=PHI-Pi;  SHA *)
(* PHI=-0.754334; *)
PHI=PHI-Pi; (* Mathematica FORTRAN ArcTan Issue *)

(* Print["R(1)=",Ray[[1]]]; *)


(*********CARTESIAN COORDINATES OF RAY POINT,ORIGIN AT TRANSMITTER  *)
XP=Re[Ray[[1]]]*STH*Cos [PHI]-XR;
YP=Re[Ray[[1]]]*STH*Sin[  PHI]-YR;
ZP=Re[Ray[[1]]]*CTH-ZR;
 (* Print["Raypath ", " XR= ", XR," YR= ",YR," ZR= ",ZR];  
  Print["Raypath ", " XP= ", XP," YP= ",YP," ZP= ",ZP];  *)

(********CARTESIAN COORDINATES OF RAY POINT,ORIGIN AT TRANSMITTER AND  ROTATED *)
EPS=XP*G1[[1,1]]+YP*G1[[1,2]]+ZP*G1[[1,3]];

ETA=XP*G1[[2,1]]+YP*G1[[2,2]]+ZP*G1[[2,3]];
ZETA=XP*G1[[3,1]]+YP*G1[[3,2]]+ZP*G1[[3,3]];
RCE2=ETA^2+ZETA^2;
RCE=Sqrt [RCE2];

(*********GROUND RANGE  Corrected based on Jones *)
RANGE=EARTHR*ArcTan [RCE,(EARTHR+EPS+XMTRH)];
 



  (*Print["RayPath Range RCE,EARTHR,EPS,XMTRH","R0= ",R0," RANGE=",RANGE," RCE- ",  RCE," EARTHR= ",EARTHR," EPS=  ",EPS," XMTRH= ",XMTRH," ArcTan [RCE,EARTHR+EPS+XMTRH]=",ArcTan [RCE,EARTHR+EPS+XMTRH]*DEGS];  *)
 
(**** ANGLE OF WAVE NORMAL WITH LOCAL HORIZONTAL *)
ELL=ArcTan[Re[Ray[[4]]],Sqrt [Re[Ray[[5]]]^2+Re[Ray[[6]]]^2]]*DEGS;

If[ELL > 90., ELL=90-(ELL-90)];
(*********STRAIGHT LINE DISTANCE IE FROM TRANSMITTER TO RAY POINT  *)
SR=Sqrt [RCE2+EPS^2];
ANGA=ArcTan [ETA,ZETA];

AZDEV=180.-Mod[540.-(AZ1-ANGA)*DEGS,360.];
(* IF (R(5).NE.0..OR.R(6).NE.0.) GO TO 34  *) 
(*********WAVE NORMAL IS VERTICAL,SO AZIMUTH DIRECTION CANNOT BE CALCULATED *)
 
 ANA=ANGA-ALPH;

       SANA=Sin [ANA];
        SPHI=SANA*STHR/STH;
    CPHI=-Cos[ANA]*Cos[(PHI-PHIR)]+SANA*Sin [PHI-PHIR]*CTHR;
       AZA=180.-Mod[540.-(ArcTan[SPHI,CPHI]-ArcTan[Re[Ray[[6]]] ,Re[Ray[[5]]]])*DEGS,360.];
      Print["uu  ",ii," V=",V," NWHY=",NWHY," H=",H," Range=",RANGE, " AZDEV*AZA=",AZDEV*AZA," El=",EL," ELL=",ELL," SR=", SR," T=",T];

If[ii >1, AppendTo[CoordinatesRay,{ RANGE,H}];];
 

] (* End RapPath *)




(* ****************************************************************************)
(* *************** PolCar   *************************************)
(***   CONVERTS SPHERICAL COORDINATES TO CARTESIAN   ***)

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

PolCarAnalysis[  IENTR_]:=
 Module [{ SINA, COSA,SINP,COSP,TEMP, Done},
Done=False;
If [IENTR==0,
(*  CONVERTS SPHERICAL COORDINATES TO CARTESIAN   *)
If [!(Re[Ray[[5]]]==0.&& Re[Ray[[6]]]==0.),  
VERT=0.;
                 SINA=Sin[Ray[[2]]];
                COSA=Sin[Pi/2-Ray[[2]]];
                SINP=Sin[Ray[[3]]];
                COSP=Sin[Pi/2-Ray[[3]]];
                 XO[[1]]=Ray[[1]]*SINA*COSP;
                 XO[[2]]=Ray[[1]]*SINA*SINP;
                 XO[[3]]=Ray[[1]]*COSA;
    XXX[[4]]=Ray[[4]]*SINA*COSP+Ray[[5]]*COSA*COSP-Ray[[6]]*SINP;
    XXX[[5]]=Ray[[4]]*SINA*SINP+Ray[[5]]*COSA*SINP+Ray[[6]]*COSP;
                 XXX[[6]]=Ray[[4]]*COSA-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]])/Ray[[1]];
                R[[5]]=(XXX[[3]]*(XXX[[1]]*XXX[[4]]+XXX[[2]]*XXX[[5]])-(XXX[[1]]^2 +XXX[[2]]^2)*XXX[[6]])/(Ray[[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[];


];

(***************************************************************************)
(***********************   RAY_PATH                **************************)
(***************************************************************************)
(*rayMatrix.wdx is a record of every ray vector in the ray trace computation*)
RayMatrix =Import["rayMatrix.wdx"];
 
lengthMatrix=Length[RayMatrix]
Print ["Length Matrix=",lengthMatrix];

(*        W.wdx is the parameter database                                    *)
WPost=Import["W.wdx"];

(* Will be  alist of height versus range points to plot                      *)
CoordinatesRay={};


(* 
 * Initialize
 *)
PIT2=2.*Pi
PID2=Pi/2.
DEGS=180./Pi
RAD=Pi/180.
CC=2.997925*10^5
K=2.81785E-15*CC^2/Pi

TLON=WPost[[5]]
TLAT=WPost[[4]]
Print ["XMTRH=",XMTRH," PLAT=",PLAT," PLON=",PLON]
Print [" TLAT=",TLAT," TLON=",TLON]
 
 EARTHR=WPost[[2]];  (* Earth Radious Km  *)
  (* lengthMatrix=362;  *) 
For[ii=1,ii< lengthMatrix, ii++,
            Ray=Re[RayMatrix[[ii]] ];


       RayPath[];




] (* End For Loop  *)
Print["Coordinates =",CoordinatesRay];

 ListPlot[CoordinatesRay,Joined->True]


sfile=OpenWrite["ray_height_versus_range.txt",FormatType->StandardForm]



For[ii=1,ii<lengthMatrix, ii++,
pair=CoordinatesRay[[ii]];
(* Print[pair[[1]]," ", pair[[2]]];*)
Write[sfile,pair[[1]]," ", pair[[2]]];

]
Close[sfile];



sfile=OpenWrite["earth_center_spherical.txt",FormatType->StandardForm]



For[ii=1,ii<lengthMatrix, ii++,
 
 Ray=RayMatrix[[ii]] ;
Write[sfile,Re[Ray[[1]]]," ",Re[Ray[[2]]]," ",Re[Ray[[3]]]];

]
Close[sfile];