f=(1+Log[#])/(1+e Cos[#])/2& (* Sztuczna trajektoria, na potrzeby rysunku *)

ϕ0=7/2

e=7/10

{xP,yP}={f[ϕ0]Cos[ϕ0],f[ϕ0]Sin[ϕ0]}

elipsa1=PolarPlot[f[ϕ],{ϕ,0,2Pi},Epilog->{
Disk[{xP,yP},1/30],
Disk[{0,0},1/10],
Arrow[{{0,0},{xP,yP}}],
Line[{{0,0},{1/(1+e),0}}],
Circle[{0,0},1/4,{0,ϕ0}],
Arrow[{{xP,yP},{xP,yP}1.2}],
Arrow[{{xP,yP},{xP,yP}-{yP,-xP}0.2}],
Text["M",{0,-1/5},BaseStyle->{FontSize->24}],
Text["m",{xP,yP},{0,-3/2},BaseStyle->{FontSize->24}],
Text["ϕ(t)",{0,1/3},BaseStyle->{FontSize->24}],
Text["r(t)",{xP/2,yP/2},{0,1},BaseStyle->{FontSize->24}],
Text["(Subscript[(OverHat[e]), r])",{xP,yP}1.2,{1,1},BaseStyle->{FontSize->24}],
Text["(Subscript[(OverHat[e]), ϕ])",{xP,yP}-{yP,-xP}0.2,{1,1},BaseStyle->{FontSize->24}]
},Axes->None,PlotRange->{{-4,2},{-2,2}},ImageSize->600]

f=1/(1+e Cos[#])&

p=1

a=(1/(1+e)+1/(1-e))/2

b=Integrate[f[t]^2/2,{t,0,2Pi}]/Pi/a

{xP,yP}={f[ϕ0]Cos[ϕ0],f[ϕ0]Sin[ϕ0]}

elipsa2=PolarPlot[f[ϕ],{ϕ,0,2Pi},PlotStyle->Thick,Epilog->{
Disk[{xP,yP},1/60],
Disk[{0,0},1/20],
Disk[{-e a,0},1/100] (* srodek elipsy, center *),
Disk[{-2 e a,0},1/100] (*drugie ognisko, 2nd focus *),
Arrow[{{0,0},{xP,yP}}],
Line[{{0,0},{0,1}}],
Line[{{0,0},{1/(1+e),0}}],
Line[{{1/(1+e),0},{1/(1+e),3/2b}}],
Line[{{-1/(1-e),0},{-1/(1-e),3/2b}}],
{Arrowheads[{-0.03,0.03}],Arrow[{{-1/(1-e),3/2b},{1/(1+e),3/2b}}]},
Text["2a",{-e a,3/2b},{0,-1},BaseStyle->{FontSize->24}], 
Line[{{-2a,b},{-a e,b}}],
Line[{{-2a,-b},{-e a,-b}}],
{Arrowheads[{-0.03,0.03}],Arrow[{{-2a,b},{-2a ,-b}}]},
Text["2b",{-2 a,0},{-1,0},BaseStyle->{FontSize->24}], 
Circle[{0,0},1/4,{0,ϕ0}],
{Dotted,Line[{{0,0},{xP,-yP},{-2e a,0}}]},
{Arrowheads[{-0.03,0.03}],Arrow[{{-2e a,0},{-e a ,0}}]},
Text["e a",{-3/2e a,0},{0,1},BaseStyle->{FontSize->24}],
Text["p=a(1-(e^2))",{0,1/2},{1,0},BaseStyle->{FontSize->24}],  
Text["M",{0,-1/5},BaseStyle->{FontSize->24}],
Text["m",{xP,yP},{0,-3/2},BaseStyle->{FontSize->24}],
Text["ϕ",{-1/4,1/4},BaseStyle->{FontSize->24}],
Text["r",{xP/2,yP/2},{0,1},BaseStyle->{FontSize->24}]},Axes->None,PlotRange->{{-2a,a},{-2b,2b}},ImageSize->600]

Eliptyczność wszystkich orbit zamkniętych w polu 1/r jest znanym m.in. z mechaniki klasycznej faktem. Poniżej pokażemy, jak powiązać geometryczne parametry elipsy z zachowanymi wielkościami fizycznymi.  Wychodzimy od lagrangianu cząstki we współrzędnych biegunowych.

L=1/2 m (r'[t]^2+r[t]^2 ϕ'[t]^2)+G M m/r[t]

D[L,r'[t]]

D[L,ϕ'[t]]

D[L,r[t]]

D[L,ϕ[t]]

rown1=D[D[L,r'[t]],t]-D[L,r[t]] (* Pierwsze równanie Lagrange-Eulera, na r(t) *)

rown2=D[L,ϕ'[t]]==J  (* Drugie równanie Lagrange-Eulera, na ϕ(t) *)

Energia=r'[t]D[L,r'[t]]+ϕ'[t]D[L,ϕ'[t]]-L//Simplify//Expand (* Dalej skorzystamy z zachowanej energii *)

Energia2 = Energia/.Solve[rown2,ϕ'[t]][[1]] (* Eliminujemy ϕ(t) kortzystając z II równania Lagrange'a-Eulera *)

Energia3 = Energia2/.{r->(1/u[#]&)} (* Znany trik dla orbit, upraszczający równania,  polega na zastąpieniu r przez odwrotność *)

Energia4 = Energia3 /.{u->(U[ϕ[#]]&)} (* Wprowadzamy nową funkcję, opisującą zależność u(ϕ). Na tablicy/papierze zwykle nadal używa sie na nią tej samej literki u, ale papier przyjmie wszystko. Mathematica zmusza nas do zrozumienia, że jest to inna funkcja, którą oznaczamy wielką literą U. *)

Energia5 = Energia4 /.Solve[rown2,ϕ'[t]][[1]] (* Ponownie korzystamy z II r. L-E *)

Energia6=Energia5/.{ϕ[t]->ϕ,r[t]->1/U[ϕ]} (* Konwertujemy równanie do postaci nie zawierającej czasu t *)

D[Energia6,ϕ]//Simplify (* Warto zwrócić uwagę nma fakt, iż pochodna powyższego równania to zwykły oscylator harmoniczny *)

oscylator=m %/J^2/.{U'[ϕ]->1}//Expand (* Oscylator harmoniczny dla sprężyny o dł. spoczynkowej G (m^2) M/(J^2) *)

DSolve[oscylator==0,U,ϕ]

Energia6/.U->(A Cos[#]+G M m^2/J^2&) (* Podstawiamy powyższe z powrotem do równania toru *)

Simplify[%]

Solve[%==E0,A]

Solve[Energia6==E0,U'[ϕ]]

fpodcalkowa = U'[ϕ]^2/.%[[-1]]//Expand

p=J^2/(G M m^2)

e=(Sqrt[1+((2 E0 (J^2))/((G^2) (m^3) (M^2)))])

fpodcalkowa -(((2 E0 m)/(J^2))+(((G^2) (m^4) (M^2))/(J^4))-(U[ϕ]-G m^2 M/J^2)^2)//Simplify

fpodcalkowa - (e^2/p^2-(U[ϕ]-1/p)^2)//Simplify

Fundamentalne wzory na wielką półoś a oraz mimośród elipsy e wyrażone poprzez fizyczne parametry układu: energia E0 oraz moment pędu J.  Dla układów związanych E0<0.

e=FullSimplify[J^2Sqrt[((2 E0 m)/(J^2))+(((G^2) (m^4) (M^2))/(J^4))]/(G M m^2),J>0&&G>0&&m>0&&M>0]^2//Simplify//Sqrt

a=1/2(p/(1+e)+p/(1-e))//FullSimplify

1+((2 E0 (J^2))/((G^2) (m^3) (M^2)))/.{J->m v R,E0->-G M m /R+1/2 m v^2} (* Test ,,sensowności'' wzoru na e w przypadku orbity kołowej - ma wyjść e=0 *)

FullSimplify[%]

%/.v->Sqrt[G M /R]