1 files changed, 173 insertions, 0 deletions
diff --git a/st/siderial_time.c b/st/siderial_time.c
new file mode 100644
index 0000000..5175511
--- /dev/null
+++ b/st/siderial_time.c
@@ -0,0 +1,173 @@
+#include <stdio.h>
+#include <math.h>
+
+typedef struct
+{
+  unsigned year;
+  unsigned is_leap;
+  unsigned jday;
+  unsigned month;
+  unsigned mday;
+  unsigned wday;
+  unsigned hour;
+  unsigned minute;
+  unsigned second;
+  unsigned nanosecond;
+} UTC;
+
+typedef struct
+{
+  unsigned hour;
+  unsigned minute;
+  unsigned second;
+  unsigned nanosecond;
+} ST;
+
+static int
+is_leap (unsigned year)
+{
+  if (year % 4)
+    return 0;
+
+  if (year % 100)
+    return 1;
+
+  if (year % 400)
+    return 0;
+
+  return 1;
+}
+
+
+double
+time_utc_to_tjd (UTC u)
+{
+  unsigned y400;
+  unsigned y100;
+  unsigned y4;
+
+  double ret, fd;
+  unsigned jd;
+
+  static int const mdays[] =
+    { 0, 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334 };
+  static int const lmdays[] =
+    { 0, 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335 };
+
+  u.is_leap = is_leap (u.year);
+
+  if (u.year < 100)
+    u.year += 2000;
+
+
+  if (!u.jday)
+    {
+      if (u.is_leap)
+        u.jday = u.mday + lmdays[u.month];
+      else
+        u.jday = u.mday + mdays[u.month];
+    }
+
+  u.year -= 1601;
+  y400 = u.year / 400;
+  u.year -= y400 * 400;
+  y100 = u.year / 100;
+  u.year -= y100 * 100;
+  y4 = u.year / 4;
+  u.year -= y4 * 4;
+
+  jd = u.jday - 1;
+  jd += u.year * 365;
+  jd += y4 * 1461;
+  jd += y100 * 36524;
+  jd += y400 * 146097;
+
+
+  jd += 2305813;
+  jd -= 2451545;
+
+
+
+  ret = (double) u.nanosecond;
+  ret /= 1000000000.;
+  ret = (double) u.second;
+  ret /= 60.;
+  ret += (double) u.minute;
+  ret /= 60.;
+  ret += (double) u.hour;
+  ret /= 24.;
+
+  ret += .5;
+
+  ret += (double) jd;
+
+  return ret;
+}
+
+double
+time_utc_to_ra (UTC u)
+{
+  double tjd = time_utc_to_tjd (u);
+  double T = tjd / 36525.;
+  double theta0 =
+    280.46061837 + (360.98564736629 * tjd) + (0.000387933 * T * T) -
+    (T * T * T / 38710000.0);
+
+  return remainder (theta0, 360.);
+}
+
+
+ST
+time_ra_to_st (double ra)
+{
+  ST ret;
+  unsigned i;
+
+  while (ra < 0.)
+    ra += 360.;
+  if (ra >= 360.0)
+    ra = remainder (ra, 360.);
+
+  ra *= 240.0;
+
+  i = (int) (floor (ra) + .5);
+  ra -= (double) i;
+  ra *= 1000000000.;
+
+  ret.nanosecond = (unsigned) (ra + .5);
+  ret.second = i % 60;
+  i /= 60;
+  ret.minute = i % 60;
+  i /= 60;
+  ret.hour = i;
+
+  return ret;
+}
+
+
+
+
+int
+main (int argc, char *argv[])
+{
+
+  UTC u = { 0 };
+  double local_lon = 0.0;
+  ST st;
+  double ra;
+
+  u.year = 2021;
+  u.month = 2;
+  u.mday = 12;
+  u.hour = 13;
+  u.minute =28;
+
+
+  ra = time_utc_to_ra (u);
+  st = time_ra_to_st (ra + local_lon);
+
+  printf ("%d:%02d:%02d\n", st.hour, st.minute, st.second);
+  return 0;
+}
+
+