kill(all)$
load(lrats)$
load(vect)$

orbit equation (after burnout,'vb,rb,β' ), version 1.0
peter.vlasschaert@gmail.com,27/08/2018.

'vb,rb,β' 

rem :  distance (I,E) = hs

 polar coordinates. (2D)

see: transformation of two massproblem:total energy.(website)

p1:('diff(r,t,2) -r*('diff(θ,t))^2)*m = -μ/r^2*m;
p2:1/r*'diff(r^2*'diff(θ,t),t)          = 0;

rem : Newton's law attraction :  F = -(m*μ)/r^2
=>  μ = G*Me                                     
       Me = Mass earth .
       m   = Mass  object around  Me.
       G   = Gravitation constant.
       r    = Distance between two masses (Me,m)

p3:part(p2,1,1);
p3a:part(p2,1,2);

rem : 1/r ≠ 0 , 'diff((r^2*('diff(θ,t,1))),t,1) = 0

p3b:p3=0;

p3c: part(p3b,1,1) = h;

p3d:solve(p3c,('diff(θ,t,1)));

p3e:p3d[1];

rem : p3c = h =  'constant angular momentum'
=======================================
rem : result from 'p1' and use of 'p3e'

p4:subst(p3e,p1)/m;

rem : insert  'p3e' into 'p1'.

subst : ' u= 1/r ' =>   'r = 1/u '

p4a:u=1/r;
p4b:solve(p4a,r)$
p4c:p4b[1];

 find : dr/dt  ?
===========

depends([u],t);

p5:'diff(r,t) = subst(p4b,'diff(r,t));

p5a:lhs(p5) =ev(rhs(p5),diff);

rem : chian rule. (u,t) , θ
====================

p5b:['diff(u,t) = 'diff(u,θ)*'diff(θ,t)];

p5c:subst(p5b,p5a);

rem : insert 'p3e'
==============

assume ( h>0);

p5d:lratsubst([p3e,p4a],p5c)$
p5e:lhs(p5d) = subst([rhs(p4a)=lhs(p4a)],rhs(p5d));

find : diff(r,t,2) ? d/dt(dr/dt)
======================

rem :chain rule

p6:'diff(r,t,2) = diff(rhs(p5e),t);

p6a:lhs(p6) =-h*'diff(('diff(u,θ,1)),θ)*'diff(θ,t);

p6b:subst(p3d,p6a);

p7b:lhs(p6b)=lratsubst([rhs(p4a)=lhs(p4a)],rhs(p6b));

find u(θ),ODE
===========

p8:lratsubst([p7b,rhs(p4a)=lhs(p4a)],p4);

1e) find coef :  'diff(u,θ,2)
--------------------------------------

p8a:ratcoeff(lhs(p8),('diff(u,θ,2)));

2e) divide :  ODE by 'p8a'
--------------------------------------

p8b:ratsimp(p8/p8a);

3e) solve : ODE : u = u(θ)
--------------------------------------

use : ode2  ' maxima '

p8c:ode2(p8b,u,θ);

p8d:subst([%k1=c1,%k2=c2],p8c);

rem : c1,c2 => unknown

rem :  C*cos(θ-θ[0]) = c1*sin(θ)+c2*cos(θ)

rem : C,θ[0] => unknow

p8e:lhs(p8d) =part(p8d,2,1)+C*cos(θ-θ[0]);

rem : u = 1/r , p= h^2/μ .
===================

p8f:subst([p4a],p8e);

p8g:solve(p8f,r)[1];

p8h:p = h^2/μ;
p8i1:part(p8g,2,1)/μ;
p8i2 :expand(part(p8g,2,2)/μ);

p9:lhs(p8g) = lhs(p8h)/p8i2;

rem : def 'ε'
=========

p9a:ε =ratcoeff(p8i2,cos(θ-θ[0]));

result :  'r'
========

p10:lratsubst([rhs(p9a)=lhs(p9a)],p9);

depends([θ],t);

calculation : Energy per unit mass (= m)
==================================

p11: E = T+U;

rem : T = kinetic energy per unit mass (= m)
           U = potential energy per unit mass(= m)

velocity  :  vb     = burnout velocity        ,' 8 a 9'  Km /s 
distance :  rb     = distance @ burnout .
                     R[e] = radius of earth
                     hs    = height above sea level.

