""" orbit_determination.py ======================= Détermination d'orbite préliminaire d'un astéroïde à partir de trois observations angulaires (RA, Dec, temps) : méthode classique de Gauss. Référence de l'algorithme : H.D. Curtis, "Orbital Mechanics for Engineering Students", chap. 5 (méthode de Gauss + Gibbs / Herrick-Gibbs), et Vallado, "Fundamentals of Astrodynamics and Applications". Limites importantes (à lire avant utilisation) ----------------------------------------------- - La méthode de Gauss ne fonctionne bien QUE pour un arc d'observation assez court (quelques heures à quelques jours) : c'est une méthode de "première orbite", pas un ajustement précis. - Elle est très sensible aux erreurs d'astrométrie ; avec seulement 3 points, il n'y a pas de moyen de contrôler la qualité intrinsèquement. Pour un résultat fiable, on utilisera plutôt : * plusieurs triplets d'observations (parmi toutes les images disponibles) et on comparera les éléments obtenus, * puis idéalement un ajustement aux moindres carrés (differential correction) ou un logiciel dédié (Find_Orb, OrbFit) pour affiner. - On suppose ici des observations optiques géocentriques/topocentriques et on calcule une orbite héliocentrique. La position de l'observateur est calculée via les éphémérides intégrées à astropy (JPL via get_body_barycentric), en négligeant si besoin l'offset topocentrique (faible pour un astéroïde lointain mais non négligeable pour un objet proche de la Terre). Prérequis --------- pip install astropy numpy --break-system-packages Utilisation ----------- from orbit_determination import Observation, gauss_preliminary_orbit obs = [ Observation(time=t1, ra_deg=.., dec_deg=..), Observation(time=t2, ra_deg=.., dec_deg=..), Observation(time=t3, ra_deg=.., dec_deg=..), ] elements = gauss_preliminary_orbit(obs) print(elements) """ from __future__ import annotations from dataclasses import dataclass import numpy as np from astropy.time import Time from astropy.coordinates import get_body_barycentric, EarthLocation, ITRS import astropy.units as u # Constante gravitationnelle du Soleil (km^3/s^2) MU_SUN = 1.32712440018e11 AU_KM = 149597870.7 # Obliquité moyenne de l'écliptique à J2000 (deg) - IAU 2006 OBLIQUITY_J2000_DEG = 23.43929111 def equatorial_to_ecliptic(vec: np.ndarray) -> np.ndarray: """Fait pivoter un vecteur (position ou vitesse) du repère équatorial ICRS/J2000 vers le repère écliptique J2000, par rotation d'angle -obliquité autour de l'axe X (qui est commun aux deux repères, dirigé vers le point vernal). C'EST UNE ÉTAPE INDISPENSABLE avant de calculer des éléments orbitaux "écliptiques" classiques : sans cette rotation, l'inclinaison calculée est mesurée par rapport à l'équateur céleste et non à l'écliptique, ce qui la fausse d'un angle proche de l'obliquité (23.44°) pour les objets proches du plan de l'écliptique. """ eps = np.deg2rad(OBLIQUITY_J2000_DEG) x, y, z = vec return np.array([ x, y * np.cos(eps) + z * np.sin(eps), -y * np.sin(eps) + z * np.cos(eps), ]) @dataclass class Observation: time: Time ra_deg: float dec_deg: float @dataclass class OrbitalElements: a_AU: float # demi-grand axe (UA) e: float # excentricité i_deg: float # inclinaison raan_deg: float # longitude du nœud ascendant (Omega) argp_deg: float # argument du périhélie (omega) nu_deg: float # anomalie vraie à l'époque M_deg: float # anomalie moyenne à l'époque epoch: Time # instant de référence (temps de l'observation centrale) period_days: float def __str__(self): return ( "Éléments orbitaux préliminaires (héliocentriques, écliptiques J2000 " "approx. -> voir note repère):\n" f" a (demi-grand axe) = {self.a_AU:.4f} UA\n" f" e (excentricité) = {self.e:.4f}\n" f" i (inclinaison) = {self.i_deg:.3f} deg\n" f" Omega (noeud ascendant) = {self.raan_deg:.3f} deg\n" f" omega (arg. périhélie) = {self.argp_deg:.3f} deg\n" f" nu (anomalie vraie) = {self.nu_deg:.3f} deg\n" f" M (anomalie moyenne) = {self.M_deg:.3f} deg\n" f" Période = {self.period_days:.2f} jours\n" f" Époque de référence = {self.epoch.isot}" ) # -------------------------------------------------------------------------- # Position héliocentrique de l'observateur # -------------------------------------------------------------------------- def observer_heliocentric_position_km(time: Time, location: EarthLocation | None = None) -> np.ndarray: """Position héliocentrique (repère équatorial ICRS) de l'observateur, en km. Si `location` (EarthLocation) est fourni, on ajoute l'offset topocentrique (site par rapport au centre de la Terre) ; sinon on prend le centre de la Terre (approximation géocentrique). """ earth_bary = get_body_barycentric("earth", time) # AU, ICRS, barycentre SSB sun_bary = get_body_barycentric("sun", time) # AU, ICRS earth_helio = (earth_bary - sun_bary).xyz.to(u.km).value # (3,) if location is not None: itrs = location.get_itrs(obstime=time) gcrs = itrs.transform_to(ITRS(obstime=time)) # déjà ITRS; conversion faite via location # offset géocentrique du site, en km, repère ITRS ~ proche de GCRS pour ce besoin site_geocentric = np.array([itrs.x.to(u.km).value, itrs.y.to(u.km).value, itrs.z.to(u.km).value]) earth_helio = earth_helio + site_geocentric return earth_helio # -------------------------------------------------------------------------- # Vecteurs unitaires de visée # -------------------------------------------------------------------------- def radec_to_unit_vector(ra_deg: float, dec_deg: float) -> np.ndarray: ra = np.deg2rad(ra_deg) dec = np.deg2rad(dec_deg) return np.array([ np.cos(dec) * np.cos(ra), np.cos(dec) * np.sin(ra), np.sin(dec), ]) # -------------------------------------------------------------------------- # Méthode de Gauss (algorithme de Curtis 5.10 / 5.11) # -------------------------------------------------------------------------- def _stumpff_like_lagrange(r2_norm, tau1, tau3, mu=MU_SUN): """Coefficients f, g approchés (ordre 4 en tau) pour Gibbs amélioré / correction de Gauss (utilisés pour affiner rho2 -> r2, r1, r3). """ f1 = 1 - 0.5 * mu * tau1**2 / r2_norm**3 f3 = 1 - 0.5 * mu * tau3**2 / r2_norm**3 g1 = tau1 - (1.0 / 6.0) * mu * tau1**3 / r2_norm**3 g3 = tau3 - (1.0 / 6.0) * mu * tau3**3 / r2_norm**3 return f1, f3, g1, g3 def gauss_preliminary_orbit( observations: list[Observation], location: EarthLocation | None = None, n_iter_refine: int = 5, ) -> OrbitalElements: """Détermine une orbite préliminaire à partir de 3 observations (RA, Dec, temps) via la méthode de Gauss, puis Gibbs pour la vitesse. Renvoie un objet OrbitalElements. Lève ValueError si aucune racine physique n'est trouvée. """ if len(observations) != 3: raise ValueError("La méthode de Gauss nécessite exactement 3 observations.") obs = sorted(observations, key=lambda o: o.time.jd) t1, t2, t3 = [o.time for o in obs] # temps en secondes, relatifs à t2 (repère TDB idéalement; UTC utilisé en approx) tau1 = (t1 - t2).sec tau3 = (t3 - t2).sec tau = tau3 - tau1 rho1_hat = radec_to_unit_vector(obs[0].ra_deg, obs[0].dec_deg) rho2_hat = radec_to_unit_vector(obs[1].ra_deg, obs[1].dec_deg) rho3_hat = radec_to_unit_vector(obs[2].ra_deg, obs[2].dec_deg) R1 = observer_heliocentric_position_km(t1, location) R2 = observer_heliocentric_position_km(t2, location) R3 = observer_heliocentric_position_km(t3, location) p1 = np.cross(rho2_hat, rho3_hat) p2 = np.cross(rho1_hat, rho3_hat) p3 = np.cross(rho1_hat, rho2_hat) D0 = np.dot(rho1_hat, p1) if abs(D0) < 1e-12: raise ValueError("Configuration géométrique dégénérée (D0 ~ 0) : " "observations trop proches ou alignées.") D = np.zeros((3, 3)) Rs = [R1, R2, R3] ps = [p1, p2, p3] for i in range(3): for j in range(3): D[i, j] = np.dot(Rs[i], ps[j]) A = (1.0 / D0) * (-D[0, 1] * (tau3 / tau) + D[1, 1] + D[2, 1] * (tau1 / tau)) B = (1.0 / (6.0 * D0)) * ( D[0, 1] * (tau3**2 - tau**2) * (tau3 / tau) + D[2, 1] * (tau**2 - tau1**2) * (tau1 / tau) ) E = np.dot(R2, rho2_hat) R2sq = np.dot(R2, R2) a_coef = -(A**2 + 2 * A * E + R2sq) b_coef = -2 * MU_SUN * B * (A + E) c_coef = -(MU_SUN**2) * (B**2) # polynôme en r2 : r2^8 + a_coef r2^6 + b_coef r2^3 + c_coef = 0 coeffs = [1, 0, a_coef, 0, 0, b_coef, 0, 0, c_coef] roots = np.roots(coeffs) real_positive = [r.real for r in roots if abs(r.imag) < 1e-6 and r.real > 0] if not real_positive: raise ValueError("Aucune racine réelle positive trouvée pour r2 : " "vérifiez les observations (temps, RA/Dec, unités).") # heuristique : on choisit la racine physiquement plausible # (typiquement entre 0.5 et 6 UA pour un astéroïde du système solaire interne/ceinture) candidates = sorted(real_positive) r2_norm = min(candidates, key=lambda r: abs(r - AU_KM * 2.0)) # a priori ~2 UA, ajustera si besoin # --- Estimation initiale (Curtis Algorithme 5.10) --- rho1 = (1.0 / D0) * ( (6 * (D[2, 0] * (tau1 / tau3) + D[1, 0] * (tau / tau3)) * r2_norm**3 + MU_SUN * D[2, 0] * (tau**2 - tau1**2) * (tau1 / tau3)) / (6 * r2_norm**3 + MU_SUN * (tau**2 - tau3**2)) - D[0, 0] ) rho2 = A + MU_SUN * B / r2_norm**3 rho3 = (1.0 / D0) * ( (6 * (D[0, 2] * (tau3 / tau1) - D[1, 2] * (tau / tau1)) * r2_norm**3 + MU_SUN * D[0, 2] * (tau**2 - tau3**2) * (tau3 / tau1)) / (6 * r2_norm**3 + MU_SUN * (tau**2 - tau1**2)) - D[2, 2] ) r1_vec = R1 + rho1 * rho1_hat r2_vec = R2 + rho2 * rho2_hat r3_vec = R3 + rho3 * rho3_hat # --- Raffinement itératif (Curtis Algorithme 5.11) : séries f,g --- for _ in range(n_iter_refine): r2_norm = np.linalg.norm(r2_vec) f1, f3, g1, g3 = _stumpff_like_lagrange(r2_norm, tau1, tau3) denom = f1 * g3 - f3 * g1 c1 = g3 / denom c3 = -g1 / denom rho1_new = (1.0 / D0) * (-D[0, 0] + D[1, 0] / c1 - c3 * D[2, 0] / c1) rho2_new = (1.0 / D0) * (-c1 * D[0, 1] + D[1, 1] - c3 * D[2, 1]) rho3_new = (1.0 / D0) * (-(c1 / c3) * D[0, 2] + D[1, 2] / c3 - D[2, 2]) r1_vec = R1 + rho1_new * rho1_hat r2_vec = R2 + rho2_new * rho2_hat r3_vec = R3 + rho3_new * rho3_hat if abs(rho2_new - rho2) < 1e-8 * abs(rho2_new): rho1, rho2, rho3 = rho1_new, rho2_new, rho3_new break rho1, rho2, rho3 = rho1_new, rho2_new, rho3_new # Vitesse au temps t2 par la méthode de Gibbs v2_vec = gibbs_velocity(r1_vec, r2_vec, r3_vec, mu=MU_SUN) # --- Garde-fou : solution dégénérée (rho2 quasi nul) --- # Si la distance topocentrique/géocentrique résolue est très petite par # rapport à la distance Terre-Soleil, r2_vec ≈ R2 : la méthode a # (numériquement) convergé vers "l'objet est confondu avec la Terre", # ce qui donne une orbite ≈ l'orbite terrestre (a≈1 UA, e≈0.0167, # période≈365 j). C'est le signe d'une géométrie d'observation trop # dégénérée pour Gauss (temps mal répartis entre les 3 observations, # arc trop court, ou bruit d'astrométrie excessif), pas d'un vrai # résultat. R2_norm = np.linalg.norm(R2) if abs(rho2) < 0.05 * R2_norm: import warnings warnings.warn( f"Solution dégénérée probable : rho2 = {rho2/AU_KM:.4f} UA (quasi nul). " "L'orbite obtenue risque de n'être qu'une réplique de l'orbite terrestre. " "Essayez d'autres triplets d'observations (mieux répartis dans le temps, " "idéalement avec tau1 et tau3 du même ordre de grandeur), ou ajoutez une " "3e nuit d'observation pour un arc plus robuste." ) # --- Conversion repère équatorial ICRS -> écliptique J2000 --- # Indispensable : sans cela, l'inclinaison calculée est mesurée par # rapport à l'équateur céleste (décalage ~obliquité, 23.44°) et non # par rapport à l'écliptique. r2_ecl = equatorial_to_ecliptic(r2_vec) v2_ecl = equatorial_to_ecliptic(v2_vec) return state_vector_to_elements(r2_ecl, v2_ecl, epoch=t2) def gibbs_velocity(r1: np.ndarray, r2: np.ndarray, r3: np.ndarray, mu: float = MU_SUN) -> np.ndarray: """Méthode de Gibbs : calcule le vecteur vitesse au temps de r2 à partir de trois vecteurs position coplanaires (r1, r2, r3), en km et km/s. """ r1n, r2n, r3n = (np.linalg.norm(v) for v in (r1, r2, r3)) Z12 = np.cross(r1, r2) Z23 = np.cross(r2, r3) Z31 = np.cross(r3, r1) # test de coplanarité (angle entre Z23 et r1) alpha_cop = 90 - np.degrees(np.arccos(np.dot(Z23, r1) / (np.linalg.norm(Z23) * r1n))) if abs(alpha_cop) > 3: import warnings warnings.warn(f"Test de coplanarité Gibbs = {alpha_cop:.2f} deg (>3 deg) : " "orbite peu fiable, envisager la méthode de Herrick-Gibbs " "pour un arc plus court, ou vérifier les observations.") N = r1n * Z23 + r2n * Z31 + r3n * Z12 Dd = Z12 + Z23 + Z31 S = (r2n - r3n) * r1 + (r3n - r1n) * r2 + (r1n - r2n) * r3 Nn = np.linalg.norm(N) Dn = np.linalg.norm(Dd) v2 = np.sqrt(mu / (Nn * Dn)) * (np.cross(Dd, r2) / r2n + S) return v2 # -------------------------------------------------------------------------- # Conversion vecteur d'état (r, v) -> éléments orbitaux classiques # -------------------------------------------------------------------------- def state_vector_to_elements(r: np.ndarray, v: np.ndarray, epoch: Time, mu: float = MU_SUN) -> OrbitalElements: r_norm = np.linalg.norm(r) v_norm = np.linalg.norm(v) h_vec = np.cross(r, v) h_norm = np.linalg.norm(h_vec) n_vec = np.cross([0, 0, 1], h_vec) n_norm = np.linalg.norm(n_vec) e_vec = (1.0 / mu) * ((v_norm**2 - mu / r_norm) * r - np.dot(r, v) * v) e = np.linalg.norm(e_vec) energy = v_norm**2 / 2 - mu / r_norm a = -mu / (2 * energy) i = np.arccos(np.clip(h_vec[2] / h_norm, -1, 1)) if n_norm > 1e-10: raan = np.arccos(np.clip(n_vec[0] / n_norm, -1, 1)) if n_vec[1] < 0: raan = 2 * np.pi - raan else: raan = 0.0 if n_norm > 1e-10 and e > 1e-10: argp = np.arccos(np.clip(np.dot(n_vec, e_vec) / (n_norm * e), -1, 1)) if e_vec[2] < 0: argp = 2 * np.pi - argp else: argp = 0.0 if e > 1e-10: nu = np.arccos(np.clip(np.dot(e_vec, r) / (e * r_norm), -1, 1)) if np.dot(r, v) < 0: nu = 2 * np.pi - nu else: nu = 0.0 E = 2 * np.arctan2(np.sqrt(1 - e) * np.sin(nu / 2), np.sqrt(1 + e) * np.cos(nu / 2)) M = E - e * np.sin(E) M = M % (2 * np.pi) period_s = 2 * np.pi * np.sqrt(a**3 / mu) return OrbitalElements( a_AU=a / AU_KM, e=e, i_deg=np.degrees(i), raan_deg=np.degrees(raan) % 360, argp_deg=np.degrees(argp) % 360, nu_deg=np.degrees(nu) % 360, M_deg=np.degrees(M) % 360, epoch=epoch, period_days=period_s / 86400.0, ) def select_observation_triplet(track: list[dict]) -> list[Observation]: """Sélectionne 3 observations pour la méthode de Gauss à partir d'une trajectoire mesurée sur plusieurs images (dicts {"time","ra","dec"}), en regroupant par nuit d'observation et en prenant un point représentatif (moyenne robuste) par nuit plutôt qu'un point brut isolé - ce qui réduit l'effet du bruit d'astrométrie image par image. - Si >= 3 nuits distinctes : prend la première, une nuit médiane, et la dernière (les plus espacées possible, pour maximiser l'arc tout en gardant une répartition raisonnable de tau1/tau3). - Si exactement 2 nuits : ATTENTION, c'est insuffisant pour une méthode de Gauss bien conditionnée (déséquilibre temporel inévitable entre "l'écart entre nuits" et "l'écart intra-nuit"). On fait de son mieux (moyenne nuit 1, premier et dernier point de la nuit 2) mais un avertissement est émis : ajoutez une 3e nuit pour un résultat fiable. - Si 1 seule nuit : lève une erreur (aucune méthode de Gauss raisonnable n'est possible sans plusieurs nuits). """ import warnings as _warnings # regroupe par date UTC (nuit d'observation) nights: dict[str, list[dict]] = {} for pt in track: night_key = pt["time"].isot[:10] # YYYY-MM-DD ; approximation simple nights.setdefault(night_key, []).append(pt) night_keys = sorted(nights.keys()) n_nights = len(night_keys) def _mean_point(pts): t_mean = Time(np.mean([p["time"].jd for p in pts]), format="jd", scale="utc") ra_mean = np.mean([p["ra"] for p in pts]) dec_mean = np.mean([p["dec"] for p in pts]) return Observation(time=t_mean, ra_deg=ra_mean, dec_deg=dec_mean) if n_nights < 2: raise ValueError( "Une seule nuit d'observation disponible : impossible de calculer " "une orbite préliminaire (il faut au moins 2, idéalement 3 nuits " "bien espacées)." ) if n_nights == 2: _warnings.warn( "Seulement 2 nuits d'observation disponibles : la méthode de Gauss " "sera mal conditionnée (déséquilibre temporel entre tau1 et tau3). " "Le résultat est à prendre avec beaucoup de précaution. Ajoutez une " "3e nuit d'observation, si possible espacée de plusieurs jours, " "pour un résultat nettement plus fiable." ) night1 = nights[night_keys[0]] night2 = nights[night_keys[1]] obs1 = _mean_point(night1) night2_sorted = sorted(night2, key=lambda p: p["time"].jd) obs2 = Observation(time=night2_sorted[0]["time"], ra_deg=night2_sorted[0]["ra"], dec_deg=night2_sorted[0]["dec"]) obs3 = Observation(time=night2_sorted[-1]["time"], ra_deg=night2_sorted[-1]["ra"], dec_deg=night2_sorted[-1]["dec"]) return [obs1, obs2, obs3] # >= 3 nuits : première, "médiane" (la plus proche du milieu temporel), dernière mid_key = night_keys[n_nights // 2] chosen_keys = [night_keys[0], mid_key, night_keys[-1]] return [_mean_point(nights[k]) for k in chosen_keys] def orbit_stability_check(track: list[dict], location: EarthLocation | None = None): """Évalue la SENSIBILITÉ du résultat de la méthode de Gauss au choix précis des observations, en essayant plusieurs triplets différents (premier / milieu / dernier point de chaque nuit disponible) plutôt qu'un seul triplet fixe. Utile en particulier avec seulement 2 nuits (cas mal conditionné, cf. avertissement de select_observation_triplet) : si les éléments orbitaux varient énormément d'un triplet à l'autre, c'est la preuve directe que le résultat n'est pas exploitable en l'état - pas besoin de deviner, la dispersion des essais le montre. Affiche un tableau des résultats et renvoie la liste des (description_triplet, OrbitalElements | None). """ import warnings as _warnings nights: dict[str, list[dict]] = {} for pt in track: night_key = pt["time"].isot[:10] nights.setdefault(night_key, []).append(pt) for k in nights: nights[k].sort(key=lambda p: p["time"].jd) night_keys = sorted(nights.keys()) if len(night_keys) < 2: raise ValueError("Il faut au moins 2 nuits pour ce test de robustesse.") def _candidate_points(pts): idxs = sorted(set([0, len(pts) // 2, len(pts) - 1])) return [pts[i] for i in idxs] if len(night_keys) == 2: n1_candidates = _candidate_points(nights[night_keys[0]]) n2_candidates = _candidate_points(nights[night_keys[1]]) combos = [] for p1 in n1_candidates: for a_idx in range(len(n2_candidates)): for b_idx in range(a_idx + 1, len(n2_candidates)): combos.append((p1, n2_candidates[a_idx], n2_candidates[b_idx])) else: # >= 3 nuits : on fait varier la nuit "du milieu" choisie combos = [] first_pts = _candidate_points(nights[night_keys[0]]) last_pts = _candidate_points(nights[night_keys[-1]]) for mid_key in night_keys[1:-1]: mid_pts = _candidate_points(nights[mid_key]) for p1 in first_pts[:1]: for pm in mid_pts: for p3 in last_pts[:1]: combos.append((p1, pm, p3)) results = [] print(f"\n{'Triplet (t1 / t2 / t3)':55s} {'a (UA)':>9s} {'e':>7s} {'i (deg)':>8s}") print("-" * 85) for p1, p2, p3 in combos: obs = [ Observation(time=p1["time"], ra_deg=p1["ra"], dec_deg=p1["dec"]), Observation(time=p2["time"], ra_deg=p2["ra"], dec_deg=p2["dec"]), Observation(time=p3["time"], ra_deg=p3["ra"], dec_deg=p3["dec"]), ] label = f"{p1['time'].isot[11:19]} / {p2['time'].isot[5:19]} / {p3['time'].isot[11:19]}" try: with _warnings.catch_warnings(): _warnings.simplefilter("ignore") el = gauss_preliminary_orbit(obs, location=location) print(f"{label:55s} {el.a_AU:9.3f} {el.e:7.3f} {el.i_deg:8.2f}") results.append((label, el)) except Exception as exc: print(f"{label:55s} échec ({exc})") results.append((label, None)) valid = [el for _, el in results if el is not None] if len(valid) >= 2: a_vals = np.array([el.a_AU for el in valid]) e_vals = np.array([el.e for el in valid]) print("-" * 85) print(f"Dispersion sur {len(valid)} triplets valides : " f"a = {a_vals.min():.2f}-{a_vals.max():.2f} UA (écart-type {a_vals.std():.2f}), " f"e = {e_vals.min():.3f}-{e_vals.max():.3f} (écart-type {e_vals.std():.3f})") if a_vals.std() / max(a_vals.mean(), 1e-6) > 0.2 or e_vals.std() > 0.1: print("-> forte dispersion : le résultat de la méthode de Gauss n'est PAS " "fiable avec les observations actuelles. Une 3e nuit d'observation, " "bien espacée dans le temps, est nécessaire pour un résultat exploitable.") return results if __name__ == "__main__": import argparse import csv parser = argparse.ArgumentParser( description="Détermination d'orbite préliminaire (méthode de Gauss) " "à partir d'un fichier CSV de 3+ observations (isot_time, ra_deg, dec_deg)." ) parser.add_argument("csv_file", help="CSV avec colonnes isot_time, ra_deg, dec_deg (ex: sortie de detect_asteroid.py)") parser.add_argument("--lon", type=float, default=None, help="longitude de l'observatoire (deg, est positif)") parser.add_argument("--lat", type=float, default=None, help="latitude de l'observatoire (deg)") parser.add_argument("--height", type=float, default=0.0, help="altitude de l'observatoire (m)") parser.add_argument("--stability-check", action="store_true", help="teste plusieurs triplets d'observations et affiche la dispersion " "des résultats (recommandé avec seulement 2 nuits)") args = parser.parse_args() rows = [] with open(args.csv_file) as f: reader = csv.DictReader(f) for row in reader: rows.append(row) if len(rows) < 3: raise SystemExit("Il faut au moins 3 observations dans le CSV.") track = [ {"time": Time(r["isot_time"], format="isot", scale="utc"), "ra": float(r["ra_deg"]), "dec": float(r["dec_deg"])} for r in rows ] observations = select_observation_triplet(track) print("Observations sélectionnées pour la méthode de Gauss :") for o in observations: print(f" {o.time.isot} RA={o.ra_deg:.6f} deg Dec={o.dec_deg:.6f} deg") location = None if args.lon is not None and args.lat is not None: location = EarthLocation(lon=args.lon * u.deg, lat=args.lat * u.deg, height=args.height * u.m) elements = gauss_preliminary_orbit(observations, location=location) print() print(elements) if args.stability_check: print("\n--- Test de robustesse : sensibilité au choix des observations ---") orbit_stability_check(track, location=location)