Coverage for python / lsst / pipe / tasks / ssp / moid.py: 14%
152 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-17 09:20 +0000
1import numpy as np
2from astropy import units as u
3from astropy.constants import au as AU_CONST, GM_sun
4from collections import namedtuple
6# Convert AU to km from astropy
7AU_KM = AU_CONST.to_value(u.km)
10EarthElements = namedtuple(
11 "EarthElements",
12 ["a_AU", "e", "inc_deg", "Omega_deg", "omega_deg"],
13)
15# Named result type
16MOIDResult = namedtuple(
17 "MOIDResult",
18 [
19 "MOID_AU",
20 "DeltaV_kms",
21 "EclipticLongitude_deg",
22 "TrueAnomaly1_deg",
23 "TrueAnomaly2_deg",
24 ],
25)
28def earth_orbit_J2000():
29 """
30 Return canonical heliocentric Keplerian orbital elements for the Earth,
31 suitable for MOID computations, in the J2000 mean ecliptic and equinox
32 reference frame.
34 The elements describe a *fixed Keplerian ellipse* that approximates the
35 Earth's mean orbit at epoch J2000. This is the standard convention used
36 in MOID literature and software (e.g., Milani & Gronchi, NEODyS, JPL)
37 when computing the MOID between small bodies and the Earth.
39 Returns
40 -------
41 EarthElements
42 A namedtuple with fields:
44 - a_AU : float
45 Semi-major axis in astronomical units (AU).
46 Canonical value: 1.00000011 AU
48 - e : float
49 Eccentricity (dimensionless).
50 Canonical value: 0.01671022
52 - inc_deg : float
53 Inclination to the ecliptic in degrees.
54 Canonical value: 0.00005 deg
56 Note: by definition the ecliptic plane is the Earth's mean orbital
57 plane, so the true inclination is ~0. A very small non-zero value
58 is used in practice to avoid singularities in rotation matrices
59 and angle definitions.
61 - Omega_deg : float
62 Longitude of the ascending node in degrees, measured in the
63 ecliptic plane from the J2000 vernal equinox.
64 Canonical value: -11.26064 deg (equivalent to 348.73936 deg)
66 - omega_deg : float
67 Argument of perihelion in degrees, measured from the ascending
68 node along the Earth's orbital plane.
69 Canonical value: 102.94719 deg
71 Notes
72 -----
73 * These values are the widely used J2000 "mean elements" of the Earth’s
74 heliocentric orbit, compatible with VSOP87-style low-order models and
75 with the orbital elements quoted in many celestial mechanics references.
77 * They are intended for *geometric* computations such as MOID, where one
78 wants a fixed Keplerian ellipse representing the Earth's orbit, rather
79 than the true, time-varying perturbed Earth orbit.
81 * If you are using the MOIDSolver defined above, you can plug this
82 directly into the solver as the Earth orbit, e.g.:
84 earth = earth_orbit_J2000()
85 result = solver.compute(asteroid_elements, earth)
87 or reverse the order depending on which way you want to label the
88 "orbit 1" quantities.
90 * The epoch associated with these elements is J2000.0; MOID is a purely
91 geometric quantity, so there is no mean anomaly / phase dependence.
93 References
94 ----------
95 These values (or very close variants) are quoted in many sources, e.g.:
97 - Explanatory Supplement to the Astronomical Almanac
98 - VSOP87/IAU-style mean orbital element lists for planets
99 - Milani, A., & Gronchi, G. F., "Theory of Orbit Determination"
100 - Online resources summarizing J2000 planetary elements
101 """
102 return EarthElements(
103 a_AU=1.00000011,
104 e=0.01671022,
105 inc_deg=0.00005,
106 Omega_deg=-11.26064,
107 omega_deg=102.94719,
108 )
111class MOIDSolver:
112 """
113 MOID solver using an adaptive 2D grid in (f1, f2).
115 Units:
116 - Semi-major axes a: AU
117 - Distances / MOID: AU
118 - Velocities: km/s
119 - μ: astropy quantity (km^3/s^2) or float
120 - Public angles: degrees
121 - Internal angles: radians
122 """
124 def __init__(
125 self,
126 mu=GM_sun,
127 n_samples=128,
128 refine_factor=5.0,
129 tol_MOID_abs=1e-7 * u.AU,
130 tol_MOID_rel=1e-6,
131 max_refine=10,
132 ):
133 # Convert mu to float in km^3/s^2
134 self.mu = mu.to_value(u.km**3 / u.s**2) if hasattr(mu, "unit") else float(mu)
136 self.n_samples = int(n_samples)
137 self.refine_factor = float(refine_factor)
139 # Absolute tolerance expressed in AU
140 self.tol_MOID_abs = tol_MOID_abs.to_value(u.AU)
141 self.tol_MOID_rel = float(tol_MOID_rel)
143 self.max_refine = int(max_refine)
145 # ------------------------------------------------------------
146 # Geometry helpers
147 # ------------------------------------------------------------
149 @staticmethod
150 def _make_rotation_matrix(inc, Om, om):
151 """Rz(Omega) * Rx(inc) * Rz(omega)."""
152 cO = np.cos(Om)
153 sO = np.sin(Om)
154 ci = np.cos(inc)
155 si = np.sin(inc)
156 cw = np.cos(om)
157 sw = np.sin(om)
159 A00 = cw
160 A01 = -sw
161 A02 = 0.0
163 A10 = ci * sw
164 A11 = ci * cw
165 A12 = -si
167 A20 = si * sw
168 A21 = si * cw
169 A22 = ci
171 Q = np.empty((3, 3), float)
173 Q[0, 0] = cO * A00 - sO * A10
174 Q[0, 1] = cO * A01 - sO * A11
175 Q[0, 2] = cO * A02 - sO * A12
177 Q[1, 0] = sO * A00 + cO * A10
178 Q[1, 1] = sO * A01 + cO * A11
179 Q[1, 2] = sO * A02 + cO * A12
181 Q[2, 0] = A20
182 Q[2, 1] = A21
183 Q[2, 2] = A22
185 return Q
187 @classmethod
188 def _make_orbit_params(cls, a_AU, e, inc, Om, om):
189 """Return p in AU, rotation matrix Q."""
190 p = a_AU * (1 - e * e) # AU
191 Q = cls._make_rotation_matrix(inc, Om, om)
192 return p, e, Q
194 @staticmethod
195 def _orbit_positions(p_AU, e, Q, f):
196 """Vectorized positions in AU."""
197 f = np.asarray(f)
198 cf = np.cos(f)
199 sf = np.sin(f)
201 r = p_AU / (1 + e * cf) # AU
203 x_pf = r * cf
204 y_pf = r * sf
206 X = Q[0, 0] * x_pf + Q[0, 1] * y_pf
207 Y = Q[1, 0] * x_pf + Q[1, 1] * y_pf
208 Z = Q[2, 0] * x_pf + Q[2, 1] * y_pf
209 return np.column_stack((X, Y, Z)) # AU
211 @staticmethod
212 def _rv_from_params(p_AU, e, Q, f, mu_km3_s2):
213 """Return (r in AU, v in km/s)."""
214 cf = np.cos(f)
215 sf = np.sin(f)
217 # Position (AU)
218 r_AU = p_AU / (1 + e * cf)
219 x_pf = r_AU * cf
220 y_pf = r_AU * sf
222 X = Q[0, 0] * x_pf + Q[0, 1] * y_pf
223 Y = Q[1, 0] * x_pf + Q[1, 1] * y_pf
224 Z = Q[2, 0] * x_pf + Q[2, 1] * y_pf
225 r_vec_AU = np.array([X, Y, Z], float)
227 # Velocity (km/s): convert p (AU) to km
228 p_km = p_AU * AU_KM
229 fac = np.sqrt(mu_km3_s2 / p_km)
231 vx_pf = -fac * sf
232 vy_pf = fac * (e + cf)
234 VX = Q[0, 0] * vx_pf + Q[0, 1] * vy_pf
235 VY = Q[1, 0] * vx_pf + Q[1, 1] * vy_pf
236 VZ = Q[2, 0] * vx_pf + Q[2, 1] * vy_pf
237 v_vec_kms = np.array([VX, VY, VZ], float)
239 return r_vec_AU, v_vec_kms
241 # ------------------------------------------------------------
242 # Adaptive 2D grid MOID
243 # ------------------------------------------------------------
245 def _moid_grid_search(self, p1, e1, Q1, p2, e2, Q2):
246 """
247 Adaptive search in AU.
248 Returns: (MOID_AU, f1_best, f2_best)
249 """
250 n = self.n_samples
251 rf = self.refine_factor
252 tol_abs = self.tol_MOID_abs
253 tol_rel = self.tol_MOID_rel
254 max_ref = self.max_refine
256 two_pi = 2 * np.pi
258 center1 = np.pi
259 center2 = np.pi
260 width1 = two_pi
261 width2 = two_pi
263 best_d2 = np.inf
264 best_f1 = 0.0
265 best_f2 = 0.0
266 prev_moid = np.inf
268 for level in range(max_ref + 1):
269 f1_lo = center1 - width1 / 2
270 f1_hi = center1 + width1 / 2
271 f2_lo = center2 - width2 / 2
272 f2_hi = center2 + width2 / 2
274 f1_grid = (f1_lo + (np.arange(n) + 0.5) * (f1_hi - f1_lo) / n) % two_pi
275 f2_grid = (f2_lo + (np.arange(n) + 0.5) * (f2_hi - f2_lo) / n) % two_pi
277 r1 = self._orbit_positions(p1, e1, Q1, f1_grid)
278 r2 = self._orbit_positions(p2, e2, Q2, f2_grid)
280 diff = r1[:, None, :] - r2[None, :, :]
281 d2 = np.einsum("ijk,ijk->ij", diff, diff)
283 i0, j0 = np.unravel_index(np.argmin(d2), d2.shape)
284 d2_local = d2[i0, j0]
285 f1_local, f2_local = f1_grid[i0], f2_grid[j0]
287 if d2_local < best_d2:
288 best_d2 = d2_local
289 best_f1 = f1_local
290 best_f2 = f2_local
292 moid = float(np.sqrt(best_d2)) # AU
294 if level > 0:
295 delta = abs(moid - prev_moid)
296 scale = max(abs(moid), abs(prev_moid), 1.0)
297 if delta <= tol_abs + tol_rel * scale:
298 break
300 prev_moid = moid
301 center1 = f1_local
302 center2 = f2_local
303 width1 /= rf
304 width2 /= rf
306 return moid, best_f1, best_f2
308 # ------------------------------------------------------------
309 # Public API
310 # ------------------------------------------------------------
312 def compute(self, el1, el2):
313 """
314 Compute MOID and related quantities.
316 el = (a_AU, e, inc_deg, Ω_deg, ω_deg)
317 Returns MOIDResult namedtuple.
318 """
319 a1, e1, i1d, O1d, w1d = el1
320 a2, e2, i2d, O2d, w2d = el2
322 i1 = np.deg2rad(i1d)
323 O1 = np.deg2rad(O1d)
324 w1 = np.deg2rad(w1d)
326 i2 = np.deg2rad(i2d)
327 O2 = np.deg2rad(O2d)
328 w2 = np.deg2rad(w2d)
330 p1, e1, Q1 = self._make_orbit_params(a1, e1, i1, O1, w1)
331 p2, e2, Q2 = self._make_orbit_params(a2, e2, i2, O2, w2)
333 MOID_AU, f1_best, f2_best = self._moid_grid_search(p1, e1, Q1, p2, e2, Q2)
335 r1, v1 = self._rv_from_params(p1, e1, Q1, f1_best, self.mu)
336 r2, v2 = self._rv_from_params(p2, e2, Q2, f2_best, self.mu)
338 dv_kms = float(np.linalg.norm(v1 - v2))
340 lon = np.degrees(np.arctan2(r1[1], r1[0])) % 360
341 f1_deg = np.degrees(f1_best) % 360
342 f2_deg = np.degrees(f2_best) % 360
344 return MOIDResult(
345 MOID_AU=MOID_AU,
346 DeltaV_kms=dv_kms,
347 EclipticLongitude_deg=lon,
348 TrueAnomaly1_deg=f1_deg,
349 TrueAnomaly2_deg=f2_deg,
350 )
353# -------------------------------------------------------------------
354# Example
355# -------------------------------------------------------------------
356if __name__ == "__main__": 356 ↛ 357line 356 didn't jump to line 357 because the condition on line 356 was never true
357 solver = MOIDSolver(
358 mu=GM_sun,
359 n_samples=128,
360 refine_factor=5.0,
361 tol_MOID_abs=1e-8 * u.AU,
362 tol_MOID_rel=1e-8,
363 max_refine=12,
364 )
366 el1 = (1.0, 0.1, 5.0, 30.0, 45.0)
367 el2 = (1.5, 0.2, 15.0, 60.0, 10.0)
369 result = solver.compute(el1, el2)
370 print(result)