(* Copyright (C)2024 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.*)


(* Orbit Position as a Function of Time for Elliptic Orbits  *)
(* Orbital Mechanics                                               *)
(* Author:Sasan Ardalan                                            *)
(* ScientificWorks.org                                             *)
(* Date:September 10,2024                                          *)

(* Paamters from:                                                        *)
(* Fundamentals of Astrodynamics and Applications, 4th Edition, Vallado  *)
(* Earth Orbit                                                           *)

mu= 3.986004415*10^5;

e=0.6;
p=20410;
a=31890;
b=25512;
rp=12756;
ra=51024;
c=19134;
earthRadius=6371;

Print["Period is:"]
period=2*Pi*Sqrt[a*a*a/mu]/3600

Show[PolarPlot[p/(1+e*Cos[theta]), {theta, 1,2.5*Pi}], 
PolarPlot[earthRadius, {theta, 1,2.5*Pi}]]
(* Plot[Sqrt[mu*(2/(p/(1+e*Cos[theta])-1/a))],{theta, 1,2.5*Pi}] *)

(* Drawing tangent line on a curve MATHEMATICA
https://community.wolfram.com/groups/-/m/t/214757  *)

SinE[theta_,e_]:=Sqrt[1-e*e]Sin[theta]/(1+e*Cos[theta])
CosE[theta_,e_]:=(e+Cos[theta])/(1+e*Cos[theta])

SinTheta[Ev_,e_]:=(Sqrt[1-e*e]*Sin[Ev])/(1-e*Cos[Ev])

CosTheta[Ev_,e_]:=(Cos[Ev]-e)/(1-e*Cos[Ev])



 

mu= 3.986004415*10^5

e=0.6
p=20410
a=31890
b=25512
rp=12756
ra=51024
c=19134
earthRadius=6371
P=2*Pi*Sqrt[a*a*a/mu]
Phours=P/3600
n=2*Pi/P
 
 
satPosTime[th_,p_,e_,n_]:=Module[{t,t0,M0,EvRoot,T,M,Ev,thetaval},t=th*3600;t0=0; T=0; M0=n*(t0-T); M=n*(t-t0)+M0;
EvRoot=Ev/. FindRoot[M==Ev-e*Sin[Ev],{Ev,M0}]; 
sinTheta=SinTheta[EvRoot,e];
cosTheta=CosTheta[EvRoot,e];
 
thetaval=ArcTan[sinTheta/cosTheta];
thetaval=If[thetaval <0, thetaval+Pi,thetaval];
  Return[thetaval]  ]
 
 


Print["Sat Positione Time:"]
satPosTime[0.1,p,e,n]*180/Pi
 
satpos[th_,p_,e_,n_]:=p/(1+e*Cos[satPosTime[th,p,e,n]]);

Print["Sat Radius Time:"]
th1=7
satPosTime[th1,p,e,n]*180/Pi
satpos[th1,p,e,n]

track={,}
path={,}
th1=0
delta=0.1;
For[i=0,i<75,i++, th1=th1+delta;
        theta=satPosTime[th1,p,e,n];
         r=satpos[th1,p,e,n];
         velocity=Sqrt[mu*(2/r-1/a)];
         (* track=Append[track,{theta,velocity}]; *)
         (*track=Append[track,{th1,velocity}];*)

         track=Append[track,{th1,  theta  }];

         path=Append[path,{theta,r}];

        (*Print[r,":",theta*180/Pi,":",velocity];*)
]
(*Print ["The Track"]
Print[track]*)
ListPlot[track,AxesLabel->{Time,theta}]
(* Print ["The Path"]
Print[path] *)
ListPolarPlot[path]