kill(all)$

rocket launch,missile flight.(freeflight)  'version 1.0'
peter.vlasschaert@gmail.com,10/06/2018

rem : we take 2D  version of the earth  'circle'.

see : fundamentals of astrodynamics (Roger R.Bate,etc....)

depends([x,y],t);

freeflight : start after burningout
initial position ,velocity.

initial position : ' above the the earth' = h0
r0 = rE + h0 
x0 = r0*cos(θ0)
y0=  r0*sin(θ0)

rem : angle be measured from equator = x-axes

p1:r_0=r_E+h_0;
p2:x_0 = r_0*cos(θ_0);
p3:y_0= r_0*sin(θ_0);

initial velocity : 'above the earth ' = h0
vx0 =V*cos(π/2-γ+θ0)
vy0= V*sin(π/2-γ+θ0)
γ = angle between horizontal  local horizon ( = θ  direction).=> V*cos(γ)
vx0 : angle between x-axes and speed = V ( => θ0+π/2-γ)     'cos'
vy0 : angle between x-axes and speed = V ( => θ0+π/2-γ)     'sin'

p5:vx_0=V*cos(%pi/2-γ+θ_0);
p6:vy_0=V*sin (%pi/2-γ+θ_0);

relation : between 'i1,j1' and  'i2,j2' 
rotation : θ  =>  θ+π/2

p8:i2=cos(θ)*i1+sin(θ)*j1;
p9:j2=cos(θ+%pi/2)*i1+sin(θ+%pi/2)*j1;

the rate change : new coordinates : differentation respect to θ
i2:i2(θ),j2:j2(θ).   'unit vectors'

depends([i2,j2],θ);

p10:diff(p8,θ,1);
p11:factor(diff(p9,θ,1));

p12:ratsubst(lhs(p9),rhs(p9),p10);
p13:ratsubst(lhs(p8),rhs(p8),p11);

rate of change : 'new coordinates' as function of time : chain rule.
reason : θ will change in time with direction.

depends([θ,r],t);

p14:'diff(i2,t) = diff(θ,t,1)*rhs(p12);
p15:'diff(j2,t) = diff(θ,t,1)*rhs(p13);

rem : R= i2*r, R = vector, r=magnitude,i2 unit vector same direction of

velocity :

p16:'diff(R,t,1) = diff(i2*r,t,1);

p17:ratsubst(rhs(p12),lhs(p12),p16);

acceleration :

p18:'diff(R,t,2) = diff(rhs(p17),t,1);

∅;

use :lratsubst;(need to load packages).

load("lrats")$

p19:lratsubst([p12,p13],p18);

above 'p12' ,two components i2,j2

p20a:coeff(rhs(p19),i2);
p20b:coeff(rhs(p19),j2);

rem : 'p20a' component  'unit vector : i2'
            i2 : same direction  ' distance to rocket' (Newton's Law of attraction)
           'p20b' component  'unit vector :j2'
            j2 :  direction orthogonal i2 not direction of r .

Newton's of Law attraction.
=======================
F = -(G.m1.m2)/r^2*i2
i2: unit vector in direction of r
m1 : m_r=mass (rocket)
m2 : M_e=mass (earth)

p21a : F=-(G*m1*m2)/r^2*i2;
p21b : m1=m_r;
p21c : m2=M_e;

p21d:sublis([p21b,p21c],p21a);

component : i2, 'p20a,p21d'

p21e:  (m_r*p20a*i2=rhs(p21d))/i2;

assume(m_r>0);

assume (m_r>0);

p21f:ratsimp(p21e)/m_r;

p21g:ratsimp(p21f);

component: j2

p22:p20b=0;

differential equations : we want solve both 'p21g,p22'
=============================================
with :                    initial conditions.(distance,velocities)

p23:p21g;
p24:p22;

We want to find : 'p24' 

p25a:r^2*diff(θ,t,1);
p25b:factor(diff(p25a,t) = rhs(p24));

check : 'p24,p25b'
reason : r≠0

p26:lhs(p24)=part(p25b,1,2);

p27:is(p26);

We want solve : both ode 2
=======================

p28:p23;
p29:p24;

because : ' p25a' = constante
=> constante of integration'initial': 
r*(r*diff(θ,t))=moment arm*tangential velocity
reason : dR/dt = dr/dt*i2+r*dθ/dt*j2

p30a :p23;
p30b :p25a = lhs(p1)*V*cos(γ);

We want to set :
p=V*r_0*cos(γ)
w = 1/r

p : ?

p31:p=rhs(p30b);

assume(p>0);

w:1/r ?

p32:w=1/r;

we want  ,  2 ODE2  => 1 ODE2  (in w)
subst: w=1/r

1e) chain rule : rem : w:w(θ)

