eqofmov.html

Stanislav Barton

Mendel University of Agriculture in Brno Karel Englis University

Zemedelska 1 Sujanovo namesti 1

613 00 Brno & 602 00 Brno

Czech Republic Czech Republic

barton@mendelu.cz

Equations of the movement - modeling of the physical process

The equations of movement represent usually the system of common differential equations of second order. This system could be formally written in form of:
, where = positional vector of the moving object, m = its mass.

The left side of equation contains the components of vector of force described by mass times acceleration in the corresponding coordinate. The right side contains function of coordinates and their first time derivatives. It is possible to find an analytic solution in the simple case, but mostly there is a very complicated system of non-linear equations. This system is soluble numerically only. Let's show the application of movement equations at several interesting cases.

1. Movement in the Atmosphere

1.1. Equations of the movement

> restart;

> with(plots):

Warning, the name changecoords has been redefined

> Fa:=-1/2*Cx*S*rho*V^2; Newton's term of the aerodynamics resistance - absolute value of the force

> Fax:=Fa*cos(alpha); Fay:=Fa*sin(alpha); Components of the vector of the aerodynamics force

> Fax:=subs(cos(alpha)=Vx/V,Fax); Fay:=subs(sin(alpha)=Vy/V,Fay); Substitution of the goniometric functions

> Vx:=diff(x(t),t); Vy:=diff(y(t),t); V:=sqrt(Vx^2+Vy^2); = [diff(x(t),t), diff(y(t),t)]

> Ax:=diff(x(t),t,t)=Fax/m; Ay:=diff(y(t),t,t)=Fay/m-g; The final equations of movement, describing aerodynamics deceleration in x and y direction.

> save Ax, Ay, `Equations.sav`; Equations will be used more times, save them into file

1.2. Free fall with the constant resistance

> Ffall:=simplify(subs(Vx=0,Ay)); In case of the free fall we shall orientate system so that y axis is antiparallel with the gravitational acceleration. in this case.

> Ffall:=subs(csgn(diff(y(t),t))=-1,Ffall); If the square root in the equation Ay can be simplified into or -diff(y(t),t) . The physically relevant is the negative value.

> Ini:=y(0)=H,D(y)(0)=0; The initial conditions

> Y:=dsolve({Ffall,Ini},y(t)); Analytical solution of the free fall trajectory as the function of the time.

> Y:=simplify(rhs(Y),symbolic); Simplification

> Su:=2^(1/2)/m^(1/2)*g^(1/2)*Cx^(1/2)*S^(1/2)*rho^(1/2)*t=T,exp(T)=ET,exp(2*T)=ET^2,m=M^2; These substitutions leads to other simplifications.

> Y1:=simplify(subs(Su,Y),symbolic); Simplification

> Bs:=ET=exp(1/m^(1/2)*2^(1/2)*g^(1/2)*Cx^(1/2)*S^(1/2)*rho^(1/2)*t),M=sqrt(m); Back substitution

> Y:=unapply(subs(Bs,Y1),Cx,t); y:=unapply(limit(Y(Cx,t),Cx=0),t); Height as the function of the Cx and t . Height of the free fall in the vacuum is the limit case for Cx = 0 .

> V:=unapply(simplify(diff(Y(Cx,t),t)),Cx,t); v:=unapply(limit(V(Cx,t),Cx=0),t); Velocity as the function of the Cx and t . Velocity of the free fall in the vacuum is the limit case for Cx = 0 .

> A:=unapply(simplify(diff(V(Cx,t),t)),Cx,t); a:=unapply(limit(A(Cx,t),Cx=0),t); Acceleration as the function of the Cx and t . Acceleration of the free fall in the vacuum is the limit case for Cx = 0 .

> Su:=g=9.81,S=evalf(Pi*0.25^2),rho=1.2,H=1200.0,m=100.0; Numerical values used for the numerical computations

