warehack.ing/pg_orrery
2
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
pg_orrery/src/de_funcs.c
Ryan Malloy dbc1f20a46 Add v0.20.0: Lagrange point SQL functions, DE variants, regression tests 
37 new SQL objects (188 → 225 total):
- Sun-planet L1-L5: heliocentric, observe, equatorial, distance (5 IMMUTABLE)
- Earth-Moon L1-L5: observe, equatorial via ELP2000-82B (2 IMMUTABLE)
- Planetary moon L1-L5: Galilean/Saturn/Uranus/Mars families (8 IMMUTABLE)
- Hill radius, zone radius, mass ratio, point name (5 IMMUTABLE)
- DE variants with VSOP87/ELP2000-82B fallback (17 STABLE)

All 31 regression tests pass. 210/210 standalone math tests pass.
2026-02-28 14:21:28 -07:00

2277 lines
72 KiB
C
Raw Blame History

/*
* de_funcs.c -- SQL-facing DE ephemeris function variants
*
* Each _de() function is a STABLE STRICT PARALLEL SAFE variant of an
* existing IMMUTABLE function. On any DE failure, falls back to the
* compiled-in VSOP87/ELP2000-82B equivalent with a NOTICE.
*
* The observation pipeline is identical:
* 1. Heliocentric ecliptic J2000 position (DE or fallback)
* 2. Geocentric ecliptic (subtract Earth's heliocentric)
* 3. observe_from_geocentric() -> topocentric az/el
*
* Constant chain of custody rule 7:
* Both target and Earth ALWAYS come from the same provider.
* If DE fails for the target, we don't use DE for Earth either.
*/
#include "postgres.h"
#include "fmgr.h"
#include "funcapi.h"
#include "catalog/pg_type.h"
#include "utils/builtins.h"
#include "utils/timestamp.h"
#include "types.h"
#include "astro_math.h"
#include "eph_provider.h"
#include "vsop87.h"
#include "elp82b.h"
#include "kepler.h"
#include "lambert.h"
#include "l12.h"
#include "tass17.h"
#include "gust86.h"
#include "marssat.h"
#include "lagrange.h"
#include <math.h>
#include <string.h>
/* Forward declarations */
PG_FUNCTION_INFO_V1(planet_heliocentric_de);
PG_FUNCTION_INFO_V1(planet_observe_de);
PG_FUNCTION_INFO_V1(sun_observe_de);
PG_FUNCTION_INFO_V1(moon_observe_de);
PG_FUNCTION_INFO_V1(lambert_transfer_de);
PG_FUNCTION_INFO_V1(lambert_c3_de);
PG_FUNCTION_INFO_V1(galilean_observe_de);
PG_FUNCTION_INFO_V1(saturn_moon_observe_de);
PG_FUNCTION_INFO_V1(uranus_moon_observe_de);
PG_FUNCTION_INFO_V1(mars_moon_observe_de);
PG_FUNCTION_INFO_V1(planet_equatorial_de);
PG_FUNCTION_INFO_V1(moon_equatorial_de);
PG_FUNCTION_INFO_V1(planet_observe_apparent_de);
PG_FUNCTION_INFO_V1(sun_observe_apparent_de);
PG_FUNCTION_INFO_V1(moon_observe_apparent_de);
PG_FUNCTION_INFO_V1(planet_equatorial_apparent_de);
PG_FUNCTION_INFO_V1(moon_equatorial_apparent_de);
PG_FUNCTION_INFO_V1(small_body_observe_apparent_de);
PG_FUNCTION_INFO_V1(galilean_equatorial_de);
PG_FUNCTION_INFO_V1(saturn_moon_equatorial_de);
PG_FUNCTION_INFO_V1(uranus_moon_equatorial_de);
PG_FUNCTION_INFO_V1(mars_moon_equatorial_de);
PG_FUNCTION_INFO_V1(pg_orrery_ephemeris_info);
/* Lagrange DE variants */
PG_FUNCTION_INFO_V1(lagrange_heliocentric_de);
PG_FUNCTION_INFO_V1(lagrange_observe_de);
PG_FUNCTION_INFO_V1(lagrange_equatorial_de);
PG_FUNCTION_INFO_V1(lagrange_distance_de);
PG_FUNCTION_INFO_V1(lagrange_distance_oe_de);
PG_FUNCTION_INFO_V1(lunar_lagrange_observe_de);
PG_FUNCTION_INFO_V1(lunar_lagrange_equatorial_de);
PG_FUNCTION_INFO_V1(galilean_lagrange_observe_de);
PG_FUNCTION_INFO_V1(galilean_lagrange_equatorial_de);
PG_FUNCTION_INFO_V1(saturn_moon_lagrange_observe_de);
PG_FUNCTION_INFO_V1(saturn_moon_lagrange_equatorial_de);
PG_FUNCTION_INFO_V1(uranus_moon_lagrange_observe_de);
PG_FUNCTION_INFO_V1(uranus_moon_lagrange_equatorial_de);
PG_FUNCTION_INFO_V1(mars_moon_lagrange_observe_de);
PG_FUNCTION_INFO_V1(mars_moon_lagrange_equatorial_de);
PG_FUNCTION_INFO_V1(hill_radius_de);
PG_FUNCTION_INFO_V1(lagrange_zone_radius_de);
/* ================================================================
* planet_heliocentric_de(body_id int, timestamptz) -> heliocentric
*
* DE variant of planet_heliocentric(). STABLE.
* Falls back to VSOP87 if DE is unavailable.
* ================================================================
*/
Datum
planet_heliocentric_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int64 ts = PG_GETARG_INT64(1);
double jd;
double xyz[6];
pg_heliocentric *result;
if (body_id == BODY_SUN)
{
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);
/* Try DE first */
if (!eph_de_planet(body_id, jd, xyz))
{
int vsop_body = body_id - 1;
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable for this query, falling back to VSOP87")));
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_de(body_id int, observer, timestamptz) -> topocentric
*
* DE variant of planet_observe(). STABLE.
* Rule 7: both planet and Earth from the same provider.
* ================================================================
*/
Datum
planet_observe_de(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];
pg_topocentric *result;
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("planet_observe_de: 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);
/* Try DE for both planet and Earth (rule 7: same provider) */
if (eph_de_planet(body_id, jd, planet_xyz) &&
eph_de_earth(jd, earth_xyz))
{
/* DE succeeded */
}
else
{
int vsop_body = body_id - 1;
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable for this query, falling back to VSOP87")));
GetVsop87Coor(jd, 2, earth_xyz);
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_de(observer, timestamptz) -> topocentric
*
* DE variant of sun_observe(). STABLE.
* ================================================================
*/
Datum
sun_observe_de(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);
/* Sun geocentric = -Earth_heliocentric */
if (eph_de_earth(jd, earth_xyz))
{
geo_ecl[0] = -earth_xyz[0];
geo_ecl[1] = -earth_xyz[1];
geo_ecl[2] = -earth_xyz[2];
}
else
{
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable, falling back to VSOP87")));
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_de(observer, timestamptz) -> topocentric
*
* DE variant of moon_observe(). STABLE.
* ================================================================
*/
Datum
moon_observe_de(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 */
if (!eph_de_moon(jd, moon_ecl))
{
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable, falling back to ELP2000-82B")));
GetElp82bCoor(jd, moon_ecl);
}
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_from_geocentric(moon_ecl, jd, obs, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* Lambert transfer functions with DE positions
* ================================================================
*/
/*
* Compute planet heliocentric velocity via numerical differentiation.
*
* use_de: if true, use DE positions; if false, use VSOP87.
* Must match the provider used for the corresponding position query
* (rule 7: same provider for position and velocity).
*/
static void
planet_velocity_de(int body_id, double jd, bool use_de, double vel[3])
{
double pos_fwd[6], pos_bwd[6];
double dt = 0.01; /* days */
if (use_de)
{
bool got_fwd = eph_de_planet(body_id, jd + dt, pos_fwd);
bool got_bwd = eph_de_planet(body_id, jd - dt, pos_bwd);
if (!got_fwd || !got_bwd)
{
/* DE boundary straddled — use VSOP87 for both to stay consistent */
int vsop_body = body_id - 1;
GetVsop87Coor(jd + dt, vsop_body, pos_fwd);
GetVsop87Coor(jd - dt, vsop_body, pos_bwd);
}
}
else
{
int vsop_body = body_id - 1;
GetVsop87Coor(jd + dt, vsop_body, pos_fwd);
GetVsop87Coor(jd - dt, vsop_body, pos_bwd);
}
vel[0] = (pos_fwd[0] - pos_bwd[0]) / (2.0 * dt);
vel[1] = (pos_fwd[1] - pos_bwd[1]) / (2.0 * dt);
vel[2] = (pos_fwd[2] - pos_bwd[2]) / (2.0 * dt);
}
Datum
lambert_transfer_de(PG_FUNCTION_ARGS)
{
int32 dep_body = PG_GETARG_INT32(0);
int32 arr_body = PG_GETARG_INT32(1);
int64 dep_ts = PG_GETARG_INT64(2);
int64 arr_ts = PG_GETARG_INT64(3);
double dep_jd, arr_jd, tof_days;
double r1[6], r2[6];
double v_planet_dep[3], v_planet_arr[3];
double v_inf_dep[3], v_inf_arr[3];
double v_inf_dep_mag, v_inf_arr_mag;
double c3_dep, c3_arr;
lambert_result lr;
TupleDesc tupdesc;
Datum values[6];
bool nulls[6];
HeapTuple tuple;
double au_per_day_to_km_per_s;
int k;
bool have_de;
if (dep_body < 1 || dep_body > 8)
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("lambert_transfer_de: dep_body_id %d must be 1-8", dep_body)));
if (arr_body < 1 || arr_body > 8)
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("lambert_transfer_de: arr_body_id %d must be 1-8", arr_body)));
if (dep_body == arr_body)
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("lambert_transfer_de: departure and arrival bodies must differ")));
dep_jd = timestamptz_to_jd(dep_ts);
arr_jd = timestamptz_to_jd(arr_ts);
tof_days = arr_jd - dep_jd;
if (tof_days <= 0.0)
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("lambert_transfer_de: arrival must be after departure")));
/* Try DE for both positions (rule 7: same provider) */
have_de = eph_de_planet(dep_body, dep_jd, r1) &&
eph_de_planet(arr_body, arr_jd, r2);
if (!have_de)
{
int dep_vsop = dep_body - 1;
int arr_vsop = arr_body - 1;
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable, falling back to VSOP87")));
GetVsop87Coor(dep_jd, dep_vsop, r1);
GetVsop87Coor(arr_jd, arr_vsop, r2);
}
if (!lambert_solve_uv(r1, r2, tof_days, GAUSS_K2, 1, &lr))
PG_RETURN_NULL();
/* Planet velocities (same provider as positions — rule 7) */
planet_velocity_de(dep_body, dep_jd, have_de, v_planet_dep);
planet_velocity_de(arr_body, arr_jd, have_de, v_planet_arr);
au_per_day_to_km_per_s = AU_KM / 86400.0;
for (k = 0; k < 3; k++) {
v_inf_dep[k] = (lr.v1[k] - v_planet_dep[k]) * au_per_day_to_km_per_s;
v_inf_arr[k] = (lr.v2[k] - v_planet_arr[k]) * au_per_day_to_km_per_s;
}
v_inf_dep_mag = sqrt(v_inf_dep[0]*v_inf_dep[0] +
v_inf_dep[1]*v_inf_dep[1] +
v_inf_dep[2]*v_inf_dep[2]);
v_inf_arr_mag = sqrt(v_inf_arr[0]*v_inf_arr[0] +
v_inf_arr[1]*v_inf_arr[1] +
v_inf_arr[2]*v_inf_arr[2]);
c3_dep = v_inf_dep_mag * v_inf_dep_mag;
c3_arr = v_inf_arr_mag * v_inf_arr_mag;
if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE)
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("function returning record called in context "
"that cannot accept type record")));
tupdesc = BlessTupleDesc(tupdesc);
memset(nulls, 0, sizeof(nulls));
values[0] = Float8GetDatum(c3_dep);
values[1] = Float8GetDatum(c3_arr);
values[2] = Float8GetDatum(v_inf_dep_mag);
values[3] = Float8GetDatum(v_inf_arr_mag);
values[4] = Float8GetDatum(tof_days);
values[5] = Float8GetDatum(lr.sma);
tuple = heap_form_tuple(tupdesc, values, nulls);
PG_RETURN_DATUM(HeapTupleGetDatum(tuple));
}
Datum
lambert_c3_de(PG_FUNCTION_ARGS)
{
int32 dep_body = PG_GETARG_INT32(0);
int32 arr_body = PG_GETARG_INT32(1);
int64 dep_ts = PG_GETARG_INT64(2);
int64 arr_ts = PG_GETARG_INT64(3);
double dep_jd, arr_jd, tof_days;
double r1[6], r2[6];
double v_planet_dep[3];
double v_inf_dep[3];
double c3_dep;
lambert_result lr;
double au_per_day_to_km_per_s;
int k;
bool have_de;
if (dep_body < 1 || dep_body > 8 || arr_body < 1 || arr_body > 8)
PG_RETURN_NULL();
if (dep_body == arr_body)
PG_RETURN_NULL();
dep_jd = timestamptz_to_jd(dep_ts);
arr_jd = timestamptz_to_jd(arr_ts);
tof_days = arr_jd - dep_jd;
if (tof_days <= 0.0)
PG_RETURN_NULL();
have_de = eph_de_planet(dep_body, dep_jd, r1) &&
eph_de_planet(arr_body, arr_jd, r2);
if (!have_de)
{
int dep_vsop = dep_body - 1;
int arr_vsop = arr_body - 1;
GetVsop87Coor(dep_jd, dep_vsop, r1);
GetVsop87Coor(arr_jd, arr_vsop, r2);
}
if (!lambert_solve_uv(r1, r2, tof_days, GAUSS_K2, 1, &lr))
PG_RETURN_NULL();
planet_velocity_de(dep_body, dep_jd, have_de, v_planet_dep);
au_per_day_to_km_per_s = AU_KM / 86400.0;
for (k = 0; k < 3; k++)
v_inf_dep[k] = (lr.v1[k] - v_planet_dep[k]) * au_per_day_to_km_per_s;
c3_dep = v_inf_dep[0]*v_inf_dep[0] +
v_inf_dep[1]*v_inf_dep[1] +
v_inf_dep[2]*v_inf_dep[2];
PG_RETURN_FLOAT8(c3_dep);
}
/* ================================================================
* Planetary moon observation with DE parent positions
*
* For each planetary moon, the moon-theory offset (L1.2, TASS17,
* GUST86, MarsSat) is computed relative to its parent planet.
* The parent's position comes from DE instead of VSOP87, giving
* sub-arcsecond accuracy for the parent while keeping the
* moon-theory accuracy for the relative offset.
* ================================================================
*/
/*
* Internal: observe a planetary moon using DE for the parent planet
* and Earth positions. Falls back to VSOP87 if DE is unavailable.
*
* moon_rel[3]: moon position relative to parent (ecliptic J2000, AU)
* parent_body_id: pg_orrery body ID of parent (5=Jupiter, 6=Saturn, etc.)
*/
static void
observe_moon_de(const double moon_rel[3], int parent_body_id,
double jd, const pg_observer *obs,
pg_topocentric *result)
{
double parent_xyz[6];
double earth_xyz[6];
double geo_ecl[3];
bool have_de;
/* Rule 7: both parent and Earth from same provider */
have_de = eph_de_planet(parent_body_id, jd, parent_xyz) &&
eph_de_earth(jd, earth_xyz);
if (!have_de)
{
int vsop_parent = parent_body_id - 1;
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable, falling back to VSOP87")));
GetVsop87Coor(jd, vsop_parent, parent_xyz);
GetVsop87Coor(jd, 2, earth_xyz); /* VSOP87 body 2 = Earth */
}
/* Moon geocentric = (parent + moon_relative) - Earth */
geo_ecl[0] = (parent_xyz[0] + moon_rel[0]) - earth_xyz[0];
geo_ecl[1] = (parent_xyz[1] + moon_rel[1]) - earth_xyz[1];
geo_ecl[2] = (parent_xyz[2] + moon_rel[2]) - earth_xyz[2];
observe_from_geocentric(geo_ecl, jd, obs, result);
}
Datum
galilean_observe_de(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 moon_xyz[3];
pg_topocentric *result;
if (body_id < L12_IO || body_id > L12_CALLISTO)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("galilean_observe_de: body_id %d must be 0-3 (Io-Callisto)",
body_id)));
jd = timestamptz_to_jd(ts);
GetL12Coor(jd, body_id, moon_xyz, NULL);
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_moon_de(moon_xyz, BODY_JUPITER, jd, obs, result);
PG_RETURN_POINTER(result);
}
Datum
saturn_moon_observe_de(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 moon_xyz[3];
pg_topocentric *result;
if (body_id < TASS17_MIMAS || body_id > TASS17_HYPERION)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("saturn_moon_observe_de: body_id %d must be 0-7 (Mimas-Hyperion)",
body_id)));
jd = timestamptz_to_jd(ts);
GetTass17Coor(jd, body_id, moon_xyz, NULL);
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_moon_de(moon_xyz, BODY_SATURN, jd, obs, result);
PG_RETURN_POINTER(result);
}
Datum
uranus_moon_observe_de(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 moon_xyz[3];
pg_topocentric *result;
if (body_id < GUST86_MIRANDA || body_id > GUST86_OBERON)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("uranus_moon_observe_de: body_id %d must be 0-4 (Miranda-Oberon)",
body_id)));
jd = timestamptz_to_jd(ts);
GetGust86Coor(jd, body_id, moon_xyz, NULL);
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_moon_de(moon_xyz, BODY_URANUS, jd, obs, result);
PG_RETURN_POINTER(result);
}
Datum
mars_moon_observe_de(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 moon_xyz[3];
pg_topocentric *result;
if (body_id < MARS_SAT_PHOBOS || body_id > MARS_SAT_DEIMOS)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("mars_moon_observe_de: body_id %d must be 0-1 (Phobos-Deimos)",
body_id)));
jd = timestamptz_to_jd(ts);
GetMarsSatCoor(jd, body_id, moon_xyz, NULL);
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_moon_de(moon_xyz, BODY_MARS, jd, obs, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* Planetary moon equatorial coordinates with DE parent positions
*
* Same DE/VSOP87 dispatch as observe_moon_de(), but returns
* geocentric RA/Dec instead of topocentric az/el.
* ================================================================
*/
static void
equatorial_moon_de(const double moon_rel[3], int parent_body_id,
double jd, pg_equatorial *result)
{
double parent_xyz[6];
double earth_xyz[6];
double geo_ecl[3];
bool have_de;
/* Rule 7: both parent and Earth from same provider */
have_de = eph_de_planet(parent_body_id, jd, parent_xyz) &&
eph_de_earth(jd, earth_xyz);
if (!have_de)
{
int vsop_parent = parent_body_id - 1;
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable, falling back to VSOP87")));
GetVsop87Coor(jd, vsop_parent, parent_xyz);
GetVsop87Coor(jd, 2, earth_xyz); /* VSOP87 body 2 = Earth */
}
/* Moon geocentric = (parent + moon_relative) - Earth */
geo_ecl[0] = (parent_xyz[0] + moon_rel[0]) - earth_xyz[0];
geo_ecl[1] = (parent_xyz[1] + moon_rel[1]) - earth_xyz[1];
geo_ecl[2] = (parent_xyz[2] + moon_rel[2]) - earth_xyz[2];
geocentric_to_equatorial(geo_ecl, jd, result);
}
Datum
galilean_equatorial_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int64 ts = PG_GETARG_INT64(1);
double jd;
double moon_xyz[3];
pg_equatorial *result;
if (body_id < L12_IO || body_id > L12_CALLISTO)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("galilean_equatorial_de: body_id %d must be 0-3 (Io-Callisto)",
body_id)));
jd = timestamptz_to_jd(ts);
GetL12Coor(jd, body_id, moon_xyz, NULL);
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
equatorial_moon_de(moon_xyz, BODY_JUPITER, jd, result);
PG_RETURN_POINTER(result);
}
Datum
saturn_moon_equatorial_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int64 ts = PG_GETARG_INT64(1);
double jd;
double moon_xyz[3];
pg_equatorial *result;
if (body_id < TASS17_MIMAS || body_id > TASS17_HYPERION)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("saturn_moon_equatorial_de: body_id %d must be 0-7 (Mimas-Hyperion)",
body_id)));
jd = timestamptz_to_jd(ts);
GetTass17Coor(jd, body_id, moon_xyz, NULL);
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
equatorial_moon_de(moon_xyz, BODY_SATURN, jd, result);
PG_RETURN_POINTER(result);
}
Datum
uranus_moon_equatorial_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int64 ts = PG_GETARG_INT64(1);
double jd;
double moon_xyz[3];
pg_equatorial *result;
if (body_id < GUST86_MIRANDA || body_id > GUST86_OBERON)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("uranus_moon_equatorial_de: body_id %d must be 0-4 (Miranda-Oberon)",
body_id)));
jd = timestamptz_to_jd(ts);
GetGust86Coor(jd, body_id, moon_xyz, NULL);
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
equatorial_moon_de(moon_xyz, BODY_URANUS, jd, result);
PG_RETURN_POINTER(result);
}
Datum
mars_moon_equatorial_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int64 ts = PG_GETARG_INT64(1);
double jd;
double moon_xyz[3];
pg_equatorial *result;
if (body_id < MARS_SAT_PHOBOS || body_id > MARS_SAT_DEIMOS)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("mars_moon_equatorial_de: body_id %d must be 0-1 (Phobos-Deimos)",
body_id)));
jd = timestamptz_to_jd(ts);
GetMarsSatCoor(jd, body_id, moon_xyz, NULL);
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
equatorial_moon_de(moon_xyz, BODY_MARS, jd, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* planet_equatorial_de(body_id int, timestamptz) -> equatorial
*
* DE variant of planet_equatorial(). STABLE.
* Rule 7: both planet and Earth from the same provider.
* ================================================================
*/
Datum
planet_equatorial_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int64 ts = PG_GETARG_INT64(1);
double jd;
double earth_xyz[6];
double planet_xyz[6];
double geo_ecl[3];
pg_equatorial *result;
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("planet_equatorial_de: 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);
/* Rule 7: both planet and Earth from same provider */
if (eph_de_planet(body_id, jd, planet_xyz) &&
eph_de_earth(jd, earth_xyz))
{
/* DE succeeded */
}
else
{
int vsop_body = body_id - 1;
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable for this query, falling back to VSOP87")));
GetVsop87Coor(jd, 2, earth_xyz);
GetVsop87Coor(jd, vsop_body, planet_xyz);
}
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_equatorial *) palloc(sizeof(pg_equatorial));
geocentric_to_equatorial(geo_ecl, jd, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* moon_equatorial_de(timestamptz) -> equatorial
*
* DE variant of moon_equatorial(). STABLE.
* ================================================================
*/
Datum
moon_equatorial_de(PG_FUNCTION_ARGS)
{
int64 ts = PG_GETARG_INT64(0);
double jd;
double moon_ecl[3];
pg_equatorial *result;
jd = timestamptz_to_jd(ts);
if (!eph_de_moon(jd, moon_ecl))
{
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable, falling back to ELP2000-82B")));
GetElp82bCoor(jd, moon_ecl);
}
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
geocentric_to_equatorial(moon_ecl, jd, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* Earth velocity via DE (numerical differentiation) or VSOP87 (analytic).
*
* DE path: central difference (earth(jd+dt) - earth(jd-dt)) / (2*dt)
* VSOP87 path: analytic velocity from earth_xyz[3..5]
*
* use_de: must match the provider used for position (rule 7).
* vel_ecl[3]: output Earth velocity in ecliptic J2000 (AU/day).
* ================================================================
*/
static void
earth_velocity_de(double jd, bool use_de, double vel_ecl[3])
{
if (use_de)
{
double pos_fwd[6], pos_bwd[6];
double dt = 0.01; /* days */
bool got_fwd = eph_de_earth(jd + dt, pos_fwd);
bool got_bwd = eph_de_earth(jd - dt, pos_bwd);
if (got_fwd && got_bwd)
{
vel_ecl[0] = (pos_fwd[0] - pos_bwd[0]) / (2.0 * dt);
vel_ecl[1] = (pos_fwd[1] - pos_bwd[1]) / (2.0 * dt);
vel_ecl[2] = (pos_fwd[2] - pos_bwd[2]) / (2.0 * dt);
return;
}
/* DE boundary straddled — fall through to VSOP87 */
}
{
double earth_xyz[6];
GetVsop87Coor(jd, 2, earth_xyz);
vel_ecl[0] = earth_xyz[3];
vel_ecl[1] = earth_xyz[4];
vel_ecl[2] = earth_xyz[5];
}
}
/* ================================================================
* planet_observe_apparent_de(body_id int, observer, timestamptz) -> topocentric
*
* DE variant of planet_observe_apparent(). STABLE.
* Light-time + annual aberration. Rule 7.
* ================================================================
*/
Datum
planet_observe_apparent_de(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];
double geo_dist, tau;
double vel_ecl[3];
pg_topocentric *result;
bool have_de;
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("planet_observe_apparent_de: 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);
/* Rule 7: both planet and Earth from same provider */
have_de = eph_de_planet(body_id, jd, planet_xyz) &&
eph_de_earth(jd, earth_xyz);
if (!have_de)
{
int vsop_body = body_id - 1;
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable for this query, falling back to VSOP87")));
GetVsop87Coor(jd, 2, earth_xyz);
GetVsop87Coor(jd, vsop_body, planet_xyz);
}
/* Geometric geocentric distance */
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];
geo_dist = sqrt(geo_ecl[0]*geo_ecl[0] + geo_ecl[1]*geo_ecl[1] + geo_ecl[2]*geo_ecl[2]);
/* Retarded planet position (same provider) */
tau = geo_dist / C_LIGHT_AU_DAY;
if (have_de)
eph_de_planet(body_id, jd - tau, planet_xyz);
else
GetVsop87Coor(jd - tau, body_id - 1, planet_xyz);
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];
/* Earth velocity for aberration */
earth_velocity_de(jd, have_de, vel_ecl);
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_from_geocentric_aberrated(geo_ecl, jd, obs, vel_ecl, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* sun_observe_apparent_de(observer, timestamptz) -> topocentric
*
* DE variant of sun_observe_apparent(). STABLE.
* ================================================================
*/
Datum
sun_observe_apparent_de(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];
double vel_ecl[3];
pg_topocentric *result;
bool have_de;
jd = timestamptz_to_jd(ts);
have_de = eph_de_earth(jd, earth_xyz);
if (!have_de)
{
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable, falling back to VSOP87")));
GetVsop87Coor(jd, 2, earth_xyz);
}
geo_ecl[0] = -earth_xyz[0];
geo_ecl[1] = -earth_xyz[1];
geo_ecl[2] = -earth_xyz[2];
earth_velocity_de(jd, have_de, vel_ecl);
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_from_geocentric_aberrated(geo_ecl, jd, obs, vel_ecl, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* moon_observe_apparent_de(observer, timestamptz) -> topocentric
*
* DE variant with light-time + aberration. STABLE.
* Moon position from DE, Earth velocity from DE (numerical) or VSOP87.
* ================================================================
*/
Datum
moon_observe_apparent_de(PG_FUNCTION_ARGS)
{
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
int64 ts = PG_GETARG_INT64(1);
double jd;
double moon_ecl[3];
double geo_dist, tau;
double vel_ecl[3];
pg_topocentric *result;
bool have_de;
jd = timestamptz_to_jd(ts);
have_de = eph_de_moon(jd, moon_ecl);
if (!have_de)
{
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable, falling back to ELP2000-82B")));
GetElp82bCoor(jd, moon_ecl);
}
/* Light-time correction */
geo_dist = sqrt(moon_ecl[0]*moon_ecl[0] + moon_ecl[1]*moon_ecl[1] + moon_ecl[2]*moon_ecl[2]);
tau = geo_dist / C_LIGHT_AU_DAY;
if (have_de)
eph_de_moon(jd - tau, moon_ecl);
else
GetElp82bCoor(jd - tau, moon_ecl);
/* Earth velocity for aberration (DE numerical or VSOP87 analytic) */
earth_velocity_de(jd, have_de, vel_ecl);
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_from_geocentric_aberrated(moon_ecl, jd, obs, vel_ecl, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* planet_equatorial_apparent_de(body_id int, timestamptz) -> equatorial
*
* DE variant with light-time + aberration. STABLE.
* ================================================================
*/
Datum
planet_equatorial_apparent_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int64 ts = PG_GETARG_INT64(1);
double jd;
double earth_xyz[6];
double planet_xyz[6];
double geo_ecl[3];
double geo_dist, tau;
double vel_ecl[3];
pg_equatorial *result;
bool have_de;
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("planet_equatorial_apparent_de: 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);
have_de = eph_de_planet(body_id, jd, planet_xyz) &&
eph_de_earth(jd, earth_xyz);
if (!have_de)
{
int vsop_body = body_id - 1;
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable for this query, falling back to VSOP87")));
GetVsop87Coor(jd, 2, earth_xyz);
GetVsop87Coor(jd, vsop_body, planet_xyz);
}
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];
geo_dist = sqrt(geo_ecl[0]*geo_ecl[0] + geo_ecl[1]*geo_ecl[1] + geo_ecl[2]*geo_ecl[2]);
tau = geo_dist / C_LIGHT_AU_DAY;
if (have_de)
eph_de_planet(body_id, jd - tau, planet_xyz);
else
GetVsop87Coor(jd - tau, body_id - 1, planet_xyz);
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];
earth_velocity_de(jd, have_de, vel_ecl);
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
geocentric_to_equatorial_aberrated(geo_ecl, jd, vel_ecl, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* moon_equatorial_apparent_de(timestamptz) -> equatorial
*
* DE variant with light-time + aberration. STABLE.
* ================================================================
*/
Datum
moon_equatorial_apparent_de(PG_FUNCTION_ARGS)
{
int64 ts = PG_GETARG_INT64(0);
double jd;
double moon_ecl[3];
double geo_dist, tau;
double vel_ecl[3];
pg_equatorial *result;
bool have_de;
jd = timestamptz_to_jd(ts);
have_de = eph_de_moon(jd, moon_ecl);
if (!have_de)
{
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable, falling back to ELP2000-82B")));
GetElp82bCoor(jd, moon_ecl);
}
geo_dist = sqrt(moon_ecl[0]*moon_ecl[0] + moon_ecl[1]*moon_ecl[1] + moon_ecl[2]*moon_ecl[2]);
tau = geo_dist / C_LIGHT_AU_DAY;
if (have_de)
eph_de_moon(jd - tau, moon_ecl);
else
GetElp82bCoor(jd - tau, moon_ecl);
earth_velocity_de(jd, have_de, vel_ecl);
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
geocentric_to_equatorial_aberrated(moon_ecl, jd, vel_ecl, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* small_body_observe_apparent_de(orbital_elements, observer, timestamptz) -> topocentric
*
* DE variant of small_body_observe_apparent(). Uses DE for Earth
* position (rule 7 satisfied: body is always Kepler, Earth from
* best available provider). STABLE.
* ================================================================
*/
Datum
small_body_observe_apparent_de(PG_FUNCTION_ARGS)
{
pg_orbital_elements *oe = (pg_orbital_elements *) PG_GETARG_POINTER(0);
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1);
int64 ts = PG_GETARG_INT64(2);
double jd;
double body_helio[3];
double earth_xyz[6];
double geo_ecl[3];
double geo_dist, tau;
double vel_ecl[3];
pg_topocentric *result;
bool have_de;
jd = timestamptz_to_jd(ts);
kepler_position(oe->q, oe->e, oe->inc, oe->arg_peri, oe->raan,
oe->tp, jd, body_helio);
have_de = eph_de_earth(jd, earth_xyz);
if (!have_de)
{
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable, falling back to VSOP87")));
GetVsop87Coor(jd, 2, earth_xyz);
}
geo_ecl[0] = body_helio[0] - earth_xyz[0];
geo_ecl[1] = body_helio[1] - earth_xyz[1];
geo_ecl[2] = body_helio[2] - earth_xyz[2];
geo_dist = sqrt(geo_ecl[0]*geo_ecl[0] + geo_ecl[1]*geo_ecl[1] + geo_ecl[2]*geo_ecl[2]);
tau = geo_dist / C_LIGHT_AU_DAY;
kepler_position(oe->q, oe->e, oe->inc, oe->arg_peri, oe->raan,
oe->tp, jd - tau, body_helio);
geo_ecl[0] = body_helio[0] - earth_xyz[0];
geo_ecl[1] = body_helio[1] - earth_xyz[1];
geo_ecl[2] = body_helio[2] - earth_xyz[2];
earth_velocity_de(jd, have_de, vel_ecl);
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_from_geocentric_aberrated(geo_ecl, jd, obs, vel_ecl, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* pg_orrery_ephemeris_info() -> RECORD
*
* Diagnostic function: returns current ephemeris provider status.
* STABLE (not STRICT — no args), PARALLEL SAFE.
* ================================================================
*/
Datum
pg_orrery_ephemeris_info(PG_FUNCTION_ARGS)
{
TupleDesc tupdesc;
Datum values[6];
bool nulls[6];
HeapTuple tuple;
if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE)
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("function returning record called in context "
"that cannot accept type record")));
tupdesc = BlessTupleDesc(tupdesc);
memset(nulls, 0, sizeof(nulls));
if (eph_de_available())
{
const char *path = eph_get_path();
values[0] = CStringGetTextDatum("JPL_DE");
values[1] = path ? CStringGetTextDatum(path) : CStringGetTextDatum("");
values[2] = Float8GetDatum(eph_de_start_jd());
values[3] = Float8GetDatum(eph_de_end_jd());
values[4] = Int32GetDatum(eph_de_version());
values[5] = Float8GetDatum(eph_de_au_km());
}
else
{
values[0] = CStringGetTextDatum("VSOP87");
values[1] = (Datum) 0;
values[2] = (Datum) 0;
values[3] = (Datum) 0;
values[4] = (Datum) 0;
nulls[1] = true; /* no file path */
nulls[2] = true; /* no start_jd */
nulls[3] = true; /* no end_jd */
nulls[4] = true; /* no version */
values[5] = Float8GetDatum((double)AU_KM);
}
tuple = heap_form_tuple(tupdesc, values, nulls);
PG_RETURN_DATUM(HeapTupleGetDatum(tuple));
}
/* ================================================================
* Lagrange point functions with DE ephemeris
*
* DE variants of the Lagrange functions in lagrange_funcs.c.
* Each uses DE for planet/Earth positions where possible,
* falling back to VSOP87/ELP2000-82B on any DE failure.
* Rule 7 always holds: both target and Earth from the same provider.
* ================================================================
*/
/*
* Validate body_id and point_id for Sun-planet Lagrange functions.
* Duplicated from lagrange_funcs.c (no cross-TU symbols).
*/
static void
validate_lagrange_args_de(int body_id, int point_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 (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("%s: point_id %d must be 1-5 (L1-L5)",
func_name, point_id)));
}
/*
* Compute Sun-planet Lagrange point using DE for planet position+velocity.
* Falls back to VSOP87 if DE unavailable (rule 7: consistent provider).
*
* Returns 0 on success, -1 on solver failure.
* Sets *used_de to indicate which provider was used.
*/
static int
compute_sun_planet_lagrange_de(int body_id, int point_id, double jd,
double result[3], bool *used_de)
{
double planet_xyz[6];
double sun[3] = {0.0, 0.0, 0.0};
double planet_vel[3];
double ratio, mu;
/* Try DE for planet position */
if (eph_de_planet(body_id, jd, planet_xyz))
{
/* Velocity via central difference (DE provides position only) */
double pos_fwd[6], pos_bwd[6];
double dt = 0.01; /* days */
bool got_fwd = eph_de_planet(body_id, jd + dt, pos_fwd);
bool got_bwd = eph_de_planet(body_id, jd - dt, pos_bwd);
if (got_fwd && got_bwd)
{
planet_vel[0] = (pos_fwd[0] - pos_bwd[0]) / (2.0 * dt);
planet_vel[1] = (pos_fwd[1] - pos_bwd[1]) / (2.0 * dt);
planet_vel[2] = (pos_fwd[2] - pos_bwd[2]) / (2.0 * dt);
*used_de = true;
}
else
{
/* DE boundary straddled — fall through to VSOP87 */
goto vsop87_fallback;
}
}
else
{
vsop87_fallback:
{
int vsop_body = body_id - 1;
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable for this query, falling back to VSOP87")));
GetVsop87Coor(jd, vsop_body, planet_xyz);
/* VSOP87 provides analytic velocity in xyz[3..5] */
planet_vel[0] = planet_xyz[3];
planet_vel[1] = planet_xyz[4];
planet_vel[2] = planet_xyz[5];
*used_de = false;
}
}
ratio = sun_planet_ratio(body_id);
mu = mu_from_ratio(ratio);
return lagrange_position(sun, planet_xyz, planet_vel, mu, point_id, result);
}
/*
* Compute Earth-Moon Lagrange point using DE for Moon position.
* Falls back to ELP2000-82B if DE unavailable.
*
* Moon velocity always via central difference (neither DE nor ELP
* provide direct velocity). Result is geocentric ecliptic J2000.
*/
static int
compute_lunar_lagrange_de(int point_id, double jd, double result_geo[3],
bool *used_de)
{
double moon_pos[3], moon_fwd[3], moon_bwd[3];
double moon_vel[3];
double earth[3] = {0.0, 0.0, 0.0};
double mu;
double dt = 0.001; /* days */
if (eph_de_moon(jd, moon_pos))
{
bool got_fwd = eph_de_moon(jd + dt, moon_fwd);
bool got_bwd = eph_de_moon(jd - dt, moon_bwd);
if (got_fwd && got_bwd)
{
*used_de = true;
}
else
{
/* DE boundary straddled */
goto elp_fallback;
}
}
else
{
elp_fallback:
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable, falling back to ELP2000-82B")));
GetElp82bCoor(jd, moon_pos);
GetElp82bCoor(jd + dt, moon_fwd);
GetElp82bCoor(jd - dt, moon_bwd);
*used_de = false;
}
moon_vel[0] = (moon_fwd[0] - moon_bwd[0]) / (2.0 * dt);
moon_vel[1] = (moon_fwd[1] - moon_bwd[1]) / (2.0 * dt);
moon_vel[2] = (moon_fwd[2] - moon_bwd[2]) / (2.0 * dt);
mu = mu_from_ratio(EARTH_MOON_EMRAT);
return lagrange_position(earth, moon_pos, moon_vel, mu, point_id,
result_geo);
}
/*
* Observe a planetary moon Lagrange point using DE for parent planet
* and Earth positions. Falls back to VSOP87 if DE unavailable.
*
* lp_rel[3]: L-point offset relative to parent (ecliptic J2000, AU)
* parent_body_id: pg_orrery body ID (5=Jupiter, 6=Saturn, etc.)
*/
static void
observe_pmoon_lagrange_de(const double lp_rel[3], int parent_body_id,
double jd, const pg_observer *obs,
pg_topocentric *result)
{
double parent_xyz[6];
double earth_xyz[6];
double geo_ecl[3];
bool have_de;
/* Rule 7: both parent and Earth from same provider */
have_de = eph_de_planet(parent_body_id, jd, parent_xyz) &&
eph_de_earth(jd, earth_xyz);
if (!have_de)
{
int vsop_parent = parent_body_id - 1;
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable, falling back to VSOP87")));
GetVsop87Coor(jd, vsop_parent, parent_xyz);
GetVsop87Coor(jd, 2, earth_xyz);
}
geo_ecl[0] = (parent_xyz[0] + lp_rel[0]) - earth_xyz[0];
geo_ecl[1] = (parent_xyz[1] + lp_rel[1]) - earth_xyz[1];
geo_ecl[2] = (parent_xyz[2] + lp_rel[2]) - earth_xyz[2];
observe_from_geocentric(geo_ecl, jd, obs, result);
}
static void
equatorial_pmoon_lagrange_de(const double lp_rel[3], int parent_body_id,
double jd, pg_equatorial *result)
{
double parent_xyz[6];
double earth_xyz[6];
double geo_ecl[3];
bool have_de;
have_de = eph_de_planet(parent_body_id, jd, parent_xyz) &&
eph_de_earth(jd, earth_xyz);
if (!have_de)
{
int vsop_parent = parent_body_id - 1;
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable, falling back to VSOP87")));
GetVsop87Coor(jd, vsop_parent, parent_xyz);
GetVsop87Coor(jd, 2, earth_xyz);
}
geo_ecl[0] = (parent_xyz[0] + lp_rel[0]) - earth_xyz[0];
geo_ecl[1] = (parent_xyz[1] + lp_rel[1]) - earth_xyz[1];
geo_ecl[2] = (parent_xyz[2] + lp_rel[2]) - earth_xyz[2];
geocentric_to_equatorial(geo_ecl, jd, result);
}
/*
* Compute planetary moon Lagrange offset (no cross-TU call).
* Duplicated from lagrange_funcs.c.
*/
static int
compute_planetary_moon_lagrange_de(const double moon_rel[3],
const double moon_vel[3],
char family, int moon_id,
int point_id,
double lp_rel[3])
{
double planet_origin[3] = {0.0, 0.0, 0.0};
double ratio, mu;
ratio = planet_moon_ratio(family, moon_id);
if (ratio < 0.0)
return -1;
mu = mu_from_ratio(ratio);
return lagrange_position(planet_origin, moon_rel, moon_vel, mu,
point_id, lp_rel);
}
/* ================================================================
* 1. lagrange_heliocentric_de(body_id, point_id, timestamptz)
* -> heliocentric
*
* DE variant of lagrange_heliocentric(). STABLE.
* ================================================================
*/
Datum
lagrange_heliocentric_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int32 point_id = PG_GETARG_INT32(1);
int64 ts = PG_GETARG_INT64(2);
double jd;
double lp[3];
bool used_de;
pg_heliocentric *result;
validate_lagrange_args_de(body_id, point_id, "lagrange_heliocentric_de");
jd = timestamptz_to_jd(ts);
if (compute_sun_planet_lagrange_de(body_id, point_id, jd, lp, &used_de) != 0)
ereport(ERROR,
(errcode(ERRCODE_INTERNAL_ERROR),
errmsg("lagrange_heliocentric_de: solver failed for body %d L%d",
body_id, point_id)));
result = (pg_heliocentric *) palloc(sizeof(pg_heliocentric));
result->x = lp[0];
result->y = lp[1];
result->z = lp[2];
PG_RETURN_POINTER(result);
}
/* ================================================================
* 2. lagrange_observe_de(body_id, point_id, observer, timestamptz)
* -> topocentric
*
* DE variant of lagrange_observe(). STABLE.
* Rule 7: L-point + Earth from same provider.
* ================================================================
*/
Datum
lagrange_observe_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int32 point_id = PG_GETARG_INT32(1);
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2);
int64 ts = PG_GETARG_INT64(3);
double jd;
double lp[3], earth_xyz[6], geo_ecl[3];
bool used_de;
pg_topocentric *result;
validate_lagrange_args_de(body_id, point_id, "lagrange_observe_de");
jd = timestamptz_to_jd(ts);
if (compute_sun_planet_lagrange_de(body_id, point_id, jd, lp, &used_de) != 0)
ereport(ERROR,
(errcode(ERRCODE_INTERNAL_ERROR),
errmsg("lagrange_observe_de: solver failed")));
/* Earth from same provider as the L-point (rule 7) */
if (used_de && eph_de_earth(jd, earth_xyz))
{
/* DE succeeded for both */
}
else
{
GetVsop87Coor(jd, 2, earth_xyz);
}
geo_ecl[0] = lp[0] - earth_xyz[0];
geo_ecl[1] = lp[1] - earth_xyz[1];
geo_ecl[2] = lp[2] - earth_xyz[2];
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_from_geocentric(geo_ecl, jd, obs, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* 3. lagrange_equatorial_de(body_id, point_id, timestamptz)
* -> equatorial
*
* DE variant of lagrange_equatorial(). STABLE.
* ================================================================
*/
Datum
lagrange_equatorial_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int32 point_id = PG_GETARG_INT32(1);
int64 ts = PG_GETARG_INT64(2);
double jd;
double lp[3], earth_xyz[6], geo_ecl[3];
bool used_de;
pg_equatorial *result;
validate_lagrange_args_de(body_id, point_id, "lagrange_equatorial_de");
jd = timestamptz_to_jd(ts);
if (compute_sun_planet_lagrange_de(body_id, point_id, jd, lp, &used_de) != 0)
ereport(ERROR,
(errcode(ERRCODE_INTERNAL_ERROR),
errmsg("lagrange_equatorial_de: solver failed")));
if (used_de && eph_de_earth(jd, earth_xyz))
{
/* DE succeeded */
}
else
{
GetVsop87Coor(jd, 2, earth_xyz);
}
geo_ecl[0] = lp[0] - earth_xyz[0];
geo_ecl[1] = lp[1] - earth_xyz[1];
geo_ecl[2] = lp[2] - earth_xyz[2];
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
geocentric_to_equatorial(geo_ecl, jd, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* 4. lagrange_distance_de(body_id, point_id, heliocentric, timestamptz)
* -> float8
*
* Distance from a heliocentric position to a Sun-planet L-point (AU).
* Uses DE for planet position. STABLE.
* ================================================================
*/
Datum
lagrange_distance_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int32 point_id = PG_GETARG_INT32(1);
pg_heliocentric *pos = (pg_heliocentric *) PG_GETARG_POINTER(2);
int64 ts = PG_GETARG_INT64(3);
double jd;
double lp[3];
double dx, dy, dz, dist;
bool used_de;
validate_lagrange_args_de(body_id, point_id, "lagrange_distance_de");
jd = timestamptz_to_jd(ts);
if (compute_sun_planet_lagrange_de(body_id, point_id, jd, lp, &used_de) != 0)
ereport(ERROR,
(errcode(ERRCODE_INTERNAL_ERROR),
errmsg("lagrange_distance_de: solver failed")));
dx = pos->x - lp[0];
dy = pos->y - lp[1];
dz = pos->z - lp[2];
dist = sqrt(dx*dx + dy*dy + dz*dz);
PG_RETURN_FLOAT8(dist);
}
/* ================================================================
* 5. lagrange_distance_oe_de(body_id, point_id, orbital_elements, timestamptz)
* -> float8
*
* Distance from an asteroid/comet to a Sun-planet L-point (AU).
* Uses DE for L-point computation. STABLE.
* ================================================================
*/
Datum
lagrange_distance_oe_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int32 point_id = PG_GETARG_INT32(1);
pg_orbital_elements *oe = (pg_orbital_elements *) PG_GETARG_POINTER(2);
int64 ts = PG_GETARG_INT64(3);
double jd;
double lp[3], body_pos[3];
double dx, dy, dz, dist;
bool used_de;
validate_lagrange_args_de(body_id, point_id, "lagrange_distance_oe_de");
jd = timestamptz_to_jd(ts);
if (compute_sun_planet_lagrange_de(body_id, point_id, jd, lp, &used_de) != 0)
ereport(ERROR,
(errcode(ERRCODE_INTERNAL_ERROR),
errmsg("lagrange_distance_oe_de: solver failed")));
kepler_position(oe->q, oe->e, oe->inc, oe->arg_peri, oe->raan,
oe->tp, jd, body_pos);
dx = body_pos[0] - lp[0];
dy = body_pos[1] - lp[1];
dz = body_pos[2] - lp[2];
dist = sqrt(dx*dx + dy*dy + dz*dz);
PG_RETURN_FLOAT8(dist);
}
/* ================================================================
* 6. lunar_lagrange_observe_de(point_id, observer, timestamptz)
* -> topocentric
*
* Earth-Moon L-point using DE Moon position. STABLE.
* ================================================================
*/
Datum
lunar_lagrange_observe_de(PG_FUNCTION_ARGS)
{
int32 point_id = PG_GETARG_INT32(0);
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1);
int64 ts = PG_GETARG_INT64(2);
double jd;
double lp_geo[3];
bool used_de;
pg_topocentric *result;
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("lunar_lagrange_observe_de: point_id %d must be 1-5",
point_id)));
jd = timestamptz_to_jd(ts);
if (compute_lunar_lagrange_de(point_id, jd, lp_geo, &used_de) != 0)
ereport(ERROR,
(errcode(ERRCODE_INTERNAL_ERROR),
errmsg("lunar_lagrange_observe_de: solver failed")));
/* lp_geo is already geocentric ecliptic J2000 (AU) */
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_from_geocentric(lp_geo, jd, obs, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* 7. lunar_lagrange_equatorial_de(point_id, timestamptz) -> equatorial
*
* Earth-Moon L-point using DE Moon position. STABLE.
* ================================================================
*/
Datum
lunar_lagrange_equatorial_de(PG_FUNCTION_ARGS)
{
int32 point_id = PG_GETARG_INT32(0);
int64 ts = PG_GETARG_INT64(1);
double jd;
double lp_geo[3];
bool used_de;
pg_equatorial *result;
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("lunar_lagrange_equatorial_de: point_id %d must be 1-5",
point_id)));
jd = timestamptz_to_jd(ts);
if (compute_lunar_lagrange_de(point_id, jd, lp_geo, &used_de) != 0)
ereport(ERROR,
(errcode(ERRCODE_INTERNAL_ERROR),
errmsg("lunar_lagrange_equatorial_de: solver failed")));
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
geocentric_to_equatorial(lp_geo, jd, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* Planetary moon Lagrange points with DE parent positions
*
* Moon-theory offset (L1.2, TASS17, GUST86, MarsSat) computed
* relative to parent, then Lagrange solver applied. Parent planet
* and Earth positions from DE (with VSOP87 fallback).
* ================================================================
*/
/* ── 8. Galilean Lagrange observe DE ─────────────────────── */
Datum
galilean_lagrange_observe_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int32 point_id = PG_GETARG_INT32(1);
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2);
int64 ts = PG_GETARG_INT64(3);
double jd;
double moon_xyz[3], moon_vel[3], lp_rel[3];
pg_topocentric *result;
if (body_id < L12_IO || body_id > L12_CALLISTO)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("galilean_lagrange_observe_de: body_id %d must be 0-3",
body_id)));
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("galilean_lagrange_observe_de: point_id %d must be 1-5",
point_id)));
jd = timestamptz_to_jd(ts);
GetL12Coor(jd, body_id, moon_xyz, moon_vel);
if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'g', body_id,
point_id, lp_rel) != 0)
ereport(ERROR,
(errcode(ERRCODE_INTERNAL_ERROR),
errmsg("galilean_lagrange_observe_de: solver failed")));
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_pmoon_lagrange_de(lp_rel, BODY_JUPITER, jd, obs, result);
PG_RETURN_POINTER(result);
}
/* ── 9. Galilean Lagrange equatorial DE ──────────────────── */
Datum
galilean_lagrange_equatorial_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int32 point_id = PG_GETARG_INT32(1);
int64 ts = PG_GETARG_INT64(2);
double jd;
double moon_xyz[3], moon_vel[3], lp_rel[3];
pg_equatorial *result;
if (body_id < L12_IO || body_id > L12_CALLISTO)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("galilean_lagrange_equatorial_de: body_id %d must be 0-3",
body_id)));
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("galilean_lagrange_equatorial_de: point_id %d must be 1-5",
point_id)));
jd = timestamptz_to_jd(ts);
GetL12Coor(jd, body_id, moon_xyz, moon_vel);
if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'g', body_id,
point_id, lp_rel) != 0)
ereport(ERROR,
(errcode(ERRCODE_INTERNAL_ERROR),
errmsg("galilean_lagrange_equatorial_de: solver failed")));
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
equatorial_pmoon_lagrange_de(lp_rel, BODY_JUPITER, jd, result);
PG_RETURN_POINTER(result);
}
/* ── 10. Saturn moon Lagrange observe DE ─────────────────── */
Datum
saturn_moon_lagrange_observe_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int32 point_id = PG_GETARG_INT32(1);
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2);
int64 ts = PG_GETARG_INT64(3);
double jd;
double moon_xyz[3], moon_vel[3], lp_rel[3];
pg_topocentric *result;
if (body_id < TASS17_MIMAS || body_id > TASS17_HYPERION)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("saturn_moon_lagrange_observe_de: body_id %d must be 0-7",
body_id)));
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("saturn_moon_lagrange_observe_de: point_id %d must be 1-5",
point_id)));
jd = timestamptz_to_jd(ts);
GetTass17Coor(jd, body_id, moon_xyz, moon_vel);
if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 's', body_id,
point_id, lp_rel) != 0)
ereport(ERROR,
(errcode(ERRCODE_INTERNAL_ERROR),
errmsg("saturn_moon_lagrange_observe_de: solver failed")));
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_pmoon_lagrange_de(lp_rel, BODY_SATURN, jd, obs, result);
PG_RETURN_POINTER(result);
}
/* ── 11. Saturn moon Lagrange equatorial DE ──────────────── */
Datum
saturn_moon_lagrange_equatorial_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int32 point_id = PG_GETARG_INT32(1);
int64 ts = PG_GETARG_INT64(2);
double jd;
double moon_xyz[3], moon_vel[3], lp_rel[3];
pg_equatorial *result;
if (body_id < TASS17_MIMAS || body_id > TASS17_HYPERION)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("saturn_moon_lagrange_equatorial_de: body_id %d must be 0-7",
body_id)));
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("saturn_moon_lagrange_equatorial_de: point_id %d must be 1-5",
point_id)));
jd = timestamptz_to_jd(ts);
GetTass17Coor(jd, body_id, moon_xyz, moon_vel);
if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 's', body_id,
point_id, lp_rel) != 0)
ereport(ERROR,
(errcode(ERRCODE_INTERNAL_ERROR),
errmsg("saturn_moon_lagrange_equatorial_de: solver failed")));
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
equatorial_pmoon_lagrange_de(lp_rel, BODY_SATURN, jd, result);
PG_RETURN_POINTER(result);
}
/* ── 12. Uranus moon Lagrange observe DE ─────────────────── */
Datum
uranus_moon_lagrange_observe_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int32 point_id = PG_GETARG_INT32(1);
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2);
int64 ts = PG_GETARG_INT64(3);
double jd;
double moon_xyz[3], moon_vel[3], lp_rel[3];
pg_topocentric *result;
if (body_id < GUST86_MIRANDA || body_id > GUST86_OBERON)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("uranus_moon_lagrange_observe_de: body_id %d must be 0-4",
body_id)));
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("uranus_moon_lagrange_observe_de: point_id %d must be 1-5",
point_id)));
jd = timestamptz_to_jd(ts);
GetGust86Coor(jd, body_id, moon_xyz, moon_vel);
if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'u', body_id,
point_id, lp_rel) != 0)
ereport(ERROR,
(errcode(ERRCODE_INTERNAL_ERROR),
errmsg("uranus_moon_lagrange_observe_de: solver failed")));
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_pmoon_lagrange_de(lp_rel, BODY_URANUS, jd, obs, result);
PG_RETURN_POINTER(result);
}
/* ── 13. Uranus moon Lagrange equatorial DE ──────────────── */
Datum
uranus_moon_lagrange_equatorial_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int32 point_id = PG_GETARG_INT32(1);
int64 ts = PG_GETARG_INT64(2);
double jd;
double moon_xyz[3], moon_vel[3], lp_rel[3];
pg_equatorial *result;
if (body_id < GUST86_MIRANDA || body_id > GUST86_OBERON)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("uranus_moon_lagrange_equatorial_de: body_id %d must be 0-4",
body_id)));
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("uranus_moon_lagrange_equatorial_de: point_id %d must be 1-5",
point_id)));
jd = timestamptz_to_jd(ts);
GetGust86Coor(jd, body_id, moon_xyz, moon_vel);
if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'u', body_id,
point_id, lp_rel) != 0)
ereport(ERROR,
(errcode(ERRCODE_INTERNAL_ERROR),
errmsg("uranus_moon_lagrange_equatorial_de: solver failed")));
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
equatorial_pmoon_lagrange_de(lp_rel, BODY_URANUS, jd, result);
PG_RETURN_POINTER(result);
}
/* ── 14. Mars moon Lagrange observe DE ───────────────────── */
Datum
mars_moon_lagrange_observe_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int32 point_id = PG_GETARG_INT32(1);
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2);
int64 ts = PG_GETARG_INT64(3);
double jd;
double moon_xyz[3], moon_vel[3], lp_rel[3];
pg_topocentric *result;
if (body_id < MARS_SAT_PHOBOS || body_id > MARS_SAT_DEIMOS)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("mars_moon_lagrange_observe_de: body_id %d must be 0-1",
body_id)));
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("mars_moon_lagrange_observe_de: point_id %d must be 1-5",
point_id)));
jd = timestamptz_to_jd(ts);
GetMarsSatCoor(jd, body_id, moon_xyz, moon_vel);
if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'm', body_id,
point_id, lp_rel) != 0)
ereport(ERROR,
(errcode(ERRCODE_INTERNAL_ERROR),
errmsg("mars_moon_lagrange_observe_de: solver failed")));
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_pmoon_lagrange_de(lp_rel, BODY_MARS, jd, obs, result);
PG_RETURN_POINTER(result);
}
/* ── 15. Mars moon Lagrange equatorial DE ────────────────── */
Datum
mars_moon_lagrange_equatorial_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int32 point_id = PG_GETARG_INT32(1);
int64 ts = PG_GETARG_INT64(2);
double jd;
double moon_xyz[3], moon_vel[3], lp_rel[3];
pg_equatorial *result;
if (body_id < MARS_SAT_PHOBOS || body_id > MARS_SAT_DEIMOS)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("mars_moon_lagrange_equatorial_de: body_id %d must be 0-1",
body_id)));
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("mars_moon_lagrange_equatorial_de: point_id %d must be 1-5",
point_id)));
jd = timestamptz_to_jd(ts);
GetMarsSatCoor(jd, body_id, moon_xyz, moon_vel);
if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'm', body_id,
point_id, lp_rel) != 0)
ereport(ERROR,
(errcode(ERRCODE_INTERNAL_ERROR),
errmsg("mars_moon_lagrange_equatorial_de: solver failed")));
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
equatorial_pmoon_lagrange_de(lp_rel, BODY_MARS, jd, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* 16. hill_radius_de(body_id, timestamptz) -> float8 (AU)
*
* Hill radius using DE for planet heliocentric distance. STABLE.
* ================================================================
*/
Datum
hill_radius_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int64 ts = PG_GETARG_INT64(1);
double jd;
double planet_xyz[6], sep, ratio, mu;
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("hill_radius_de: body_id %d must be 1-8", body_id)));
jd = timestamptz_to_jd(ts);
if (!eph_de_planet(body_id, jd, planet_xyz))
{
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable, falling back to VSOP87")));
GetVsop87Coor(jd, body_id - 1, planet_xyz);
}
sep = sqrt(planet_xyz[0]*planet_xyz[0] +
planet_xyz[1]*planet_xyz[1] +
planet_xyz[2]*planet_xyz[2]);
ratio = sun_planet_ratio(body_id);
mu = mu_from_ratio(ratio);
PG_RETURN_FLOAT8(lagrange_hill_radius(sep, mu));
}
/* ================================================================
* 17. lagrange_zone_radius_de(body_id, point_id, timestamptz)
* -> float8 (AU)
*
* Libration zone radius using DE for planet distance. STABLE.
* ================================================================
*/
Datum
lagrange_zone_radius_de(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int32 point_id = PG_GETARG_INT32(1);
int64 ts = PG_GETARG_INT64(2);
double jd;
double planet_xyz[6], sep, ratio, mu, zr;
validate_lagrange_args_de(body_id, point_id, "lagrange_zone_radius_de");
jd = timestamptz_to_jd(ts);
if (!eph_de_planet(body_id, jd, planet_xyz))
{
if (eph_de_available())
ereport(NOTICE,
(errmsg("DE ephemeris unavailable, falling back to VSOP87")));
GetVsop87Coor(jd, body_id - 1, planet_xyz);
}
sep = sqrt(planet_xyz[0]*planet_xyz[0] +
planet_xyz[1]*planet_xyz[1] +
planet_xyz[2]*planet_xyz[2]);
ratio = sun_planet_ratio(body_id);
mu = mu_from_ratio(ratio);
zr = lagrange_zone_radius(sep, mu, point_id);
PG_RETURN_FLOAT8(zr);
}
Reference in New Issue View Git Blame Copy Permalink