11277630968

Committed 10 Oct 2024 03:54PM UTC coverage: 73.794% (-0.2%) from 73.99%

Build # 11277630968

Build Type

Pull #166

github

Committed by

web-flow

Commit Message

Merge 4ac6449af into 3172ac203

Pull Request Pull Request #166: Swap to pdm

Run Details

213 of 288 new or added lines in 37 files covered. (73.96%)

2768 of 3751 relevant lines covered (73.79%)

0.74 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

70.49

/src/thor/orbits/gauss.py

import numpy as np
from adam_core.constants import Constants as c
from adam_core.coordinates import (
    CartesianCoordinates,
    Origin,
    SphericalCoordinates,
    transform_coordinates,
)
from adam_core.orbits import Orbits
from adam_core.time import Timestamp
from numba import jit
from numpy import roots

from .gibbs import calcGibbs
from .herrick_gibbs import calcHerrickGibbs
from .iterators import iterateStateTransition

__all__ = [
    "_calcV",
    "_calcA",
    "_calcB",
    "_calcLambdas",
    "_calcRhos",
    "approxLangrangeCoeffs",
    "calcGauss",
    "gaussIOD",
]

MU = c.MU
C = c.C


@jit("f8(f8[:], f8[:], f8[:])", nopython=True, cache=True)
def _calcV(rho1_hat, rho2_hat, rho3_hat):
    # Vector triple product that gives the area of
    # the "volume of the parallelepiped" or according to
    # to Milani et al. 2008: 3x volume of the pyramid with vertices q, r1, r2, r3.
    # Note that vector triple product rules apply here.
    return np.dot(np.cross(rho1_hat, rho2_hat), rho3_hat)


@jit("f8(f8[:], f8[:], f8[:], f8[:], f8[:], f8, f8, f8)", nopython=True, cache=True)
def _calcA(q1, q2, q3, rho1_hat, rho3_hat, t31, t32, t21):
    # Equation 21 from Milani et al. 2008
    return np.linalg.norm(q2) ** 3 * np.dot(np.cross(rho1_hat, rho3_hat), (t32 * q1 - t31 * q2 + t21 * q3))


@jit("f8(f8[:], f8[:], f8[:], f8[:], f8, f8, f8, f8)", nopython=True, cache=True)
def _calcB(q1, q3, rho1_hat, rho3_hat, t31, t32, t21, mu=MU):
    # Equation 19 from Milani et al. 2008
    return mu / 6 * t32 * t21 * np.dot(np.cross(rho1_hat, rho3_hat), ((t31 + t32) * q1 + (t31 + t21) * q3))


@jit("UniTuple(f8, 2)(f8, f8, f8, f8, f8)", nopython=True, cache=True)
def _calcLambdas(r2_mag, t31, t32, t21, mu=MU):
    # Equations 16 and 17 from Milani et al. 2008
    lambda1 = t32 / t31 * (1 + mu / (6 * r2_mag**3) * (t31**2 - t32**2))
    lambda3 = t21 / t31 * (1 + mu / (6 * r2_mag**3) * (t31**2 - t21**2))
    return lambda1, lambda3


@jit(
    "UniTuple(f8[:], 3)(f8, f8, f8[:], f8[:], f8[:], f8[:], f8[:], f8[:], f8)",
    nopython=True,
    cache=True,
)
def _calcRhos(lambda1, lambda3, q1, q2, q3, rho1_hat, rho2_hat, rho3_hat, V):
    # This can be derived by taking a series of scalar products of the coplanarity condition equation
    # with cross products of unit vectors in the direction of the observer, in particular, see Chapter 9 in
    # Milani's book on the theory of orbit determination
    numerator = -lambda1 * q1 + q2 - lambda3 * q3
    rho1_mag = np.dot(numerator, np.cross(rho2_hat, rho3_hat)) / (lambda1 * V)
    rho2_mag = np.dot(numerator, np.cross(rho1_hat, rho3_hat)) / V
    rho3_mag = np.dot(numerator, np.cross(rho1_hat, rho2_hat)) / (lambda3 * V)
    rho1 = rho1_mag * rho1_hat
    rho2 = rho2_mag * rho2_hat
    rho3 = rho3_mag * rho3_hat
    return rho1, rho2, rho3


