THE SPHERICAL PENDULUM WITH NON-ZERO ANGULAR MOMENTUM 1
K.6.3.3.b
Franz Michel, Bernhard Schnizer

Institut fuer Theoretische Physik
Technische Universitaet Graz
A-8010  Graz, Austria
Projekt Multimedial Lehre des
Zentralen Informatikdienstes der TU Graz

2002-01-31

After the second section (Proof of Analytical Relations) has been run, any of the subsections 
of the section Static Plots runs independently.

Introduction

Proof of Analytical Relations

Clear["@", R]
Needs["Graphics`Graphics`"]
Needs["Graphics`Colors`"]

(Subscript[ω, 0]) =  (Sqrt[(g/R)]);                                                                 (*  Angular Frequency for small oscillations *)

T0=PowerExpand[((2 π)/(Subscript[ω, 0]))]                                  (*   Period for small oscillations  *)

Comparison of coefficients for the two representations of polynomial P(u)

In this section eqs.(4) and (5)  are proved.

pu  =  (u  -  u1) (u  -  u2) (u  -  
              u3);
pl  =  Collect[pu // Expand, u];  
    Print[pl // TraditionalForm]

pg  =  
      u^3  -  ϵ u^2  -  
        u  +  ϵ  - λ^2/2 ;  
    Print[pg // TraditionalForm]

cu  =  
      Table[Coefficient[u pl // Expand, u^k], {k, Length[pl]}]
cg  =  
      Table[Coefficient[u pg // Expand, u^k], {k, 
          Length[pg] - 1}]

sel  =  
      Solve[cg  ==  cu,  {ϵ,  λ}]

sel  =  sel[[2]]

k  =  2;
su3  =  
      Solve[cg[[k]]  ==  cu[[k]],  u3] // 
        Flatten

Expressing  ϵ  and  λ  by   (Subscript[u, 1]) and  (Subscript[u, 2]).

fe[u1_, u2_]   =  (ϵ  /.  sel)  /.  su3

fl[u1_, u2_]   =  (λ  /.  sel)  /.  su3  // 
        FullSimplify

Transforming the Equation of Motion for  u  =  cos θ  to Standard Elliptic Integral

ps  =  (u2 - u1)^2 (u3   -  u1) 
        Sin[ψ]^2 Cos[ψ]^2 (1  -  
            kq  Sin[ψ]^2)

uu  =  u1  +  (u2   -  u1) Sin[ψ]^2

duu  =  D[uu, ψ]

PowerExpand[duu/Sqrt[ps]]

Implementation of the Expression for the Azimuth in Eq.(11)

Numeric Implementation

Azimuth  ϕ  expressed as function of auxiliary variable  ψ. Function  ϕ(ψ) found by numerical integration of integral  (11).

phn[ψ_]  :=  λ/
          su  NIntegrate[1/
              Sqrt[1  -  
                  kq Sin[ψψ]^2]/(1  -  (u1  +  (u2  -  u1) Sin[ψψ]^2 )^2), {ψψ, 0, ψ}]

Analytic Implementation

Azimuth  ϕ  expressed as function of auxiliary variable  ψ. Function  ϕ(ψ) expressed by elliptic integrals of third kind, eq. (11).

n1  =  (u1  -  u2)/(u1  -  1);
n2  =  (u1  -  u2)/(u1  +  1) ;

php[ψ_] :=  (((Sqrt[1/2]) fl[u1, u2])/(Sqrt[u3  -   u1]))((1/(1   -    u1))EllipticPi[n1,  ψ,  kq]+(1/(1   +    u1))EllipticPi[n2,  ψ,  
                    kq]) /. su3 ;

Proof of the Decomposition for Elliptic Integrals of 3rd kind in Eq.(11)

Check of Identity (16)

au=(1/(1-((((Subscript[u, 1])+((Subscript[u, 2])-(Subscript[u, 1]))(sin^2)β))^2)));

α(Subscript[q, 1]) = (((Subscript[u, 1]) - (Subscript[u, 2]))/((Subscript[u, 1]) - 1));

α(Subscript[q, 2]) = (((Subscript[u, 1]) - (Subscript[u, 2]))/((Subscript[u, 1]) + 1));

au == ((1/2)/(1   -    (Subscript[u, 1])))(1/((1  -     α(Subscript[q, 1])  (sin^2) β)))+ ((1/2)/(1   +    (Subscript[u, 1])))(1/((1  -   α(Subscript[q, 2])  (sin^2) β)))//Simplify//ExpandAll

Perfoming the decomposition in  (16)

an  =  Apart[au]

cf  =  1/2;
nn  =   {an[[1]] // Factor, an[[2]]} /
        cf

vf = {1 - (Subscript[u, 1]), 1  +  (Subscript[u, 1])}^-1;

cc  =  Collect[Denominator[nn], (sin^2) β] vf

aq  =  - 
            Coefficient[Collect[cc, (sin^2) β], (sin^2) β] // 
          ExpandAll // Together

cf vf  . (1  -   aq (sin^2) β)^-1

Static Plots

Small Oscillation  (0 < 2° ≤  δ  ≤  6°)

Input Parameters. Choice of latitudes

Resulting Parameters

\(u1\  = \ \ Cos[θ1] // N\ ; \ 
    Print[\*""\!\(u\_1\) = "", u1]\)
\(u2\  = \ \ Cos[θ2] // N; \ \ Print[\*""\!\(u\_2\) = "", 
      u2]\)
\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Print[\*""(Subscript[u, 3]) = "", 
      u3\  //. \ su3\ ]\)

ϵ     =  (ϵ   /. sel) /.  
        su3 ;   Print["ϵ = ", ϵ]
λ     =  (λ  /. sel) /.  
        su3;    Print["λ = ", λ]
s1  =  Sin[θ1] // Chop ;   Print["r1 = ", 
      s1, "  = Relative radius of lower circle of latitude "]
s2  =  Sin[θ2] // N;            Print["r2 = ", 
      s2, "  = Relative radius of higher circle of latitude "]

\(kq\  = \ \((u2 - u1)\)/\((u3 - u1)\)\  /. \ su3; \ 
    Print[\*""m = \!\(k\^2\) = "", 
      kq, \ "  = Parameter of elliptic integrals"]\)

su  =  Sqrt[(u3  -  u1)/2]  /.  su3

Period

tt =  
      4  EllipticK[kq]/
            su /(2 π);     Print["T/T0 = ", 
      tt]

u  =  cos θ   expressed by auxiliary variable  ψ.

fu[ψ_]  =  uu

Choosing an Expression for the Azimuth

As a check the two definitions of  ϕ(ψ)  are compared for some numerical values:

Table[{phn[ψ], php[ψ]}, {ψ, 
        0, π/2, π/10}]

The computing times of both approaches are compared:

tn  =  Timing[
          Table[phn[ψ], {ψ, 0, 
              10  π, π/20}]][[1]]

tp  =  Timing[
          Table[php[ψ], {ψ, 0, 
              10  π, π/20}]][[1]]

So it is more expedient to use the analytic formula for  ϕ(ψ) :

ph[ψ_] := php[ψ]

Curves for θ(ψ),   (Subscript[ω, 0]) t(ψ)  and  ϕ(ψ).

ppth = Plot[( π-ArcCos[fu[ψ]]) 180/π,{ψ, 0, 6  π},
                      \(AxesLabel\  -> \ {"ψ", "θ [°]"}\), TextStyle  ->  {FontSize  ->   12},
                      Ticks  ->  {PiScale, 
                Automatic},PlotRange  ->  {1, 7},
		         DisplayFunction  ->  Identity];

\(\(ppt\  = \ 
        Plot[{ψ, EllipticF[ψ, kq]/su}, {ψ, 0, 
            6  π}, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle\  -> {GrayLevel[0], {Thickness[ .01], 
                Hue[0]}}, 
		DisplayFunction\  -> \ 
            Identity, \ 
          Ticks\  -> \ 
            PiScale, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel\  -> \ {"ψ", \*""(Subscript[ω, 0])t""}, \ 
          TextStyle\  -> \ {FontSize\  -> \ \ 12}];\)\)

\(\(pph\  = \ 
        Plot[{ψ, ph[ψ]}, {ψ, 0, 
            6  π}, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle\  -> {GrayLevel[0], {Thickness[ .01], 
                Hue[0]}}, 
		\ DisplayFunction\  -> \ Identity, \ \ Ticks\  -> \ 
            PiScale, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel\  -> \ {"ψ", "ϕ"}, \ 
          TextStyle\  -> \ {FontSize\  -> \ \ 12}];\)\)

Show[GraphicsArray[{{ppth}, {ppt, pph}}],  
        ImageSize  ->  500];

The Projection of the Trajectory on the x,y-plane

R  =   .1;

\(\(pr\  = \ 
        ParametricPlot[
          R\ Sqrt[1\  - \ fu[ψ]^2] {Sin[ph[ψ]], 
              Cos[ph[ψ]]}, {ψ, 0, 2  π}, 
	
          AspectRatio\  -> \ Automatic, \ 
          ImageSize\  -> \ 450, 
	
          PlotStyle\  -> {{Thickness[ .003], 
                Hue[0]}}, 
	
          Epilog\  -> \ {Circle[{0, 0}, R\ s1], Circle[{0, 0}, R\ s2], 
              Circle[{0, 0}, R\ ], 
	Hue[0], 
              Thickness[ .003], 
              Line[{{0, 0}, 
                  R\ s1 {Sin[ph[0]], Cos[ph[0]]}}], 
	
              Hue[ .7], 
              Line[{{0, 0}, 
                  0.5\ R\ s1 {Sin[ph[2  π]], 
                      Cos[ph[2  π]]}}]}];\)\)

Fig. 6.9

Show[pr, 
        PlotRange  ->  
          R {{-1, 
                1},  .15 {-1, 
                  1}},  
        Ticks  ->  
          None,  ImageSize  ->  600];

The Trajectory on the Sphere

pd  =  
        ParametricPlot3D[{Sqrt[1  -  fu[ψ]^2] Sin[ph[ψ]], 
            Sqrt[1  -  fu[ψ]^2] Cos[ph[ψ]],  
            fu[ψ], {Thickness[ .005], Hue[0]}}, {ψ, 0, 2  π},  
          BoxRatios  ->  {1, 1,  .1}];

sph  =  
        ParametricPlot3D[{Sin[th] Cos[pp], Sin[th] Sin[pp], Cos[th]}, {th, 
            0, π}, {pp, 0, 2  π}];

Show[pd, sph,  ViewPoint  -> {0, 0.1, -2},  
        ImageSize  ->  450];

Large Oscillation 1  (0 < 20° ≤  δ  ≤  70°)

Input Parameters. Choice of latitudes

Resulting Parameters

\(u1\  = \ \ Cos[θ1] // N\ ; \ 
    Print[\*""\!\(u\_1\) = "", u1]\)
\(u2\  = \ \ Cos[θ2] // N; \ \ Print[\*""\!\(u\_2\) = "", 
      u2]\)
\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Print[\*""(Subscript[u, 3]) = "", 
      u3\  //. \ su3\ ]\)

ϵ     =  (ϵ   /. sel) /.  
        su3 ;   Print["ϵ = ", ϵ]
λ     =  (λ  /. sel) /.  
        su3;    Print["λ = ", λ]
s1  =  Sin[θ1] // Chop ;   Print["r1 = ", 
      s1, "  = Relative radius of lower circle of latitude "]
s2  =  Sin[θ2] // N;            Print["r2 = ", 
      s2, "  = Relative radius of higher circle of latitude "]

\(kq\  = \ \((u2 - u1)\)/\((u3 - u1)\)\  /. \ su3; \ 
    Print[\*""m = \!\(k\^2\) = "", 
      kq, \ "  = Parameter of elliptic integrals"]\)

su  =  Sqrt[(u3  -  u1)/2]  /.  su3

tt =  
      4  EllipticK[kq]/
            su /(2 π);     Print["T/T0 = ", 
      tt]

u  =  cos θ   expressed by auxiliary variable  ψ.

fu[ψ_]  =  uu

ph[ψ_] := php[ψ]

Curves for θ(ψ),   (Subscript[ω, 0]) t(ψ)  and  ϕ(ψ). Fig.6.13

ppth = Plot[( π-ArcCos[fu[ψ]]) 180/π,{ψ, 0, 6  π},
                      \(AxesLabel\  -> \ {"ψ", "θ [°]"}\), TextStyle  ->  {FontSize  ->   12},
                      Ticks  ->  {PiScale, 
                Automatic},PlotRange  ->  {10, 70},
		         DisplayFunction  ->  Identity];

\(\(ppt\  = \ 
        Plot[{ψ, EllipticF[ψ, kq]/su}, {ψ, 0, 
            6  π}, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle\  -> {GrayLevel[0], {Thickness[ .01], 
                Hue[0]}}, 
		DisplayFunction\  -> \ 
            Identity, \ 
          Ticks\  -> \ 
            PiScale, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel\  -> \ {"ψ", \*""(Subscript[ω, 0])t""}, \ 
          TextStyle\  -> \ {FontSize\  -> \ \ 12}, 
		Epilog\  -> \ {Text[\*""(Subscript[ω, 0])t = ψ"", π {5, 3}], 
              Hue[0], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[
                Subscript[\*""(Subscript[ω, 0])t"", 
                  ex], π {2.5, 5}]}];\)\)

\(\(pph\  = \ 
        Plot[{ψ, ph[ψ]}, {ψ, 0, 
            6  π}, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle\  -> {GrayLevel[0], {Thickness[ .01], 
                Hue[0]}}, 
		\ DisplayFunction\  -> \ Identity, \ \ Ticks\  -> \ 
            PiScale, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel\  -> \ {"ψ", "ϕ"}, \ 
          TextStyle\  -> \ {FontSize\  -> \ \ 12}, 
		Epilog\  -> \ {Text["ϕ = ψ", π {5, 3}], Hue[
                0], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[
                Subscript["ϕ(ψ)", 
                  ex], π {2.5, 
                    5}], 
				\ Hue[ .6], 
              Line[π {{2.15, 2.29}, {3, 1.6}}], \ 
              Text["Δϕ", π {3.4, 
                    1.5}], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Thickness[ .01], 
              Line[{π {2, 2}, {2  π, ph[2  π]}}]}];\)\)

Show[GraphicsArray[{{ppth}, {ppt, pph}}],  
        ImageSize  ->  500];

Projection of Trajectory on  x,y-plane

R  =   .1;

\(\(pr\  = \ 
        ParametricPlot[
          R\ Sqrt[1\  - \ fu[ψ]^2] {Sin[ph[ψ]], 
              Cos[ph[ψ]]}, {ψ, 0, 
            2  π}, 
\ \ \ \ \ \ \ \ \ \ \ AspectRatio\  -> \ Automatic, 
          ImageSize\  -> \ 
            450, 
\ \ \ \ \ \ \ \ \ \ \ \ PlotStyle\  -> {{Thickness[ .005], 
                Hue[0]}}, \ 
	\ Epilog\  -> \ {Circle[{0, 0}, R\ s1], Circle[{0, 0}, R\ s2], 
              Circle[{0, 0}, R\ ], 
	Thickness[ .005], 
              Hue[0], Line[{{0, 0}, 
                  R\ s1 {Sin[ph[0]], Cos[ph[0]]}}], 
	
              Hue[ .7], 
              Line[{{0, 0}, 
                  R\ s1 {Sin[ph[2  π]], 
                      Cos[ph[2  π]]}}]}];\)\)

ps  =  
        ParametricPlot[
          Evaluate[
            R Sqrt[1  -  fu[ψ]^2] {Sin[ph[ψ]], 
                Cos[ph[ψ]]}], {ψ, 2  π, 
            4  π}, 
AspectRatio  ->  Automatic, 
          PlotStyle  -> {{Thickness[ .005], Hue[0.7],  
                Dashing[{ .03}]}}];

Graphics Fig. 6.14

dg  =  
        Graphics[{Circle[{0, 0}, 
              R s1 0.4, {-ph[2  π] + 5  π/2, π/
                  2}],   Circle[{0, 0}, 
              R s1 0.43, {-ph[2  π] + 5  π/2, π/
                  2}]}];
dh  =  
        Graphics[
          Text["Δϕ", 
            R s1 { .4,  .55}]];
vbxz  =  
        Show[pr, ps, dg, dh, 
          PlotRange  ->  
            R {{-1, 1}, {-1, 
                  1}},  
Ticks  ->  None,  
          TextStyle  ->  {FontSize  ->  18}];

The Trajectory in Space

lt  =  Thickness[0.007];
pd  =  
        ParametricPlot3D[ {R Sqrt[1  -  fu[ψ]^2] Sin[ph[ψ]], 
            R Sqrt[1  -  fu[ψ]^2] Cos[ph[ψ]],  
            R fu[ψ], {lt, Hue[0]}}, {ψ, 0, 2  π},  
          BoxRatios  ->  {1, 1, 1}];

lt  =  Thickness[0.004];
pe1  =  
        ParametricPlot3D[{R Sqrt[1  -  fu[ψ]^2] Sin[ph[ψ]], 
            R Sqrt[1  -  fu[ψ]^2] Cos[ph[ψ]], 
            R fu[ψ],  {Dashing[{ .005,  .5}], lt, Hue[0.7]}}, {ψ, 
            2  π, 3  π},  
          BoxRatios  ->  {1, 1, 1}];

pe2  =  
        ParametricPlot3D[{R Sqrt[1  -  fu[ψ]^2] Sin[ph[ψ]], 
            R Sqrt[1  -  fu[ψ]^2] Cos[ph[ψ]], 
            R fu[ψ],  {Dashing[{ .005,  .5}], lt, Hue[0.7]}}, {ψ, 
            3  π, 4  π},  
          BoxRatios  ->  {1, 1, 1}];

pf  =  
        ParametricPlot3D[{R Sqrt[1  -  fu[ψ]^2] Sin[ph[ψ]], 
            R Sqrt[1  -  fu[ψ]^2] Cos[ph[ψ]], 
            R  fu[ψ], {lt, Hue[0]}}, {ψ, 2  π, 4  π},  
          BoxRatios  ->  {1, 1, 1}];

lt  =  Thickness[ .004];

k1  =  
        ParametricPlot3D[{R s1 Sin[ϕ], R s1 Cos[ϕ], R u1, 
            lt}, {ϕ, 0, 2  π},  BoxRatios  ->  {1, 1, 1},  
          DisplayFunction  ->  Identity];
k2  =  
        ParametricPlot3D[{R s2 Sin[ϕ], R s2 Cos[ϕ], R u2, 
            lt}, {ϕ, 0, 2  π},  BoxRatios  ->  {1, 1, 1},  
          DisplayFunction  ->  Identity];
ϕ3  =    π/4  -  0.1;
k3  =  
        ParametricPlot3D[{R Sin[θ] Sin[ϕ3], 
            R Sin[θ]  Cos[ϕ3], R Cos[θ], 
            lt}, {θ, θ1, θ2},  
          BoxRatios  ->  {1, 1, 1},  
          DisplayFunction  ->  Identity];
ϕ4  =    ϕ3  +  π;
k4  =  
        ParametricPlot3D[{R Sin[θ] Sin[ϕ4], 
            R Sin[θ]  Cos[ϕ4], R Cos[θ], 
            lt}, {θ, θ1, θ2},  
          BoxRatios  ->  {1, 1, 1},  
          DisplayFunction  ->  Identity];

lh  =   .005;  vh  = -R;
fp  =  
        Graphics3D[{Line[{{lh, 0, vh}, {-lh, 0, vh}}], 
            Line[{{0, lh, vh}, {0, -lh, vh}}]}];

vb  =  
        Show[pd, pe1, pe2, k1, k2, k3, k4, fp, 
          BoxRatios  ->  {1, 1,  .5},  Axes  ->  False,  
          ImageSize  ->  500];

fvb  = 
        Show[vb, PlotRange  ->  {- .1, - .015}, 
          Axes  ->  False, 
          ViewPoint -> {-4.275,  -0.643,  1.8}];

po1  =  
        ParametricPlot3D[{Sin[th] Cos[pp], Sin[th] Sin[pp], Cos[th]}, {th, 
            0, ArcCos[u2]}, {pp, 0, 2  Pi}];

Show[ke, po1];

po2  =  
        ParametricPlot3D[{Sin[th] Cos[pp], Sin[th] Sin[pp], 
            Cos[th]}, {th, π, ArcCos[u1]}, {pp, 0, 2  Pi}];

re  =  1.005;
ke =  
        ParametricPlot3D[{re Sin[ϕ], re   Cos[ϕ], 
            0, {Thickness[ .008], GrayLevel[0]}}, {ϕ, 0, 2  π},  
          BoxRatios  ->  {1, 1, 1}];

sw1  =  
        Show[po1  /.  Polygon[x___]  ->  {EdgeForm[], Polygon[x]},  
          ke];

sw2  =  
        Show[po2  /.  
            Polygon[x___]  ->  {EdgeForm[], Polygon[x]}];

po  =  
        ParametricPlot3D[{Sin[th] Cos[pp], Sin[th] Sin[pp], Cos[th]}, {th, 
            ArcCos[u1], ArcCos[u2]}, {pp, 0, 2  Pi}, 
 
          ViewPoint  ->  {5, 5, -1},  Boxed  ->  False,  
          Axes  ->  None];

pc  =  
        Show[sw1, sw2, po,   ViewPoint  ->  {5, 5, -1.9},  
          Boxed  ->  False,  Axes  ->  None];

The Trajectory in Space and the Permitted Zone on the Sphere, Fig.6.6

Show[GraphicsArray[{pc, fvb}], GraphicsSpacing  -> { .1, 0},  
        ImageSize  ->  500];

Further Representation of the Trajectory in Space

sph  =  
        ParametricPlot3D[
          R {Sin[th] Cos[pp], Sin[th] Sin[pp], 
              Cos[th]}, 
{th, 0, π}, {pp, 0, 
            2  π},  DisplayFunction  ->  Identity];
sn  =  
        Show[sph  /.  Polygon[x___]  ->  {EdgeForm[], Polygon[x]},  
          DisplayFunction  ->  Identity];

Show[pd, pe, sn, k1, k2,  ViewPoint  -> {0, 0.1, -2},  
        ImageSize  ->   500];

Show[pd, pe, k1, k2, sn,  ViewPoint  -> {0, 2,  .2},  
        ImageSize  ->  500];

Show[pd, pe, k1, k2, sn,  ViewPoint  -> {0, 2, -1.8},  
        ImageSize  ->  500];

Large Oscillation 2  (0 < 15° ≤  δ  ≤  80°)

Input Parameters. Choice of latitudes

Feel free to use your own values for the two latitudes. But not all values are possible from the physics point of view.
In particular a case with both latitudes in the upper hemisphere is not admissible.

Resulting Parameters

\(u1\  = \ \ Cos[θ1] // N\ ; \ 
    Print[\*""\!\(u\_1\) = "", u1]\)
\(u2\  = \ \ Cos[θ2] // N; \ \ Print[\*""\!\(u\_2\) = "", 
      u2]\)
\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Print[\*""(Subscript[u, 3]) = "", 
      u3\  //. \ su3\ ]\)

ϵ     =  (ϵ   /. sel) /.  
        su3 ;   Print["ϵ = ", ϵ]
λ     =  (λ  /. sel) /.  
        su3;    Print["λ = ", λ]
s1  =  Sin[θ1] // Chop ;   Print["r1 = ", 
      s1, "  = Relative radius of lower circle of latitude "]
s2  =  Sin[θ2] // N;            Print["r2 = ", 
      s2, "  = Relative radius of higher circle of latitude "]

