/*
 * star_funcs.c -- Star and fixed-position object observation
 *
 * Takes J2000 catalog coordinates (RA in hours, Dec in degrees),
 * applies IAU 1976 precession to date of observation, computes
 * local hour angle, and converts to topocentric azimuth/elevation.
 *
 * Range and range_rate are zero -- stars are effectively at infinity.
 *
 * star_observe / star_observe_safe: catalog J2000 coords only.
 * star_observe_pm / star_equatorial_pm: proper motion, parallax, RV.
 * star_equatorial: catalog J2000 to apparent equatorial of date.
 */

#include "postgres.h"
#include "fmgr.h"
#include "utils/timestamp.h"
#include "types.h"
#include "astro_math.h"
#include "vsop87.h"

PG_FUNCTION_INFO_V1(star_observe);
PG_FUNCTION_INFO_V1(star_observe_safe);
PG_FUNCTION_INFO_V1(star_observe_pm);
PG_FUNCTION_INFO_V1(star_equatorial_pm);
PG_FUNCTION_INFO_V1(star_equatorial);

/*
 * star_observe(ra_hours, dec_degrees, observer, timestamptz) -> topocentric
 *
 * Compute az/el of a fixed celestial object from an observer at a time.
 * Uses IAU 1976 precession (~1 arcsecond accuracy for centuries near J2000).
 */
Datum
star_observe(PG_FUNCTION_ARGS)
{
    double       ra_hours = PG_GETARG_FLOAT8(0);
    double       dec_deg  = PG_GETARG_FLOAT8(1);
    pg_observer *obs      = (pg_observer *) PG_GETARG_POINTER(2);
    int64        ts       = PG_GETARG_INT64(3);

    double          jd;
    double          ra_j2000, dec_j2000;
    double          ra_date, dec_date;
    double          gmst, lst, ha;
    double          az, el;
    pg_topocentric *result;

    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.")));

    jd = timestamptz_to_jd(ts);

    ra_j2000  = ra_hours * (M_PI / 12.0);
    dec_j2000 = dec_deg * DEG_TO_RAD;

    precess_j2000_to_date(jd, ra_j2000, dec_j2000, &ra_date, &dec_date);

    gmst = gmst_from_jd(jd);
    lst  = gmst + obs->lon;
    ha   = lst - ra_date;

    equatorial_to_horizontal(ha, dec_date, obs->lat, &az, &el);

    result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
    result->azimuth    = az;
    result->elevation  = el;
    result->range_km   = 0.0;
    result->range_rate = 0.0;

    PG_RETURN_POINTER(result);
}


/*
 * star_observe_safe -- returns NULL if inputs are out of range.
 * For batch queries over star catalogs.
 */
Datum
star_observe_safe(PG_FUNCTION_ARGS)
{
    double       ra_hours = PG_GETARG_FLOAT8(0);
    double       dec_deg  = PG_GETARG_FLOAT8(1);
    pg_observer *obs      = (pg_observer *) PG_GETARG_POINTER(2);
    int64        ts       = PG_GETARG_INT64(3);

    double          jd;
    double          ra_j2000, dec_j2000;
    double          ra_date, dec_date;
    double          gmst, lst, ha;
    double          az, el;
    pg_topocentric *result;

    if (ra_hours < 0.0 || ra_hours >= 24.0)
        PG_RETURN_NULL();
    if (dec_deg < -90.0 || dec_deg > 90.0)
        PG_RETURN_NULL();

    jd = timestamptz_to_jd(ts);

    ra_j2000  = ra_hours * (M_PI / 12.0);
    dec_j2000 = dec_deg * DEG_TO_RAD;

    precess_j2000_to_date(jd, ra_j2000, dec_j2000, &ra_date, &dec_date);

    gmst = gmst_from_jd(jd);
    lst  = gmst + obs->lon;
    ha   = lst - ra_date;

    equatorial_to_horizontal(ha, dec_date, obs->lat, &az, &el);

    result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
    result->azimuth    = az;
    result->elevation  = el;
    result->range_km   = 0.0;
    result->range_rate = 0.0;

    PG_RETURN_POINTER(result);
}


/*
 * star_observe_pm(ra_hours, dec_deg, pm_ra_masyr, pm_dec_masyr,
 *                 parallax_mas, rv_kms, observer, timestamptz) -> topocentric
 *
 * Full star observation with proper motion correction applied before
 * IAU 1976 precession.  pm_ra is mu_alpha*cos(delta) in mas/yr
 * (Hipparcos/Gaia convention).  Parallax and radial velocity accepted
 * for API symmetry with star_equatorial_pm but not yet applied to
 * the topocentric result (would require Earth's heliocentric position).
 */
