deck_in_a_dash/engine/src/astro.c

#define _POSIX_C_SOURCE 200809L /* for gmtime_r */

#include "astro.h"
#include "i18n.h"
#include "../third_party/astronomy/astronomy.h"

#include <math.h>
#include <stdbool.h>
#include <string.h>

/* DEG2RAD / RAD2DEG come from astronomy.h. */

static const astro_body_t k_astro_body[NUM_BODIES] = {
  [PLANET_SUN] = BODY_SUN,
  [PLANET_MOON] = BODY_MOON,
  [PLANET_MERCURY] = BODY_MERCURY,
  [PLANET_VENUS] = BODY_VENUS,
  [PLANET_MARS] = BODY_MARS,
  [PLANET_JUPITER] = BODY_JUPITER,
  [PLANET_SATURN] = BODY_SATURN,
  [PLANET_URANUS] = BODY_URANUS,
  [PLANET_NEPTUNE] = BODY_NEPTUNE,
  [PLANET_PLUTO] = BODY_PLUTO,
};

static double normalize_degrees(double deg) {
  double d = fmod(deg, 360.0);
  return d < 0.0 ? d + 360.0 : d;
}

static void longitude_to_position(double longitude, PlanetPosition *out) {
  out->ecliptic_longitude = normalize_degrees(longitude);
  out->sign = (ZodiacSign)((int)(out->ecliptic_longitude / 30.0) % 12);
  out->degree_in_sign = fmod(out->ecliptic_longitude, 30.0);
}

/* Mean obliquity of the ecliptic (IAU low-precision polynomial), matching
 * what Astronomy Engine computes internally in its (non-exported)
 * mean_obliq(). T is Julian centuries since J2000.0. */
static double mean_obliquity_deg(double t_centuries) {
  double t = t_centuries;
  return 23.4392911 - 0.0130042 * t - 0.00000016 * t * t + 0.000000504 * t * t * t;
}

/* Ascendant (rising ecliptic degree), via the standard RAMC/obliquity/
 * latitude identity (Duffett-Smith & Zwart, "Practical Astronomy with
 * your Calculator or Spreadsheet"). Astronomy Engine has no built-in
 * Ascendant function. */
static double compute_ascendant(astro_time_t *time, double latitude_deg,
                                 double longitude_deg) {
  double gast_hours = Astronomy_SiderealTime(time);
  double ramc_deg = normalize_degrees(gast_hours * 15.0 + longitude_deg);
  double eps_deg = mean_obliquity_deg(time->tt / 36525.0);

  double ramc = ramc_deg * DEG2RAD;
  double eps = eps_deg * DEG2RAD;
  double lat = latitude_deg * DEG2RAD;

  double y = -cos(ramc);
  double x = sin(ramc) * cos(eps) + tan(lat) * sin(eps);
  return normalize_degrees(atan2(y, x) * RAD2DEG);
}

static astro_time_t time_from_birth(const BirthData *birth) {
  astro_time_t local = Astronomy_MakeTime(birth->year, birth->month, birth->day,
                                           birth->hour, birth->minute, 0.0);
  return Astronomy_AddDays(local, -birth->utc_offset_hours / 24.0);
}

static astro_time_t time_from_unix(time_t utc_moment) {
  struct tm tm_utc;
  gmtime_r(&utc_moment, &tm_utc);
  return Astronomy_MakeTime(tm_utc.tm_year + 1900, tm_utc.tm_mon + 1, tm_utc.tm_mday,
                             tm_utc.tm_hour, tm_utc.tm_min, (double)tm_utc.tm_sec);
}

static void compute_body_positions(astro_time_t time, PlanetPosition out[NUM_BODIES]) {
  /* Astronomy_EclipticLongitude() computes *heliocentric* longitude (and
   * outright rejects BODY_SUN) - wrong for astrology, which needs the
   * apparent geocentric position. Astronomy_GeoVector() + Astronomy_Ecliptic()
   * gives that for every body, Sun included. */
  for (int b = 0; b < NUM_BODIES; b++) {
    astro_vector_t geo = Astronomy_GeoVector(k_astro_body[b], time, ABERRATION);
    astro_ecliptic_t eclip = Astronomy_Ecliptic(geo);
    longitude_to_position(eclip.elon, &out[b]);
  }
}

void astro_compute_natal_chart(const BirthData *birth, NatalChart *out) {
  astro_time_t time = time_from_birth(birth);

  compute_body_positions(time, out->bodies);
  out->ascendant_longitude = compute_ascendant(&time, birth->latitude, birth->longitude);

  int ascendant_sign = (int)(out->ascendant_longitude / 30.0) % 12;
  for (int house = 0; house < 12; house++) {
    out->houses[house] = (ZodiacSign)((ascendant_sign + house) % 12);
  }
}