\(kq\  = \ \((u2 - u1)\)/\((u3 - u1)\)\  /. \ su3; \ 
    Print[\*""m = \!\(k\^2\) = "", 
      kq, \ "  = Parameter of elliptic integrals"]\)

su  =  Sqrt[(u3  -  u1)/2]  /.  su3

tt =  
      4  EllipticK[kq]/
            su /(2 π);     Print["T/T0 = ", 
      tt]

u  =  cos θ   expressed by auxiliary variable  ψ.

fu[ψ_]  =  uu

ph[ψ_] := php[ψ]

Curves for θ(ψ),   (Subscript[ω, 0]) t(ψ)  and  ϕ(ψ).

ppth = Plot[( π-ArcCos[fu[ψ]]) 180/π,{ψ, 0, 6  π},
                      \(AxesLabel\  -> \ {"ψ", "θ [°]"}\), TextStyle  ->  {FontSize  ->   12},
                      Ticks  ->  {PiScale, 
                Automatic},PlotRange  ->  {30, 80},
		         DisplayFunction  ->  Identity];

\(\(ppt\  = \ 
        Plot[{ψ, EllipticF[ψ, kq]/su}, {ψ, 0, 
            6  π}, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle\  -> {GrayLevel[0], {Thickness[ .01], 
                Hue[0]}}, 
		DisplayFunction\  -> \ 
            Identity, \ 
          Ticks\  -> \ 
            PiScale, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel\  -> \ {"ψ", \*""(Subscript[ω, 0])t""}, \ 
          TextStyle\  -> \ {FontSize\  -> \ \ 12}, 
		Epilog\  -> \ {Text[\*""(Subscript[ω, 0])t = ψ"", π {5, 3}], 
              Hue[0], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[
                Subscript[\*""(Subscript[ω, 0])t"", 
                  ex], π {2.5, 5}]}];\)\)

