/*
 * test_lagrange.c -- Standalone unit test for the Lagrange solver
 *
 * Verifies quintic solutions, L4/L5 geometry, Hill radius,
 * zone radius, and co-rotating to physical frame transform.
 *
 * No PostgreSQL dependency.
 *
 * Build:  cc -Wall -Werror -Isrc -o test/test_lagrange \
 *             test/test_lagrange.c -lm
 * Run:    ./test/test_lagrange
 */

#include "lagrange.h"

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* ── Test harness ───────────────────────────────────────── */

static int n_run, n_pass;

#define RUN(cond, msg) do { \
    n_run++; \
    if (!(cond)) \
        fprintf(stderr, "FAIL: %s  [line %d]\n", (msg), __LINE__); \
    else { n_pass++; fprintf(stderr, "  ok: %s\n", (msg)); } \
} while (0)

#define CLOSE(a, b, tol, msg) do { \
    n_run++; \
    double _a = (a), _b = (b); \
    if (fabs(_a - _b) > (tol)) \
        fprintf(stderr, "FAIL: %s: %.15g vs %.15g (diff %.3e)  [line %d]\n", \
                (msg), _a, _b, fabs(_a - _b), __LINE__); \
    else { n_pass++; fprintf(stderr, "  ok: %s\n", (msg)); } \
} while (0)

/* ── Tests ─────────────────────────────────────────────── */

/*
 * Verify equilibrium: at a Lagrange point, the net force in the
 * co-rotating frame should vanish. We check the effective potential
 * gradient by evaluating the quintic polynomial.
 */
static void
test_equilibrium_check(double mu, int point_id, const char *label)
{
    double x, y;
    int    rc;
    char   buf[128];

    rc = lagrange_corotating(mu, point_id, &x, &y);
    snprintf(buf, sizeof(buf), "%s: convergence", label);
    RUN(rc == 0, buf);

    if (rc != 0)
        return;

    if (point_id <= LAGRANGE_L3)
    {
        /*
         * For collinear points, verify equilibrium directly.
         * At equilibrium on the x-axis:
         *   x - (1-mu)*(x+mu)/|x+mu|^3 - mu*(x-1+mu)/|x-1+mu|^3 = 0
         */
        double dx1 = x + mu;       /* distance from primary (at -mu) */
        double dx2 = x - 1.0 + mu; /* distance from secondary (at 1-mu) */
        double r1  = fabs(dx1);
        double r2  = fabs(dx2);
        double residual;

        residual = x - (1.0 - mu) * dx1 / (r1 * r1 * r1)
                     - mu * dx2 / (r2 * r2 * r2);

        snprintf(buf, sizeof(buf), "%s: equilibrium residual", label);
        CLOSE(residual, 0.0, 1e-12, buf);
    }
    else
    {
        /* L4/L5: equidistant from both primaries at unit distance */
        double r1 = sqrt((x + mu) * (x + mu) + y * y);
        double r2 = sqrt((x - 1.0 + mu) * (x - 1.0 + mu) + y * y);

        snprintf(buf, sizeof(buf), "%s: distance to primary", label);
        CLOSE(r1, 1.0, 1e-14, buf);

        snprintf(buf, sizeof(buf), "%s: distance to secondary", label);
        CLOSE(r2, 1.0, 1e-14, buf);
    }
}


static void
test_sun_earth(void)
{
    double mu = mu_from_ratio(SUN_EARTH_RATIO);
    double x, y;
    int    rc;

    fprintf(stderr, "\n── Sun-Earth system (mu = %.6e) ──\n", mu);

    /* L1: between Sun and Earth, ~0.01 AU from Earth */
    rc = lagrange_corotating(mu, LAGRANGE_L1, &x, &y);
    RUN(rc == 0, "L1 converges");
    /* L1 should be between barycenter and secondary */
    RUN(x > -mu && x < 1.0 - mu, "L1 between primaries");

    /* Distance from secondary (Earth at 1-mu) */
    {
        double d_from_earth = (1.0 - mu) - x;
        CLOSE(d_from_earth, 0.01, 0.002, "L1 ~0.01 AU from Earth");
    }

    /* L2: beyond Earth, also ~0.01 AU */
    rc = lagrange_corotating(mu, LAGRANGE_L2, &x, &y);
    RUN(rc == 0, "L2 converges");
    {
        double d_from_earth = x - (1.0 - mu);
        CLOSE(d_from_earth, 0.01, 0.002, "L2 ~0.01 AU from Earth");
    }

    /* L3: opposite side from Earth */
    rc = lagrange_corotating(mu, LAGRANGE_L3, &x, &y);
    RUN(rc == 0, "L3 converges");
    RUN(x < -mu, "L3 beyond primary (opposite side)");

    test_equilibrium_check(mu, LAGRANGE_L1, "Sun-Earth L1");
    test_equilibrium_check(mu, LAGRANGE_L2, "Sun-Earth L2");
    test_equilibrium_check(mu, LAGRANGE_L3, "Sun-Earth L3");
    test_equilibrium_check(mu, LAGRANGE_L4, "Sun-Earth L4");
    test_equilibrium_check(mu, LAGRANGE_L5, "Sun-Earth L5");
}


