pg_orrery/src/radio_funcs.c

/*
 * radio_funcs.c -- Jupiter decametric radio burst prediction
 *
 * Jupiter emits intense radio bursts at 4-39.5 MHz, modulated by:
 *   1. Io's orbital phase around Jupiter (Io-related storms)
 *   2. Jupiter's Central Meridian Longitude (CML, System III)
 *
 * Radio JOVE operators use CML-vs-Io-phase lookup charts to predict
 * burst probability. This module computes both parameters from SQL,
 * enabling batch prediction queries over time series.
 *
 * References:
 *   Carr, Desch & Alexander (1983), "Phenomenology of Magnetospheric
 *   Radio Emissions", in Physics of the Jovian Magnetosphere.
 *   Radio JOVE project: https://radiojove.gsfc.nasa.gov/
 */

#include "postgres.h"
#include "fmgr.h"
#include "funcapi.h"
#include "utils/timestamp.h"
#include "types.h"
#include "astro_math.h"
#include "vsop87.h"
#include "l12.h"
#include <math.h>

PG_FUNCTION_INFO_V1(io_phase_angle);
PG_FUNCTION_INFO_V1(jupiter_cml);
PG_FUNCTION_INFO_V1(jupiter_burst_probability);


/* ================================================================
 * io_phase_angle(timestamptz) -> float8
 *
 * Returns Io's orbital phase angle in degrees [0, 360).
 * Measured from superior geocentric conjunction:
 *   0 = Io behind Jupiter (as seen from Earth)
 *   90 = Io east of Jupiter (western elongation)
 *   180 = Io in front of Jupiter (inferior conjunction)
 *   270 = Io west of Jupiter (eastern elongation)
 *
 * This is the standard convention used by Radio JOVE charts.
 * ================================================================
 */
Datum
io_phase_angle(PG_FUNCTION_ARGS)
{
    int64   ts = PG_GETARG_INT64(0);
    double  jd;
    double  io_xyz[3];
    double  jup_xyz[6], earth_xyz[6];
    double  jup_geo[3];
    double  io_rel_geo[3];
    double  phase;
    double  proj_x, proj_y;
    double  jup_dist;
    double  east_x, east_y, east_z;
    double  north_x, north_y, north_z;
    double  io_east, io_north;

    jd = timestamptz_to_jd(ts);

    /* Io position relative to Jupiter (VSOP87 ecliptic J2000, AU) */
    GetL12Coor(jd, L12_IO, io_xyz, NULL);

    /* Jupiter and Earth heliocentric positions */
    GetVsop87Coor(jd, 4, jup_xyz);   /* VSOP87 body 4 = Jupiter */
    GetVsop87Coor(jd, 2, earth_xyz); /* VSOP87 body 2 = Earth */

    /* Jupiter geocentric direction */
    jup_geo[0] = jup_xyz[0] - earth_xyz[0];
    jup_geo[1] = jup_xyz[1] - earth_xyz[1];
    jup_geo[2] = jup_xyz[2] - earth_xyz[2];

    jup_dist = sqrt(jup_geo[0] * jup_geo[0] +
                    jup_geo[1] * jup_geo[1] +
                    jup_geo[2] * jup_geo[2]);

    /* Normalize Jupiter geocentric direction */
    jup_geo[0] /= jup_dist;
    jup_geo[1] /= jup_dist;
    jup_geo[2] /= jup_dist;

    /*
     * Project Io's position onto the plane perpendicular to
     * the Earth-Jupiter line of sight. We need two orthogonal
     * axes in this plane: "east" (position angle) and "north"
     * (toward north ecliptic pole).
     *
     * North ecliptic pole = (0, 0, 1) in ecliptic J2000.
     * East = North x LineOfSight (cross product)
     */
    north_x = -(jup_geo[2] * jup_geo[0]);  /* north - (north.los)*los */
    north_y = -(jup_geo[2] * jup_geo[1]);
    north_z = 1.0 - jup_geo[2] * jup_geo[2];

    {
        double nmag = sqrt(north_x * north_x + north_y * north_y + north_z * north_z);
        if (nmag > 1e-10) {
            north_x /= nmag;
            north_y /= nmag;
            north_z /= nmag;
        }
    }

    /* East = North x LOS */
    east_x = north_y * jup_geo[2] - north_z * jup_geo[1];
    east_y = north_z * jup_geo[0] - north_x * jup_geo[2];
    east_z = north_x * jup_geo[1] - north_y * jup_geo[0];

    /* Project Io relative position onto east/north */
    io_east  = io_xyz[0] * east_x  + io_xyz[1] * east_y  + io_xyz[2] * east_z;
    io_north = io_xyz[0] * north_x + io_xyz[1] * north_y + io_xyz[2] * north_z;

    /* Along LOS component (positive = behind Jupiter from Earth) */
    proj_x = io_xyz[0] * jup_geo[0] + io_xyz[1] * jup_geo[1] + io_xyz[2] * jup_geo[2];
    proj_y = io_east;

    /*
     * Phase angle convention (Radio JOVE):
     *   0 = superior conjunction (Io behind Jupiter)
     *   Measured counter-clockwise as seen from north ecliptic pole
     */
    phase = atan2(-proj_y, proj_x) * RAD_TO_DEG;
    if (phase < 0.0)
        phase += 360.0;

    PG_RETURN_FLOAT8(phase);
}