> Ys:=unapply(subs(Su,Y(Cx,t)),Cx,t); ys:=unapply(subs(Su,y(t)),t);
Vs:=unapply(subs(Su,V(Cx,t)),Cx,t); vs:=unapply(subs(Su,v(t)),t);
As:=unapply(subs(Su,A(Cx,t)),Cx,t); as:=unapply(subs(Su,a(t)),t); Substituted functions

> CX:=1/100*[$1..100]: Vector of the aerodynamics coefficients

> Tcx:=map(u->fsolve(subs(Cx=u,Ys(Cx,t)),t=14..20),CX): Duration of the free fall for different Cx, Cx<>0!

> Vcx:=zip((u,v)->evalf(Vs(u,v)),CX,Tcx): Landing velocity as function of Cx

> Acx:=zip((u,v)->evalf(As(u,v)),CX,Tcx): Acceleration as function of Cx

> Tcx:=[fsolve(ys(t),t=14..20),Tcx[]]: Duration of the free fall for Cx = 0 is added as the first value into Tcx

> Vcx:=[vs(Tcx[1]),Vcx[]]: Landing velocity for Cx = 0 is added

> Acx:=[as(Tcx[1]),Acx[]]: Landing acceleration for Cx = 0 is added

> CX:=[0,CX[]]: Value Cx = 0 is added

>    display({plot(ys(t),t=0..20,color=red,thickness=3),
         plot({seq(Ys(Cx,t),Cx=0.05*[$0..20])},t=0..20,color=blue)},
         title="Y(t,Cx)",labels=["t","Y"],view=[0..20,0..1200]); Plot of the height as the function of Cx and t . Cx=0 .

>    display({plot({seq(-Vs(Cx,t),Cx=0.05*[$1..20])},t=0..20,color=blue),
         plot(zip((u,v)->[u,-v],Tcx,Vcx),color=black,thickness=3),
         plot(-vs(t),t=0..20,color=red,thickness=3)},
        labels=["t","V"],title="V(t,Cx)"); Velocity as the function of the Cx and time. Cx=0 . Black curve indicates landing velocity and duration of the fall from the height H=1200m.

>    display({plot({seq(-As(Cx,t),Cx=0.05*[$1..20])},t=0..20,color=blue),
         plot(zip((u,v)->[u,-v],Tcx,Acx),color=black,thickness=3),
         plot(-as(t),t=0..20,color=red,thickness=3)},
         labels=["t","A"],title="A(t,Cx)"); Acceleration as the function of the Cx and time. Cx=0 . Black curve indicates landing acceleration and duration of the fall from the height H=1200m..

>    display({plot({seq([Ys(Cx,t),-Vs(Cx,t),t=0..20],Cx=0.05*[$1..20])},color=blue),
         plot([ys(t),-vs(t),t=0..20],color=red,thickness=3)},
         view=[0..1200,0..160],labels=["h","V"],title="V(h,Cx)"); Velocity as the function of the Cx and height. Cx=0 .

>    display({plot({seq([Ys(Cx,t),-As(Cx,t),t=0..20],Cx=0.05*[$1..20])},color=blue),
         plot([ys(t),-as(t),t=0..20],color=red,thickness=3)},
         view=[0..1200,0..10],labels=["h","A"],title="A(h,Cx)"); Acceleration as the function of the Cx and height. Cx=0 .

>    display({plot({seq([-Vs(Cx,t),-As(Cx,t),t=0..20],Cx=0.05*[$1..20])},color=blue),
         plot([-vs(t),-as(t),t=0..20],color=red,thickness=3)},
         plot(zip((u,v)->[-u,-v],Vcx,Acx),color=black,thickness=3),
         view=[0..160,0..10],labels=["V","A"],title="A(V,Cx)"); Acceleration as the function of the Cx and velocity. Cx=0 . Black curve indicates landing acceleration and velocity of the fall from the height H=1200m..

