/*
 * magnitude_funcs.c -- Planet magnitude, solar elongation, phase fraction
 *
 * Uses the Mallama & Hilton (2018) magnitude model with
 * VSOP87 positions for distances and phase angles.
 *
 * Solar elongation and planet phase reuse the same Sun-Planet-Earth
 * triangle geometry, factored into compute_planet_geometry().
 *
 * Reference: Mallama & Hilton, "Computing Apparent Planetary
 * Magnitudes for The Astronomical Almanac", A&C vol. 25, 2018.
 */

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

PG_FUNCTION_INFO_V1(planet_magnitude);
PG_FUNCTION_INFO_V1(solar_elongation);
PG_FUNCTION_INFO_V1(planet_phase);
PG_FUNCTION_INFO_V1(saturn_ring_tilt);


/*
 * Per-planet phase correction -- Mallama & Hilton (2018).
 *
 * Mercury uses a 6th-order polynomial (their Eq. 1).
 * Venus and Mars are piecewise with different coefficients for
 * small vs large phase angles.  Jupiter is piecewise at 12 deg.
 * Saturn, Uranus, Neptune use simpler models.
 *
 * Saturn: globe model (Eq. 11/12) plus ring tilt correction (Eq. 10)
 * using IAU 2000 Saturn pole direction for sub-observer latitude B'.
 */

static double
phase_correction(int body_id, double i)
{
	double i2 = i * i;

	switch (body_id)
	{
		case 1:		/* Mercury: 6th-order polynomial */
			return i * (6.3280e-02
				+ i * (-1.6336e-03
				+ i * (3.3644e-05
				+ i * (-3.4265e-07
				+ i * (1.6893e-09
				+ i * (-3.0334e-12))))));

		case 2:		/* Venus: piecewise at 163.7 deg */
			if (i < 163.7)
				return i * (-1.044e-03
					+ i * (3.687e-04
					+ i * (-2.814e-06
					+ i * 8.938e-09)));
			else
				return (236.05828 + 4.384) + i * (-2.81914e+00
					+ i * 8.39034e-03);

		case 4:		/* Mars: piecewise at 50 deg */
			if (i <= 50.0)
				return i * (2.267e-02 + i * (-1.302e-04));
			else
				return (-1.601 + 0.367) + i * (-0.02573 + i * 0.0003445);

		case 5:		/* Jupiter: piecewise at 12 deg */
			if (i <= 12.0)
				return i * (6.16e-04 * i - 3.7e-04);
			else
			{
				double a = i / 180.0;
				return (-9.428 + 9.395) + (-2.5)
					* log10(1.0 - 1.507 * a - 0.363 * a * a
						- 0.062 * a * a * a
						+ 2.809 * a * a * a * a
						- 1.876 * a * a * a * a * a);
			}

		case 6:		/* Saturn: globe phase (Eq. 11/12), ring tilt added in planet_magnitude() */
			if (i <= 6.5)
				return -3.7e-04 * i + 6.16e-04 * i2;
			else
				return 2.446e-04 * i + 2.672e-04 * i2
					- 1.506e-06 * i2 * i + 4.767e-09 * i2 * i2;

		case 7:		/* Uranus */
			if (i <= 3.1)
				return 0.0;
			return i * (6.587e-03 + i * 1.045e-04);

		case 8:		/* Neptune */
			if (i <= 1.9)
				return 0.0;
			return i * (7.944e-03 + i * 9.617e-05);

		default:
			return 0.0;
	}
}


/*
 * V(1,0) per planet -- absolute magnitude at unit distances, zero phase.
 * Mercury through Neptune.  Mars piecewise handled in phase_correction().
 */
static const double planet_v10[] = {
	[0] =  0.0,		/* Sun: unused */
	[1] = -0.613,		/* Mercury */
	[2] = -4.384,		/* Venus */
	[3] =  0.0,		/* Earth: unused */
	[4] = -1.601,		/* Mars (i <= 50; piecewise shifts in phase_correction) */
	[5] = -9.395,		/* Jupiter (i <= 12; piecewise shifts in phase_correction) */
	[6] = -8.95,		/* Saturn (globe + ring) */
	[7] = -7.110,		/* Uranus */
	[8] = -7.00,		/* Neptune */
};


/*
 * Shared Sun-Planet-Earth geometry.
 *
 * Computes the three distances (r, delta, R) and the phase angle
 * (Sun-Planet-Earth angle) from VSOP87 positions.  Used by
 * magnitude, elongation, and phase functions.
 */
typedef struct
{
	double r;       /* Sun-Planet distance (AU) */
	double delta;   /* Earth-Planet distance (AU) */
	double R;       /* Sun-Earth distance (AU) */
	double i_deg;   /* Phase angle, degrees (Sun-Planet-Earth vertex) */
	double gv[3];   /* geocentric ecliptic J2000 (AU) — for Saturn ring tilt */
} planet_geometry;