depends([w],θ);

p33:'diff(w,t,1)  = diff(θ,t,1)*diff(w,θ,1);

2e)  'p' must into 'p33'

p34:lhs(p30b) = lhs(p31);

p35:solve(p34,diff(θ,t,1));

p36:lratsubst(p35,p33);

3e) chain rule : rem : w:w(r)

depends([w],r);

p37:'diff(w,t)=diff(w,r)*diff(r,t);

p38:lhs(p37) =sublis([p32],rhs(p37));

p39:lhs(p38)=ev(rhs(p38),diff);

4e) combine 'p36 = p39'

p40:rhs(p36) =rhs(p39);

How to solve : ODE2 ?

p41:solve(p40,diff(r,t))[1];

set  :  z: diff(r,t)

p41a: z=lhs(p41);

1e) chain rule : z  changes with respect 't'

depends([z],θ);

p42:'diff(z,t) = diff(θ,t)*diff(z,θ);

depends([w],θ);

p43:lhs(p41a) =rhs(p41);

p44:lhs(p42) = sublis([p43],rhs(p42));

p45:lratsubst(p35,p44);

p46:sublis([p41a],p45);

p47:lhs(p46)=-p^2/r^2*diff(w,θ,2);

We want to solve: ODE2:below with the above.

p48:p30a;

p49:lratsubst([p47,p35[1]],p48);

p50:expand(p49)*r^2;

p51:expand(p50);

p52:factor(p51)/p^2;

p53:expand(p52);

p54:lratsubst([rhs(p32)=lhs(p32)],p53);

p55:factor(p54)/(-1);

set :  Δ =r0*V^2/(G*Me) ? 

p56:sublis([p31],p55);

p57:factor(args(p56)[2]);

p58: Δ =r_0*V^2/(G*M_e);
p59:part(p57,1);
p60:part(p57,2);

p61:expand(p58*p59)*r_0;

p62:lratsubst([rhs(p61)=lhs(p61)],p56);

try : ODE2
r(0)        = r_0
dr(0)/dt = V*sin(γ)  :vertical horizon

pa:'r(0) =r_0;
pb:'diff(r(0),t) = V*sin(γ);

p63:ode2(p62,w,θ);

How to find : %k1,%k2

p64:diff(p63,θ,1);

p65:solve(p41,lhs(p64))[1];

p66:lratsubst([p31],p65);

p67:lratsubst([lhs(p41)=rhs(pb)],p66);

p67a:sublis([θ=0],rhs(p67)=rhs(p64));
p67b:solve(p67a,%k1)[1];

w(0) = 1/r(o),

p68:[p32];
p69:sublis(p68,p63);
p70:sublis([r=r_0,θ=0],p69);
p71:solve(p70,%k2)[1];

'p63' we put 'p67b,p71'

p72:sublis([p67b,p71],p63);
p73:ratsimp(p72);

p74:part(p73,2,1);
p75:-1*part(p73,2,1,1);
p76:part(p73,2,1,2);

p77:expand(p75);

coeff : = Δ after  / 'p76'

p78:coeff(p77,Δ)/p76;

 p78a:trigreduce(p78);

coeff : ≠ Δ after / 'p76'

p79:factor(ev(p77,Δ=0))/p76;

result :  ODE2 after use 'initial conditions'

p80:w=p78*Δ+p79;
p81:trigreduce(p80);

p82:sublis([p32],p80);
p83:lratsubst([sec(γ)=1/cos(γ)],p82);
p84:solve(p83,r)[1];