>    display({plot3d(-As(Cx,t),t=0..20,Cx=0..1,grid=[40,40]),
         spacecurve([seq([Tcx[j],CX[j],-Acx[j]],j=1..nops(CX))],color=black,thickness=8)},
         axes=boxed,title="A(t,Cx)",orientation=[40,50],labels=["t","Cx","A"]); Plot of the acceleration as a function of the Cx and t . Black curve indicates landing acceleration from H = 1200m ,

>    display({plot3d([Ys(Cx,t),-Vs(Cx,t),-As(Cx,t)],t=0..20,Cx=0..1,grid=[40,40],style=patchnogrid,shading=XY),
                 spacecurve({seq([Ys(Cx,t),-Vs(Cx,t),-As(Cx,t)],Cx=0.1*[$1..10])},t=0..20,color=black),
                 spacecurve({seq([Ys(Cx,t),-Vs(Cx,t),-As(Cx,t)],t=2*[$0..10])},Cx=0..1,color=red)},
                 axes=boxed,orientation=[-150,70],
                 view=[0..1200,0..160,0..10],labels=["h(Cx,t)","V(Cx,t)","A(Cx,t)"]); All variables h(Cx,t), V(Cx,t) and A(Cx,t), describing free fall can be determined using this 3d plot

Constan t time t=2s , Constant resistance Cx=0.1

1.3. Dynamics of the Skydiver

Let us assume skydiver falling from the initial height 3000 m. He is falling approximately free for =30 seconds, later he opens his parachute. Let us assume coefficient of his aerodynamics resistance Cx0=0.2 , cross section of his body S0=0.2 , if he opens parachute changes into Cxp=1.5 and cross section of the parachute will be SP=45 .

Weight of the skydiver is m=100 kg and g=9.81 kg/ is the gravitational acceleration.

> tau:=30.: Cx0:=0.2: CxP:=0.6: S0:=0.2: SP:=45: C:=1.5: m:=100.: g:=9.81: y:='y':

> Cx:=Cx0+(CxP-Cx0)*((t-tau)/sqrt((t-tau)^2+C^2)/2+1/2);
S:=S0+(SP-S0)*((t-tau)/sqrt((t-tau)^2+C^2)/2+1/2); Function describing increasing of Cx and S as a function of time. The variable C determines slope of the function = describes velocity of opening of the parachute - increasing of the aerodynamic coefficient and the cross section surface.

> plot(Cx*S,t=0..40,labels=["t [s]","Cx*S"],title="Opening of parachute");

> Skydiver:=subs(m=100.,g=9.81,rho=1.2,Ffall); Equation of skydiver motion - nonlinear differential equation of the 2nd order. Can be solved only numerically.

> Ini:=y(0)=3000.,D(y)(0)=0; Initial conditions

> Y:=dsolve({Skydiver,Ini},y(t),numeric); Numerical solution

> plots[odeplot](Y,[t,y(t)],0..80,labels=["t [s]","h [m]"],title="H(t)"); Graph of the height as the function of time. Landing time can be determined using Newton's numerical procedure. This procedure computes corrections dT for T solving Y(T)=0, where dT = -Y(T)/diff(Y(T),T) . As we can see, it is the simplest to use Maple's numerical procedure Y to compute values of Y(T) and .

> T:=70.0; dT:=1; Approximate value of T leading to Y(T)=0 . dT = initial correction of T .

> Y(T); Values returned by the procedure Y .

> while abs(dT)>0.000001 do; dT:=subs(Y(T),-y(t)/diff(y(t),t)); T:=T+dT; od; Loop computing landing time of the skydiver, is running till dT is greater than or equal as desired accuracy.

>    display([odeplot(Y,[t,y(t)],0..T,numpoints=1000,color=black,legend="Height [m]"),
         odeplot(Y,[t,-20*diff(y(t),t)],0..T,numpoints=1000,color=red,legend="Velocity [*20m/s]"),
         odeplot(Y,[t,100*rhs(Skydiver)],0..T,numpoints=1000,color=blue,legend="Acceleration [*100m/s^2]")],
         thickness=2,title="Dynamics of the skydiver",labels=["",""]); Height, Velocity and Acceleration of the Skydiver as a functions of time.

