doc(Kepler problem): #83 add mathematical-physical descriptions (#93)

* doc(kepler): #83 tau_of_u (elliptic)

* chore(test): #83 more strict test

* chore(markdown): #83 refactor

* chore(markdown): #83 fix yml

* chore: #83 fix equations

* doc(kepler): #83 tau_of_u (elliptic)

* doc(kepler): #83 the parabolic case

* doc(kepler): #83 the hyperbolic case

* chore: #83 fixes

* chore: #83 fixes

* doc(kepler): #83 derivatives

* doc(kepler): #83 hyperbolic approximations

* doc(kepler): #83 model

* doc(kepler): #83 finish

docs/references/maths/trigonometrics.md

@@ -0,0 +1,3 @@
+# Trigonometrics
+
+::: hamilflow.maths.trigonometrics

docs/references/models/kepler_problem/dynamics.md

@@ -0,0 +1,3 @@
+# Dynamics
+
+::: hamilflow.models.kepler_problem.dynamics

docs/references/models/kepler_problem/model.md

@@ -0,0 +1,3 @@
+# Model
+
+::: hamilflow.models.kepler_problem.model

docs/references/models/kepler_problem/numerics.md

@@ -0,0 +1,3 @@
+# Numerics
+
+::: hamilflow.models.kepler_problem.numerics

hamilflow/models/kepler_problem/__init__.py

@@ -1,5 +1,5 @@
 """Package for the Kepler problem."""
 
-from .model import Kepler2D, Kepler2DIoM, Kepler2DSystem
+from .model import Kepler2D, Kepler2DFI, Kepler2DSystem
 
-__all__ = ["Kepler2D", "Kepler2DIoM", "Kepler2DSystem"]
+__all__ = ["Kepler2D", "Kepler2DFI", "Kepler2DSystem"]

hamilflow/models/kepler_problem/dynamics.py

@@ -12,28 +12,55 @@
 _1_3 = 1 / 3
 
 
-def _tau_of_u_exact_elliptic(
+def tau_of_u_exact_elliptic(
     ecc: float,
     u: "npt.NDArray[np.float64]",
 ) -> "npt.NDArray[np.float64]":
+    r"""Exact solution for tau of u in the elliptic case.
+
+    For $-e \le u \le e$,
+    $$ \tau = -\frac{\sqrt{e^2-u^2}}{(1-e^2)(1+u)}
+    + \frac{\frac{\pi}{2} - \arctan\frac{e^2+u}{\sqrt{(1-e^2)(e^2-u^2)}}}{(1-e^2)^{\frac{3}{2}}}\,. $$
+    """
     cosqr, eusqrt = 1 - ecc**2, np.sqrt(ecc**2 - u**2)
     trig_numer = np.pi / 2 - np.arctan((ecc**2 + u) / np.sqrt(cosqr) / eusqrt)
     return -eusqrt / cosqr / (1 + u) + trig_numer / cosqr**1.5
 
 
-def _tau_of_e_plus_u_elliptic(
+def tau_of_e_plus_u_elliptic(
     ecc: float,
     u: "npt.NDArray[np.float64]",
 ) -> "npt.NDArray[np.float64]":
+    r"""Expansion for tau of u in the elliptic case at $u = -e+0$.
+
+    The exact solution has a removable singularity at $u = -e$, hence this
+    expansion helps with numerics.
+
+    Let $\epsilon = \sqrt{\frac{2(e+u)}{e}}$,
+    $$ \tau = \frac{\pi}{(1-e^2)^\frac{3}{2}}
+    - \frac{1}{(1-e)^2}\epsilon
+    - \frac{1-9e}{24(1-e)^3}\epsilon^3
+    + O\left(\epsilon^5\right)\,. $$
+    """
     epu = np.sqrt(2 * (ecc + u) / ecc)
     const = np.pi / (1 - ecc**2) ** 1.5
-    return const - epu / (ecc - 1) ** 2 + epu**3 * (1 - 9 * ecc) / 24 / (1 - ecc) ** 3
+    return const - epu / (ecc - 1) ** 2 - epu**3 * (1 - 9 * ecc) / 24 / (1 - ecc) ** 3
 
 
-def _tau_of_e_minus_u_elliptic(
+def tau_of_e_minus_u_elliptic(
     ecc: float,
     u: "npt.NDArray[np.float64]",
 ) -> "npt.NDArray[np.float64]":
+    r"""Expansion for tau of u in the elliptic case at $u = +e-0$.
+
+    The exact solution has a removable singularity at $u = +e$, hence this
+    expansion helps with numerics.
+
+    Let $\epsilon = \sqrt{\frac{2(e-u)}{e}}$,
+    $$ \tau = \frac{1}{(1+e)^2}\epsilon
+    - \frac{1+9e}{24(1+e)^3}\epsilon^3
+    + O\left(\epsilon^5\right)\,. $$
+    """
     emu = np.sqrt(2 * (ecc - u) / ecc)
     return emu / (1 + ecc) ** 2 - emu**3 * (1 + 9 * ecc) / 24 / (1 + ecc) ** 3
 
@@ -48,14 +75,17 @@ def tau_of_u_elliptic(ecc: float, u: "npt.ArrayLike") -> "npt.NDArray[np.float64
     return _approximate_at_termina(
         ecc,
         u,
-        _tau_of_u_exact_elliptic,
-        _tau_of_e_plus_u_elliptic,
-        _tau_of_e_minus_u_elliptic,
+        tau_of_u_exact_elliptic,
+        tau_of_e_plus_u_elliptic,
+        tau_of_e_minus_u_elliptic,
     )
 
 
 def tau_of_u_parabolic(ecc: float, u: "npt.ArrayLike") -> "npt.NDArray[np.float64]":
-    """Calculate the scaled time tau from u in the parabolic case.
+    r"""Calculate the scaled time tau from u in the parabolic case.
+
+    For $-1 < u \le 1$,
+    $$ \tau = \frac{\sqrt{1-u}(2+u)}{3(1+u)^\frac{3}{2}}\,. $$
 
     :param ecc: eccentricity, ecc == 1 (unchecked, unused)
     :param u: convenient radial inverse
@@ -65,31 +95,57 @@ def tau_of_u_parabolic(ecc: float, u: "npt.ArrayLike") -> "npt.NDArray[np.float6
     return np.sqrt(1 - u) * (2 + u) / 3 / (1 + u) ** 1.5
 
 
-def _tau_of_u_exact_hyperbolic(
+def tau_of_u_exact_hyperbolic(
     ecc: float,
     u: "npt.ArrayLike",
 ) -> "npt.NDArray[np.float64]":
+    r"""Exact solution for tau of u in the hyperbolic case.
+
+    For $-1 < u \le e$,
+    $$ \tau = \frac{\sqrt{e^2-u^2}}{(e^2-1)(1+u)}
+    - \frac{\mathrm{artanh}\frac{e^2+u}{\sqrt{(e^2-1)(e^2-u^2)}}}{(e^2-1)^{\frac{3}{2}}}\,. $$
+    """
     u = np.asarray(u)
     cosqr, eusqrt = ecc**2 - 1, np.sqrt(ecc**2 - u**2)
     trig_numer = np.arctanh(np.sqrt(cosqr) * eusqrt / (ecc**2 + u))
     return eusqrt / cosqr / (1 + u) - trig_numer / cosqr**1.5
 
 
-def _tau_of_1_plus_u_hyperbolic(
+def tau_of_1_plus_u_hyperbolic(
     ecc: float,
     u: "npt.NDArray[np.float64]",
 ) -> "npt.NDArray[np.float64]":
+    r"""Expansion for tau of u in the hyperbolic case at $u = -1+0$.
+
+    The exact solution has a logarithmic singularity at $u = -1$, hence this
+    expansion helps with numerics.
+
+    Let $\epsilon = \frac{e(1+u)}{2(e^2-1)}$,
+    $$ \tau = \ln \epsilon + \frac{e}{2\epsilon}
+    + 1 - \frac{e^2+2}{e}\epsilon + \left(3+\frac{2}{e^2}\right)\epsilon^2
+    + O(\epsilon^3)\,. $$
+    """
     cosqr = ecc**2 - 1
     up1 = ecc * (1 + u) / 2 / cosqr
     diverge = np.log(up1) + ecc / 2 / up1
     regular = 1 - (2 + ecc**2) / ecc * up1 + (3 + 2 / ecc**2) * up1**2
     return (diverge + regular) / cosqr**1.5
 
 
-def _tau_of_e_minus_u_hyperbolic(
+def tau_of_e_minus_u_hyperbolic(
     ecc: float,
     u: "npt.NDArray[np.float64]",
 ) -> "npt.NDArray[np.float64]":
+    r"""Expansion for tau of u in the elliptic case at $u = +e-0$.
+
+    The exact solution has a removable singularity at $u = +e$, hence this
+    expansion helps with numerics.
+
+    Let $\epsilon = \sqrt{\frac{2(e-u)}{e}}$,
+    $$ \tau = \frac{1}{(1+e)^2}\epsilon
+    + \frac{1+9e}{24(1+e)^3}\epsilon^3
+    + O\left(\epsilon^5\right)\,. $$
+    """
     emu = np.sqrt(2 * (ecc - u) / ecc)
     return emu / (1 + ecc) ** 2 + emu**3 * (1 + 9 * ecc) / 24 / (1 + ecc) ** 3
 
@@ -104,14 +160,16 @@ def tau_of_u_hyperbolic(ecc: float, u: "npt.ArrayLike") -> "npt.NDArray[np.float
     return _approximate_at_termina(
         ecc,
         u,
-        _tau_of_u_exact_hyperbolic,
-        _tau_of_1_plus_u_hyperbolic,
-        _tau_of_e_minus_u_hyperbolic,
+        tau_of_u_exact_hyperbolic,
+        tau_of_1_plus_u_hyperbolic,
+        tau_of_e_minus_u_hyperbolic,
     )
 
 
 def tau_of_u_prime(ecc: float, u: "npt.ArrayLike") -> "npt.NDArray[np.float64]":
-    """Calculate the first derivative of scaled time tau with respect to u.
+    r"""Calculate the first derivative of scaled time tau with respect to u.
+
+    $$ \tau'(u) = -\frac{1}{(1+u)^2\sqrt{e^2-u^2}}\,. $$
 
     :param ecc: eccentricity, ecc >= 0 (unchecked)
     :param u: convenient radial inverse
@@ -122,7 +180,9 @@ def tau_of_u_prime(ecc: float, u: "npt.ArrayLike") -> "npt.NDArray[np.float64]":
 
 
 def tau_of_u_prime2(ecc: float, u: "npt.ArrayLike") -> "npt.NDArray[np.float64]":
-    """Calculate the second derivative of scaled time tau with respect to u.
+    r"""Calculate the second derivative of scaled time tau with respect to u.
+
+    $$ \tau''(u) = \frac{2e^2 - u - 3u^2}{(1+u)^3 (e^2-u^2)^\frac{3}{2}} $$
 
     :param ecc: eccentricity, ecc >= 0 (unchecked)
     :param u: convenient radial inverse
@@ -137,7 +197,12 @@ def esolve_u_from_tau_parabolic(
     ecc: float,
     tau: "npt.ArrayLike",
 ) -> "npt.NDArray[np.float64]":
-    """Calculate the convenient radial inverse u from tau in the parabolic case, using the exact solution.
+    r"""Calculate the convenient radial inverse u from tau in the parabolic case, using the exact solution.
+
+    Let $\kappa = 1+9\tau^2$,
+    $$ u = -1
+    + \left(\frac{3\tau}{\kappa^\frac{3}{2}} + \frac{1}{\kappa}\right)^\frac{1}{3}
+    + \left(\frac{1}{3\tau \kappa^\frac{3}{2}+\kappa^2}\right)^\frac{1}{3}\,. $$
 
     :param ecc: eccentricity, ecc = 0 (unchecked, unused)
     :param tau: scaled time