(*  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 Ray Points in 3D with Origin at Earth Center                        *)

(*  Uses Plot Range  to Show Ray Paths                                         *)
(***************************************************************************)
(***********************   RAY_PATH Ray Vector 3D Plot      *****************)
(***************************************************************************)


RayMatrix =Import["rayMatrix.wdx"];
 
lengthMatrix=Length[RayMatrix]
Print ["Length Matrix=",lengthMatrix];

WPost=Import["W.wdx"];


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  *)



(* Store in column format in text file *)

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];

a={0,0,0};
arrowArray={}
AppendTo[arrowArray,{Green,Ball[a,6370]}];
Graphics3D[arrowArray,Axes->True,Boxed->True,AxesLabel->{x,y,z} ];
Print["lengthMatrix=",lengthMatrix];
For[ii=3,ii<lengthMatrix, ii++,

   

Ray=RayMatrix[[ii]] ;
spherical={Re[Ray[[1]]],Re[Ray[[2]]],Re[Ray[[3]]]};

(* Print["Spherical+",spherical]; *)
(* FromSphericalCoordinates *)
  b=FromSphericalCoordinates[spherical];  

(* bb={spherical[[1]]*Cos[spherical[[3]]]*Sin[spherical[[2]]],spherical[[1]]*Sin[spherical[[3]]]*Sin[spherical[[2]]],spherical[[1]]*Cos[spherical[[2]]]}; *)
R=Sqrt[b[[1]]^2+b[[2]]^2+b[[3]]^2];
  Print["ii=",ii," b=",b, "  magnitude b=",R , " (R - EARTHR)=  ",R-EARTHR ];  
 
 AppendTo[arrowArray,{Red,Point[{a,b} ]}];  
(* AppendTo[arrowArray,{Red,Arrow[{a,b} ]}]; *)


];

(*Print ["arrowArray=",arrowArray];*) 

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

Graphics3D[arrowArray,Axes->True,Boxed->True,AxesLabel->{x,y,z},PlotRange->{{3000,3610},{-3340,-3302},{4000,4700}}];  (* 45 Degrees *)

(* Graphics3D[arrowArray,Axes->True,Boxed->True,AxesLabel->{x,y,z},PlotRange->{{-3020,-2900},{3100,3300},{4700,4900}} ];  90 Degrees *)

%566

%534