1 files changed, 127 insertions, 6 deletions
diff --git a/app/time_fn.c b/app/time_fn.c
index fc3f829..64a92ff 100644
--- a/app/time_fn.c
+++ b/app/time_fn.c
@@ -15,6 +15,39 @@ static int is_leap (unsigned year)
   return 1;
 }
 
+EPOCH time_epoch_sub (EPOCH a, EPOCH b)
+{
+  EPOCH ret;
+
+  ret.ns = a.ns - b.ns;
+  ret.s = a.s - b.s;
+
+  if (!ret.s) return ret;
+
+  if (ret.ns < 0) {
+    ret.s--;
+    ret.ns += 1000000000;
+  }
+
+  return ret;
+}
+
+
+EPOCH time_epoch_add (EPOCH a, EPOCH b)
+{
+  EPOCH ret;
+
+  ret.ns = a.ns + b.ns;
+  ret.s = a.s + b.s;
+
+  while (ret.ns > 1000000000) {
+    ret.s++;
+    ret.ns -= 1000000000;
+  }
+
+  return ret;
+}
+
 
 UTC time_epoch_to_utc (EPOCH epoch)
 {
@@ -293,6 +326,20 @@ time_utc_to_tjd (UTC u)
   return ret;
 }
 
+
+double ra_normalize (double ra)
+{
+
+  while (ra < 0.)
+    ra += 360.;
+
+  if (ra >= 360.0)
+    ra = remainder (ra, 360.);
+
+  return ra;
+}
+
+
 double
 time_utc_to_ra (UTC u)
 {
@@ -302,9 +349,31 @@ time_utc_to_ra (UTC u)
     280.46061837 + (360.98564736629 * tjd) + (0.000387933 * T * T) -
     (T * T * T / 38710000.0);
 
-  return remainder (theta0, 360.);
+
+  return ra_normalize (remainder (theta0, 360.));
 }
 
+double time_st_to_ra (ST st)
+{
+  double ra;
+  unsigned i;
+
+  i = st.hour;
+  i *= 60;
+  i += st.minute;
+  i *= 60;
+  i += st.second;
+
+  ra = (double) st.nanosecond;
+  ra /= 1000000000.;
+  ra += (double) i;
+
+  ra /= 240.;
+
+  return ra;
+}
+
+
 
 ST
 time_ra_to_st (double ra)
@@ -312,11 +381,7 @@ time_ra_to_st (double ra)
   ST ret;
   unsigned i;
 
-  while (ra < 0.)
-    ra += 360.;
-
-  if (ra >= 360.0)
-    ra = remainder (ra, 360.);
+  ra = ra_normalize (ra);
 
   ra *= 240.0;
 
@@ -344,4 +409,60 @@ ST time_utc_to_lst (UTC u, double lon)
   return time_ra_to_st (ra + lon);
 }
 
+int time_ra_cmp (double a, double b)
+{
+  return (ra_normalize (a - b) < 180.) ? -1 : 1;;
+}
+
+
+/*Find the next occuring time for RA tra */
+
+EPOCH time_ra_to_next_epoch (EPOCH l, double tra)
+{
+  EPOCH r, m;
+  unsigned n;
+  UTC u;
+  double ra;
+
+  printf ("time_ra_to_next_epoch %.9f\n", tra);
+
+  /* XXX: we should use newton raphson, but */
+  /* 1) we're on team hooke */
+  /* 2) we want to limit the domain of solutions */
+  /* So IB works better */
+
+
+  r = l;
+  r.s += 86400.;
+
+  for (n = 0; n < 64; ++n) {
+
+    m = time_epoch_sub (r, l);
+
+    if (m.s < 0) {
+      m.s = 0;
+      m.ns = 0;
+      break;
+    }
+
+    if (m.s & 1)
+      m.ns += 1000000000;
+
+    m.s /= 2;
+    m.ns /= 2;
+
+    m = time_epoch_add (l, m);
+
+    u = time_epoch_to_utc (m);
+    ra = time_utc_to_ra (u);
+
+    if (time_ra_cmp (ra, tra) < 0)
+      r = m;
+
+    else
+      l = m;
+  }
+
+  return m;
+}