static void
test_sun_jupiter(void)
{
    double mu = mu_from_ratio(SUN_JUPITER_RATIO);
    double x, y;
    int    rc;

    fprintf(stderr, "\n── Sun-Jupiter system (mu = %.6e) ──\n", mu);

    /* L4/L5: should be at 60 degrees from Jupiter */
    rc = lagrange_corotating(mu, LAGRANGE_L4, &x, &y);
    RUN(rc == 0, "L4 converges");
    {
        /* Angle from secondary: atan2(y, x - (1-mu)) */
        double angle = atan2(y, x - (1.0 - mu));
        double angle_deg = angle * 180.0 / M_PI;
        /* L4 leads secondary by ~60 degrees (but angle from barycenter) */
        /* Actually check equilateral property */
        double d_prim = sqrt((x + mu) * (x + mu) + y * y);
        double d_sec  = sqrt((x - 1.0 + mu) * (x - 1.0 + mu) + y * y);
        CLOSE(d_prim, 1.0, 1e-14, "L4 unit distance from primary");
        CLOSE(d_sec,  1.0, 1e-14, "L4 unit distance from secondary");
        RUN(y > 0.0, "L4 above x-axis (leading)");
        (void)angle_deg;  /* used implicitly via assertions */
    }

    test_equilibrium_check(mu, LAGRANGE_L1, "Sun-Jupiter L1");
    test_equilibrium_check(mu, LAGRANGE_L2, "Sun-Jupiter L2");
    test_equilibrium_check(mu, LAGRANGE_L3, "Sun-Jupiter L3");
    test_equilibrium_check(mu, LAGRANGE_L4, "Sun-Jupiter L4");
    test_equilibrium_check(mu, LAGRANGE_L5, "Sun-Jupiter L5");
}


static void
test_earth_moon(void)
{
    double mu = mu_from_ratio(EARTH_MOON_EMRAT);

    fprintf(stderr, "\n── Earth-Moon system (mu = %.6e) ──\n", mu);

    test_equilibrium_check(mu, LAGRANGE_L1, "Earth-Moon L1");
    test_equilibrium_check(mu, LAGRANGE_L2, "Earth-Moon L2");
    test_equilibrium_check(mu, LAGRANGE_L3, "Earth-Moon L3");
    test_equilibrium_check(mu, LAGRANGE_L4, "Earth-Moon L4");
    test_equilibrium_check(mu, LAGRANGE_L5, "Earth-Moon L5");

    /* Earth-Moon L1 should be ~326,000 km from Earth (~84.7% of separation) */
    {
        double x, y;
        int    rc;
        double earth_moon_km = 384400.0;  /* mean distance */

        rc = lagrange_corotating(mu, LAGRANGE_L1, &x, &y);
        RUN(rc == 0, "Earth-Moon L1 converges");

        /* In co-rotating frame, Earth is at -mu, Moon at 1-mu.
         * L1 is between them. Distance from Earth = x + mu. */
        {
            double frac = (x + mu);  /* fraction of separation from Earth */
            double km_from_earth = frac * earth_moon_km;
            CLOSE(km_from_earth, 326000.0, 5000.0,
                  "E-M L1 ~326,000 km from Earth");
        }
    }
}


static void
test_l4_l5_symmetry(void)
{
    double mu = mu_from_ratio(SUN_JUPITER_RATIO);
    double x4, y4, x5, y5;
    int    rc;

    fprintf(stderr, "\n── L4/L5 symmetry ──\n");

    rc = lagrange_corotating(mu, LAGRANGE_L4, &x4, &y4);
    RUN(rc == 0, "L4 converges");
    rc = lagrange_corotating(mu, LAGRANGE_L5, &x5, &y5);
    RUN(rc == 0, "L5 converges");

    CLOSE(x4, x5, 1e-15, "L4 and L5 same x-coordinate");
    CLOSE(y4, -y5, 1e-15, "L4 and L5 mirror in y");
}


