pg_orrery/src/equatorial_funcs.c

/*
 * equatorial_funcs.c -- Equatorial coordinate type and observation functions
 *
 * Provides the equatorial PostgreSQL type (RA/Dec/distance) and functions
 * to compute equatorial coordinates for satellites, planets, Sun, Moon,
 * small bodies, and stars.
 *
 * The type stores 3 doubles (24 bytes): ra (radians), dec (radians),
 * distance (km).  RA is displayed in hours [0,24), dec in degrees [-90,90].
 *
 * Two satellite paths:
 *   - Topocentric (with observer): TEME -> ECEF -> subtract observer -> back
 *     to TEME -> RA/Dec from the range vector.
 *   - Geocentric (no observer): RA/Dec is the direction of the TEME position.
 *
 * Solar system path: VSOP87/ELP82B/Kepler heliocentric ecliptic J2000 ->
 * subtract Earth -> ecliptic-to-equatorial -> precess to date.
 */

#include "postgres.h"
#include "fmgr.h"
#include "utils/timestamp.h"
#include "utils/builtins.h"
#include "libpq/pqformat.h"
#include "types.h"
#include "astro_math.h"
#include "vsop87.h"
#include "elp82b.h"
#include <math.h>

/* Type I/O */
PG_FUNCTION_INFO_V1(equatorial_in);
PG_FUNCTION_INFO_V1(equatorial_out);
PG_FUNCTION_INFO_V1(equatorial_recv);
PG_FUNCTION_INFO_V1(equatorial_send);

/* Accessors */
PG_FUNCTION_INFO_V1(eq_ra);
PG_FUNCTION_INFO_V1(eq_dec);
PG_FUNCTION_INFO_V1(eq_distance);

/* Computation functions */
PG_FUNCTION_INFO_V1(eci_to_equatorial);
PG_FUNCTION_INFO_V1(eci_to_equatorial_geo);
PG_FUNCTION_INFO_V1(planet_equatorial);
PG_FUNCTION_INFO_V1(sun_equatorial);
PG_FUNCTION_INFO_V1(moon_equatorial);

/* Angular distance and cone search */
PG_FUNCTION_INFO_V1(eq_angular_distance);
PG_FUNCTION_INFO_V1(eq_within_cone);


/* ----------------------------------------------------------------
 * Static helper -- observer geodetic to ECEF.
 *
 * Duplicated from pass_funcs.c / coord_funcs.c because both files
 * define it as static.  Too small to warrant a shared module, and
 * keeping it static preserves the no-cross-TU-symbol convention.
 * ----------------------------------------------------------------
 */
static void
observer_to_ecef(const pg_observer *obs, double *ecef)
{
    double  sinlat = sin(obs->lat);
    double  coslat = cos(obs->lat);
    double  sinlon = sin(obs->lon);
    double  coslon = cos(obs->lon);
    double  N;      /* radius of curvature in the prime vertical */
    double  alt_km = obs->alt_m / 1000.0;

    N = WGS84_A / sqrt(1.0 - WGS84_E2 * sinlat * sinlat);

    ecef[0] = (N + alt_km) * coslat * coslon;
    ecef[1] = (N + alt_km) * coslat * sinlon;
    ecef[2] = (N * (1.0 - WGS84_E2) + alt_km) * sinlat;
}


/* ================================================================
 * Type I/O
 *
 * Text format: (ra_hours, dec_degrees, distance_km)
 *
 * RA displayed in hours [0,24), stored in radians [0, 2*pi).
 * Dec displayed in degrees [-90,90], stored in radians.
 * Distance in km.
 * ================================================================
 */

Datum
equatorial_in(PG_FUNCTION_ARGS)
{
    char            *str = PG_GETARG_CSTRING(0);
    pg_equatorial   *result;
    double           ra_hours, dec_deg, distance;
    int              nfields;

    result = (pg_equatorial *) palloc(sizeof(pg_equatorial));

    nfields = sscanf(str, " ( %lf , %lf , %lf )",
                     &ra_hours, &dec_deg, &distance);

    if (nfields != 3)
        ereport(ERROR,
                (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
                 errmsg("invalid input syntax for type equatorial: \"%s\"", str),
                 errhint("Expected (ra_hours,dec_degrees,distance_km).")));

    if (ra_hours < 0.0 || ra_hours >= 24.0)
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("right ascension out of range: %.6f", ra_hours),
                 errhint("RA must be in [0, 24) hours.")));

    if (dec_deg < -90.0 || dec_deg > 90.0)
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("declination out of range: %.6f", dec_deg),
                 errhint("Declination must be between -90 and +90 degrees.")));

    result->ra       = ra_hours * (M_PI / 12.0);
    result->dec      = dec_deg * DEG_TO_RAD;
    result->distance = distance;

    PG_RETURN_POINTER(result);
}

