fix: use WGS-84 ellipsoid for observer positioning instead of spherical Earth

Observer ECEF/ECI conversion treated geodetic lat/lon as geocentric,
placing the observer up to ~21.7 km off at mid-latitudes. This caused
up to ~1° azimuth error on east/west LEO passes and ~2.3% slant range
error. Add shared geodetic helpers (geodeticToEcef, geodeticToEci) using
WGS-84 prime vertical radius of curvature and replace spherical formula
in az-el, doppler, and magnitude modules.

Author: sandrwich

src/astro/az-el.ts

@@ -2,52 +2,49 @@
  * Compute azimuth and elevation from a ground observer to a satellite.
  * Pure math — no THREE.js or DOM dependencies. Usable in Web Workers.
  *
- * Uses spherical Earth model (R = 6371 km).
+ * Uses WGS-84 ellipsoid for observer positioning.
  *
  * @param eciX  Standard ECI X position (km) — toward vernal equinox
  * @param eciY  Standard ECI Y position (km) — completes right-hand system
  * @param eciZ  Standard ECI Z position (km) — toward north pole
  * @param gmstRad  Greenwich Mean Sidereal Time in radians
- * @param obsLatDeg  Observer latitude in degrees (-90 to 90)
- * @param obsLonDeg  Observer longitude in degrees (-180 to 180)
- * @param obsAltM  Observer altitude in meters above sea level
+ * @param obsLatDeg  Observer geodetic latitude in degrees (-90 to 90)
+ * @param obsLonDeg  Observer geodetic longitude in degrees (-180 to 180)
+ * @param obsAltM  Observer altitude in meters above WGS-84 ellipsoid
  * @returns { az, el } in degrees. Az: 0=N, 90=E, 180=S, 270=W. El: 0=horizon, 90=zenith.
  */
