© Kamal K. Barley, Andreas Ruffing, and Sergei K. Suslov
“Old Quantum Mechanics by Bohr and Sommerfeld from Modern Perspective”

© К. Барли, А. Руффинг, С.К. Суслов
“Старая квантовая механика Бора и Зоммерфельда с современной точки зрения”

Успехи физических наук, Январь 2026 г., Том 196, № 1, с. 83-103

(Last modified/executed on December 4, 2026;  4:15 PM Arizona time.)

Animation of the relativistic Kepler motion

(* ::Package::*)(* ::Title::*)(*Relativistic Kepler Ellipse:Animate*)(* ::Section::*)(*Setup& Derived Parameters*)ClearAll["Global`*"];

(*Physical/notation defaults (dimensionless unless you set a0 below)*)
a0:=1;                       (*set a0=1 for dimensionless;or Bohr radius if you want SI scaling*)
αDefault=1/137.036;          (*fine-structure constant*)
ZDefault=79;                 (*gold,matching your example*)
ntDefault=1;nrDefault=1;  (*quantum numbers used in your notebook*)

(*Derived parameters from the notes:ω=Sqrt[nt^2-(Z α)^2]/nt ε=Sqrt[nr (nr+2 Sqrt[nt^2-(Z α)^2])]/(nr+Sqrt[nt^2-(Z α)^2]) a=(a0/Z) (nr+Sqrt[nt^2-(Z α)^2]) Sqrt[(Z α)^2+(nr+Sqrt[nt^2-(Z α)^2])^2]*)

ω[Z_,α_,nt_]:=Sqrt[nt^2-(Z α)^2]/nt;

ε[Z_,α_,nt_,nr_]:=Module[{Δ=Sqrt[nt^2-(Z α)^2]},Sqrt[nr (nr+2 Δ)]/(nr+Δ)];

aParam[Z_,α_,nt_,nr_,a0_:1]:=Module[{Δ=Sqrt[nt^2-(Z α)^2]},(a0/Z) (nr+Δ) Sqrt[(Z α)^2+(nr+Δ)^2]];

(*Polar radius r(θ) using the relativistic ellipse 1/r=(1+ε Cos[ω θ])/(a(1-ε^2))*)
rOf[θ_,a_,eps_,om_]:=a (1-eps^2)/(1+eps Cos[om θ]);

(*Cartesian param from polar (helpful for a moving point)*)
xy[θ_,a_,eps_,om_]:=With[{r=rOf[θ,a,eps,om]},{r Cos[θ],r Sin[θ]}];

(*A helper to guard parameter ranges (nt must exceed Z α).*)
validParamsQ[Z_,α_,nt_]:=nt>Z α


(* ::Section::*)
(*Interactive Controls (Manipulate)*)

