1 files changed, 101 insertions, 0 deletions
diff --git a/stm32/app/events.c b/stm32/app/events.c
new file mode 100644
index 0000000..2ad9b47
--- /dev/null
+++ b/stm32/app/events.c
@@ -0,0 +1,101 @@
+#include "project.h"
+
+#define D2R(a) ((a)*(M_PI/180.))
+
+
+
+
+// shamelessly stolen from Meeus Astronmical Algorithms Chapter 27 table 27.B
+// chapter 25 is probably more tractable, but this is easier to code up
+
+static double mean_vernal_equinox (unsigned year)
+{
+  double y = year;
+  y -= 2000;
+  y *= 0.001;
+  return 2451623.80984 + 365242.37404 * y + 0.05169 * (y * y) - 0.00411 * (y * y * y) - 0.00057 * (y * y * y * y);
+}
+
+static double mean_summer_solstice (unsigned year)
+{
+  double y = year;
+  y -= 2000;
+  y *= 0.001;
+  return 2451716.56767 + 365241.62603 * y + 0.00325 * (y * y) + 0.00888 * (y * y * y) - 0.00030 * (y * y * y * y);
+}
+
+static double mean_autumnal_equinox (unsigned year)
+{
+  double y = year;
+  y -= 2000;
+  y *= 0.001;
+  return 2451810.21715 + 365242.01767 * y - 0.11575 * (y * y) + 0.00337 * (y * y * y) + 0.00078 * (y * y * y * y);
+}
+
+static double mean_winter_solstice (unsigned year)
+{
+  double y = year;
+  y -= 2000;
+  y *= 0.001;
+  return 2451900.05952 + 365242.74049 * y - 0.06223 * (y * y) - 0.00823 * (y * y * y) + 0.00032 * (y * y * y * y);
+}
+
+
+static double orbital_periodic_terms (double t)
+{
+#define N 24
+  const double A[N] = {485, 203, 199, 182, 156, 136, 77, 74, 70, 58, 52, 50, 45,
+                       44, 29, 18, 17, 16, 14, 12, 12, 12, 9, 8
+                      };
+  const double B[N] = {D2R (324.96), D2R (337.23), D2R (342.08), D2R (27.85),
+                       D2R (73.14), D2R (171.52), D2R (222.54), D2R (296.72),
+                       D2R (243.58), D2R (119.81), D2R (297.17), D2R (21.02),
+                       D2R (247.54), D2R (325.15), D2R (60.93), D2R (155.12),
+                       D2R (288.79), D2R (198.04), D2R (199.76), D2R (95.39),
+                       D2R (287.11), D2R (320.81), D2R (227.73), D2R (15.45)
+                      };
+  const double C[N] = {D2R (1934.136), D2R (32964.467), D2R (20.186),
+                       D2R (445267.112), D2R (45036.886), D2R (22518.443),
+                       D2R (65928.934), D2R (3034.906), D2R (9037.513),
+                       D2R (33718.147), D2R (150.678), D2R (2281.226),
+                       D2R (29929.562), D2R (31555.956), D2R (4443.417),
+                       D2R (67555.328), D2R (4562.452), D2R (62894.029),
+                       D2R (31436.921), D2R (14577.848), D2R (31931.756),
+                       D2R (34777.259), D2R (1222.114), D2R (16859.074)
+                      };
+  double s = 0;
+  unsigned i;
+
+  for (i = 0; i < N; ++i)
+    s += A[i] * cos (B[i] + (C[i] * t));
+
+  return s;
+}
+
+
+static double mean_to_real (double j0)
+{
+
+  double t = (j0 - 2451545.) / 36525.;
+  double w = D2R ((35999.373 * t) - 2.47);
+  double dl = 1 + 0.0334 * cos (w) + 0.0007 * cos (2. * w);
+  double s = orbital_periodic_terms (t);
+
+#if 0
+  printf ("j0=%.6f\r\n", j0);
+  printf ("t=%.6f\r\n", t);
+  printf ("w=%.6f\r\n", w);
+  printf ("dl=%.6f\r\n", dl);
+  printf ("s=%.6f\r\n", s);
+#endif
+
+  return j0 + ((0.00001 * s) / dl);
+}
+
+
+
+double autumnal_equinox (unsigned y)
+{
+  return mean_to_real (mean_autumnal_equinox (y));
+  //  return mean_to_real (mean_summer_solstice (y));
+}