Datum
star_observe_pm(PG_FUNCTION_ARGS)
{
    double       ra_hours     = PG_GETARG_FLOAT8(0);
    double       dec_deg      = PG_GETARG_FLOAT8(1);
    double       pm_ra_masyr  = PG_GETARG_FLOAT8(2);
    double       pm_dec_masyr = PG_GETARG_FLOAT8(3);
    double       parallax_mas = PG_GETARG_FLOAT8(4);
    double       rv_kms       = PG_GETARG_FLOAT8(5);
    pg_observer *obs          = (pg_observer *) PG_GETARG_POINTER(6);
    int64        ts           = PG_GETARG_INT64(7);

    double       jd, dt_years;
    double       ra_j2000, dec_j2000;
    double       cos_dec, ra_corrected, dec_corrected;
    double       ra_date, dec_date;
    double       gmst, lst, ha;
    double       az, el;
    pg_topocentric *result;

    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.")));

    jd = timestamptz_to_jd(ts);
    dt_years = (jd - J2000_JD) / 365.25;

    if (fabs(dt_years) > 200.0)
        ereport(NOTICE,
                (errmsg("proper motion extrapolation %.0f years from J2000 — accuracy degrades beyond ~200 years", dt_years)));

    ra_j2000  = ra_hours * (M_PI / 12.0);
    dec_j2000 = dec_deg * DEG_TO_RAD;

    /* Apply proper motion (linear in J2000 frame).
     * pm_ra is mu_alpha*cos(delta), so divide by cos(dec) to get dRA. */
    cos_dec = cos(dec_j2000);
    if (fabs(cos_dec) < cos(89.99 * DEG_TO_RAD))
        cos_dec = (cos_dec >= 0.0 ? 1.0 : -1.0) * cos(89.99 * DEG_TO_RAD);

    ra_corrected  = ra_j2000  + (pm_ra_masyr / 3.6e6) * DEG_TO_RAD / cos_dec * dt_years;
    dec_corrected = dec_j2000 + (pm_dec_masyr / 3.6e6) * DEG_TO_RAD * dt_years;

    if (dec_corrected > M_PI / 2.0)
        dec_corrected = M_PI / 2.0;
    if (dec_corrected < -M_PI / 2.0)
        dec_corrected = -M_PI / 2.0;

    ra_corrected = fmod(ra_corrected, 2.0 * M_PI);
    if (ra_corrected < 0.0)
        ra_corrected += 2.0 * M_PI;

    /* Annual parallax: displace star position by Earth's heliocentric
     * position projected onto the star direction.
     * Green (1985), Eq. 11.3: d_RA and d_Dec from parallax. */
    if (parallax_mas > 0.0)
    {
        double earth_xyz[6], earth_equ[3];
        double p_rad = (parallax_mas / 1000.0) * ARCSEC_TO_RAD;
        double sin_ra_c = sin(ra_corrected);
        double cos_ra_c = cos(ra_corrected);
        double sin_dec_c = sin(dec_corrected);
        double cos_dec_c = cos(dec_corrected);

        GetVsop87Coor(jd, 2, earth_xyz);
        ecliptic_to_equatorial(earth_xyz, earth_equ);

        /* Parallax displacement (Green 1985 Eq. 11.3) */
        if (fabs(cos_dec_c) > 1e-12)
        {
            ra_corrected  += p_rad * (earth_equ[0] * sin_ra_c
                                    - earth_equ[1] * cos_ra_c) / cos_dec_c;
            dec_corrected += p_rad * (earth_equ[0] * cos_ra_c * sin_dec_c
                                    + earth_equ[1] * sin_ra_c * sin_dec_c
                                    - earth_equ[2] * cos_dec_c);
        }

        ra_corrected = fmod(ra_corrected, 2.0 * M_PI);
        if (ra_corrected < 0.0)
            ra_corrected += 2.0 * M_PI;
    }
    (void) rv_kms;

    precess_j2000_to_date(jd, ra_corrected, dec_corrected, &ra_date, &dec_date);

    gmst = gmst_from_jd(jd);
    lst  = gmst + obs->lon;
    ha   = lst - ra_date;

    equatorial_to_horizontal(ha, dec_date, obs->lat, &az, &el);

    result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
    result->azimuth    = az;
    result->elevation  = el;
    result->range_km   = 0.0;
    result->range_rate = 0.0;

    PG_RETURN_POINTER(result);
}


/*
 * star_equatorial_pm(ra_hours, dec_deg, pm_ra_masyr, pm_dec_masyr,
 *                    parallax_mas, rv_kms, timestamptz) -> equatorial
 *
 * Proper-motion-corrected apparent equatorial coordinates of date.
 * No observer needed (geocentric).  If parallax > 0, distance is
 * derived as 1000/parallax_mas parsecs converted to km.
 */