Datum
equatorial_out(PG_FUNCTION_ARGS)
{
    pg_equatorial *eq = (pg_equatorial *) PG_GETARG_POINTER(0);

    PG_RETURN_CSTRING(psprintf("(%.8f,%.8f,%.3f)",
                               eq->ra * (12.0 / M_PI),
                               eq->dec * RAD_TO_DEG,
                               eq->distance));
}

Datum
equatorial_recv(PG_FUNCTION_ARGS)
{
    StringInfo      buf = (StringInfo) PG_GETARG_POINTER(0);
    pg_equatorial  *result;

    result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
    result->ra       = pq_getmsgfloat8(buf);
    result->dec      = pq_getmsgfloat8(buf);
    result->distance = pq_getmsgfloat8(buf);

    PG_RETURN_POINTER(result);
}

Datum
equatorial_send(PG_FUNCTION_ARGS)
{
    pg_equatorial  *eq = (pg_equatorial *) PG_GETARG_POINTER(0);
    StringInfoData  buf;

    pq_begintypsend(&buf);
    pq_sendfloat8(&buf, eq->ra);
    pq_sendfloat8(&buf, eq->dec);
    pq_sendfloat8(&buf, eq->distance);

    PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
}


/* ================================================================
 * Accessor functions
 *
 * RA returns hours, Dec returns degrees (matching display convention).
 * ================================================================
 */

Datum
eq_ra(PG_FUNCTION_ARGS)
{
    pg_equatorial *eq = (pg_equatorial *) PG_GETARG_POINTER(0);
    PG_RETURN_FLOAT8(eq->ra * (12.0 / M_PI));
}

Datum
eq_dec(PG_FUNCTION_ARGS)
{
    pg_equatorial *eq = (pg_equatorial *) PG_GETARG_POINTER(0);
    PG_RETURN_FLOAT8(eq->dec * RAD_TO_DEG);
}

Datum
eq_distance(PG_FUNCTION_ARGS)
{
    pg_equatorial *eq = (pg_equatorial *) PG_GETARG_POINTER(0);
    PG_RETURN_FLOAT8(eq->distance);
}


/* ================================================================
 * eci_to_equatorial(eci_position, observer, timestamptz) -> equatorial
 *
 * Topocentric satellite RA/Dec.  Subtracts observer parallax via
 * ECEF round-trip, then extracts RA/Dec from the range vector.
 * For LEO, topocentric vs geocentric parallax is ~1 degree.
 * ================================================================
 */
Datum
eci_to_equatorial(PG_FUNCTION_ARGS)
{
    pg_eci         *eci = (pg_eci *) PG_GETARG_POINTER(0);
    pg_observer    *obs = (pg_observer *) PG_GETARG_POINTER(1);
    int64           ts  = PG_GETARG_INT64(2);
    double          jd;
    double          pos_teme[3];
    double          obs_ecef[3];
    double          gmst;
    pg_equatorial  *result;

    jd = timestamptz_to_jd(ts);
    gmst = gmst_from_jd(jd);

    pos_teme[0] = eci->x;
    pos_teme[1] = eci->y;
    pos_teme[2] = eci->z;

    observer_to_ecef(obs, obs_ecef);

    result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
    teme_to_equatorial_topo(pos_teme, gmst, obs_ecef, result);

    PG_RETURN_POINTER(result);
}


/* ================================================================
 * eci_to_equatorial_geo(eci_position, timestamptz) -> equatorial
 *
 * Geocentric satellite RA/Dec.  No observer subtraction -- the
 * position vector direction in TEME is the RA/Dec.
 * ================================================================
 */
Datum
eci_to_equatorial_geo(PG_FUNCTION_ARGS)
{
    pg_eci         *eci = (pg_eci *) PG_GETARG_POINTER(0);
    int64           ts  = PG_GETARG_INT64(1);
    double          pos_teme[3];
    pg_equatorial  *result;

    /* ts used to establish the equatorial frame; teme_to_equatorial_geo()
     * doesn't need jd because TEME is already ~equatorial of date. */
    (void) ts;

    pos_teme[0] = eci->x;
    pos_teme[1] = eci->y;
    pos_teme[2] = eci->z;

    result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
    teme_to_equatorial_geo(pos_teme, result);

    PG_RETURN_POINTER(result);
}


/* ================================================================
 * planet_equatorial(body_id int, timestamptz) -> equatorial
 *
 * Geocentric RA/Dec of a planet via VSOP87.
 * Body IDs: 1=Mercury, ..., 8=Neptune (not Sun, Earth, or Moon).
 * ================================================================
 */