/* ================================================================
 * jupiter_cml(observer, timestamptz) -> float8
 *
 * Jupiter's Central Meridian Longitude in System III (1965.0).
 * Returns degrees [0, 360).
 *
 * System III is defined by Jupiter's magnetic field rotation,
 * with period 9h 55m 29.711s (IAU).
 *
 * The CML is the System III longitude of the sub-observer point.
 * For Earth-based observation, this is approximately the
 * sub-Earth longitude.
 *
 * Reference: IAU (1976) System III longitude definition.
 * ================================================================
 */
Datum
jupiter_cml(PG_FUNCTION_ARGS)
{
    pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
    int64        ts  = PG_GETARG_INT64(1);
    double       jd;
    double       jup_xyz[6], earth_xyz[6];
    double       dx, dy, dz, geo_dist;
    double       light_time_days;
    double       jd_lt;
    double       d;
    double       cml;

    (void)obs;  /* observer not needed for CML (geocentric is sufficient) */

    jd = timestamptz_to_jd(ts);

    /* Jupiter and Earth positions */
    GetVsop87Coor(jd, 4, jup_xyz);
    GetVsop87Coor(jd, 2, earth_xyz);

    /* Geocentric distance for light-time correction */
    dx = jup_xyz[0] - earth_xyz[0];
    dy = jup_xyz[1] - earth_xyz[1];
    dz = jup_xyz[2] - earth_xyz[2];
    geo_dist = sqrt(dx * dx + dy * dy + dz * dz);

    /* Light-time correction (speed of light = 173.1446 AU/day) */
    light_time_days = geo_dist / 173.1446327;
    jd_lt = jd - light_time_days;

    /*
     * System III longitude of central meridian.
     *
     * Reference epoch: J1965.0 = JD 2438761.5
     * Initial longitude: 305.38 deg (IAU 1976 convention)
     * Rotation rate: 870.536 deg/day (System III, 9h 55m 29.711s period)
     *
     * The sub-Earth longitude decreases with time (Jupiter rotates
     * prograde), so CML = initial + rate * elapsed_days.
     */
    d = jd_lt - 2438761.5;  /* Days since J1965.0 */
    cml = 305.38 + 870.536 * d;

    /* Normalize to [0, 360) */
    cml = fmod(cml, 360.0);
    if (cml < 0.0)
        cml += 360.0;

    PG_RETURN_FLOAT8(cml);
}


/* ================================================================
 * jupiter_burst_probability(io_phase float8, cml float8) -> float8
 *
 * Estimated probability of Jupiter decametric radio burst
 * based on Io phase angle and CML (System III).
 *
 * Returns a value 0.0-1.0 representing relative likelihood.
 *
 * Based on the empirical source regions from Carr et al. (1983):
 *   Source A: CML 200-260, Io phase 195-265 (Io-A)
 *   Source B: CML 90-200,  Io phase 70-110  (Io-B)
 *   Source C: CML 290-10,  Io phase 215-260 (Io-C)
 *   Source D: CML 0-50,    no Io dependence  (non-Io-D)
 *
 * This is a simplified model. Real-world burst probability
 * depends on additional factors (frequency, solar activity,
 * antenna orientation).
 * ================================================================
 */
Datum
jupiter_burst_probability(PG_FUNCTION_ARGS)
{
    double io_phase = PG_GETARG_FLOAT8(0);
    double cml      = PG_GETARG_FLOAT8(1);
    double p;

    /* Normalize inputs */
    io_phase = fmod(io_phase, 360.0);
    if (io_phase < 0.0) io_phase += 360.0;
    cml = fmod(cml, 360.0);
    if (cml < 0.0) cml += 360.0;

    p = 0.0;

    /* Source A (Io-A): strongest, most predictable */
    if (cml >= 200.0 && cml <= 260.0 && io_phase >= 195.0 && io_phase <= 265.0)
        p = 0.8;

    /* Source B (Io-B): second strongest */
    if (cml >= 90.0 && cml <= 200.0 && io_phase >= 70.0 && io_phase <= 110.0)
    {
        double pb = 0.6;
        if (pb > p) p = pb;
    }

    /* Source C (Io-C): wraps around 360/0 boundary */
    if ((cml >= 290.0 || cml <= 10.0) && io_phase >= 215.0 && io_phase <= 260.0)
    {
        double pc = 0.5;
        if (pc > p) p = pc;
    }

    /* Non-Io-D: CML-dependent only, weaker */
    if (cml >= 0.0 && cml <= 50.0)
    {
        double pd = 0.15;
        if (pd > p) p = pd;
    }

    PG_RETURN_FLOAT8(p);
}