p11a:μ = G*Me;
p11b: rb = R[e] +hs;

p11c: T = vb^2/2;
p11d: U= -1*μ/rb;

p12 : subst([p11c,p11d],p11)$
p12a:subst([p11a,p11b],p12);

end of rocket burn : see E , 'p12a'

find relation between : ε = ε (E)
==========================

math : v = dr/dt *e[r] +r*dθ/dt*e[θ]
===========================
rem : e[r] perpendicular to e[θ] 
===========================
 => v = 'vector'
===========================

calculation :  ' T ' per unit mass,   1/2*v^2
==========

1e)  dr/dt from ,  'p10'
--------------------------------

p13: 'diff(r,t) =diff(rhs(p10),t)$
p13a:lhs(p13) =subst(['diff(θ[0],t,1)=0],rhs(p13));

p13b : T = 1/2*(lhs(p13a)^2+(r*'diff(θ,t))^2);

p13c:subst([p13a],p13b);

p13d:subst([p3e],p13c);

p13e:[(rhs(p10))^2=(lhs(p10))^2];

p13f1:part(p13d,2,1,1);
p13f11:lratsubst(p13e,p13f1);
p13f111:factor(p13f11);
p13f2:part(p13d,2,1,2);
p13f3:part(p13d,2,2);

 result : T
========

p14:lhs(p13d) = (p13f111+p13f2)/p13f3;

p14a:solve(p10,p)$
p14b:factor(p14a)[1]$
p14c:(rhs(p14b)=lhs(p14b))^2;
p14d:[p14c];
p14e:subst(p14d,p13f111);
p14f:subst([p8h],p14e);

p15:1/p13e[1];

p15a: subst([rhs(p15) = lhs(p15)],p13f2);

insert  T:

rhs : T

p16:1/p13f3*(p14e+p15a);
p16a:factor(p16);
p16b:trigsimp(p16a);
p16c:factor(p16b);

final result :  T
==========

p17:lhs(p14) = p16c;

We want to find : U
================

p18:U = - μ/r;

p18a: [1/p10];

p18b:subst(p18a,p18);

p18c:solve(p8h,μ);

p18d:subst(p18c,p18b);

total energy : E  ?
===============

p19:lhs(p11) = rhs(p17)+rhs(p18d);

p19a:expand(p19);

p19b:solve(p19a,ε);

p19c:p19b[2];

What about : h ?
==============
--------------------------------------------------------------------------------------------------------------
h = r*V[θ] = rb*vb*cos(angle between local horizon and burnout velocity)
angle of attack = angle between local horizon and burnout velocity = β
============================================================
rem : μ = G*Me = 3.986 *10^14 m^3/s^2
============================================================

p20:h=rb*vb*cos(β);
p20a: p8h;
p20b:subst([p20,p11a,p11b],p20a);

p20c:subst([p20b,p20],p19c);

p20d:subst([p12a],p20c);

assume(G>0,R[e]>0,Me>0);

p20e:ratsimp(p20d);

How to find  : θ[0] ?
================
-----------------------------
for : r=rb => θ = 0
-----------------------------

p21:p10;

p21a:[r=rb,θ=0];

p21b:lratsubst(p21a,p21);

p21c:solve(p21b,θ[0]);

p21d:subst([p20b,p20e],p21c);

p21e:p21d[1];

relation between  E and ε
=====================

Ellipse        : ε < 1, E < 0
Parabola    : ε = 1, E = 0
Hyperbola :  ε > 1, E> 0
rem : circle  :ε =0

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

elliptical orbit
============

 How to find : (ap)ogee , (pe)rigee.
==========================

p22:2*a = r[pe]+r[ap];

p22a: θ - θ[0] = α;
p22b:subst([p22a],p10);

see :  p, ε  , see above

 perigee
=======
r=r[pe] = rmin  , α = 0

p22c:subst([r=r[pe],α =0],p22b);

 apogee
=======
r=r[ap] = rmax  , α = π

p22d:subst([r=r[ap],α =π],p22b);

find  a : ? 'semi major axis'
======================

p23:solve(p22,a)[1];

p23a:subst([p22c,p22d],p23);

p23b:expand(p23a);

p23c:ratsimp(p23b);

p23d:combine(solve(p23c,p));

find b : ? 'semi minor axis'
======================

p24:b=a*(1-ε^2)^(1/2);

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