//                                meridian orrery
//                                                     by illusionmanager 2026
//
// this version should work with or without the upgrade with a touch sensor
#include <Wire.h>
#include <RTClib.h>
#include <math.h>
#include <stdio.h>
#include <string.h>

// --- hardware + motion constants ---
// MICRO_STEPS is fixed at 8
#define MICRO_STEPS 8

// STEPS_PER_DEGREE is a measured value for motor labeled 1:236
// which turns out to be about 1:236.14               
// pick the one that you have.
// #define STEPS_PER_DEGREE 13.11944 // motor ratio labeled 1:236 = 236.14
#define STEPS_PER_DEGREE 33.52681 // motor ratio labeled 1:603 = 603.48

// PULSE_TIME determines the speed at which the motor runs
// A lower value makes it run faster, but below some value,
// the motor wouldn't run at all, or not reliable.
#define PULSE_TIME 58

// interchange these two if the orrery doesn't run clockwise during reset
// this depends on how you wired the motor.
#define CLOCKWISE HIGH
#define COUNTERCLOCKWISE LOW

// use the ZODIAC_OFFSET (0-359) if Neptune isn't exactly at Aries/Pisces boundary after reset.
#define ZODIAC_OFFSET 2 

#define DIR_PIN GPIO_NUM_1
#define STEP_PIN GPIO_NUM_2
#define ENABLE_PIN GPIO_NUM_3
#define TOUCH_PIN GPIO_NUM_4 // for the upgraded version

#define SDA_PIN GPIO_NUM_8
#define SCL_PIN GPIO_NUM_9

#define MAGNET_PIN GPIO_NUM_10

#define EXTRA_HOME 40 

// These values depend on printing accuracy, you shouldn't need to change them.
// But if you want the way to do it is to connect the device to the computer
// send the 'r' command and then enter 2490 followed by a few 1's until Neptune 
// moves a tiny bit, add those numbers together and subtract 1. That is the value
// for Neptune. Next enter -2130 followed by a few -1's until Uranus starts to move, 
// continue until you have done all the planets. (you could go in smaller steps like 0.5 degree)
static const double engage[8] = {
  0.0,    // Mercury
  356.0,  // Venus
  710.5,  // Earth
  1066.0, // Mars
  1423.0, // Jupiter
  1780.5, // Saturn
  2139.0, // Uranus
  2497.0  // Neptune
};

RTC_DS3231 rtc;

static int prev_day = 0;

struct Poly4 {
  double a0, a1, a2, a3;
};

struct MeeusElemDate {
  struct Poly4 L;
  struct Poly4 a_au;
  struct Poly4 e;
  struct Poly4 i;
  struct Poly4 omega;
  struct Poly4 pi;
};

static const char *planet_names[8] = {
  "Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune"
};

static float start_offset[8];

