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

flight over rotating spherical earth(part2:orbitplaneunit),version 1.0
peter.vlasschaert@gmail.com,13/09/2018

orbital plane (orthogonal set  P,Q,W )

definition : rotation from part 1
=========================

definition : x 'rotation with angle : α '
-----------------------------------------------------

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

definition : y 'rotation with angle : α '
------------------------------------------------------

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

definition : z 'rotation with angle : α '
-----------------------------------------------------

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

find: orthogonal set  ? {P,Q,W}
==========================

 first rotation : ΩO    (about z-axes)
=============================

p4a:rotz(ΩO);

second rotation : iO      (about x-axes)
================================

p4b:rotx(iO);

third rotation : ω    (about z-axes)
============================

p4c:rotz(ωO);

rem : relation 'zyz' , (are eulerian angles)

use :stacked principal for rotation .

p5:tot[rotation] = p4c.p4b.p4a;

modification  needed because :  =>   Q xP = W  => angle : + iO (see figure1)
                                                                      from figure follows 
                                                                      PxQ = -W  => angle : - iO
see this : angle ( = given from equator to the orbital plane),why we need to use negative angle.

p5a:rotx(-iO);

p5b:corrected1l[rotation] = p4c.p5a.p4a;

next modification needed because :  P : negative  x direction  => (see figure1)
                                                                      Q : negative  y direction => (see figure1)
----------------------------------------------------------------------------------------------------------------------
see this : when orbital plane have no inclination ( = iO) :result =>  : ωO  => - ωO
                                                                                                                    =>  : ΩO  =>  - ΩO
----------------------------------------------------------------------------------------------------------------------

p6a:rotz(-ΩO);
p6b:rotz(-ωO);

p6c:corrected2[rotation] = p6b.p5a.p6a;

  relation between both orthogonal sets : {P,Q,W} and {e1,e2,e3}.
--------------------------------------------------------------------------------------------------
relation : cross product =>   'QxP = W', 'e1 x e2 = e3'
--------------------------------------------------------------------------------------------------

p6d:rhs(p6c);

find term from ' matrix' ,'p6d'
------------------------------------------

p6e:p6d[3,3];

relation between two sets
--------------------------------------

p7a:P=p6d[1,1]*e1+p6d[1,2]*e2+p6d[1,3]*e3;
p7b:Q=p6d[2,1]*e1+p6d[2,2]*e2+p6d[2,3]*e3;
p7c:W=p6d[3,1]*e1+p6d[3,2]*e2+p6d[3,3]*e3;

find matrix coef : from 'p7a','p7b','p7c'
----------------------------------------------------------
rem : 'p6d' =(-1)*coefmatrix()
----------------------------------------------------------

p7d:coefmatrix([p7a,p7b,p7c],[e1,e2,e3])*(-1);

relation : vector P = P(px,py,pz)
                  vector Q = Q(qx,qy,z) 
                  vector W = W(wx,wy,wz)

p8a:P=px*e1+py*e2+pz*e3;
p8b:Q=qx*e1+qy*e2+qz*e3;
p8c:W=wx*e1+wy*e2+wz*e3;

p8d:coefmatrix([p8a,p8b,p8c],[e1,e2,e3])*(-1);

result : 'p7d' and 'p8d'
---------------------------------

p8e:p8d=p7d;

see result ,p,q,w , components vectors ,P,Q,W

We find relation between : P=>xw,Q=>yw,W=> zw and e1 =>x,e2 => y,e3 => z
-------------------------------------------------------------------------------------------------------------------

p9:[P = xw,Q = yw, W = zw,e1=x,e2=y,e3=z];

p9a:subst(p9,P=px*e1+py*e2+pz*e3);
p9b:subst(p9,Q=qx*e1+qy*e2+qz*e3);
p9c:subst(p9,W=wx*e1+wy*e2+wz*e3);

rem :coef matrix : 'p8d'

find {x,y,z} => {xw,yw,zw}
----------------------------------------------------------------
Ax=b  => A^(-1)*A*x=A^(-1)*b => x = A^(-1)*b
rotation matrix : A^(-1) = transpose(A) ,orthogonal

p10:[xw,yw,zw];
p101:transpose(p10);
p10a:[x,y,z];
p10b:transpose(p8d);

p10c : transpose(p10a) = p10b.'' p101;

single vector equation  : R = R(x,y,z)
                  vector P = P(px,py,pz)
                  vector Q = Q(qx,qy,z) 
                  vector W = W(wx,wy,wz)         
================================                                  

p10d: R=xw*P+yw*Q+zw*W;

because  zw = 0 ( motion always in (xw,yw) plane,kepler assumption 'two body problem')
===========================================================================

p10e:subst([zw=0],p10d);

rem : R,P,Q are vectors.

relation between : elliptic elements (geometry)
---------------------------------------------------------------------
a = semi-major axis
c = a/e  =>  c =a*e
---------------------------------------------------------------------
what about : vO 
---------------------------------------------------------------------

angles  E, vO

from figure 2  : projction from r to   xw = M'M
                            projection from r to yw = B[1]M'
=======================================

p11:  xw = r*cos(vO);
p11a:yw = r*sin(vO);

use E : angle from circumscribing circle ( Eccentric Anomaly)
===================================================
from figure2 : =>  Oorb  to N      = a*cos(E) .
                          =>  Oorb  to B[1]  = a*e          .
                          =>  B[1]    to  N      = xw          .
===================================================

p12:xw =a*cos(E) -ae;
p12a: r =a*(1-e*cos(E));

from :pythagoras
--------------------------

p12b:r^2=(xw)^2+(yw)^2;

rw,xw,yw : ?
-------------------

p12c:rw=a*(1-e*cos(E));
p12d : rhs(p11) = rhs(p12);
p12e : rhs(p11a) = a*sqrt(1-e^2)*sin(E);

depends([vO,E,r],t);

differentation :  'p12a','p12d','p12e'
----------------------------------------------------

p13:diff(p12a,t);
p13a:diff(p12d,t);
p13b:diff(p12e,t);

relation : d(xw)/dt ,d(yw)/dt.
-----------------------------------------

p14   : 'diff(xw,t) = lhs(p13a);
p14a : 'diff(xw,t) = rhs(p13a);
p15   : 'diff(yw,t) = lhs(p13b);
p15a : 'diff(yw,t) = rhs(p13b);
p15b:p13;

p16:solve(p11,vO);
p16a:solve(p12a,E);

p17:diff(p16[1],t);
p17a:diff(p16a[1],t);

other : relation possible , see the problem, for differentation

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