@jit("UniTuple(f8, 2)(f8, f8, f8)", nopython=True, cache=True)
def approxLangrangeCoeffs(r_mag, dt, mu=MU):
    # Calculate the Lagrange coefficients (Gauss f and g series)
    f = 1 - (1 / 2) * mu / r_mag**3 * dt**2
    g = dt - (1 / 6) * mu / r_mag**3 * dt**2
    return f, g


def calcGauss(r1, r2, r3, t1, t2, t3):
    """
    Calculates the velocity vector at the location of the second position vector (r2) with Gauss's
    original method.

    .. math::
        f_1 \approx 1 - \frac{1}{2}\frac{\mu}{r_2^3} (t_1 - t_2)^2

        f_3 \approx 1 - \frac{1}{2}\frac{\mu}{r_2^3} (t_3 - t_2)^2

        g_1 \approx (t_1 - t_2) - \frac{1}{6}\frac{\mu}{r_2^3} (t_1 - t_2)^2

        g_3 \approx (t_3 - t_2) - \frac{1}{6}\frac{\mu}{r_2^3} (t_3 - t_2)^2

        \vec{v}_2 = \frac{1}{f_1 g_3 - f_3 g_1} (-f_3 \vec{r}_1 + f_1 \vec{r}_3)


    For more details on theory see Chapter 5 in Howard D. Curtis' "Orbital Mechanics for
    Engineering Students".

    Parameters
    ----------
    r1 : `~numpy.ndarray` (3)
        Heliocentric position vector at time 1 in cartesian coordinates in units
        of AU.
    r2 : `~numpy.ndarray` (3)
        Heliocentric position vector at time 2 in cartesian coordinates in units
        of AU.
    r3 : `~numpy.ndarray` (3)
        Heliocentric position vector at time 3 in cartesian coordinates in units
        of AU.
    t1 : float
        Time at r1. Units of MJD or JD work or any decimal time format (one day is 1.00) as
        long as all times are given in the same format.
    t2 : float
        Time at r2. Units of MJD or JD work or any decimal time format (one day is 1.00) as
        long as all times are given in the same format.
    t3 : float
        Time at r3. Units of MJD or JD work or any decimal time format (one day is 1.00) as
        long as all times are given in the same format.

    Returns
    -------
    v2 : `~numpy.ndarray` (3)
        Velocity of object at position r2 at time t2 in units of AU per day.
    """
    t12 = t1 - t2
    t32 = t3 - t2
    r2_mag = np.linalg.norm(r2)
    f1, g1 = approxLangrangeCoeffs(r2_mag, t12)
    f3, g3 = approxLangrangeCoeffs(r2_mag, t32)
    v2 = (1 / (f1 * g3 - f3 * g1)) * (-f3 * r1 + f1 * r3)
    return v2


