Combining Gravity Turn and Orbit Models

Question

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] == 0, y[0] == r, x'[0] == v, y'[0] == 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[0] == 1, γ[0] == 90 Degree , x[0] == 0, h[0] == 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[0] == 1, γ[0] == 90 Degree , x[0] == 0, y[0] == 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.

Answer 1

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[0]==v0,
γ[0]==γ0,
r[0]==r0,
φ[0]==φ0}

Or more readably:

$${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.