# How to Solve these equations by Mathematica and plot the results' curve?

I'm working on an Apollo re-entry program(from 160 Km to 60 Km) and would like to solve the equations of the trajectory. I'm using these equations of motion to calculate the trajectory of the spacecraft.

V̇=-(D/m)-g Sin[θ]

OverDot[φ]=-(V Cos[θ]/(r+h))

OverDOt[θ]=(L/m V)-(g Cos[θ]/V)

ḣ=V Sin[θ]

Drag and lift are calculated using Newtonian Theory:

cd=2 Sin[α]^3

cl=2 Sin[α]^2 Cos[α]

D=(1/2) cd ρ V^2 s

L=(1/2) cl ρ V^2 s

Density and gravity:

g=g0*(r/r+h)

ρ=ρ0*Exp[h0-h/H]

Constants:

s=Quantity[12.0687,"Meters"^2];(*surface area of the capsule*)

r=Quantity[6371,"Kilometers"];(*radius of the Earth*)

g0=Quantity[9.80612,"Meters"*"Seconds"^-2];(*gravity at 160 km*)

ρ0=Quantity[1.20654,"Kilograms"*"Meters"^-3];(*density at 160 km*)

m=Quantity[5500,"Kilograms"];(*mass of capsule*)

h0=Quantity[60,"Kilometers"];(*end of sensible atmosphere*)

H=Quantity[7,"Kilometers"];(*scaled height*)

Inputs:

V=Quantity[7162.8,"Meters"*"Seconds"^-1]; (*velocity of capsule*)

φ=10.556°; (*angular range*)

θ=-2.0°; (*flight path angle*)

α=53°; (*attack angle*)

h=Quantity[160,"Kilometers"]; (*altitude*)

I am just beginning on working on Mathematica and would like to ask how can I solve these equations and get as a result a curve of altitude(km) vs time(s)

• The air density profile is really given by the equation ρ=ρ0*Exp[h0-h/H]? Cause this expression do not reproduce the expected behavior of the atmosphere. Is the attack angle considered constant over the entire reentrance? Jun 1, 2017 at 16:14

You can solve the above system of coupled differential equations using NDSolve Mathematica function numerically.

\begin{align} \frac{\mathrm dV}{\mathrm dt} &= -\frac{D}{m}-g\sin\theta\\ \frac{\mathrm d\theta}{\mathrm dt} &= \frac{1}{V}\left(\frac{L}{m} -g\cos\theta\right)\\ \frac{\mathrm d\varphi}{\mathrm dt} &= -\frac{V\cos\theta}{r+h}\\ \frac{\mathrm dh}{\mathrm dt} &= V\sin\theta \end{align}

First one needs to assign correctly the equations, initial values and its variables. The expressions for g and $$\rho$$ are incorrect and have to be modified. A little care must to be taken in respect to different units used (km,m). The drag force D symbol cannot be used because MMA will recognize this as the derivative function. Here I changed it to Fd. Trigonometric functions entries must be in radians instead degrees and need to be changed with N[angle Degrees]. So we have

s = 12.0687;                       (*surface area of the capsule m^2*)
r = 6371000;                       (*radius of the Earth m *)
g0 = 9.80612;                      (*gravity at surface m/s^2 *)
\[Rho]0 = 1.20654;                 (*density at surface kg/m^3 *)
m = 5500;                          (*mass of capsule kg *)
h0 = 60000;                        (*end of sensible atmosphere m *)
H = 7000;                          (*scaled height m *)

V0 = 7162.8;                       (*velocity of capsule  m/s*)
h0 = 160000;                       (*altitude m *)
\[CurlyPhi]0 = N[10.556 Degree];   (*angular range*)
\[Theta]0 = N[-2.0 Degree];        (*flight path angle*)
\[Alpha] =  N[53 Degree];          (*attack angle*)

g[z_] := g0*(r^2/(r + z)^2);
\[Rho][z_] := \[Rho]0*Exp[-z/H];

cd = 2 (Sin[\[Alpha]])^3;
cl = 2 (Sin[\[Alpha]])^2 Cos[\[Alpha]];

Fd[v_, z_, t_] := (1/2) cd \[Rho][z] v^2 s;
L[v_, z_, t_] := (1/2) cl \[Rho][z] v^2 s;

