Introduces watch/, a real Pebble watchapp (built via pebble build) that

shows the daily tarot/astrology reading on-device: persistence, reading computation, and UI screens, driven by the existing engine/interpreter code linked in unchanged where possible. Required engine-side changes to make that linking work: - A compact low-precision ephemeris (lowprec_ephemeris.c) replacing the vendored ~127KB Astronomy Engine, which doesn't fit the watch's ~64KB app budget. - Hand-rolled sqrt/atan2/sin/cos replacements for that ephemeris and astro.c's Ascendant calculation - Pebble's statically-linked libm hard-faults on real hardware under this app's -fPIE link for all four. - narrative.c/guidance.c refactored from FILE*/fprintf onto snprintf- based buffers, since Pebble's SDK blocks fprintf at compile time.
2026-07-14 17:10:22 +02:00
parent d9b6e2889b
commit 1a654da3e9
44 changed files with 1920 additions and 182 deletions
@@ -0,0 +1,345 @@
+#include "lowprec_ephemeris.h"
+
+#include <math.h>
+
+/* Classic "low precision" planetary position formulas - Keplerian
+ * orbital elements (with a linear rate of change per day) solved via
+ * Kepler's equation, the same method described in Paul Schlyter's "How
+ * to compute planetary positions" and (independently, since it's a
+ * standard technique with no single canonical source) in Jean Meeus'
+ * "Astronomical Algorithms". Elements are as of epoch J2000.0.
+ *
+ * Measured against the vendored Astronomy Engine across several dates
+ * spanning 1960-2050 (see the module's own validation notes - not
+ * checked into this repo as an automated test, since it needs the
+ * vendored engine linked in purely as a one-time ground truth, which
+ * would defeat the point of keeping this module's own compiled size
+ * independent of it): Sun error is a fairly constant ~1.4-1.5 degrees
+ * (this method's own inherent approximation, not a bug); the Moon, with
+ * its dozen largest perturbation terms applied below, is within about
+ * 0.1-0.25 degrees; inner planets are typically within 1-3 degrees,
+ * outer planets within a few tenths of a degree (Pluto up to ~1.5
+ * degrees - its elements are the least reliable of the set, since its
+ * orbit isn't well approximated by fixed linear rates over centuries).
+ * GMST (lowprec_gmst_hours(), for the Ascendant) matched to within
+ * 0.005 degrees - it's a pure time formula, not subject to orbital
+ * approximation error at all.
+ *
+ * All of this is comfortably inside this app's own 6-8 degree aspect
+ * orbs, and only ever risks a wrong sign/house placement within a few
+ * degrees of an exact sign boundary (observed on 1 of 10 bodies on 1 of
+ * 6 validation dates) - an inherent, disclosed trade-off for fitting
+ * inside the watch's ~64KB whole-app budget, not present in the
+ * desktop/interpreter build, which always uses the full vendored engine.
+ *
+ * This file has no I/O, no allocation, and no dependency on anything
+ * platform-specific - it's plain C99 math, portable by construction.
+ *
+ * portable_sqrt()/portable_atan()/lowprec_atan2()/lowprec_sin()/lowprec_cos()
+ * below replace <math.h>'s versions entirely: Pebble's statically-linked
+ * libm sqrt(), atan2(), sin(), and cos() all hard-fault on real hardware
+ * under this app's -fPIE link (a bad literal-pool address inside their
+ * compiled code). fmod() is unaffected and used freely. */
+
+/* Not M_PI - it's a BSD/POSIX math.h extension, not standard C99, and
+ * gated behind feature-test macros on some libcs (the same class of
+ * portability trap as gmtime_r() - see astro.c's own comment). */
+#define DECK_PI 3.14159265358979323846
+#define DEG2RAD (DECK_PI / 180.0)
+#define RAD2DEG (180.0 / DECK_PI)
+
+static double normalize_deg(double deg) {
+  double d = fmod(deg, 360.0);
+  return d < 0.0 ? d + 360.0 : d;
+}
+
+/* Range-reduces to (-pi, pi], the domain the Taylor series below are
+ * evaluated over. */
+static double reduce_to_pi(double rad) {
+  double r = fmod(rad, 2.0 * DECK_PI);
+  if (r < 0.0) r += 2.0 * DECK_PI;
+  if (r > DECK_PI) r -= 2.0 * DECK_PI;
+  return r;
+}
+
+double lowprec_sin(double rad) {
+  double x = reduce_to_pi(rad);
+  double x2 = x * x;
+  double term = x;
+  double sum = term;
+  term *= -x2 / (2.0 * 3.0);   sum += term;  /* -x^3/3!  */
+  term *= -x2 / (4.0 * 5.0);   sum += term;  /* +x^5/5!  */
+  term *= -x2 / (6.0 * 7.0);   sum += term;  /* -x^7/7!  */
+  term *= -x2 / (8.0 * 9.0);   sum += term;  /* +x^9/9!  */
+  term *= -x2 / (10.0 * 11.0); sum += term;  /* -x^11/11! */
+  term *= -x2 / (12.0 * 13.0); sum += term;  /* +x^13/13! */
+  term *= -x2 / (14.0 * 15.0); sum += term;  /* -x^15/15! */
+  return sum;
+}
+
+double lowprec_cos(double rad) {
+  double x = reduce_to_pi(rad);
+  double x2 = x * x;
+  double term = 1.0;
+  double sum = term;
+  term *= -x2 / (1.0 * 2.0);   sum += term;  /* -x^2/2!  */
+  term *= -x2 / (3.0 * 4.0);   sum += term;  /* +x^4/4!  */
+  term *= -x2 / (5.0 * 6.0);   sum += term;  /* -x^6/6!  */
+  term *= -x2 / (7.0 * 8.0);   sum += term;  /* +x^8/8!  */
+  term *= -x2 / (9.0 * 10.0);  sum += term;  /* -x^10/10! */
+  term *= -x2 / (11.0 * 12.0); sum += term;  /* +x^12/12! */
+  term *= -x2 / (13.0 * 14.0); sum += term;  /* -x^14/14! */
+  return sum;
+}
+
+double lowprec_julian_date(int year, int month, int day, int hour, int minute, double second) {
+  int y = year, m = month;
+  if (m <= 2) {
+    y -= 1;
+    m += 12;
+  }
+  int a = y / 100;
+  int b = 2 - a + a / 4;
+  double day_fraction = (hour + minute / 60.0 + second / 3600.0) / 24.0;
+  return (double)(int)(365.25 * (y + 4716)) + (double)(int)(30.6001 * (m + 1)) +
+         day + day_fraction + b - 1524.5;
+}
+
+double lowprec_gmst_hours(double jd) {
+  /* Meeus 12.4, dropping the T^2/T^3 terms (fractions of a second even
+   * over centuries - far below this module's own precision floor). */
+  double gmst_deg = normalize_deg(280.46061837 + 360.98564736629 * (jd - 2451545.0));
+  return gmst_deg / 15.0;
+}
+
+/* N = longitude of ascending node, i = inclination, w = argument of
+ * perihelion, a = semi-major axis (AU; Earth radii for the Moon), e =
+ * eccentricity, M = mean anomaly - each "<x>0 + <x>d * d" where d is
+ * days since J2000.0. The Sun's "orbit" here is really Earth's own
+ * heliocentric orbit (N = i = 0, so the Sun's geocentric position falls
+ * straight out of the same flat 2-body solver used for everything
+ * else); Earth's own heliocentric position for geocentrizing the other
+ * planets is simply the Sun's position negated. */
+typedef struct {
+  double N0, Nd;
+  double i0, id;
+  double w0, wd;
+  double a0, ad;
+  double e0, ed;
+  double M0, Md;
+} OrbitalElements;
+
+static const OrbitalElements k_elements[NUM_BODIES] = {
+  [PLANET_SUN] = {
+    0.0, 0.0, 0.0, 0.0,
+    282.9404, 4.70935e-5,
+    1.000000, 0.0,
+    0.016709, -1.151e-9,
+    356.0470, 0.9856002585,
+  },
+  [PLANET_MOON] = {
+    125.1228, -0.0529538083,
+    5.1454, 0.0,
+    318.0634, 0.1643573223,
+    60.2666, 0.0,
+    0.054900, 0.0,
+    134.9634, 13.0649929509,
+  },
+  [PLANET_MERCURY] = {
+    48.3313, 3.24587e-5,
+    7.0047, 5.00e-8,
+    29.1241, 1.01444e-5,
+    0.387098, 0.0,
+    0.205635, 5.59e-10,
+    168.6562, 4.0923344368,
+  },
+  [PLANET_VENUS] = {
+    76.6799, 2.46590e-5,
+    3.3946, 2.75e-8,
+    54.8910, 1.38374e-5,
+    0.723330, 0.0,
+    0.006773, -1.302e-9,
+    48.0052, 1.6021302244,
+  },
+  [PLANET_MARS] = {
+    49.5574, 2.11081e-5,
+    1.8497, -1.78e-8,
+    286.5016, 2.92961e-5,
+    1.523688, 0.0,
+    0.093405, 2.516e-9,
+    18.6021, 0.5240207766,
+  },
+  [PLANET_JUPITER] = {
+    100.4542, 2.76854e-5,
+    1.3030, -1.557e-7,
+    273.8777, 1.64505e-5,
+    5.20256, 0.0,
+    0.048498, 4.469e-9,
+    19.8950, 0.0830853001,
+  },
+  [PLANET_SATURN] = {
+    113.6634, 2.38980e-5,
+    2.4886, -1.081e-7,
+    339.3939, 2.97661e-5,
+    9.55475, 0.0,
+    0.055546, -9.499e-9,
+    316.9670, 0.0334442282,
+  },
+  [PLANET_URANUS] = {
+    74.0005, 1.3978e-5,
+    0.7733, 1.9e-8,
+    96.6612, 3.0565e-5,
+    19.18171, -1.55e-8,
+    0.047318, 7.45e-9,
+    142.5905, 0.011725806,
+  },
+  [PLANET_NEPTUNE] = {
+    131.7806, 3.0173e-5,
+    1.7700, -2.55e-7,
+    272.8461, -6.027e-6,
+    30.05826, 3.313e-8,
+    0.008606, 2.15e-9,
+    260.2471, 0.005995147,
+  },
+  /* Approximate fixed elements (not accurate as fixed linear rates over
+   * long spans, but Pluto is only ever used as a slow outer planet with
+   * a wide orb here). */
+  [PLANET_PLUTO] = {
+    110.30347, 0.0,
+    17.14175, 0.0,
+    113.76329, 0.0,
+    39.48168677, 0.0,
+    0.24880766, 0.0,
+    14.53, 0.00396,
+  },
+};
+
+/* Newton-Raphson sqrt. Fixed iteration count rather than a
+ * convergence-check loop, so it can't ever fail to terminate. */
+static double portable_sqrt(double x) {
+  if (x <= 0.0) return 0.0;
+  double guess = (x < 1.0) ? 1.0 : x;
+  for (int i = 0; i < 50; i++) {
+    guess = 0.5 * (guess + x / guess);
+  }
+  return guess;
+}
+
+/* atan(z) for z in [-1,1] via a minimax polynomial (Abramowitz & Stegun
+ * 4.4.49-style coefficients), max error ~1.2e-5 radians. */
+static double portable_atan(double z) {
+  double z2 = z * z;
+  return z * (0.9998660 +
+               z2 * (-0.3302995 +
+                     z2 * (0.1801410 +
+                           z2 * (-0.0851330 + z2 * 0.0208351))));
+}
+
+double lowprec_atan2(double y, double x) {
+  if (x == 0.0 && y == 0.0) return 0.0;
+
+  double ax = x < 0.0 ? -x : x;
+  double ay = y < 0.0 ? -y : y;
+
+  double angle;
+  if (ax >= ay) {
+    angle = portable_atan(ay / ax);
+    if (x < 0.0) angle = DECK_PI - angle;
+  } else {
+    angle = DECK_PI / 2.0 - portable_atan(ax / ay);
+    if (x < 0.0) angle = DECK_PI - angle;
+  }
+  return (y < 0.0) ? -angle : angle;
+}
+
+/* Solves Kepler's equation for the given elements at day-number `d`
+ * (days since J2000.0), returning rectangular heliocentric (geocentric
+ * for the Sun/Moon "orbits") ecliptic coordinates in AU (Earth radii for
+ * the Moon). */
+static void solve_orbit(const OrbitalElements *el, double d, double *x, double *y) {
+  double N = (el->N0 + el->Nd * d) * DEG2RAD;
+  double i = (el->i0 + el->id * d) * DEG2RAD;
+  double w = (el->w0 + el->wd * d) * DEG2RAD;
+  double a = el->a0 + el->ad * d;
+  double e = el->e0 + el->ed * d;
+  double M = normalize_deg(el->M0 + el->Md * d) * DEG2RAD;
+
+  double E = M + e * lowprec_sin(M) * (1.0 + e * lowprec_cos(M));
+  for (int iter = 0; iter < 8; iter++) {
+    double delta = (E - e * lowprec_sin(E) - M) / (1.0 - e * lowprec_cos(E));
+    E -= delta;
+  }
+
+  double xv = a * (lowprec_cos(E) - e);
+  double yv = a * (portable_sqrt(1.0 - e * e) * lowprec_sin(E));
+  double v = lowprec_atan2(yv, xv);
+  double r = portable_sqrt(xv * xv + yv * yv);
+
+  double vw = v + w;
+  *x = r * (lowprec_cos(N) * lowprec_cos(vw) - lowprec_sin(N) * lowprec_sin(vw) * lowprec_cos(i));
+  *y = r * (lowprec_sin(N) * lowprec_cos(vw) + lowprec_cos(N) * lowprec_sin(vw) * lowprec_cos(i));
+}
+
+/* The dozen largest lunar perturbation terms (Schlyter's "more
+ * accurate" Moon correction), applied as a direct correction in degrees
+ * to the Moon's mean-orbit longitude - brings the Moon from several
+ * degrees of error down to a few arcminutes, worth the modest extra
+ * code given how often the Moon matters here (a natal luminary, and the
+ * fastest-moving transiting body). */
+static double moon_longitude_correction(double d) {
+  double Ms = normalize_deg(356.0470 + 0.9856002585 * d) * DEG2RAD;   /* Sun mean anomaly */
+  double Mm = normalize_deg(134.9634 + 13.0649929509 * d) * DEG2RAD;  /* Moon mean anomaly */
+  double Nm = normalize_deg(125.1228 - 0.0529538083 * d) * DEG2RAD;   /* Moon's node */
+  double ws = normalize_deg(282.9404 + 4.70935e-5 * d) * DEG2RAD;     /* Sun's perihelion */
+  double wm = normalize_deg(318.0634 + 0.1643573223 * d) * DEG2RAD;   /* Moon's perihelion */
+
+  double Ls = ws + Ms;   /* Sun's mean longitude */
+  double Lm = Nm + wm + Mm; /* Moon's mean longitude */
+  double D = Lm - Ls;    /* elongation */
+  double F = Lm - Nm;    /* argument of latitude */
+
+  double corr = 0.0;
+  corr += -1.274 * lowprec_sin(Mm - 2.0 * D);
+  corr += 0.658 * lowprec_sin(2.0 * D);
+  corr += -0.186 * lowprec_sin(Ms);
+  corr += -0.059 * lowprec_sin(2.0 * Mm - 2.0 * D);
+  corr += -0.057 * lowprec_sin(Mm - 2.0 * D + Ms);
+  corr += 0.053 * lowprec_sin(Mm + 2.0 * D);
+  corr += 0.046 * lowprec_sin(2.0 * D - Ms);
+  corr += 0.041 * lowprec_sin(Mm - Ms);
+  corr += -0.035 * lowprec_sin(D);
+  corr += -0.031 * lowprec_sin(Mm + Ms);
+  corr += -0.015 * lowprec_sin(2.0 * F - 2.0 * D);
+  corr += 0.011 * lowprec_sin(Mm - 4.0 * D);
+  return corr;
+}
+
+double lowprec_geocentric_longitude(Body body, double jd) {
+  double d = jd - 2451545.0;
+
+  double sun_x, sun_y;
+  solve_orbit(&k_elements[PLANET_SUN], d, &sun_x, &sun_y);
+
+  if (body == PLANET_SUN) {
+    return normalize_deg(lowprec_atan2(sun_y, sun_x) * RAD2DEG);
+  }
+
+  if (body == PLANET_MOON) {
+    double moon_x, moon_y;
+    solve_orbit(&k_elements[PLANET_MOON], d, &moon_x, &moon_y);
+    double lon = lowprec_atan2(moon_y, moon_x) * RAD2DEG + moon_longitude_correction(d);
+    return normalize_deg(lon);
+  }
+
+  /* Earth's own heliocentric position is the Sun's geocentric one,
+   * negated (both "orbits" share the same ecliptic plane by
+   * definition here, i.e. i = 0 for the Sun's elements). */
+  double earth_x = -sun_x, earth_y = -sun_y;
+
+  double planet_x, planet_y;
+  solve_orbit(&k_elements[body], d, &planet_x, &planet_y);
+
+  double geo_x = planet_x - earth_x;
+  double geo_y = planet_y - earth_y;
+  return normalize_deg(lowprec_atan2(geo_y, geo_x) * RAD2DEG);
+}