\(\(pph\  = \ 
        Plot[{ψ, ph[ψ]}, {ψ, 0, 
            6  π}, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle\  -> {GrayLevel[0], {Thickness[ .01], 
                Hue[0]}}, 
		\ DisplayFunction\  -> \ Identity, \ \ Ticks\  -> \ 
            PiScale, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel\  -> \ {"ψ", "ϕ"}, \ 
          TextStyle\  -> \ {FontSize\  -> \ \ 12}, 
		Epilog\  -> \ {Text["ϕ = ψ", π {5, 3}], Hue[
                0], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[
                Subscript["ϕ(ψ)", 
                  ex], π {2.5, 
                    5}], 
				\ Hue[ .6], 
              Line[π {{2.15, 2.29}, {3, 1.6}}], \ 
              Text["Δϕ", π {3.4, 
                    1.5}], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Thickness[ .01], 
              Line[{π {2, 2}, {2  π, ph[2  π]}}]}];\)\)

Show[GraphicsArray[{{ppth}, {ppt, pph}}],  
        ImageSize  ->  500];

Projection of Trajectory on  x,y-plane

R  =   .1;

\(\(pr\  = \ 
        ParametricPlot[
          R\ Sqrt[1\  - \ fu[ψ]^2] {Sin[ph[ψ]], 
              Cos[ph[ψ]]}, {ψ, 0, 
            2  π}, 
\ \ \ \ \ \ \ \ \ \ \ AspectRatio\  -> \ Automatic, 
          ImageSize\  -> \ 
            450, 
\ \ \ \ \ \ \ \ \ \ \ \ PlotStyle\  -> {{Thickness[ .005], 
                Hue[0]}}, \ 
	\ Epilog\  -> \ {Circle[{0, 0}, R\ s1], Circle[{0, 0}, R\ s2], 
              Circle[{0, 0}, R\ ], 
	Thickness[ .005], 
              Hue[0], Line[{{0, 0}, 
                  R\ s1 {Sin[ph[0]], Cos[ph[0]]}}], 
	
              Hue[ .7], 
              Line[{{0, 0}, 
                  R\ s1 {Sin[ph[2  π]], 
                      Cos[ph[2  π]]}}]}];\)\)