expand(Δ*cos(γ)*sin(γ)*sin(θ)+(1-Δ*cos(γ)^2)*cos(θ)-1);

p85a:part(p84,2,0);
p85b:part(p84,2,1,1);
p85c:part(p84,2,1);
p85d:part(p85c,1);
p85e:expand(-1*part(p85c,2));
;

coeff :Δ (=)

-----------------------------------------------------------------------------------------

p86:coeff(p85e,Δ);
p86a:trigreduce(factor(p86));

------------------------------------------------------------------------------------------

coeff*Δ

p86b:Δ*p86;

coeff: Δ(≠)

p87:ev(p85e,Δ=0);

p88:p87+p86b;

p89:expand(p88);

collectterms : cos(θ),sin(θ).

p90:collectterms(p89,cos(θ),sin(θ));

p91:args(p90);

find :coeff  'cos(θ)', 'p90'

p92:coeff(p90,cos(θ));

find :coeff  'sin(θ)', 'p90'

p93:coeff(p90,sin(θ));

     ? Ellipse 'polar coordinates'
----------------------------------------------------
--------------------------------------------------------------------------------------------------------
        r=(a*(1-e^2))/(1-e*cos(θ-ω))
--------------------------------------------------------------------------------------------------------
e  :eccentricity, e=0 :circle,0<e<1:ellipse ,e=1:parabola,e>1:hyperbola
a  : semimajor axis
Δ : velocity
ω : argument of apogee
θ : central angle
--------------------------------------------------------------------------------------------------------

γ : free-path angle  (angle between local horizon and velocity vector)
rem :Δ = 1       :circle motion
          0 < Δ < 2:elliptical motion
          Δ = 2       :parabolic motion
          Δ >2        :hyperbolic motion

p94:rr=(a*(1-e^2))/(1-e*cos(θ-ω));

p95:args(p94);

p95a:p95[2];

p95b:part(p95a,1);
p95c:part(p95a,2);

find : relation between  'e' and 'Δ,γ'
----------------------------------------------------

1e) general elliptic equation 'from above'

p95d:expand(trigexpand(p95c));

p95d1:collectterms(p95d,cos(θ),sin(θ));

p95d2:coeff(p95d1,cos(θ));
p95d3:coeff(p95d1,sin(θ));

2e)  ODEs : find other equation

p96a:p92;
p96b:p93;

relation : Δ,e,γ

p97a:p95d2=p92;
p97b:p95d3=p93;

1e)  take power 2 ,for both above.

p98a:(p97a)^2;
p98b:(p97b)^2;

p99a:expand(trigsimp(p98a+p98b));

p99b:sqrt(p99a);

rem :missile :  0<Δ<1 : 'projectile'

Impact : Equation  r = r_E, φ     'missile'
---------------------------------------------------

We use : φ : the total central angle .
s=distance :  'along the earth '
r_E               : 'radius of the earth'
φ                  : 'angle for impact '
elementray mathematics    :   s=r*θ

1e) use :  'p84'
============

p100a:p84;

2e) r=r_E, θ = φ   'impact'.
=====================

p100b:sublis([r=r_E,θ=φ],p100a);

3e)  Δ.
=====
Δ =    ( v_ob*r_ob) /(G*M_e),
 b=burnout the fuel to start the free flight .

p100c:p58;

4e)  velocity  V ?
==============

p100d:solve(p100b,r_0);

p100e:sublis(p100d,p100c);

p100f:solve(p100e,V);

p100g:p100f[2];

Flight time : 'free flight time ' t=t_free.
===============================

1e) use :  'p101a'; see above 'p30b', ODE2 , 'second equation'

p101a:p30b;

2e) use : 'p100a',  r => r^2

 p102a : (p100a)^2;

3e) integration  θ : 0 to φ , V 'see above'
to start the free flight .
==================================

1e) p103a: integrate(rhs(p102a), θ, 0, %pi);
2e) p103b: t_free =1/(V*r_0*cos(γ))*p103a;

assume(φ>0,Δ<2,Δ>0.3,γ<%pi/2);

                                                          ==========================================