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

equations of motion spherical coordinates,'version 1.0'
peter.vlasschaert@gmail.com,15/09/2018

1e) relation between : cartesian coordinates <=> spherical coordinates
===========================================================
(x,y,z)  <=> (r,θ,φ)
r : radius       : range  = 0 ≤ r < infinity
θ : longitude : range  = 0 ≤ θ ≤ 2*π
φ :tatitude     : range  = - π/2 ≤ φ ≤ π/2  : φ = 0 @ equator (=evenaar)
===========================================================

p1   : x=r*cos(φ)*cos(θ);
p1a : y=r*cos(φ)*sin(θ);
p1b : z=r*sin(φ);

2e)  T = Kinetic energy  '(x,y,z)  <=> (r,θ,φ)'
===================================

a)  T = 1/2*m*((dx/dt)^2+(dy/dt)^2+(dz/dt)^2) : start
---------------------------------------------------------------------------

depends([x,y,z],t);

p2: T= 1/2*m*((diff(x,t))^2+(diff(y,t))^2+(diff(z,t))^2);

b) differentation , x,y,z with respect to 't' (use 'p1,p1a,p1b')
-----------------------------------------------------------------------------------------

depends([r,θ,φ],t);

1e) dx/dt,dy/dt,dz/dt
------------------------------

p2a:diff(p1,t);
p2b:diff(p1a,t);
p2c:diff(p1b,t);

2e) (dx/dt,dy/dt,dz/dt )^2
--------------------------------------

p3a:(p2a)^2;
p3b:(p2b)^2;
p3c:(p2c)^2;

(dx/dt)^2
--------------

p3a1:expand(p3a);
p3a2:trigsimp(p3a1);

(dy/dt)^2
--------------

p3b1:expand(p3b);
p3b2:trigsimp(p3b1);

(dz/dt)^2
--------------

p3c1:expand(p3c);
p3c2:trigsimp(p3c1);

result : T(x,y,z)  <=> T(r,θ,φ)
=======================

p4:lhs(p2) =1/2*m*(rhs(p3a2)+rhs(p3b2)+rhs(p3c2));

p4a:lhs(p4)=trigsimp(rhs(p4));

3e)  U = potential energy  '(x,y,z)  <=> (r,θ,φ)'
====================================

central gravitational field : U = U(r)
----------------------------------------------------

depends([U],[r,θ,φ]);

rem : relation between : force and U :  F = - grad(U)

4e )  Lagrange equation = Lag ( 'lagrangian')
=====================================

q        = generalized coordinate    =>  maxima = q
dq/dt = generalized velocities      =>  maxima = qt

= >    d/dt(∂L/qt) -dL/dq  = 0

p5:  lag(L,q,qt) := diff(diff(L,qt),t) -diff(L,q);

spherical coordinates 
-----------------------------------------------
q        = r,θ,φ
qt       = dq/dt = dr/dt,dθ/dt,dφ/dt
------------------------------------------------

a)  L = T-U
----------------

p6:L = rhs(p4a) -U;

b)  Lagrange equation : ' r,θ,φ '
----------------------------------------------

lagrangian : q = 'r' ,dq/dt = dr/dt  (use 'p5')
----------------------------------------------------------------

p7a:lag(rhs(p6),r,diff(r,t)) = 0;
p7b:expand(p7a);
p7c: solve(p7b,diff(U,r))[1];
p7d:(rhs(p7c) = lhs(p7c))*(-1)/m;
p7e:expand(p7d);

lagrangian : q = 'θ' ,dq/dt = dθ/dt  (use 'p5')
----------------------------------------------------------------

p8a:lag(rhs(p6),θ,diff(θ,t)) = 0;
p8b:expand(p8a);
p8c: solve(p8b,diff(U,θ))[1];
p8d:(rhs(p8c) = lhs(p8c))*(-1)/m;
p8e:expand(p8d);

lagrangian : q = 'φ' ,dq/dt = dφ/dt  (use 'p5')
----------------------------------------------------------------

p9a:lag(rhs(p6),φ,diff(φ,t)) = 0;
p9b:expand(p9a);
p9c: solve(p9b,diff(U,φ))[1];
p9d:(rhs(p9c) = lhs(p9c))*(-1)/m;
p9e:expand(p9d);

result : three equations of motion
============================

p10:p7e;
p11:p8e;
p12:p9e;

some simplifications :   'p11,p12'
-------------------------------------------------

for   equation   'p12'
================

p13 :'diff((r^2*diff(φ,t)),t)=diff((r^2*diff(φ,t)),t);

so final equation ' p12'.
===================

p13a:part(lhs(p12),1);
p13b:part(lhs(p12),2);
p13c:part(lhs(p12),3);

p13d:lhs(p13)+p13c = rhs(p12);

for   equation   'p11'
================

p14:'diff((r^2*(cos(φ)^2)*('diff(θ,t))),t) = diff((r^2*(cos(φ)^2)*('diff(θ,t))),t);

rem : check  , is ,' command maxima ', result : true or false

p14a:rhs(p14) = lhs(p8e);

p14b:is(p14a);

p14c:lhs(p14) = rhs(p8e);

final solution : 'simplify'
====================

p15:p7e;
p16:p13d;
p17:p14c;

central gravitational field ; conservative field  : U = U(r)
=============================================

p18:['diff(U,φ,1)=0,'diff(U,θ,1)=0]$

p19a:subst(p18,p15);
p19b:subst(p18,p16);
p19c:subst(p18,p17);

rem : kepler first law from lagrange equations.(above)
           =>  φ = 0  so dφ/dt = 0

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

-  coordinates  :  'p1,p1a,p1b'
   ------------------------------------------

p20:[φ = 0];
p20a:subst(p20,p1);
p20b:subst(p20,p1a);
p20c:subst(p20,p1b);

rem : polar coordinates : definition

- equations motion polar coordinates.
   ------------------------------------------------------

p21a:subst(p20,p19a);
p21b:subst(p20,p19b);
p21c:subst(p20,p19c);

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

rem : p21c

p22:ode2(p21c,θ,t);

rem : df/dt = 0  follows f = constant

p22a : part(p21c,1,1) = constant ;

rem : r^2*('diff(θ,t,1)) = r*(r*dθ/dt)= r*v[θ].
------------------------------------------------------------------------------
rem : v = dr/dt *e[t] + r*dθ/dt*e[n] = v[r]*e[r]+v[θ]*e[n]
t  : tangential component
n : normal      component
e[] : unit vector
------------------------------------------------------------------------------

rem : p22a 'first Law of Kepler'

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