static const struct MeeusElemDate meeus_elem[8] = {
  // Mercury
  { { 252.250906, 149474.0722491, 0.00030350, 0.000000018 },
    { 0.387098310, 0.0, 0.0, 0.0 },
    { 0.20563175, 0.000020407, -0.0000000283, -0.00000000018 },
    { 7.004986, 0.0018215, -0.00001810, 0.000000056 },
    { 48.330893, 1.1861883, 0.00017542, 0.000000215 },
    { 77.456119, 1.5564776, 0.00029544, 0.000000009 } },

  // Venus
  { { 181.979801, 58519.2130302, 0.00031014, 0.000000015 },
    { 0.723329820, 0.0, 0.0, 0.0 },
    { 0.00677192, -0.000047765, 0.0000000981, 0.00000000046 },
    { 3.394662, 0.0010037, -0.00000088, -0.000000007 },
    { 76.679920, 0.9011206, 0.00040618, -0.000000093 },
    { 131.563703, 1.4022288, -0.00107618, -0.000005678 } },

  // Earth
  { { 100.466457, 36000.7698278, 0.00030322, 0.000000020 },
    { 1.000001018, 0.0, 0.0, 0.0 },
    { 0.01670863, -0.000042037, -0.0000001267, 0.00000000014 },
    { 0.0, 0.0, 0.0, 0.0 },
    { 0.0, 0.0, 0.0, 0.0 },
    { 102.937348, 1.7195366, 0.00045688, -0.000000018 } },

  // Mars
  { { 355.433000, 19141.6964471, 0.00031052, 0.000000016 },
    { 1.523679342, 0.0, 0.0, 0.0 },
    { 0.09340065, 0.000090484, -0.0000000806, -0.00000000025 },
    { 1.849726, -0.0006011, 0.00001276, -0.000000007 },
    { 49.558093, 0.7720959, 0.00001557, 0.000002267 },
    { 336.060234, 1.8410449, 0.00013477, 0.000000536 } },

  // Jupiter
  { { 34.351519, 3036.3027748, 0.00022330, 0.000000037 },
    { 5.202603209, 0.0000001913, 0.0, 0.0 },
    { 0.04849793, 0.000163225, -0.0000004714, -0.00000000201 },
    { 1.303267, -0.0054965, 0.00000466, -0.000000002 },
    { 100.464407, 1.0209774, 0.00040315, 0.000000404 },
    { 14.331207, 1.6126352, 0.00103042, -0.000004464 } },

  // Saturn
  { { 50.077444, 1223.5110686, 0.00051908, -0.000000030 },
    { 9.554909192, -0.0000021390, 0.000000004, 0.0 },
    { 0.05554814, -0.000346641, -0.0000006436, 0.00000000340 },
    { 2.488879, -0.0037362, -0.00001519, 0.000000087 },
    { 113.665503, 0.8770880, -0.00012176, -0.000002249 },
    { 93.057237, 1.9637613, 0.00083753, 0.000004928 } },

  // Uranus
  { { 314.055005, 429.8640561, 0.00030390, 0.000000026 },
    { 19.218446062, -0.0000000372, 0.00000000098, 0.0 },
    { 0.04638122, -0.000027293, 0.0000000789, 0.00000000024 },
    { 0.773197, 0.0007744, 0.00003749, -0.000000092 },
    { 74.005957, 0.5211278, 0.00133947, 0.000018484 },
    { 173.005291, 1.4863790, 0.00021406, 0.000000434 } },

  // Neptune
  { { 304.348665, 219.8833092, 0.00030882, 0.000000018 },
    { 30.110386869, -0.0000001663, 0.00000000069, 0.0 },
    { 0.00945575, 0.000006033, 0.0, -0.00000000005 },
    { 1.769953, -0.0093082, -0.00000708, 0.000000027 },
    { 131.784057, 1.1022039, 0.00025952, -0.000000637 },
    { 48.120276, 1.4262957, 0.00038434, 0.000000020 } }
};

volatile bool magnet_seen = false;

char cmd_buf[64];
int cmd_pos = 0;

static double deg2rad(double x) {
  return x * (PI / 180.0);
}

static double rad2deg(double x) {
  return x * (180.0 / PI);
}

static double wrap360(double x) {
  while (x < 0.0) x += 360.0;
  while (x >= 360.0) x -= 360.0;
  return x;
}

static double poly4_eval(const struct Poly4 &p, double T) {
  return ((p.a3 * T + p.a2) * T + p.a1) * T + p.a0;
}

static double solve_kepler_E(double M, double e) {
  double E = M + e * sin(M) * (1.0 + e * cos(M));
  for (int i = 0; i < 10; i++) {
    double f = E - e * sin(E) - M;
    double fp = 1.0 - e * cos(E);
    double dE = f / fp;
    E -= dE;
    if (fabs(dE) < 1e-12) break;
  }
  return E;
}

static double heliocentric_longitude(double N, double inc, double w, double a_au, double e, double M) {
  N = deg2rad(wrap360(N));
  inc = deg2rad(inc);
  w = deg2rad(wrap360(w));
  M = deg2rad(wrap360(M));

  double Ec = solve_kepler_E(M, e);

  double xv = a_au * (cos(Ec) - e);
  double yv = a_au * (sqrt(1.0 - e * e) * sin(Ec));

  double v = atan2(yv, xv);
  double r = sqrt(xv * xv + yv * yv);

  double xh = r * (cos(N) * cos(v + w) - sin(N) * sin(v + w) * cos(inc));
  double yh = r * (sin(N) * cos(v + w) + cos(N) * sin(v + w) * cos(inc));

  return wrap360(rad2deg(atan2(yh, xh)));
}

static double julian_day_from_datetime(const DateTime &dt) {
  int y = dt.year();
  int m = dt.month();
  int d = dt.day();

  if (m <= 2) {
    y -= 1;
    m += 12;
  }

  int A = y / 100;
  int B = 2 - A + (A / 4);

  double day_fraction =
      (dt.hour() +
      (dt.minute() +
      dt.second() / 60.0) / 60.0) / 24.0;

  double jd =
      floor(365.25 * (y + 4716))
    + floor(30.6001 * (m + 1))
    + d + day_fraction + B - 1524.5;

  return jd;
}

