kill(all)$
load(lrats)$

Central Force Motion ( = Force along distance between two masses ),version 1.0
peter.vlasschaert@gmail.com,09/08/2018.

total energy for orbit, μ = reduced mass.
==================================

p1:E=1/2*μ*(('diff(r,t))^2+r^2*('diff(θ,t))^2) +U(r);

from definition l =μ*r^2*'diff(θ,t).
==========================

p2:l =μ*r^2*'diff(θ,t);

p2a:solve(p2,('diff(θ,t,1)));

use : definition of l into 'p1' = E
==========================

p3:subst(p2a,p1);

p4:expand(p3);

We want velocity = dr/dt  from 'p4' = E 
===============================

Equation of Motion.
================

part 1: dr/dt ?
--------------------

p5:solve(p4,('diff(r,t,1)));

p5a:p5[1];
p5b:p5[2];

part  2:θ = θ (r )  ?
-------------------------

rem :  polar equation velocity :  v^2 = r^2*('diff(θ,t,1))^2+('diff(r,t,1))^2
------------------------------------------------------------------------------------------------------
=>   θ = θ (r )   in general : Central Force Motion.
------------------------------------------------------------------------------------------------------

dθ = (dθ/dt)*(dt/dr)*dr 
rem : r , variable for integration for  θ = θ (r )   =  'dr'
dθ = (dθ/dt)/(dr/dt)*dr 
rem : r , variable for integration for  θ = θ (r )   =  'dr'
dθ/dr =(dθ/dt)/(dr/dt)

p6:'diff(θ,r) = 'diff(θ,t)/'diff(r,t);

the use , definition l ,'see above'
===========================

p6a:subst(p2a,p6);

p6b:ode2(p6a,θ,r);

assume(μ>0);

p6c:subst([p5b],p6b);

let :E-U(r)=m, ' can't use subst => lratsubst '
-----------------------------------------------------------------

p6d:lratsubst([E-U(r)=m],p6c);

rem : V= U(r)  = (m1*m2*G)/r = const/r= k/r
---------------------------------------------------------------
values : m1,m2,G = gravitational constant
---------------------------------------------------------------

p6e:subst([p6a],p6c);

p6f:p6e-%c;

part 3: U(r) = k/r
-----------------------

p6g:U(r) = k/r;

assume(k>0);

p6h:subst([p6g],p6f);

menu : calculus => change variable :
===============================
Integral /sum :p6c
old   variable : r
new variable : u
equation         : u=1/r
===============================

p6i:changevar(p6h, u=1/r, u, r);

p6j:expand(p6i);

p6k:subst([abs(u)=u],p6j);

rem : integrate : 1/(a*x^2+b*x+c),x.
============================
a: l^2 >0      =>  a < 0.
b: k>0,μ>0  =>  b < 0
c: sign (E)   =>  c>0,c=0,c<0
============================

p6l:p6k;

find : parts from integral  'p6l';

p6m:lhs(p6l);
p6m1:rhs(p6l);

p6m1a:part(p6m1,0);
p6m1b:part(p6m1,1,1);

parts : integral (use command :part)
==============================

p6m2:part(p6m1,1,2,0);

p6m3:part(p6m1,1,2,1);

p6m4:part(p6m1,1,2,1,2,1);

s: integrate : 1/(a*x^2+b*x+c),x = -1/sqrt(-a)*asin((2*a*x+b)/sqrt(b*b-4*a*c));
=============================================================
from : table 'from books'
=============================================================
=>  a <0, b^2> = 4a*c :'real'  otherwise 'complex'
=============================================================

find : a,b,c 

a ?
==

p7: a=coeff(p6m4,u^2);

b ?
==

p7a:b=coeff(p6m4,u);

c ?
==

p7b:c=ev(p6m4,u=0);

p8:lhs(p6l)=-1/sqrt(-a)*asin((2*a*x+b)/sqrt(b*b-4*a*c));

p8a:[p7,p7a,p7b,x=-u];

rem : abs(u) = u , abs(u) =-u,why use :  'p6k'

p8b:subst(p8a,p8);

p8c: lhs(p8b) =rhs(p8b)*(-l);

p8d:subst([abs(l)=l],p8c);

part 4: find  %c
-----------------------

p9:map(sin,p8d);

rem :  % c = -π/2

p9a:sin(θ+π/2);

p9b:subst([%c=-π/2],p9);

subst : back u=1/r
---------------------------

p9c:subst([u=1/r],p9b);

minimum value of  r , when θ = 0
------------------------------------------------

p9b1:part(rhs(p9c),1)/(2*k*μ);
p9b2:part(rhs(p9c),2)/(2*k*μ);

p9c1:expand(p9b1);
p9c2:factor(p9b2);

p10:p9b2=sqrt(1+w);

assume(k^2*μ+2*E*l^2>0);

p10a:solve(p10,w);

planetary motion,orbit motion ......( kepler problem)
---------------------------------------------------------------------------

general polar equation. ( ellipse,....)
=============================

p11:r=p/(1+ε*cos(θ));

from : 'p9c'
----------------

p11a:lhs(p9c) =subst(p10a,rhs(p9c1/p10));

from: 'p11a'
----------------

p11b:part(p11a,2,1);
p11c: part(p11a,2,2);

p11d:p11a*p11c;

p11e:subst([p11c=ε],p11d);

p11f:p11e+1;

solve : r

p11g:factor(solve(p11f,r))[1];

find p ? : from polar equation .(see above 'p11')
=======================================

p11h: rhs(p11) = rhs(p11g);

final result :kepler problem.
=======================

p11h1:p11;

p11i:solve(p11h,p)[1];

p11j:ε = subst(p10a,rhs(p10));

p11k:p11a;

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