Manipulate[Module[{ok=validParamsQ[Z,α,nt],om,eps,a,rmin,rmax,θmax,basePlotRange},If[!ok,Grid[{{Style["Constraint: nt must be greater than Z α.",14,Red]}}],(*else compute and draw*)om=ω[Z,α,nt];
eps=ε[Z,α,nt,nr];
a=aParam[Z,α,nt,nr,a0Scale];
rmin=a (1-eps);(*perihelion*)rmax=a (1+eps);(*aphelion*)θmax=2 π turns;(*total sweep angle*)basePlotRange={{-1.1 rmax,1.1 rmax},{-1.1 rmax,1.1 rmax}};
Column[{Row[{Style["Derived:  ",12,GrayLevel[0.3]],Style[Row[{"ω = ",NumberForm[om,{6,6}],"   ε = ",NumberForm[eps,{6,6}],"   a = ",NumberForm[a,{6,6}],"   rmin = ",NumberForm[rmin,{6,6}],"   rmax = ",NumberForm[rmax,{6,6}]}],12]}],(*Static polar curve over requested number of turns to show precession*)Show[ParametricPlot[xy[θ,a,eps,om],{θ,0,θmax},PlotRange->basePlotRange,PlotPoints->300,MaxRecursion->2,AspectRatio->1,Axes->True,AxesLabel->{"x","y"},PlotLabel->Row[{"Orbit (",turns," turn",If[turns==1,"","s"],")"}],ImageSize->Large],Graphics[{{Red,AbsolutePointSize[6],Point[{0,0}]},(*nucleus at origin*)If[showPeriAp,{Directive[GrayLevel[0.3],Dashed],Circle[{0,0},rmin],Circle[{0,0},rmax]},{}]}]]}]]],{{Z,ZDefault,"Z"},1,92,1,Appearance->"Labeled"},{{α,αDefault,"α"},1/150,1/120,Appearance->"Labeled"},Delimiter,{{nt,ntDefault,"nθ (azimuthal)"},1,10,1,Appearance->"Labeled"},{{nr,nrDefault,"nr (radial)"},0,10,1,Appearance->"Labeled"},Delimiter,{{turns,4.25,"turns to draw"},0.5,12,0.25,Appearance->"Labeled"},{{a0Scale,1,"a0 scale"},{1->"dimensionless (a0=1)",0.52921*10^-8->"Bohr radius (cm)"}},{{showPeriAp,True,"show rmin/rmax circles"},{True,False}},TrackedSymbols:>True,SynchronousUpdating->False]


(* ::Section::*)
(*Pure Animate of the Moving Point (with optional trail)*)

With[{Z=ZDefault,α=αDefault,nt=ntDefault,nr=nrDefault,a0Scale=1},If[!validParamsQ[Z,α,nt],Style["Adjust nt so that nt > Z α",14,Red],Module[{om=ω[Z,α,nt],eps=ε[Z,α,nt,nr],a=aParam[Z,α,nt,nr,a0Scale]},With[{rmin=a (1-eps),rmax=a (1+eps),baseRange=1.2},Animate[Module[{pt=xy[θ,a,eps,om],trailSteps=500,θ0=Max[0,θ-4 π]},Show[{(*trail*)ParametricPlot[xy[ϕ,a,eps,om],{ϕ,θ0,θ},PlotRange->{{-baseRange rmax,baseRange rmax},{-baseRange rmax,baseRange rmax}},PlotPoints->100,MaxRecursion->1,PlotStyle->{Directive[Thickness[0.002]]},PerformanceGoal->"Quality"],(*current point*)Graphics[{{Red,AbsolutePointSize[6],Point[{0,0}]},(*nucleus*){Black,AbsolutePointSize[8],Point[pt]}}]},Axes->True,AspectRatio->1,ImageSize->Large]],{θ,0,10 2 π},(*sweep angle for animation;increase for more precession shown*)AnimationRate->1.0,AnimationRepetitions->1]]]]];


(* ::Section::*)
(*Optional:Export helpers*)

exportOrbit[path_String:"relativistic-ellipse.mp4",θmax_:(10 2 π)]:=Module[{Z=ZDefault,α=αDefault,nt=ntDefault,nr=nrDefault,a0Scale=1,om,eps,a,rmin,rmax,baseRange=1.25,frames=240},If[!validParamsQ[Z,α,nt],Return[$Failed]];
om=ω[Z,α,nt];eps=ε[Z,α,nt,nr];a=aParam[Z,α,nt,nr,a0Scale];
rmin=a (1-eps);rmax=a (1+eps);
Export[path,Table[Show[ParametricPlot[xy[ϕ,a,eps,om],{ϕ,0,θ},PlotRange->{{-baseRange rmax,baseRange rmax},{-baseRange rmax,baseRange rmax}},PlotPoints->100,MaxRecursion->1,PlotStyle->Thickness[0.002]],Graphics[{Red,AbsolutePointSize[6],Point[{0,0}]}],Axes->True,AspectRatio->1,ImageSize->800],{θ,Subdivide[0,θmax,frames]}],"MP4","VideoEncoding"->"H264","FrameRate"->30]];

Ellipses Au