static double planet_helio_lon(double jd_ut, int p) {
  double T = (jd_ut - 2451545.0) / 36525.0; // centuries
  const struct MeeusElemDate &E = meeus_elem[p];

  double L = wrap360(poly4_eval(E.L, T));
  double a = poly4_eval(E.a_au, T);
  double ecc = poly4_eval(E.e, T);
  double inc = poly4_eval(E.i, T);
  double Om = wrap360(poly4_eval(E.omega, T));
  double piL = wrap360(poly4_eval(E.pi, T));

  double M = wrap360(L - piL);
  double w = wrap360(piL - Om);

  return wrap360(heliocentric_longitude(Om, inc, w, a, ecc, M));
}

static double moon_phase_angle(double jd) {
  double T = (jd - 2451545.0) / 36525.0;

  double L0 = 218.3164477 + 481267.88123421 * T;
  double M  = 357.5291092 + 35999.0502909 * T;
  double Mp = 134.9633964 + 477198.8675055 * T;
  double D  = 297.8501921 + 445267.1114034 * T;
  double F  = 93.2720950 + 483202.0175233 * T;

  double Mr  = deg2rad(M);
  double Mpr = deg2rad(Mp);
  double Dr  = deg2rad(D);
  double Fr  = deg2rad(F);

  double Ls =
      280.46646 + 36000.76983 * T
    + 1.914602 * sin(Mr)
    + 0.019993 * sin(2.0 * Mr);

  double Lm =
      L0
    + 6.288774 * sin(Mpr)
    + 1.274027 * sin(2.0 * Dr - Mpr)
    + 0.658314 * sin(2.0 * Dr)
    + 0.213618 * sin(2.0 * Mpr)
    - 0.185116 * sin(Mr)
    - 0.114332 * sin(2.0 * Fr);

  return wrap360(Lm - Ls);
}

static void print_angles() {
  DateTime now = rtc.now();
  double jd = julian_day_from_datetime(now);

  Serial.printf("Angles for %02d-%02d-%04d %02d:%02d:%02d\n",
    now.day(), now.month(), now.year(),
    now.hour(), now.minute(), now.second());

  for (int p = 0; p < 8; p++) {
    double lon = planet_helio_lon(jd, p);
    Serial.printf("%-8s  %8.3f deg\n", planet_names[p], lon);
  }

  Serial.printf("%-8s  %8.3f deg\n", "Moon", moon_phase_angle(jd));
}

// --- motor helpers ---

void enable_motor() {
  digitalWrite(ENABLE_PIN, LOW);
}

void disable_motor() {
  digitalWrite(ENABLE_PIN, HIGH);
}

// the rotate use complex code to decelerate nicely
void rotate(float deg) {
  if (deg < 0) digitalWrite(DIR_PIN, CLOCKWISE);
  else         digitalWrite(DIR_PIN, COUNTERCLOCKWISE);

  long steps_to_do = lround(fabs(deg) * STEPS_PER_DEGREE * MICRO_STEPS);

 
  const float decel_degrees = 30.0;
  const float decel_factor = 3.0;

  const long decel_steps = decel_degrees * STEPS_PER_DEGREE * MICRO_STEPS;
  while (steps_to_do > 0) {

    int this_pulse_time = PULSE_TIME;

    if (steps_to_do < decel_steps && decel_steps > 1) {
      float t = 1.0f - (float)(steps_to_do - 1) / (float)(decel_steps - 1);

      // S-curve easing from 0 to 1
      float u = t * t * (3.0f - 2.0f * t);

      // speed factor goes from 1.0 down to 1.0 / decel_factor
      float min_speed_factor = 1.0f / decel_factor;
      float speed_factor = 1.0f - (1.0f - min_speed_factor) * u;

      // pulse time is inverse of speed
      this_pulse_time = (int)lroundf((float)PULSE_TIME / speed_factor);
    }

    digitalWrite(STEP_PIN, HIGH);
    delayMicroseconds(this_pulse_time);
    digitalWrite(STEP_PIN, LOW);
    delayMicroseconds(this_pulse_time);

    steps_to_do--;
  }
}
// --- homing ---

void IRAM_ATTR magnet_isr() {
  magnet_seen = true;
}