ps  =  
        ParametricPlot[
          R Sqrt[1  -  fu[ψ]^2] {Sin[ph[ψ]], 
              Cos[ph[ψ]]}, {ψ, 2  π, 
            4  π}, 
AspectRatio  ->  Automatic, 
          PlotStyle  -> {{Thickness[ .005], Hue[0.7]}}];

dg  =  
        Graphics[{Circle[{0, 0}, 
              R s1 0.4, {-ph[2  π] + 5  π/2, π/
                  2}],   Circle[{0, 0}, 
              R s1 0.43, {-ph[2  π] + 5  π/2, π/
                  2}]}];
dh  =  
        Graphics[
          Text["Δϕ", 
            R s1 { .4,  .55}]];
vbxz  =  
        Show[pr, ps, dg, dh, 
          PlotRange  ->  
            R {{-1, 1}, {-1, 
                  1}},  
Ticks  ->  None,  
          TextStyle  ->  {FontSize  ->  18}];

The Trajectory in Space

pd  =  
        ParametricPlot3D[ {R Sqrt[1  -  fu[ψ]^2] Sin[ph[ψ]], 
            R Sqrt[1  -  fu[ψ]^2] Cos[ph[ψ]],  
            R fu[ψ], {Thickness[ .007], Hue[0]}}, {ψ, 0, 
            2  π},  BoxRatios  ->  {1, 1, 1}];