static void
compute_planet_geometry(int body_id, double jd, planet_geometry *geo)
{
	double	earth_xyz[6], planet_xyz[6];
	double	gv[3];
	double	cos_i;
	int		vsop_body = body_id - 1;	/* pg_orrery 1-based -> VSOP87 0-based */

	GetVsop87Coor(jd, 2, earth_xyz);			/* Earth (VSOP87 body 2) */
	GetVsop87Coor(jd, vsop_body, planet_xyz);	/* target planet */

	/* Heliocentric distance to planet */
	geo->r = sqrt(planet_xyz[0] * planet_xyz[0] +
				  planet_xyz[1] * planet_xyz[1] +
				  planet_xyz[2] * planet_xyz[2]);

	/* Geocentric vector and distance */
	gv[0] = planet_xyz[0] - earth_xyz[0];
	gv[1] = planet_xyz[1] - earth_xyz[1];
	gv[2] = planet_xyz[2] - earth_xyz[2];
	geo->delta = sqrt(gv[0] * gv[0] + gv[1] * gv[1] + gv[2] * gv[2]);
	geo->gv[0] = gv[0];
	geo->gv[1] = gv[1];
	geo->gv[2] = gv[2];

	/* Sun-Earth distance */
	geo->R = sqrt(earth_xyz[0] * earth_xyz[0] +
				  earth_xyz[1] * earth_xyz[1] +
				  earth_xyz[2] * earth_xyz[2]);

	/* Phase angle via law of cosines: vertex at planet */
	cos_i = (geo->r * geo->r + geo->delta * geo->delta - geo->R * geo->R)
			/ (2.0 * geo->r * geo->delta);
	if (cos_i > 1.0) cos_i = 1.0;
	if (cos_i < -1.0) cos_i = -1.0;
	geo->i_deg = acos(cos_i) * RAD_TO_DEG;
}


/*
 * Saturn pole direction in ecliptic J2000.
 *
 * IAU 2000 pole: RA0 = 40.589 deg, Dec0 = 83.537 deg (equatorial J2000).
 * Converted to ecliptic J2000 via rotation by obliquity.
 *
 *   ecl_x = cos(dec)*cos(ra)
 *   ecl_y = cos(dec)*sin(ra)*cos(eps) + sin(dec)*sin(eps)
 *   ecl_z = -cos(dec)*sin(ra)*sin(eps) + sin(dec)*cos(eps)
 *
 * Pre-computed unit vector (constant across timescales relevant here).
 */
static const double saturn_pole_ecl[3] = {
	0.08547883,     /* x: cos(83.537)*cos(40.589) */
	0.46244181,     /* y: cos(83.537)*sin(40.589)*cos(23.4393) + sin(83.537)*sin(23.4393) */
	0.88251965      /* z: -cos(83.537)*sin(40.589)*sin(23.4393) + sin(83.537)*cos(23.4393) */
};

/*
 * Compute sub-observer latitude of Earth relative to Saturn's ring plane.
 *
 * B' = arcsin(dot(geocentric_unit_vector, saturn_pole_ecl))
 *
 * When |B'| is large, rings are maximally tilted toward Earth (brighter).
 * When B' ~ 0, rings are edge-on (dimmest, nearly invisible).
 * Range: [-27, +27] deg (Saturn's axial tilt is 26.73 deg).
 */
static double
compute_ring_tilt(const double gv[3], double delta)
{
	double gv_unit[3];
	double dot;

	gv_unit[0] = gv[0] / delta;
	gv_unit[1] = gv[1] / delta;
	gv_unit[2] = gv[2] / delta;

	dot = gv_unit[0] * saturn_pole_ecl[0] +
		  gv_unit[1] * saturn_pole_ecl[1] +
		  gv_unit[2] * saturn_pole_ecl[2];

	if (dot > 1.0) dot = 1.0;
	if (dot < -1.0) dot = -1.0;

	return asin(dot);  /* radians */
}


/*
 * Validate planet body_id for magnitude/elongation/phase.
 * Must be 1-8 (Mercury-Neptune), not 3 (Earth).
 */
static void
validate_planet_body_id(int body_id, const char *func_name)
{
	if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
		ereport(ERROR,
				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
				 errmsg("%s: body_id %d must be 1-8 (Mercury-Neptune)",
						func_name, body_id)));

	if (body_id == BODY_EARTH)
		ereport(ERROR,
				(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
				 errmsg("%s: cannot compute for Earth from Earth",
						func_name)));
}


