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

the plane of the orbit from  3d to 2d ( = coordinate transformation) , version 1.0
peter.vlasschaert@gmail.com,21/09/2018.

rem : appendix : 'flight over a spherical earth '

rotation matrix , righ hand rule 'use'.
------------------------------------------------------

rotx(α):=matrix(
		[1,	0,	0],
		[0,	cos(α),	-sin(α)],
		[0,	sin(α),	cos(α)]
	);

roty(α):=matrix(
		[cos(α),	0,	sin(α)],
		[0,	1,	0],
		[-sin(α),	0,	cos(α)]
	);

rotz(α):=matrix(
		[cos(α),	-sin(α),	0],
		[sin(α),	cos(α),	0],
		[0,	0,	1]
	);

 use : fixed system X,Y,Z with  e=[e1,e2,e3] , unit vector.
----------------------------------------------------------------------------------

cross product :  C= R x dR/dt.   , C=[cx,cy,cz], kepelerian orbits  
----------------------------------------------------------------------------------------------
 =>  vector C : perpendicular to R and dR/dt
 =>  vector C : constant
----------------------------------------------------------------------------------------------
result :stays one plane (per unit mass), see website
----------------------------------------------------------------------------------------------

use : menu => algebra => enter matrix
===============================

p1:C= matrix(
 [e1,e2,e3], 
 [x,y,z], 
 ['diff(x,t),'diff(y,t),'diff(z,t)]
);

determinant  'C'
------------------------

p1a:determinant(rhs(p1));

p2a: [e1,e2,e3];
p2b: [cx,cy,cz];
p2c: lhs(p1) =sum(p2a[i]*p2b[i],i,1,3);

find  coef 'e1' ,coef 'e2',coef ' e3
-------------------------------------------------------

p3a:ratcoeff(rhs(p2c) = p1a,e1);
p3b:ratcoeff(rhs(p2c) = p1a,e2);
p3c:ratcoeff(rhs(p2c) = p1a,e3);

Let us now use the fact that the motion is planar and re-orient our coordinates in such
a way that  'one of the body’s coordinate ' = 0 
----------------------------------------------------------------------------------------------------------------------------------
This means that the body :
remains in the plane z = 0 of the coordinate system
suppose  : z=0,cz=c,x=ζ,y=η
---------------------------------------------------------------------------------------------------------------------------------
rem   : ζ  = zeta ,ξ = xi

p4:[z=0,cz=c,x=ζ,y=η];
p4a:subst(p4,p3a);
p4b:subst(p4,p3b);
p4c:subst(p4,p3c);

  d/dt(dR/dt) = -κ^2*R/r^3
=====================
κ = m1*m2*G = m[e]*m*G
-------------------------------------------------------
rem : m[e] = mass earth
          G       = gravity constant.
         m       = mass around big mass.
-------------------------------------------------------
rem : R : vector , r : radius

p5:'diff(R,t,2) = -κ^2*R/r^3;

insert : New coordinate system : (ζ,η) , r 
=>  (x,y)  = (ζ,η)
polar coordinates
=>  (r*cos(θ),r*sin(θ)) =(r*cos(u),r*sin(u))

old :(x,y)  = > polar coordinates
-----------------------------------------------

p5a:x=r*cos(θ);
p5b:y=r*sin(θ);

new :(ζ,η)  = > polar coordinates
------------------------------------------------

p5c:ζ=r*cos(u);
p5d:η=r*sin(u);

depends([ζ,η,r,u],t);

differentation , new(ζ,η) = (dζ/dt,dη/dt)
---------------------------------------------------------

p6a:diff(p5c,t);
p6b:diff(p5d,t);

p6c:(p5c)^2+(p5d)^2;
p6d:trigsimp(p6c);
p6e:rhs(p6d) = lhs(p6d);
p6f:solve(p6e,r);
p6f1:p6f[1];
p6f2:p6f[2];
p6f3:diff(p6f2,t);

p6f4:subst([rhs(p6f2)=lhs(p6f2)],p6f3);

p6f5:p6f4*2*r;

calculation : c (new coordinates , when z= 0)
------------------------------------------------------------------

p7:  subst([p6a,p6b],p4c);
p7a:subst([p6f3],p7);
p7b:subst([rhs(p6f5)=lhs(p6f5)],p7a);
p7c:subst([rhs(p6f2)=lhs(p6f2)],p7b);
p7d:subst([p5c,p5d],p7c);
p7e:trigsimp(p7d);

c : per unit mass  c = r*v[u] = r*(r*du/dt)]
-----------------------------------------------------------
analog :  r^2*dθ/dt = r.v[θ] = r*(r*d(θ)/dt)