def gaussIOD(
    coords,
    observation_times,
    coords_obs,
    velocity_method="gibbs",
    light_time=True,
    iterate=True,
    iterator="state transition",
    mu=MU,
    max_iter=10,
    tol=1e-15,
):
    """
    Compute up to three intial orbits using three observations in angular equatorial
    coordinates.

    Parameters
    ----------
    coords : `~numpy.ndarray` (3, 2)
        RA and Dec of three observations in units of degrees.
    observation_times : `~numpy.ndarray` (3)
        Times of the three observations in units of decimal days (MJD or JD for example).
    coords_obs : `~numpy.ndarray` (3, 2)
        Heliocentric position vector of the observer at times t in units of AU.
    velocity_method : {'gauss', gibbs', 'herrick+gibbs'}, optional
        Which method to use for calculating the velocity at the second observation.
        [Default = 'gibbs']
    light_time : bool, optional
        Correct for light travel time.
        [Default = True]
    iterate : bool, optional
        Iterate initial orbit using universal anomaly to better approximate the
        Lagrange coefficients.
    mu : float, optional
        Gravitational parameter (GM) of the attracting body in units of
        AU**3 / d**2.
    max_iter : int, optional
        Maximum number of iterations over which to converge to solution.
    tol : float, optional
        Numerical tolerance to which to compute chi using the Newtown-Raphson
        method.

    Returns
    -------
    epochs : `~numpy.ndarry` (<3)
        Epochs in units of decimal days at which preliminary orbits are
        defined.
    orbits : `~numpy.ndarray` ((<3, 6) or (0))
        Up to three preliminary orbits (as cartesian state vectors).
    """
    coords = transform_coordinates(
        SphericalCoordinates.from_kwargs(
            lon=coords[:, 0],
            lat=coords[:, 1],
            origin=Origin.from_kwargs(code=np.full(len(coords), "Unknown", dtype="object")),
            frame="equatorial",
        ).to_unit_sphere(),
        representation_out=CartesianCoordinates,
        frame_out="ecliptic",
    )
    rho = coords.values[:, :3]

    rho1_hat = rho[0, :]
    rho2_hat = rho[1, :]
    rho3_hat = rho[2, :]
    q1 = coords_obs[0, :]
    q2 = coords_obs[1, :]
    q3 = coords_obs[2, :]
    q2_mag = np.linalg.norm(q2)

    # Make sure rhohats are unit vectors
    rho1_hat = rho1_hat / np.linalg.norm(rho1_hat)
    rho2_hat = rho2_hat / np.linalg.norm(rho2_hat)
    rho3_hat = rho3_hat / np.linalg.norm(rho3_hat)

    t1 = observation_times[0]
    t2 = observation_times[1]
    t3 = observation_times[2]
    t31 = t3 - t1
    t21 = t2 - t1
    t32 = t3 - t2

    A = _calcA(q1, q2, q3, rho1_hat, rho3_hat, t31, t32, t21)
    B = _calcB(q1, q3, rho1_hat, rho3_hat, t31, t32, t21, mu=mu)
    V = _calcV(rho1_hat, rho2_hat, rho3_hat)
    coseps2 = np.dot(q2, rho2_hat) / q2_mag
    C0 = V * t31 * q2_mag**4 / B
    h0 = -A / B

    if np.isnan(C0) or np.isnan(h0):
        return np.array([])

    # Find roots to eighth order polynomial
    all_roots = roots(
        [
            C0**2,
            0,
            -(q2_mag**2) * (h0**2 + 2 * C0 * h0 * coseps2 + C0**2),
            0,
            0,
            2 * q2_mag**5 * (h0 + C0 * coseps2),
            0,
            0,
            -(q2_mag**8),
        ]
    )

    # Keep only positive real roots (which should at most be 3)
    r2_mags = np.real(all_roots[np.isreal(all_roots) & (np.real(all_roots) >= 0)])
    num_solutions = len(r2_mags)
    if num_solutions == 0:
        return np.array([])

    orbits = []
    epochs = []
    for r2_mag in r2_mags:
        lambda1, lambda3 = _calcLambdas(r2_mag, t31, t32, t21, mu=mu)
        rho1, rho2, rho3 = _calcRhos(lambda1, lambda3, q1, q2, q3, rho1_hat, rho2_hat, rho3_hat, V)

        if np.dot(rho2, rho2_hat) < 0:
            continue

        # Test if we get the same rho2 as using equation 22 in Milani et al. 2008
        rho2_mag = (h0 - q2_mag**3 / r2_mag**3) * q2_mag / C0
        # np.testing.assert_almost_equal(np.dot(rho2_mag, rho2_hat), rho2)

        r1 = q1 + rho1
        r2 = q2 + rho2
        r3 = q3 + rho3

        if velocity_method == "gauss":
            v2 = calcGauss(r1, r2, r3, t1, t2, t3)
        elif velocity_method == "gibbs":
            v2 = calcGibbs(r1, r2, r3)
        elif velocity_method == "herrick+gibbs":
            v2 = calcHerrickGibbs(r1, r2, r3, t1, t2, t3)
        else:
            raise ValueError("velocity_method should be one of {'gauss', 'gibbs', 'herrick+gibbs'}")

        epoch = t2
        orbit = np.concatenate([r2, v2])

        if iterate is True:
            if iterator == "state transition":
                orbit = iterateStateTransition(
                    orbit,
                    t21,
                    t32,
                    q1,
                    q2,
                    q3,
                    rho1,
                    rho2,
                    rho3,
                    light_time=light_time,
                    mu=mu,
                    max_iter=max_iter,
                    tol=tol,
                )

        if light_time is True:
            rho2_mag = np.linalg.norm(orbit[:3] - q2)
            lt = rho2_mag / C
            epoch -= lt

        if np.linalg.norm(orbit[3:]) >= C:
            continue

        if (np.linalg.norm(orbit[:3]) > 300.0) or (np.linalg.norm(orbit[3:]) > 1.0):
            # Orbits that crash PYOORB:
            # 58366.84446725786 : 9.5544354809296721e+01  1.4093228616761269e+01 -6.6700146960148423e+00 -6.2618123281073522e+01 -9.4167879481188717e+00  4.4421501034359023e+0
            continue

        orbits.append(orbit)
        epochs.append(epoch)

    epochs = np.array(epochs)
    orbits = np.array(orbits)
    if len(orbits) > 0:
        epochs = epochs[~np.isnan(orbits).any(axis=1)]
        orbits = orbits[~np.isnan(orbits).any(axis=1)]

        return Orbits.from_kwargs(
            coordinates=CartesianCoordinates.from_kwargs(
                x=orbits[:, 0],
                y=orbits[:, 1],
                z=orbits[:, 2],
                vx=orbits[:, 3],
                vy=orbits[:, 4],
                vz=orbits[:, 5],
                time=Timestamp.from_mjd(epochs, scale="utc"),
                origin=Origin.from_kwargs(code=np.full(len(orbits), "SUN", dtype="object")),
                frame="ecliptic",
            )
        )

    else:
        return Orbits.empty()