void go_home(float extra_after_home) {
  Serial.println("going home");
  enable_motor();

  if (digitalRead(MAGNET_PIN)) {
    Serial.println("magnet active, clearing it first");

    digitalWrite(DIR_PIN, CLOCKWISE);

    long quiet_steps = lround(20.0 * STEPS_PER_DEGREE * MICRO_STEPS);
    long steps_since_seen = 0;

    while (steps_since_seen < quiet_steps) {
      digitalWrite(STEP_PIN, HIGH);
      delayMicroseconds(PULSE_TIME);
      digitalWrite(STEP_PIN, LOW);
      delayMicroseconds(PULSE_TIME);

      if (digitalRead(MAGNET_PIN)) {
        steps_since_seen = 0;
      } else {
        steps_since_seen++;
      }
    }
  }

  magnet_seen = false;
  attachInterrupt(digitalPinToInterrupt(MAGNET_PIN), magnet_isr, RISING);

  digitalWrite(DIR_PIN, CLOCKWISE);

  for (long i = 0; i < (long)(STEPS_PER_DEGREE * MICRO_STEPS * 360.0 * 6.0); i++) {
    digitalWrite(STEP_PIN, HIGH);
    delayMicroseconds(PULSE_TIME);
    digitalWrite(STEP_PIN, LOW);
    delayMicroseconds(PULSE_TIME);

    if (magnet_seen) {
      Serial.println("zero detected");
      break;
    }
  }

  detachInterrupt(digitalPinToInterrupt(MAGNET_PIN));

  rotate(-(360.0 * 3.0) - EXTRA_HOME - extra_after_home);

  disable_motor();
  Serial.println("home done");
}
// --- derived geometry ---

void compute_start_offsets() {
  float tab_loss[8];

  for (int i = 0; i < 8; i++) {
    float turns = roundf((float)engage[i] / 360.0f);
    tab_loss[i] = turns * 360.0f - (float)engage[i];
  }

  start_offset[7] = 0.0f;
  start_offset[6] = tab_loss[6] - tab_loss[7];
  start_offset[5] = 0.0f;
  start_offset[4] = tab_loss[4] - tab_loss[5];
  start_offset[3] = 0.0f;
  start_offset[2] = tab_loss[2] - tab_loss[3];
  start_offset[1] = 0.0f;
  start_offset[0] = tab_loss[0] - tab_loss[1];

}

// --- placement ---

static void place_targets(const float target[8]) {
 
  enable_motor();

  float prev_target = -ZODIAC_OFFSET;
  int dir = 1;

  for (int i = 7; i >= 0; i--) {
    delay(250);
    float current_angle = wrap360(prev_target + start_offset[i]);
    float extra = dir * wrap360(dir * (target[i] - current_angle));
    float move = dir * engage[i] + extra;

    Serial.printf("Rotating to %.1f deg for %s\n", target[i], planet_names[i]);

    rotate(move);

    prev_target = target[i];
    dir = -dir;
  }
  disable_motor();
}

static void do_solar_system() {
  float target[8];
  DateTime now = rtc.now();
  print_date_time();
  double jd = julian_day_from_datetime(now);

  for (int p = 0; p < 8; p++) {
    target[p] = (float)planet_helio_lon(jd, p);
  }

  float moon_phase = (float)moon_phase_angle(jd);
  float moon_correction = wrap360(moon_phase + 11.0f * target[2] - ZODIAC_OFFSET) / 11.0f;

  go_home(moon_correction);
  Serial.printf("Moon phase %.1f\n", moon_phase);
  place_targets(target);
  Serial.println("Solar System done.");
  print_date_time();
}

// --- serial commands ---

void print_help() {
  Serial.println("Commands:");
  Serial.println("  h or ?                  help");
  Serial.println("  d                       print date and time");
  Serial.println("  dDD-MM-YYYY             set date");
  Serial.println("  dHH:MM                  set time");
  Serial.println("  dDD-MM-YYYY HH:MM       set date and time");
  Serial.println("  <number>                rotate motor in degrees for example 5 or -30");
  Serial.println("  r                       reset / go home");
  Serial.println("  p                       print planet angles");
  Serial.println("  s                       set all planets and moon");
}

void print_date_time() {
  DateTime now = rtc.now();
  Serial.printf("%02d-%02d-%04d %02d:%02d:%02d\n",
    now.day(), now.month(), now.year(),
    now.hour(), now.minute(), now.second());
}

