# Combining Gravity Turn and Orbit Models

I have a mathematical model for the motion of an orbiting spacecraft about Earth:

G = 6.672*10^-11; (*Gravitational Constant*)
M = 5.97219*10^24 ; (*Mass of Earth*)
R = 6.378 *10^6;(*Radius of Earth*)
r = R + 150000; (*Orbital radius*)
tmax = 5500;  (*Simulation time*)
v = Sqrt[(G M)/R]; (*Circular orbital velocity*)

orbit = NDSolve[{
x''[t] == -((G M x[t])/(x[t]^2 + y[t]^2)^(3/2)),
y''[t] == -((G M y[t])/(x[t]^2 + y[t]^2)^(3/2)),
x == 0, y == r, x' == v, y' == 0}, {x[t], y[t]}, {t, 0,
tmax} , MaxSteps -> 1000000, Method -> "StiffnessSwitching"]

ParametricPlot[Evaluate[{x[t], y[t]} /. orbit], {t, 0, tmax},
AxesLabel -> {x, y}, PlotStyle -> Automatic, PlotRange -> Full,
ImageSize -> Large] As well as a model for a spacecraft's surface launch using a gravity turn:

Remove["Global*"]
Unprotect[D]; (*Using symbol D for drag*)

G = 6.672*10^-11; (*Gravitational Constant*)
M = 5.97219*10^24 ; (*Mass of Earth*)
R = 6.378 *10^6;(*Radius of Earth*)
g0 = 9.81; (*Sea level gravitational acceleration*)
g = g0/(1 + h[t]/R)^2; (*Gravitational acceleration w.r.t. height*)

d = 5;  (*Diameter*)
A = (π d^2)/4;  (*Area*)
Subscript[C, D] = 0.5; (*Drag coefficient*)
Subscript[ρ, 0] = 1.225; (*Sea level air density*)
Subscript[h, 0] = 7500; (*Height scale*)
ρ = Subscript[ρ, 0]
Exp[-h[t]/Subscript[h, 0]]; (*Atmospheric air density*)
D = 1/2 ρ v[t]^2 A Subscript[C, D]; (*Drag*)
tburn = 260; (*Engine burn time*)
T = If[t <= tburn, 800000, 0];(*Thrust*)
m0 = 68000; (*Initial mass*)
mdot = If[t <= tburn, 244.1,  0];(*Engine mass flow rate*)
m = m0 - mdot*t; (*Mass of rocket w.r.t. time*)
tmax = 260; (*Simulation running time*)

traj = NDSolve[{
v'[t] == T/m - D/m - (g - v[t]^2/(R + h[t])) Sin[γ[t]],
γ'[t] == -(1/v[t]) (g - v[t]^2/(R + h[t])) Cos[γ[t]],
h'[t] == v[t] Sin[γ[t]],
x'[t] == v[t] Cos[γ[t]],
v == 1, γ == 90 Degree , x == 0, h == 0,
WhenEvent[h[t] == 1000, γ[t] -> 89.6 Degree]},
{v[t], γ[t], x[t], h[t]}, {t, 0, tmax} ]

ParametricPlot[{x[t], h[t]} /. traj, {t, 0, tmax}, AxesLabel -> {x, y}] What I'm trying to do, though, is somehow combine the two models so that I can simulate a gravity turn launch and subsequent orbit about Earth at the same time. As it stands, it seems like the gravity turn code unfortunately only considers a flat Earth surface. Does anyone know how I might go about combining the two and allowing for a curved surface in the gravity turn code?

EDIT: I've edited the above gravity turn code to take into account a curved Earth surface, but once the engine stops its burn something strange happens: The spacecraft carries on gaining speed as if its thrust is still going.

Remove["Global*"]
Unprotect[D]; (*Using symbol D for drag*)

G = 6.672*10^-11; (*Gravitational Constant*)
M = 5.97219*10^24 ; (*Mass of Earth*)
R = 6.378 *10^6;(*Radius of Earth*)

d = 5;  (*Diameter*)
A = (π d^2)/4;  (*Area*)
Subscript[C, D] = 0.5; (*Drag coefficient*)
Subscript[ρ, 0] = 1.225; (*Sea level air density*)
Subscript[y, 0] = 7500; (*Height scale*)
ρ = Subscript[ρ, 0]
Exp[-y[t]/Subscript[y, 0]]; (*Atmospheric air density*)
D = 1/2 ρ v[t]^2 A Subscript[C, D]; (*Drag*)
tburn = 260; (*Engine burn time*)
T = If[t <= tburn, 800000, 0];(*Thrust*)
m0 = 68000; (*Initial mass*)
mdot = If[t <= tburn, 244.1,  0];(*Engine mass flow rate*)
m[t] = m0 - mdot*t; (*Mass of rocket w.r.t. time*)
tmax = tburn; (*Simulation running time*)

