Coverage for python / lsst / pipe / tasks / ssp / moid.py: 14%
155 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-07 08:39 +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(epoch_mjd=51544.5):
29 """Return Earth’s heliocentric Keplerian elements at a given epoch.
31 Evaluates a linear secular model for Earth’s orbital elements
32 in the J2000 ecliptic and equinox frame. The model coefficients
33 are from JPL’s "Keplerian Elements for Approximate Positions of
34 the Major Planets" (Standish 1992, Table 1), valid for
35 3000 BC -- 3000 AD.
37 Parameters
38 ----------
39 epoch_mjd : float, optional
40 Epoch as Modified Julian Date. Default is 51544.5
41 (J2000.0 = 2000 Jan 1.5 TDB).
43 Returns
44 -------
45 EarthElements
46 A namedtuple with fields:
48 ``a_AU``
49 Semi-major axis in AU.
50 ``e``
51 Eccentricity.
52 ``inc_deg``
53 Inclination to the ecliptic in degrees. Set to a small
54 nonzero value (0.00005°) to avoid singularities in
55 rotation matrices; the true value is ~0 by definition
56 since the ecliptic *is* Earth’s mean orbital plane.
57 ``Omega_deg``
58 Longitude of the ascending node in degrees. Zero by
59 definition in the ecliptic frame.
60 ``omega_deg``
61 Argument of perihelion in degrees. Since Omega = 0 in
62 the ecliptic frame, this equals the longitude of
63 perihelion (ϖ).
65 Notes
66 -----
67 The linear model for each element is::
69 element(T) = element_0 + element_dot * T
71 where T is Julian centuries from J2000.0. The coefficients are:
73 ====== =============== =================
74 Param Value at J2000 Rate (per century)
75 ====== =============== =================
76 a 1.00000261 AU +0.00000562 AU
77 e 0.01671123 -0.00004392
78 ϖ 102.93768193° +0.32327364°
79 ====== =============== =================
81 Inclination and Omega are fixed at ~0 (ecliptic frame).
83 MOID is a purely geometric quantity (no mean anomaly / phase
84 dependence), but the *shape and orientation* of Earth’s orbit
85 do evolve slowly. Evaluating at the asteroid’s osculating
86 element epoch ensures the two orbits are compared at a
87 consistent time.
89 References
90 ----------
91 Standish, E.M. (1992). "Keplerian Elements for Approximate
92 Positions of the Major Planets." Solar System Dynamics Group,
93 JPL. https://ssd.jpl.nasa.gov/planets/approx_pos.html
94 """
95 T = (epoch_mjd - 51544.5) / 36525.0 # Julian centuries from J2000
96 return EarthElements(
97 a_AU=1.00000261 + 0.00000562 * T,
98 e=0.01671123 - 0.00004392 * T,
99 inc_deg=0.00005,
100 Omega_deg=0.0,
101 omega_deg=102.93768193 + 0.32327364 * T,
102 )
105def earth_orbit_J2000():
106 """Return Earth’s elements at J2000.0. Deprecated; use
107 `earth_orbit(epoch_mjd)` instead.
108 """
109 return earth_orbit(51544.5)
112class MOIDSolver:
113 """
114 MOID solver using an adaptive 2D grid in (f1, f2).
116 Units:
117 - Semi-major axes a: AU
118 - Distances / MOID: AU
119 - Velocities: km/s
120 - μ: astropy quantity (km^3/s^2) or float
121 - Public angles: degrees
122 - Internal angles: radians
123 """
125 def __init__(
126 self,
127 mu=GM_sun,
128 n_samples=128,
129 refine_factor=5.0,
130 tol_MOID_abs=1e-7 * u.AU,
131 tol_MOID_rel=1e-6,
132 max_refine=10,
133 ):
134 # Convert mu to float in km^3/s^2
135 self.mu = mu.to_value(u.km**3 / u.s**2) if hasattr(mu, "unit") else float(mu)
137 self.n_samples = int(n_samples)
138 self.refine_factor = float(refine_factor)
140 # Absolute tolerance expressed in AU
141 self.tol_MOID_abs = tol_MOID_abs.to_value(u.AU)
142 self.tol_MOID_rel = float(tol_MOID_rel)
144 self.max_refine = int(max_refine)
146 # ------------------------------------------------------------
147 # Geometry helpers
148 # ------------------------------------------------------------
150 @staticmethod
151 def _make_rotation_matrix(inc, Om, om):
152 """Rz(Omega) * Rx(inc) * Rz(omega)."""
153 cO = np.cos(Om)
154 sO = np.sin(Om)
155 ci = np.cos(inc)
156 si = np.sin(inc)
157 cw = np.cos(om)
158 sw = np.sin(om)
160 A00 = cw
161 A01 = -sw
162 A02 = 0.0
164 A10 = ci * sw
165 A11 = ci * cw
166 A12 = -si
168 A20 = si * sw
169 A21 = si * cw
170 A22 = ci
172 Q = np.empty((3, 3), float)
174 Q[0, 0] = cO * A00 - sO * A10
175 Q[0, 1] = cO * A01 - sO * A11
176 Q[0, 2] = cO * A02 - sO * A12
178 Q[1, 0] = sO * A00 + cO * A10
179 Q[1, 1] = sO * A01 + cO * A11
180 Q[1, 2] = sO * A02 + cO * A12
182 Q[2, 0] = A20
183 Q[2, 1] = A21
184 Q[2, 2] = A22
186 return Q
188 @classmethod
189 def _make_orbit_params(cls, a_AU, e, inc, Om, om):
190 """Return p in AU, rotation matrix Q."""
191 p = a_AU * (1 - e * e) # AU
192 Q = cls._make_rotation_matrix(inc, Om, om)
193 return p, e, Q
195 @staticmethod
196 def _orbit_positions(p_AU, e, Q, f):
197 """Vectorized positions in AU."""
198 f = np.asarray(f)
199 cf = np.cos(f)
200 sf = np.sin(f)
202 r = p_AU / (1 + e * cf) # AU
204 x_pf = r * cf
205 y_pf = r * sf
207 X = Q[0, 0] * x_pf + Q[0, 1] * y_pf
208 Y = Q[1, 0] * x_pf + Q[1, 1] * y_pf
209 Z = Q[2, 0] * x_pf + Q[2, 1] * y_pf
210 return np.column_stack((X, Y, Z)) # AU
212 @staticmethod
213 def _rv_from_params(p_AU, e, Q, f, mu_km3_s2):
214 """Return (r in AU, v in km/s)."""
215 cf = np.cos(f)
216 sf = np.sin(f)
218 # Position (AU)
219 r_AU = p_AU / (1 + e * cf)
220 x_pf = r_AU * cf
221 y_pf = r_AU * sf
223 X = Q[0, 0] * x_pf + Q[0, 1] * y_pf
224 Y = Q[1, 0] * x_pf + Q[1, 1] * y_pf
225 Z = Q[2, 0] * x_pf + Q[2, 1] * y_pf
226 r_vec_AU = np.array([X, Y, Z], float)
228 # Velocity (km/s): convert p (AU) to km
229 p_km = p_AU * AU_KM
230 fac = np.sqrt(mu_km3_s2 / p_km)
232 vx_pf = -fac * sf
233 vy_pf = fac * (e + cf)
235 VX = Q[0, 0] * vx_pf + Q[0, 1] * vy_pf
236 VY = Q[1, 0] * vx_pf + Q[1, 1] * vy_pf
237 VZ = Q[2, 0] * vx_pf + Q[2, 1] * vy_pf
238 v_vec_kms = np.array([VX, VY, VZ], float)
240 return r_vec_AU, v_vec_kms
242 # ------------------------------------------------------------
243 # Adaptive 2D grid MOID
244 # ------------------------------------------------------------
246 def _moid_grid_search(self, p1, e1, Q1, p2, e2, Q2):
247 """
248 Adaptive search in AU.
249 Returns: (MOID_AU, f1_best, f2_best)
250 """
251 n = self.n_samples
252 rf = self.refine_factor
253 tol_abs = self.tol_MOID_abs
254 tol_rel = self.tol_MOID_rel
255 max_ref = self.max_refine
257 two_pi = 2 * np.pi
259 center1 = np.pi
260 center2 = np.pi
261 width1 = two_pi
262 width2 = two_pi
264 best_d2 = np.inf
265 best_f1 = 0.0
266 best_f2 = 0.0
267 prev_moid = np.inf
269 for level in range(max_ref + 1):
270 f1_lo = center1 - width1 / 2
271 f1_hi = center1 + width1 / 2
272 f2_lo = center2 - width2 / 2
273 f2_hi = center2 + width2 / 2
275 f1_grid = (f1_lo + (np.arange(n) + 0.5) * (f1_hi - f1_lo) / n) % two_pi
276 f2_grid = (f2_lo + (np.arange(n) + 0.5) * (f2_hi - f2_lo) / n) % two_pi
278 r1 = self._orbit_positions(p1, e1, Q1, f1_grid)
279 r2 = self._orbit_positions(p2, e2, Q2, f2_grid)
281 diff = r1[:, None, :] - r2[None, :, :]
282 d2 = np.einsum("ijk,ijk->ij", diff, diff)
284 i0, j0 = np.unravel_index(np.argmin(d2), d2.shape)
285 d2_local = d2[i0, j0]
286 f1_local, f2_local = f1_grid[i0], f2_grid[j0]
288 if d2_local < best_d2:
289 best_d2 = d2_local
290 best_f1 = f1_local
291 best_f2 = f2_local
293 moid = float(np.sqrt(best_d2)) # AU
295 if level > 0:
296 delta = abs(moid - prev_moid)
297 scale = max(abs(moid), abs(prev_moid), 1.0)
298 if delta <= tol_abs + tol_rel * scale:
299 break
301 prev_moid = moid
302 center1 = f1_local
303 center2 = f2_local
304 width1 /= rf
305 width2 /= rf
307 return moid, best_f1, best_f2
309 # ------------------------------------------------------------
310 # Public API
311 # ------------------------------------------------------------
313 def compute(self, el1, el2):
314 """
315 Compute MOID and related quantities.
317 el = (a_AU, e, inc_deg, Ω_deg, ω_deg)
318 Returns MOIDResult namedtuple.
319 """
320 a1, e1, i1d, O1d, w1d = el1
321 a2, e2, i2d, O2d, w2d = el2
323 i1 = np.deg2rad(i1d)
324 O1 = np.deg2rad(O1d)
325 w1 = np.deg2rad(w1d)
327 i2 = np.deg2rad(i2d)
328 O2 = np.deg2rad(O2d)
329 w2 = np.deg2rad(w2d)
331 p1, e1, Q1 = self._make_orbit_params(a1, e1, i1, O1, w1)
332 p2, e2, Q2 = self._make_orbit_params(a2, e2, i2, O2, w2)
334 MOID_AU, f1_best, f2_best = self._moid_grid_search(p1, e1, Q1, p2, e2, Q2)
336 r1, v1 = self._rv_from_params(p1, e1, Q1, f1_best, self.mu)
337 r2, v2 = self._rv_from_params(p2, e2, Q2, f2_best, self.mu)
339 dv_kms = float(np.linalg.norm(v1 - v2))
341 lon = np.degrees(np.arctan2(r1[1], r1[0])) % 360
342 f1_deg = np.degrees(f1_best) % 360
343 f2_deg = np.degrees(f2_best) % 360
345 return MOIDResult(
346 MOID_AU=MOID_AU,
347 DeltaV_kms=dv_kms,
348 EclipticLongitude_deg=lon,
349 TrueAnomaly1_deg=f1_deg,
350 TrueAnomaly2_deg=f2_deg,
351 )
354# -------------------------------------------------------------------
355# Example
356# -------------------------------------------------------------------
357if __name__ == "__main__": 357 ↛ 358line 357 didn't jump to line 358 because the condition on line 357 was never true
358 solver = MOIDSolver(
359 mu=GM_sun,
360 n_samples=128,
361 refine_factor=5.0,
362 tol_MOID_abs=1e-8 * u.AU,
363 tol_MOID_rel=1e-8,
364 max_refine=12,
365 )
367 el1 = (1.0, 0.1, 5.0, 30.0, 45.0)
368 el2 = (1.5, 0.2, 15.0, 60.0, 10.0)
370 result = solver.compute(el1, el2)
371 print(result)