"""Linear covariance (Fisher) forecast of non-gravitational parameters from real astrometry epochs. Two dynamical models are provided. numeric_fisher : numerical integration of the equations of motion and of the variational equations (state-transition matrix + sensitivities to A1, A2, A3) in the field of the Sun, Earth and Moon (optionally all planets), with [revision] the Earth's J2 inside 0.002 au and the Sun's 1PN term. Used for every arc <= 180 d, every 1-10 m object and the validation set. kepler_fisher : heliocentric two-body motion; a radial 1/r^2 acceleration A1 is *exactly* a rescaling of the solar GM, so its partial is obtained from Kepler propagation with mu(1-beta). Used for multi-year arcs. Both return the normal matrix of [x0(6), A1, A2, A3] for optical (and radar) astrometry with the assumed data weights. No real orbit fit is performed: the trajectory comes from JPL's own solution; only the *information content* of the observation epochs/geometry is evaluated. """ import math import numpy as np import erfa from scipy.integrate import odeint from common import MU_SUN, C_AUD ARCSEC = np.pi / 648000.0 R_EARTH_AU = 6378.137 / 1.495978707e8 OBLIQ = 84381.448 * ARCSEC # IAU76 obliquity (JPL ecliptic J2000 elements) # Sun/body mass ratios (IAU / DE4xx values, rounded) J2_EARTH = 1.0826359e-3 C2_AUD = C_AUD * C_AUD MU_EARTH = MU_SUN / 332946.05 MU_MOON = MU_SUN / (332946.05 * 81.30057) PLANETS = {1: 6023600.0, 2: 408523.71, 4: 3098708.0, 5: 1047.3486, 6: 3497.898, 7: 22902.98, 8: 19412.24} _LEAP = np.array([41317, 41499, 41683, 42048, 42413, 42778, 43144, 43509, 43874, 44239, 44786, 45151, 45516, 46247, 47161, 47892, 48257, 48804, 49169, 49534, 50083, 50630, 51179, 53736, 54832, 56109, 57204, 57754], float) def utc_to_tt(mjd_utc): """MJD(UTC) -> MJD(TT); TDB-TT (<2 ms) is ignored.""" mjd_utc = np.asarray(mjd_utc, float) tai_utc = 9.0 + np.searchsorted(_LEAP, mjd_utc, side="right") return mjd_utc + (tai_utc + 32.184) / 86400.0 class _Tables: """Heliocentric positions/velocities of Earth, Moon and the planets from the ERFA analytic ephemerides (epv00, moon98, plan94), tabulated once on a 0.25-day grid and interpolated with cubic Hermite polynomials (Earth and Moon are reproduced to ~1e-10 au).""" T0, T1, H = 11000.0, 61700.0, 0.25 # MJD 1888 - 2027 P = V = None @classmethod def build(cls): if cls.P is not None: return tg = cls.T0 + cls.H * np.arange(int((cls.T1 - cls.T0) / cls.H) + 1) pvh, _ = erfa.epv00(2400000.5, tg) pm = erfa.moon98(2400000.5, tg) pp = [pvh["p"], pvh["p"] + pm["p"]]; vv = [pvh["v"], pvh["v"] + pm["v"]] for k in PLANETS: pv = erfa.plan94(2400000.5, tg, k) pp.append(pv["p"]); vv.append(pv["v"]) cls.P = np.ascontiguousarray(np.concatenate(pp, axis=1)) # (n, 27) cls.V = np.ascontiguousarray(np.concatenate(vv, axis=1)) * cls.H cls.NMAX = tg.size - 2 @classmethod def interp(cls, t, cols=slice(0, 3)): """Vectorised Hermite interpolation (positions) at times t.""" cls.build() x = (np.asarray(t, float) - cls.T0) / cls.H i = np.clip(x.astype(np.int64), 0, cls.NMAX) s = (x - i)[:, None] s2 = s * s; s3 = s2 * s P = cls.P[:, cols]; V = cls.V[:, cols] return ((2 * s3 - 3 * s2 + 1) * P[i] + (s3 - 2 * s2 + s) * V[i] + (-2 * s3 + 3 * s2) * P[i + 1] + (s3 - s2) * V[i + 1]) def earth_helio(mjd_tt): """Heliocentric Earth position and velocity [au, au/d], BCRS axes.""" t = np.atleast_1d(np.asarray(mjd_tt, float)) if t.size and t.min() > _Tables.T0 + 1 and t.max() < _Tables.T1 - 1: p = _Tables.interp(t) h = 0.01 v = (_Tables.interp(t + h) - _Tables.interp(t - h)) / (2 * h) return p, v pvh, _ = erfa.epv00(2400000.5, t) return pvh["p"], pvh["v"] def moon_geo(mjd_tt): pv = erfa.moon98(2400000.5, mjd_tt) return pv["p"], pv["v"] def observer_offset(mjd_utc, lon_deg, rho_cos, rho_sin): """Geocentric position of ground stations in the GCRS [au].""" mjd_utc = np.atleast_1d(np.asarray(mjd_utc, float)) tt = utc_to_tt(mjd_utc) gast = erfa.gst00b(2400000.5, mjd_utc) # UT1 ~ UTC th = gast + np.radians(lon_deg) tod = np.stack([rho_cos * np.cos(th), rho_cos * np.sin(th), rho_sin * np.ones_like(th)], axis=-1) bpn = erfa.pnm00b(2400000.5, tt) # GCRS -> true of date return np.einsum("nji,nj->ni", bpn, tod) * R_EARTH_AU def elements_to_state(a, e, inc, om, w, ma, mu=MU_SUN): """Heliocentric ecliptic-J2000 Keplerian elements (deg) -> ICRF state.""" inc, om, w, ma = np.radians([inc, om, w, ma]) if e < 1.0: E = ma for _ in range(60): dE = (E - e * np.sin(E) - ma) / (1.0 - e * np.cos(E)) E -= dE if abs(dE) < 1e-14: break xp = a * (np.cos(E) - e); yp = a * np.sqrt(1 - e * e) * np.sin(E) r = a * (1 - e * np.cos(E)); n = np.sqrt(mu / a**3) vxp = -a * n * np.sin(E) * a / r; vyp = a * n * np.sqrt(1 - e * e) * np.cos(E) * a / r else: raise ValueError("hyperbolic elements: supply a state vector instead") cO, sO, cw, sw, ci, si = np.cos(om), np.sin(om), np.cos(w), np.sin(w), np.cos(inc), np.sin(inc) P = np.array([cO * cw - sO * sw * ci, sO * cw + cO * sw * ci, sw * si]) Q = np.array([-cO * sw - sO * cw * ci, -sO * sw + cO * cw * ci, cw * si]) r_ecl = xp * P + yp * Q; v_ecl = vxp * P + vyp * Q ce, se = np.cos(OBLIQ), np.sin(OBLIQ) rot = np.array([[1, 0, 0], [0, ce, -se], [0, se, ce]]) return np.concatenate([rot @ r_ecl, rot @ v_ecl]) # --------------------------------------------------------------- numeric model class Ephem: """Perturber positions for the force model: Earth, Moon and (optionally) Mercury-Neptune, read from the global tables.""" def __init__(self, tmin, tmax, planets=True): _Tables.build() if tmin < _Tables.T0 + 1 or tmax > _Tables.T1 - 1: raise ValueError("epoch outside the tabulated ephemeris range") self.ncol = 27 if planets else 6 self.mu_list = [MU_EARTH, MU_MOON] + ([MU_SUN / r for r in PLANETS.values()] if planets else []) def pos_flat(self, t): """Flat python list [x0,y0,z0, x1,...] of the perturbers' positions.""" x = (t - _Tables.T0) / _Tables.H i = int(x) s = x - i s2 = s * s; s3 = s2 * s P = _Tables.P; V = _Tables.V; n = self.ncol return ((2 * s3 - 3 * s2 + 1) * P[i, :n] + (s3 - 2 * s2 + s) * V[i, :n] + (-2 * s3 + 3 * s2) * P[i + 1, :n] + (s3 - s2) * V[i + 1, :n]).tolist() class _Rhs: """Equations of motion + variational equations. y = [r, v, Z(6x9)] with Z = d(r,v)/d(r0, v0, A1, A2, A3). The force evaluation is written with python scalars (much faster than numpy for 3-vectors).""" def __init__(self, ng, eph): self.a1, self.a2, self.a3 = [float(x) for x in ng] self.eph = eph self.mu = eph.mu_list self.nb = len(self.mu) def __call__(self, t, y): x, yy, z, vx, vy, vz = y[0:6].tolist() pb = self.eph.pos_flat(t) r2 = x * x + yy * yy + z * z r1 = math.sqrt(r2) s3 = MU_SUN / (r2 * r1) ax = -s3 * x; ay = -s3 * yy; az = -s3 * z q = 3.0 * s3 / r2 gxx = -s3 + q * x * x; gyy = -s3 + q * yy * yy; gzz = -s3 + q * z * z gxy = q * x * yy; gxz = q * x * z; gyz = q * yy * z mu = self.mu for k in range(self.nb): bx = pb[3 * k]; by = pb[3 * k + 1]; bz = pb[3 * k + 2] dx = x - bx; dy_ = yy - by; dz = z - bz d2 = dx * dx + dy_ * dy_ + dz * dz w3 = mu[k] / (d2 * math.sqrt(d2)) b2 = bx * bx + by * by + bz * bz wb = mu[k] / (b2 * math.sqrt(b2)) ax -= w3 * dx + wb * bx; ay -= w3 * dy_ + wb * by; az -= w3 * dz + wb * bz if k == 0 and d2 < 4e-6: # [revision] Earth oblateness (J2) inside 0.002 au; the pole is taken along the ICRF # z-axis (precession since J2000 < 0.2 deg). Acceleration only: its contribution # to the variational equations is < 1e-3 of the point-mass gradient. cj = 1.5 * J2_EARTH * mu[0] * R_EARTH_AU * R_EARTH_AU / (d2 * d2 * math.sqrt(d2)) zz = 5.0 * dz * dz / d2 ax -= cj * dx * (1.0 - zz); ay -= cj * dy_ * (1.0 - zz); az -= cj * dz * (3.0 - zz) q = 3.0 * w3 / d2 gxx += -w3 + q * dx * dx; gyy += -w3 + q * dy_ * dy_; gzz += -w3 + q * dz * dz gxy += q * dx * dy_; gxz += q * dx * dz; gyz += q * dy_ * dz # [revision] Schwarzschild term of the Sun (first post-Newtonian order), acceleration only v2 = vx * vx + vy * vy + vz * vz rv = x * vx + yy * vy + z * vz fr_ = s3 / C2_AUD g1 = 4.0 * MU_SUN / r1 - v2 ax += fr_ * (g1 * x + 4.0 * rv * vx); ay += fr_ * (g1 * yy + 4.0 * rv * vy); az += fr_ * (g1 * z + 4.0 * rv * vz) # radial / transverse / normal unit vectors, scaled by (1 au / r)^2 hx = yy * vz - z * vy; hy = z * vx - x * vz; hz = x * vy - yy * vx hn = math.sqrt(hx * hx + hy * hy + hz * hz) rx = x / r1; ry = yy / r1; rz = z / r1 nx = hx / hn; ny = hy / hn; nz = hz / hn tx = ny * rz - nz * ry; ty = nz * rx - nx * rz; tz = nx * ry - ny * rx i2 = 1.0 / r2 a1, a2, a3 = self.a1, self.a2, self.a3 ax += (a1 * rx + a2 * tx + a3 * nx) * i2 ay += (a1 * ry + a2 * ty + a3 * ny) * i2 az += (a1 * rz + a2 * tz + a3 * nz) * i2 dy = np.empty(60) dy[0] = vx; dy[1] = vy; dy[2] = vz; dy[3] = ax; dy[4] = ay; dy[5] = az Z = y[6:60].reshape(6, 9) dZ = np.empty((6, 9)) dZ[0:3] = Z[3:6] dZ[3:6] = np.array(((gxx, gxy, gxz), (gxy, gyy, gyz), (gxz, gyz, gzz))) @ Z[0:3] dZ[3:6, 6:9] += np.array(((rx * i2, tx * i2, nx * i2), (ry * i2, ty * i2, ny * i2), (rz * i2, tz * i2, nz * i2))) dy[6:60] = dZ.ravel() return dy def propagate(state0, t0, times, ng=(0, 0, 0), planets=True, rtol=1e-10): """Integrate trajectory + variational equations from t0 to every epoch in `times` (any order). Returns Y (60, n) and the number of force evaluations.""" times = np.asarray(times, float) eph = Ephem(min(times.min(), t0) - 1.0, max(times.max(), t0) + 1.0, planets) rhs = _Rhs(ng, eph) Z0 = np.zeros((6, 9)); Z0[:, 0:6] = np.eye(6) y0 = np.concatenate([state0, Z0.ravel()]) Y = np.empty((60, times.size)) nfev = 0 ok = True for sign in (-1, 1): m = (times < t0) if sign < 0 else (times >= t0) if not m.any(): continue tu, inv = np.unique(times[m], return_inverse=True) if sign < 0: tu = tu[::-1]; inv = tu.size - 1 - inv tt = np.concatenate([[t0], tu]) if tu[0] != t0 else tu sol, info = odeint(rhs, y0, tt, tfirst=True, rtol=rtol, atol=1e-14, mxstep=200000, full_output=True, hmax=20.0) if tu[0] != t0: sol = sol[1:] Y[:, np.flatnonzero(m)] = sol[inv].T nfev += int(info["nfe"][-1]) ok = ok and (info["message"] == "Integration successful.") return Y, nfev, ok def _design_rows(dvec, dist, dr_dtheta): """dvec (N,3) observer->object, dist (N,), dr_dtheta (N,3,P) -> two (N,P) blocks.""" u = dvec / dist[:, None] ra = np.arctan2(u[:, 1], u[:, 0]); dec = np.arcsin(np.clip(u[:, 2], -1, 1)) ea = np.stack([-np.sin(ra), np.cos(ra), np.zeros_like(ra)], axis=1) ed = np.stack([-np.sin(dec) * np.cos(ra), -np.sin(dec) * np.sin(ra), np.cos(dec)], axis=1) Ha = np.einsum("ni,nip->np", ea, dr_dtheta) / dist[:, None] Hd = np.einsum("ni,nip->np", ed, dr_dtheta) / dist[:, None] return ra, dec, Ha, Hd def numeric_fisher(state0, mjd_ref_tt, obs, radar=None, ng=(0, 0, 0), planets=True, rtol=1e-10): """obs: dict with mjd_utc, ra, dec (deg), sig (arcsec), off (N,3) topocentric offset [au]. radar: optional dict with mjd_utc, kind ('d'/'D'), sig, freq. Returns dict with normal matrix N (9x9: x0, A1, A2, A3), residuals etc. Light time is applied to first order (r(t-lt) = r(t) - v*lt).""" tt = utc_to_tt(obs["mjd_utc"]) nopt = tt.size tr = utc_to_tt(radar["mjd_utc"]) if radar is not None and len(radar["mjd_utc"]) else np.array([]) Yall, nfev, ok = propagate(np.asarray(state0, float), mjd_ref_tt, np.concatenate([tt, tr]), ng, planets, rtol) y = Yall[:, :nopt] ep, ev = earth_helio(tt) obspos = ep + obs["off"] r = y[0:3].T; v = y[3:6].T dist = np.linalg.norm(r - obspos, axis=1) dvec = r - v * (dist / C_AUD)[:, None] - obspos dist = np.linalg.norm(dvec, axis=1) dr = y[6:60].T.reshape(-1, 6, 9)[:, 0:3, :] # (N,3,9) ra, dec, Ha, Hd = _design_rows(dvec, dist, dr) w = 1.0 / (obs["sig"] * ARCSEC) ** 2 N = (Ha * w[:, None]).T @ Ha + (Hd * w[:, None]).T @ Hd dra = (np.radians(obs["ra"]) - ra + np.pi) % (2 * np.pi) - np.pi res_a = dra * np.cos(dec) / ARCSEC res_d = (np.radians(obs["dec"]) - dec) / ARCSEC out = dict(N_opt=N.copy(), res_ra=res_a, res_dec=res_d, dist=dist, nfev=nfev, ok=ok, Ha=Ha, Hd=Hd) if tr.size: y = Yall[:, nopt:] ep, ev = earth_helio(tr) r = y[0:3].T; v = y[3:6].T d0 = np.linalg.norm(r - ep, axis=1) dvec = r - v * (d0 / C_AUD)[:, None] - ep d0 = np.linalg.norm(dvec, axis=1); u = dvec / d0[:, None] Zr = y[6:60].T.reshape(-1, 6, 9) dr = Zr[:, 0:3, :]; dv = Zr[:, 3:6, :] vrel = v - ev vperp = vrel - np.sum(vrel * u, axis=1)[:, None] * u d_range = np.einsum("ni,nip->np", u, dr) # au per unit theta d_rrate = np.einsum("ni,nip->np", u, dv) + np.einsum("ni,nip->np", vperp, dr) / d0[:, None] Hr = np.where((radar["kind"] == "d")[:, None], 2.0 * d_range / C_AUD * 86400e6, # delay, microseconds -2.0 * radar["freq"][:, None] * 1e6 * d_rrate / C_AUD) # Doppler, Hz wr = 1.0 / radar["sig"] ** 2 N = N + (Hr * wr[:, None]).T @ Hr out["N"] = N return out # ---------------------------------------------------------------- Kepler model def _kepler_rv(state, mu, dt): """Propagate elliptical two-body states. state (V,6), mu (V,), dt (T,) -> r (V,T,3), v (V,T,3).""" r0 = state[:, 0:3]; v0 = state[:, 3:6] r0n = np.linalg.norm(r0, axis=1); v2 = np.sum(v0 * v0, axis=1) a = 1.0 / (2.0 / r0n - v2 / mu) n = np.sqrt(mu / a**3) ecosE0 = 1.0 - r0n / a esinE0 = np.sum(r0 * v0, axis=1) / np.sqrt(mu * a) e = np.hypot(ecosE0, esinE0) E0 = np.arctan2(esinE0, ecosE0) M = (E0 - esinE0)[:, None] + n[:, None] * dt[None, :] k = np.round(M / (2 * np.pi)); Mr = M - 2 * np.pi * k ee = e[:, None] E = Mr + ee * np.sin(Mr) / (1.0 - np.sin(Mr + ee) + np.sin(Mr)) for _ in range(12): f = E - ee * np.sin(E) - Mr fp = 1.0 - ee * np.cos(E) fpp = ee * np.sin(E) E = E - f / (fp - 0.5 * f * fpp / fp) dE = (E - E0[:, None]) + 2 * np.pi * k cdE = np.cos(dE); sdE = np.sin(dE) aa = a[:, None]; nn = n[:, None]; rr0 = r0n[:, None] fF = 1.0 - aa / rr0 * (1.0 - cdE) gG = dt[None, :] + (sdE - dE) / nn r = fF[:, :, None] * r0[:, None, :] + gG[:, :, None] * v0[:, None, :] rn = np.linalg.norm(r, axis=2) fd = -aa * aa * nn / (rn * rr0) * sdE gd = 1.0 - aa / rn * (1.0 - cdE) v = fd[:, :, None] * r0[:, None, :] + gd[:, :, None] * v0[:, None, :] return r, v def _apparition_shifts(state0, mjd_ref_tt, tt, uobs, obspos, gap=60.0, nsel=5): """Two-body propagation of one osculating orbit over decades drifts along the orbit (planetary perturbations change the mean motion). For each apparition (observations separated by < `gap` days) find the time shift that makes the two-body position match the *observed* directions; this restores the correct observer-object geometry for the partials.""" st = state0[None, :]; mu = np.array([MU_SUN]) r0n = np.linalg.norm(state0[:3]); v2 = state0[3:] @ state0[3:] a = 1.0 / (2.0 / r0n - v2 / MU_SUN) P = 2 * np.pi * np.sqrt(a**3 / MU_SUN) order = np.argsort(tt) ts = tt[order] breaks = np.flatnonzero(np.diff(ts) > gap) + 1 shift = np.zeros(tt.size) sep_out = np.zeros(tt.size) for seg in np.split(order, breaks): k = seg[np.unique(np.linspace(0, seg.size - 1, min(nsel, seg.size)).astype(int))] tk = tt[k]; uk = uobs[k]; ok = obspos[k] def cost(deltas): dt = (tk[None, :] + deltas[:, None] - mjd_ref_tt).ravel() r, _ = _kepler_rv(st, mu, dt) d = r[0].reshape(deltas.size, tk.size, 3) - ok[None, :, :] d /= np.linalg.norm(d, axis=2)[:, :, None] return np.sum(1.0 - np.sum(d * uk[None, :, :], axis=2), axis=1) # ~ sum sep^2 / 2 best = 0.0 c0 = cost(np.array([0.0]))[0] if c0 > 0.5 * tk.size * (0.3 * np.pi / 180.0) ** 2: # only when off by > 0.3 deg width = 0.5 * P; n = 361 for _ in range(4): grid = best + np.linspace(-width, width, n) c = cost(grid) j = int(np.argmin(c)) best = grid[j] width = 2.0 * width / (n - 1); n = 41 c0 = c[j] shift[seg] = best sep_out[seg] = np.sqrt(2.0 * c0 / tk.size) return shift, sep_out def kepler_fisher(state0, mjd_ref_tt, obs, radar=None, phase_fix=True): """Two-body version. Parameters: x0(6), A1 (exact: mu -> mu(1-A1/mu)), A2 (secular along-track approximation), A3 not modelled (column of zeros is dropped by the caller).""" state0 = np.asarray(state0, float) dx = np.array([1e-7, 1e-7, 1e-7, 2e-9, 2e-9, 2e-9]) * 0.5 dA1 = 1e-10 V = np.tile(state0, (15, 1)); mu = np.full(15, MU_SUN) for j in range(6): V[1 + 2 * j, j] += dx[j]; V[2 + 2 * j, j] -= dx[j] mu[13] -= dA1; mu[14] += dA1 # +A1 lowers the effective GM (r0 = 1 au) r0n = np.linalg.norm(state0[:3]); v2 = state0[3:] @ state0[3:] a = 1.0 / (2.0 / r0n - v2 / MU_SUN); n = np.sqrt(MU_SUN / a**3) h = np.cross(state0[:3], state0[3:]); ome2 = (h @ h) / MU_SUN / a def partials(tt, shift): dt = tt + shift - mjd_ref_tt r, v = _kepler_rv(V, mu, dt) dr = np.empty((tt.size, 3, 8)); dv = np.empty((tt.size, 3, 8)) for j in range(6): dr[:, :, j] = (r[1 + 2 * j] - r[2 + 2 * j]) / (2 * dx[j]) dv[:, :, j] = (v[1 + 2 * j] - v[2 + 2 * j]) / (2 * dx[j]) dr[:, :, 6] = (r[13] - r[14]) / (2 * dA1) dv[:, :, 6] = (v[13] - v[14]) / (2 * dA1) # transverse A2: da/dt = 2 A2 /(n a^2 (1-e^2)) -> along-track time shift dtau = -1.5 * dt**2 / (n * a**3 * ome2) dr[:, :, 7] = v[0] * dtau[:, None] acc = -MU_SUN * r[0] / np.linalg.norm(r[0], axis=1)[:, None] ** 3 dv[:, :, 7] = acc * dtau[:, None] return r[0], v[0], dr, dv tt = utc_to_tt(obs["mjd_utc"]) ep, ev = earth_helio(tt) obspos = ep + obs["off"] ra_o = np.radians(obs["ra"]); de_o = np.radians(obs["dec"]) uobs = np.stack([np.cos(de_o) * np.cos(ra_o), np.cos(de_o) * np.sin(ra_o), np.sin(de_o)], axis=1) shift = np.zeros(tt.size) if phase_fix: r0_, _ = _kepler_rv(V[:1], mu[:1], tt - mjd_ref_tt) d0_ = r0_[0] - obspos d0_ /= np.linalg.norm(d0_, axis=1)[:, None] sep = np.degrees(np.arccos(np.clip(np.sum(d0_ * uobs, axis=1), -1, 1))) out_sep0 = float(np.median(sep)) if out_sep0 > 1.0: # two-body geometry off by > 1 deg: re-phase shift, _ = _apparition_shifts(state0, mjd_ref_tt, tt, uobs, obspos) r, v, dr, dv = partials(tt, shift) dist = np.linalg.norm(r - obspos, axis=1) dvec = r - v * (dist / C_AUD)[:, None] - obspos dist = np.linalg.norm(dvec, axis=1) ra, dec, Ha, Hd = _design_rows(dvec, dist, dr) w = 1.0 / (obs["sig"] * ARCSEC) ** 2 N = (Ha * w[:, None]).T @ Ha + (Hd * w[:, None]).T @ Hd dra = (ra_o - ra + np.pi) % (2 * np.pi) - np.pi out = dict(N_opt=N.copy(), res_ra=dra * np.cos(dec) / ARCSEC, res_dec=(de_o - dec) / ARCSEC, dist=dist, nfev=0, ok=True, max_shift=float(np.abs(shift).max()), Ha=Ha, Hd=Hd) if radar is not None and len(radar["mjd_utc"]): tr = utc_to_tt(radar["mjd_utc"]) ep, ev = earth_helio(tr) # radar epochs inherit the shift of the nearest optical observation sh = shift[np.argmin(np.abs(tt[None, :] - tr[:, None]), axis=1)] r, v, dr, dv = partials(tr, sh) d0 = np.linalg.norm(r - ep, axis=1) dvec = r - v * (d0 / C_AUD)[:, None] - ep; d0 = np.linalg.norm(dvec, axis=1); u = dvec / d0[:, None] vrel = v - ev vperp = vrel - np.sum(vrel * u, axis=1)[:, None] * u d_range = np.einsum("ni,nip->np", u, dr) d_rrate = np.einsum("ni,nip->np", u, dv) + np.einsum("ni,nip->np", vperp, dr) / d0[:, None] Hr = np.where((radar["kind"] == "d")[:, None], 2.0 * d_range / C_AUD * 86400e6, -2.0 * radar["freq"][:, None] * 1e6 * d_rrate / C_AUD) N = N + (Hr * (1.0 / radar["sig"] ** 2)[:, None]).T @ Hr out["N"] = N return out # ------------------------------------------------- measurement of A1 (revision) def _solve_normal(N, y): """Solve N x = y for a symmetric positive (semi-)definite N with column scaling. Returns (x, covariance diagonal) or (None, None) when N is numerically singular.""" d = np.sqrt(np.clip(np.diag(N), 1e-300, None)) Ns = N / np.outer(d, d) w, V = np.linalg.eigh(Ns) if w.min() <= 1e-13 * w.max(): return None, None Ninv = (V / w) @ V.T return (Ninv @ (y / d)) / d, np.diag(Ninv) / d**2 LIN_MARGIN = 0.75 # the linearised estimate of a LARGE signal can be low by up to ~23% (noisy # injection test, s05_injection.main_measurement: recovered/injected >= 0.77 # in 160 trials); upper bounds therefore use max(A1, 0) / LIN_MARGIN + 3 sigma RESOLVE_Z = 4.0 # a linearised estimate may clear a flagged object only where it is itself not # significant (A1/sigma < RESOLVE_Z): for significant values the linear step is biased def upper_bound(a1, sig, lin=LIN_MARGIN): """3-sigma upper bound on an outward A1 from a linearised estimate (sunward values count as 0).""" return np.maximum(a1, 0.0) / lin + 3.0 * sig def fit_a1(Ha, Hd, res_a, res_d, sig, groups=None, clip=3.0, max_iter=6, npar=7, max_groups=400, sig_d=None, cols=None, a1_ref=0.0): """Linearised estimate of the radial parameter A1 from *real* astrometric residuals. Ha, Hd : (n, >=7) design rows [rad per unit parameter] of the 6 state components and A1 along the reference trajectory; res_a, res_d: observed-minus-computed [arcsec]; sig: assumed 1-sigma [arcsec]. The residuals are fitted by weighted least squares with (a) 6 state corrections and (b) 6 state corrections + A1. Outliers are rejected iteratively on the residuals of fit (b), never of fit (a), so that a real radial signal is not clipped away. groups : optional integer label per observation (station-night). If given, the fit is repeated leaving each group out in turn; the smallest resulting 3-sigma upper bound (upper_bound(), on the total A1 = a1_ref + correction) is returned together with its estimate and sigma (used only to vet candidates, see s04_limits). a1_ref : A1 of the reference trajectory (only used for the leave-one-out bound). sig_d : optional separate sigma for declination (default: sig for both coordinates). cols : optional list of design columns to solve for; the first six must be the state and the LAST one is the reported parameter. cols=[0,1,2,3,4,5,7,6] solves the transverse parameter A2 as a nuisance parameter together with A1 (needed on multi-year arcs, where a Yarkovsky drift otherwise leaks into A1). Returned A1 values are *corrections* to the A1 of the reference trajectory [au/d^2].""" n = len(sig) if sig_d is None: sig_d = sig if cols is None: cols = list(range(npar)) npar = len(cols) A = np.vstack([Ha[:, cols] / (sig * ARCSEC)[:, None], Hd[:, cols] / (sig_d * ARCSEC)[:, None]]) b = np.concatenate([res_a / sig, res_d / sig_d]) out = dict(m_n_all=n) def solve(keep): k2 = np.concatenate([keep, keep]) Ak, bk = A[k2], b[k2] N = Ak.T @ Ak; y = Ak.T @ bk x7, c7 = _solve_normal(N, y) if x7 is None: return None x6, _ = _solve_normal(N[:6, :6], y[:6]) r7 = b - A @ x7 # all observations, also rejected ones chi7 = float(np.sum(r7[k2] ** 2)) chi6 = float(np.sum((bk - Ak[:, :6] @ x6) ** 2)) if x6 is not None else np.nan return dict(x=x7, sig=float(np.sqrt(c7[npar - 1])), r7=r7, chi6=chi6, chi7=chi7, N=N, y=y, bb=float(bk @ bk), nk=int(keep.sum())) keep = np.ones(n, bool) f = solve(keep) if f is None: return out out.update(m_A1_all=float(f["x"][npar - 1]), m_sig_all=f["sig"], m_rms6_all=float(np.sqrt(f["chi6"] / max(2 * n - 6, 1))), m_rms7_all=float(np.sqrt(f["chi7"] / max(2 * n - npar, 1)))) for _ in range(max_iter): bad = (np.abs(f["r7"][:n]) > clip) | (np.abs(f["r7"][n:]) > clip) new = ~bad if new.sum() < 5 or np.array_equal(new, keep): break g = solve(new) if g is None: break keep, f = new, g nk = f["nk"] rms7 = float(np.sqrt(f["chi7"] / max(2 * nk - npar, 1))) out.update(m_n=nk, m_A1=float(f["x"][npar - 1]), m_sig=f["sig"], m_rms6=float(np.sqrt(f["chi6"] / max(2 * nk - 6, 1))), m_rms7=rms7, m_dchi2=f["chi6"] - f["chi7"]) if groups is not None: lab = np.asarray(groups)[keep] ug, inv = np.unique(lab, return_inverse=True) out["m_ngrp"] = int(ug.size) if 2 <= ug.size <= max_groups: k2 = np.concatenate([keep, keep]) Ak, bk = A[k2], b[k2] inv2 = np.concatenate([inv, inv]) best_ul, best_z, best_g, best_n, best_x, best_s = np.inf, np.inf, -1, 0, np.nan, np.nan for j in range(ug.size): m = inv2 == j nj = nk - int(m.sum() // 2) if nj < 5: continue Aj, bj = Ak[m], bk[m] N = f["N"] - Aj.T @ Aj; y = f["y"] - Aj.T @ bj x7, c7 = _solve_normal(N, y) if x7 is None: continue chi = f["bb"] - float(bj @ bj) - float(x7 @ y) sj = float(np.sqrt(c7[npar - 1])) * max(1.0, np.sqrt(max(chi, 0.0) / max(2 * nj - npar, 1))) xt = float(x7[npar - 1]) + a1_ref ul = float(upper_bound(xt, sj)) if ul < best_ul: best_ul, best_g, best_n, best_x, best_s = ul, j, int(m.sum() // 2), xt, sj best_z = min(best_z, xt / sj) if np.isfinite(best_ul): out.update(m_lono_ul=best_ul, m_lono_A1=best_x, m_lono_sig=best_s, m_lono_zmin=best_z, m_lono_n=best_n) return out def marginal_sigma(N, idx, target): """1-sigma of parameter `target` when the parameters `idx` are solved for.""" idx = list(idx) M = N[np.ix_(idx, idx)] d = np.sqrt(np.clip(np.diag(M), 1e-300, None)) Ms = M / np.outer(d, d) w, Vv = np.linalg.eigh(Ms) if w.min() <= 1e-13 * w.max(): return np.inf j = idx.index(target) var = np.sum(Vv[j, :] ** 2 / w) return float(np.sqrt(var) / d[j])