traj = NDSolve[{
v'[t] ==
T/m[t] - D/
m[t] - ((G M)/(Sqrt[x[t]^2 + y[t]^2])^2 - v[t]^2/Sqrt[
x[t]^2 + y[t]^2]) Sin[γ[t]],
γ'[
t] == -(1/
v[t]) ((G M)/(Sqrt[x[t]^2 + y[t]^2])^2 - v[t]^2/Sqrt[
x[t]^2 + y[t]^2]) Cos[γ[t]],
y'[t] == v[t] Sin[γ[t]],
x'[t] == v[t] Cos[γ[t]],
v == 1, γ == 90 Degree , x == 0, y == R,
WhenEvent[y[t] == R + 1000, γ[t] -> 88.85 Degree]}, {v[
t], γ[t], x[t], y[t]}, {t, 0, tmax} ]

ParametricPlot[{x[t], y[t] - R} /. traj, {t, 0, tmax},
AxesLabel -> {x, y}]
Plot[{v[t]} /. traj, {t, 0, tmax}, AxesLabel -> {t, v}]
Plot[{γ[t]} /. traj, {t, 0, tmax}, AxesLabel -> {t, γ}]


This is the trajectory after the 260 second engine burn: And this is the trajectory after 5000 seconds: As can be seen, something is definitely awry.

• Hi ! Are you sure this is not a physics question ? I might be still sleepy, but I think you are having trouble with the underlying physics/mathematics. Nov 6 '14 at 11:05
• The Drexel University Mathematica Forum might make a better place to post this question. As I remember, the forum would entertain (even celebrate) any question either about (1.) how to do things in Mathematica or (2.) any problem that someone uses Mathematica to solve. This forum has a narrower focus, but I think it would benefit from expanding it to include questions like this one. Nov 6 '14 at 13:08
• Hi Sektor, although this is a physics based question, I put it in the Mathematica section because of the amount of code I posted. I'm hoping Mathematica users with a background in physcis/mathematics will be able to help. Nov 7 '14 at 14:14
• Hi Jagra, what is the focus of this Mathematica forum? Is it only used for language specific coding issues? Nov 7 '14 at 14:15
• Unprotect[D]; (*Using symbol D for drag*) <-- this is a really really bad idea ... it might also stop some people from reading the code because there's always the suspicion that it's causing something to go wrong (and it's so easy to fix). Nov 7 '14 at 14:39

This question really belongs on physics.stackexchange.com. You can post your propagation equations there using TeX. TeXForm[] will do the conversion for you if you don't know TeX.

For a 2-D case around a non-rotating spherical planet as you appear to be attempting, I use these propagation equations with NDSolve[]:

{v'[t]==-Sin[γ[t]] μ/r[t]^2-D/m,
γ'[t]==Cos[γ[t]] 1/r[t] (v[t]-μ/(v[t]r[t]))+L/(m v[t]),
r'[t]==Sin[γ[t]] v[t],
φ'[t]==Cos[γ[t]] 1/r[t] v[t],
v==v0,
γ==γ0,
r==r0,
φ==φ0}


$${dv\over dt}=-\frac{\mu}{r^2}\sin\gamma-\frac{D}{m}$$

$${d\gamma\over dt}={1\over r}{\left(v-\frac{\mu }{r v}\right)\cos\gamma}+\frac{L}{m v}$$

$${dr\over dt}=v\sin\gamma$$

$${d\phi\over dt}=\frac{v}{r}\cos\gamma$$

$\mu$ is $G M$. $r$ is the radius from the center of the body and $\phi$ is the central angle around the body. $x$ and $y$ can be readily computed from $r$ and $\phi$.

I use these for entries as opposed to launches ($D$ is drag and $L$ is lift), but you can replace $D$ and $L$ with whatever you like, including thrust along the velocity direction ($-D$) and orthogonal to the velocity direction ($\pm L$). Note that when using these I never set $D$ (or $L$) — instead I use /. D->....

If $D$ and $L$ are replaced with zeros, you get simple Kepler propagation, resulting in ellipses or hyperbolas about the body. So you can continue to use the same equations to propagate after the launch.

• Hi Mark, thanks very much. I had a go using these equations but had limited success. They worked great for simulating a rocket already in orbit, but I couldn't figure out the correct initial conditions for a successful gravity turn surface launch. Am I correct in assuming that with an initial radius = 6.378*10^6 and an initial central angle phi = 90 degrees, we'd have the rocket sitting on the "north pole". Furthermore, do you have any links showing the derivation of the above equations, as I'd be very interested in learning more and haven't ever seen equations of motion in this form before. Nov 7 '14 at 22:58
• The initial $\phi$ doesn't matter. You should set it to zero so you can easily see how much central angle is traversed. There is no "pole" here. All points on the surface of a non-rotating sphere are equivalent. Nov 7 '14 at 23:50
• The initial conditions would be $r$ at the surface and $v$ and $\phi$ at zero. Since $v$ is zero, the initial $\gamma$ doesn't matter. However if you will be applying thrust in the $\gamma$ direction, then you should set it to $\pi/2$ for straight up. Assuming you want to go straight up. Nov 7 '14 at 23:52
• Note that your thrust to weight ratio had better be more than one. Nov 7 '14 at 23:53
• To get the gravity turn, you will need to give it an initial kick angle at some point after launch to get it off vertical. Nov 8 '14 at 0:01