pe  =  
        ParametricPlot3D[{R Sqrt[1  -  fu[ψ]^2] Sin[ph[ψ]], 
            R Sqrt[1  -  fu[ψ]^2] Cos[ph[ψ]], 
            R fu[ψ], {Thickness[ .007], Hue[0.7]}}, {ψ, 2  π, 
            4  π},  BoxRatios  ->  {1, 1, 1}];

pf  =  
        ParametricPlot3D[{R Sqrt[1  -  fu[ψ]^2] Sin[ph[ψ]], 
            R Sqrt[1  -  fu[ψ]^2] Cos[ph[ψ]], 
            R  fu[ψ], {Thickness[ .007], Hue[0]}}, {ψ, 2  π, 
            4  π},  BoxRatios  ->  {1, 1, 1}];

k1  =  
        ParametricPlot3D[{R s1 Sin[ϕ], R s1 Cos[ϕ], 
            R u1, {Thickness[ .005]}}, {ϕ, 0, 2  π},  
          BoxRatios  ->  {1, 1, 1},  
          DisplayFunction  ->  Identity];
k2  =  
        ParametricPlot3D[{R s2 Sin[ϕ], R s2 Cos[ϕ], 
            R u2, {Thickness[ .005]}}, {ϕ, 0, 2  π},  
          BoxRatios  ->  {1, 1, 1},  
          DisplayFunction  ->  Identity];
ϕ3  =    π/4  +   .1;
k3  =  
        ParametricPlot3D[{R Sin[θ] Sin[ϕ3], 
            R Sin[θ]  Cos[ϕ3], 
            R Cos[θ], {Thickness[ .005]}}, {θ, θ1, θ2},  BoxRatios  ->  {1, 1, 1},  
          DisplayFunction  ->  Identity];
ϕ4  =    ϕ3  +  π;
k4  =  
        ParametricPlot3D[{R Sin[θ] Sin[ϕ4], 
            R Sin[θ]  Cos[ϕ4], 
            R Cos[θ], {Thickness[ .005]}}, {θ, θ1, θ2},  BoxRatios  ->  {1, 1, 1},  
          DisplayFunction  ->  Identity];

Show[pd, pe, k1, k2, k3, k4, BoxRatios  ->  {1, 1,  .5}, 
        Axes  ->  False, ImageSize  ->  500];

sph  =  
        ParametricPlot3D[
          R {Sin[th] Cos[pp], Sin[th] Sin[pp], 
              Cos[th]}, 
{th, 0, π}, {pp, 0, 
            2  π},  DisplayFunction  ->  Identity];
sn  =  
        Show[sph  /.  Polygon[x___]  ->  {EdgeForm[], Polygon[x]},  
          DisplayFunction  ->  Identity];

Show[pd, pe, sn, k1, k2,  ViewPoint  -> {0, 0.1, -2},  
        ImageSize  ->   500];

Show[pd, pe, k1, k2, sn,  ViewPoint  -> {0, 2,  .2},  
        ImageSize  ->  500];

Show[pd, pe, k1, k2, sn,  ViewPoint  -> {0, 2, -1.8},  
        ImageSize  ->  500];

Large Oscillation 3  (0 < 35° ≤  δ  ≤  130°)

Input Parameters. Choice of latitudes

Resulting Parameters

\(u1\  = \ \ Cos[θ1] // N\ ; \ 
    Print[\*""\!\(u\_1\) = "", u1]\)
\(u2\  = \ \ Cos[θ2] // N; \ \ Print[\*""\!\(u\_2\) = "", 
      u2]\)
\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Print[\*""(Subscript[u, 3]) = "", 
      u3\  //. \ su3\ ]\)

ϵ     =  (ϵ   /. sel) /.  
        su3 ;   Print["ϵ = ", ϵ]
λ     =  (λ  /. sel) /.  
        su3;    Print["λ = ", λ]
s1  =  Sin[θ1] // Chop ;   Print["r1 = ", 
      s1, "  = Relative radius of lower circle of latitude "]
s2  =  Sin[θ2] // N;            Print["r2 = ", 
      s2, "  = Relative radius of higher circle of latitude "]

\(kq\  = \ \((u2 - u1)\)/\((u3 - u1)\)\  /. \ su3; \ 
    Print[\*""m = \!\(k\^2\) = "", 
      kq, \ "  = Parameter of elliptic integrals"]\)

su  =  Sqrt[(u3  -  u1)/2]  /.  su3

tt =  
      4  EllipticK[kq]/
            su /(2 π);     Print["T/T0 = ", 
      tt]

u  =  cos θ   expressed by auxiliary variable  ψ.

fu[ψ_]  =  uu

ph[ψ_] := php[ψ]

Curves for θ(ψ),   (Subscript[ω, 0]) t(ψ)  and  ϕ(ψ).

ppth = Plot[( π-ArcCos[fu[ψ]]) 180/π,{ψ, 0, 6  π},
                      \(AxesLabel\  -> \ {"ψ", "θ [°]"}\), TextStyle  ->  {FontSize  ->   12},
                      Ticks  ->  {PiScale, 
                Automatic},PlotRange  ->  {0, 140},
		         DisplayFunction  ->  Identity];

\(\(ppt\  = \ 
        Plot[{ψ, EllipticF[ψ, kq]/su}, {ψ, 0, 
            6  π}, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle\  -> {GrayLevel[0], {Thickness[ .01], 
                Hue[0]}}, 
		DisplayFunction\  -> \ 
            Identity, \ 
          Ticks\  -> \ 
            PiScale, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel\  -> \ {"ψ", \*""(Subscript[ω, 0])t""}, \ 
          TextStyle\  -> \ {FontSize\  -> \ \ 12}, 
		Epilog\  -> \ {Text[\*""(Subscript[ω, 0])t = ψ"", π {5, 3}], 
              Hue[0], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[
                Subscript[\*""(Subscript[ω, 0])t"", 
                  ex], π {2.5, 5}]}];\)\)

\(\(pph\  = \ 
        Plot[{ψ, ph[ψ]}, {ψ, 0, 
            6  π}, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle\  -> {GrayLevel[0], {Thickness[ .01], 
                Hue[0]}}, 
		\ DisplayFunction\  -> \ Identity, \ \ Ticks\  -> \ 
            PiScale, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel\  -> \ {"ψ", "ϕ"}, \ 
          TextStyle\  -> \ {FontSize\  -> \ \ 12}, 
		Epilog\  -> \ {Text["ϕ = ψ", π {5, 3}], Hue[
                0], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[
                Subscript["ϕ(ψ)", 
                  ex], π {2.5, 
                    5}], 
				\ Hue[ .6], 
              Line[π {{2.15, 2.29}, {3, 1.6}}], \ 
              Text["Δϕ", π {3.4, 
                    1.5}], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Thickness[ .01], 
              Line[{π {2, 2}, {2  π, ph[2  π]}}]}];\)\)

Show[GraphicsArray[{{ppth}, {ppt, pph}}],  
        ImageSize  ->  500];

Projection of Trajectory on  x,y-plane

