Ryan Malloy
0544a78276
Phase 1 — Stars, comets, Keplerian propagation: - star_observe() / star_observe_safe(): fixed star alt/az via IAU 1976 precession, equatorial-to-horizontal transform - kepler_propagate(): two-body Keplerian orbit propagation for elliptic, parabolic, and hyperbolic orbits - comet_observe(): observe comets/asteroids from orbital elements - heliocentric type: ecliptic J2000 position (x, y, z in AU) Phase 2 — VSOP87 planets, ELP82B Moon, Sun: - planet_heliocentric(): VSOP87 heliocentric ecliptic J2000 positions for Mercury through Neptune (Bretagnon & Francou, MIT) - planet_observe(): full observation pipeline for any planet - sun_observe(): Sun position from negated Earth VSOP87 - moon_observe(): ELP2000-82B lunar position (Chapront-Touzé, MIT) - Clean-room precession (IAU 2006) and sidereal time (IERS 2010) - elliptic_to_rectangular utility (Stellarium, MIT) All Stellarium extractions are MIT-licensed, thread-safe (static caching removed for PARALLEL SAFE), zero external data files. All 9 regression tests pass (90ms total).
237 lines
7.2 KiB
C
237 lines
7.2 KiB
C
/*
|
|
* planet_funcs.c -- VSOP87 planet and ELP82B Moon observation
|
|
*
|
|
* SQL functions for planetary/Sun/Moon position and observation.
|
|
* Wraps VSOP87 (Bretagnon & Francou) for 8 planets + Sun, and
|
|
* ELP2000-82B (Chapront-Touze) for the Moon.
|
|
*
|
|
* The observation pipeline:
|
|
* 1. Heliocentric ecliptic J2000 position (VSOP87 or ELP82B)
|
|
* 2. Geocentric ecliptic (subtract Earth's heliocentric)
|
|
* 3. Ecliptic -> equatorial J2000 (obliquity rotation)
|
|
* 4. Spherical coordinates (RA, Dec, distance)
|
|
* 5. Precession J2000 -> date (IAU 1976, Lieske)
|
|
* 6. Sidereal time -> hour angle -> az/el
|
|
* 7. Return topocentric result
|
|
*/
|
|
|
|
#include "postgres.h"
|
|
#include "fmgr.h"
|
|
#include "funcapi.h"
|
|
#include "utils/timestamp.h"
|
|
#include "types.h"
|
|
#include "astro_math.h"
|
|
#include "vsop87.h"
|
|
#include "elp82b.h"
|
|
#include <math.h>
|
|
|
|
PG_FUNCTION_INFO_V1(planet_heliocentric);
|
|
PG_FUNCTION_INFO_V1(planet_observe);
|
|
PG_FUNCTION_INFO_V1(sun_observe);
|
|
PG_FUNCTION_INFO_V1(moon_observe);
|
|
|
|
|
|
/* ================================================================
|
|
* Internal: geocentric observation pipeline
|
|
*
|
|
* Takes geocentric ecliptic J2000 position, observer location,
|
|
* and Julian date. Converts through equatorial, precesses to
|
|
* date, and computes topocentric az/el.
|
|
*
|
|
* geo_ecl_au[3] = geocentric ecliptic J2000 position in AU
|
|
* ================================================================
|
|
*/
|
|
static 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 -> date */
|
|
precess_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 */
|
|
}
|
|
|
|
|
|
/* ================================================================
|
|
* planet_heliocentric(body_id int, timestamptz) -> heliocentric
|
|
*
|
|
* Returns the heliocentric ecliptic J2000 position of a planet
|
|
* (or the Sun, as the barycenter proxy at (0,0,0)).
|
|
*
|
|
* Body IDs: 0=Sun, 1=Mercury, ..., 8=Neptune
|
|
* ================================================================
|
|
*/
|
|
Datum
|
|
planet_heliocentric(PG_FUNCTION_ARGS)
|
|
{
|
|
int32 body_id = PG_GETARG_INT32(0);
|
|
int64 ts = PG_GETARG_INT64(1);
|
|
double jd;
|
|
double xyz[6];
|
|
int vsop_body;
|
|
pg_heliocentric *result;
|
|
|
|
if (body_id == BODY_SUN)
|
|
{
|
|
/* Sun is at the origin of heliocentric coordinates */
|
|
result = (pg_heliocentric *) palloc(sizeof(pg_heliocentric));
|
|
result->x = 0.0;
|
|
result->y = 0.0;
|
|
result->z = 0.0;
|
|
PG_RETURN_POINTER(result);
|
|
}
|
|
|
|
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
|
|
ereport(ERROR,
|
|
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
|
errmsg("invalid body_id %d: must be 0 (Sun) or 1-8 (Mercury-Neptune)",
|
|
body_id)));
|
|
|
|
jd = timestamptz_to_jd(ts);
|
|
vsop_body = body_id - 1; /* pg_orbit 1-based -> VSOP87 0-based */
|
|
GetVsop87Coor(jd, vsop_body, xyz);
|
|
|
|
result = (pg_heliocentric *) palloc(sizeof(pg_heliocentric));
|
|
result->x = xyz[0];
|
|
result->y = xyz[1];
|
|
result->z = xyz[2];
|
|
|
|
PG_RETURN_POINTER(result);
|
|
}
|
|
|
|
|
|
/* ================================================================
|
|
* planet_observe(body_id int, observer, timestamptz) -> topocentric
|
|
*
|
|
* Full pipeline: VSOP87 position -> geocentric -> equatorial ->
|
|
* precession -> topocentric az/el.
|
|
*
|
|
* Body IDs: 1=Mercury, ..., 8=Neptune (not Sun or Moon)
|
|
* ================================================================
|
|
*/
|
|
Datum
|
|
planet_observe(PG_FUNCTION_ARGS)
|
|
{
|
|
int32 body_id = PG_GETARG_INT32(0);
|
|
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1);
|
|
int64 ts = PG_GETARG_INT64(2);
|
|
double jd;
|
|
double earth_xyz[6];
|
|
double planet_xyz[6];
|
|
double geo_ecl[3];
|
|
int vsop_body;
|
|
pg_topocentric *result;
|
|
|
|
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
|
|
ereport(ERROR,
|
|
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
|
errmsg("planet_observe: body_id %d must be 1-8 (Mercury-Neptune)",
|
|
body_id)));
|
|
|
|
if (body_id == BODY_EARTH)
|
|
ereport(ERROR,
|
|
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
|
|
errmsg("cannot observe Earth from Earth")));
|
|
|
|
jd = timestamptz_to_jd(ts);
|
|
|
|
/* Earth's heliocentric position */
|
|
GetVsop87Coor(jd, 2, earth_xyz); /* VSOP87 body 2 = Earth */
|
|
|
|
/* Target planet heliocentric position */
|
|
vsop_body = body_id - 1;
|
|
GetVsop87Coor(jd, vsop_body, planet_xyz);
|
|
|
|
/* Geocentric ecliptic = planet - Earth */
|
|
geo_ecl[0] = planet_xyz[0] - earth_xyz[0];
|
|
geo_ecl[1] = planet_xyz[1] - earth_xyz[1];
|
|
geo_ecl[2] = planet_xyz[2] - earth_xyz[2];
|
|
|
|
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
|
observe_from_geocentric(geo_ecl, jd, obs, result);
|
|
|
|
PG_RETURN_POINTER(result);
|
|
}
|
|
|
|
|
|
/* ================================================================
|
|
* sun_observe(observer, timestamptz) -> topocentric
|
|
*
|
|
* The Sun's geocentric position is the negation of Earth's
|
|
* heliocentric VSOP87 position.
|
|
* ================================================================
|
|
*/
|
|
Datum
|
|
sun_observe(PG_FUNCTION_ARGS)
|
|
{
|
|
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
|
|
int64 ts = PG_GETARG_INT64(1);
|
|
double jd;
|
|
double earth_xyz[6];
|
|
double geo_ecl[3];
|
|
pg_topocentric *result;
|
|
|
|
jd = timestamptz_to_jd(ts);
|
|
|
|
/* Earth heliocentric -> Sun geocentric by negation */
|
|
GetVsop87Coor(jd, 2, earth_xyz);
|
|
geo_ecl[0] = -earth_xyz[0];
|
|
geo_ecl[1] = -earth_xyz[1];
|
|
geo_ecl[2] = -earth_xyz[2];
|
|
|
|
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
|
observe_from_geocentric(geo_ecl, jd, obs, result);
|
|
|
|
PG_RETURN_POINTER(result);
|
|
}
|
|
|
|
|
|
/* ================================================================
|
|
* moon_observe(observer, timestamptz) -> topocentric
|
|
*
|
|
* ELP2000-82B returns geocentric ecliptic J2000 Moon position.
|
|
* No Earth subtraction needed.
|
|
* ================================================================
|
|
*/
|
|
Datum
|
|
moon_observe(PG_FUNCTION_ARGS)
|
|
{
|
|
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
|
|
int64 ts = PG_GETARG_INT64(1);
|
|
double jd;
|
|
double moon_ecl[3];
|
|
pg_topocentric *result;
|
|
|
|
jd = timestamptz_to_jd(ts);
|
|
|
|
/* Moon geocentric ecliptic J2000 in AU */
|
|
GetElp82bCoor(jd, moon_ecl);
|
|
|
|
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
|
observe_from_geocentric(moon_ecl, jd, obs, result);
|
|
|
|
PG_RETURN_POINTER(result);
|
|
}
|