pg_orrery/src/rise_set_funcs.c
Ryan Malloy 8ca4383b2e v0.14.0: refracted planet/moon rise/set, constellation identification
Add 4 refracted rise/set functions completing the rise/set feature set:
- planet_next_rise/set_refracted: -0.569 deg threshold (refraction only,
  point source — even Jupiter at opposition is only 24 arcsec)
- moon_next_rise/set_refracted: -0.833 deg threshold (refraction +
  mean semidiameter, same as Sun)

Add IAU constellation identification from Roman (1987) CDS VI/42:
- 357 boundary segments covering all 88 constellations
- Precesses J2000 coordinates to B1875.0 epoch for lookup
- Two overloads: constellation(equatorial) and constellation(float8, float8)
- IMMUTABLE (compiled-in static data)

141 -> 147 SQL objects. 24 -> 25 regression suites. All 25 pass.
2026-02-25 17:02:08 -07:00

560 lines
18 KiB
C

/*
* rise_set_funcs.c -- Rise/set prediction for solar system bodies
*
* Adapts the satellite pass prediction bisection algorithm from
* pass_funcs.c for planets, Sun, and Moon. The core difference:
* elevation is computed via VSOP87/ELP82B -> observe_from_geocentric()
* instead of SGP4 propagation.
*
* Coarse scan at 60-second steps (planets move slowly compared to LEO
* satellites at 30s), then bisection to 0.1-second precision.
*
* Returns NULL if the body doesn't rise/set within the search window
* (circumpolar or perpetually below horizon at observer latitude).
*/
#include "postgres.h"
#include "fmgr.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_next_rise);
PG_FUNCTION_INFO_V1(planet_next_set);
PG_FUNCTION_INFO_V1(sun_next_rise);
PG_FUNCTION_INFO_V1(sun_next_set);
PG_FUNCTION_INFO_V1(moon_next_rise);
PG_FUNCTION_INFO_V1(moon_next_set);
PG_FUNCTION_INFO_V1(sun_next_rise_refracted);
PG_FUNCTION_INFO_V1(sun_next_set_refracted);
PG_FUNCTION_INFO_V1(planet_next_rise_refracted);
PG_FUNCTION_INFO_V1(planet_next_set_refracted);
PG_FUNCTION_INFO_V1(moon_next_rise_refracted);
PG_FUNCTION_INFO_V1(moon_next_set_refracted);
#define COARSE_STEP_JD (60.0 / 86400.0) /* 60 seconds */
#define BISECT_TOL_JD (0.1 / 86400.0) /* 0.1 second */
#define DEFAULT_WINDOW_DAYS 7.0
/* body_type encoding for the elevation helper */
#define BTYPE_PLANET 0
#define BTYPE_SUN 1
#define BTYPE_MOON 2
/*
* Standard almanac refraction correction for rise/set of Sun and Moon.
* The Sun/Moon are considered to rise/set when their geometric center
* is 0.833 degrees below the geometric horizon:
* 0.569 deg = atmospheric refraction at horizon (Bennett 1982)
* 0.266 deg = mean solar/lunar semidiameter
*/
#define SUN_MOON_REFRACTED_HORIZON_RAD (-0.01454) /* -0.833 deg */
/*
* Refraction-only horizon for point sources (planets).
* No semidiameter correction needed — even Jupiter at opposition
* subtends only ~24" (0.4 arcmin), negligible against 34' refraction.
* Error from treating planets as point sources: <1 second in time.
*/
#define REFRACTION_ONLY_HORIZON_RAD (-0.00993) /* -0.569 deg */
/* ----------------------------------------------------------------
* elevation_at_jd_body -- compute topocentric elevation for a body
*
* Returns geometric elevation in radians. No error return path --
* VSOP87/ELP82B always succeed for reasonable dates.
* ----------------------------------------------------------------
*/
static double
elevation_at_jd_body(int body_type, int body_id,
const pg_observer *obs, double jd)
{
double earth_xyz[6];
double target_xyz[6];
double geo_ecl[3];
pg_topocentric topo;
switch (body_type)
{
case BTYPE_PLANET:
{
int vsop_body = body_id - 1;
GetVsop87Coor(jd, 2, earth_xyz); /* Earth */
GetVsop87Coor(jd, vsop_body, target_xyz);
geo_ecl[0] = target_xyz[0] - earth_xyz[0];
geo_ecl[1] = target_xyz[1] - earth_xyz[1];
geo_ecl[2] = target_xyz[2] - earth_xyz[2];
break;
}
case BTYPE_SUN:
{
GetVsop87Coor(jd, 2, earth_xyz);
geo_ecl[0] = -earth_xyz[0];
geo_ecl[1] = -earth_xyz[1];
geo_ecl[2] = -earth_xyz[2];
break;
}
case BTYPE_MOON:
{
double moon_ecl[3];
GetElp82bCoor(jd, moon_ecl);
geo_ecl[0] = moon_ecl[0];
geo_ecl[1] = moon_ecl[1];
geo_ecl[2] = moon_ecl[2];
break;
}
default:
return -M_PI; /* unreachable */
}
observe_from_geocentric(geo_ecl, jd, obs, &topo);
return topo.elevation;
}
/* ----------------------------------------------------------------
* find_next_crossing -- coarse scan + bisection for horizon crossing
*
* Scans from start_jd to stop_jd looking for the next rising or
* setting event. Returns the Julian date of the crossing, or -1
* if no crossing is found within the window.
*
* rising=true: find where elevation crosses threshold upward
* rising=false: find where elevation crosses threshold downward
* ----------------------------------------------------------------
*/
static double
find_next_crossing(int body_type, int body_id,
const pg_observer *obs,
double start_jd, double stop_jd,
double threshold_rad,
bool rising)
{
double jd = start_jd;
double prev_el, curr_el;
prev_el = elevation_at_jd_body(body_type, body_id, obs, jd);
while (jd < stop_jd)
{
jd += COARSE_STEP_JD;
if (jd > stop_jd)
jd = stop_jd;
curr_el = elevation_at_jd_body(body_type, body_id, obs, jd);
if (rising)
{
/* Rising: was below threshold, now above */
if (prev_el <= threshold_rad && curr_el > threshold_rad)
{
double lo = jd - COARSE_STEP_JD;
double hi = jd;
while (hi - lo > BISECT_TOL_JD)
{
double mid = (lo + hi) / 2.0;
if (elevation_at_jd_body(body_type, body_id, obs, mid) > threshold_rad)
hi = mid;
else
lo = mid;
}
return (lo + hi) / 2.0;
}
}
else
{
/* Setting: was above threshold, now below */
if (prev_el > threshold_rad && curr_el <= threshold_rad)
{
double lo = jd - COARSE_STEP_JD;
double hi = jd;
while (hi - lo > BISECT_TOL_JD)
{
double mid = (lo + hi) / 2.0;
if (elevation_at_jd_body(body_type, body_id, obs, mid) > threshold_rad)
lo = mid;
else
hi = mid;
}
return (lo + hi) / 2.0;
}
}
prev_el = curr_el;
}
return -1.0; /* no crossing found */
}
/* ================================================================
* planet_next_rise(body_id, observer, timestamptz) -> timestamptz
*
* Returns the next time a planet rises above the geometric horizon.
* NULL if the planet doesn't rise within 7 days (circumpolar or
* perpetually below horizon).
*
* Body IDs: 1=Mercury, ..., 8=Neptune (not Sun, Earth, or Moon)
* ================================================================
*/
Datum
planet_next_rise(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 start_jd, stop_jd, result_jd;
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("planet_next_rise: 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")));
start_jd = timestamptz_to_jd(ts);
stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
result_jd = find_next_crossing(BTYPE_PLANET, body_id, obs,
start_jd, stop_jd, 0.0, true);
if (result_jd < 0.0)
PG_RETURN_NULL();
PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
}
/* ================================================================
* planet_next_set(body_id, observer, timestamptz) -> timestamptz
* ================================================================
*/
Datum
planet_next_set(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 start_jd, stop_jd, result_jd;
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("planet_next_set: 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")));
start_jd = timestamptz_to_jd(ts);
stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
result_jd = find_next_crossing(BTYPE_PLANET, body_id, obs,
start_jd, stop_jd, 0.0, false);
if (result_jd < 0.0)
PG_RETURN_NULL();
PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
}
/* ================================================================
* sun_next_rise(observer, timestamptz) -> timestamptz
* ================================================================
*/
Datum
sun_next_rise(PG_FUNCTION_ARGS)
{
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
int64 ts = PG_GETARG_INT64(1);
double start_jd, stop_jd, result_jd;
start_jd = timestamptz_to_jd(ts);
stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
result_jd = find_next_crossing(BTYPE_SUN, 0, obs,
start_jd, stop_jd, 0.0, true);
if (result_jd < 0.0)
PG_RETURN_NULL();
PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
}
/* ================================================================
* sun_next_set(observer, timestamptz) -> timestamptz
* ================================================================
*/
Datum
sun_next_set(PG_FUNCTION_ARGS)
{
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
int64 ts = PG_GETARG_INT64(1);
double start_jd, stop_jd, result_jd;
start_jd = timestamptz_to_jd(ts);
stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
result_jd = find_next_crossing(BTYPE_SUN, 0, obs,
start_jd, stop_jd, 0.0, false);
if (result_jd < 0.0)
PG_RETURN_NULL();
PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
}
/* ================================================================
* moon_next_rise(observer, timestamptz) -> timestamptz
* ================================================================
*/
Datum
moon_next_rise(PG_FUNCTION_ARGS)
{
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
int64 ts = PG_GETARG_INT64(1);
double start_jd, stop_jd, result_jd;
start_jd = timestamptz_to_jd(ts);
stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
result_jd = find_next_crossing(BTYPE_MOON, 0, obs,
start_jd, stop_jd, 0.0, true);
if (result_jd < 0.0)
PG_RETURN_NULL();
PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
}
/* ================================================================
* moon_next_set(observer, timestamptz) -> timestamptz
* ================================================================
*/
Datum
moon_next_set(PG_FUNCTION_ARGS)
{
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
int64 ts = PG_GETARG_INT64(1);
double start_jd, stop_jd, result_jd;
start_jd = timestamptz_to_jd(ts);
stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
result_jd = find_next_crossing(BTYPE_MOON, 0, obs,
start_jd, stop_jd, 0.0, false);
if (result_jd < 0.0)
PG_RETURN_NULL();
PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
}
/* ================================================================
* sun_next_rise_refracted(observer, timestamptz) -> timestamptz
*
* Uses -0.833 degree threshold (standard almanac: 0.569 deg refraction
* at horizon + 0.266 deg solar semidiameter). Refracted sunrise is
* earlier than geometric by ~4 minutes at mid-latitudes.
* ================================================================
*/
Datum
sun_next_rise_refracted(PG_FUNCTION_ARGS)
{
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
int64 ts = PG_GETARG_INT64(1);
double start_jd, stop_jd, result_jd;
start_jd = timestamptz_to_jd(ts);
stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
result_jd = find_next_crossing(BTYPE_SUN, 0, obs,
start_jd, stop_jd,
SUN_MOON_REFRACTED_HORIZON_RAD, true);
if (result_jd < 0.0)
PG_RETURN_NULL();
PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
}
/* ================================================================
* sun_next_set_refracted(observer, timestamptz) -> timestamptz
*
* Refracted sunset is later than geometric by ~4 minutes.
* ================================================================
*/
Datum
sun_next_set_refracted(PG_FUNCTION_ARGS)
{
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
int64 ts = PG_GETARG_INT64(1);
double start_jd, stop_jd, result_jd;
start_jd = timestamptz_to_jd(ts);
stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
result_jd = find_next_crossing(BTYPE_SUN, 0, obs,
start_jd, stop_jd,
SUN_MOON_REFRACTED_HORIZON_RAD, false);
if (result_jd < 0.0)
PG_RETURN_NULL();
PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
}
/* ================================================================
* planet_next_rise_refracted(body_id, observer, timestamptz) -> timestamptz
*
* Uses -0.569 degree threshold (refraction only, point source).
* Planets are too small for semidiameter to matter — Jupiter at
* opposition is 24 arcseconds, <1 second of time error.
* ================================================================
*/
Datum
planet_next_rise_refracted(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 start_jd, stop_jd, result_jd;
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("planet_next_rise_refracted: 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")));
start_jd = timestamptz_to_jd(ts);
stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
result_jd = find_next_crossing(BTYPE_PLANET, body_id, obs,
start_jd, stop_jd,
REFRACTION_ONLY_HORIZON_RAD, true);
if (result_jd < 0.0)
PG_RETURN_NULL();
PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
}
/* ================================================================
* planet_next_set_refracted(body_id, observer, timestamptz) -> timestamptz
*
* Refracted planet set is later than geometric.
* ================================================================
*/
Datum
planet_next_set_refracted(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 start_jd, stop_jd, result_jd;
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("planet_next_set_refracted: 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")));
start_jd = timestamptz_to_jd(ts);
stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
result_jd = find_next_crossing(BTYPE_PLANET, body_id, obs,
start_jd, stop_jd,
REFRACTION_ONLY_HORIZON_RAD, false);
if (result_jd < 0.0)
PG_RETURN_NULL();
PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
}
/* ================================================================
* moon_next_rise_refracted(observer, timestamptz) -> timestamptz
*
* Uses -0.833 degree threshold (same as Sun: 0.569 deg refraction +
* 0.264 deg mean lunar semidiameter). Moon semidiameter varies
* 14.7'-16.7'; mean value error is ~1 arcmin → ~15 seconds in time.
* ================================================================
*/
Datum
moon_next_rise_refracted(PG_FUNCTION_ARGS)
{
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
int64 ts = PG_GETARG_INT64(1);
double start_jd, stop_jd, result_jd;
start_jd = timestamptz_to_jd(ts);
stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
result_jd = find_next_crossing(BTYPE_MOON, 0, obs,
start_jd, stop_jd,
SUN_MOON_REFRACTED_HORIZON_RAD, true);
if (result_jd < 0.0)
PG_RETURN_NULL();
PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
}
/* ================================================================
* moon_next_set_refracted(observer, timestamptz) -> timestamptz
*
* Refracted moonset is later than geometric.
* ================================================================
*/
Datum
moon_next_set_refracted(PG_FUNCTION_ARGS)
{
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
int64 ts = PG_GETARG_INT64(1);
double start_jd, stop_jd, result_jd;
start_jd = timestamptz_to_jd(ts);
stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
result_jd = find_next_crossing(BTYPE_MOON, 0, obs,
start_jd, stop_jd,
SUN_MOON_REFRACTED_HORIZON_RAD, false);
if (result_jd < 0.0)
PG_RETURN_NULL();
PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
}