kill(all)$
load(lrats)$

Hohmann transfer , version 1.0
peter.vlasschaert@gmail.com,17/08/2018.

info : https://en.wikipedia.org/wiki/Hohmann_transfer_orbit

 Hohmann Transfer .

Legend : Figure 1
===============

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

vector : u  =  v1+Δv1     , velocity,(circle ,ellipse)
                v  =  v2 +(Δv2) ,velocity,(circle ,ellipse)
                w                        ,velocity,(transfer orbit,'e')
point   :  F = peripoint.
                G = apopoint
                M = point on the transfer orbit 'e'
                C = point on the big        circle 'd'
                D = point on  the small   circle 'c'
                E = point on  the small   circle 'c'
circle     : r1 ( = radius )= small circle , c
                  r2 ( = radius) = big     circle , d
ellipse  : e

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

circular orbit velocity. v[cir]

rem : angle between ( R, v)  = π/2
           =======================

simple math : ?
=============

p1:θ = ω.t;
p2:   x = R*cos(rhs(p1));
p2a: y = R*sin(rhs(p1));

rem : figure 2, (g,f)  = (x,y).

p3:'diff(x,t)=diff(rhs(p2),t);
p3a:'diff(y,t) =diff(rhs(p2a),t);

rem : figure 2, (p3,p3a)  = (vx,vy).
===========================
           =>        v=sqrt(vx^2+vy^2).
-------------------------------------------------

p4:(p3)^2;
p4a:(p3a)^2;

p5:v^2=rhs(p4)+rhs(p4a);

assume(R>0,ω>0);

p6:sqrt(trigsimp(p5));

p7:'diff(lhs(p3),t) =diff(rhs(p3),t);
p7a:'diff(lhs(p3a),t) =diff(rhs(p3a),t);

rem : figure 2, (p7,p7a)  = (ax,ay).
===========================
           =>        a=sqrt(ax^2+ay^2).
------------------------------------------------

p8:a^2=(rhs(p7))^2+(rhs(p7a))^2;

p9:sqrt(trigsimp(p8));

acceleration : a

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

p10a:lhs(p10)=lratsubst([rhs(p6)=lhs(p6)],rhs(p10)*R)/R;

? find : v[cir].

1e)  second law newton =  m*m[e]*G/R^2
2e)  m*v^2/R
  =>  1e) = 2e)  condition ,  ' circular orbit '

p11:F = m*m[e]*G/R^2;
p11a: rhs(p10a)*m=rhs(p11);

p12:subst([v=v[cir]],p11a);

p13:solve(p12,v[cir]);

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

m[e] = mass earth.
m      = mass  object( = sattelite,planet).
R       = distance (earth to objct).
G       = gravitation constant.

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

p13a :p13[2];

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

find : Δv1,Δv2 , for Hohmann Transfer.  ?
--------------------------------------------------------------

=========================================
=>  solve the problem we need '  two equations '.
=========================================

rem : Hohmann Transfer from circular orbit to elliptical orbit.

What are the assumption to do the calculation . ' below'
------------------------------------------------------------------------------------

1e) orbits : conics.
2e) total : angular momentum and total energy is conserved.

see figure 1 :

peripoint  ( = r1 ) , vp = v1+Δv1 ( peri -speed)
apopoint  ( = r2 ) , va = v2           ( apo -speed)

p14:  vp =v1+Δv1;
p14a : va = v2;

angular momentum: 'first equation'
-----------------------------------------------------
-----------------------------------------------------

angular momentum : simple because ' angle (r,v) ' =π/2

p15:m*r1*(rhs(p14)) *sin(π/2)=m* r2*(rhs(p14a))*sin(π/2);

p15a:solve(p15,v2);

total energy 'E' :  'kinetic energy + potential energy = constant '
----------------------------------------------------------------------------------------------
-----------------------------------------------------------------------------------------------

E = 1/2*m *v^2 -G*m*m[e]/r 
-----------------------------------------------------------
E < 0 , elliptical      orbit  or circular orbit.
E = 0 , parablic      orbit ( = escape  orbit)
E > 0 , hyperbolic orbit  ( = escape orbit)

 at peripoint : ' E ' = energy
----------------------------------------

p16:1/2*m*(rhs(p14))^2-G*m[E]*m/r1;

at apopoint : ' E ' = energy
----------------------------------------

use : ' p15a ' for v2

p17: 1/2*m*v2^2-G*m[e]*m/r2;
p17a:subst(p15a,p17);

E = constant  ?

p18:p16=p17a;

solve for : Δv1

p19:solve(p18,Δv1);

p19a:p19[1];

rem :  v1= v[circ] , R= r1
            v2 =v[circ] , R= r2
====================
see above : ' p13a '

p20: subst([R=r1,v[cir]=v1],p13a);

p20a:subst([p20],p19a);

Δv1 : extra thruster burner of rocket to leave orbit ( r1).

=>  we need extra : Δv2 , so we don't fall back to circular orbit  ( r1 )

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