R  =   .1;

\(\(pr\  = \ 
        ParametricPlot[
          R\ Sqrt[1\  - \ fu[ψ]^2] {Sin[ph[ψ]], 
              Cos[ph[ψ]]}, {ψ, 0, 
            2  π}, 
\ \ \ \ \ \ \ \ \ \ \ AspectRatio\  -> \ Automatic, 
          ImageSize\  -> \ 
            450, 
\ \ \ \ \ \ \ \ \ \ \ \ PlotStyle\  -> {{Thickness[ .005], 
                Hue[0]}}, \ 
	\ Epilog\  -> \ {Circle[{0, 0}, R\ s1], Circle[{0, 0}, R\ s2], 
              Circle[{0, 0}, R\ ], 
	Thickness[ .005], 
              Hue[0], Line[{{0, 0}, 
                  R\ s1 {Sin[ph[0]], Cos[ph[0]]}}], 
	
              Hue[ .7], 
              Line[{{0, 0}, 
                  R\ s1 {Sin[ph[2  π]], 
                      Cos[ph[2  π]]}}]}];\)\)

ps  =  
        ParametricPlot[
          R Sqrt[1  -  fu[ψ]^2] {Sin[ph[ψ]], 
              Cos[ph[ψ]]}, {ψ, 2  π, 
            4  π}, 
AspectRatio  ->  Automatic, 
          PlotStyle  -> {{Thickness[ .005], Hue[0.7]}}];

dg  =  
        Graphics[{Circle[{0, 0}, 
              R s1 0.4, {-ph[2  π] + 5  π/2, π/
                  2}],   Circle[{0, 0}, 
              R s1 0.43, {-ph[2  π] + 5  π/2, π/
                  2}]}];
dh  =  
        Graphics[
          Text["Δϕ", 
            R s1 { .4,  .55}]];
Show[pr, ps, dg, dh, 
        PlotRange  ->  
          R {{-1, 1}, {-1, 
                1}},  
Ticks  ->  None,  
        TextStyle  ->  {FontSize  ->  14}];

Trajectory in space

pd  =  
        ParametricPlot3D[ {R Sqrt[1  -  fu[ψ]^2] Sin[ph[ψ]], 
            R Sqrt[1  -  fu[ψ]^2] Cos[ph[ψ]],  
            R fu[ψ], {Thickness[ .007], Hue[0]}}, {ψ, 0, 
            2  π},  BoxRatios  ->  {1, 1, 1}];

pe  =  
        ParametricPlot3D[{R Sqrt[1  -  fu[ψ]^2] Sin[ph[ψ]], 
            R Sqrt[1  -  fu[ψ]^2] Cos[ph[ψ]], 
            R fu[ψ], {Thickness[ .007], Hue[0.7]}}, {ψ, 2  π, 
            4  π},  BoxRatios  ->  {1, 1, 1}];

pf  =  
        ParametricPlot3D[{R Sqrt[1  -  fu[ψ]^2] Sin[ph[ψ]], 
            R Sqrt[1  -  fu[ψ]^2] Cos[ph[ψ]], 
            R  fu[ψ], {Thickness[ .007], Hue[0]}}, {ψ, 2  π, 
            4  π},  BoxRatios  ->  {1, 1, 1}];

k1  =  
        ParametricPlot3D[{R s1 Sin[ϕ], R s1 Cos[ϕ], 
            R u1, {Thickness[ .007]}}, {ϕ, 0, 2  π},  
          BoxRatios  ->  {1, 1, 1},  
          DisplayFunction  ->  Identity];
k2  =  
        ParametricPlot3D[{R s2 Sin[ϕ], R s2 Cos[ϕ], 
            R u2, {Thickness[ .007]}}, {ϕ, 0, 2  π},  
          BoxRatios  ->  {1, 1, 1},  
          DisplayFunction  ->  Identity];
ϕ3  =    π/4  +   .1;
k3  =  
        ParametricPlot3D[{R Sin[θ] Sin[ϕ3], 
            R Sin[θ]  Cos[ϕ3], 
            R Cos[θ], {Thickness[ .007]}}, {θ, θ1, θ2},  BoxRatios  ->  {1, 1, 1},  
          DisplayFunction  ->  Identity];
ϕ4  =    ϕ3  +  π;
k4  =  
        ParametricPlot3D[{R Sin[θ] Sin[ϕ4], 
            R Sin[θ]  Cos[ϕ4], 
            R Cos[θ], {Thickness[ .007]}}, {θ, θ1, θ2},  BoxRatios  ->  {1, 1, 1},  
          DisplayFunction  ->  Identity];

Show[pd, pe, k1, k2, k3, k4, BoxRatios  ->  {1, 1,  .75}, 
        Axes  ->  False, ImageSize  ->  500];

sph  =  
        ParametricPlot3D[
          R {Sin[th] Cos[pp], Sin[th] Sin[pp], 
              Cos[th]}, 
{th, 0, π}, {pp, 0, 
            2  π},  DisplayFunction  ->  Identity];
sn  =  
        Show[sph  /.  Polygon[x___]  ->  {EdgeForm[], Polygon[x]},  
          DisplayFunction  ->  Identity];

Show[pd, pe, sn, k1, k2,  ViewPoint  -> {0, 0.1, -2},  
        ImageSize  ->   500];

Show[pd, pe, k1, k2, sn,  ViewPoint  -> {0, 2,  .2},  
        ImageSize  ->  500];

Show[pd, pe, k1, k2, sn,  ViewPoint  -> {0, 2, -1.8},  
        ImageSize  ->  500];

Large Oscillation 4  (0 < 5° ≤  δ  ≤  130°)

Input Parameters. Choice of latitudes

Resulting Parameters

\(u1\  = \ \ Cos[θ1] // N\ ; \ 
    Print[\*""\!\(u\_1\) = "", u1]\)
\(u2\  = \ \ Cos[θ2] // N; \ \ Print[\*""\!\(u\_2\) = "", 
      u2]\)
\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Print[\*""(Subscript[u, 3]) = "", 
      u3\  //. \ su3\ ]\)

ϵ     =  (ϵ   /. sel) /.  
        su3 ;   Print["ϵ = ", ϵ]
λ     =  (λ  /. sel) /.  
        su3;    Print["λ = ", λ]
s1  =  Sin[θ1] // Chop ;   Print["r1 = ", 
      s1, "  = Relative radius of lower circle of latitude "]
s2  =  Sin[θ2] // N;            Print["r2 = ", 
      s2, "  = Relative radius of higher circle of latitude "]

\(kq\  = \ \((u2 - u1)\)/\((u3 - u1)\)\  /. \ su3; \ 
    Print[\*""m = \!\(k\^2\) = "", 
      kq, \ "  = Parameter of elliptic integrals"]\)

su  =  Sqrt[(u3  -  u1)/2]  /.  su3

tt =  
      4  EllipticK[kq]/
            su /(2 π);     Print["T/T0 = ", 
      tt]

u  =  cos θ   expressed by auxiliary variable  ψ.

fu[ψ_]  =  uu

ph[ψ_] := php[ψ]

Curves for θ(ψ),   (Subscript[ω, 0]) t(ψ)  and  ϕ(ψ).