1	import numpy as np	1✔
2	from adam_core.constants import Constants as c	1✔
3	from adam_core.coordinates import (	1✔
4	CartesianCoordinates,
5	Origin,
6	SphericalCoordinates,
7	transform_coordinates,
8	)
9	from adam_core.orbits import Orbits	1✔
10	from adam_core.time import Timestamp	1✔
11	from numba import jit	1✔
12	from numpy import roots	1✔
13
14	from .gibbs import calcGibbs	1✔
15	from .herrick_gibbs import calcHerrickGibbs	1✔
16	from .iterators import iterateStateTransition	1✔
17
18	__all__ = [	1✔
19	"_calcV",
20	"_calcA",
21	"_calcB",
22	"_calcLambdas",
23	"_calcRhos",
24	"approxLangrangeCoeffs",
25	"calcGauss",
26	"gaussIOD",
27	]
28
29	MU = c.MU	1✔
30	C = c.C	1✔
31
32
33	@jit("f8(f8[:], f8[:], f8[:])", nopython=True, cache=True)	1✔
34	def _calcV(rho1_hat, rho2_hat, rho3_hat):	1✔
35	# Vector triple product that gives the area of
36	# the "volume of the parallelepiped" or according to
37	# to Milani et al. 2008: 3x volume of the pyramid with vertices q, r1, r2, r3.
38	# Note that vector triple product rules apply here.
39	return np.dot(np.cross(rho1_hat, rho2_hat), rho3_hat)	×
40
41
42	@jit("f8(f8[:], f8[:], f8[:], f8[:], f8[:], f8, f8, f8)", nopython=True, cache=True)	1✔
43	def _calcA(q1, q2, q3, rho1_hat, rho3_hat, t31, t32, t21):	1✔
44	# Equation 21 from Milani et al. 2008
NEW 45	return np.linalg.norm(q2) ** 3 * np.dot(np.cross(rho1_hat, rho3_hat), (t32 * q1 - t31 * q2 + t21 * q3))	×
46
47
48	@jit("f8(f8[:], f8[:], f8[:], f8[:], f8, f8, f8, f8)", nopython=True, cache=True)	1✔
49	def _calcB(q1, q3, rho1_hat, rho3_hat, t31, t32, t21, mu=MU):	1✔
50	# Equation 19 from Milani et al. 2008
NEW 51	return mu / 6 * t32 * t21 * np.dot(np.cross(rho1_hat, rho3_hat), ((t31 + t32) * q1 + (t31 + t21) * q3))	×
52
53
54	@jit("UniTuple(f8, 2)(f8, f8, f8, f8, f8)", nopython=True, cache=True)	1✔
55	def _calcLambdas(r2_mag, t31, t32, t21, mu=MU):	1✔
56	# Equations 16 and 17 from Milani et al. 2008
57	lambda1 = t32 / t31 * (1 + mu / (6 * r2_mag*3) (t312 - t322))	×
58	lambda3 = t21 / t31 * (1 + mu / (6 * r2_mag*3) (t312 - t212))	×
59	return lambda1, lambda3	×
60
61
62	@jit(	1✔
63	"UniTuple(f8[:], 3)(f8, f8, f8[:], f8[:], f8[:], f8[:], f8[:], f8[:], f8)",
64	nopython=True,
65	cache=True,
66	)
67	def _calcRhos(lambda1, lambda3, q1, q2, q3, rho1_hat, rho2_hat, rho3_hat, V):	1✔
68	# This can be derived by taking a series of scalar products of the coplanarity condition equation
69	# with cross products of unit vectors in the direction of the observer, in particular, see Chapter 9 in
70	# Milani's book on the theory of orbit determination
71	numerator = -lambda1 * q1 + q2 - lambda3 * q3	×
72	rho1_mag = np.dot(numerator, np.cross(rho2_hat, rho3_hat)) / (lambda1 * V)	×
73	rho2_mag = np.dot(numerator, np.cross(rho1_hat, rho3_hat)) / V	×
74	rho3_mag = np.dot(numerator, np.cross(rho1_hat, rho2_hat)) / (lambda3 * V)	×
75	rho1 = rho1_mag * rho1_hat	×
76	rho2 = rho2_mag * rho2_hat	×
77	rho3 = rho3_mag * rho3_hat	×
78	return rho1, rho2, rho3	×
79
80
81	@jit("UniTuple(f8, 2)(f8, f8, f8)", nopython=True, cache=True)	1✔
82	def approxLangrangeCoeffs(r_mag, dt, mu=MU):	1✔
83	# Calculate the Lagrange coefficients (Gauss f and g series)
84	f = 1 - (1 / 2) * mu / r_mag*3 dt**2	×
85	g = dt - (1 / 6) * mu / r_mag*3 dt**2	×
86	return f, g	×
87
88
89	def calcGauss(r1, r2, r3, t1, t2, t3):	1✔
90	"""
91	Calculates the velocity vector at the location of the second position vector (r2) with Gauss's
92	original method.
93
94	.. math::
95	f_1 \approx 1 - \frac{1}{2}\frac{\mu}{r_2^3} (t_1 - t_2)^2
96
97	f_3 \approx 1 - \frac{1}{2}\frac{\mu}{r_2^3} (t_3 - t_2)^2
98
99	g_1 \approx (t_1 - t_2) - \frac{1}{6}\frac{\mu}{r_2^3} (t_1 - t_2)^2
100
101	g_3 \approx (t_3 - t_2) - \frac{1}{6}\frac{\mu}{r_2^3} (t_3 - t_2)^2
102
103	\vec{v}_2 = \frac{1}{f_1 g_3 - f_3 g_1} (-f_3 \vec{r}_1 + f_1 \vec{r}_3)
104
105
106	For more details on theory see Chapter 5 in Howard D. Curtis' "Orbital Mechanics for
107	Engineering Students".
108
109	Parameters
110	----------
111	r1 : `~numpy.ndarray` (3)
112	Heliocentric position vector at time 1 in cartesian coordinates in units
113	of AU.
114	r2 : `~numpy.ndarray` (3)
115	Heliocentric position vector at time 2 in cartesian coordinates in units
116	of AU.
117	r3 : `~numpy.ndarray` (3)
118	Heliocentric position vector at time 3 in cartesian coordinates in units
119	of AU.
120	t1 : float
121	Time at r1. Units of MJD or JD work or any decimal time format (one day is 1.00) as
122	long as all times are given in the same format.
123	t2 : float
124	Time at r2. Units of MJD or JD work or any decimal time format (one day is 1.00) as
125	long as all times are given in the same format.
126	t3 : float
127	Time at r3. Units of MJD or JD work or any decimal time format (one day is 1.00) as
128	long as all times are given in the same format.
129
130	Returns
131	-------
132	v2 : `~numpy.ndarray` (3)
133	Velocity of object at position r2 at time t2 in units of AU per day.
134	"""
135	t12 = t1 - t2	×
136	t32 = t3 - t2	×
137	r2_mag = np.linalg.norm(r2)	×
138	f1, g1 = approxLangrangeCoeffs(r2_mag, t12)	×
139	f3, g3 = approxLangrangeCoeffs(r2_mag, t32)	×
140	v2 = (1 / (f1 * g3 - f3 * g1)) * (-f3 * r1 + f1 * r3)	×
141	return v2	×
142
143
144	def gaussIOD(	1✔
145	coords,
146	observation_times,
147	coords_obs,
148	velocity_method="gibbs",
149	light_time=True,
150	iterate=True,
151	iterator="state transition",
152	mu=MU,
153	max_iter=10,
154	tol=1e-15,
155	):
156	"""
157	Compute up to three intial orbits using three observations in angular equatorial
158	coordinates.
159
160	Parameters
161	----------
162	coords : `~numpy.ndarray` (3, 2)
163	RA and Dec of three observations in units of degrees.
164	observation_times : `~numpy.ndarray` (3)
165	Times of the three observations in units of decimal days (MJD or JD for example).
166	coords_obs : `~numpy.ndarray` (3, 2)
167	Heliocentric position vector of the observer at times t in units of AU.
168	velocity_method : {'gauss', gibbs', 'herrick+gibbs'}, optional
169	Which method to use for calculating the velocity at the second observation.
170	[Default = 'gibbs']
171	light_time : bool, optional
172	Correct for light travel time.
173	[Default = True]
174	iterate : bool, optional
175	Iterate initial orbit using universal anomaly to better approximate the
176	Lagrange coefficients.
177	mu : float, optional
178	Gravitational parameter (GM) of the attracting body in units of
179	AU3 / d2.
180	max_iter : int, optional
181	Maximum number of iterations over which to converge to solution.
182	tol : float, optional
183	Numerical tolerance to which to compute chi using the Newtown-Raphson
184	method.
185
186	Returns
187	-------
188	epochs : `~numpy.ndarry` (<3)
189	Epochs in units of decimal days at which preliminary orbits are
190	defined.
191	orbits : `~numpy.ndarray` ((<3, 6) or (0))
192	Up to three preliminary orbits (as cartesian state vectors).
193	"""
194	coords = transform_coordinates(	1✔
195	SphericalCoordinates.from_kwargs(
196	lon=coords[:, 0],
197	lat=coords[:, 1],
198	origin=Origin.from_kwargs(code=np.full(len(coords), "Unknown", dtype="object")),
199	frame="equatorial",
200	).to_unit_sphere(),
201	representation_out=CartesianCoordinates,
202	frame_out="ecliptic",
203	)
204	rho = coords.values[:, :3]	1✔
205
206	rho1_hat = rho[0, :]	1✔
207	rho2_hat = rho[1, :]	1✔
208	rho3_hat = rho[2, :]	1✔
209	q1 = coords_obs[0, :]	1✔
210	q2 = coords_obs[1, :]	1✔
211	q3 = coords_obs[2, :]	1✔
212	q2_mag = np.linalg.norm(q2)	1✔
213
214	# Make sure rhohats are unit vectors
215	rho1_hat = rho1_hat / np.linalg.norm(rho1_hat)	1✔
216	rho2_hat = rho2_hat / np.linalg.norm(rho2_hat)	1✔
217	rho3_hat = rho3_hat / np.linalg.norm(rho3_hat)	1✔
218
219	t1 = observation_times[0]	1✔
220	t2 = observation_times[1]	1✔
221	t3 = observation_times[2]	1✔
222	t31 = t3 - t1	1✔
223	t21 = t2 - t1	1✔
224	t32 = t3 - t2	1✔
225
226	A = _calcA(q1, q2, q3, rho1_hat, rho3_hat, t31, t32, t21)	1✔
227	B = _calcB(q1, q3, rho1_hat, rho3_hat, t31, t32, t21, mu=mu)	1✔
228	V = _calcV(rho1_hat, rho2_hat, rho3_hat)	1✔
229	coseps2 = np.dot(q2, rho2_hat) / q2_mag	1✔
230	C0 = V * t31 * q2_mag**4 / B	1✔
231	h0 = -A / B	1✔
232
233	if np.isnan(C0) or np.isnan(h0):	1✔
234	return np.array([])	×
235
236	# Find roots to eighth order polynomial
237	all_roots = roots(	1✔
238	[
239	C0**2,
240	0,
241	-(q2_mag*2) (h0*2 + 2 C0 * h0 * coseps2 + C0**2),
242	0,
243	0,
244	2 * q2_mag*5 (h0 + C0 * coseps2),
245	0,
246	0,
247	-(q2_mag**8),
248	]
249	)
250
251	# Keep only positive real roots (which should at most be 3)
252	r2_mags = np.real(all_roots[np.isreal(all_roots) & (np.real(all_roots) >= 0)])	1✔
253	num_solutions = len(r2_mags)	1✔
254	if num_solutions == 0:	1✔
255	return np.array([])	×
256
257	orbits = []	1✔
258	epochs = []	1✔
259	for r2_mag in r2_mags:	1✔
260	lambda1, lambda3 = _calcLambdas(r2_mag, t31, t32, t21, mu=mu)	1✔
261	rho1, rho2, rho3 = _calcRhos(lambda1, lambda3, q1, q2, q3, rho1_hat, rho2_hat, rho3_hat, V)	1✔
262
263	if np.dot(rho2, rho2_hat) < 0:	1✔
264	continue	×
265
266	# Test if we get the same rho2 as using equation 22 in Milani et al. 2008
267	rho2_mag = (h0 - q2_mag3 / r2_mag3) * q2_mag / C0	1✔
268	# np.testing.assert_almost_equal(np.dot(rho2_mag, rho2_hat), rho2)
269
270	r1 = q1 + rho1	1✔
271	r2 = q2 + rho2	1✔
272	r3 = q3 + rho3	1✔
273
274	if velocity_method == "gauss":	1✔
275	v2 = calcGauss(r1, r2, r3, t1, t2, t3)	×
276	elif velocity_method == "gibbs":	1✔
277	v2 = calcGibbs(r1, r2, r3)	1✔
278	elif velocity_method == "herrick+gibbs":	×
279	v2 = calcHerrickGibbs(r1, r2, r3, t1, t2, t3)	×
280	else:
NEW 281	raise ValueError("velocity_method should be one of {'gauss', 'gibbs', 'herrick+gibbs'}")	×
282
283	epoch = t2	1✔
284	orbit = np.concatenate([r2, v2])	1✔
285
286	if iterate is True:	1✔
287	if iterator == "state transition":	×
288	orbit = iterateStateTransition(	×
289	orbit,
290	t21,
291	t32,
292	q1,
293	q2,
294	q3,
295	rho1,
296	rho2,
297	rho3,
298	light_time=light_time,
299	mu=mu,
300	max_iter=max_iter,
301	tol=tol,
302	)
303
304	if light_time is True:	1✔
305	rho2_mag = np.linalg.norm(orbit[:3] - q2)	1✔
306	lt = rho2_mag / C	1✔
307	epoch -= lt	1✔
308
309	if np.linalg.norm(orbit[3:]) >= C:	1✔
310	continue	×
311
312	if (np.linalg.norm(orbit[:3]) > 300.0) or (np.linalg.norm(orbit[3:]) > 1.0):	1✔
313	# Orbits that crash PYOORB:
314	# 58366.84446725786 : 9.5544354809296721e+01 1.4093228616761269e+01 -6.6700146960148423e+00 -6.2618123281073522e+01 -9.4167879481188717e+00 4.4421501034359023e+0
315	continue	×
316
317	orbits.append(orbit)	1✔
318	epochs.append(epoch)	1✔
319
320	epochs = np.array(epochs)	1✔
321	orbits = np.array(orbits)	1✔
322	if len(orbits) > 0:	1✔
323	epochs = epochs[~np.isnan(orbits).any(axis=1)]	1✔
324	orbits = orbits[~np.isnan(orbits).any(axis=1)]	1✔
325
326	return Orbits.from_kwargs(	1✔
327	coordinates=CartesianCoordinates.from_kwargs(
328	x=orbits[:, 0],
329	y=orbits[:, 1],
330	z=orbits[:, 2],
331	vx=orbits[:, 3],
332	vy=orbits[:, 4],
333	vz=orbits[:, 5],
334	time=Timestamp.from_mjd(epochs, scale="utc"),
335	origin=Origin.from_kwargs(code=np.full(len(orbits), "SUN", dtype="object")),
336	frame="ecliptic",
337	)
338	)
339
340	else:
341	return Orbits.empty()	×

moeyensj / thor / 11277630968

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous