Skip to content

Numerics¤

Numerics for the Kepler problem.

nsolve_u_from_tau_bisect(ecc, tau) ¤

Calculate the convenient radial inverse u from tau in the elliptic or parabolic case, using the bisect method.

Parameters:

Name Type Description Default
ecc float

eccentricity, ecc > 0, ecc != 1

 required
tau ArrayLike

scaled time

 required

Returns:

Type Description
list[OptimizeResult]

numeric OptimizeResult from scipy
Source code in hamilflow/models/kepler_problem/numerics.py 
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
def nsolve_u_from_tau_bisect(ecc: float, tau: "npt.ArrayLike") -> list[OptimizeResult]:
    """Calculate the convenient radial inverse u from tau in the elliptic or parabolic case, using the bisect method.

    :param ecc: eccentricity, ecc > 0, ecc != 1
    :param tau: scaled time
    :return: numeric OptimizeResult from scipy
    """
    tau_s = np.asarray(tau).reshape(-1)
    if 0 < ecc < 1:
        tau_of_u = tau_of_u_elliptic
    elif ecc > 1:
        tau_of_u = tau_of_u_hyperbolic
    else:
        msg = f"Expect ecc > 0, ecc != 1, got {ecc}"
        raise ValueError(msg)

    def f(u: float, tau: float) -> np.float64:
        return (
            np.finfo(np.float64).max if u == -1 else np.float64(tau_of_u(ecc, u) - tau)
        )

    return [toms748(f, max(-1, -ecc), ecc, (ta,), 2, full_output=True) for ta in tau_s]

nsolve_u_from_tau_newton(ecc, tau) ¤

Calculate the convenient radial inverse u from tau in the elliptic or parabolic case, using the Newton method.

Parameters:

Name Type Description Default
ecc float

eccentricity, ecc > 0, ecc != 1

 required
tau ArrayLike

scaled time

 required

Returns:

Type Description
OptimizeResult

numeric OptimizeResult from scipy

Raises:

Type Description
ValueError

when ecc is invalid
Source code in hamilflow/models/kepler_problem/numerics.py 
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
def nsolve_u_from_tau_newton(ecc: float, tau: "npt.ArrayLike") -> OptimizeResult:
    """Calculate the convenient radial inverse u from tau in the elliptic or parabolic case, using the Newton method.

    :param ecc: eccentricity, ecc > 0, ecc != 1
    :param tau: scaled time
    :raises ValueError: when `ecc` is invalid
    :return: numeric OptimizeResult from scipy
    """
    tau = np.asarray(tau)
    if 0 < ecc < 1:
        tau_of_u = tau_of_u_elliptic
        u0 = _u0_elliptic(ecc, tau)
    elif ecc > 1:
        tau_of_u = tau_of_u_hyperbolic
        u0 = _u0_hyperbolic(ecc, tau)
    else:
        msg = f"Expect ecc > 0, ecc != 1, got {ecc}"
        raise ValueError(msg)

    def f(u: float, tau: "npt.NDArray[np.float64]") -> "npt.NDArray[np.float64]":
        return tau_of_u(ecc, u) - tau

    def fprime(
        u: float,
        tau: "npt.ArrayLike",  # noqa: ARG001
    ) -> "npt.NDArray[np.float64]":
        return tau_of_u_prime(ecc, u)

    def fprime2(
        u: float,
        tau: "npt.ArrayLike",  # noqa: ARG001
    ) -> "npt.NDArray[np.float64]":
        return tau_of_u_prime2(ecc, u)

    return newton(f, u0, fprime, (tau,), fprime2=fprime2, full_output=True, disp=False)

u_of_tau(ecc, tau, method='bisect') ¤

Calculate the convenient radial inverse u from tau, using numeric methods.

Parameters:

Name Type Description Default
ecc float

eccentricity, ecc >= 0

 required
tau ArrayLike

scaled time

 required
method Literal['bisect', 'newton']

"newton" or "bisect" numeric methods, defaults to "bisect"

 'bisect'

Returns:

Type Description
NDArray[float64]

convenient radial inverse u

Raises:

Type Description
ValueError

when method is invalid
ValueError

when ecc is invalid
Source code in hamilflow/models/kepler_problem/numerics.py 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
def u_of_tau(
    ecc: float,
    tau: "npt.ArrayLike",
    method: Literal["bisect", "newton"] = "bisect",
) -> "npt.NDArray[np.float64]":
    """Calculate the convenient radial inverse u from tau, using numeric methods.

    :param ecc: eccentricity, ecc >= 0
    :param tau: scaled time
    :param method: "newton" or "bisect" numeric methods, defaults to "bisect"
    :raises ValueError: when `method` is invalid
    :raises ValueError: when `ecc` is invalid
    :return: convenient radial inverse u
    """
    tau = np.asarray(tau)
    if ecc == 0:
        return np.zeros(tau.shape)
    elif ecc == 1:
        return esolve_u_from_tau_parabolic(ecc, tau)
    elif ecc > 0:
        match method:
            case "bisect":
                return np.array([s[0] for s in nsolve_u_from_tau_bisect(ecc, tau)])
            case "newton":
                return nsolve_u_from_tau_newton(ecc, tau)[0]
            case _:
                msg = f"Expect 'bisect' or 'newton', got {method}"
                raise ValueError(msg)
    else:
        msg = f"Expect ecc >= 0, got {ecc}"
        raise ValueError(msg)