ppth = Plot[( π-ArcCos[fu[ψ]]) 180/π,{ψ, 0, 6  π},
                      \(AxesLabel\  -> \ {"ψ", "θ [°]"}\), TextStyle  ->  {FontSize  ->   12},
                      Ticks  ->  {PiScale, 
                Automatic},PlotRange  ->  {0, 140},
		         DisplayFunction  ->  Identity];

\(\(ppt\  = \ 
        Plot[{ψ, EllipticF[ψ, kq]/su}, {ψ, 0, 
            6  π}, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle\  -> {GrayLevel[0], {Thickness[ .01], 
                Hue[0]}}, 
		DisplayFunction\  -> \ 
            Identity, \ 
          Ticks\  -> \ 
            PiScale, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel\  -> \ {"ψ", \*""(Subscript[ω, 0])t""}, \ 
          TextStyle\  -> \ {FontSize\  -> \ \ 12}, 
		Epilog\  -> \ {Text[\*""(Subscript[ω, 0])t = ψ"", π {5, 3}], 
              Hue[0], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[
                Subscript[\*""(Subscript[ω, 0])t"", 
                  ex], π {2.5, 5}]}];\)\)

\(\(pph\  = \ 
        Plot[{ψ, ph[ψ]}, {ψ, 0, 
            6  π}, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle\  -> {GrayLevel[0], {Thickness[ .01], 
                Hue[0]}}, 
		\ DisplayFunction\  -> \ Identity, \ \ Ticks\  -> \ 
            PiScale, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel\  -> \ {"ψ", "ϕ"}, \ 
          TextStyle\  -> \ {FontSize\  -> \ \ 12}, 
		Epilog\  -> \ {Text["ϕ = ψ", π {5, 3}], Hue[
                0], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[
                Subscript["ϕ(ψ)", 
                  ex], π {2.5, 
                    5}], 
				\ Hue[ .6], 
              Line[π {{2.15, 2.29}, {3, 1.6}}], \ 
              Text["Δϕ", π {3.4, 
                    1.5}], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Thickness[ .01], 
              Line[{π {2, 2}, {2  π, ph[2  π]}}]}];\)\)

Show[GraphicsArray[{{ppth}, {ppt, pph}}],  
        ImageSize  ->  500];

Projection of Trajectory on  x,y-plane

R  =   .1;

\(\(pr\  = \ 
        ParametricPlot[
          R\ Sqrt[1\  - \ fu[ψ]^2] {Sin[ph[ψ]], 
              Cos[ph[ψ]]}, {ψ, 0, 
            2  π}, 
\ \ \ \ \ \ \ \ \ \ \ AspectRatio\  -> \ Automatic, 
          ImageSize\  -> \ 
            450, 
\ \ \ \ \ \ \ \ \ \ \ \ PlotStyle\  -> {{Thickness[ .005], 
                Hue[0]}}, \ 
	\ Epilog\  -> \ {Circle[{0, 0}, R\ s1], Circle[{0, 0}, R\ s2], 
              Circle[{0, 0}, R\ ], 
	Thickness[ .005], 
              Hue[0], Line[{{0, 0}, 
                  R\ s1 {Sin[ph[0]], Cos[ph[0]]}}], 
	
              Hue[ .7], 
              Line[{{0, 0}, 
                  R\ s1 {Sin[ph[2  π]], 
                      Cos[ph[2  π]]}}]}];\)\)

ps  =  
        ParametricPlot[
          R Sqrt[1  -  fu[ψ]^2] {Sin[ph[ψ]], 
              Cos[ph[ψ]]}, {ψ, 2  π, 
            4  π}, 
AspectRatio  ->  Automatic, 
          PlotStyle  -> {{Thickness[ .005], Hue[0.7]}}];

dg  =  
        Graphics[{Circle[{0, 0}, 
              R s1 0.4, {-ph[2  π] + 5  π/2, π/
                  2}],   Circle[{0, 0}, 
              R s1 0.43, {-ph[2  π] + 5  π/2, π/
                  2}]}];
dh  =  
        Graphics[
          Text["Δϕ", 
            R s1 { .4,  .55}]];
Show[pr, ps, dg, dh, 
        PlotRange  ->  
          R {{-1, 1}, {-1, 
                1}},  
Ticks  ->  None,  
        TextStyle  ->  {FontSize  ->  14}];

Trajectory in space

pd  =  
        ParametricPlot3D[ {R Sqrt[1  -  fu[ψ]^2] Sin[ph[ψ]], 
            R Sqrt[1  -  fu[ψ]^2] Cos[ph[ψ]],  
            R fu[ψ], {Thickness[ .007], Hue[0]}}, {ψ, 0, 
            2  π},  BoxRatios  ->  {1, 1, 1}];

pe  =  
        ParametricPlot3D[{R Sqrt[1  -  fu[ψ]^2] Sin[ph[ψ]], 
            R Sqrt[1  -  fu[ψ]^2] Cos[ph[ψ]], 
            R fu[ψ], {Thickness[ .007], Hue[0.7]}}, {ψ, 2  π, 
            4  π},  BoxRatios  ->  {1, 1, 1}];

pf  =  
        ParametricPlot3D[{R Sqrt[1  -  fu[ψ]^2] Sin[ph[ψ]], 
            R Sqrt[1  -  fu[ψ]^2] Cos[ph[ψ]], 
            R  fu[ψ], {Thickness[ .007], Hue[0]}}, {ψ, 2  π, 
            4  π},  BoxRatios  ->  {1, 1, 1}];

k1  =  
        ParametricPlot3D[{R s1 Sin[ϕ], R s1 Cos[ϕ], 
            R u1, {Thickness[ .007]}}, {ϕ, 0, 2  π},  
          BoxRatios  ->  {1, 1, 1},  
          DisplayFunction  ->  Identity];
k2  =  
        ParametricPlot3D[{R s2 Sin[ϕ], R s2 Cos[ϕ], 
            R u2, {Thickness[ .007]}}, {ϕ, 0, 2  π},  
          BoxRatios  ->  {1, 1, 1},  
          DisplayFunction  ->  Identity];
ϕ3  =    π/4  +   .1;
k3  =  
        ParametricPlot3D[{R Sin[θ] Sin[ϕ3], 
            R Sin[θ]  Cos[ϕ3], 
            R Cos[θ], {Thickness[ .007]}}, {θ, θ1, θ2},  BoxRatios  ->  {1, 1, 1},  
          DisplayFunction  ->  Identity];
ϕ4  =    ϕ3  +  π;
k4  =  
        ParametricPlot3D[{R Sin[θ] Sin[ϕ4], 
            R Sin[θ]  Cos[ϕ4], 
            R Cos[θ], {Thickness[ .007]}}, {θ, θ1, θ2},  BoxRatios  ->  {1, 1, 1},  
          DisplayFunction  ->  Identity];

Show[pd, pe, k1, k2, k3, k4, BoxRatios  ->  {1, 1,  .75}, 
        Axes  ->  False, ImageSize  ->  500];

sph  =  
        ParametricPlot3D[
          R {Sin[th] Cos[pp], Sin[th] Sin[pp], 
              Cos[th]}, 
{th, 0, π}, {pp, 0, 
            2  π},  DisplayFunction  ->  Identity];
sn  =  
        Show[sph  /.  Polygon[x___]  ->  {EdgeForm[], Polygon[x]},  
          DisplayFunction  ->  Identity];

