Zabawa z tzw. problemem Keplera dotyczącym ruchu w polu siły centralnej

Zaczynamy od parametrów takich, jak dla ruchu Ziemi wokół Słońca po orbicie kołowej po to, by określić iloczyn kappa=G*M*m w rozsądnych jednostkach
G - stała grawitacji
M - masa Słońca
m - masa Ziemi 
Korzystamy z warunku że siła dośrodkowa równa się sile grawitacyjnej
    Fg= G*M*m/r^2 = kappa/r^2
    Fd = m*omega^2*r
czyli
    kappa/r^2= m*omega^2*r,
wiec
     kappa= m*omega^2*r^3, 
 gdzie 
      omega = 2*pi/T 
  to prędkość kątowa.
  T to okres obiegu Ziemi dookoła Słońca, a  dist to odległość Ziemia-Słońce.

dist=149.597870691; (* w milionach kilometrow *)
omega=(2.*Pi)/365.0; (* rad/rok*)
kappa=omega*omega*dist*dist*dist;
m=1; (* w masach Ziemi *)

Parametry krzywej stożkowej wyrażone przez energię całkowitą (Etot) oraz całkowity moment pędu (Ltot).  Rozwiązaniem równań ruchu jest krzywa
   r(phi) = p/(1+e*cos(phi),
gdzie
    p = L^2/(m*kappa)
    e  = sqrt (1 + 2*E*L^2/(m*kappa^2))

Ltot=m*dist*omega*dist;

Etot=-kappa/dist +m*(dist*omega)^2/2;

p=Ltot^2/(m*kappa) ;

e=Sqrt[1+2*Etot*Ltot^2/(m*kappa^2)];

rho[fi_]:=p/(1+e*Cos[fi])

Ltot

Etot

p

e

Dla znikającego mimośrodu e mamy okrąg. Przy ustalonym Ltot, Etot dla okręgu jest minimalną dopuszczalną wartością. Inne wartości Etot mogą być tylko większe !   
Ta minimalna wartosc energii moze byc znaleziona ze wzoru
e=0 -> 1+2*Etot*Ltot^2/(m*kappa^2)=0 -> Emin = m*kappa^2/(2*Ltot)
(Etot < 0)

PolarPlot[rho[fi],{fi,0,2*Pi}, AxesLabel->{"x","y"}]

Etot=.;e=.;

Teraz weźmy ujemną wartość Etot, ale dwukrotnie bliższą zera

Etot=(-kappa/dist +m*(dist*omega)^2/2)/2

e=Sqrt[1+2*Etot*Ltot^2/(m*kappa^2)]

PolarPlot[rho[fi],{fi,0,2*Pi}, AxesLabel->{"x","y"}]

Etot=.;e=.;

Weźmy ujemną wartość Etot, ale jeszcze bardziej bliską zera

Etot=(-kappa/dist +m*(dist*omega)^2/2)/3

e=Sqrt[1+2*Etot*Ltot^2/(m*kappa^2)]

PolarPlot[rho[fi],{fi,0,2*Pi}, AxesLabel->{"x","y"}]

Etot=.;e=.;

Weźmy Etot dokładnie równą zero

Etot=0

e=Sqrt[1+2*Etot*Ltot^2/(m*kappa^2)]

PolarPlot[rho[fi],{fi,-Pi/1.5,Pi/1.5}, AxesLabel->{"x","y"}]

Etot=.;e=.;

Weźmy małą dodatnią wartość Etot

Etot=1

e=Sqrt[1+2*Etot*Ltot^2/(m*kappa^2)]

PolarPlot[rho[fi],{fi,-Pi/1.5,Pi/1.5}, AxesLabel->{"x","y"}]

Etot=.;e=.;

Weźmy nieco większą wartość Etot

Etot=2

e=Sqrt[1+2*Etot*Ltot^2/(m*kappa^2)]

PolarPlot[rho[fi],{fi,-Pi/1.5,Pi/1.5}, AxesLabel->{"x","y"}]

Etot=.;e=.;

... i jeszcze większą wartość Etot

Etot=20

e=Sqrt[1+2*Etot*Ltot^2/(m*kappa^2)]

PolarPlot[rho[fi],{fi,-Pi/1.7,Pi/1.7}, AxesLabel->{"x","y"}]