static MoonPhaseName moon_phase_from_angle(double angle_deg) {
  int bucket = ((int)((angle_deg + 22.5) / 45.0)) % 8;
  return (MoonPhaseName)bucket;
}

static const double k_aspect_angle[5] = {
  [ASPECT_CONJUNCTION] = 0.0,
  [ASPECT_SEXTILE] = 60.0,
  [ASPECT_SQUARE] = 90.0,
  [ASPECT_TRINE] = 120.0,
  [ASPECT_OPPOSITION] = 180.0,
};

static double angular_separation(double a, double b) {
  double diff = fabs(normalize_degrees(a) - normalize_degrees(b));
  return diff > 180.0 ? 360.0 - diff : diff;
}

static double orb_for_pair(Body transiting, Body natal) {
  /* Wider orb when a luminary (Sun or Moon) is involved, per standard
   * convention. */
  bool luminary = transiting == PLANET_SUN || transiting == PLANET_MOON ||
                  natal == PLANET_SUN || natal == PLANET_MOON;
  return luminary ? 8.0 : 6.0;
}

void astro_compute_daily_transits(time_t utc_moment, const NatalChart *natal,
                                   DailyTransits *out) {
  astro_time_t time = time_from_unix(utc_moment);

  compute_body_positions(time, out->bodies);

  astro_angle_result_t phase = Astronomy_MoonPhase(time);
  out->moon_phase = moon_phase_from_angle(phase.angle);

  out->aspect_count = 0;
  for (int t = 0; t < NUM_BODIES; t++) {
    for (int n = 0; n < NUM_BODIES; n++) {
      double separation = angular_separation(out->bodies[t].ecliptic_longitude,
                                              natal->bodies[n].ecliptic_longitude);
      double orb_limit = orb_for_pair((Body)t, (Body)n);

      for (int a = 0; a < 5; a++) {
        double orb = fabs(separation - k_aspect_angle[a]);
        if (orb <= orb_limit && out->aspect_count < MAX_ASPECTS) {
          Aspect *aspect = &out->aspects[out->aspect_count++];
          aspect->transiting_planet = (Body)t;
          aspect->natal_planet = (Body)n;
          aspect->type = (AspectType)a;
          aspect->orb = orb;
        }
      }
    }
  }
}

static const char *const k_body_names[NUM_BODIES] = {
  "Sun", "Moon", "Mercury", "Venus", "Mars",
  "Jupiter", "Saturn", "Uranus", "Neptune", "Pluto",
};

static const char *const k_sign_names[12] = {
  "Aries", "Taurus", "Gemini", "Cancer", "Leo", "Virgo",
  "Libra", "Scorpio", "Sagittarius", "Capricorn", "Aquarius", "Pisces",
};

static const char *const k_moon_phase_names[8] = {
  "New Moon", "Waxing Crescent", "First Quarter", "Waxing Gibbous",
  "Full Moon", "Waning Gibbous", "Last Quarter", "Waning Crescent",
};

static const char *const k_aspect_names[5] = {
  "Conjunction", "Sextile", "Square", "Trine", "Opposition",
};

/* Slugs used to build i18n.c lookup keys - kept separate from the C enum
 * names so a translation file's keys don't depend on identifiers that
 * might get renamed. */
static const char *const k_body_slug[NUM_BODIES] = {
  "sun", "moon", "mercury", "venus", "mars",
  "jupiter", "saturn", "uranus", "neptune", "pluto",
};

static const char *const k_sign_slug[12] = {
  "aries", "taurus", "gemini", "cancer", "leo", "virgo",
  "libra", "scorpio", "sagittarius", "capricorn", "aquarius", "pisces",
};

static const char *const k_moon_phase_slug[8] = {
  "new", "waxing_crescent", "first_quarter", "waxing_gibbous",
  "full", "waning_gibbous", "last_quarter", "waning_crescent",
};

static const char *const k_aspect_slug[5] = {
  "conjunction", "sextile", "square", "trine", "opposition",
};

static const char *lookup(const char *prefix, const char *slug, const char *fallback) {
  char key[64];
  strcpy(key, prefix);
  strcat(key, slug);
  return i18n_get(key, fallback);
}

const char *astro_body_name(Body body) {
  return lookup("body.", k_body_slug[body], k_body_names[body]);
}
const char *astro_sign_name(ZodiacSign sign) {
  return lookup("sign.", k_sign_slug[sign], k_sign_names[sign]);
}
const char *astro_moon_phase_name(MoonPhaseName phase) {
  return lookup("moonphase.", k_moon_phase_slug[phase], k_moon_phase_names[phase]);
}
const char *astro_aspect_name(AspectType type) {
  return lookup("aspect.", k_aspect_slug[type], k_aspect_names[type]);
}