Datum
planet_equatorial(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];
    int             vsop_body;
    pg_equatorial  *result;

    if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("planet_equatorial: 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_equatorial *) palloc(sizeof(pg_equatorial));
    geocentric_to_equatorial(geo_ecl, jd, result);

    PG_RETURN_POINTER(result);
}


/* ================================================================
 * sun_equatorial(timestamptz) -> equatorial
 *
 * Geocentric RA/Dec of the Sun.  The Sun's geocentric position is
 * the negation of Earth's heliocentric VSOP87 position.
 * ================================================================
 */
Datum
sun_equatorial(PG_FUNCTION_ARGS)
{
    int64           ts = PG_GETARG_INT64(0);
    double          jd;
    double          earth_xyz[6];
    double          geo_ecl[3];
    pg_equatorial  *result;

    jd = timestamptz_to_jd(ts);

    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_equatorial *) palloc(sizeof(pg_equatorial));
    geocentric_to_equatorial(geo_ecl, jd, result);

    PG_RETURN_POINTER(result);
}


/* ================================================================
 * moon_equatorial(timestamptz) -> equatorial
 *
 * Geocentric RA/Dec of the Moon via ELP2000-82B.
 * ELP returns geocentric ecliptic J2000 in AU -- no Earth subtraction.
 * ================================================================
 */
Datum
moon_equatorial(PG_FUNCTION_ARGS)
{
    int64           ts = PG_GETARG_INT64(0);
    double          jd;
    double          moon_ecl[3];
    pg_equatorial  *result;

    jd = timestamptz_to_jd(ts);

    GetElp82bCoor(jd, moon_ecl);

    result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
    geocentric_to_equatorial(moon_ecl, jd, result);

    PG_RETURN_POINTER(result);
}


/* ================================================================
 * eq_angular_distance(equatorial, equatorial) -> float8
 *
 * Angular separation in degrees between two equatorial positions.
 * Uses the Vincenty formula (numerically stable at all separations
 * including 0 and 180 degrees, unlike the haversine or dot product).
 *
 * Vincenty: atan2(num, den) where
 *   num = sqrt((cos d2 sin dRA)^2 + (cos d1 sin d2 - sin d1 cos d2 cos dRA)^2)
 *   den = sin d1 sin d2 + cos d1 cos d2 cos dRA
 * ================================================================
 */
Datum
eq_angular_distance(PG_FUNCTION_ARGS)
{
    pg_equatorial *a = (pg_equatorial *) PG_GETARG_POINTER(0);
    pg_equatorial *b = (pg_equatorial *) PG_GETARG_POINTER(1);
    double d_ra, cos_d_ra, sin_d_ra;
    double sin_d1, cos_d1, sin_d2, cos_d2;
    double num1, num2, num, den;

    d_ra     = b->ra - a->ra;
    cos_d_ra = cos(d_ra);
    sin_d_ra = sin(d_ra);

    sin_d1 = sin(a->dec);
    cos_d1 = cos(a->dec);
    sin_d2 = sin(b->dec);
    cos_d2 = cos(b->dec);

    num1 = cos_d2 * sin_d_ra;
    num2 = cos_d1 * sin_d2 - sin_d1 * cos_d2 * cos_d_ra;
    num  = sqrt(num1 * num1 + num2 * num2);
    den  = sin_d1 * sin_d2 + cos_d1 * cos_d2 * cos_d_ra;

    PG_RETURN_FLOAT8(atan2(num, den) * RAD_TO_DEG);
}


/* ================================================================
 * eq_within_cone(equatorial, equatorial, float8) -> bool
 *
 * Returns true if the first position is within `radius_deg` degrees
 * of the second position.  Cosine shortcut avoids the expensive
 * atan2 in the Vincenty formula for most reject cases.
 * ================================================================
 */
Datum
eq_within_cone(PG_FUNCTION_ARGS)
{
    pg_equatorial *a      = (pg_equatorial *) PG_GETARG_POINTER(0);
    pg_equatorial *b      = (pg_equatorial *) PG_GETARG_POINTER(1);
    double         radius = PG_GETARG_FLOAT8(2);
    double cos_r, d_ra, cos_sep;

    if (radius < 0.0)
        PG_RETURN_BOOL(false);

    cos_r = cos(radius * DEG_TO_RAD);

    d_ra = b->ra - a->ra;
    cos_sep = sin(a->dec) * sin(b->dec)
            + cos(a->dec) * cos(b->dec) * cos(d_ra);

    PG_RETURN_BOOL(cos_sep >= cos_r);
}