Kinetic energy :  T : old coordinates./per unit mass
===========================================

p8:v[r] = diff(r,t);

cartesian coordinates (x,y)
---------------------------------

p9a:T = 1/2*(('diff(x,t))^2+('diff(y,t))^2);

polar coordinates   (r,θ)
---------------------------

p9b: v[θ] = r*'diff(θ,t);

p9c:lhs(p9a)=1/2*((rhs(p8))^2+(rhs(p9b))^2);

new coordinates      (ζ,η)
-------------------------------------

p9d:T = 1/2*(('diff(ζ,t))^2+('diff(η,t))^2);

polar coordinates   (r,u)
---------------------------

p9e: v[u] = r*'diff(u,t);
p9f:lhs(p9a)=1/2*((rhs(p8))^2+(rhs(p9e))^2);

use relation,  ' p5', orbit equation.

p10:p5;

R : vector  , R[x,y]  => R[ζ,η]

vector component : ODE's '(ζ,η)'
-----------------------------------------

p10a:subst([R=ζ],p10);
p10b:subst([R=η],p10);

from p6f5 , to do modification

p10c:p6f5;
p10c1:rhs(p10c)=lhs(p10c);

find coef 'ζ','η'
----------------------

p10d:ratcoeff(rhs(p10c),ζ);
p10e:ratcoeff(rhs(p10c),η);

multiply : 'p10a' by 'p10d'
                  'p10b' by 'p10e'

p11:p10a*p10d;
p11a:p10b*p10e;

sum(p11+p11a)

p11b:p11+p11a;

use : 'p10c1'

p11c:lratsubst([p10c1],p11b);

we will rewrite : second differential from 'ζ','η'

p11d1:('diff(ζ,t,1))=f1;
p11d11:('diff(ζ,t,2)) = 'diff(rhs(p11d1),t);
p11d2:('diff(η,t,1))=f2;
p11d22:('diff(η,t,2)) = 'diff(rhs(p11d2),t);

use this for 'p11d..'

p11e:subst([p11d1,p11d11,p11d2,p11d22],p11c );

let(f2*2*('diff(f2,t,1)),'diff(f2^2,t));
let(f1*2*('diff(f1,t,1)),'diff(f1^2,t));

p11f:letsimp(p11e)*dt;

maxima  :code will result into integration with respect 't'
                  =======================================

arg :  'r' ,  'rhs'
---------------------

p11g1: ratcoeff(rhs(p11f),dt*('diff(r,t,1)));

p11g2:integrate(p11g1,r);

result : kinetic energy
---------------------------------

p12:f1^2+f2^2 =p11g2+ constant;

p12a:p12/2;
p12b:lhs(p12a)=expand(rhs(p12a));

p12c1:subst([rhs(p11d1)=lhs(p11d1)],p12b);
p12c2:subst([rhs(p11d2)=lhs(p11d2)],p12c1);

p12c3:subst([rhs(p9d)=lhs(p9d)],p12c2);

rem : found out that   T+V = constant

p12c4:lhs(p12c3) = part(rhs(p12c3),1)+E;

p12c5:p12c4 - part(rhs(p12c3),1);

E =energy of orbit
==================
T = kinetic energy
V = potential energy
==================

p12c6: -1*part(rhs(p12c3),1) = V;

total result : Energy equation .
=========================

p13:subst([p12c6],p12c5);

polar coordinates :
-----------------------------

p13a:subst([rhs(p12c6)=lhs(p12c6)],p13) $
p13b:subst([p9f],p13a);

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