Show[pd, pe, sn, k1, k2,  ViewPoint  -> {0, 0.1, -2},  
        ImageSize  ->   500];

Show[pd, pe, k1, k2, sn,  ViewPoint  -> {0, 2,  .2},  
        ImageSize  ->  500];

Show[pd, pe, k1, k2, sn,  ViewPoint  -> {0, 2, -1.8},  
        ImageSize  ->  500];

Regular Precession: Motion on a Circle of Latitude  (0 <  δ  =  70°)

Input Parameters. Choice of latitudes

Resulting Parameters

\(u1\  = \ \ Cos[θ1] // N\ ; \ 
    Print[\*""\!\(u\_1\) = "", u1]\)
\(u2\  = \ \ Cos[θ2] // N; \ \ Print[\*""\!\(u\_2\) = "", 
      u2]\)
\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Print[\*""(Subscript[u, 3]) = "", 
      u3\  //. \ su3\ ]\)

ϵ     =  (ϵ   /. sel) /.  
        su3 ;   Print["ϵ = ", ϵ]
λ     =  (λ  /. sel) /.  
        su3;    Print["λ = ", λ]
s1  =  Sin[θ1] // Chop ;   Print["r1 = ", 
      s1, "  = Relative radius of lower circle of latitude "]
s2  =  Sin[θ2] // N;            Print["r2 = ", 
      s2, "  = Relative radius of higher circle of latitude "]

\(kq\  = \ \((u2 - u1)\)/\((u3 - u1)\)\  /. \ su3; \ 
    Print[\*""m = \!\(k\^2\) = "", 
      kq, \ "  = Parameter of elliptic integrals"]\)

su  =  Sqrt[(u3  -  u1)/2]  /.  su3

tt =  
      4  EllipticK[kq]/
            su /(2 π);     Print["T/T0 = ", 
      tt]

u  =  cos θ   expressed by auxiliary variable  ψ.

fu[ψ_]  =  uu

ph[ψ_] := php[ψ]

Curves for θ(ψ),   (Subscript[ω, 0]) t(ψ)  and  ϕ(ψ).

ppth = Plot[( π-ArcCos[fu[ψ]]) 180/π,{ψ, 0, 6  π},
                      \(AxesLabel\  -> \ {"ψ", "θ [°]"}\), TextStyle  ->  {FontSize  ->   12},
                      Ticks  ->  {PiScale, 
                Automatic},PlotRange  ->  {0, 140},
		         DisplayFunction  ->  Identity];

\(\(ppt\  = \ 
        Plot[{ψ, EllipticF[ψ, kq]/su}, {ψ, 0, 
            6  π}, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle\  -> {GrayLevel[0], {Thickness[ .01], 
                Hue[0]}}, 
		DisplayFunction\  -> \ 
            Identity, \ 
          Ticks\  -> \ 
            PiScale, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel\  -> \ {"ψ", \*""(Subscript[ω, 0])t""}, \ 
          TextStyle\  -> \ {FontSize\  -> \ \ 12}, 
		Epilog\  -> \ {Text[\*""(Subscript[ω, 0])t = ψ"", π {5, 3}], 
              Hue[0], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[
                Subscript[\*""(Subscript[ω, 0])t"", 
                  ex], π {2.5, 5}]}];\)\)

\(\(pph\  = \ 
        Plot[{ψ, ph[ψ]}, {ψ, 0, 
            6  π}, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle\  -> {GrayLevel[0], {Thickness[ .01], 
                Hue[0]}}, 
		\ DisplayFunction\  -> \ Identity, \ \ Ticks\  -> \ 
            PiScale, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel\  -> \ {"ψ", "ϕ"}, \ 
          TextStyle\  -> \ {FontSize\  -> \ \ 12}, 
		Epilog\  -> \ {Text["ϕ = ψ", π {5, 3}], Hue[
                0], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Text[
                Subscript["ϕ(ψ)", 
                  ex], π {2.5, 
                    5}], 
				\ Hue[ .6], 
              Line[π {{2.15, 2.29}, {3, 1.6}}], \ 
              Text["Δϕ", π {3.4, 
                    1.5}], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Thickness[ .01], 
              Line[{π {2, 2}, {2  π, ph[2  π]}}]}];\)\)

Show[GraphicsArray[{{ppth}, {ppt, pph}}],  
        ImageSize  ->  500];

Projection of Trajectory on  x,y-plane

R  =   .1;

\(\(pr\  = \ 
        ParametricPlot[
          R\ Sqrt[1\  - \ fu[ψ]^2] {Sin[ph[ψ]], 
              Cos[ph[ψ]]}, {ψ, 0, 
            2  π}, 
\ \ \ \ \ \ \ \ \ \ \ AspectRatio\  -> \ Automatic, 
          ImageSize\  -> \ 
            450, 
\ \ \ \ \ \ \ \ \ \ \ \ PlotStyle\  -> {{Thickness[ .005], 
                Hue[0]}}, \ 
	\ Epilog\  -> \ {Circle[{0, 0}, R\ s1], Circle[{0, 0}, R\ s2], 
              Circle[{0, 0}, R\ ], 
	Thickness[ .005], 
              Hue[0], Line[{{0, 0}, 
                  R\ s1 {Sin[ph[0]], Cos[ph[0]]}}], 
	
              Hue[ .7], 
              Line[{{0, 0}, 
                  R\ s1 {Sin[ph[2  π]], 
                      Cos[ph[2  π]]}}]}];\)\)

dg  =  
        Graphics[{Circle[{0, 0}, 
              R s1 0.4, {-ph[2  π] + 5  π/2, π/
                  2}],   Circle[{0, 0}, 
              R s1 0.43, {-ph[2  π] + 5  π/2, π/
                  2}]}];
dh  =  
        Graphics[
          Text["Δϕ", 
            R s1 { .4,  .55}]];
Show[pr, dg, dh, 
        PlotRange  ->  
          R {{-1, 1}, {-1, 
                1}},  
Ticks  ->  None,  
        TextStyle  ->  {FontSize  ->  14}];

Trajectory in space

pd  =  
        ParametricPlot3D[ {R Sqrt[1  -  fu[ψ]^2] Sin[ph[ψ]], 
            R Sqrt[1  -  fu[ψ]^2] Cos[ph[ψ]],  
            R fu[ψ], {Thickness[ .007], Hue[0]}}, {ψ, 0, 
            2  π},  BoxRatios  ->  {1, 1, 1}];

Show[pd, BoxRatios  ->  {1, 1,  .75}, Axes  ->  False, 
        ImageSize  ->  500];

sph  =  
        ParametricPlot3D[
          R {Sin[th] Cos[pp], Sin[th] Sin[pp], 
              Cos[th]}, 
{th, 0, π}, {pp, 0, 
            2  π},  DisplayFunction  ->  Identity];
sn  =  
        Show[sph  /.  Polygon[x___]  ->  {EdgeForm[], Polygon[x]},  
          DisplayFunction  ->  Identity];

Show[pd, sn,  ViewPoint  -> {0, 0.1, -2},  
        ImageSize  ->   500];

Show[pd, sn,  ViewPoint  -> {0, 2,  .2},  
        ImageSize  ->  500];

Show[pd, sn,  ViewPoint  -> {0, 2, -1.8},  
        ImageSize  ->  500];