Fix planet_magnitude: use full Mallama & Hilton (2018) polynomials

The original implementation only used c1+c2 coefficients from the
simplified model. Mercury's 6th-order polynomial diverges badly
beyond ~60 deg phase angle with only 2 terms — returning -23 mag
at 130 deg (should be +1.1).

Now uses the complete piecewise models from the paper:
- Mercury: full 6th-order polynomial (Eq. 1)
- Venus: piecewise at 163.7 deg (Eq. 2/3)
- Mars: piecewise at 50 deg (Eq. 5/6)
- Jupiter: piecewise at 12 deg with log term (Eq. 7/8)
- Saturn: globe-only model (Eq. 11/12), ring tilt still not modeled
- Uranus: phase threshold at 3.1 deg (Eq. 14)
- Neptune: phase threshold at 1.9 deg (Eq. 17)

Bug found by astrolock: Mercury at superior conjunction (i=130.6 deg)
returned -23.3 mag instead of +1.1.

src/magnitude_funcs.c

@@ -20,35 +20,99 @@ PG_FUNCTION_INFO_V1(planet_magnitude);
 
 
 /*
- * Planet magnitude parameters -- Mallama & Hilton (2018), simplified.
+ * Per-planet phase correction -- Mallama & Hilton (2018).
  *
- * V(1,0) = absolute magnitude at r=1 AU, delta=1 AU, i=0 deg
- * Phase corrections are polynomial fits to i (phase angle in degrees).
+ * 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.
  *
- * We use the linear+quadratic terms which are sufficient for
- * phase angles encountered from Earth (typically <180 deg).
- *
- * Saturn caveat: visual magnitude depends heavily on ring tilt
- * (can vary by ~1.5 mag).  The simplified model here uses a fixed
- * V(1,0) without ring correction.
+ * Saturn caveat: ring tilt contribution (their Eq. 10) requires
+ * saturnicentric sub-observer latitude, which we don't compute.
+ * We use the globe-only model (Eq. 11/12) — error up to ~1.5 mag.
  */
-typedef struct {
-	double	v10;	/* V(1,0) */
-	double	c1;		/* coefficient for i */
-	double	c2;		/* coefficient for i^2 */
-	double	c3;		/* coefficient for i^3 (0 if unused) */
-} mag_params;
 
-static const mag_params planet_mag[] = {
-	[0] = { 0, 0, 0, 0 },								/* Sun: unused placeholder */
-	[1] = { -0.613, 6.328e-2, -1.6336e-3, 0 },		/* Mercury */
-	[2] = { -4.384, 1.044e-3, 3.687e-4, 0 },			/* Venus */
-	[3] = { 0, 0, 0, 0 },								/* Earth: unused */
-	[4] = { -1.601, 2.267e-2, -1.302e-4, 0 },			/* Mars */
-	[5] = { -9.395, 3.7e-4, 0, 0 },					/* Jupiter */
-	[6] = { -8.95, 0, 0, 0 },							/* Saturn (ring tilt NOT modeled) */
-	[7] = { -7.110, 0, 0, 0 },							/* Uranus */
-	[8] = { -7.00, 0, 0, 0 },							/* Neptune */
+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-only (Eq. 11), no ring tilt */
+			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-only) */
+	[7] = -7.110,		/* Uranus */
+	[8] = -7.00,		/* Neptune */
 };
 
 
@@ -65,7 +129,6 @@ compute_planet_magnitude(int body_id, double jd)
 	double	geo[3];
 	double	r, delta, R;
 	double	cos_i, i_deg;
-	const mag_params *p;
 	double	V;
 	int		vsop_body = body_id - 1;	/* pg_orrery 1-based -> VSOP87 0-based */
 
@@ -94,13 +157,10 @@ compute_planet_magnitude(int body_id, double jd)
 	if (cos_i < -1.0) cos_i = -1.0;
 	i_deg = acos(cos_i) * RAD_TO_DEG;
 
-	/* Mallama & Hilton (2018) magnitude formula */
-	p = &planet_mag[body_id];
-	V = p->v10
+	/* Mallama & Hilton (2018) magnitude with full phase correction */
+	V = planet_v10[body_id]
 		+ 5.0 * log10(r * delta)
-		+ p->c1 * i_deg
-		+ p->c2 * i_deg * i_deg
-		+ p->c3 * i_deg * i_deg * i_deg;
+		+ phase_correction(body_id, i_deg);
 
 	return V;
 }