3Pi/Quantity[1,"GravitationalConstant"]/Quantity[11,"Minutes"]^2//UnitConvert

Block[{R=1,z=0},Plot[{-1/Abs@r,-2/Pi EllipticK[-4 r R/((r-R)^2+z^2)]/Sqrt[(r-R)^2+z^2]},{r,-4,4},PlotRange->{-2,0}]]

Block[{R=1},ContourPlot[-2/Pi EllipticK[-4 r R/((r-R)^2+z^2)]/Sqrt[(r-R)^2+z^2],{r,-4,4},{z,-4,4},Contours->{0,-0.2,-0.4,-0.6,-0.8,-1},PlotPoints->128,PlotRange->{-2,0}]]

Block[{R=1},ContourPlot[-2/Pi EllipticK[-4 r R/((r-R)^2+z^2)]/Sqrt[(r-R)^2+z^2],{r,0,2},{z,-1,1},Contours->{0,-0.2,-0.4,-0.6,-0.8,-1,-1.2,-1.5,-2},PlotPoints->128,PlotRange->{-2,0}]]

Block[{R=1},ContourPlot[-2/Pi EllipticK[-4 r R/((r-R)^2+z^2)]/Sqrt[(r-R)^2+z^2],{r,0.75,1.25},{z,-0.25,0.25},Contours->{0,-0.2,-0.4,-0.6,-0.8,-1,-1.2,-1.5,-2,-2.2},PlotPoints->128,PlotRange->{-2,0},Axes->True,AxesOrigin->{1,0}]]

R=1(* Promień okregu *)

m=1 (* Masa próbna; nieistotna - skraca się *)

M=1 (* Masa okręgu *)

G=1 (* Stała Grawitacji *)

R0=16 (* Początkowa odległość *)

ω=Sqrt[G M/R0^3] (* Częstość orbity kołowej *)

v0=2/3ω R0 (* Prędkość początkowa *)

T=2Pi/ω (* Oszacowanie okresu orbitanego, przybliżone *)

Tmax=300 (* czas całkowania w okresach orbitalnych, UWAGA : może braknąć pamięci: dla 2 GB Tmax=200, 8 GB Tmax=2000, 32 GB Tmax>10000 *)

r=Sqrt[x[t]^2+y[t]^2]

V=-2 G m M/Pi EllipticK[-4 r R/((r-R)^2+z[t]^2)]/Sqrt[(r-R)^2+z[t]^2] (* Energia potencjalna grawitacyjna jednorodnego okręgu o promieniu R *)

(* Siła = gradient energii potencjalnej *)

Fx=- D[V,x[t]]//Simplify;
Fy=- D[V,y[t]]//Simplify;
Fz=-D[V,z[t]]//Simplify;

(* Równania Newtona *)

eqX=m x''[t]==Fx;
eqY=m y''[t]==Fy;
eqZ=m z''[t]==Fz;

(* Rozwiązanie numeryczne *)

start=Now