sol = NDSolve[{
V'[t] == -(Fd[V[t], h[t], t]/m) - g[h[t]] Sin[\[Theta][t]],
\[Theta]'[t] == L[V[t], h[t], t]/(m V[t]) - g[h[t]] Cos[\[Theta][t]]/V[t],
\[CurlyPhi]'[t] == -(V[t] Cos[\[Theta][t]]/(r + h[t])),
h'[t] == V[t] Sin[\[Theta][t]],
h[0] == h0,
V[0] == V0,
\[Theta][0] == \[Theta]0,
\[CurlyPhi][0] == \[CurlyPhi]0
}, {V, \[Theta], \[CurlyPhi], h}, {t, 0, 800}];

pl1 = Plot[Evaluate[V[t] /. sol], {t, 0, 800}, PlotRange -> All,
FrameLabel -> {StyleForm["t (s)", FontSize -> 13],
StyleForm["v (km/s)", FontSize -> 13]},
Frame -> True, PlotStyle -> {Directive[Red, Thick]},
ImageSize -> 300];
pl2 = Plot[Evaluate[(180*\[Theta][t]/\[Pi]) /. sol], {t, 0, 800},
PlotRange -> All,
FrameLabel -> {StyleForm["t (s)", FontSize -> 13],
StyleForm["\[Theta] (degrees)", FontSize -> 13]},
Frame -> True, PlotStyle -> {Directive[Green, Thick]},
ImageSize -> 300];
pl3 = Plot[Evaluate[(180*\[CurlyPhi][t]/\[Pi]) /. sol], {t, 0, 800},
PlotRange -> All,
FrameLabel -> {StyleForm["t (s)", FontSize -> 13],
StyleForm["\[CurlyPhi] (degrees)", FontSize -> 13]},
Frame -> True, PlotStyle -> {Directive[Orange, Thick]},
ImageSize -> 300];
pl4 = Plot[Evaluate[(h[t]/1000) /. sol], {t, 0, 800},
PlotRange -> All,
FrameLabel -> {StyleForm["t (s)", FontSize -> 13],
StyleForm["h (km)", FontSize -> 13]},
Frame -> True, PlotStyle -> {Directive[Blue, Thick]},
ImageSize -> 300];
GraphicsGrid[{{pl1, pl2}, {pl3, pl4}}]

ApoloPath =
Table[{\[CurlyPhi][t], ((r + h[t])/1000)}, {t, 0, 800}] /. sol;
pl5 = ListPolarPlot[ApoloPath, PlotStyle -> Red,
PolarAxes -> Automatic,
PolarGridLines -> {{0, Pi/2, Pi, 3 Pi/2}, {0.2*r/1000, 0.4*r/1000,
0.6*r/1000, 0.8*r/1000}},
PolarTicks -> {"Degrees", Automatic}];
pl6 = Graphics[{Hue[0.6, 0.3, 0.7, .7],
Disk[{0, 0}, r/1000, {-0.3, 0.3}]}];
Show[pl5, pl6,
PlotRange -> {{0, r/1000}, {-7000 Sin[0.3], 7000 Sin[0.3]}}]
ApoloPathTime =

Table[{r (\[CurlyPhi][t] + 0.15336)/
1852., (h[t]/(1000*0.3048))}, {t, 0, 800}] /. sol;

TimeTicks =
Table[{r (\[CurlyPhi][t] + 0.15336)/
1852., (h[t]/(1000*0.3048))}, {t, 0, 800, 30}] /. sol;

pl7 = ListLinePlot[ApoloPathTime,
ScalingFunctions -> {"Reverse", Identity}, PlotRange -> All,
FrameLabel -> {StyleForm["Range (Nautical Miles)", FontSize -> 13],
StyleForm["h (1000 ft)", FontSize -> 13]}, Frame -> True,
PlotStyle -> {Directive[Black, Thick]}, ImageSize -> 300];

pl8 = ListPlot[TimeTicks, ScalingFunctions -> {"Reverse", Identity},
PlotMarkers -> {"|"}, PlotRange -> All,
FrameLabel -> {StyleForm["Range (Nautical Miles)", FontSize -> 13],
StyleForm["h (1000 ft)", FontSize -> 13]}, Frame -> True,
PlotStyle -> {Directive[Red, Thick]}, ImageSize -> 300];

Show[pl7, pl8]


In this manner the output functions will be

The above behavior shows us that the landing must occur in a time around 800 seconds. Also, one can plot the the trajectory of the Apollo spacecraft parametrically ($$\varphi(t)$$,$$r+h(t)$$) in polar coordinates adjusting correctly the units to meters and degrees.

These results can be compared with the next plot