(*  Ray Trace Through the IONOSPHERE  Mathematica Notebook                    *)
(*  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 22,2018                                                     *)
(*  http://www.radiocalc.com                                                  *)
(*                                                                            *)
(*  The Mathematica Code 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                                                               *)
 
(* The Code was scanned and OCR'd (back in 2002).                             *) 
(* The corrections took a major effort. Debugged more in December 2018.       *)
(* Obviously a very error prone procedure. The code in the OT Report 75-76    *)
(*  had some issues since it was coded for the CDC Mainframe.                 *)
(* Also the CDC FORTRAN allowed for some non standard variable persistance    *)
(* between subroutine calls.                                                  *)
(* Moreover, some assignments were fixed especially for DRDT(1).              *)
(* A debug option with extensive reporting on computations  was added.        *)
(*                                                                            *)
(* The Mathematica code eliminates all the FORTRAN GO TO statements.          *)
(* Also one gets a better view of global variables.                           *)
(* Mathematica uses color coding on variables (local global).                 *)
(* The Mathematica Modules hopefully clearup what variables are global.       *)
(* Be nice to elliminate most gobal variables                                 *)
(*     (except the W Array) from the code                                     *)
(* The routine POL CAR was modified.                                          *)
(*                                                                            *)
(* The Mathematica Code is a great basis for porting to Java  and other       *)
(* Programming languages ( we have eliminated the Go TO's as stated before    *)
(* plus extensive work for structured programming.                            *)
(*                                                                            *)
(* The effort is to capture the setup used by Jones and Stephenson in the     *)
(* tabulated results.                                                         *)
(*                                                                            *)
(* With Mathematica   vast 3D Graphing capabilties are available as well      *)
(*  as 2-D Graphing.                                                          *)
(* Also easy to access Matrices and Arrays for analysis.                      *)
(*                                                                            *)
(* At any rate getting the excellent work by Jones and Stephenson             *)
(* into Mathematica  opens up a huge opportunity for Research and Development *)
(*  in Ray Tracing through the IONOSPHERE.                                    *)
(*                                                                            *)

(* The GPL Copyright applies to the Mathematica Code.                         *)
(*                                                                            *)
(*          THIS IS AN ALPHA RELEASE                                          *)
(*         Single and Double HOP.                                             *)
(*  Debugging Required for 3rd HOP or even 2nd HOP                             *)
(*                                                                            *)
(*                                                                            *)
(*        Usage Guidance                                                      *)
(* The code is written so one case is analysed.                               *)
(* Set the Elevation and Azimuth and the Transmitter Longitude and Lattitude  *)
(* All   key parameters including above  read through the test file:  ray.dat *)
(* This is the same file used by Jones and Stephenson.                        *)
(* Earth Radius set in the code:                                              *)
(*    EARTHR=6370. ;  W[[2]]=EARTHR;                                          *)
(* The Mathematica code uses the W array as in FORTRAN                        *)
(*   as a database of parameters just as in Jones and Stephenson              *)
(* Set the maximum iterations in the Trace routine loop:                      *)
(*                MAX_ITERATIONS                                              *)
(* Multi Hop different than Jones and Stephenson since we only use the        *)
(* Iterations loop.                                                           *)
(* The results are store in the matrix: RayMatrix                             *)
(* which is stored ina file using Mathematica:                                 *)
(*              Export["RayMatrix.wdx",RayMatrix];                            *)
(* This way the ray information (state vector) is available for Post Processing *)
(* See the WEB Site https://www.radiocalc.com  for update on                 *)
(*   Post Processing notebooks including 3-D Graphics.                        *)
(* We have also translated the Subroutine PRINTR to Mathematica code          *)
(* There are some issus with Fortran ATAN2 and Mathematica ArcTan           *)
(* so this is a work in progress.                                             *)
(*                                                                            *)
(*  The Mathematica Notebook implements one  case out of many possible models *)
(* as explained in Jones and Stephenson.                                      *)
(* Hopefully, by following the Mathematica code, one can substitute various   *)
(* models as outlined in Jones and Stephenson.                                *)


 
(* ****************************************************************************)

 


(* Hamiltonian and Derivatives *)

(* Ray Vector  *)









(* ****************************************************************************)
(* ********************      START              ********************************)
(* ****************************************************************************)
(* ****************************************************************************)



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

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

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


];


(* ****************************************************************************)
(*  CHAPMAN LAYER WITH TILTS,RIPPLES9 AND GRADIENTS                             *)

     (*  W(104)=0.5 FOR AN ALPHA-CHAPMAN LAYER                                       *)
(*        =1.0 FOR A BETA-CHAPMAN LAYER                                         *)
(* ****************************************************************************)
(* Globals:
W Array   R Array  
  CC 
pxpr,pxpth
 
*)


Electx[   ]:=
 Module [{EARTHR, FC,F,ALPHA,A,CC,PERT,HM,HMAX, THETA,THETA2,H,SH, D,B,temp,PID2, E,EXZ},
               
(* Print["Electx Enter"]; *)
PID2=Pi/2;
                 EARTHR=W[[2]];
                  FC=W[[101]];
                   F=W[[7]];
                   ALPHA=W[[104]];
                    A=W[[105]];
                    CC=W[[107]];

                    PERT=W[[150]];


                THETA=R[[2]];
                E=W[[108]]; 
                THETA2=THETA-PID2;
                HM=W[[102]];
                HMAX=HM+EARTHR*E*THETA2;
                H=R[[1]]-EARTHR;
                SH=W[[103]];   (* e.g. 62.0 scale height,km *)
                Z=(H-HMAX)/SH;
                 D=0.;
                 B=W[[106]];
(* Print["Electx Befor Ib B  X=",X]; *)
                 If[B ≠ 0., D=PID2/B];


                  temp=1.+A*Sin[D*THETA2]+CC*THETA2;
                  EXZ=1.-Exp[-Z];
                   X=(FC/F)^2*temp*Exp[ALPHA*(EXZ-Z)];


(*Print["Electx X=",X]; *)
                    pxpr=-ALPHA*X*EXZ/SH;
                    pxpth =X*(D*A*Sin[PID2-D*THETA2]+CC)/temp-pxpr*EARTHR*E;
                    (* IF (PERT.NE.0.) CALL ELECT1   *)

                    (*  Print["Electx 7777 X=",X,"  Z= ",Z, "R (1) =", R[[1]]," H=",H]; *)
                     (* Print["Electx 3333 pxpr=",pxpr,"  pxpth= ",pxpth]; *)

]
(* ****************************************************************************)
(* *************** Program DIPOLY MAGY   *************************************)
(* ****************************************************************************)

(* Globals:
W Array   R Array  
EARTHR, 
    pyrpr,pyrpt,pytpr,pytpt,pypr,pypth

YR, YTH
*)


MAGY[    ]:=
 Module [{EARTHR,term9,PID2,SINTH, COSTH,T1,F,FH   },
               PID2=Pi/2;
(* Print["Entering MAGY "]; *)
F=W[[7]];
FH=W[[201]];
EARTHR=W[[2]];

               SINTH=Sin[R[[2]]];

                   COSTH=Sin[PID2-R[[2]]];
term9=Sqrt[1.+3.*COSTH^2];
T1=FH*(EARTHR/R[[1]])^3/F;



Y=T1*term9;

(* Print["MAGY FH=",FH," Y=",Y, "TERM9=",term9,"T1=",T1]; *)

YR=2.*T1*COSTH;
YTH=T1*SINTH;
     pyrpr=-3.*YR/R[[1]];
    pyrpt =-2.*YTH;
    pytpr=-3.*YTH/R[[1]];
   pytpt=.5*YR;
  pypr =-3.*Y/R[[1]];
 pypth=-3.*Y*SINTH*COSTH/term9^2;


]

(* ****************************************************************************)
(* *************** Program EXPZ2  ENTRY COLFRZ *************************************)
(* ****************************************************************************)
(*  COLLISION FREQUENCY PROFILE FROM TWO EXPONENTIALS   ****************)
(* Globals:
W Array   R Array  
  
   Z,pzpr
*)


COLFRZ[    ]:=
 Module [{EARTHR,PIT2,NU1, H1,T1,H2,F,A1,NU2, A2,H,EXP1,EXP2 },
                
PIT2=2*Pi;
EARTHR=W[[2]];
F=W[[7]];
NU1=N[W[[251]]];
H1=W[[252]];
H2=W[[255]];
A1=W[[253]];
NU2=W[[254]];
A2=W[[256]];



H=R[[1]]-EARTHR;
     EXP1=NU1*Exp[-A1*(H-H1)];
     EXP2=NU2*Exp[-A2*(H-H2)];
     Z=(EXP1+EXP2)/(PIT2*F*1.*10^6);

     pzpr=(-A1*EXP1-A2*EXP2)/(PIT2*F*1*10^6);
(* Print["COLFRZ In H and H1 H2 ",H,"   ", H1, "   ",H2, "  NU1=   ",NU1," Nu2= ",NU2]; *)
(* Print["COLFRZ In Z=",Z," EXP1= ",EXP1, " EXP2= ",EXP2," A1 and A2 ", A1, A2, "  pzpr=",pzpr]; *)


]




(* ****************************************************************************)
(* *************** Program AHWFWC RINDEX   *************************************)
(* ****************************************************************************)
(* **   CALCULATES THE REFRACTIVE INDEX AND ITS GRADIENT USING THE  *************)
(* APPLETON-HARTREE FORMULA WITH FIELD WITH COLLISIONS                         *)
(* Globals:
W Array   R Array  
EARTHR,
   CC 
    pyrpr,pyrpt,pytpr,pytpt,pypr,pypth
Computed in MAGY:  YR, YTH
    COLL,FIELD
    SPACE
    KAY2,H

phpt,phpr,phpth, phpph,phpom,phpkr,phpkth,phpkph,kphpk,POLAR,LPOLAR
 pypr,pypth,pyppa,pyrpr,pyrpt,pyrpp,pytpr,
      pytpt,pytpp,YPH,pyppr,pyppt,pyppp
 pxpr,pxpth,pxpph,pxpt,

   X,Y,Z
  YR,YTH,HMAX
   N2 
*)
RINDEX[    ]:=
 Module [{RAY,OM,OM2,VR,C2,  K2,V2, F,VDOTY,YLV ,YL2 ,YT2,YT4 
     ,UX,UX2,U,D,D2, RAD,scale,GAM,nnp},
   (* Print["Entering RINDEX"]; *)

RAY=W[[1]];
COLL=1.;
       FIELD=1.;
 
KAY2I=0.0;
pxpr=0.;
pxpth=0.;
pxpph= 0.;
pxpt=0.;
pypr=0.;
pypth=0.;
pypph=0.;
pyrpr=0.;
pyrpt=0.;
pyrpp=0.;
pytpr=0.;
pytpt=0.;
pytpt=0.;
pytpt=0.;
pytpp=0.;
pyppr=0.;
pyppt=0.;
pyppp=0.;
pzpr= 0.0;
 pzpth=0.;
pzpph=0.;
ABSLIM=1.*10^(-5);

YPH=0;

 


F=W[[7]];
OM=2*Pi*1.*10^6*F;
 

     C2=CC*CC;
      K2=KR*KR+KTH*KTH+KPH*KPH;

     OM2=OM*OM;
    VR=CC/OM*KR;
    VTH=CC/OM*KTH;
    VPH=CC/OM*KPH;

(* Print["RINDEX F,CC,C2,OM,K2,VR,VTH,VPH:",F,", ",CC,", ",C2", ",OM", ",K2", ",VR", ",VTH", ",VPH];
*)

Electx[];
 
               MAGY[];
 

(* Print["RINDEX after Electx Z=",Z];
*)

V2=VR^2+VTH^2+VPH^2;

(*Print["RINDEX  V2 ... ",V2];
*)

     VDOTY=VR*YR+VTH*YTH+VPH*YPH;
     YLV=VDOTY/V2;
     YL2=VDOTY^2/V2;
      YT2=Y^2-YL2;
      YT4=YT2*YT2;
(* Print["RINDEX  V2 ... ",V2,", ",VDOTY,", ",YLV,", ",YL2,", ",YT2,", ",YT4 ]; *)

 COLFRZ[];
(* 
   Print["COLFRZ"];
*)
(* Note Z modified in COLFRZ *)
(* pzpr calculated in COLFRZ *)
(*Print["RINDEX000 X  =",X," Y= ",Y," Z= ",Z ];
Print["RINDEX Z=",Z];
*)
      U= (1.-Z *I);
     UX=U-X;
     UX2=UX*UX;
     RAD=RAY*Sqrt[YT4+4.*YL2*UX2];

      D=2.*U*UX-YT2+RAD;
       D2=D*D;
       N2=1.-2.*X*UX/D;

(* Print["RINDEX U,UX,UX2,RAD,D,N2",U,", ",UX,", ",UX2,", ",RAD,", ",D,",N2= ",N2 ];

Print["RINDEX X Y Z Y=",X,"  ",Y,"  ",Z,"  ",Y];
*)

pnpps=2.*X*UX*(-1.+(YT2-2.*UX2)/RAD)/D2;
ppspr=YL2/Y*pypr-(VR*pyrpr+VTH*pytpr+VPH*pyppr)*YLV;
ppspth=YL2/Y*pypth-(VR*pyrpt+VTH*pytpt+VPH*pyppt)*YLV;
ppspph=YL2/Y*pypph-(VR*pyrpp+VTH*pytpp+VPH*pyppp)*YLV;


(* Print["RINDEX X222 RAD=",RAD," D2=",D2, " V2=",V2, " UX =",UX];
*)

pnpx=-(2.*U*UX2-YT2*(U-2.*X)+(YT4*(U-2.*X)+5 4.*YL2*UX*UX2)/RAD)/D2;
pnpy=2.*X*UX*(-YT2+(YT4+2.*YL2*UX2)/RAD)/(D2*Y);
pnpz=I*X*(-2.*UX2-YT2+YT4/RAD)/D2;
pnpr=pnpx*pxpr+pnpy*pypr+pnpz*pzpr+pnpps*ppspr;
pnpth=pnpx*pxpth +pnpy*pypth+pnpz*pzpth+pnpps*ppspth;
pnpph=pnpx*pxpph+pnpy*pypph+pnpz*pzpph+pnpps*ppspph;
pnpvr=pnpps*(VR*YL2/V2-YLV*YR);
pnpvth=pnpps*(VTH*YL2/V2-YLV*YTH);
pnpvph=pnpps*(VPH*YL2/V2-YLV*YPH);
nnp=N2-(2.*X*pnpx+Y*pnpy+Z*pnpz);
 
pnpt=pnpx*pxpt;

TT1=(Re[N2]==1.);
TT2=(Abs[Im[N2]]< ABSLIM);

 
(* Print["RINDEX 3333000  N2=",N2," TT1=",TT1, "  TT2=",TT2]; *)
 
     SPACE=(Re[N2]==1. &&  Abs[Im[N2]]< ABSLIM);
 
(* Print["RINDEX 33330111  SPACE=",SPACE]; *)
(* Print["RINDEX 3333 CC=",CC," VDOTY=",VDOTY,"  UX=",UX]; *)
 

POLAR=-I*Sqrt[V2]*(-YT2+RAD)/(2.*VDOTY*UX);
     GAM=(-YT2+RAD)/(2.*UX);
     LPOLAR=I*X*Sqrt(YT2)/(UX*(U+GAM));
      KAY2=OM2/C2*N2;
(*
Print["RINDEX 4444 UX=",UX,"  OM2=", OM2, "  K2=",K2,"  KAY2= ",KAY2];
*)
      (*   IF[RSTART.EQ.0.] GO TO 1   *)
If[RSTART!=0.,
(* Print["RINDEX 66666    K2=",K2,"  KAY2= ",KAY2]; *)
      scale=Sqrt [Re[KAY2]/K2];
(* Print["RINDEX 77777    SCALE=",scale ]; *)

      KR=scale*KR;
     KTH=scale*KTH;
     KPH=scale*KPH;
];
(*  1 CONTINUE  *)
(* *********CALCULATES A HAMILTONIAN H   *)
      Hamiltonian= (C2*K2/OM2-N2)*0.5;
(**********AND ITS PARTTAL DERIVATIVES WITH RESPECT TO  *)
(********TIME,R,THETA,PHI,OMEGA,KR,KTHETA,AND KPHI.     *)
             phpt=-pnpt;
             phpr=-pnpr;
             phpth=-pnpth;
              phpph=-pnpph;
               phpom=-nnp/OM;

phpkr=C2/OM2*KR-CC/OM*pnpvr;

phpkth=C2/OM2*KTH-CC/OM*pnpvth;
     phpkph=C2/OM2*KPH-CC/OM*pnpvph;
     kphpk=N2;


(* Print ["RINDEX pxpr  pxpth  pxpph  pxpt  pypr  pypth  pypph  pyrpr  pyrpt  pyrpp  pytpr  pytpt  pytpt  pytpt  pytpp  pyppr  pyppt  pyppp  pzpr  pzpth  pzpph ",pxpr,pxpth,pxpph,pxpt,pypr,pypth,pypph,pyrpr,pyrpt,pyrpp,pytpr,pytpt,pytpt,pytpt,pytpp,pyppr,pyppt,pyppp,pzpr,pzpth,pzpph]; *)

(* Print["RINDEX Returning"]; *)


]

(*****************************************************************************)
(* ****************************************************************************)
(* *************** Procedure  RKAMExit Utility   *******************************)
(* ****************************************************************************)
RKAMExit[ ]:=
 Module [{ SSE,ii,KK,K,Z,DEL,E   },
          (***************** GOTO 42   EXIT   ***************************)
(* Print["RKAM  GOTO 42 EEXXXIIITTTT Procedure "]; *)
           LL=2;
                        MM=4;
           For [K=1, K≤ 3,K++,
                        XV[[K]]=XV[[K+1]];
                       For[ii=1, ii≤ NN,ii++,
 
                                                                 FV[[K,ii]]=FV[[K+1,ii]];
                                                                 YU [[K,ii]]=YU [[K+1,ii]];
                      ];
          ];
                                     XV[[4]]=XV[[5]];
            For[ii=1, ii≤ NN,ii++,
 
                                                     FV[[4,ii]]=DRDT[[ii]];
                                                     YU[[4,ii]]=YU[[5,ii]];
          ];
           (* Print["RKAM exit 111"]; *)
           If[MODE  ≤ 2, Return[]];
                                   E=Abs[XV[[4]]-ALPHA];
                 (* Print["RKAM exit 222"]; *)
           If [ E ≤ EPM, 
                       Print [" RKAM: GOTO 2000 horror need reentry rewritten code "] 
         ];
(* Print["RKAM exit 333"]; *)
             EPM=E;
              Return[];
      

];
(*****************************************************************************)

(* ****************************************************************************)
(* *************** Program RKAM   *************************************)
(* ****************************************************************************)
(*******  NUMERICAL INTEGRATION OF DIFFERENTIAL EQUATIONS           ****)
(*   GLOBAL
    NN  See NITIAL 
     Y List 20 elements
     T
     STEP
     DRDT  List 20 elements
     MODE
      
     SPACE
     E1MAX,E1MIN,E2MAX,E2MIN
     FACT


*)
 
RKAM[ ]:=
 Module [{ SSE,ii,KKFlag,flag1,flag2000,K,Z,DEL,E,EPSIL   },
(* Print["Entered RKAM"]; *)
flag1=False;
flag2000=False;
If[RSTART == 1,
(* Print["Entered RSTART 1 NN=",NN]; *)
flag1=True;
      LL=1;
                            MM=1;
                           If [MODE==1, MM=4];
                          ALPHA=T;
                          EPM=0.0;
                          BETT[[2]]=0.5;
                           BETT[[1]]=BETT[[2]];
                          BETT[[3]]=1.0;
                         BETT[[4]]=0.0;

(* Print["RKAM RESTART BET =", BETT[[1]]," ",BETT[[2]], "  ", BETT[[3]], "  ", BETT[[4]]]; *)
                         (* Was STEP=SPACE in FORTRAN which is rediculous *) 

                   STEP=1.0;
                  If[SPACE== False , STEP=0];
                  STEP=1.0;

                         RR=19.0/270.0;
                          XV[[MM]]=T;
                          If[E1MIN < 0., E1MIN=E1MAX/55. ];
                          If[FACT ≤ 0., FACT=0.5 ]


                            HAMLTN[];

      For [ii=1,ii<=NN,ii++,
            FV[[MM,ii]]=DRDT[[ii]];
                                     YU[[MM,ii]]=R[[ii]];
(* Print["RKAM RESTART ii=",ii," Y(ii)=",R[[ii]]," FV[[MM,ii]]=",FV[[MM,ii]]]; *)

      ];

                  (*  Print["RKAM RESTART 2"]; *)
                               RSTART=0.;
(* 1000 IF(MODE.NE.1) GO TO 2000  *) (* FIX THIS *)
];(* RSTART==0  *)

If[flag1==False,

 If[MODE ≠ 1,
          (*  Print["RKAM GOTO 2000"]; *)
              flag2000=True;
     ];
];
(* ELSE RESTART==1  *)
(* All important 1000   *)
While[True,

If[flag2000== False ,
 (*
     * RUNGE-KUTTA
  *)

(* 1001 *)
       For[  K=1,K ≤ 4   ,K++,(* 1034*)
           For [ii=1,ii≤ NN,ii++,(* 1350*)
           
                                               DELY[[K,ii]]=STEP*FV[[MM,ii]];
                                                Z=YU[[MM,ii]];
(*1350*)                             R[[ii]]=Z+BETT[[K]]*DELY[[K,ii]];
                                 (* Print["RKAM 1350 R[[ii]]=",R[[ii]]," Z= ", Z, " DELY[[K,ii]]=",DELY[[K,ii]], " FV[[MM,ii]]=",FV[[MM,ii]], " STEP=",STEP]; *)
            ];
                                       T=BETT[[K]]*STEP+XV[[MM]];
 
               HAMLTN[ ];
               For[ii=1, ii≤NN,ii++,
                                                 FV[[MM,ii]]=DRDT[[ii]];
             ];
        ]; (* 1034*)
 (* Print["RKAM 1034"]; *)

                                For[ii=1, ii≤ NN,ii++,
                    DEL=(DELY[[1,ii]]+2.0*DELY[[2,ii]]+2.0*DELY[[3,ii]]+DELY[[4,ii]])/6.0;
                                            YU[[MM+1,ii]]=YU[[MM,ii]]+DEL;
          ]; (*1039*)
 (* Print["RKAM 1039"]; *)
                                  MM=MM+1;
                                  XV[[MM]]=XV[[MM-1]]+STEP;
            For[ii=1, ii≤ NN,ii++,
                     R[[ii]]=YU[[MM,ii]];
                    ];
             T=XV[[MM]];
              HAMLTN[ ];
(*   IF(MODE.EQ.1) GO TO 42   *) (* FIX THIS *)
              If[ MODE==1,
                     (***************** GOTO 42   EXIT   ***************************)
                        RKAMExit[];
                         Return[];
                ];

                                          (*****************   END EXIT   ***********************)
For[ii=1, ii≤ NN,ii++,
  FV[[MM,ii]]=DRDT[[ii]];
];

(* Print["RKAM 55555 MM=",MM]; *)



(*   END WHILE    If[Not[MM <= 3],Break[]]  *)
]; (* flag2000  *)
flag1=False;
If [MM > 3  || flag2000==True ,
(* Print["RKAM MM >3 MM=",MM]; *)
If [flag2000 == True, abcd=1(* Print["RKAM GOTO 2000 landed "];*) ];
flag1=True;
 (*********                         ******)
           (*********   ADAMS-MOULTON         ******)
 (*********                         ******)
 (* 2000 *)

For[ii=1, ii≤ NN,ii++,
       DEL=STEP*(55.*FV[[4,ii]]-59.0*FV[[3,ii]]+37.*FV[[2,ii]]-9.*FV[[1,ii]])/24.;
                                R[[ii]]=YU[[4,ii]]+DEL;
                                 (*2048*)    DELY[[1,ii]]=R[[ii]]
];
 T=XV[[4]]+STEP;
 HAMLTN[ ];
XV[[5]]=T;

         For[ii=1, ii≤ NN,ii++,
                                  DEL=STEP*(9.*DRDT[[ii]]+19.*FV[[4,ii]]-5.*FV[[3,ii]]+FV[[2,ii]])/24.;
                                 YU[[5,ii]]=YU[[4,ii]]+DEL;
                                  R[[ii]]=YU[[5,ii]];

 
];


HAMLTN[ ];
(* Print["RKAM 66666   "]; *)
 

If[MODE ≤ 2,
(* Print["RKAM 66666 MODE ≤ 2 "]; *)
                          RKAMExit[];
                         Return[];
            ];

];  (* End If [MM > 3   *)
If[flag1 == True,
(* Print["RKAM ERROR ANALYSIS flag1 true   "]; *)
(*
    * ERROR ANALYSIS
 *)
                        SSE=0.0;
                        For[ii=1, ii≤ NN,ii++,
                                   EPSIL=RR*Abs[R[[ii]]-DELY [[1,ii]]];
           If[ MODE== 3 && R[[ii]] ≠ 0.0 ,
                       EPSIL=EPSIL/Abs[R[[ii]]] ;
          ];
(* Print["RKAM RKAM ERROR ANALYSIS ii=",ii,"  R[[ii]=",R[[ii]],"  DELY [[1,ii]]= ",DELY [[1,ii]],"    EPSIL=", EPSIL ," SSE= ",SSE, " RR=", RR];*)
           If [SSE< EPSIL , SSE=EPSIL];

           
         ];
             
  
(* 3033*)
(* Print["RKAM RKAM ERROR ANALYSIS STEP=",STEP, "  SSE= ",SSE," EPSIL= ", EPSIL, " MODE= ", MODE]; *)
KKFlag=False;
If [E1MAX ≤ SSE,
 (* Print["RKAM 8888 Error Analysis STEP=",STEP,"E2MIN= ", E2MIN ]; *)
If[Abs[STEP] ≤ E2MIN, 
      (*****************GOTO 42    EXIT   ***************************)
                 RKAMExit[];
                         Return[];
 ];

     (* Print["RKAM 776677 "]; *)
              (*****************   END EXIT   ***********************)

KKFlag=True;
LL=1;
MM=1;
STEP=STEP*FACT;

];


(* GOTO 1001 *)
 (* Print["RKAM  GOTO 1001 "]; *)
If [KKFlag== False, 
    If [LL≤ 1|| Abs[SSE ]≥ E1MIN ||E2MAX <=Abs[STEP] ,
                   (*GOTO 42 *)
          (***************** GOTO 42   EXIT   ***************************)
                         RKAMExit[];
                         Return[];
  ];
     

        (*****************   END EXIT   ***********************)
             
 
  
   LL=2;
                      MM=3;
                     XV[[2]]=XV[[3]];
                      XV[[3]]=XV[[5]];
                      For[ii=1, ii≤ NN,ii++,

                                    FV[[2,ii]]=FV[[3,ii]];
                                    FV [[3,ii]]=DRDT[[ii]];
                                    YU[[2,ii]]=YU[[3,ii]];
                                     YU[[3,ii]]=YU[[5,ii]];
                        ];
     STEP=2.0*STEP;
]   (*  KKFlag Flag    *)

]; (* flag1 *)
 (* Print["RKAM  End While goto 1000 "]; *)
] (* END WHILE   *)   (* was goto 1000 *)
] (* End RKAM  *)

(* ****************************************************************************)
(* *************** Program HAMLTN   *************************************)
(*********   ********CALCULATES HAMILTONS EQUATIONS FOR RAY TRACING   ********)
(* Globals:
W Array R Array
RAD, K, CC
COLL,FIELD,SPACE,KAY2,KAY2I,H,HI,
phpt,phpti,phpr,phpri,phpth,phpthi,phpph,phpphi,phpom,phpomi,phpkr,phpkri,phpkth,phpkti,phpkph,phpkpi,kphpk,kphpk
POLAR,POLARI,LPOLAR,LPOLRI
R(20),T,STP,DRDT(20)


*)
(* ****************************************************************************)
HAMLTN[   ]:=
 Module [{TH,PH, DTHDT,DPHDT,DKRDT,DKTHDT,DKPHDT,F,
    OM,STH, CTH,RSTH,RCTH,NR,DRDT1},
TH=R[[2]];
PH=R[[3]];
KR=R[[4]];
KTH=R[[5]];
KPH= R[[6]];
DTHDT=DRDT[[2]];
DPHDT=DRDT[[3]];
DKRDT=DRDT[[4]];
DKTHDT=DRDT[[5]];
DKPHDT=DRDT[[6]];
F=W[[7]];

OM=Pi*2*10^6*F;
STH=Sin[TH];
     CTH=Sin[Pi/2-TH];
      RSTH=R[[1]]*STH;
     RCTH=R[[1]]*CTH;
 
RINDEX[  ];

DRDT1=- phpkr/( phpom  *CC);
DRDT[[1]]=DRDT1;

DTHDT=-phpkth/(phpom*R[[1]]*CC);
     DPHDT=-phpkph/(phpom*RSTH*CC);
     DKRDT=phpr/(phpom*CC)+KTH*DTHDT+KPH*STH*DPHDT;
     DKTHDT=(phpth/(phpom*CC)-KTH*DRDT1+KPH*RCTH*DPHDT)/R[[1]];
      DKPHDT=(phpph/(phpom*CC)-KPH*STH*DRDT1-KPH*RCTH*DTHDT)/RSTH;
      NR=6;
      (***   PHASE PATH  ***)
       If [W[[57]] ≠  0.,
NR=NR+1;
DRDT[[NR]]=-kphpk/phpom/OM;

] ;
(***   ABSORPTION  ***)
If [W[[58]] ≠  0.,
      NR=NR+1;
             DRDT[[NR]]=10./Log[10]*kphpk*KAY2I/(KR*KR+KTH*KTH+KPH*KPH)/phpom/CC;
     ];
(***   DOPPLER SHIFT  ***)
If [W[[59]]≠ 0, 
    NR=NR+1;
           DRDT[[NR]]=-phpt/phpom/CC/(2*Pi)
];
(***   GEOMETRICAL PATH LENGTH  ***)
NR=NR+1;
     DRDT[[NR]]=-Sqrt(phpkr^2+phpkth^2+phpkph^2)/phpom/CC;


 DRDT[[2]]=DTHDT;
 DRDT[[3]]=DPHDT;
 DRDT[[4]]=DKRDT;
 DRDT[[5]]=DKTHDT;
 DRDT[[6]]=DKPHDT;


(***   OTHER CALCULATIONS  ***)
(* Print["HAMLTN Done"];*)
]

(* ****************************************************************************)
(* *************** Program PRINTR   *************************************)
(* ****************************************************************************)

PRINTR [ theText_,Val_]:=
 Module [{nnn},

Print["RAY TRACE(444444): ",theText];

]

(* ****************************************************************************)
(* *************** Program REACH   *************************************)
(* ****************************************************************************)

(* CALCULATES THE STRAIGHT LINE.RAY PATH BETWEEN THE EARTH      *)

(* AND THE IONOSPHERE OR BETWEEN IONOSPHERIC LAYERS             *)
(* 
* Globals
*)
(*
     /RK/  NN,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART
    /TRAC/    GROUND,PERIGE,THERE,MINDIS,NEWRAY,SMT
     /COORD/    S
    /RIN/     COLL,FIELD,SPACE,N2,N2I,PNP(10),POLAR,LPOLAR
     /XX/ 
   X,pxpr,pxpth,pxpph,pxpt,HMAX
      R(20),T,STP,DRDT(20)/WW/ID(10),W0,W(400)
*)
 
(* Reach25[ S1_,EARTHR_ ,NN_ ,CROSSG_,CROSS_,CROSSR_,SG_,SR_]:= *)
Reach25[]:=
 Module [{nnn,V   },

(* Print[" Entered Reach25"]; *)
 



If[CROSS ,   S=Min [S,S1]; 
];
(*** CONVERT POSITION AND DIRECTION TO SPHERICAL POLAR,COORDINATES  *)
(*  AT A DISTANCE S ALONG THE RAY  *)
PolCar[1];
R[[1]]=Max[Re[R[[1]]],EARTHR];
If[NN≥ 7,
    For [nnn=7, nnn ≤NN, nnn++,
            R[[nnn]]=R[[nnn]]+S*DRDT[[nnn]];
];
];
T=T+S;
RINDEX[   ];
(* AT A PERIGEE  *)
PERIGE=(S==(-UP));
(*  CORRECT MINOR,ERRORS   **)
If [PERIGE, R[[4]]=0.
];
(*   ****KEEP CONSISTENCY AFTER CORRECTING MINOR ERRORS  *)
DRDT[[1]]=R[[4]];
(***ON THE GROUND   ***)
GROUND=(S==SG && CROSSG);
 
Print["Reach S=",S," SG=", SG, " CROSSG =",CROSSG, " GROUND=", GROUND];
Print["Reach S1=",S1," SR= ", SR, " CROSS=", CROSS, " SPACE=",SPACE];
(* AT THE RECEIVER HEIGHT  *)
THERE==(S==SR &&  CROSSR && !PERIGE);
(* **AT A CLOSEST APPROACH TO THE RECEIVER HEIGHT  *)
MINDIS=(PERIGE && S==SR && CROSSR);
       RSPACE=SPACE;
        V=Sqrt[N2/(R[[4]]^2+R[[5]]^2+R[[6]]^2)];
(*** RENORMALIZE THE WAVE NORMAL DIRECTION TO SQPT(REAL(N**2))**)

R[[4]]=R[[4]]*V;
       R[[5]]=R[[5]]*V;
       R[[6]]=R[[6]]*V;

If [ !SPACE,
         PRINTR["ENTR_ION",0.];
];

        RSTART=1;
  
(* Print["Reach25 Space=",SPACE," PERIGE=",PERIGE, " GROUND=",GROUND," THERE=",THERE," MINDIS=",MINDIS," RSPACE=", RSPACE]; *)


];

REACH[ ]:=
 (*Module [{EARTHR,XMTRH,RCVRH,NSTEP,H,V,RADG,DISTG,RADR,DISTR,S1,SR,SG,
   CROSSG,CROSSR,CROSS },*)
  Module [{EARTHR,XMTRH,RCVRH,NSTEP,H,V,RADG,DISTG,RADR,DISTR },

EARTHR=W[[2]];
 XMTRH=W[[3]];
RCVRH=W[[20]];
 NSTEP=500;
VERT=0.0;
(* Print["REACH  Enter"]; *)
HAMLTN[ ];
(* Print["REACH  Post HAMLTN  RSPACE=",RSPACE]; *)
H=Re[R[[1]]]-EARTHR;
(****   IF (.NOT.NEWRAY.AND..NOT.RSPACE) CALL PRINTR(8HEXIT_ION,0.)  **)
If [!NEWRAY&&!RSPACE,
       PRINTR["EXIT_ION",0.];
];

NEWRAY= False;

V=Sqrt[Re[R[[4]]]^2+Re[R[[5]]]^2+Re[R[[6]]]^2];
(**********NORMALIZE THE WAVE NORMAL DIRECTION TO ONE **)
R[[4]]=Re[R[[4]]]/V;
     R[[5]]=Re[R[[5]]]/V;
     R[[6]]=Re[R[[6]]]/V;
(*  NEGATIVE OF DISTANCE ALONG RAY TO CLOSEST APPROACH TO CENTER OF EARTH *)
UP=Re[R[[1]]]*Re[R[[4]]];

     RADG=EARTHR^2-Re[R[[1]]]^2*(Re[R[[5]]]^2+Re[R[[6]]]^2);
(* Print["REACH 111 RADG=",RADG, " R1 = ", R[[1]]]; *)
RADG=Re[RADG];
(* Print["REACH RADG=",RADG, " R1 = ", R[[1]]]; *)

     DISTG=Sqrt[ Max[0., RADG ]];
 
(* DISTANCE ALONG RAY TO FIRST INTERSECTION WITH OR CLOSEST APPROACH TO THE EARTH  *)
SG= -UP-DISTG;
(* CROSSING IS TRUE IF THE RAY WILL INTERSECT OR TOUCH THE EARTH  *)
CROSSG=UP<0  && RADG ≥  0.;

(* Print["REACH 1111 "]; *)

RADR=(EARTHR+RCVRH)^2-R[[1]]^2*(R[[5]]^2+R[[6]]^2);
RADR=Re[RADR];
      DISTR=Sqrt [Max[0., RADR ]];
(*  DISTANCE ALONG RAY TO THE FIRST INTERSECTION WITH   CLOSEST APPROACH TO THE RECEIVER,HEIGHT  *)

 SR=DISTR-UP;
If[UP <0.0 && DISTR < -UP && Re[R[[1]]] ≠ EARTHR+RCVRH,
        SR= -DISTR-UP; 
];
(* Print["REACH 2222 "]; *)
(* CROSSR IS TRUE IF THE RAY WILL INTERSECT WITH OR MAKE A CLOSEST APPROACH TO THE RECEIVER HEIGHT *) 
CROSSR=Re[R[[4]]]< 0. || Re[R[[1]] ]<(EARTHR+RCVRH);
CROSS=CROSSG | CROSSR;
(*  MAXIMUM DISTANCE IN WHICH TO LOOK FOR THE IONOSPHERE *)
S1=Min[SR,SG];
If[ !CROSSG, S1=SR];
If[ UP < 0.0 ,
CROSS=True;

(* IF RAY IS GOING DOWN,S1 IS AT MOST THE DISTANCE TO A PERIGEE   *)
S1=Min[Re[S1],-UP];

];
(*  CONVERT THE POSITION AND DIRECTION OF THE RAY TO CARTESIAN   COORDINATES *)


VERT=0.;
PolCar[0];

SSTEP=100.;
        S=SSTEP;
(* Print["REACH NSTEP=",NSTEP]; *)
 
For [ii=1,  ii≤ NSTEP, ii++,
    (* Print["REACH **********************    ii=",ii,"  S=", S, "  SSTEP=", SSTEP, " SPACE=",SPACE]; *)

    If[(S-SSTEP) >S1 && CROSS,
                                  (* Print["REACH 1111111 Reach25 "]; *)
           (* Reach25[S1,S,EARTHR,NN,CROSSG,CROSS,CROSSR,SG,SR];*)
Reach25[];
            RSTART=1.;

            Return[];
           (* GOTO 25 *)
  ];
PolCar[1];
Electx[];
If[X ≠ 0 ,
     RINDEX[  ];
    (*    ****EFFECTIVELY FREE SPACE  *****     *)
    If[ ! SPACE,
           If [(SSTEP  <= 0.5*10^(-4)),   (* original 0.5*10^(-4)  *)
                    (* Print["REACH Reach25 2222222"]; *)
                     (* Reach25[S1, EARTHR,NN,CROSSG,CROSS,CROSSR,SG,SR];*)
                          Reach25[];
                        RSTART=1.;

(* Print["REACH Reach25 4444444"];*)
                     Return[];
         (* GOTO 25 *)
          ];
             (*  ***RAY IN THE IONOSPHERE.STEP BACK OUT   *)
             S=S-SSTEP;
              (* DECREASE STEP SIZE *)
              SSTEP=SSTEP/10.;
                          
];
   
];
S=S+SSTEP;

]; (* For loop *)

Print[" EXCEEDED,",NSTEP,"   STEPS IN SUBROUTINE REACH"];	

   Abort[];      




];

(* ****************************************************************************)
(* *************** Program BackUp   *************************************)

(*
 * Globals
* /RK/N,STEP,MODE,E1MAX,EIMIN,E2MAX,E2MIN,FACT,RSTART
 * COMMON/FLG/NTYP,NEWWR,NEWWP,PENET,LINES,IHOP,HPUNCH
* /TRAC/GROUND,PERIGE,THERE,MINDIS,NEWRAY,SMT* /RIN/MODRIN(3),COLL,FIELD,SPACE,N2,PNP(10),POLAR,LPOLAR
* /XX/MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX
*R(20),T,STP,DRDT(20)/WW/ID(10),W0,W(400)
*
*)
(* ****************************************************************************)
BackUp[ EARTHR_  ]:=
Module [{H, 
     RSTART},
(*   For RKAM   *)

LL=1;
      MM=1;
      ALPHA=0;
       EPM=0;
        RR=0;
(* End  For RKAM   *)
If[IGRAZE ≠ 1,
   HOME=DRDT[[1]]*(R[[1]]-EARTHR-HS) ≥ 0.;
   If[HS> 0.&& !HOME|| HS==0.&& DRDT[[1]]> 0.,
                 Print["Backup GOTO 30"];
];
(*****FIND NEAREST INTERSECTION OF RAY WITH THE HEIGHT HS  **)
For[ ii=1, ii≤10,ii++,
     STEP=-(R[[1]]-EARTHR-HS)/DRDT[[1]];

                           STEP= Min[Abs[STP],Abs[STEP] ]*Sign[STEP]; (*** CHECK THIS ***)
If [Abs(R[[1]]-EARTHR-HS)< .5*10^(-4)  && Abs[STEP] <1.,
             Print["Backup GOTO 60"];
];
MODE=1;
RSTART=1.;

(* RKAM[DELY,BET,XV,FV,YU,LL,MM,ALPHA,EPM,RR,RSTART]; *)
RKAM[ ];
RSTART=0.;
]  (*10*)

];(* 55 Continue   *)
(*********FIND NEAREST CLOSEST APPROACH OF RAY TO THE HEIGHT HS ***)
THERE=False;
If[!(SMT> Abs[R [[1]]-EARTHR-HS]), 

For [ii=1, ii≤10,ii++,       (*20*)
STEP=-R[[4]]/DRDT[[4]];
 
STEP= Min[Abs[STP],Abs[STEP] ]*Sign[STEP]; (*** CHECK THIS ***)
 If [Abs[R[[4]]] ≤ 1.*10^(-6)  && Abs[STEP]< 1.,
  (* GO TO 60 *)
       MODE=IntegrationType;
      STEP=STEP1;
      Return[];
];
(* GRAZE *)
MODE=1;
            RSTART=1.;

(*RKAM[DELY,BET,XV,FV,YU,LL,MM,ALPHA,EPM,RR,RSTART]; *)
RKAM[ ];
RSTART=1.;
 If [DRDT[[4]]*(R[[1]]-EARTHR-HS)< 0.,
     (*  GO TO 30 *)
   Break; (* Break out of for loop *)
];
If[R[[5]]==0.&& R[[6]]==0.,
  (* GO TO 60 *)
       MODE=IntegrationType;
      STEP=STEP1;
      Return[];
];
];   (* 20 CONTINUE  *)

]; (* GO TO 30 *) 
(* BACKUP *)
MODE=1;
(* *******ESTIMATE DISTANCE TO NEAREST INTERSECTION OF RAY WITH HEIGHT  *)
(* *******HS BEHIND THE PRESENT RAY POINT  *)
H=R[[1]]-EARTHR;
STEP=(-R[[4]]-Sqrt[R[[4]]^2-2.*(H-HS)*R[[4]]/DRDT[[1]]*DRDT[[4]]])/DRDT[[4]];
RSTART=1.;
(*RKAM[DELY,BET,XV,FV,YU,LL,MM,ALPHA,EPM,RR,RSTART]; *)
        RKAM[ ];
        RSTART=1.;
(***** ****FIND NEAREST INTERSECTION OF RAY WITH HEIGHT HS     ***)
For[ ii=1, ii≤10,ii++,  (*40*)
       STEP=-(R[[1]]-EARTHR-HS)/DRDT[[1]];
       STEP= Min[Abs[STP],Abs[STEP] ]*Sign[STEP]; (*** CHECK THIS ***)
       If [Abs[R[[1]]-EARTHR-HS]<.5*10^(-4) && Abs[STEP]< 1., 
              (* GO TO 60 *)
               MODE=IntegrationType;
                STEP=STEP1;
                Return[];

     ];

(* BACK UP *)
RSTART=1.;
(* RKAM[DELY,BET,XV,FV,YU,LL,MM,ALPHA,EPM,RR,RSTART]; *)
        RKAM[ ];
        RSTART=1.;



]; (* 40*)
THERE=True;
    MODE=IntegrationType;
                STEP=STEP1;
                Return[];
];

(* ****************************************************************************)
(* *************** Program Trace   *************************************)

(*
 * Globals
* /RK/N,STEP,MODE,E1MAX,EIMIN,E2MAX,E2MIN,FACT,RSTART
 * COMMON/FLG/NTYP,NEWWR,NEWWP,PENET,LINES,IHOP,HPUNCH
* /TRAC/GROUND,PERIGE,THERE,MINDIS,NEWRAY,SMT* /RIN/MODRIN(3),COLL,FIELD,SPACE,N2,PNP(10),POLAR,LPOLAR
* /XX/MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX
*R(20),T,STP,DRDT(20)/WW/ID(10),W0,W(400)
*
*)
(* ****************************************************************************)
TraceRT[   ]:=
 Module [{EARTHR,RCVRH,HOP,MAXSTP,SKIP,RAYSET,PLT,NHOP,MAX,NSKIP,
   IHOP,HTMAX,APHT,H,j,ll,TOLD,WAS,WASNT,TMP,
     skipThese },

(* DIMENSION DELY(4,20),BET(4),XV(5),FV(4,20),YU(5,20)
DOUBLE PRECISION YU
*)
(*   For RKAM   *)
XV= Table[ 0,5];
BETT= Table[ 0,4];
DELY= Table[0,4,20];
FV= Table[0,4,20];
YU= Table[0,5,20];
LL=1;
      MM=1;
      ALPHA=0;
       EPM=0;
        RR=0;
(* End  For RKAM   *)

RSPACE=False;


EARTHR=W[[2]];
RCVRH=W[[20]];
HOP=W[[22]];
MAXSTP=W[[23]];
 SKIP=W[[71]];
RAYSET=W[[72]];
PLT=W[[81]];

NHOP=HOP;
     MAX=MAXSTP;
    NSKIP=SKIP;
    RSTART=1.;

HAMLTN[ ];
HOME=Re[DRDT[[1]]]*(Re[R[[1]]]-EARTHR-RCVRH) ≥ 0.;
(*  IHOP=O TELLS PRINTR TO PRINT HEADING AND PUNCH A TRANSMITTER *)
(******RAYSET AND TELLS RAYPLT TO START A NEW RAY  *)
IHOP=0;
PRINTR ["XMTR____",0.];
HTMAX=0.;
    NEWRAY=True;
THERE=(Re[R[[1]]]-EARTHR)==RCVRH;
(*10*)  IHOP=IHOP+1;
If [IHOP> NHOP, Return[]];
PENET=False;
APHT=RCVRH;

(* Print["TraceRT MAX=",MAX," SPACE+",SPACE]; *)
(******************************************************)
(*   ****LOOP ON MAXIMUM NUMBER OF STEPS PER HOP       *)
(******************************************************)


(* For [j=1, j≤MAX, j++,*) 
(* For [j=1, j≤362, j++,  *)
For [j=1, j≤603, j++, 
 Print["TraceRT j=",j];
AppendTo[RayMatrix,R];



skipThese=False;
H=Re[R[[1]]-EARTHR];
If [Abs[H-RCVRH] > Abs[APHT-RCVRH], APHT=H];
HTMAX=Max[H,HTMAX];
If [SPACE ,   (* if ! space goto 12 *)

(* Print["TraceRT SPACE=",SPACE]; *)
REACH[];
RSTART=1.;
H=Re[R[[1]]-EARTHR];
If [Abs[H-RCVRH] > Abs[APHT-RCVRH],  APHT=H];
(* Print["TraceRT SPACE 1111 "]; *)


HTMAX=Max[H,HTMAX];
If [SPACE,  (* if ! space goto 12 *)
Print["TraceRT Enter SPACE 2222 "];
If [PERIG,   PRINTR ["PERIGEE_",0.] ];
If [THERE, Print["TraceRT GOTO 51"],
If  [MINDIS, Print["TraceRT  GO TO 40"],
If [GROUND, 
    Print["TraceRT  GO TO 60"]; 
    PRINTR["GRND_REF",RAYSET];
    (* 60*) R[[1]]=EARTHR;
                                    If [Abs[RCVRH] > Abs[APHT-RCVRH] , APHT=0.];
                                    R[[4]]=Abs [Re[R[[4]]]];



                                   DRDT[[1]]=Abs [Re[DRDT[[1]]]];
                                   RSTART=1.;
                                  HPUNCH=HTMAX;
 
                                  HTMAX=0.;
            skipThese=True;
            HOME=True;
             THERE=False;
            H=0;
              PENET=False;
              APHT=RCVRH;

 




 
If[PERIGE,  Print ["TraceRT GO TO 79"],
];
];
];
  Print["TraceRT SPACE 33333 "];
];
Print["TraceRT SPACE 444444 "];
]; 
 Print["TraceRT  GO TO 12 SPACE  ",SPACE];
 ]; (* GO TO 12 *)

 If[skipThese == False,


  (* Print["TraceRT 666666  SPACE=  " , SPACE]; *)
(*12*) For [ll=1,ll≤ NN,ll++,
    ROLD[[ll]]=R[[ll]];
   DROLD[[ll]]=DRDT[[ll]];
];
TOLD=T;
 WAS=THERE;

(*RKAM[DELY,BET,XV,FV,YU,LL,MM,ALPHA,EPM,RR,RSTART]; *)
 
(* Print["TRACE Calling RKAM"]; *)
RKAM[];
 

For[ ll=1, ll≤  NN, ll++,
(* Print["TraceRT After RKAM Check ROLD[[ll]]=",  ROLD[[ll]], " R[[ll]]=",R[[ll]]]; *)
abcd=1;

];


H=Re[R[[1]]]-EARTHR;
        THERE=False;

WASNT=!HOME;

        HOME=(Re[DRDT[[1]]]*(H-RCVRH)) ≥ 0.0;
TMP=(Re[DRDT[[1]]-DROLD[[1]]])*(T-TOLD);

       SMT=0.;
      If [TMP ≠ 0., SMT=Re[0.5*(R[[1]]-ROLD[[1]]+0.5*TMP)^2/Abs[TMP]]
];
If [(((H-RCVRH)*(Re[ROLD[[1]]]-EARTHR-RCVRH)< 0.) && !WAS)||  (WAS&& (Re[DRDT[[1]]]*Re[ROLD[[1]]]<0.) && HOME)  ,Print["TraceRT GO TO 50  XXXX"];
];
       If[HOME&& WASNT, Print["TraceRT  GO TO 30"]
];
If [H<0.||Re[DRDT[[1]]]>0.&& Re[DROLD[[1]]]<0.&& SMT>H, Print ["TraceRT GO TO 20"]
];
 
If [Re[DROLD[[1]]] ≤ 0.&& Re[DRDT[[1]]]> 0.,  PRINTR["PERIGEE_",0.]];

If [Re[DROLD[[1]]]> 0.&& Re[DRDT[[1]]]< 0.,   PRINTR["APOGEE__",0.]];

If [Re[DROLD[[2]]]*Re[DRDT[[2]]]< 0.,   PRINTR["MAX_LAT_",0.]];
       If[Re[DROLD[[3]]]*Re[DRDT[[3]]]<0.,     PRINTR["MAX_LONG",0.]];

For[ ii=4, ii≤ 6, ii++,
    (* If [ROLD[[ii]]*R[[ii]] <0.,   PRINTR["WAVE_REV",0.]]; *)
If [Re[ROLD[[ii]]]*Re[R[[ii]]] <0.,   PRINTR["WAVE_REV",0.]];
(* Print["TRACE WAVE REV ROLD[ii]=",ROLD[[ii]], " R[[ii]] =",R[[ii]]];*)
(* If [ROLD[[ii]]*R[[ii]] <0.,   PRINTR["WAVE_REV",0.]]; *)
             abcd=1;
             ];
(*75 *)  
        If [H> HMAX && H> RCVRH && Re[DRDT[[1]]]>0.,
(*  GO TO 90 *) 
        PENET=True;  PRINTR["PENETRAT",RAYSET]; Return[]
]; 

(* If [Mod[j,NSKIP]== 0,  PRINTR["_____Mod[j,NSKIP]== 0 ",0.]]; *)
 
]; (* skip these  *)
]; (* end for loop  steps  79   *)
(**** *******EXCEEDED MAXIMUM NUMBER OF STEPS   ******)
 
  PRINTR["STEP_MAX",RAYSET];

 

 Export["RayMatrix.wdx",RayMatrix];  

(* Print[RayMatrix]; *)


Abort[];
   Return[];




] (* end module *)

(* ****************************************************************************)
(* *************** Program START nitial  ********************************************)
(* ****************************************************************************)
(*   *********CONSTANTS   ****)
(* PI=3.1415926536  *) (* Original *)
PIT2=2.*Pi
PID2=Pi/2.
DEGS=180./Pi
RAD=Pi/180.
CC=2.997925*10^5
K=2.81785E-15*CC^2/Pi
LOGTEN=Log[10.]
(**********INITIALIZE SOME VARIABLES     ***)
 
PLON=0.
PLAT=PID2
 
(**********INITIALIZE SOME VARIABLES IN THE W ARRAY   ***)

ClearAll;

WW=Import["~/Silicon_DSP_Repositories/IONOSPHERIC_RESEARCH/RayTraceMathematica/ray.csv","CSV"]
wel= WW[[27]]
wel[[3]]
wlength=Length[WW]
 
W= Table[ 0,256]
For[i=1, i<=wlength, i++,
    W[[i]]=0.0
]
For[i=1, i< wlength, i++,
    wel= WW[[i]]; 
Print [i];
    index=wel[[1]];
value=wel[[2]];
Print[ "Index=",index," Value=",value];
W[[index]]=value;
 
    

] 
W

For[i=1, i< 43, i++,
    wel= WW[[i]]; 
Print [i];
    index=wel[[1]];
value=wel[[2]];
Print[ "Index=",index," Value=",value];
W[[index]]=value;
 
    

] 
W
For [ ii=1, ii<=256, ii++,
 
If[W[[ii]] ≠ 0.0 , 
             Print["W[",ii,"]=",W[[ii]]];
            ]
]

F=0.
BETA=0.
AZ1=0.

(*  CHECK THIS 
MODE=-1  (*EXTRAORDINARY  1= "ORDINARY" 0="NO FIELD"  *)
*)

(**********INITIALIZE SOME VARIABLES     ***)
 
 

PLON=W[[25]];
PLAT=W[[24]];

PLON=PLON/DEGS;
PLAT=PLAT/DEGS;


EARTHR=6370. ;  W[[2]]=EARTHR;  (* Earth Radious Km  *)
IntegrationType=3.
W[[41]]=IntegrationType
MAXERR=1.0*10^(-4);    W[[42]]=MAXERR;
ERATIO=50.;  W[[43]]=ERATIO;
STEP1=1. ;   W[[44]]=STEP1;
STPMAX=100.;  W[[45]]=STEPMAX;
STPMIN=1.*10^(-8);  W[[46]]=STPMIN;
FACTR=0.5;   W[[47]]=FACTR;
MAXSTP=1000.;  W[[23]]=MAXSTP;
(*
 * Globale DRDT Vector
 *)
DRDT= Table[ 0,20]
(*
 * Globale R Vector
 *)
R= Table[ 0,20]

(*
 * Globale DRDT Vector
 *)
DROLD= Table[ 0,20]
(*
 * Globale ROLD Vector
 *)
ROLD= Table[ 0,20]

(*
 * Globale RO  Vector  PolCar
 *)
RO= Table[ 0,20]

(*
 * Globale XO  Vector  PolCar
 *)
XO= Table[ 0,20]

(*
 * Globale XXX  Vector  PolCar
 *)
XXX= Table[ 0,20]

(*
 * Globale Y  Vector Was Y in Fortran but conflict with scalar Y. So use YY  
 *)
YY= Table[ 0,20]

F=W[[7]];

For [ ii=0, ii<256, ii++,
 
If[W[[ii]] ≠ 0.0 , 
             Print["W[",ii,"]=",W[[ii]]];
            ]
]

HOP=1.
(**********LET SUBROUTINES PRINTR AND RAYPLT KNOW THERE IS A NEW W ARRAY *)
NEWWP=True
NEWWR=True

(**********INITIALIZE PARAMETERS FOR INTEGRATION SUBROUTINE RKAM **)
NN=6
For[NR=7, NR ≤ 20, NR++,
 If [W[[50+NR]]≠ 0.,
       NN=NN+1
]
 

 ]

MODE=IntegrationType
STEP=STEP1
E1MAX=MAXERR
E1MIN=MAXERR/ERATIO
E2MAX=STPMAX
E2MIN=STPMIN
FACT=FACTR


XMTRH=W[[3]]
 (*  W[[24]]=PLAT
  W[[25]]=PLON *)

TLON=W[[5]]
TLAT=W[[4]]
Print ["XMTRH=",XMTRH," PLAT=",PLAT," PLON=",PLON]
Print [" TLAT=",TLAT," TLON=",TLON]
 

Print ["XMTRH=",XMTRH," LAT=",PLAT," LON=",PLON]


(***** STORE W Array   *******)
Export["W.wdx",W];







(*  *********DETERMINE TRANSMITTER LOCATION IN COMPUTATIONAL COORDINATES  *)
(*  *********SYSTEM (GEOMAGNETIC COORDINATES IF DIPOLE FIELD IS USED)   *)
R0=EARTHR+XMTRH
SP=Sin [PLAT]
CP=Cos [PLAT]
SDPH=Sin [TLON/DEGS-PLON]
CDPH=Cos [ (TLON/DEGS-PLON)]
SL=Sin [TLAT/DEGS]
CL=Cos [ TLAT/DEGS]
ALPHA=ArcTan[-SDPH*CP,-CDPH*CP*SL+SP*CL]
TH0=ArcCos [CDPH*CP*CL+SP*SL]
  PH0=ArcTan[SDPH*CL,CDPH*SP*CL-CP*SL] 
 

Print ["TH0=",TH0*DEGS, " PH0=",PH0*DEGS]

(*    Azimuth   *)
J=1
AZBEG=W[[11]]
AZSTEP= W[[13]] 
AZ1=AZBEG+(J-1)*AZSTEP
AZA=AZ1
AZ1=AZ1/DEGS;

GAMMA=Pi-AZ1+ALPHA
SGAMMA=Sin [GAMMA]
CGAMMA=Sin [PID2-GAMMA]

Print["Gamma SGAMMA CGAMMA ",GAMMA,SGAMMA,CGAMMA]

(*  Elevation  BETA  *)
ELSTEP=W[[17]]
ELBEG = W[[15]]

(* ELBEG=45  *)
III=1  (* For Loop for later *)
BETA=ELBEG-(III-1)*ELSTEP 

BETA= BETA/DEGS
EL=BETA*DEGS

Print["EL=",EL]

CBETA=Sin[PID2-BETA]

 W[[14]]=BETA;

R[[1]]=R0
R[[2]]=TH0
R[[3]]=PH0
R[[4]]=Sin[ BETA]
R[[5]]=CBETA*CGAMMA
R[[6]]=CBETA*SGAMMA

Print ["R=",R]

T=0.
 
 

(* Matrix of Rays ****)
RayMatrix ={}
AppendTo[RayMatrix,R];

kz=R[[4]]*Cos[R[[6]]]
     kxy=R[[4]]*Sin[R[[6]]]
    kx= kxy*Sin[R[[5]]]
    ky=kxy*Cos[R[[5]]]

kz=R[[4]] 
 
    kx=  R[[5]]
    ky= R[[6]]

a={0,0,0};
b={kx,ky,kz};
Print["b=",b]
Graphics3D[{Red,Arrow[{a,b}]},Axes->True,Boxed->True,AxesLabel->{x,y,z}]

a={0,0,0};
zz=R[[1]]*Cos[R[[2]]]
xx=R[[1]]*Cos[R[[2]]]*Cos[R[[3]]]
yy=R[[1]]*Cos[R[[2]]]*Cos[R[[3]]]
b={xx,yy,zz} 

Graphics3D[{Red,Arrow[{a,b}]},Axes->True,Boxed->True,AxesLabel->{x,y,z}]

T=0.;
RSTART=1.;
(**** SGN=l-(NEED FOR RAY TRACING IN COMPLEX SPACE-)                  ***)

(*****   ALLOW IONOSPHERIC MODEL SUBROUTINES TO READ AND PRINT DATA ***)
(* KRm KTH and KPH Undefined in FORTRAN at this point  *)

KR=1;
KTH=0.1;
KPH =0.1;
RINDEX[  ];
Print["FREQUENCY=",F,"MHZ ,AZIMUTH ANGLE OF TRANSMISION = ",AZA, " DEG"] 
Print ["ELEVATION ANGLE OF TRANSMISSION =",EL,"Degrees"]
Print["INITIAL N2=",N2];
(* N2 is Complex *)
If[Abs[N2] ≤ 0,
Electx[] ;
  FN=  (Sqrt [Abs [X]]*F)*Sign[X]  ;
Print[ "TRANSMITTER IN EVANESCENT REGION, TRANSMISSION IMPOSSBLE PLASMA FREGUENCY =",FN];
      Return[];
]
MU=Sqrt [N2/(R[[4]]^2+R[[5]]^2+R[[6]]^2)]
R[[4]]=R[[4]]*MU
R[[5]]=R[[5]]*MU
R[[6]]=R[[6]]*MU

For [kk=7, kk≤ NN, kk++,
       R[[kk]]=0.0;
]

TraceRT[]
(* THIS RAY CALCULATION TOOK  Calculate time in Mathematica to do calculations and print*)

(*Print RayMatrix *)

{0,0,0}

3.65*10^4