Ryan Malloy
22b272fd0c
162 → 174 SQL objects, 27 → 28 test suites, 3 new C source files. Features: - solar_elongation(body_id, ts): Sun-Earth-Planet angle [0,180] degrees - planet_phase(body_id, ts): illuminated disk fraction [0,1] - satellite_is_eclipsed/next_eclipse_entry/exit/eclipse_fraction: cylindrical shadow model (Vallado §5.3) for Earth shadow prediction - observing_night_quality(observer, ts): composite PL/pgSQL scoring based on astronomical darkness duration and Moon interference - moon_libration_longitude/latitude/position_angle/libration/subsolar_longitude: optical libration from Meeus (1998) Ch. 53 Refactored magnitude_funcs.c to extract shared compute_planet_geometry() used by magnitude, elongation, and phase — single VSOP87 evaluation per call. All 28 regression suites pass. Zero compiler warnings.
299 lines
8.2 KiB
C
299 lines
8.2 KiB
C
/*
|
|
* magnitude_funcs.c -- Planet magnitude, solar elongation, phase fraction
|
|
*
|
|
* Uses the Mallama & Hilton (2018) magnitude model with
|
|
* VSOP87 positions for distances and phase angles.
|
|
*
|
|
* Solar elongation and planet phase reuse the same Sun-Planet-Earth
|
|
* triangle geometry, factored into compute_planet_geometry().
|
|
*
|
|
* Reference: Mallama & Hilton, "Computing Apparent Planetary
|
|
* Magnitudes for The Astronomical Almanac", A&C vol. 25, 2018.
|
|
*/
|
|
|
|
#include "postgres.h"
|
|
#include "fmgr.h"
|
|
#include "utils/timestamp.h"
|
|
#include "types.h"
|
|
#include "astro_math.h"
|
|
#include "vsop87.h"
|
|
#include <math.h>
|
|
|
|
PG_FUNCTION_INFO_V1(planet_magnitude);
|
|
PG_FUNCTION_INFO_V1(solar_elongation);
|
|
PG_FUNCTION_INFO_V1(planet_phase);
|
|
|
|
|
|
/*
|
|
* Per-planet phase correction -- Mallama & Hilton (2018).
|
|
*
|
|
* 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.
|
|
*
|
|
* 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.
|
|
*/
|
|
|
|
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 */
|
|
};
|
|
|
|
|
|
/*
|
|
* Shared Sun-Planet-Earth geometry.
|
|
*
|
|
* Computes the three distances (r, delta, R) and the phase angle
|
|
* (Sun-Planet-Earth angle) from VSOP87 positions. Used by
|
|
* magnitude, elongation, and phase functions.
|
|
*/
|
|
typedef struct
|
|
{
|
|
double r; /* Sun-Planet distance (AU) */
|
|
double delta; /* Earth-Planet distance (AU) */
|
|
double R; /* Sun-Earth distance (AU) */
|
|
double i_deg; /* Phase angle, degrees (Sun-Planet-Earth vertex) */
|
|
} planet_geometry;
|
|
|
|
static void
|
|
compute_planet_geometry(int body_id, double jd, planet_geometry *geo)
|
|
{
|
|
double earth_xyz[6], planet_xyz[6];
|
|
double gv[3];
|
|
double cos_i;
|
|
int vsop_body = body_id - 1; /* pg_orrery 1-based -> VSOP87 0-based */
|
|
|
|
GetVsop87Coor(jd, 2, earth_xyz); /* Earth (VSOP87 body 2) */
|
|
GetVsop87Coor(jd, vsop_body, planet_xyz); /* target planet */
|
|
|
|
/* Heliocentric distance to planet */
|
|
geo->r = sqrt(planet_xyz[0] * planet_xyz[0] +
|
|
planet_xyz[1] * planet_xyz[1] +
|
|
planet_xyz[2] * planet_xyz[2]);
|
|
|
|
/* Geocentric vector and distance */
|
|
gv[0] = planet_xyz[0] - earth_xyz[0];
|
|
gv[1] = planet_xyz[1] - earth_xyz[1];
|
|
gv[2] = planet_xyz[2] - earth_xyz[2];
|
|
geo->delta = sqrt(gv[0] * gv[0] + gv[1] * gv[1] + gv[2] * gv[2]);
|
|
|
|
/* Sun-Earth distance */
|
|
geo->R = sqrt(earth_xyz[0] * earth_xyz[0] +
|
|
earth_xyz[1] * earth_xyz[1] +
|
|
earth_xyz[2] * earth_xyz[2]);
|
|
|
|
/* Phase angle via law of cosines: vertex at planet */
|
|
cos_i = (geo->r * geo->r + geo->delta * geo->delta - geo->R * geo->R)
|
|
/ (2.0 * geo->r * geo->delta);
|
|
if (cos_i > 1.0) cos_i = 1.0;
|
|
if (cos_i < -1.0) cos_i = -1.0;
|
|
geo->i_deg = acos(cos_i) * RAD_TO_DEG;
|
|
}
|
|
|
|
|
|
/*
|
|
* Validate planet body_id for magnitude/elongation/phase.
|
|
* Must be 1-8 (Mercury-Neptune), not 3 (Earth).
|
|
*/
|
|
static void
|
|
validate_planet_body_id(int body_id, const char *func_name)
|
|
{
|
|
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
|
|
ereport(ERROR,
|
|
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
|
errmsg("%s: body_id %d must be 1-8 (Mercury-Neptune)",
|
|
func_name, body_id)));
|
|
|
|
if (body_id == BODY_EARTH)
|
|
ereport(ERROR,
|
|
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
|
|
errmsg("%s: cannot compute for Earth from Earth",
|
|
func_name)));
|
|
}
|
|
|
|
|
|
/* ================================================================
|
|
* planet_magnitude(body_id int4, timestamptz) -> float8
|
|
*
|
|
* Returns the apparent visual magnitude of a planet as seen from
|
|
* Earth. Uses Mallama & Hilton (2018) magnitude model.
|
|
*
|
|
* Body IDs: 1=Mercury, ..., 8=Neptune (not Sun 0, Earth 3, or Moon 10)
|
|
*
|
|
* NOTE: Saturn magnitude does not account for ring tilt, which
|
|
* can vary the apparent magnitude by ~1.5 mag. The returned value
|
|
* is approximate for Saturn.
|
|
* ================================================================
|
|
*/
|
|
Datum
|
|
planet_magnitude(PG_FUNCTION_ARGS)
|
|
{
|
|
int32 body_id = PG_GETARG_INT32(0);
|
|
int64 ts = PG_GETARG_INT64(1);
|
|
double jd;
|
|
planet_geometry geo;
|
|
double V;
|
|
|
|
validate_planet_body_id(body_id, "planet_magnitude");
|
|
|
|
jd = timestamptz_to_jd(ts);
|
|
compute_planet_geometry(body_id, jd, &geo);
|
|
|
|
V = planet_v10[body_id]
|
|
+ 5.0 * log10(geo.r * geo.delta)
|
|
+ phase_correction(body_id, geo.i_deg);
|
|
|
|
PG_RETURN_FLOAT8(V);
|
|
}
|
|
|
|
|
|
/* ================================================================
|
|
* solar_elongation(body_id int4, timestamptz) -> float8
|
|
*
|
|
* Sun-Earth-Planet angle in degrees [0, 180].
|
|
* How far a planet appears from the Sun in the sky.
|
|
*
|
|
* Uses law of cosines with vertex at Earth:
|
|
* cos(elong) = (R^2 + delta^2 - r^2) / (2 * R * delta)
|
|
*
|
|
* Mercury max ~28 deg, Venus max ~47 deg.
|
|
* Superior planets can reach ~180 deg (opposition).
|
|
* ================================================================
|
|
*/
|
|
Datum
|
|
solar_elongation(PG_FUNCTION_ARGS)
|
|
{
|
|
int32 body_id = PG_GETARG_INT32(0);
|
|
int64 ts = PG_GETARG_INT64(1);
|
|
double jd;
|
|
planet_geometry geo;
|
|
double cos_elong, elong_deg;
|
|
|
|
validate_planet_body_id(body_id, "solar_elongation");
|
|
|
|
jd = timestamptz_to_jd(ts);
|
|
compute_planet_geometry(body_id, jd, &geo);
|
|
|
|
/* Law of cosines, vertex at Earth */
|
|
cos_elong = (geo.R * geo.R + geo.delta * geo.delta - geo.r * geo.r)
|
|
/ (2.0 * geo.R * geo.delta);
|
|
if (cos_elong > 1.0) cos_elong = 1.0;
|
|
if (cos_elong < -1.0) cos_elong = -1.0;
|
|
elong_deg = acos(cos_elong) * RAD_TO_DEG;
|
|
|
|
PG_RETURN_FLOAT8(elong_deg);
|
|
}
|
|
|
|
|
|
/* ================================================================
|
|
* planet_phase(body_id int4, timestamptz) -> float8
|
|
*
|
|
* Illuminated fraction of a planet's disk as seen from Earth [0, 1].
|
|
* k = (1 + cos(i)) / 2
|
|
* where i is the phase angle (Sun-Planet-Earth).
|
|
*
|
|
* Inner planets vary dramatically (Venus crescent).
|
|
* Outer planets are always near 1.0.
|
|
* ================================================================
|
|
*/
|
|
Datum
|
|
planet_phase(PG_FUNCTION_ARGS)
|
|
{
|
|
int32 body_id = PG_GETARG_INT32(0);
|
|
int64 ts = PG_GETARG_INT64(1);
|
|
double jd;
|
|
planet_geometry geo;
|
|
double i_rad, k;
|
|
|
|
validate_planet_body_id(body_id, "planet_phase");
|
|
|
|
jd = timestamptz_to_jd(ts);
|
|
compute_planet_geometry(body_id, jd, &geo);
|
|
|
|
i_rad = geo.i_deg * DEG_TO_RAD;
|
|
k = (1.0 + cos(i_rad)) / 2.0;
|
|
|
|
PG_RETURN_FLOAT8(k);
|
|
}
|