Datum
star_equatorial_pm(PG_FUNCTION_ARGS)
{
    double        ra_hours     = PG_GETARG_FLOAT8(0);
    double        dec_deg      = PG_GETARG_FLOAT8(1);
    double        pm_ra_masyr  = PG_GETARG_FLOAT8(2);
    double        pm_dec_masyr = PG_GETARG_FLOAT8(3);
    double        parallax_mas = PG_GETARG_FLOAT8(4);
    double        rv_kms       = PG_GETARG_FLOAT8(5);
    int64         ts           = PG_GETARG_INT64(6);

    double        jd, dt_years;
    double        ra_j2000, dec_j2000;
    double        cos_dec, ra_corrected, dec_corrected;
    double        ra_date, dec_date;
    pg_equatorial *result;

    (void) rv_kms;

    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.")));

    jd = timestamptz_to_jd(ts);
    dt_years = (jd - J2000_JD) / 365.25;

    if (fabs(dt_years) > 200.0)
        ereport(NOTICE,
                (errmsg("proper motion extrapolation %.0f years from J2000 — accuracy degrades beyond ~200 years", dt_years)));

    ra_j2000  = ra_hours * (M_PI / 12.0);
    dec_j2000 = dec_deg * DEG_TO_RAD;

    cos_dec = cos(dec_j2000);
    if (fabs(cos_dec) < cos(89.99 * DEG_TO_RAD))
        cos_dec = (cos_dec >= 0.0 ? 1.0 : -1.0) * cos(89.99 * DEG_TO_RAD);

    ra_corrected  = ra_j2000  + (pm_ra_masyr / 3.6e6) * DEG_TO_RAD / cos_dec * dt_years;
    dec_corrected = dec_j2000 + (pm_dec_masyr / 3.6e6) * DEG_TO_RAD * dt_years;

    if (dec_corrected > M_PI / 2.0)
        dec_corrected = M_PI / 2.0;
    if (dec_corrected < -M_PI / 2.0)
        dec_corrected = -M_PI / 2.0;

    ra_corrected = fmod(ra_corrected, 2.0 * M_PI);
    if (ra_corrected < 0.0)
        ra_corrected += 2.0 * M_PI;

    /* Annual parallax displacement (Green 1985 Eq. 11.3) */
    if (parallax_mas > 0.0)
    {
        double earth_xyz[6], earth_equ[3];
        double p_rad = (parallax_mas / 1000.0) * ARCSEC_TO_RAD;
        double sin_ra_c = sin(ra_corrected);
        double cos_ra_c = cos(ra_corrected);
        double sin_dec_c = sin(dec_corrected);
        double cos_dec_c = cos(dec_corrected);

        GetVsop87Coor(jd, 2, earth_xyz);
        ecliptic_to_equatorial(earth_xyz, earth_equ);

        if (fabs(cos_dec_c) > 1e-12)
        {
            ra_corrected  += p_rad * (earth_equ[0] * sin_ra_c
                                    - earth_equ[1] * cos_ra_c) / cos_dec_c;
            dec_corrected += p_rad * (earth_equ[0] * cos_ra_c * sin_dec_c
                                    + earth_equ[1] * sin_ra_c * sin_dec_c
                                    - earth_equ[2] * cos_dec_c);
        }

        ra_corrected = fmod(ra_corrected, 2.0 * M_PI);
        if (ra_corrected < 0.0)
            ra_corrected += 2.0 * M_PI;
    }

    precess_j2000_to_date(jd, ra_corrected, dec_corrected, &ra_date, &dec_date);

    result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
    result->ra  = ra_date;
    result->dec = dec_date;

    /* Distance from parallax: d_pc = 1000/parallax_mas, d_AU = d_pc * 206265 */
    if (parallax_mas > 0.0)
        result->distance = (1000.0 / parallax_mas) * 206265.0 * AU_KM;
    else
        result->distance = 0.0;

    PG_RETURN_POINTER(result);
}


/*
 * star_equatorial(ra_hours, dec_deg, timestamptz) -> equatorial
 *
 * Precess J2000 catalog coordinates to apparent equatorial of date.
 * Distance is zero (no parallax information).
 */
Datum
star_equatorial(PG_FUNCTION_ARGS)
{
    double        ra_hours = PG_GETARG_FLOAT8(0);
    double        dec_deg  = PG_GETARG_FLOAT8(1);
    int64         ts       = PG_GETARG_INT64(2);

    double        jd, ra_j2000, dec_j2000, ra_date, dec_date;
    pg_equatorial *result;

    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.")));

    jd = timestamptz_to_jd(ts);

    ra_j2000  = ra_hours * (M_PI / 12.0);
    dec_j2000 = dec_deg * DEG_TO_RAD;

    precess_j2000_to_date(jd, ra_j2000, dec_j2000, &ra_date, &dec_date);

    result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
    result->ra       = ra_date;
    result->dec      = dec_date;
    result->distance = 0.0;

    PG_RETURN_POINTER(result);
}