From f659c173bf1dda9e15285c5780fcccb2dd0fc6a2 Mon Sep 17 00:00:00 2001
From: Daniel Warner <contact@danrw.com>
Date: Sun, 11 Dec 2011 22:33:57 +0000
Subject: [PATCH] Correct date set in SGP4. Added original GMST calculation to
 Julian.cpp

---
 Eci.cpp    | 11 +++++++----
 Julian.cpp | 22 +++++++++++++++++-----
 SGP4.cpp   |  5 +++--
 3 files changed, 27 insertions(+), 11 deletions(-)

diff --git a/Eci.cpp b/Eci.cpp
index 79ce479..2aeff9e 100644
--- a/Eci.cpp
+++ b/Eci.cpp
@@ -1,5 +1,7 @@
 #include "Eci.h"
 
+#include "Globals.h"
+
 Eci::Eci(const Julian &date, const CoordGeodetic &geo)
 : date_(date) {
 
@@ -59,11 +61,12 @@ Eci::~Eci(void) {
 CoordGeodetic Eci::ToGeodetic() const {
 
     const double theta = AcTan(position_.y, position_.x);
-    /*
-     * changes lon to 0>= and <360
-     * const double lon = Globals::Fmod2p(theta - date_.ToGreenwichSiderealTime());
-     */
+
+    // 0 >= lon < 360
+    // const double lon = Fmod2p(theta - date_.ToGreenwichSiderealTime());
+    // 180 >= lon < 180
     const double lon = fmod(theta - date_.ToGreenwichSiderealTime(), kTWOPI);
+
     const double r = sqrt((position_.x * position_.x) + (position_.y * position_.y));
     static const double e2 = kF * (2.0 - kF);
 
diff --git a/Julian.cpp b/Julian.cpp
index 4b1e4cd..6f7214b 100644
--- a/Julian.cpp
+++ b/Julian.cpp
@@ -247,15 +247,27 @@ time_t Julian::ToTime() const {
 double Julian::ToGreenwichSiderealTime() const {
 
 #if 0
-    double theta;
-    double tut1;
+    const double UT = fmod(jul + 0.5, 1.0);
+    const double TU = (jul - 2451545.0 - UT) / 36525.0;
 
+    double GMST = 24110.54841 + TU *
+        (8640184.812866 + TU * (0.093104 - TU * 6.2e-06));
+
+    GMST = fmod(GMST + SEC_PER_DAY * OMEGA_E * UT, SEC_PER_DAY);
+
+    if (GMST < 0.0)
+        GMST += SEC_PER_DAY;  // "wrap" negative modulo value
+
+    return  (TWOPI * (GMST / SEC_PER_DAY));
+#endif
+
+#if 0
     // tut1 = Julian centuries from 2000 Jan. 1 12h UT1 (since J2000 which is 2451545.0)
     // a Julian century is 36525 days
-    tut1 = (date_ - 2451545.0) / 36525.0;
+    const double tut1 = (date_ - 2451545.0) / 36525.0;
 
     // Rotation angle in arcseconds
-    theta = 67310.54841 + (876600.0 * 3600.0 + 8640184.812866) * tut1 + 0.093104 * pow(tut1, 2.0) - 0.0000062 * pow(tut1, 3.0);
+    double theta = 67310.54841 + (876600.0 * 3600.0 + 8640184.812866) * tut1 + 0.093104 * pow(tut1, 2.0) - 0.0000062 * pow(tut1, 3.0);
 
     // 360.0 / 86400.0 = 1.0 / 240.0
     theta = fmod(DegreesToRadians(theta / 240.0), kTWOPI);
@@ -285,7 +297,7 @@ double Julian::ToGreenwichSiderealTime() const {
     const double c1p2p = C1 + kTWOPI;
     double gsto = fmod(THGR70 + C1 * ds70 + c1p2p * tfrac + ts70 * ts70 * FK5R, kTWOPI);
     if (gsto < 0.0)
-        gsto = gsto + kTWOPI;
+        gsto += kTWOPI;
 
     return gsto;
 }
diff --git a/SGP4.cpp b/SGP4.cpp
index fc3600b..72da522 100644
--- a/SGP4.cpp
+++ b/SGP4.cpp
@@ -519,8 +519,9 @@ void SGP4::CalculateFinalPositionVelocity(Eci* eci, const double& tsince, const
     const double zdot = (rdotk * uz + rfdotk * vz) * kXKMPER / 60.0;
     Vector velocity(xdot, ydot, zdot);
 
-    Timespan diff(tsince);
-    Julian julian = elements_.Epoch() + diff;
+    Julian julian(elements_.Epoch());
+    julian.AddMin(tsince);
+
     (*eci) = Eci(julian, position, velocity);
 }