Solution

First, load the necessary constants.

In[31]:=

<<Miscellaneous`PhysicalConstants`

R = EarthRadius/Meter//N ;

g = AccelerationDueToGravity/(Meter/Second^2)//N ;

Then, the parameters

In[34]:=

m0 = 1000 ; f0 = 28000 ; c = 1 ; p = 600 ;

Define the time-dependent parameters, i.e. mass and thrust force.

In[35]:=

m[t_] := If[c t<p, m0 - c t, m0 - p] ;

fm[t_] := If[c t<p, f0, 0] ;

Define the position and velocity vectors. All forces will be presented in terms of these.

In[37]:=

Clear[r, t, x, y, f, v] ;

r[t_] = {x[t], y[t]} ;

v[t_] = {x '[t], y '[t]} ;

f[t_] = fm[t] v[t]/Norm[v[t]] ;

grav[t_] = g R^2r[t]/Norm[r[t]]^3 ;

A differential equation for both directions.

In[42]:=

dex = m[t] x''[t] Part[f[t], 1] - m[t] Part[grav[t], 1] ;

dey = m[t] y''[t] Part[f[t], 2] - m[t] Part[grav[t], 2] ;

A fancy way to determine the initial conditions.

In[44]:=

latitude = 90 ;

x0 = R Cos[latitude * Pi/180] ;

y0 = R Sin[latitude * Pi/180] ;

vx0 = 0.001 ;

vy0 = 0.1 ;

(*vyinit = vinit Sin[latitude * Pi/180] ; *)

Now, commence with the solution.

In[49]:=

T = 2000 ;

Sol = Part[NDSolve[{dex, dey, x[0] x0, x '[0] vx0, y[0] y0, y '[0] vy0}, {x[t], y[t]}, {t, 0, T}], 1] ;

flight = ParametricPlot[{x[t]/.Sol, y[t]/.Sol}, {t, 0, T}, DisplayFunctionIdentity] ;

Draw a picture of the Earth to get the scale of the trajectory visible.

In[52]:=

Earth = Graphics[Circle[{0, 0}, R]] ;

In[53]:=

Show[flight, Earth, DisplayFunction$DisplayFunction, PlotRangeAll, TicksNone] ;

Plot[Norm[r[t]/.Sol], {t, 0, T}, TicksNone]

[Graphics:../HTMLFiles/ex08_solutions_109.gif]

[Graphics:../HTMLFiles/ex08_solutions_110.gif]

Out[54]=

⁃Graphics⁃

Created by Mathematica  (April 10, 2007) Valid XHTML 1.1!