diff --git a/src/magnitude_funcs.c b/src/magnitude_funcs.c index 77cbc07..59a9eda 100644 --- a/src/magnitude_funcs.c +++ b/src/magnitude_funcs.c @@ -20,35 +20,99 @@ PG_FUNCTION_INFO_V1(planet_magnitude); /* - * Planet magnitude parameters -- Mallama & Hilton (2018), simplified. + * Per-planet phase correction -- Mallama & Hilton (2018). * - * V(1,0) = absolute magnitude at r=1 AU, delta=1 AU, i=0 deg - * Phase corrections are polynomial fits to i (phase angle in degrees). + * Mercury uses a 6th-order polynomial (their Eq. 1). + * Venus and Mars are piecewise with different coefficients for + * small vs large phase angles. Jupiter is piecewise at 12 deg. + * Saturn, Uranus, Neptune use simpler models. * - * We use the linear+quadratic terms which are sufficient for - * phase angles encountered from Earth (typically <180 deg). - * - * Saturn caveat: visual magnitude depends heavily on ring tilt - * (can vary by ~1.5 mag). The simplified model here uses a fixed - * V(1,0) without ring correction. + * Saturn caveat: ring tilt contribution (their Eq. 10) requires + * saturnicentric sub-observer latitude, which we don't compute. + * We use the globe-only model (Eq. 11/12) — error up to ~1.5 mag. */ -typedef struct { - double v10; /* V(1,0) */ - double c1; /* coefficient for i */ - double c2; /* coefficient for i^2 */ - double c3; /* coefficient for i^3 (0 if unused) */ -} mag_params; -static const mag_params planet_mag[] = { - [0] = { 0, 0, 0, 0 }, /* Sun: unused placeholder */ - [1] = { -0.613, 6.328e-2, -1.6336e-3, 0 }, /* Mercury */ - [2] = { -4.384, 1.044e-3, 3.687e-4, 0 }, /* Venus */ - [3] = { 0, 0, 0, 0 }, /* Earth: unused */ - [4] = { -1.601, 2.267e-2, -1.302e-4, 0 }, /* Mars */ - [5] = { -9.395, 3.7e-4, 0, 0 }, /* Jupiter */ - [6] = { -8.95, 0, 0, 0 }, /* Saturn (ring tilt NOT modeled) */ - [7] = { -7.110, 0, 0, 0 }, /* Uranus */ - [8] = { -7.00, 0, 0, 0 }, /* Neptune */ +static double +phase_correction(int body_id, double i) +{ + double i2 = i * i; + + switch (body_id) + { + case 1: /* Mercury: 6th-order polynomial */ + return i * (6.3280e-02 + + i * (-1.6336e-03 + + i * (3.3644e-05 + + i * (-3.4265e-07 + + i * (1.6893e-09 + + i * (-3.0334e-12)))))); + + case 2: /* Venus: piecewise at 163.7 deg */ + if (i < 163.7) + return i * (-1.044e-03 + + i * (3.687e-04 + + i * (-2.814e-06 + + i * 8.938e-09))); + else + return (236.05828 + 4.384) + i * (-2.81914e+00 + + i * 8.39034e-03); + + case 4: /* Mars: piecewise at 50 deg */ + if (i <= 50.0) + return i * (2.267e-02 + i * (-1.302e-04)); + else + return (-1.601 + 0.367) + i * (-0.02573 + i * 0.0003445); + + case 5: /* Jupiter: piecewise at 12 deg */ + if (i <= 12.0) + return i * (6.16e-04 * i - 3.7e-04); + else + { + double a = i / 180.0; + return (-9.428 + 9.395) + (-2.5) + * log10(1.0 - 1.507 * a - 0.363 * a * a + - 0.062 * a * a * a + + 2.809 * a * a * a * a + - 1.876 * a * a * a * a * a); + } + + case 6: /* Saturn: globe-only (Eq. 11), no ring tilt */ + if (i <= 6.5) + return -3.7e-04 * i + 6.16e-04 * i2; + else + return 2.446e-04 * i + 2.672e-04 * i2 + - 1.506e-06 * i2 * i + 4.767e-09 * i2 * i2; + + case 7: /* Uranus */ + if (i <= 3.1) + return 0.0; + return i * (6.587e-03 + i * 1.045e-04); + + case 8: /* Neptune */ + if (i <= 1.9) + return 0.0; + return i * (7.944e-03 + i * 9.617e-05); + + default: + return 0.0; + } +} + + +/* + * V(1,0) per planet -- absolute magnitude at unit distances, zero phase. + * Mercury through Neptune. Mars piecewise handled in phase_correction(). + */ +static const double planet_v10[] = { + [0] = 0.0, /* Sun: unused */ + [1] = -0.613, /* Mercury */ + [2] = -4.384, /* Venus */ + [3] = 0.0, /* Earth: unused */ + [4] = -1.601, /* Mars (i <= 50; piecewise shifts in phase_correction) */ + [5] = -9.395, /* Jupiter (i <= 12; piecewise shifts in phase_correction) */ + [6] = -8.95, /* Saturn (globe-only) */ + [7] = -7.110, /* Uranus */ + [8] = -7.00, /* Neptune */ }; @@ -65,7 +129,6 @@ compute_planet_magnitude(int body_id, double jd) double geo[3]; double r, delta, R; double cos_i, i_deg; - const mag_params *p; double V; int vsop_body = body_id - 1; /* pg_orrery 1-based -> VSOP87 0-based */ @@ -94,13 +157,10 @@ compute_planet_magnitude(int body_id, double jd) if (cos_i < -1.0) cos_i = -1.0; i_deg = acos(cos_i) * RAD_TO_DEG; - /* Mallama & Hilton (2018) magnitude formula */ - p = &planet_mag[body_id]; - V = p->v10 + /* Mallama & Hilton (2018) magnitude with full phase correction */ + V = planet_v10[body_id] + 5.0 * log10(r * delta) - + p->c1 * i_deg - + p->c2 * i_deg * i_deg - + p->c3 * i_deg * i_deg * i_deg; + + phase_correction(body_id, i_deg); return V; }