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}]]

(*Export["H:\\Dydaktyka\\AA\\Okrag_potencjal.png",%]*)

R=0.25(* Promień okregu *)

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

M=1 (* Masa okręgu *)

G=1 (* Stała Grawitacji *)

R0=5 (* 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=200 (* 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 *)

rozw=NDSolve[{eqX,eqY,eqZ,x[0]==R0,y[0]==0,x'[0]==v0,y'[0]==v0,z[0]==0,z'[0]==v0},{x,y,z},{t,-Tmax T,Tmax T},Method->"ExplicitRungeKutta"][[1]]

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

X=x/.rozw;

Y=y/.rozw;

Z = z/.rozw;

(* Podstawowe testy *)

X[1]

Plot[Y[t],{t,-2T,2T}]

ParametricPlot[{X[t],Y[t]},{t,-12T,12T},
PlotRange->4{{-R0,R0},{-R0,R0}}]

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

(* 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},10RungeLenz/.t->k T}]],{k,-Tmax,Tmax-1,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

(* 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);

Plot[Etot,{t,0,100T},PlotRange->All]

(* Sprawdzenie czy wielka półoś elipsy istotnie się 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,Norm@J},{t,-Tmax T,Tmax T},PlotRange->All,PlotStyle->{Red,Blue,Green,Black}]

(* 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,{t,-Tmax T,Tmax T},PlotRange->All,PlotStyle->{Red,Blue,Green}]

(* Sprawdzenie prostopadłości *)

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

(* 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}]]

(* Mimośród orbity *)

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