> odeplot(Y,[y(t),diff(y(t),t),rhs(Skydiver)],0..T,numpoints=1000,color=red,thickness=3,labels=["y","y'","y''"],axes=boxed,title="Dynamics of the skydiver"); The same result as above, as a spacecurve.

1.4. Ballistics of the Gun Projectile - Czech self-propelled howitzer DANA

> restart; with(plots):

Warning, the name changecoords has been redefined

> read `Equations.sav`; Equations of movement of the projectile - see subsection 1.1. This system can be solved only numerically. We shall try to find the elevation angle leading to the maximal range of the projectile.

Czech self-propelled howitzer DANA

> g:=9.81: v:=1500.0: Cx:=0.1: d:=0.15: L:=0.6: Rho:=2500: S:=evalf(Pi*d^2/4): m:=S*L*Rho: rho:=1.22: g = gravitational acceleration, v = projectile initial velocity, Cx = coefficient of the projectile's aerodynamic resistance, d = projectile's diameter, Rho = projectile's mass density, S = projectile's cross section, m = mass of the projectile, rho = density of the atmosphere

> Ini:=x(0)=0,y(0)=0,D(x)(0)=v*cos(alpha),D(y)(0)=v*sin(alpha); Initial conditions. = elevation angle

> alpha:=evalf(convert(15*degrees,radians)); da:=2*alpha; da = elevation angle increment

> D1:=0: D2:=1: tau:=60: D1, D2 = variables controlling numerical procedure computing maximal range, = initial value of time

>    DMS:=proc(r) local rf, d, m, s;
  rf:=evalf(r*180/Pi);
  d:=floor(rf);
  rf:=frac(rf)*60;
  m:=floor(rf);
  s:=frac(rf)*60;
  cat(``,`Elevation = `,convert(d,string),`� `,
      convert(m,string),`' `,convert(s,string),`'' `);
end; Procedure converting in radians into degrees, minutes and seconds

