From adc1c36a5b18a4ff83103df2674a307c14e08aec Mon Sep 17 00:00:00 2001
From: Daniel Warner <contact@danrw.com>
Date: Mon, 27 Aug 2018 12:29:44 +0100
Subject: [PATCH] Create function to compute perturbed values.

---
 libsgp4/SGP4.cc | 92 ++++++++++++++++++++++++++++---------------------
 libsgp4/SGP4.h  |  8 +++++
 2 files changed, 61 insertions(+), 39 deletions(-)

diff --git a/libsgp4/SGP4.cc b/libsgp4/SGP4.cc
index 2df4676..9dfe842 100644
--- a/libsgp4/SGP4.cc
+++ b/libsgp4/SGP4.cc
@@ -61,10 +61,16 @@ void SGP4::Initialise()
         throw SatelliteException("Inclination out of range");
     }
 
-    common_consts_.cosio = cos(elements_.Inclination());
-    common_consts_.sinio = sin(elements_.Inclination());
+    RecomputeConstants(elements_.Inclination(),
+                       common_consts_.sinio,
+                       common_consts_.cosio,
+                       common_consts_.x3thm1,
+                       common_consts_.x1mth2,
+                       common_consts_.x7thm1,
+                       common_consts_.xlcof,
+                       common_consts_.aycof);
+
     const double theta2 = common_consts_.cosio * common_consts_.cosio;
-    common_consts_.x3thm1 = 3.0 * theta2 - 1.0;
     const double eosq = elements_.Eccentricity() * elements_.Eccentricity();
     const double betao2 = 1.0 - eosq;
     const double betao = sqrt(betao2);
@@ -127,8 +133,6 @@ void SGP4::Initialise()
         + 0.75 * kCK2 * tsi / psisq * common_consts_.x3thm1
         * (8.0 + 3.0 * etasq * (8.0 + etasq)));
     common_consts_.c1 = elements_.BStar() * c2;
-    common_consts_.a3ovk2 = -kXJ3 / kCK2 * kAE * kAE * kAE;
-    common_consts_.x1mth2 = 1.0 - theta2;
     common_consts_.c4 = 2.0 * elements_.RecoveredMeanMotion()
         * coef1 * elements_.RecoveredSemiMajorAxis() * betao2
         * (common_consts_.eta * (2.0 + 0.5 * etasq) + elements_.Eccentricity()
@@ -155,18 +159,6 @@ void SGP4::Initialise()
     common_consts_.xnodcf = 3.5 * betao2 * xhdot1 * common_consts_.c1;
     common_consts_.t2cof = 1.5 * common_consts_.c1;
 
-    if (fabs(common_consts_.cosio + 1.0) > 1.5e-12)
-    {
-        common_consts_.xlcof = 0.125 * common_consts_.a3ovk2 * common_consts_.sinio * (3.0 + 5.0 * common_consts_.cosio) / (1.0 + common_consts_.cosio);
-    }
-    else
-    {
-        common_consts_.xlcof = 0.125 * common_consts_.a3ovk2 * common_consts_.sinio * (3.0 + 5.0 * common_consts_.cosio) / 1.5e-12;
-    }
-
-    common_consts_.aycof = 0.25 * common_consts_.a3ovk2 * common_consts_.sinio;
-    common_consts_.x7thm1 = 7.0 * theta2 - 1.0;
-
     if (use_deep_space_)
     {
         deepspace_consts_.gsto = elements_.Epoch().ToGreenwichSiderealTime();
@@ -316,29 +308,21 @@ Eci SGP4::FindPositionSDP4(double tsince) const
     /*
      * re-compute the perturbed values
      */
-    const double perturbed_sinio = sin(xincl);
-    const double perturbed_cosio = cos(xincl);
-
-    const double perturbed_theta2 = perturbed_cosio * perturbed_cosio;
-
-    const double perturbed_x3thm1 = 3.0 * perturbed_theta2 - 1.0;
-    const double perturbed_x1mth2 = 1.0 - perturbed_theta2;
-    const double perturbed_x7thm1 = 7.0 * perturbed_theta2 - 1.0;
-
+    double perturbed_sinio;
+    double perturbed_cosio;
+    double perturbed_x3thm1;
+    double perturbed_x1mth2;
+    double perturbed_x7thm1;
     double perturbed_xlcof;
-    if (fabs(perturbed_cosio + 1.0) > 1.5e-12)
-    {
-        perturbed_xlcof = 0.125 * common_consts_.a3ovk2 * perturbed_sinio
-            * (3.0 + 5.0 * perturbed_cosio) / (1.0 + perturbed_cosio);
-    }
-    else
-    {
-        perturbed_xlcof = 0.125 * common_consts_.a3ovk2 * perturbed_sinio
-            * (3.0 + 5.0 * perturbed_cosio) / 1.5e-12;
-    }
-
-    const double perturbed_aycof = 0.25 * common_consts_.a3ovk2
-        * perturbed_sinio;
+    double perturbed_aycof;
+    RecomputeConstants(xincl,
+                       perturbed_sinio,
+                       perturbed_cosio,
+                       perturbed_x3thm1,
+                       perturbed_x1mth2,
+                       perturbed_x7thm1,
+                       perturbed_xlcof,
+                       perturbed_aycof);
 
     /*
      * using calculated values, find position and velocity
@@ -359,6 +343,36 @@ Eci SGP4::FindPositionSDP4(double tsince) const
                                           perturbed_sinio);
 }
 
+void SGP4::RecomputeConstants(const double xincl,
+                              double& sinio,
+                              double& cosio,
+                              double& x3thm1,
+                              double& x1mth2,
+                              double& x7thm1,
+                              double& xlcof,
+                              double& aycof)
+{
+    sinio = sin(xincl);
+    cosio = cos(xincl);
+
+    const double theta2 = cosio * cosio;
+
+    x3thm1 = 3.0 * theta2 - 1.0;
+    x1mth2 = 1.0 - theta2;
+    x7thm1 = 7.0 * theta2 - 1.0;
+
+    if (fabs(cosio + 1.0) > 1.5e-12)
+    {
+        xlcof = 0.125 * kA3OVK2 * sinio * (3.0 + 5.0 * cosio) / (1.0 + cosio);
+    }
+    else
+    {
+        xlcof = 0.125 * kA3OVK2 * sinio * (3.0 + 5.0 * cosio) / 1.5e-12;
+    }
+
+    aycof = 0.25 * kA3OVK2 * sinio;
+}
+
 Eci SGP4::FindPositionSGP4(double tsince) const
 {
     /*
diff --git a/libsgp4/SGP4.h b/libsgp4/SGP4.h
index 7899624..fc181a4 100644
--- a/libsgp4/SGP4.h
+++ b/libsgp4/SGP4.h
@@ -189,6 +189,14 @@ private:
     };
     
     void Initialise();
+    static void RecomputeConstants(const double xincl,
+                                   double& sinio,
+                                   double& cosio,
+                                   double& x3thm1,
+                                   double& x1mth2,
+                                   double& x7thm1,
+                                   double& xlcof,
+                                   double& aycof);
     Eci FindPositionSDP4(const double tsince) const;
     Eci FindPositionSGP4(double tsince) const;
     Eci CalculateFinalPositionVelocity(