rozw=NDSolve[{eqX,eqY,eqZ,x[0]==R0,y[0]==0,x'[0]==v0/3,y'[0]==2/3v0,z[0]==0,z'[0]==6/8v0},{x,y,z},{t,- Tmax T,Tmax T}(*,Method->"ImplicitRungeKutta",WorkingPrecision->32,MaxSteps->10^6*)][[1]]

finish=Now

DateDifference[start,finish,"Minutes"]

MaxMemoryUsed[]

(* Odczytanie funkcji z rozwiązania powyżej, sprawa techniczna *)

X=x/.rozw;
Y=y/.rozw;
Z = z/.rozw;

(* Podstawowe testy *)

X[1]

Plot[X[t],{t,-1T,1T}]

ParametricPlot3D[{X[t],Y[t],Z[t]},{t,-T,T},
PlotRange->4{{-R0,R0},{-R0,R0},{-R0,R0}},AxesOrigin->{0,0,0},PlotStyle->{Thin,Black},Axes->True,AxesLabel->{"X","Y","Z"} ,BaseStyle->{FontSize->16},PlotPoints->4*128]

(* Test zachowania energii  *)

Etot=-G M/2/Pi 4 EllipticK[-4 Sqrt[X[t]^2+Y[t]^2] R/((Sqrt[X[t]^2+Y[t]^2]-R)^2+Z[t]^2)]/Sqrt[((Sqrt[X[t]^2+Y[t]^2]-R)^2+Z[t]^2)]+1/2 m (X'[t]^2+Y'[t]^2+Z'[t]^2);

Etot0=Etot/.t->0

Plot[1-Etot/Etot0,{t,-Tmax T,Tmax T},PlotRange->All]

(* Sprawdzenie czy wielka półoś elipsy istotnie się nie zmienia *)

Plot[-G m M/Etot/2,{t,-Tmax T,Tmax T},PlotRange->{0,2 R0}]

(* Moment pędu *)

J=m Cross[{X[t],Y[t],Z[t]},{X'[t],Y'[t],Z'[t]}];

(* Składowe wektora J i jego długość *)

Plot[{J[[1]],J[[2]],J[[3]],Norm@J(*,J.J + R^2 Z'[t]^2+2 R Z[t] *((2 EllipticK[-((4 (Sqrt[((X[t])^2)+((Y[t])^2)]))/((((-1+(Sqrt[((X[t])^2)+((Y[t])^2)])))^2)+((Z[t])^2)))])/(π (Sqrt[(((-1+(Sqrt[((X[t])^2)+((Y[t])^2)])))^2)+((Z[t])^2)])))*)},{t,-Tmax T,Tmax T},PlotRange->All,PlotStyle->{Red,Blue,Green,Black,Purple}]

(* Wektor Rungego-Lenza *)

RungeLenz=Cross[ m {X'[t],Y'[t],Z'[t]},
Cross[{X[t],Y[t],Z[t]},m{X'[t],Y'[t],Z'[t]}]]-G m^2 M{X[t],Y[t],Z[t]}/Norm[{X[t],Y[t],Z[t]}] ;

Plot[{RungeLenz[[1]],RungeLenz[[2]],RungeLenz[[3]],Norm@RungeLenz},{t,-Tmax T,Tmax T},PlotRange->All,PlotStyle->{Red,Blue,Green,Black}]

(* Sprawdzenie prostopadłości *)

Plot[Dot[RungeLenz,J],{t,-Tmax T,Tmax T},PlotRange->All,PlotStyle->{Red,Blue,Green}]

(* Ewolucja orbity eliptycznej *)

Manipulate[Show[ParametricPlot3D[{X[t],Y[t],Z[t]},{t,k T,(k+2)T},
PlotRange->3{{-R0,R0},{-R0,R0},{-R0,R0}},
AxesOrigin->{0,0,0},
PlotStyle->{Thin,Black},Axes->True,BaseStyle->{FontSize->12},PlotPoints->128],Graphics3D@Arrow[{{0,0,0},10J/.t->k T}],Graphics3D@Arrow[{{0,0,0},30RungeLenz/.t->k T}]],{k,-Tmax,Tmax-2,1}]

Table[ParametricPlot3D[{{R Cos[t],R Sin[t],0},{X[t],Y[t],Z[t]}},{t,k T,(k+2)T},
PlotRange->3{{-R0,R0},{-R0,R0},{-R0,R0}},PlotStyle->{{Red,Thick},{Thin,Black}},Axes->True,AxesLabel->{"X","Y","Z"}, BaseStyle->{FontSize->12},PlotPoints->128],{k,-Tmax,Tmax-1,Tmax/10}]//Show

(* Krzywe zakreślane przez wektory J i Rungego-Lenza *)

Show[
Graphics3D@Arrow[{{0,0,0},J/.t->0}],
Graphics3D[{Opacity[0.1],Sphere[{0,0,0},Norm@J/.t->0]}],
Graphics3D[{Opacity[0.1],Sphere[{0,0,0},Norm@RungeLenz/.t->0]}],
Graphics3D@Arrow[{{0,0,0},RungeLenz/.t->0}],ParametricPlot3D[{J,RungeLenz},{t,-Tmax T,Tmax T},PlotRange->{{-R0,R0},{-R0,R0},{-R0,R0}},PlotStyle->{Red,Blue,Black},PlotPoints->256]]

(* Mimośród orbity *)

Plot[Norm@RungeLenz/(G m^2 M),{t,-Tmax T,Tmax T},PlotRange->{0,1}]

$Assumptions=p∈Vectors[3]&&r∈Vectors[3]

J=Cross[r,p]

e=Cross[p,J]/(G M m^2)-r/Sqrt[r.r]

J.e//TensorReduce