/* * astro_math.h -- Shared astronomical math for pg_orrery * * Static inline functions used by star_funcs.c, kepler_funcs.c, * and future planet/moon observation code. * * Using static inline preserves the project convention of no * cross-translation-unit symbol coupling. */ #ifndef PG_ORRERY_ASTRO_MATH_H #define PG_ORRERY_ASTRO_MATH_H #include #include "types.h" #include "precession.h" #define DEG_TO_RAD (M_PI / 180.0) #define RAD_TO_DEG (180.0 / M_PI) #define ARCSEC_TO_RAD (M_PI / (180.0 * 3600.0)) /* Pre-computed obliquity trig (IAU 1976, OBLIQUITY_J2000 = 0.40909280422232897 rad). * Avoids recomputing cos/sin on every frame rotation call. */ #define COS_OBLIQUITY_J2000 0.91748206206918181 #define SIN_OBLIQUITY_J2000 0.39777715593191371 /* * Greenwich Mean Sidereal Time from Julian date. * Vallado, "Fundamentals of Astrodynamics", Eq. 3-47. * Returns angle in radians. */ static inline double gmst_from_jd(double jd) { double t_ut1 = (jd - 2451545.0) / 36525.0; double gmst = 67310.54841 + (876600.0 * 3600.0 + 8640184.812866) * t_ut1 + 0.093104 * t_ut1 * t_ut1 - 6.2e-6 * t_ut1 * t_ut1 * t_ut1; gmst = fmod(gmst * M_PI / 43200.0, 2.0 * M_PI); if (gmst < 0.0) gmst += 2.0 * M_PI; return gmst; } /* * IAU 1976 precession: rotate J2000 equatorial to mean equatorial of date. * * Source: Lieske (1979), A&A 73, 282. * Three angles zeta_A, z_A, theta_A define the precession rotation. * Valid within ~1 arcsecond for several centuries around J2000. */ static inline void precess_j2000_to_date(double jd, double ra_j2000, double dec_j2000, double *ra_date, double *dec_date) { double T, T2, T3; double zeta_A, z_A, theta_A; double cos_dec, sin_dec, cos_ra_zeta, sin_ra_zeta; double cos_theta, sin_theta; double A, B, C; T = (jd - J2000_JD) / 36525.0; T2 = T * T; T3 = T2 * T; /* Precession angles in arcseconds (Lieske 1979) */ zeta_A = (2306.2181 + 1.39656 * T - 0.000139 * T2) * T + (0.30188 - 0.000344 * T) * T2 + 0.017998 * T3; z_A = (2306.2181 + 1.39656 * T - 0.000139 * T2) * T + (1.09468 + 0.000066 * T) * T2 + 0.018203 * T3; theta_A = (2004.3109 - 0.85330 * T - 0.000217 * T2) * T - (0.42665 + 0.000217 * T) * T2 - 0.041833 * T3; /* Arcseconds to radians */ zeta_A *= ARCSEC_TO_RAD; z_A *= ARCSEC_TO_RAD; theta_A *= ARCSEC_TO_RAD; /* Direct formula: R3(-z_A) R2(theta_A) R3(-zeta_A) applied to (ra, dec) */ cos_dec = cos(dec_j2000); sin_dec = sin(dec_j2000); cos_ra_zeta = cos(ra_j2000 + zeta_A); sin_ra_zeta = sin(ra_j2000 + zeta_A); cos_theta = cos(theta_A); sin_theta = sin(theta_A); A = cos_dec * sin_ra_zeta; B = cos_theta * cos_dec * cos_ra_zeta - sin_theta * sin_dec; C = sin_theta * cos_dec * cos_ra_zeta + cos_theta * sin_dec; *dec_date = asin(C); *ra_date = atan2(A, B) + z_A; if (*ra_date < 0.0) *ra_date += 2.0 * M_PI; if (*ra_date >= 2.0 * M_PI) *ra_date -= 2.0 * M_PI; } /* * IAU 1976 precession + IAU 2000B nutation: J2000 -> true equatorial of date. * * Precesses to mean-of-date, then applies the dominant 4-term nutation * correction (Meeus 1998, Eq. 23.1). The combined correction reaches * ~17.2 arcsec in longitude (18.6-year period from the Moon's node). * * This is the correct transform for solar system and star observation * pipelines. Do NOT use for satellites in the TEME frame — SGP4 * already includes a simplified nutation model. */ static inline void precess_and_nutate_j2000_to_date(double jd, double ra_j2000, double dec_j2000, double *ra_date, double *dec_date) { double ra_mean, dec_mean; double dpsi, deps; /* nutation angles, arcseconds */ double eps_A, chi_A, omega_A, psi_A; /* precession angles, arcseconds */ double eps_rad; /* mean obliquity, radians */ double sin_eps, cos_eps; double sin_ra, cos_ra, tan_dec; /* Step 1: precess J2000 -> mean of date */ precess_j2000_to_date(jd, ra_j2000, dec_j2000, &ra_mean, &dec_mean); /* Step 2: compute nutation angles */ get_nutation_angles_iau2000b(jd, &dpsi, &deps); /* Step 3: mean obliquity of date (arcseconds -> radians) */ get_precession_angles_vondrak(jd, &eps_A, &chi_A, &omega_A, &psi_A); eps_rad = eps_A * ARCSEC_TO_RAD; /* Step 4: nutation correction to RA/Dec (Meeus 1998, Eq. 23.1) * Δα = (cos ε + sin ε sin α tan δ) Δψ − cos α tan δ Δε * Δδ = sin ε cos α Δψ + sin α Δε * dpsi, deps are in arcseconds — convert to radians for the shift. */ sin_eps = sin(eps_rad); cos_eps = cos(eps_rad); sin_ra = sin(ra_mean); cos_ra = cos(ra_mean); tan_dec = tan(dec_mean); *ra_date = ra_mean + (cos_eps + sin_eps * sin_ra * tan_dec) * dpsi * ARCSEC_TO_RAD - cos_ra * tan_dec * deps * ARCSEC_TO_RAD; *dec_date = dec_mean + sin_eps * cos_ra * dpsi * ARCSEC_TO_RAD + sin_ra * deps * ARCSEC_TO_RAD; /* Normalize RA to [0, 2pi) */ if (*ra_date < 0.0) *ra_date += 2.0 * M_PI; if (*ra_date >= 2.0 * M_PI) *ra_date -= 2.0 * M_PI; } /* * Equatorial (hour angle, declination) to horizontal (azimuth, elevation). * All angles in radians. * Azimuth convention: 0=N, pi/2=E, pi=S, 3*pi/2=W. */ static inline void equatorial_to_horizontal(double ha, double dec, double lat, double *az, double *el) { double sin_lat = sin(lat); double cos_lat = cos(lat); double sin_dec = sin(dec); double cos_dec = cos(dec); double cos_ha = cos(ha); double y, x; *el = asin(sin_lat * sin_dec + cos_lat * cos_dec * cos_ha); y = -cos_dec * sin(ha); x = cos_lat * sin_dec - sin_lat * cos_dec * cos_ha; *az = atan2(y, x); if (*az < 0.0) *az += 2.0 * M_PI; } /* * Ecliptic J2000 to equatorial J2000. * Simple rotation around X-axis by -obliquity. * * NOT safe for aliased (in-place) calls: ecl and equ must not overlap. * Writing equ[1] before reading ecl[1] for equ[2] produces wrong results * when ecl == equ. The vendored sgp4/sdp4.c has separate in-place versions * that use a temp variable; do not confuse the two. */ static inline void ecliptic_to_equatorial(const double *ecl, double *equ) { equ[0] = ecl[0]; equ[1] = ecl[1] * COS_OBLIQUITY_J2000 - ecl[2] * SIN_OBLIQUITY_J2000; equ[2] = ecl[1] * SIN_OBLIQUITY_J2000 + ecl[2] * COS_OBLIQUITY_J2000; } /* * Equatorial J2000 to ecliptic J2000. * Rotation around X-axis by +obliquity. * * NOT safe for aliased (in-place) calls: equ and ecl must not overlap. * Same ordering hazard as ecliptic_to_equatorial() above. */ static inline void equatorial_to_ecliptic(const double *equ, double *ecl) { ecl[0] = equ[0]; ecl[1] = equ[1] * COS_OBLIQUITY_J2000 + equ[2] * SIN_OBLIQUITY_J2000; ecl[2] = -equ[1] * SIN_OBLIQUITY_J2000 + equ[2] * COS_OBLIQUITY_J2000; } /* * Cartesian to spherical: (x, y, z) -> (ra, dec, dist). * ra in [0, 2*pi), dec in [-pi/2, pi/2], dist in same units as input. */ static inline void cartesian_to_spherical(const double *xyz, double *ra, double *dec, double *dist) { *dist = sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1] + xyz[2] * xyz[2]); *dec = asin(xyz[2] / *dist); *ra = atan2(xyz[1], xyz[0]); if (*ra < 0.0) *ra += 2.0 * M_PI; } /* * Geocentric observation pipeline (shared by all observation functions). * * Takes geocentric ecliptic J2000 position in AU, observer location, * and Julian date. Converts through equatorial, precesses to date, * and computes topocentric az/el. * * This is the canonical path: * ecliptic J2000 -> equatorial J2000 -> precess to date -> * sidereal time -> hour angle -> az/el */ static inline void observe_from_geocentric(const double geo_ecl_au[3], double jd, const pg_observer *obs, pg_topocentric *result) { double geo_equ[3]; double ra_j2000, dec_j2000, geo_dist; double ra_date, dec_date; double gmst_val, lst, ha; double az, el; /* Ecliptic J2000 -> equatorial J2000 */ ecliptic_to_equatorial(geo_ecl_au, geo_equ); /* Cartesian -> spherical */ cartesian_to_spherical(geo_equ, &ra_j2000, &dec_j2000, &geo_dist); /* Precess J2000 -> true of date (precession + nutation) */ precess_and_nutate_j2000_to_date(jd, ra_j2000, dec_j2000, &ra_date, &dec_date); /* Hour angle and az/el */ gmst_val = gmst_from_jd(jd); lst = gmst_val + obs->lon; ha = lst - ra_date; equatorial_to_horizontal(ha, dec_date, obs->lat, &az, &el); result->azimuth = az; result->elevation = el; result->range_km = geo_dist * AU_KM; result->range_rate = 0.0; /* no velocity computation yet */ } /* * Geocentric ecliptic J2000 (AU) -> equatorial RA/Dec of date. * * Captures the intermediate result that observe_from_geocentric() computes * and discards. Stops after precession -- no hour angle or az/el. * Distance is converted to km (AU_KM) to match topocentric.range_km. */ static inline void geocentric_to_equatorial(const double geo_ecl_au[3], double jd, pg_equatorial *result) { double geo_equ[3]; double ra_j2000, dec_j2000, geo_dist; double ra_date, dec_date; ecliptic_to_equatorial(geo_ecl_au, geo_equ); cartesian_to_spherical(geo_equ, &ra_j2000, &dec_j2000, &geo_dist); precess_and_nutate_j2000_to_date(jd, ra_j2000, dec_j2000, &ra_date, &dec_date); result->ra = ra_date; result->dec = dec_date; result->distance = geo_dist * AU_KM; } /* * Apply classical annual aberration to J2000 equatorial RA/Dec. * * The apparent direction of an object is displaced toward the apex of * Earth's motion by v/c (~20.5 arcseconds maximum). * * earth_vel_equ: Earth's velocity in equatorial J2000 (AU/day) * ra, dec: pointers to J2000 RA and Dec (radians), modified in-place * * Reference: Green (1985) "Spherical Astronomy", Eq. 10.22 */ static inline void apply_annual_aberration(const double earth_vel_equ[3], double *ra, double *dec) { double cos_ra = cos(*ra); double sin_ra = sin(*ra); double cos_dec = cos(*dec); double sin_dec = sin(*dec); double d_ra, d_dec; /* Near celestial poles, cos(dec) -> 0, dRA diverges. Skip. */ if (fabs(cos_dec) < 1e-12) return; d_ra = (-earth_vel_equ[0] * sin_ra + earth_vel_equ[1] * cos_ra) / (C_LIGHT_AU_DAY * cos_dec); d_dec = (-earth_vel_equ[0] * cos_ra * sin_dec - earth_vel_equ[1] * sin_ra * sin_dec + earth_vel_equ[2] * cos_dec) / C_LIGHT_AU_DAY; *ra += d_ra; *dec += d_dec; if (*ra < 0.0) *ra += 2.0 * M_PI; if (*ra >= 2.0 * M_PI) *ra -= 2.0 * M_PI; } /* * Geocentric observation pipeline with annual aberration. * * Same as observe_from_geocentric() plus aberration correction applied * in J2000 equatorial coordinates before precessing to date. * * earth_vel_ecl: Earth's velocity in ecliptic J2000 (AU/day) * — gets rotated to equatorial internally. */ static inline void observe_from_geocentric_aberrated(const double geo_ecl_au[3], double jd, const pg_observer *obs, const double earth_vel_ecl[3], pg_topocentric *result) { double geo_equ[3]; double vel_equ[3]; double ra_j2000, dec_j2000, geo_dist; double ra_date, dec_date; double gmst_val, lst, ha; double az, el; ecliptic_to_equatorial(geo_ecl_au, geo_equ); ecliptic_to_equatorial(earth_vel_ecl, vel_equ); cartesian_to_spherical(geo_equ, &ra_j2000, &dec_j2000, &geo_dist); apply_annual_aberration(vel_equ, &ra_j2000, &dec_j2000); precess_and_nutate_j2000_to_date(jd, ra_j2000, dec_j2000, &ra_date, &dec_date); gmst_val = gmst_from_jd(jd); lst = gmst_val + obs->lon; ha = lst - ra_date; equatorial_to_horizontal(ha, dec_date, obs->lat, &az, &el); result->azimuth = az; result->elevation = el; result->range_km = geo_dist * AU_KM; result->range_rate = 0.0; } /* * Geocentric ecliptic J2000 (AU) -> equatorial RA/Dec of date, with aberration. * * Same as geocentric_to_equatorial() plus annual aberration before precession. * earth_vel_ecl: Earth's velocity in ecliptic J2000 (AU/day). */ static inline void geocentric_to_equatorial_aberrated(const double geo_ecl_au[3], double jd, const double earth_vel_ecl[3], pg_equatorial *result) { double geo_equ[3]; double vel_equ[3]; double ra_j2000, dec_j2000, geo_dist; double ra_date, dec_date; ecliptic_to_equatorial(geo_ecl_au, geo_equ); ecliptic_to_equatorial(earth_vel_ecl, vel_equ); cartesian_to_spherical(geo_equ, &ra_j2000, &dec_j2000, &geo_dist); apply_annual_aberration(vel_equ, &ra_j2000, &dec_j2000); precess_and_nutate_j2000_to_date(jd, ra_j2000, dec_j2000, &ra_date, &dec_date); result->ra = ra_date; result->dec = dec_date; result->distance = geo_dist * AU_KM; } /* * TEME position (km) to equatorial RA/Dec, geocentric. * * TEME is "True Equator, Mean Equinox" -- approximately the mean * equatorial frame of date. For geocentric (no observer parallax), * RA/Dec is just the direction of the position vector in TEME. * Accuracy vs true-of-date: ~arcsecond (nutation residual). */ static inline void teme_to_equatorial_geo(const double pos_teme[3], pg_equatorial *result) { double rm = sqrt(pos_teme[0]*pos_teme[0] + pos_teme[1]*pos_teme[1] + pos_teme[2]*pos_teme[2]); result->dec = asin(pos_teme[2] / rm); result->ra = atan2(pos_teme[1], pos_teme[0]); if (result->ra < 0.0) result->ra += 2.0 * M_PI; result->distance = rm; } /* * TEME position (km) to equatorial RA/Dec, topocentric (observer-corrected). * * Subtracts observer position (via ECEF) to get the observer-relative * range vector in TEME, then computes RA/Dec from that vector. * For LEO satellites, the topocentric vs geocentric parallax is ~1 deg. * * Lifted from od_math.c:od_teme_to_radec() logic. */ static inline void teme_to_equatorial_topo(const double pos_teme[3], double gmst, const double obs_ecef[3], pg_equatorial *result) { double cg = cos(gmst), sg = sin(gmst); double pos_ecef[3], range_ecef[3], range_teme[3], rm; /* TEME -> ECEF */ pos_ecef[0] = cg * pos_teme[0] + sg * pos_teme[1]; pos_ecef[1] = -sg * pos_teme[0] + cg * pos_teme[1]; pos_ecef[2] = pos_teme[2]; /* Observer-relative range in ECEF */ range_ecef[0] = pos_ecef[0] - obs_ecef[0]; range_ecef[1] = pos_ecef[1] - obs_ecef[1]; range_ecef[2] = pos_ecef[2] - obs_ecef[2]; /* Back to TEME (inertial) for RA/Dec */ range_teme[0] = cg * range_ecef[0] - sg * range_ecef[1]; range_teme[1] = sg * range_ecef[0] + cg * range_ecef[1]; range_teme[2] = range_ecef[2]; rm = sqrt(range_teme[0]*range_teme[0] + range_teme[1]*range_teme[1] + range_teme[2]*range_teme[2]); result->dec = asin(range_teme[2] / rm); result->ra = atan2(range_teme[1], range_teme[0]); if (result->ra < 0.0) result->ra += 2.0 * M_PI; result->distance = rm; } #endif /* PG_ORRERY_ASTRO_MATH_H */