static void
test_l1_l2_ordering(void)
{
    double mu = mu_from_ratio(SUN_EARTH_RATIO);
    double x1, y1, x2, y2, x3, y3;
    int    rc;

    fprintf(stderr, "\n── L1/L2/L3 ordering ──\n");

    rc = lagrange_corotating(mu, LAGRANGE_L1, &x1, &y1);
    RUN(rc == 0, "L1 converges");
    rc = lagrange_corotating(mu, LAGRANGE_L2, &x2, &y2);
    RUN(rc == 0, "L2 converges");
    rc = lagrange_corotating(mu, LAGRANGE_L3, &x3, &y3);
    RUN(rc == 0, "L3 converges");

    /* Ordering: L3 < primary < L1 < secondary < L2 */
    RUN(x3 < -mu, "L3 < primary");
    RUN(x1 > -mu && x1 < 1.0 - mu, "L1 between primaries");
    RUN(x2 > 1.0 - mu, "L2 beyond secondary");
}


static void
test_hill_radius(void)
{
    double mu_jup, mu_earth;
    double hill_jup, hill_earth;

    fprintf(stderr, "\n── Hill radius ──\n");

    mu_jup   = mu_from_ratio(SUN_JUPITER_RATIO);
    mu_earth = mu_from_ratio(SUN_EARTH_RATIO);

    /* Jupiter at ~5.2 AU */
    hill_jup = lagrange_hill_radius(5.2, mu_jup);
    CLOSE(hill_jup, 0.355, 0.02, "Jupiter Hill radius ~0.35 AU");

    /* Earth at ~1.0 AU */
    hill_earth = lagrange_hill_radius(1.0, mu_earth);
    CLOSE(hill_earth, 0.01, 0.002, "Earth Hill radius ~0.01 AU");
}


static void
test_zone_radius(void)
{
    double mu = mu_from_ratio(SUN_JUPITER_RATIO);
    double zr;

    fprintf(stderr, "\n── Zone radius ──\n");

    zr = lagrange_zone_radius(5.2, mu, LAGRANGE_L1);
    RUN(zr > 0.0, "L1 zone radius positive");

    zr = lagrange_zone_radius(5.2, mu, LAGRANGE_L4);
    RUN(zr > 0.0, "L4 zone radius positive");

    zr = lagrange_zone_radius(5.2, mu, 99);
    RUN(zr < 0.0, "invalid point_id returns -1");
}


static void
test_physical_transform(void)
{
    double primary[3]   = {0.0, 0.0, 0.0};            /* Sun at origin */
    double secondary[3] = {1.0, 0.0, 0.0};            /* "planet" at 1 AU on x-axis */
    double sec_vel[3]   = {0.0, 0.01720209895, 0.0};  /* ~Gauss constant, circular */
    double mu = 0.001;                                 /* ~Jupiter-like */
    double result[3];
    int    rc;

    fprintf(stderr, "\n── Physical frame transform ──\n");

    /* L1: should be between Sun and planet, on x-axis */
    rc = lagrange_position(primary, secondary, sec_vel, mu, LAGRANGE_L1, result);
    RUN(rc == 0, "L1 transform succeeds");
    RUN(result[0] > 0.0 && result[0] < 1.0, "L1 between Sun and planet on x-axis");
    CLOSE(result[1], 0.0, 1e-10, "L1 y-component ~0");
    CLOSE(result[2], 0.0, 1e-10, "L1 z-component ~0");

    /* L2: beyond planet */
    rc = lagrange_position(primary, secondary, sec_vel, mu, LAGRANGE_L2, result);
    RUN(rc == 0, "L2 transform succeeds");
    RUN(result[0] > 1.0, "L2 beyond planet");
    CLOSE(result[1], 0.0, 1e-10, "L2 y-component ~0");

    /* L4: 60 degrees ahead, above x-axis in ecliptic plane */
    rc = lagrange_position(primary, secondary, sec_vel, mu, LAGRANGE_L4, result);
    RUN(rc == 0, "L4 transform succeeds");
    {
        double dist = sqrt(result[0]*result[0] + result[1]*result[1] + result[2]*result[2]);
        /* L4 should be ~1 AU from Sun (equilateral triangle) */
        CLOSE(dist, 1.0, 0.01, "L4 ~1 AU from Sun");
        RUN(result[1] > 0.0, "L4 positive y (leading)");
    }

    /* L5: symmetric with L4 */
    rc = lagrange_position(primary, secondary, sec_vel, mu, LAGRANGE_L5, result);
    RUN(rc == 0, "L5 transform succeeds");
    RUN(result[1] < 0.0, "L5 negative y (trailing)");
}