+import { geodeticToEcef } from './geodetic';
+
 export function getAzEl(
   eciX: number, eciY: number, eciZ: number,
   gmstRad: number,
   obsLatDeg: number, obsLonDeg: number, obsAltM: number,
 ): { az: number; el: number } {
   const DEG2RAD = Math.PI / 180;
-  const EARTH_R = 6371.0; // km
 
   const satR = Math.sqrt(eciX * eciX + eciY * eciY + eciZ * eciZ);
   if (satR === 0) return { az: 0, el: -90 };
 
-  // Satellite geodetic from ECI
+  // Satellite geocentric lat/lon from ECI, then to ECEF
   const satLat = Math.asin(eciZ / satR);
   const satLonEci = Math.atan2(eciY, eciX);
   const satLonEcef = satLonEci - gmstRad;
 
-  // Satellite ECEF position
   const sx = satR * Math.cos(satLat) * Math.cos(satLonEcef);
   const sy = satR * Math.cos(satLat) * Math.sin(satLonEcef);
   const sz = satR * Math.sin(satLat);
 
-  // Observer ECEF position
-  const latRad = obsLatDeg * DEG2RAD;
-  const lonRad = obsLonDeg * DEG2RAD;
-  const obsR = EARTH_R + obsAltM / 1000;
-  const ox = obsR * Math.cos(latRad) * Math.cos(lonRad);
-  const oy = obsR * Math.cos(latRad) * Math.sin(lonRad);
-  const oz = obsR * Math.sin(latRad);
+  // Observer ECEF position (WGS-84 ellipsoid)
+  const obs = geodeticToEcef(obsLatDeg, obsLonDeg, obsAltM);
 
   // Range vector (ECEF)
-  const dx = sx - ox;
-  const dy = sy - oy;
-  const dz = sz - oz;
+  const dx = sx - obs.x;
+  const dy = sy - obs.y;
+  const dz = sz - obs.z;
 
-  // Rotate to topocentric East-North-Up
+  // Rotate to topocentric East-North-Up (geodetic lat/lon defines the local frame)
+  const latRad = obsLatDeg * DEG2RAD;
+  const lonRad = obsLonDeg * DEG2RAD;
   const clat = Math.cos(latRad), slat = Math.sin(latRad);
   const clon = Math.cos(lonRad), slon = Math.sin(lonRad);
 

src/astro/doppler.ts

@@ -7,7 +7,8 @@
 import { twoline2satrec, json2satrec, propagate } from 'satellite.js';
 import type { SatRec } from 'satellite.js';
 import { epochToUnix, epochToGmst } from './epoch';
-import { EARTH_RADIUS_KM, DEG2RAD } from '../constants';
+import { DEG2RAD } from '../constants';
+import { geodeticToEcef } from './geodetic';
 
 const C_KM_S = 299792.458;        // speed of light (km/s)
 const OMEGA_E = 7.2921159e-5;     // Earth rotation rate (rad/s)
@@ -18,18 +19,6 @@ export interface DopplerResult {
   rangeRateKmS: number;   // range rate (km/s, positive = receding)
 }
 
-/** Observer ECEF position from geodetic (spherical Earth). */
-function observerEcef(latDeg: number, lonDeg: number, altM: number) {
-  const lat = latDeg * DEG2RAD;
-  const lon = lonDeg * DEG2RAD;
-  const R = EARTH_RADIUS_KM + altM / 1000;
-  return {
-    x: R * Math.cos(lat) * Math.cos(lon),
-    y: R * Math.cos(lat) * Math.sin(lon),
-    z: R * Math.sin(lat),
-  };
-}
-
 /**
  * Rotate a standard-ECI vector to ECEF via Rz(-GMST).
  * Standard ECI: x=vernal equinox, y=90° ahead, z=north pole.
@@ -83,8 +72,8 @@ export function calculateDopplerShift(
     z: velRot.z,
   };
 
-  // Observer (stationary in ECEF)
-  const obs = observerEcef(obsLatDeg, obsLonDeg, obsAltM);
+  // Observer (stationary in ECEF, WGS-84 ellipsoid)
+  const obs = geodeticToEcef(obsLatDeg, obsLonDeg, obsAltM);
 
   // Range vector (observer → satellite)
   const dx = satPos.x - obs.x;

src/astro/geodetic.ts

@@ -0,0 +1,57 @@
+/**
+ * WGS-84 geodetic coordinate conversions.
+ * Pure math — no THREE.js dependency. Usable in Web Workers.
+ */
+
+import { DEG2RAD, WGS84_A, WGS84_E2 } from '../constants';
+
+/**
+ * Convert geodetic coordinates to ECEF (Earth-Centered, Earth-Fixed) in km.
+ * Uses the WGS-84 ellipsoid for accurate positioning.
+ *
+ * @param latDeg  Geodetic latitude (degrees)
+ * @param lonDeg  Geodetic longitude (degrees)
+ * @param altM    Altitude above ellipsoid (meters)
+ */
+export function geodeticToEcef(
+  latDeg: number, lonDeg: number, altM: number,
+): { x: number; y: number; z: number } {
+  const latRad = latDeg * DEG2RAD;
+  const lonRad = lonDeg * DEG2RAD;
+  const sinLat = Math.sin(latRad);
+  const cosLat = Math.cos(latRad);
+  const altKm = altM / 1000;
+  const N = WGS84_A / Math.sqrt(1 - WGS84_E2 * sinLat * sinLat);
+  return {
+    x: (N + altKm) * cosLat * Math.cos(lonRad),
+    y: (N + altKm) * cosLat * Math.sin(lonRad),
+    z: (N * (1 - WGS84_E2) + altKm) * sinLat,
+  };
+}
+
+/**
+ * Convert geodetic coordinates to standard ECI (Earth-Centered Inertial) in km.
+ * Rotates from ECEF by GMST to align with the vernal equinox frame.
+ *
+ * Standard ECI: x = vernal equinox, y = 90° ahead, z = north pole.
+ *
+ * @param latDeg   Geodetic latitude (degrees)
+ * @param lonDeg   Geodetic longitude (degrees)
+ * @param altM     Altitude above ellipsoid (meters)
+ * @param gmstRad  Greenwich Mean Sidereal Time (radians)
+ */
+export function geodeticToEci(
+  latDeg: number, lonDeg: number, altM: number, gmstRad: number,
+): { x: number; y: number; z: number } {
+  const latRad = latDeg * DEG2RAD;
+  const sinLat = Math.sin(latRad);
+  const cosLat = Math.cos(latRad);
+  const altKm = altM / 1000;
+  const N = WGS84_A / Math.sqrt(1 - WGS84_E2 * sinLat * sinLat);
+  const theta = gmstRad + lonDeg * DEG2RAD;
+  return {
+    x: (N + altKm) * cosLat * Math.cos(theta),
+    y: (N + altKm) * cosLat * Math.sin(theta),
+    z: (N * (1 - WGS84_E2) + altKm) * sinLat,
+  };
+}

src/astro/magnitude.ts

@@ -13,9 +13,10 @@
  * and cached in localStorage. See .docs/magnitude.md for details.
  */
 
+import { geodeticToEci } from './geodetic';
+
 const DEG2RAD = Math.PI / 180.0;
 const RAD2DEG = 180.0 / Math.PI;
-const EARTH_RADIUS_KM = 6371.0;
 
 /**
  * Lambertian (diffuse) sphere phase function.
@@ -74,27 +75,17 @@ export function computePhaseAngle(
 /**
  * Compute observer position in **standard ECI** (km) from geodetic coordinates.
  * Returns standard ECI — convert to render coords (x, z, -y) before mixing with sat.currentPos.
- * Simplified: assumes spherical Earth (good enough for magnitude estimation).
+ * Uses WGS-84 ellipsoid for accurate positioning.
  *
  * @param latDeg  Geodetic latitude (degrees)
  * @param lonDeg  Geodetic longitude (degrees)
- * @param altM    Altitude above sea level (meters)
+ * @param altM    Altitude above WGS-84 ellipsoid (meters)
  * @param gmstRad Greenwich Mean Sidereal Time (radians)
  */
 export function observerEci(
   latDeg: number, lonDeg: number, altM: number, gmstRad: number,
 ): { x: number; y: number; z: number } {
-  const latRad = latDeg * DEG2RAD;
-  const lonRad = lonDeg * DEG2RAD;
-  const r = EARTH_RADIUS_KM + altM / 1000.0;
-  const cosLat = Math.cos(latRad);
-  const sinLat = Math.sin(latRad);
-  const theta = gmstRad + lonRad;  // sidereal angle
-  return {
-    x: r * cosLat * Math.cos(theta),
-    y: r * cosLat * Math.sin(theta),
-    z: r * sinLat,
-  };
+  return geodeticToEci(latDeg, lonDeg, altM, gmstRad);
 }
 
 /**

src/constants.ts

@@ -50,5 +50,10 @@ export function satColorGl(index: number): readonly [number, number, number] {
 // J2 perturbation constants
 export const J2 = 1.08263e-3;                  // Earth's J2 zonal harmonic
 export const EARTH_RADIUS_EQ_KM = 6378.137;    // WGS-84 equatorial radius (km)
+
+// WGS-84 ellipsoid constants
+export const WGS84_A = 6378.137;               // semi-major axis (km)
+export const WGS84_F = 1 / 298.257223563;      // flattening
+export const WGS84_E2 = 2 * WGS84_F - WGS84_F * WGS84_F; // first eccentricity squared
 export const ORBIT_RECOMPUTE_INTERVAL_S = 900;  // recompute orbits every 15 sim-minutes
 export const MOBILE_BREAKPOINT = 768;