$DMS := proc (r) local rf, d, m, s; rf := evalf(180*r/Pi); d := floor(rf); rf := 60*frac(rf); m := floor(rf); s := 60*frac(rf); cat(``,`Elevation = `,convert(d,string),`� `,convert(m,string),`' `,conver...$
$DMS := proc (r) local rf, d, m, s; rf := evalf(180*r/Pi); d := floor(rf); rf := 60*frac(rf); m := floor(rf); s := 60*frac(rf); cat(``,`Elevation = `,convert(d,string),`� `,convert(m,string),`' `,conver...$
$DMS := proc (r) local rf, d, m, s; rf := evalf(180*r/Pi); d := floor(rf); rf := 60*frac(rf); m := floor(rf); s := 60*frac(rf); cat(``,`Elevation = `,convert(d,string),`� `,convert(m,string),`' `,conver...$
$DMS := proc (r) local rf, d, m, s; rf := evalf(180*r/Pi); d := floor(rf); rf := 60*frac(rf); m := floor(rf); s := 60*frac(rf); cat(``,`Elevation = `,convert(d,string),`� `,convert(m,string),`' `,conver...$
$DMS := proc (r) local rf, d, m, s; rf := evalf(180*r/Pi); d := floor(rf); rf := 60*frac(rf); m := floor(rf); s := 60*frac(rf); cat(``,`Elevation = `,convert(d,string),`� `,convert(m,string),`' `,conver...$
$DMS := proc (r) local rf, d, m, s; rf := evalf(180*r/Pi); d := floor(rf); rf := 60*frac(rf); m := floor(rf); s := 60*frac(rf); cat(``,`Elevation = `,convert(d,string),`� `,convert(m,string),`' `,conver...$
$DMS := proc (r) local rf, d, m, s; rf := evalf(180*r/Pi); d := floor(rf); rf := 60*frac(rf); m := floor(rf); s := 60*frac(rf); cat(``,`Elevation = `,convert(d,string),`� `,convert(m,string),`' `,conver...$
$DMS := proc (r) local rf, d, m, s; rf := evalf(180*r/Pi); d := floor(rf); rf := 60*frac(rf); m := floor(rf); s := 60*frac(rf); cat(``,`Elevation = `,convert(d,string),`� `,convert(m,string),`' `,conver...$
$DMS := proc (r) local rf, d, m, s; rf := evalf(180*r/Pi); d := floor(rf); rf := 60*frac(rf); m := floor(rf); s := 60*frac(rf); cat(``,`Elevation = `,convert(d,string),`� `,convert(m,string),`' `,conver...$

>    while da>10^(-Digits-1) do;
  D1:=0; D2:=1;
  da:=da*0.3;
  while D2>D1 do;
    D1:=D2;
    alpha:=alpha+da;
    Ns:=dsolve({Ini,Ax,Ay},{x(t),y(t)},numeric);
    dt:=1:
    while abs(dt)>0.000001 do; dt:=subs(Ns(tau),-y(t)/diff(y(t),t)); tau:=tau+dt; od;
    D2:=subs(Ns(tau),x(t));
  od;
  Alpha:=alpha-da;
  NS:=Ns;
  print(DMS(Alpha),Range=D1,Time=tau);
  alpha:=alpha-2.25*da;
od: Loop computing maximal range. Range is computed for each elevation angle by the very similar procedure as in subsection 1.3. Thereafter elevation angle increases by da and new range is computed until new range is shorter than or equal as old one. In the case that increment is scaled down and new range is compared with the old one, until increment is greater than or equal as desired precision.

> plots[odeplot](NS,[x(t),y(t)],0..tau,scaling=constrained,
title="Trajectory",legend="rho=Const");P1:=%: Ballistic curve of the maximal range for the constant atmospheric density.

Atmospheric density is height dependent

> rho:=rho0*exp(-M*g*y(t)/R/T); Function describing atmospheric density as the function of height.

> M:=(28*4/5+32/5)*0.001: R:=8.315: T:=293.15: rho0:=1.3: M = atmospheric molar mass = 1/4 of and 4/5 of . R = universal gas constant, T = atmospheric temperature in Kelvin, rho0 = atmospheric density at zero height

> alpha:=evalf(convert(15*degrees,radians)); da:=2*alpha;

> P2:=plots[odeplot](Ns,[x(t),y(t)],0..tau,scaling=constrained,
color=blue,legend="rho=rho(h)"): P2; Ballistic curve of the maximal range for the height dependent atmospheric density.

No atmospheric resistance = vacuum

> Ax:=diff(x(t),t,t)=0; Ay:=diff(y(t),t,t)=-G; Equations of movement in vacuum

> alpha:='alpha':

> Ini:=x(0)=0,D(x)(0)=V*cos(alpha),y(0)=0,D(y)(0)=V*sin(alpha); Initial condition

> Sol:=subs(x(t)=x,y(t)=y,dsolve({Ax,Ay,Ini},{x(t),y(t)})); Analytical solution

$Sol := {x = V*cos(alpha)*t, y = -1/2*G*t^2+V*sin(alpha)*t}$

> Trajectory:=subs(solve(select(has,Sol,x),t),select(has,Sol,y))[]; Equation of the trajectory - can be obtained by the elimination of the parameter t .

> solve(subs(y=0,Trajectory),x); The y=0 if x = range.

> Range:=combine(select(has,{%},alpha)[],sincos); Select and simplify nontrivial solution

> Elevation:=solve(diff(Range,alpha)=0,alpha); Find elevation for the maximal possible range.

> Range:=value(subs(alpha=Elevation,V=v,G=g,Range)); Numerical value

> Trajectory:=value(subs(alpha=Elevation,V=v,G=g,rhs(Trajectory))); Trajectory equation.

> P3:=plot(Trajectory,x=0..Range,color=black): P3; Trajectory in vacuum.

> plots[display]([P1,P2,P3],labels=["x [m]","y [m]"]); Graphical comparison of all cases.

2. Dynamics of the space ship

2.1. Space ship on the high orbit

High orbit = no atmospheric resistance

> restart; with(plots):

Warning, the name changecoords has been redefined

> A:=kappa*M/R^2; Newton's law of the gravity. M = mass of the homogenous sphere, or of the sphere where mass density depends only on the distance from its center , later mass of the Earth, = Newton's gravitational constant, A = gravitational acceleration, Rz = Earth radius

> R:=sqrt(x(t)^2+y(t)^2);

> ax:=A*x(t)/R; ay:=A*y(t)/R; Components of the vector A in x and y axis

> Ax:=diff(x(t),t,t)=-ax; Ay:=diff(y(t),t,t)=-ay; The equations of the movement for the Earth satellite moving in the vacuum - no resistance. This set of equations has no analytical solution

> Values:=M=6.0*10^24,kappa=6.67*10^(-11),Rz=6378000.0;
assign(Values); Numerical values in SI units

> Ini:=x(0)=Rz+300000.0,D(x)(0)=0,y(0)=0,D(y)(0)=8500.0; Initial conditions for the numerical solution

> Ns:=dsolve({Ini,Ax,Ay},{x(t),y(t)},numeric);

> Ns(7500); Numerical solution returns components of the positional and velocity vectors

> Z:=plots[polarplot](Rz,phi=0..2*Pi,color=blue):

> plots[odeplot](Ns,[x(t),y(t)],0..7500,
labels=["x [m]","y [m]"],title="Trajectory"); Plot trajectory for 7500 seconds - as we can see this value is a shorter than one orbit period

> tau:=7500.0; supporting variables

>    Tau:=proc(T) local dt, tau;
  dt:=1.0;
  tau:=T;
  while abs(dt)>0.0001 do;
    dt:=subs(Ns(tau),-y(t)/diff(y(t),t));
    tau:=tau+dt;
  od;
  tau;
end;   Procedure computing duration of one orbit period. It is based on Newton's iteration, very similar to the loop used in 1.2

> tau:=Tau(tau);

> Dx:=subs(Ns(0),x(t))-subs(Ns(tau),x(t));
Dy:=subs(Ns(0),y(t))-subs(Ns(tau),y(t));
DR:=sqrt(Dx^2+Dy^2);
DVx:=subs(Ns(0),diff(x(t),t))-subs(Ns(tau),diff(x(t),t));
DVy:=subs(Ns(0),diff(y(t),t))-subs(Ns(tau),diff(y(t),t));
DV:=sqrt(DVx^2+DVy^2); As we can see, result is very precise, because difference in position after one orbit is lower than 24m and difference in velocity is lower than 2.5 cm/s .

> Tr:=odeplot(Ns,[x(t),y(t)],0..tau):

> display({Z,Tr},scaling=constrained,
labels=["x [m]","y [m]"],title="One orbit"); Trajectory and the Earh

> Pv:=[x(t),y(t),0];
Vv:=[diff(x(t),t),diff(y(t),t),0]; Components of the positional and velocity vectors in space, the coordinate z = 0 is added.

> Sv:=linalg[crossprod](Pv,Vv/2); Momentum of the velocity = Surface velocity - surface covered by the positional vector during one second. This is exact formulation for the term introduced by Johannes Keppler at his famous second law = the surface velocity is constant

> SV:=subs(Ns(0),[t,Sv[3]][2]); The initial value of the surface velocity of the given satellite

> odeplot(Ns,[t,Sv[3]],0..tau,axes=boxed,
labels=["t [s]","x*y'-x'*y"],title="Momentum=const?"); Was Keppler right? This is not constant!

> odeplot(Ns,[t,Sv[3]],0..tau,view=[0..tau,0.9*SV..1.1*SV],axes=boxed,
labels=["t [s]","x*y'-x'*y"],title="Momentum=const!"); Yes, it is constant, the preceding plot shows only rounding errors!

2.2. Space ship on the low orbit

Low orbit = it is necessary to assume impact of the aerodynamic resistance

> alias(He=Heaviside);
plot(He(x),x=-1..1,axes=boxed); Abbreviation for the Heaviside function

> h:=sqrt(x(t)^2+y(t)^2)-Rz; h = height above Earth surface

> Values:=S=evalf(Pi*2.5^2),Cx=0.6,m=2000.0;
assign(Values); Numerical values, S = cross section of the satellite, Cx = coefficient of the aerodynamics resistance, m = mass of the satellite, all in SI units

> rho:=unapply(exp((.6392146930e-10*u^2-.1447577359e-3*u+.3316213319)*He(104413.9387-u)+
(-35.02764072*u+1412587.921)/(u+54946.54461)*He(-104413.9387+u)),u); Function describing = density of the Earth atmosphere as a function of the height above Earth surface. Fubction's coefficients were computed very similarly to http://adept.maplesoft.com/categories/data_analysis_stats/regression/html/nonlinearpiecewise.html , using data from CRC Handbook of Chemistry and Physics, D.R.Linde editor, CRC Press 72nd edition, tables 14-13 and 14-14, ISBN-0-8493-0565-9

> Ax:=diff(x(t),t,t)=-ax-
1/2*S*rho(h)*Cx/m*sqrt(diff(x(t),t)^2+diff(y(t),t)^2)*diff(x(t),t);
Ay:=diff(y(t),t,t)=-ay-
1/2*S*rho(h)*Cx/m*sqrt(diff(x(t),t)^2+diff(y(t),t)^2)*diff(y(t),t): Equations of the movement of the Earth satellite moving on the low orbit. Equations are too long, so output of Ay is not presented - it is very similar to Ax . These equations have no analytical solution.

> A:=sqrt(rhs(Ax)^2+rhs(Ay)^2): The absolute value of the acceleration of the satellite.

> Ag:=A/(kappa*M/(x(t)^2+y(t)^2))-1: Acceleration recomputed into g units

> V:=sqrt(diff(x(t),t)^2+diff(y(t),t)^2); The absolute value of the velocity vector.

> Ini:=x(0)=Rz+1000000.0,D(x)(0)=0,y(0)=0,D(y)(0)=7100.0; Initial conditions

> Ns:=dsolve({Ini,Ax,Ay},{x(t),y(t)},numeric); The numerical solution

> subs(Ns(3000),h); subs(Ns(2000),h); Looking for initial value for the iteration of the landing time. is sufficient.

>    Tau:=proc(T) local dt, tau;
  dt:=1.0;
  tau:=T;
  while abs(dt)>0.0001 do;
    dt:=subs(Ns(tau),-h/diff(h,t));
    tau:=tau+dt;
  od;
  tau;
end;    Procedure computing landing time . If spaceship lands its height h = 0 .. It is based on Newton's iteration, very similar to the procedure used in 2.1.

> tau:=Tau(3000.);

> Z:=polarplot(1,phi=0..2*Pi,color=blue,numpoints=300):

> P1:=display({Z,odeplot(Ns,[x(t)/Rz,y(t)/Rz],0..tau,numpoints=1000)},scaling=constrained,
labels=["x","y"],labels=["x [Rz]","y [Rz]"],title="Trajectory"): P1; Landing trajectory

> display(P1,view=[-0.95..-0.88,0.34..0.45],scaling=constrained); Detail of the landing trajectory

> odeplot(Ns,[t,h/1000.],0..tau,labels=["t [s]","h [km]"],title="h(t)"); Height as a function of time

> odeplot(Ns,[t,h/1000.],2400..tau,labels=["t [s]","h [km]"],title="h(t)"); Height as a function of time - detail

> odeplot(Ns,[t,V/1000],0..tau,numpoints=1000,
labels=["t [s]","V [km/s]"],title="V(t)"); Velocity as a function of time

> odeplot(Ns,[h/1000,V/1000],0..tau,numpoints=1000,
labels=["h(t) [km]","V [km/s]"],title="V(h)"); Velocity as a function of height

> odeplot(Ns,[t,Ag],2000..tau,axes=framed,numpoints=1000,
labels=["t [s]","A [g]"],title="A(t)"); Acceleration as a function of time

> odeplot(Ns,[V/1000,Ag],2000..tau,numpoints=1000,axes=framed,
labels=["V [km/s]","A [g]"],title="A(V)"); Acceleration as a function of velocity

> odeplot(Ns,[h/1000,Ag],2000..tau,axes=framed,numpoints=1000,
labels=["h [km]","A [g]"],title="A(h)"); Acceleration as a function of height

> odeplot(Ns,[h/1000,V/1000,Ag],2300..tau,axes=framed,color=black,thickness=2,orientation=[150,70],
labels=["h [km]","V [km/s]","A [g]"],numpoints=1000,title="Dynamics of landing"); Plot all three variables h(t), V(t) , and A(t) in one 3D plot - as a spacecurve

> save Ax, Ay, Ag, V, h, rho, Rz, Tau, `spaceship.sav`; Save these variables for later use

2.3. Aerodynamics retarding of 2nd cosmic velocity

> restart; with(plots):

Warning, the name changecoords has been redefined

> read `spaceship.sav`: Read variables and procedures created in the preceding part

> alpha1:=evalf(convert(113.26*degrees,radians)); W:=12200.;

> Ini:=x(0)=Rz+1000000.,y(0)=0,D(x)(0)=W*cos(alpha1),D(y)(0)=W*sin(alpha1); Initial conditions

> Ns:=dsolve({Ini,Ax,Ay},{x(t),y(t)},numeric); Numerical solution

> Z:=polarplot(1,phi=0..2*Pi,color=blue,numpoints=400): Earth

> P1:=odeplot(Ns,[x(t)/Rz,y(t)/Rz],-2000..97000,numpoints=1000,color=red,thickness=3): Trajectory

> display({P1,Z},labels=["x [Rz]","y [Rz]"],title="Trajectory",scaling=constrained); Plot of the Earth and trajectory

> display({Z,P1},view=[0.35..1.25,-0.2..1],scaling=constrained,labels=["x [Rz]","y [Rz]"]); Detail of the trajectory

> display({Z,P1},view=[0.4..0.6,0.8..0.95],scaling=constrained,labels=["x [Rz]","y [Rz]"]); Detail of the landing zone

> subs(Ns(92000),h); subs(Ns(93000),h); Looking for the initial value for iteration of the landing time

> tau:=Tau(93000); Landing time

> plot([seq(evalf(subs(Ns(T),[t,Ag])),T=[$250..500])],axes=framed,
labels=["t [s]","A [g]"],title="First braking, A(t)"); First braking

> plot([seq(evalf(subs(Ns(T),[t,Ag])),T=[$tau-1000..tau])],axes=framed,
labels=["t [s]","A [g]"],title="Second braking, A(t)"); Second braking - after one orbit

>    display({plot([seq(evalf(subs(Ns(T),[V/1000,Ag])),T=[$250..500])],color=blue),
         plot([seq(evalf(subs(Ns(T),[V/1000,Ag])),T=[$tau-1000..tau])],color=red)},
         labels=["V [km/s]","A [g]"],title="Braking, A(V)",axes=framed); Deceleration as a function of the velocity

>    display({plot([seq(evalf(subs(Ns(T),[h/1000,Ag])),T=[$250..500])],color=blue),
         plot([seq(evalf(subs(Ns(T),[h/1000,Ag])),T=[$tau-1000..tau])],color=red)},
         labels=["h [km]","A [g]"],title="Braking, A(h)",view=[0..100,-1..7.5],axes=framed); Braking as a function of the height, first braking , second braking

>    display({odeplot(Ns,[h/1000,V/1000],150..550,color=blue),
         odeplot(Ns,[h/1000,V/1000],tau-1000..tau,color=red)},
         labels=["h [km]","V [km/s]"],title="V(h)"); Velocity as a function of the height, first braking , secon