/* * 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 /* 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); /* ---------------------------------------------------------------- * 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); }