kill(all)$
load(lrats)$

Central Force of  motion (=diff eq),version 1.0a
peter.vlasschaert@gmail.com, 09/08/2018.

rem : find F(r)  ?, orbit is known.
------------------------------------------------

from Lagrangain  equation :  r-coordinate.
===================================

p1:μ*('diff(r,t,2)-r*('diff(θ,t))^2) = F(r);

build u : diff eq
============= .

p2:u=1/r;

part 1:du/dθ : ?
-------------------

u(θ) => u(r)*r(θ).r(t).θ(t)

1e)

du/dθ =du/dr*dr/dθ = du/dr*dr/dθ = du/dr*dr/dt*dt/dθ.

p3:'diff(u,θ) = 'diff(u,r)*'diff(r,t)/'diff(θ,t);

p3a: 'diff(u,r) = subst([p2],'diff(u,r));

p3b:lhs(p3a) = ev(rhs(p3a),diff);

p3c:subst([p3b],p3);

2e)

use definition : l
------------------------

find : dθ/dt

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

p4a:solve(p4,('diff(θ,t,1)));

insert 'p4a' into 'p3c'

p4b:lhs(p3c)  = subst(p4a,rhs(p3c));

part 2:diff(θ,u,2) : ?
------------------------

p5:'diff(u,θ,2) = 'diff(rhs(p4b),θ);

rem : d/dθ  = d/dt*dt/dθ = d/dt (1/dθ/dt)
===============================

p5a:lhs(p5) = 1/'diff(θ,t)*'diff(rhs(p4b),t);

use definition : l , use 'p4a'

p5b:lhs(p5a)=subst(p4a,rhs(p5a));

ev : diff

depends(r,t);

depends must be used otherwise : rhs(p5b)  =  0
-------------------------------------------------------------------------

p5c:lhs(p5b)=ev(rhs(p5b),diff);

p5d:solve(p5c,('diff(r,t,2)));

part 3:  ' use above ' fill in  'p1' :?
----------------------------------------------

1e) fill in:
=======

p6:subst(p4a,p1)$

p6a:subst(p5d,p6);

p6b:expand(p6a);

2e) find coef : 'diff(u,θ,2), lhs
=======================

p6c:ratcoeff(lhs(p6b),'diff(u,θ,2));

p6c1:lratsubst([rhs(p2)=lhs(p2)],p6c);

3e) find coef : without 'diff(u,θ,2),lhs
==============================

p6d:ev(lhs(p6b),('diff(u,θ,2))=0);

p6d1:lratsubst([rhs(p2)=lhs(p2)],p6d);

4e)  into  lhs(p1)
==============

p6e:p6c1*'diff(u,θ,2)+p6d1 = rhs(p1);

p6f:1/p6c1*p6e;

p6g:ratsimp(p6f);

part  4:  solve 'p6g'
-----------------------------

1e)     F(r) = -div(V) = -k/r^2
====================
V = potential = V(r), divergence = div
------------------------------------------------------

p7:subst([F(r)=-k*u^2],p6g);

p7a:ratsimp(p7);

2e)   ode2, u
==========

1e) u = u(θ)
=========

p7b:ode2(p7a,u,θ);

p7c:subst([p2],p7b);

rem :  %k1*sin(θ)+%k2*cos(θ) = d*cos(θ-θo);
d,θo = find out  ? 

take axis so that : θo =0
------------------------------------

p7b:lhs(p7c)=rhs(p7a)+d*cos(θ);

p7c:solve(p7b,r)[1];

p7c1:part(rhs(p7c),1);
p7c2:part(rhs(p7c),2);

p7d1: p7c1/(k*μ);
p7d2:p7c2/(k*μ);

p7e1:(p7d1);
p7e2:expand(p7d2);

p7f:lhs(p7c)= p7e1/p7e2;

p = semi-latus rectum = 'p7e1'

p7g:lhs(p7f) = p/p7e2;

ε = (d*l^2)/(k*μ)

p7h:ε = (d*l^2)/(k*μ);

p7i: lratsubst([rhs(p7h)=lhs(p7h)],p7g);

rem : polar equation of conic section( ellipse,...) = 'p7i'

rem :  find d ? uniquely determines , given E, angular momentum.
--------------------------------------------------------------------------------------------------

We can use : 'p6g'  when  r=r(θ), known orbit.
use :  u = 1/r    
=>  Find force law.

Example :  r= r(θ)
=>   r=k*exp(θ)  
Find : F(r)

start :
=====

p8:p6g;

p8a:subst([p2],p8);

use :  equation of orbit. (into p8a) :
=============================

p8b:subst([r=k*exp(θ)],p8a);

p8b1:k*exp(θ) = r;

p8c:ev(p8b,diff);

p8d:solve(p8c,F(k*%e^θ))[1];

Force Law : F(r)
=============

p8e:subst([p8b1],lhs(p8d)) = rhs(p8d);

rem : removes θ,from rhs 'p8e'
---------------------------------------------

p8f:solve(p8b1,θ);

use  ' p8f ' into 'p8e'
================

p8g:subst(p8f,p8e);

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