Satellite orbits

updated version 2012 for Mathematica 8
code by jph, www.hakenberg.de

Instructions: Evaluation -> Evaluate notebook

General settings and macros:

Off[General::spell]; 
Off[General::spell1];
𝓈:=FullSimplify[#1,assume]&;
Τ★[v_,n_]:=Array[(Subscript[v, ##1])&,n];
X_∘ϕ_:=X/.Τ★[x,Length[ϕ]]→ϕ;
Μd[X_List,r_Integer,n_Integer]:=Transpose[Array[D[X,(Subscript[x, #1])]&,{n}],Flatten[{r+1,Range[r]}]];
Μd[X_List,r_Integer]:=Μd[X,r,Length[X]];
Μd[X_,n_Integer]:=Array[D[X,(Subscript[x, #1])]&,n];
x_Symbol→y_List:=Array[(Subscript[x, #1])->y[[#1]]&,Length[y]];
x_List→y_List:=Array[x[[#1]]->y[[#1]]&,Min[Length[x],Length[y]]];
x_List→0:=Array[x[[#1]]->0&,Length[x]];

Unfortunately, we are using a global variables h and T.
And why did we manually implement a Runge-Kutta scheme, instead of using Mathematica’s own function to integrate?

Clear[h,T];
SRungeKutta[X_]:=Module[{t=X[[1]],x=X[[2]],k1,k2,k3,k4},
k1=Φ[t,x];
k2=Φ[t+.5 h,x+.5 h k1];
k3=Φ[t+.5 h,x+.5 h k2];
k4=Φ[t+h,x+h k3];
{t+h,x+h {(1/6),(1/3),(1/3),(1/6)}.{k1,k2,k3,k4}}];

Our computation is set to integrate 1000 ‘steps’. But this is set
The “init” function converts a TLE into Euler coordinates - best look at the formula inside to check the units.

SIntegrate[X_,T_]:=Module[{res},h=T/Ν;
res=NestList[SRungeKutta,{0,X},Ν];
Interpolation[Transpose[{res[[All,1]],#1}]][t]&/@Transpose[res[[All,2,{1,2,3}]]]];
STarget[X_,T_]:=Module[{res},h=T/Ν;
Nest[SRungeKutta,{0,X},Ν]];
SInit[μ_,a_,ϵ_,ω_,i_,Ω_,n_]:=Module[{
RΩ=MatrixExp[(({{0, -(Subscript[x, 3]), (Subscript[x, 2])}, {(Subscript[x, 3]), 0, -(Subscript[x, 1])}, {-(Subscript[x, 2]), (Subscript[x, 1]), 0}}))∘{0,0,Ω(π/180)}],
Ri=MatrixExp[(({{0, -(Subscript[x, 3]), (Subscript[x, 2])}, {(Subscript[x, 3]), 0, -(Subscript[x, 1])}, {-(Subscript[x, 2]), (Subscript[x, 1]), 0}}))∘{i(π/180),0,0}],
Rω=MatrixExp[(({{0, -(Subscript[x, 3]), (Subscript[x, 2])}, {(Subscript[x, 3]), 0, -(Subscript[x, 1])}, {-(Subscript[x, 2]), (Subscript[x, 1]), 0}}))∘{0,0,ω(π/180)}],
rp=a (1-ϵ),mat},mat=RΩ.Ri.Rω;
{Flatten[{mat.{rp,0,0},mat.{0,(Sqrt[((μ (1+ϵ))/(a (1-ϵ)))]),0}}],T=n 2 π (Sqrt[((a^3)/μ)])}]

Analytical transformation from spherical to Euler coordinates.
Also, their analytical differentials.

assume={0<(Subscript[x, 1])<ε,0<(Subscript[x, 2])<ε,0<(Subscript[x, 3])<ε};
ϕ={(Sqrt[(Subscript[x, 1]^2)+(Subscript[x, 2]^2)+(Subscript[x, 3]^2)]),ArcTan[(Subscript[x, 1]),(Subscript[x, 2])],ArcTan[(Sqrt[(Subscript[x, 1]^2)+(Subscript[x, 2]^2)]),(Subscript[x, 3])]};
φ=(Subscript[x, 1]){Cos[(Subscript[x, 2])] Cos[(Subscript[x, 3])],Sin[(Subscript[x, 2])] Cos[(Subscript[x, 3])],Sin[(Subscript[x, 3])]};
dϕ=Μd[ϕ,1]//𝓈;dϕ//MatrixForm
dφ=Μd[φ,1]//𝓈;dφ//MatrixForm

It seems we are using km instead of SI.

(Subscript[r, E])=6378.14;
(Subscript[μ, E])=398900.;
(Subscript[μ, S])=1.327 (10^11);
(Subscript[ι, E])=23.452 π/180;
(Subscript[σ, E])=86164.11605509756;
(Subscript[d, E])=86400;

υJ={0,1082.6*^-6,-2.53*^-6,-1.62*^-6,-0.23*^-6,0.54*^-6};
υC={{},{0,1.57*^-6},{-2.19*^-6,0.31*^-6,-0.1*^-6},{0,0,0,0},{0,0,0,0,0},{0,0,0,0,0,0}};
υS={{},{0,-0.9*^-6},-{0.27*^-6,0.21*^-6,0.2*^-6},{0,0,0,0},{0,0,0,0,0},{0,0,0,0,0,0}};

Υ[ρ_,θ_,λ_]:=-((Subscript[μ, E])/ρ)+((Subscript[μ, E])/ρ)Sum[(((((Subscript[r, E])/ρ)))^n)υJ[[n]] LegendreP[n,0,-Sin[θ]]+ 
Sum[(((((Subscript[r, E])/ρ)))^n)(υC[[n,m]] Cos[m λ]+υS[[n,m]] Sin[m λ]) LegendreP[n,m,-Sin[θ]],{m,n}],
{n,2,6}];
U=Υ[(Subscript[x, 1]),(Subscript[x, 3]),(Subscript[x, 2])]//Simplify;
dU=Simplify[Μd[U,3]∘ϕ].dϕ;

In Mathematica 5.2, it was not necessary to use Re[] at this point.
SATELLITE denotes the function that Runge-Kutta uses for integration.

hU=Re[dU∘{x,y,z}];
cU=Compile[{x,y,z},Evaluate[hU]];
SATELLITE[t_,X_]:=Join[X[[{4,5,6}]],-cU@@X[[{1,2,3}]]];

ℙ[𝔱_]:=With[{X=ϕ∘(pos/.t->𝔱)},{Mod[X[[2]]-((2 π 𝔱)/(Subscript[σ, E])),2π,-π],X[[3]]}];
drwfnc:=Line/@Module[{c1,pnto,pntn,lst={},cur={},h=T/100,t=0.,c,stp},
stp={h,h/4,h/16,h/64};
While[t<T-h,
pnto=ℙ[t];
pntn=ℙ[t+#1]&/@stp;
c=Count[Norm[#1-pnto]-.1&/@pntn,_?Positive];
If[c<Length[stp],
AppendTo[cur,pntn[[c+1]]];
t+=stp[[c+1]];,
If[Sign[pnto[[1]] pntn[[1,1]]]==-1,
AppendTo[cur,pntn[[-1]]+pnto[[1]] {2π,0}];
AppendTo[lst,cur];
cur={pnto-pnto[[1]] {2π,0}};,
AppendTo[lst,cur];cur={};
];
t+=stp[[-1]];
];
];
AppendTo[lst,cur];lst
];

An example for a satellite with 
• altitude 23222 km
• eccentricity 0.3
• tracking over 75% of the analytical orital period 2 π (Sqrt[((a^3)/μ)])
• 1000 integration steps

Φ=SATELLITE;
Ν=1000;
Shome=SInit[(Subscript[μ, E]),(Subscript[r, E])+23222,0.3,10,55,13,0.75];
pos=SIntegrate@@Shome;
pP=ParametricPlot3D[pos,{t,0,T},DisplayFunction->Identity,PlotPoints->200];
pS=ParametricPlot3D[Evaluate[φ/.(Subscript[x, 1])->(Subscript[r, E])],{(Subscript[x, 2]),-π,π},{(Subscript[x, 3]),-π/2,π/2},DisplayFunction->Identity];
Show[pP,pS,DisplayFunction->$DisplayFunction,PlotRange->All]

We combine the ground track of the satellite with Mathematica’s (offline) WorldData.
I am not sure if the scale factor 1000π is correct!

<<WorldPlot`
Nation[L_]:={#1[[2]],#1[[1]]}/(1000π)&/@WorldData[L][[1]];
The=Nation/@World;
map=Graphics[Flatten[{
RGBColor[0.75,0.75,0.75],
Line/@The}]];

Show[{map,Graphics[
Flatten[{
RGBColor[0.5,0,0.5],
drwfnc,
Table[{PointSize[300/(ϕ∘pos)[[1]]],RGBColor[0,1-.7 t/T,0],Point[ℙ[t]]},{t,0,T-1,60 60}]
}]]},PlotRange->{{-π,π},{-π/2,π/2}}]