static void
test_extreme_mass_ratios(void)
{
    double x, y;
    int    rc;

    fprintf(stderr, "\n── Extreme mass ratios ──\n");

    /* Very small mu (like Mercury around the Sun) */
    {
        double mu = mu_from_ratio(SUN_MERCURY_RATIO);  /* ~1.66e-7 */
        rc = lagrange_corotating(mu, LAGRANGE_L1, &x, &y);
        RUN(rc == 0, "tiny mu L1 converges");
        test_equilibrium_check(mu, LAGRANGE_L1, "Mercury L1");
    }

    /* Moderately large mu */
    {
        double mu = 0.1;
        rc = lagrange_corotating(mu, LAGRANGE_L1, &x, &y);
        RUN(rc == 0, "mu=0.1 L1 converges");
        test_equilibrium_check(mu, LAGRANGE_L1, "mu=0.1 L1");
        test_equilibrium_check(mu, LAGRANGE_L2, "mu=0.1 L2");
        test_equilibrium_check(mu, LAGRANGE_L3, "mu=0.1 L3");
    }

    /* Equal mass (mu = 0.5, maximum) */
    {
        double mu = 0.5;
        rc = lagrange_corotating(mu, LAGRANGE_L1, &x, &y);
        RUN(rc == 0, "mu=0.5 L1 converges");
        test_equilibrium_check(mu, LAGRANGE_L1, "mu=0.5 L1");
        test_equilibrium_check(mu, LAGRANGE_L2, "mu=0.5 L2");
        test_equilibrium_check(mu, LAGRANGE_L3, "mu=0.5 L3");
        /* L4/L5 at (0, +-sqrt(3)/2) for equal mass */
        rc = lagrange_corotating(mu, LAGRANGE_L4, &x, &y);
        RUN(rc == 0, "mu=0.5 L4 converges");
        CLOSE(x, 0.0, 1e-15, "mu=0.5 L4 x=0");
        CLOSE(y, sqrt(3.0)/2.0, 1e-15, "mu=0.5 L4 y=sqrt(3)/2");
    }
}


static void
test_error_cases(void)
{
    double x, y;
    int    rc;

    fprintf(stderr, "\n── Error cases ──\n");

    rc = lagrange_corotating(0.0, LAGRANGE_L1, &x, &y);
    RUN(rc != 0, "mu=0 rejected");

    rc = lagrange_corotating(-0.1, LAGRANGE_L1, &x, &y);
    RUN(rc != 0, "negative mu rejected");

    rc = lagrange_corotating(0.6, LAGRANGE_L1, &x, &y);
    RUN(rc != 0, "mu>0.5 rejected");

    rc = lagrange_corotating(0.01, 0, &x, &y);
    RUN(rc != 0, "point_id=0 rejected");

    rc = lagrange_corotating(0.01, 6, &x, &y);
    RUN(rc != 0, "point_id=6 rejected");

    /* Mass ratio lookups */
    RUN(sun_planet_ratio(1) > 0.0, "Mercury ratio valid");
    RUN(sun_planet_ratio(8) > 0.0, "Neptune ratio valid");
    RUN(sun_planet_ratio(0) < 0.0, "Sun ratio invalid");
    RUN(sun_planet_ratio(9) < 0.0, "body 9 invalid");

    RUN(planet_moon_ratio('g', 0) > 0.0, "Io ratio valid");
    RUN(planet_moon_ratio('g', 4) < 0.0, "Galilean moon 4 invalid");
    RUN(planet_moon_ratio('s', 7) > 0.0, "Hyperion ratio valid");
    RUN(planet_moon_ratio('s', 8) < 0.0, "Saturn moon 8 invalid");
    RUN(planet_moon_ratio('u', 4) > 0.0, "Oberon ratio valid");
    RUN(planet_moon_ratio('u', 5) < 0.0, "Uranus moon 5 invalid");
    RUN(planet_moon_ratio('m', 1) > 0.0, "Deimos ratio valid");
    RUN(planet_moon_ratio('m', 2) < 0.0, "Mars moon 2 invalid");
    RUN(planet_moon_ratio('x', 0) < 0.0, "unknown family invalid");
}


static void
test_all_planets(void)
{
    int body;

    fprintf(stderr, "\n── All planets equilibrium ──\n");

    for (body = 1; body <= 8; body++)
    {
        double ratio = sun_planet_ratio(body);
        double mu    = mu_from_ratio(ratio);
        char   label[64];
        int    pt;

        for (pt = LAGRANGE_L1; pt <= LAGRANGE_L5; pt++)
        {
            snprintf(label, sizeof(label), "body %d L%d", body, pt);
            test_equilibrium_check(mu, pt, label);
        }
    }
}


/* ── Main ──────────────────────────────────────────────── */

int
main(void)
{
    fprintf(stderr, "Lagrange solver unit test\n");
    fprintf(stderr, "========================\n");

    test_sun_earth();
    test_sun_jupiter();
    test_earth_moon();
    test_l4_l5_symmetry();
    test_l1_l2_ordering();
    test_hill_radius();
    test_zone_radius();
    test_physical_transform();
    test_extreme_mass_ratios();
    test_error_cases();
    test_all_planets();

    fprintf(stderr, "\n%d/%d tests passed\n", n_pass, n_run);

    return (n_pass == n_run) ? 0 : 1;
}