/* ================================================================
 * planet_magnitude(body_id int4, timestamptz) -> float8
 *
 * Returns the apparent visual magnitude of a planet as seen from
 * Earth.  Uses Mallama & Hilton (2018) magnitude model.
 *
 * Body IDs: 1=Mercury, ..., 8=Neptune (not Sun 0, Earth 3, or Moon 10)
 *
 * Saturn includes ring tilt correction (Eq. 10) using the IAU 2000
 * pole direction and VSOP87 geometry.
 * ================================================================
 */
Datum
planet_magnitude(PG_FUNCTION_ARGS)
{
	int32			body_id = PG_GETARG_INT32(0);
	int64			ts		= PG_GETARG_INT64(1);
	double			jd;
	planet_geometry	geo;
	double			V;

	validate_planet_body_id(body_id, "planet_magnitude");

	jd = timestamptz_to_jd(ts);
	compute_planet_geometry(body_id, jd, &geo);

	V = planet_v10[body_id]
		+ 5.0 * log10(geo.r * geo.delta)
		+ phase_correction(body_id, geo.i_deg);

	/* Saturn ring tilt correction -- Mallama & Hilton (2018) Eq. 10 */
	if (body_id == BODY_SATURN)
	{
		double Bp = compute_ring_tilt(geo.gv, geo.delta);
		double sin_Bp = fabs(sin(Bp));
		double sin2_Bp = sin(Bp) * sin(Bp);
		V += -2.60 * sin_Bp + 1.25 * sin2_Bp;
	}

	PG_RETURN_FLOAT8(V);
}


/* ================================================================
 * solar_elongation(body_id int4, timestamptz) -> float8
 *
 * Sun-Earth-Planet angle in degrees [0, 180].
 * How far a planet appears from the Sun in the sky.
 *
 * Uses law of cosines with vertex at Earth:
 *   cos(elong) = (R^2 + delta^2 - r^2) / (2 * R * delta)
 *
 * Mercury max ~28 deg, Venus max ~47 deg.
 * Superior planets can reach ~180 deg (opposition).
 * ================================================================
 */
Datum
solar_elongation(PG_FUNCTION_ARGS)
{
	int32			body_id = PG_GETARG_INT32(0);
	int64			ts		= PG_GETARG_INT64(1);
	double			jd;
	planet_geometry	geo;
	double			cos_elong, elong_deg;

	validate_planet_body_id(body_id, "solar_elongation");

	jd = timestamptz_to_jd(ts);
	compute_planet_geometry(body_id, jd, &geo);

	/* Law of cosines, vertex at Earth */
	cos_elong = (geo.R * geo.R + geo.delta * geo.delta - geo.r * geo.r)
				/ (2.0 * geo.R * geo.delta);
	if (cos_elong > 1.0) cos_elong = 1.0;
	if (cos_elong < -1.0) cos_elong = -1.0;
	elong_deg = acos(cos_elong) * RAD_TO_DEG;

	PG_RETURN_FLOAT8(elong_deg);
}


/* ================================================================
 * planet_phase(body_id int4, timestamptz) -> float8
 *
 * Illuminated fraction of a planet's disk as seen from Earth [0, 1].
 *   k = (1 + cos(i)) / 2
 * where i is the phase angle (Sun-Planet-Earth).
 *
 * Inner planets vary dramatically (Venus crescent).
 * Outer planets are always near 1.0.
 * ================================================================
 */
Datum
planet_phase(PG_FUNCTION_ARGS)
{
	int32			body_id = PG_GETARG_INT32(0);
	int64			ts		= PG_GETARG_INT64(1);
	double			jd;
	planet_geometry	geo;
	double			i_rad, k;

	validate_planet_body_id(body_id, "planet_phase");

	jd = timestamptz_to_jd(ts);
	compute_planet_geometry(body_id, jd, &geo);

	i_rad = geo.i_deg * DEG_TO_RAD;
	k = (1.0 + cos(i_rad)) / 2.0;

	PG_RETURN_FLOAT8(k);
}


/* ================================================================
 * saturn_ring_tilt(timestamptz) -> float8
 *
 * Sub-observer latitude of Earth relative to Saturn's ring plane,
 * in degrees [-27, +27].  Indicates how much the rings are tilted
 * toward Earth at the given time.
 *
 * Near 0: rings edge-on (ring crossing events, e.g. 2025 March).
 * Near +/-27: rings maximally open (brightest configuration).
 *
 * Uses IAU 2000 Saturn pole and VSOP87 Earth-Saturn geometry.
 * ================================================================
 */
Datum
saturn_ring_tilt(PG_FUNCTION_ARGS)
{
	int64			ts = PG_GETARG_INT64(0);
	double			jd;
	planet_geometry	geo;

	jd = timestamptz_to_jd(ts);
	compute_planet_geometry(BODY_SATURN, jd, &geo);

	PG_RETURN_FLOAT8(compute_ring_tilt(geo.gv, geo.delta) * RAD_TO_DEG);
}