void handle_date_command(const char *s) {
  int day, month, year, hour, minute;
  DateTime now = rtc.now();

  if (strcmp(s, "d") == 0) {
    print_date_time();
    return;
  }

  if (sscanf(s, "d%2d-%2d-%4d %2d:%2d", &day, &month, &year, &hour, &minute) == 5) {
    rtc.adjust(DateTime(year, month, day, hour, minute, 0));
    Serial.println("Date and time set");
    print_date_time();
    return;
  }

  if (sscanf(s, "d%2d-%2d-%4d", &day, &month, &year) == 3) {
    rtc.adjust(DateTime(year, month, day, now.hour(), now.minute(), now.second()));
    Serial.println("Date set");
    print_date_time();
    return;
  }

  if (sscanf(s, "d%2d:%2d", &hour, &minute) == 2) {
    rtc.adjust(DateTime(now.year(), now.month(), now.day(), hour, minute, 0));
    Serial.println("Time set");
    print_date_time();
    return;
  }

  Serial.println("Bad date/time format. Use h or ?");
}

void handle_rotate_command(const char *s) {
  float deg;
  char extra;

  if (sscanf(s, "%f %c", &deg, &extra) == 1) {
    enable_motor();
    rotate(deg);
    disable_motor();
    Serial.printf("Rotated %.1f degrees\n", deg);
    Serial.flush();
  } else {
    Serial.println("Bad rotation value");
  }
}

void handle_command(const char *s) {
  if (s[0] == 0) return;

  if (strcmp(s, "h") == 0 || strcmp(s, "?") == 0) {
    print_help();
    return;
  }

  if (strcmp(s, "r") == 0) {
    go_home(0);
    return;
  }

  if (strcmp(s, "p") == 0) {
    print_angles();
    return;
  }

  if (strcmp(s, "s") == 0) {
    do_solar_system();
    return;
  }

  if (s[0] == 'd') {
    handle_date_command(s);
    return;
  }

  if ((s[0] >= '0' && s[0] <= '9') || s[0] == '-' || s[0] == '+') {
    handle_rotate_command(s);
    return;
  }

  Serial.println("Unknown command. Use h or ?");
}

void setup() {
  Serial.begin(115200);
  delay(1500);
  Serial.println("Meridian Orrery, by illusionmanager (2026)");
  
  Wire.begin(SDA_PIN, SCL_PIN);

  pinMode(STEP_PIN, OUTPUT);
  pinMode(DIR_PIN, OUTPUT);
  pinMode(ENABLE_PIN, OUTPUT);
  pinMode(MAGNET_PIN, INPUT_PULLDOWN);
  pinMode(TOUCH_PIN, INPUT_PULLDOWN); // for the upgraded version

  disable_motor();

  if (!rtc.begin()) {
    Serial.println("Couldn't find RTC");
    Serial.flush();
    while (1) delay(10);
  }

  if (rtc.lostPower()) {
    Serial.println("RTC lost power, setting compile time");
    rtc.adjust(DateTime(F(__DATE__), F(__TIME__)));
  }

  compute_start_offsets();


  Serial.println("Press h for help");
  DateTime now = rtc.now();
  prev_day = now.day();
  do_solar_system();
}

void loop() {
  static uint32_t touch_start_time = 0;
  static bool touch_active = false;
  
  DateTime now = rtc.now();
  // run the solar system once a day at this hour
  if (now.hour() == 12 && now.day() != prev_day) {
    do_solar_system();
    prev_day = now.day();
  }
  if (Serial.available()) {
    char c = Serial.read();

    if (c == '\r' || c == '\n') {
      if (cmd_pos > 0) {
        cmd_buf[cmd_pos] = 0;
        handle_command(cmd_buf);
        cmd_pos = 0;
      }
    } else {
      if (cmd_pos < (int)sizeof(cmd_buf) - 1) {
        cmd_buf[cmd_pos++] = c;
      }
    }
  }
  // next part is for when you installed the upgrade
  bool pin_is_touched = (digitalRead(TOUCH_PIN) == HIGH);

  if (pin_is_touched && !touch_active) {
    touch_start_time = millis();
    touch_active = true;
  }

  if (!pin_is_touched && touch_active) {
    uint32_t touch_duration = millis() - touch_start_time;

    if (touch_duration < 500) {
      do_solar_system();
      prev_day = now.day();
      delay(1000);
    }
    touch_active = false;
  }
}