From ae8383d31628d8703f5d6a9cbeecace0790286bb Mon Sep 17 00:00:00 2001
From: Daniel Warner <contact@danrw.com>
Date: Sat, 2 Apr 2011 23:37:01 +0100
Subject: [PATCH] Updated THDT constant. Calculate real tsince

---
 Julian.cpp | 17 ++++++++++-------
 SGDP4.cpp  | 16 ++++++++++------
 main.cpp   |  6 +++---
 3 files changed, 23 insertions(+), 16 deletions(-)

diff --git a/Julian.cpp b/Julian.cpp
index 41f4d76..4f6a62f 100644
--- a/Julian.cpp
+++ b/Julian.cpp
@@ -212,7 +212,8 @@ time_t Julian::ToTime() const {
  * Greenwich Mean Sidereal Time
  */
 double Julian::ToGMST() const {
-    double value;
+
+    double theta;
     double tut1;
 
     // tut1 = Julian centuries from 2000 Jan. 1 12h UT1 (since J2000 which is 2451545.0)
@@ -220,16 +221,18 @@ double Julian::ToGMST() const {
     tut1 = (date_ - 2451545.0) / 36525.0;
 
     // Rotation angle in arcseconds
-    value = 67310.54841 + (876600.0 * 3600.0 + 8640184.812866) * tut1 + 0.093104 * pow(tut1, 2.0) - 0.0000062 * pow(tut1, 3.0);
+    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
-    value = fmod(Globals::Deg2Rad(value / 240.0), Globals::TWOPI());
+    theta = fmod(Globals::Deg2Rad(theta / 240.0), Globals::TWOPI());
 
-    // check quadrants
-    if (value < 0.0)
-        value += Globals::TWOPI();
+    /*
+     * check quadrants
+     */
+    if (theta < 0.0)
+        theta += Globals::TWOPI();
 
-    return value;
+    return theta;
 }
 
 /*
diff --git a/SGDP4.cpp b/SGDP4.cpp
index fb02d34..e7af933 100644
--- a/SGDP4.cpp
+++ b/SGDP4.cpp
@@ -247,10 +247,14 @@ void SGDP4::Initialize(const double& theta2, const double& betao2, const double&
 
 void SGDP4::FindPosition(Eci& eci, double tsince) {
 
+    Julian tsince_epoch = Epoch();
+    tsince_epoch.AddMin(tsince);
+    double actual_tsince = tsince_epoch.SpanMin(Epoch());
+
     if (i_use_deep_space_)
-        FindPositionSDP4(eci, tsince);
+        FindPositionSDP4(eci, actual_tsince);
     else
-        FindPositionSGP4(eci, tsince);
+        FindPositionSGP4(eci, actual_tsince);
 }
 
 void SGDP4::FindPositionSDP4(Eci& eci, double tsince) {
@@ -622,7 +626,7 @@ void SGDP4::DeepSpaceInitialize(const double& eosq, const double& sinio, const d
     static const double ROOT44 = 7.3636953E-9;
     static const double ROOT52 = 1.1428639E-7;
     static const double ROOT54 = 2.1765803E-9;
-    static const double THDT = 4.3752691E-3;
+    static const double THDT = 4.37526908801129966e-3;
 
     const double aqnv = 1.0 / RecoveredSemiMajorAxis();
     const double xpidot = omgdot + xnodot;
@@ -812,7 +816,7 @@ void SGDP4::DeepSpaceInitialize(const double& eosq, const double& sinio, const d
         d_fasx4_ = 2.8843198;
         d_fasx6_ = 0.37448087;
 
-        d_xlamo_ = fmod(MeanAnomoly() + AscendingNode() + ArgumentPerigee() - d_gsto_, Globals::TWOPI());
+        d_xlamo_ = MeanAnomoly() + AscendingNode() + ArgumentPerigee() - d_gsto_;
         bfact = xmdot + xpidot - THDT;
         bfact += d_ssl_ + d_ssg_ + d_ssh_;
 
@@ -909,7 +913,7 @@ void SGDP4::DeepSpaceInitialize(const double& eosq, const double& sinio, const d
         d_d5421_ = temp * f542 * g521;
         d_d5433_ = temp * f543 * g533;
 
-        d_xlamo_ = fmod(MeanAnomoly() + AscendingNode() + AscendingNode() - d_gsto_ - d_gsto_, Globals::TWOPI());
+        d_xlamo_ = MeanAnomoly() + AscendingNode() + AscendingNode() - d_gsto_ - d_gsto_;
         bfact = xmdot + xnodot + xnodot - THDT - THDT;
         bfact = bfact + d_ssl_ + d_ssh_ + d_ssh_;
     }
@@ -1012,7 +1016,7 @@ void SGDP4::DeepSpacePeriodics(const double& t, double& em,
         const double sinis = sin(xinc);
         const double cosis = cos(xinc);
 
-        if (xinc >= 0.2) {
+        if (Inclination() >= 0.2) {
             /*
              * apply periodics directly
              */
diff --git a/main.cpp b/main.cpp
index 5b1fd89..e6bcb02 100644
--- a/main.cpp
+++ b/main.cpp
@@ -144,7 +144,7 @@ void RunTle(Tle tle, double start, double end, double inc) {
 
 void RunTest() {
 
-#if 0
+//#if 0
 
     /*
       #   ------------------ Verification test cases ----------------------
@@ -177,7 +177,7 @@ void RunTest() {
             "1 06251U 62025E   06176.82412014  .00008885  00000-0  12808-3 0  3985",
             "2 06251  58.0579  54.0425 0030035 139.1568 221.1854 15.56387291  6774"), 0.0, 2880.0, 120.0);
 
-#endif
+//#endif
     /*
     #   MOLNIYA 2-14              # 12h resonant ecc in 0.65 to 0.7 range
     1 08195U 75081A   06176.33215444  .00000099  00000-0  11873-3 0   813
@@ -187,7 +187,7 @@ void RunTest() {
             "1 08195U 75081A   06176.33215444  .00000099  00000-0  11873-3 0   813",
             "2 08195  64.1586 279.0717 6877146 264.7651  20.2257  2.00491383225656"), 0.0, 2880.0, 120.0);
 
-    return;
+    //return;
     /*
     #   MOLNIYA 1-36              ## fig 12h resonant ecc in 0.7 to 0.715 range
     1 09880U 77021A   06176.56157475  .00000421  00000-0  10000-3 0  9814