G=1; (* Stała grawitacyjna, musi być 1 dla kompatybilności z  NBodySimulation *)

m=E; (* Masa 1 ciała *)

M=Pi; (* Masa 2 ciała *)

d=Catalan;

T=2Pi/Sqrt[G(m+M)/d^3]

twoBody =NBodySimulation["InverseSquare",
{
<|"Mass"->m,"Position"->{0,0,0},"Velocity"->{0,   0.5,-1.0}|>,
<|"Mass"->M,"Position"->{d,0,0},"Velocity"->{0,-0.5,   0.0}|>
},8T]

ParametricPlot3D[Evaluate[twoBody[All,"Position",t]],{t,0,4T}]

srodekMasy=twoBody[All,"Mass"].twoBody[All,"Position",t]/Total[twoBody[All,"Mass"]];

ParametricPlot3D[Evaluate[twoBody[All,"Position",t]-srodekMasy],{t,0,4T},PlotLabel->"Ruch ciał względem środka masy"]

ParametricPlot3D[Evaluate[twoBody[2,"Position",t]-twoBody[1,"Position",t]],{t,0,4T},PlotLabel->"Ruch ciała 2 względem ciała 1"]

Standardowa redukcja problemu 2 ciał do ruchu 1 ciała w polu centralnym.

r=twoBody[2,"Position",t]-twoBody[1,"Position",t]; (* Pozycja względna *)

(*r=twoBody[2,"Position",t]-srodekMasy;*)

μ=(1/m+1/M)^-1 ;(* Masa zredukowana *)

p=μ D[r,t];

J=Cross[r,p];

Plot[{J[[1]],J[[2]],J[[3]]},{t,0,6T},PlotStyle->{Red,Green,{Blue,Dashed}},PlotLabel->"Zachowanie składowych momentu pędu"]

A=Cross[p,J]-G  μ^2 (m+M) r/Norm[r];

μ^2 (m+M)==μ m M (* Sprawdzenie, czy reguła redukcji problemu 2 ciał do ruchu 1 ciała w polu centralnym: { m->μ, M->m+M }, daje podręcznikowy wynik. *)

Simplify[%]

Plot[{A[[1]],A[[2]],A[[3]]},{t,0,6T},PlotStyle->{Red,Green,Blue},PlotLabel->"Zachowanie składowych wektora Rungego-Lenza" ]

J0=J/.t->0;
A0=A/.t->0;

e=A0/(G m M μ)//Norm (* Ekscentryczność elipsy wyznaczamy z wektora Rungego-Lenza *)

Show[ParametricPlot3D[Evaluate[twoBody[2,"Position",t]-twoBody[1,"Position",t]],{t,0,6T}],
Graphics3D@Arrow[{{0,0,0},J0/m}],
Graphics3D@Arrow[{{0,0,0},-A0 J0.J0/(G m M μ)^2/(1-e)/e}]
,PlotRange->All,PlotLabel->"Elipsa z zaznaczonymi J (prostopadły) i A (w płaszczyźnie)"]

Plot[{Norm[r],A0 J0.J0/(G m M μ)^2/(1-e)/e//Norm,A0 J0.J0/(G m M μ)^2/(1+e)/e//Norm},{t,0,T},PlotRange->{0,All}]

Plot[{VectorAngle[r,A0],Pi},{t,0,T},PlotRange->{0,All}]

Plot[{A0.r/Norm[A0]/Norm[r],J0.(A0\[Cross]r)/Norm[J0]/Norm[A0]/Norm[r]},{t,0,T},PlotLegends->{"Cos[φ]","Sin[φ]"}]

Plot[{ArcTan[A0.r/Norm[A0]/Norm[r],J0.(A0\[Cross]r)/Norm[J0]/Norm[A0]/Norm[r]],VectorAngle[r,A0]},{t,0,T},PlotLegends->{"φ - anomalia prawdziwa","φ - kąt między wektorami"},PlotStyle->{Blue,{Red,Dashed}}]

a=(
(A0 J0.J0/(G m M μ)^2/(1-e)/e//Norm)
+
(A0 J0.J0/(G m M μ)^2/(1+e)/e//Norm)
)/2 (* Ostrzeżenie: nadużywanie notacji postfixowej może prowadzić do bardzo trudno wykrywalnych błędów ! *)

a={A0 J0.J0/(G m M μ)^2/(1-e)/e//Norm,A0 J0.J0/(G m M μ)^2/(1+e)/e//Norm}//Total//(#/2&)

a=J0.J0/(G μ^2 (m+M))/(1-e^2)

Ttrue = 2Pi/ Sqrt[G(m+M)/a^3] (* Prawdziwy okres *)

T/Ttrue

Plot[ArcTan[A0.r/Norm[A0]/Norm[r],J0.(A0\[Cross]r)/Norm[J0]/Norm[A0]/Norm[r]],{t,0,2Ttrue},GridLines->{None,{-Pi,Pi}}]

Plot[ArcTan[A0.r/Norm[A0]/Norm[r],J0.(A0\[Cross]r)/Norm[J0]/Norm[A0]/Norm[r]]+2Pi Ceiling[t/Ttrue],{t,0,2.99Ttrue},GridLines->{None,{-Pi,Pi}}]