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

flight over rotating spherical earth (part3:coordinate system ),version 1.0
peter.vlasschaert@gmail.com, 14/09/2018.

coordinate system 

detail vehicle point P 

spherical coordinates : (r,θ,φ)  ,Φ= φ
===============================
θ  : longitude
r   : radius
Φ : latitude
===============================

right triangle equations : see figure 3
===============================
point M :coordinates (x,y,z)
x = OP1
y = OP2
z = OP3

1e) angle : θ
------------------

p1:cos(θ) = x/rr;
p2:cos(π/2-θ) =y/rr;

find x,y ?

p1a:solve(p1,x);
p2a:solve(p2,y);

2e) angle : Φ
-------------------

p3:cos(φ) =rr/r;
p4:cos(π/2-φ) = z/r;

find z ?

p4a:solve(p4,z);

we need eliminate rr ?

p3a:solve(p3,rr);

cartesian coordinates : x,y,z  as function r,θ,φ
=======================================

spherical coordinates (r,θ,φ)
========================

p5:subst(p3a,p1a[1]);
p5a:subst(p3a,p2a[1]);
p5b:p4a[1];

spherical coordinates (V,γ,ψ)
=========================

rem : see fig 1 also , system : Oxyz
rem : P to cc              ( direction y axis)
rem : P along  vec r ( direction x axis)

because : angle between(purple line = -------- ,V) = γ = gamma

projection : V on the x-axis
=======================

p6: cos(π/2-γ) = x/V;

p6a:solve(p6,x)[1];

rem : from above  z => x  same analog y,z see above.

projection : V on the y-axis : analog from the above
=======================

p7: y = a1*V*cos(γ);

projection : V on the z-axis : analog from the above
=======================

p7a: z= a2*V*cos(γ);

How to find : a1,a2  
================
=>  V  below  local horizon 
=>  Definition ψ : angle ()
result  :  rotation of angle ψ  from y  to  local horizon . ( a1=cos(ψ)       )
                                                                                                     ( a2=cos(π/2-ψ))

result :  spherical coordinates 
=========================

p8:[a1=cos(ψ),a2=cos(π/2-ψ)];
p8a:p6a;
p8b:subst(p8,p7);
p8c:subst(p8,p7a);

set  :vec V = VV
------------------------

p8d:VV =rhs(p8a)*ex+rhs(p8b)*ey+rhs(p8c)*ez;

vec ω = ωω   , vector direction  of the north pole.
=========================================

from fig 1 : ωω ,components are in the (xz plane).

 constant ω : rotation of the earth

projection ω vector  to  x
-------------------------------------

p9: x=ω*cos(π/2-Φ);

projection ω vector  to  z
-------------------------------------

p9a:z=ω*cos(Φ);

vector ω = ωω
----------------------

p9b:ωω = x*ex+y*ey+z*ez;
p9c:subst([p9,p9a,y=0],p9b);

cross product = ωω x  VV = cross1
=============================

p10: ωω =[rhs(p9),0,rhs(p9a)];
p10a: VV = [rhs(p8a),rhs(p8b),rhs(p8c)];

p10b:cross1 = express(rhs(p10) ~rhs(p10a));

rem : use load(vect) , vector product use the operator = ~

components of  cross 1:
------------------------------------

p10b1:x=rhs(p10b)[1];
p10b2:y=rhs(p10b)[2];
p10b3:z=rhs(p10b)[3];

vector notation cross 1:
------------------------------------

p10b4:lhs(p10b) = x*ex+y*ey+z*ez;
p10b5:subst([p10b1,p10b2,p10b3],p10b4);

cross product = ωω x (ωω x rr) =cross3
===============================================

1e)  vector r = rr,reason axis from center of earth aligned is along the  x direction.
-------------------------------------------------------------------------------------------------------------------------

p11:rr=r*ex;

vector notation : rr

p11a:rr =[r,0,0];

2e)  cross2  =  ωω x r
---------------------------------

p11b:cross2= express(rhs(p10)~rhs(p11a));

3e)  cross3  =  ωω x cross2
-----------------------------------------

p11c:cross3=express(rhs(p10)~rhs(p11b));

final result : cross3
-----------------------------

components of  cross 3:
------------------------------------

p11c1:x=rhs(p11c)[1];
p11c2:y=rhs(p11c)[2];
p11c3:z=rhs(p11c)[3];

vector notation cross 3:
------------------------------------

p11c4:lhs(p11c) = x*ex+y*ey+z*ez;
p11c5:subst([p11c1,p11c2,p11c3],p11c4);

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