kill(all)$
load(lrats)$

kepler : second and thrid laws  ' version 1.0'
peter.vlasschaert@gmail.com,17/08/2018.

orbital area : swept out by radius vector = (r) , moves trough small angle = (dθ)
rem : point C and B are very close .

def orbital area:
 ΔΑ = 1/2*s*r , s = r*dθ   (1)

Show : areal velocity of orbit  : constant  ?
=============================================
=>  ' second law ' : (kepler)
=============================================
the radius vector from the sun to a planet 'sweeps out'
equal areas in equal time intervals.
=============================================

p1:' diff(A,t) = 1/2*r^2*'diff(θ,t);

from : 
=================================================================
 new: 10/06/2018:motion stays in one plane (angular momentum) 'version 1.0'
 =>  angular momentum = constant  = h = r x v = > h = r*v[θ], angle(r,v) = π/2
=================================================================

p2: h = (r)*(r*'diff(θ,t));

p3:subst([rhs(p2)=lhs(p2)],p1);

period   T : elliptical orbit
==============================
T = area/arealvelocity 
-------------------------------------------------------
rem  :  area = total area of full ellipse.
-------------------------------------------------------

p4: T= A/(rhs(p3));

rem : total area : circle area = π*r^2 ,ellipse area = π*a*b

p5: A = π*a.b;
p6: subst([p5],p4);

rem : b = a*sqrt(1-e^2)
ellipse : e > 1
-------------------------------------------------
rem : follows geometry of ellipse.
--------------------------------------------------
rem  : e = c/a, a^2 = b^2+c^2
--------------------------------------------------

assume ( μ>0,a>0,b>0,e<1);

p6a : b = a*sqrt(1-e^2);

p6b :subst([p6a],p6);

polar equation : conics
====================

p7: r= p/(1+e*cos(θ));
p7a : p = h^2/μ;
p7b : p= a*(1-e^2);

rem : μ = G *( M(sun)+m(planet)) ,solar system
           μ = G.( M(earth ) +m(object)) , orbit around the earth.

p8:rhs(p7a) = rhs(p7b);

p8a:solve(p8,h);

p8b:(p8a)[2];

we want remove :  h  from 'p6b'.
==========================

Find : thrid law of kepler : ?
=============================

T^2/a^3 =constant , see below  'μ'

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

use : ' p6 b' and 'p8b'

p9:subst([p8b],p6b);

assume ( μ>0,a>0,e<1);

p9a:fullratsimp(p9);

p9b:(p9a)^2;

p9c:lhs(p9b)=factor(rhs(p9b));

p9d:lhs(p9c) = 4*%pi^2/μ*a^3;

rem  T^2/a^3, is not constant for all planets. (see : μ ' above')

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