/[PAMELA software]/YodaProfiler/src/sgp4.cpp

Diff of /YodaProfiler/src/sgp4.cpp

Parent Directory | Revision Log | View Patch Patch

-revision 1.1 by mocchiut,
Fri Oct 20 11:39:40 2006 UTC
+revision 1.2 by mocchiut,
Tue Jan 23 17:04:13 2007 UTC
 Line 6
  //////////////////////////////////////////////////////////////////////////////
  double sqr(const double x)
  {
-    return (x * x);
+   return (x * x);
  }
  //////////////////////////////////////////////////////////////////////////////
  double Fmod2p(const double arg)
  {
-    double modu = fmod(arg, TWOPI);
+   double modu = fmod(arg, TWOPI);
-    if (modu < 0.0)
+   if (modu < 0.0)
-       modu += TWOPI;
+     modu += TWOPI;
-    return modu;
+   return modu;
  }
  //////////////////////////////////////////////////////////////////////////////
 Line 26 
 double Fmod2p(const double arg)
  // is that it returns the correct quadrant of the angle.
  double AcTan(const double sinx, const double cosx)
  {
-    double ret;
+   double ret;
-    if (cosx == 0.0)
+   if (cosx == 0.0)
-    {
+     {
        if (sinx > 0.0)
-          ret = PI / 2.0;
+         ret = PI / 2.0;
        else
-          ret = 3.0 * PI / 2.0;
+         ret = 3.0 * PI / 2.0;
-    }
+     }
-    else
+   else
-    {
+     {
        if (cosx > 0.0)
-          ret = atan(sinx / cosx);
+         ret = atan(sinx / cosx);
        else
-          ret = PI + atan(sinx / cosx);
+         ret = PI + atan(sinx / cosx);
-    }
+     }
-    return ret;
+   return ret;
  }
  //////////////////////////////////////////////////////////////////////////////
  double rad2deg(const double r)
  {
-    const double DEG_PER_RAD = 180.0 / PI;
+   const double DEG_PER_RAD = 180.0 / PI;
-    return r * DEG_PER_RAD;
+   return r * DEG_PER_RAD;
  }
  //////////////////////////////////////////////////////////////////////////////
  double deg2rad(const double d)
  {
-    const double RAD_PER_DEG = PI / 180.0;
+   const double RAD_PER_DEG = PI / 180.0;
-    return d * RAD_PER_DEG;
+   return d * RAD_PER_DEG;
  }
  //
 Line 72 
 double deg2rad(const double d)
  cCoordGeo::cCoordGeo()
  {
-    m_Lat = 0.0;
+   m_Lat = 0.0;
-    m_Lon = 0.0;
+   m_Lon = 0.0;
-    m_Alt = 0.0;
+   m_Alt = 0.0;
  }
  //////////////////////////////////////////////////////////////////////
 Line 83 
 cCoordGeo::cCoordGeo()
  cCoordTopo::cCoordTopo()
  {
-    m_Az = 0.0;
+   m_Az = 0.0;
-    m_El = 0.0;
+   m_El = 0.0;
-    m_Range = 0.0;
+   m_Range = 0.0;
-    m_RangeRate = 0.0;
+   m_RangeRate = 0.0;
  }
 Line 102 
 cCoordTopo::cCoordTopo()
  //*****************************************************************************
  void cVector::Mul(double factor)
  {
-    m_x *= factor;
+   m_x *= factor;
-    m_y *= factor;
+   m_y *= factor;
-    m_z *= factor;
+   m_z *= factor;
-    m_w *= fabs(factor);
+   m_w *= fabs(factor);
  }
  //*****************************************************************************
 Line 113 
 void cVector::Mul(double factor)
  //*****************************************************************************
  void cVector::Sub(const cVector& vec)
  {
-    m_x -= vec.m_x;
+   m_x -= vec.m_x;
-    m_y -= vec.m_y;
+   m_y -= vec.m_y;
-    m_z -= vec.m_z;
+   m_z -= vec.m_z;
-    m_w -= vec.m_w;
+   m_w -= vec.m_w;
  }
  //*****************************************************************************
 Line 142 
 double cVector::Magnitude() const
  //*****************************************************************************
  double cVector::Dot(const cVector& vec) const
  {
-    return (m_x * vec.m_x) +
+   return (m_x * vec.m_x) +
-           (m_y * vec.m_y) +
+     (m_y * vec.m_y) +
-           (m_z * vec.m_z);
+     (m_z * vec.m_z);
  }
  //
  // cJulian.cpp
 Line 173 
 double cVector::Dot(const cVector& vec)
  // number of seconds since midnight UTC January 1, 1970.
  cJulian::cJulian(time_t time)
  {
-    struct tm *ptm = gmtime(&time);
+   struct tm *ptm = gmtime(&time);
-    assert(ptm);
+   assert(ptm);
-    int    year = ptm->tm_year + 1900;
+   int    year = ptm->tm_year + 1900;
-    double day  = ptm->tm_yday + 1 +
+   double day  = ptm->tm_yday + 1 +
-                  (ptm->tm_hour +
+     (ptm->tm_hour +
-                   ((ptm->tm_min +
+      ((ptm->tm_min +
-                    (ptm->tm_sec / 60.0)) / 60.0)) / 24.0;
+        (ptm->tm_sec / 60.0)) / 60.0)) / 24.0;
-    Initialize(year, day);
+   Initialize(year, day);
  }
  //////////////////////////////////////////////////////////////////////////////
 Line 190 
 cJulian::cJulian(time_t time)
  // Example parameters: year = 2001, day = 1.5 (Jan 1 12h)
  cJulian::cJulian(int year, double day)
  {
-    Initialize(year, day);
+   Initialize(year, day);
  }
  //////////////////////////////////////////////////////////////////////////////
 Line 203 
 cJulian::cJulian(int year,
                   double sec /* = 0.0 */) // 0..(59.999999...)
  {
-    // Calculate N, the day of the year (1..366)
+   // Calculate N, the day of the year (1..366)
-    int N;
+   int N;
-    int F1 = (int)((275.0 * mon) / 9.0);
+   int F1 = (int)((275.0 * mon) / 9.0);
-    int F2 = (int)((mon + 9.0) / 12.0);
+   int F2 = (int)((mon + 9.0) / 12.0);
-    if (IsLeapYear(year))
+   if (IsLeapYear(year))
-    {
+     {
        // Leap year
        N = F1 - F2 + day - 30;
-    }
+     }
-    else
+   else
-    {
+     {
        // Common year
        N = F1 - (2 * F2) + day - 30;
-    }
+     }
-    double dblDay = N + (hour + (min + (sec / 60.0)) / 60.0) / 24.0;
+   double dblDay = N + (hour + (min + (sec / 60.0)) / 60.0) / 24.0;
-    Initialize(year, dblDay);
+   Initialize(year, dblDay);
  }
  //////////////////////////////////////////////////////////////////////////////
  void cJulian::Initialize(int year, double day)
  {
-    // 1582 A.D.: 10 days removed from calendar
+   // 1582 A.D.: 10 days removed from calendar
-    // 3000 A.D.: Arbitrary error checking limit
+   // 3000 A.D.: Arbitrary error checking limit
-    assert((year > 1582) && (year < 3000));
+   assert((year > 1582) && (year < 3000));
-    assert((day >= 0.0) && (day <= 366.5));
+   assert((day >= 0.0) && (day <= 366.5));
-    // Now calculate Julian date
+   // Now calculate Julian date
-    year--;
+   year--;
-    // Centuries are not leap years unless they divide by 400
+   // Centuries are not leap years unless they divide by 400
-    int A = (year / 100);
+   int A = (year / 100);
-    int B = 2 - A + (A / 4);
+   int B = 2 - A + (A / 4);
-    double NewYears = (int)(365.25 * year) +
+   double NewYears = (int)(365.25 * year) +
-                      (int)(30.6001 * 14)  +
+     (int)(30.6001 * 14)  +
-                      1720994.5 + B;  // 1720994.5 = Oct 30, year -1
+     1720994.5 + B;  // 1720994.5 = Oct 30, year -1
-    m_Date = NewYears + day;
+   m_Date = NewYears + day;
  }
  //////////////////////////////////////////////////////////////////////////////
 Line 257 
 void cJulian::getComponent(int    *pYear
                             int    *pMon  /* = NULL */,
                             double *pDOM  /* = NULL */) const
  {
-    assert(pYear != NULL);
+   assert(pYear != NULL);
-    double jdAdj = getDate() + 0.5;
+   double jdAdj = getDate() + 0.5;
-    int    Z     = (int)jdAdj;  // integer part
+   int    Z     = (int)jdAdj;  // integer part
-    double F     = jdAdj - Z;   // fractional part
+   double F     = jdAdj - Z;   // fractional part
-    double alpha = (int)((Z - 1867216.25) / 36524.25);
+   double alpha = (int)((Z - 1867216.25) / 36524.25);
-    double A     = Z + 1 + alpha - (int)(alpha / 4.0);
+   double A     = Z + 1 + alpha - (int)(alpha / 4.0);
-    double B     = A + 1524.0;
+   double B     = A + 1524.0;
-    int    C     = (int)((B - 122.1) / 365.25);
+   int    C     = (int)((B - 122.1) / 365.25);
-    int    D     = (int)(C * 365.25);
+   int    D     = (int)(C * 365.25);
-    int    E     = (int)((B - D) / 30.6001);
+   int    E     = (int)((B - D) / 30.6001);
-    double DOM   = B - D - (int)(E * 30.6001) + F;
+   double DOM   = B - D - (int)(E * 30.6001) + F;
-    int    month = (E < 13.5) ? (E - 1) : (E - 13);
+   int    month = (E < 13.5) ? (E - 1) : (E - 13);
-    int    year  = (month > 2.5) ? (C - 4716) : (C - 4715);
+   int    year  = (month > 2.5) ? (C - 4716) : (C - 4715);
-    *pYear = year;
+   *pYear = year;
-    if (pMon != NULL)
+   if (pMon != NULL)
-       *pMon = month;
+     *pMon = month;
-    if (pDOM != NULL)
+   if (pDOM != NULL)
-       *pDOM = DOM;
+     *pDOM = DOM;
  }
  //////////////////////////////////////////////////////////////////////////////
 Line 295 
 void cJulian::getComponent(int    *pYear
  //       Nov/Dec 1995
  double cJulian::toGMST() const
  {
-    const double UT = fmod(m_Date + 0.5, 1.0);
+   const double UT = fmod(m_Date + 0.5, 1.0);
-    const double TU = (FromJan1_12h_2000() - UT) / 36525.0;
+   const double TU = (FromJan1_12h_2000() - UT) / 36525.0;
-    double GMST = 24110.54841 + TU *
+   double GMST = 24110.54841 + TU *
-                  (8640184.812866 + TU * (0.093104 - TU * 6.2e-06));
+     (8640184.812866 + TU * (0.093104 - TU * 6.2e-06));
-    GMST = fmod(GMST + SEC_PER_DAY * OMEGA_E * UT, SEC_PER_DAY);
+   GMST = fmod(GMST + SEC_PER_DAY * OMEGA_E * UT, SEC_PER_DAY);
-    if (GMST < 0.0)
+   if (GMST < 0.0)
-       GMST += SEC_PER_DAY;  // "wrap" negative modulo value
+     GMST += SEC_PER_DAY;  // "wrap" negative modulo value
-    return  (TWOPI * (GMST / SEC_PER_DAY));
+   return  (TWOPI * (GMST / SEC_PER_DAY));
-  }
+ }
  //////////////////////////////////////////////////////////////////////////////
  // toLMST()
 Line 317 
 double cJulian::toGMST() const
  // Vernal Equinox to the given longitude.
  double cJulian::toLMST(double lon) const
  {
-    return fmod(toGMST() + lon, TWOPI);
+   return fmod(toGMST() + lon, TWOPI);
  }
  //////////////////////////////////////////////////////////////////////////////
 Line 327 
 double cJulian::toLMST(double lon) const
  // time component.
  time_t cJulian::toTime() const
  {
-    int    nYear;
+   int    nYear;
-    int    nMonth;
+   int    nMonth;
-    double dblDay;
+   double dblDay;
-    getComponent(&nYear, &nMonth, &dblDay);
+   getComponent(&nYear, &nMonth, &dblDay);
-    // dblDay is the fractional Julian Day (i.e., 29.5577).
+   // dblDay is the fractional Julian Day (i.e., 29.5577).
-    // Save the whole number day in nDOM and convert dblDay to
+   // Save the whole number day in nDOM and convert dblDay to
-    // the fractional portion of day.
+   // the fractional portion of day.
-    int nDOM = (int)dblDay;
+   int nDOM = (int)dblDay;
-    dblDay -= nDOM;
+   dblDay -= nDOM;
-    const int SEC_PER_MIN = 60;
+   const int SEC_PER_MIN = 60;
-    const int SEC_PER_HR  = 60 * SEC_PER_MIN;
+   const int SEC_PER_HR  = 60 * SEC_PER_MIN;
-    const int SEC_PER_DAY = 24 * SEC_PER_HR;
+   const int SEC_PER_DAY = 24 * SEC_PER_HR;
-    int secs = (int)((dblDay * SEC_PER_DAY) + 0.5);
+   int secs = (int)((dblDay * SEC_PER_DAY) + 0.5);
-    // Create a "struct tm" type.
+   // Create a "struct tm" type.
-    // NOTE:
+   // NOTE:
-    // The "struct tm" type has a 1-second resolution. Any fractional
+   // The "struct tm" type has a 1-second resolution. Any fractional
-    // component of the "seconds" time value is discarded.
+   // component of the "seconds" time value is discarded.
-    struct tm tGMT;
+   struct tm tGMT;
-    memset(&tGMT, 0, sizeof(tGMT));
+   memset(&tGMT, 0, sizeof(tGMT));
-    tGMT.tm_year = nYear - 1900;  // 2001 is 101
+   tGMT.tm_year = nYear - 1900;  // 2001 is 101
-    tGMT.tm_mon  = nMonth - 1;    // January is 0
+   tGMT.tm_mon  = nMonth - 1;    // January is 0
-    tGMT.tm_mday = nDOM;          // First day is 1
+   tGMT.tm_mday = nDOM;          // First day is 1
-    tGMT.tm_hour = secs / SEC_PER_HR;
+   tGMT.tm_hour = secs / SEC_PER_HR;
-    tGMT.tm_min  = (secs % SEC_PER_HR) / SEC_PER_MIN;
+   tGMT.tm_min  = (secs % SEC_PER_HR) / SEC_PER_MIN;
-    tGMT.tm_sec  = (secs % SEC_PER_HR) % SEC_PER_MIN;
+   tGMT.tm_sec  = (secs % SEC_PER_HR) % SEC_PER_MIN;
-    tGMT.tm_isdst = 0; // No conversion desired
+   tGMT.tm_isdst = 0; // No conversion desired
-    time_t tEpoch = mktime(&tGMT);
+   time_t tEpoch = mktime(&tGMT);
-    if (tEpoch != -1)
+   if (tEpoch != -1)
-    {
+     {
        // Valid time_t value returned from mktime().
        // mktime() expects a local time which means that tEpoch now needs
        // to be adjusted by the difference between this time zone and GMT.
        tEpoch -= timezone;
-    }
+     }
-    return tEpoch;
+   return tEpoch;
  }
  //
  // cTle.cpp
 Line 410 
 const int TLE2_COL_REVATEPOCH    = 63; c
  /////////////////////////////////////////////////////////////////////////////
  cTle::cTle(string& strName, string& strLine1, string& strLine2)
  {
-    m_strName  = strName;
+   m_strName  = strName;
-    m_strLine1 = strLine1;
+   m_strLine1 = strLine1;
-    m_strLine2 = strLine2;
+   m_strLine2 = strLine2;
-    Initialize();
+   Initialize();
  }
  /////////////////////////////////////////////////////////////////////////////
  cTle::cTle(const cTle &tle)
  {
-    m_strName  = tle.m_strName;
+   m_strName  = tle.m_strName;
-    m_strLine1 = tle.m_strLine1;
+   m_strLine1 = tle.m_strLine1;
-    m_strLine2 = tle.m_strLine2;
+   m_strLine2 = tle.m_strLine2;
-    for (int fld = FLD_FIRST; fld < FLD_LAST; fld++)
+   for (int fld = FLD_FIRST; fld < FLD_LAST; fld++)
-    {
+     {
        m_Field[fld] = tle.m_Field[fld];
-    }
+     }
-    m_mapCache = tle.m_mapCache;
+   m_mapCache = tle.m_mapCache;
  }
  /////////////////////////////////////////////////////////////////////////////
 Line 450 
 double cTle::getField(eField   fld,
                        string  *pstr      /* = NULL     */,
                        bool     bStrUnits /* = false    */) const
  {
-    assert((FLD_FIRST <= fld) && (fld < FLD_LAST));
+   assert((FLD_FIRST <= fld) && (fld < FLD_LAST));
-    assert((U_FIRST <= units) && (units < U_LAST));
+   assert((U_FIRST <= units) && (units < U_LAST));
-    if (pstr)
+   if (pstr)
-    {
+     {
        // Return requested field in string form.
        *pstr = m_Field[fld];
        if (bStrUnits)
-          *pstr += getUnits(fld);
+         *pstr += getUnits(fld);
        return 0.0;
-    }
+     }
-    else
+   else
-    {
+     {
        // Return requested field in floating-point form.
        // Return cache contents if it exists, else populate cache
        FldKey key = Key(units, fld);
        if (m_mapCache.find(key) == m_mapCache.end())
-       {
+         {
-          // Value not in cache; add it
+           // Value not in cache; add it
-          double valNative = atof(m_Field[fld].c_str());
+           double valNative = atof(m_Field[fld].c_str());
-          double valConv   = ConvertUnits(valNative, fld, units);
+           double valConv   = ConvertUnits(valNative, fld, units);
-          m_mapCache[key]  = valConv;
+           m_mapCache[key]  = valConv;
-          return valConv;
+           return valConv;
-       }
+         }
        else
-       {
+         {
-          // return cached value
+           // return cached value
-          return m_mapCache[key];
+           return m_mapCache[key];
-       }
+         }
-    }
+     }
  }
  //////////////////////////////////////////////////////////////////////////////
 Line 493 
 double cTle::ConvertUnits(double valNati
                            eField fld,       // what field the value is
                            eUnits units)     // what units to convert to
  {
-    switch (fld)
+   switch (fld)
-    {
+     {
-    case FLD_I:
+     case FLD_I:
-    case FLD_RAAN:
+     case FLD_RAAN:
-    case FLD_ARGPER:
+     case FLD_ARGPER:
-    case FLD_M:
+     case FLD_M:
-      {
+       {
-        // The native TLE format is DEGREES
+         // The native TLE format is DEGREES
-        if (units == U_RAD)
+         if (units == U_RAD)
-          return valNative * RADS_PER_DEG;
+           return valNative * RADS_PER_DEG;
-      }
+       }
-    case FLD_NORADNUM:
+     case FLD_NORADNUM:
-    case FLD_INTLDESC:
+     case FLD_INTLDESC:
-    case FLD_SET:
+     case FLD_SET:
-    case FLD_EPOCHYEAR:
+     case FLD_EPOCHYEAR:
-    case FLD_EPOCHDAY:
+     case FLD_EPOCHDAY:
-    case FLD_ORBITNUM:
+     case FLD_ORBITNUM:
-    case FLD_E:
+     case FLD_E:
-    case FLD_MMOTION:
+     case FLD_MMOTION:
-    case FLD_MMOTIONDT:
+     case FLD_MMOTIONDT:
-    case FLD_MMOTIONDT2:
+     case FLD_MMOTIONDT2:
-    case FLD_BSTAR:
+     case FLD_BSTAR:
-    case FLD_LAST:
+     case FLD_LAST:
-      { // do nothing
+       { // do nothing
-      }
+       }
-    }
+     }
-    return valNative; // return value in unconverted native format
+   return valNative; // return value in unconverted native format
  }
  //////////////////////////////////////////////////////////////////////////////
  string cTle::getUnits(eField fld) const
  {
-    static const string strDegrees    = " degrees";
+   static const string strDegrees    = " degrees";
-    static const string strRevsPerDay = " revs / day";
+   static const string strRevsPerDay = " revs / day";
-    static const string strNull;
+   static const string strNull;
-    switch (fld)
+   switch (fld)
-    {
+     {
-       case FLD_I:
+     case FLD_I:
-       case FLD_RAAN:
+     case FLD_RAAN:
-       case FLD_ARGPER:
+     case FLD_ARGPER:
-       case FLD_M:
+     case FLD_M:
-          return strDegrees;
+       return strDegrees;
-       case FLD_MMOTION:
+     case FLD_MMOTION:
-          return strRevsPerDay;
+       return strRevsPerDay;
-       default:
+     default:
-          return strNull;
+       return strNull;
-    }
+     }
  }
  /////////////////////////////////////////////////////////////////////////////
 Line 559 
 string cTle::getUnits(eField fld) const
  //       " 40436+1" =  4.0436
  string cTle::ExpToDecimal(const string &str)
  {
-    const int COL_EXP_SIGN = 6;
+   const int COL_EXP_SIGN = 6;
-    const int LEN_EXP      = 2;
+   const int LEN_EXP      = 2;
-    const int LEN_BUFREAL  = 32;   // max length of buffer to hold floating point
+   const int LEN_BUFREAL  = 32;   // max length of buffer to hold floating point
                                   // representation of input string.
-    int nMan;
+   int nMan;
-    int nExp;
+   int nExp;
-    // sscanf(%d) will read up to the exponent sign
+   // sscanf(%d) will read up to the exponent sign
-    sscanf(str.c_str(), "%d", &nMan);
+   sscanf(str.c_str(), "%d", &nMan);
-    double dblMan = nMan;
+   double dblMan = nMan;
-    bool  bNeg    = (nMan < 0);
+   bool  bNeg    = (nMan < 0);
-    if (bNeg)
+   if (bNeg)
-       dblMan *= -1;
+     dblMan *= -1;
-    // Move decimal place to left of first digit
+   // Move decimal place to left of first digit
-    while (dblMan >= 1.0)
+   while (dblMan >= 1.0)
-       dblMan /= 10.0;
+     dblMan /= 10.0;
-    if (bNeg)
+   if (bNeg)
-       dblMan *= -1;
+     dblMan *= -1;
-    // now read exponent
+   // now read exponent
-    sscanf(str.substr(COL_EXP_SIGN, LEN_EXP).c_str(), "%d", &nExp);
+   sscanf(str.substr(COL_EXP_SIGN, LEN_EXP).c_str(), "%d", &nExp);
-    double dblVal = dblMan * pow(10.0, nExp);
+   double dblVal = dblMan * pow(10.0, nExp);
-    char szVal[LEN_BUFREAL];
+   char szVal[LEN_BUFREAL];
-    snprintf(szVal, sizeof(szVal), "%.9f", dblVal);
+   snprintf(szVal, sizeof(szVal), "%.9f", dblVal);
-    string strVal = szVal;
+   string strVal = szVal;
-    return strVal;
+   return strVal;
  } // ExpToDecimal()
 Line 602 
 string cTle::ExpToDecimal(const string &
  // Initialize the string array.
  void cTle::Initialize()
  {
-    // Have we already been initialized?
+   // Have we already been initialized?
-    if (m_Field[FLD_NORADNUM].size())
+   if (m_Field[FLD_NORADNUM].size())
-       return;
+     return;
-    assert(!m_strName.empty());
+   assert(!m_strName.empty());
-    assert(!m_strLine1.empty());
+   assert(!m_strLine1.empty());
-    assert(!m_strLine2.empty());
+   assert(!m_strLine2.empty());
-    m_Field[FLD_NORADNUM] = m_strLine1.substr(TLE1_COL_SATNUM, TLE1_LEN_SATNUM);
+   m_Field[FLD_NORADNUM] = m_strLine1.substr(TLE1_COL_SATNUM, TLE1_LEN_SATNUM);
-    m_Field[FLD_INTLDESC] = m_strLine1.substr(TLE1_COL_INTLDESC_A,
+   m_Field[FLD_INTLDESC] = m_strLine1.substr(TLE1_COL_INTLDESC_A,
-                                              TLE1_LEN_INTLDESC_A +
+                                             TLE1_LEN_INTLDESC_A +
-                                              TLE1_LEN_INTLDESC_B +
+                                             TLE1_LEN_INTLDESC_B +
-                                              TLE1_LEN_INTLDESC_C);
+                                             TLE1_LEN_INTLDESC_C);
-    m_Field[FLD_EPOCHYEAR] =
+   m_Field[FLD_EPOCHYEAR] =
-          m_strLine1.substr(TLE1_COL_EPOCH_A, TLE1_LEN_EPOCH_A);
+     m_strLine1.substr(TLE1_COL_EPOCH_A, TLE1_LEN_EPOCH_A);
-    m_Field[FLD_EPOCHDAY] =
+   m_Field[FLD_EPOCHDAY] =
-          m_strLine1.substr(TLE1_COL_EPOCH_B, TLE1_LEN_EPOCH_B);
+     m_strLine1.substr(TLE1_COL_EPOCH_B, TLE1_LEN_EPOCH_B);
-    if (m_strLine1[TLE1_COL_MEANMOTIONDT] == '-')
+   if (m_strLine1[TLE1_COL_MEANMOTIONDT] == '-')
-    {
+     {
        // value is negative
        m_Field[FLD_MMOTIONDT] = "-0";
-    }
+     }
-    else
+   else
-       m_Field[FLD_MMOTIONDT] = "0";
+     m_Field[FLD_MMOTIONDT] = "0";
-    m_Field[FLD_MMOTIONDT] += m_strLine1.substr(TLE1_COL_MEANMOTIONDT + 1,
+   m_Field[FLD_MMOTIONDT] += m_strLine1.substr(TLE1_COL_MEANMOTIONDT + 1,
-                                                TLE1_LEN_MEANMOTIONDT);
+                                               TLE1_LEN_MEANMOTIONDT);
-    // decimal point assumed; exponential notation
+   // decimal point assumed; exponential notation
-    m_Field[FLD_MMOTIONDT2] = ExpToDecimal(
+   m_Field[FLD_MMOTIONDT2] = ExpToDecimal(
-                                  m_strLine1.substr(TLE1_COL_MEANMOTIONDT2,
+                                          m_strLine1.substr(TLE1_COL_MEANMOTIONDT2,
-                                                    TLE1_LEN_MEANMOTIONDT2));
+                                                            TLE1_LEN_MEANMOTIONDT2));
-    // decimal point assumed; exponential notation
+   // decimal point assumed; exponential notation
-    m_Field[FLD_BSTAR] = ExpToDecimal(m_strLine1.substr(TLE1_COL_BSTAR,
+   m_Field[FLD_BSTAR] = ExpToDecimal(m_strLine1.substr(TLE1_COL_BSTAR,
-                                                        TLE1_LEN_BSTAR));
+                                                       TLE1_LEN_BSTAR));
-    //TLE1_COL_EPHEMTYPE
+   //TLE1_COL_EPHEMTYPE
-    //TLE1_LEN_EPHEMTYPE
+   //TLE1_LEN_EPHEMTYPE
-    m_Field[FLD_SET] = m_strLine1.substr(TLE1_COL_ELNUM, TLE1_LEN_ELNUM);
+   m_Field[FLD_SET] = m_strLine1.substr(TLE1_COL_ELNUM, TLE1_LEN_ELNUM);
-    TrimLeft(m_Field[FLD_SET]);
+   TrimLeft(m_Field[FLD_SET]);
-    //TLE2_COL_SATNUM
+   //TLE2_COL_SATNUM
-    //TLE2_LEN_SATNUM
+   //TLE2_LEN_SATNUM
-    m_Field[FLD_I] = m_strLine2.substr(TLE2_COL_INCLINATION,
+   m_Field[FLD_I] = m_strLine2.substr(TLE2_COL_INCLINATION,
-                                       TLE2_LEN_INCLINATION);
+                                      TLE2_LEN_INCLINATION);
-    TrimLeft(m_Field[FLD_I]);
+   TrimLeft(m_Field[FLD_I]);
-    m_Field[FLD_RAAN] = m_strLine2.substr(TLE2_COL_RAASCENDNODE,
+   m_Field[FLD_RAAN] = m_strLine2.substr(TLE2_COL_RAASCENDNODE,
-                                          TLE2_LEN_RAASCENDNODE);
+                                         TLE2_LEN_RAASCENDNODE);
-    TrimLeft(m_Field[FLD_RAAN]);
+   TrimLeft(m_Field[FLD_RAAN]);
-    // decimal point is assumed
+   // decimal point is assumed
-    m_Field[FLD_E]   = "0.";
+   m_Field[FLD_E]   = "0.";
-    m_Field[FLD_E]   += m_strLine2.substr(TLE2_COL_ECCENTRICITY,
+   m_Field[FLD_E]   += m_strLine2.substr(TLE2_COL_ECCENTRICITY,
-                                        TLE2_LEN_ECCENTRICITY);
+                                         TLE2_LEN_ECCENTRICITY);
-    m_Field[FLD_ARGPER] = m_strLine2.substr(TLE2_COL_ARGPERIGEE,
+   m_Field[FLD_ARGPER] = m_strLine2.substr(TLE2_COL_ARGPERIGEE,
-                                            TLE2_LEN_ARGPERIGEE);
+                                           TLE2_LEN_ARGPERIGEE);
-    TrimLeft(m_Field[FLD_ARGPER]);
+   TrimLeft(m_Field[FLD_ARGPER]);
-    m_Field[FLD_M]   = m_strLine2.substr(TLE2_COL_MEANANOMALY,
+   m_Field[FLD_M]   = m_strLine2.substr(TLE2_COL_MEANANOMALY,
-                                       TLE2_LEN_MEANANOMALY);
+                                        TLE2_LEN_MEANANOMALY);
-    TrimLeft(m_Field[FLD_M]);
+   TrimLeft(m_Field[FLD_M]);
-    m_Field[FLD_MMOTION]   = m_strLine2.substr(TLE2_COL_MEANMOTION,
+   m_Field[FLD_MMOTION]   = m_strLine2.substr(TLE2_COL_MEANMOTION,
-                                             TLE2_LEN_MEANMOTION);
+                                              TLE2_LEN_MEANMOTION);
-    TrimLeft(m_Field[FLD_MMOTION]);
+   TrimLeft(m_Field[FLD_MMOTION]);
-    m_Field[FLD_ORBITNUM] = m_strLine2.substr(TLE2_COL_REVATEPOCH,
+   m_Field[FLD_ORBITNUM] = m_strLine2.substr(TLE2_COL_REVATEPOCH,
-                                              TLE2_LEN_REVATEPOCH);
+                                             TLE2_LEN_REVATEPOCH);
-    TrimLeft(m_Field[FLD_ORBITNUM]);
+   TrimLeft(m_Field[FLD_ORBITNUM]);
  } // InitStrVars()
 Line 692 
 void cTle::Initialize()
  //
  bool cTle::IsValidLine(string& str, eTleLine line)
  {
-    TrimLeft(str);
+   TrimLeft(str);
-    TrimRight(str);
+   TrimRight(str);
-    size_t nLen = str.size();
+   size_t nLen = str.size();
-    if (nLen != (uint)TLE_LEN_LINE_DATA)
+   if (nLen != (uint)TLE_LEN_LINE_DATA)
-       return false;
+     return false;
-    // First char in string must be line number
+   // First char in string must be line number
-    if ((str[0] - '0') != line)
+   if ((str[0] - '0') != line)
-       return false;
+     return false;
-    // Second char in string must be blank
+   // Second char in string must be blank
-    if (str[1] != ' ')
+   if (str[1] != ' ')
-       return false;
+     return false;
-    /*
+   /*
-       NOTE: 12/96
+     NOTE: 12/96
-       The requirement that the last char in the line data must be a valid
+     The requirement that the last char in the line data must be a valid
-       checksum is too restrictive.
+     checksum is too restrictive.
-       // Last char in string must be checksum
+     // Last char in string must be checksum
-       int nSum = CheckSum(str);
+     int nSum = CheckSum(str);
-       if (nSum != (str[TLE_LEN_LINE_DATA - 1] - '0'))
+     if (nSum != (str[TLE_LEN_LINE_DATA - 1] - '0'))
-          return false;
+     return false;
-    */
+   */
-    return true;
+   return true;
  } // IsTleFormat()
 Line 735 
 bool cTle::IsValidLine(string& str, eTle
  //
  int cTle::CheckSum(const string& str)
  {
-    // The length is "- 1" because we don't include the current (existing)
+   // The length is "- 1" because we don't include the current (existing)
-    // checksum character in the checksum calculation.
+   // checksum character in the checksum calculation.
-    size_t len = str.size() - 1;
+   size_t len = str.size() - 1;
-    int xsum = 0;
+   int xsum = 0;
-    for (size_t i = 0; i < len; i++)
+   for (size_t i = 0; i < len; i++)
-    {
+     {
        char ch = str[i];
        if (isdigit(ch))
-          xsum += (ch - '0');
+         xsum += (ch - '0');
        else if (ch == '-')
-          xsum++;
+         xsum++;
-    }
+     }
-    return (xsum % 10);
+   return (xsum % 10);
  } // CheckSum()
  /////////////////////////////////////////////////////////////////////////////
  void cTle::TrimLeft(string& s)
  {
-    while (s[0] == ' ')
+   while (s[0] == ' ')
-       s.erase(0, 1);
+     s.erase(0, 1);
  }
  /////////////////////////////////////////////////////////////////////////////
  void cTle::TrimRight(string& s)
  {
-    while (s[s.size() - 1] == ' ')
+   while (s[s.size() - 1] == ' ')
-       s.erase(s.size() - 1);
+     s.erase(s.size() - 1);
  }
  //
 Line 780 
 cEci::cEci(const cVector &pos,
             const cJulian &date,
             bool  IsAeUnits /* = true */)
  {
-    m_pos      = pos;
+   m_pos      = pos;
-    m_vel      = vel;
+   m_vel      = vel;
-    m_date     = date;
+   m_date     = date;
-    m_VecUnits = (IsAeUnits ? UNITS_AE : UNITS_NONE);
+   m_VecUnits = (IsAeUnits ? UNITS_AE : UNITS_NONE);
  }
  //////////////////////////////////////////////////////////////////////
 Line 795 
 cEci::cEci(const cVector &pos,
  // Reference: www.celestrak.com (Dr. TS Kelso)
  cEci::cEci(const cCoordGeo &geo, const cJulian &date)
  {
-    m_VecUnits = UNITS_KM;
+   m_VecUnits = UNITS_KM;
-    double mfactor = TWOPI * (OMEGA_E / SEC_PER_DAY);
+   double mfactor = TWOPI * (OMEGA_E / SEC_PER_DAY);
-    double lat = geo.m_Lat;
+   double lat = geo.m_Lat;
-    double lon = geo.m_Lon;
+   double lon = geo.m_Lon;
-    double alt = geo.m_Alt;
+   double alt = geo.m_Alt;
-    // Calculate Local Mean Sidereal Time (theta)
+   // Calculate Local Mean Sidereal Time (theta)
-    double theta = date.toLMST(lon);
+   double theta = date.toLMST(lon);
-    double c = 1.0 / sqrt(1.0 + F * (F - 2.0) * sqr(sin(lat)));
+   double c = 1.0 / sqrt(1.0 + F * (F - 2.0) * sqr(sin(lat)));
-    double s = sqr(1.0 - F) * c;
+   double s = sqr(1.0 - F) * c;
-    double achcp = (XKMPER_WGS72 * c + alt) * cos(lat);
+   double achcp = (XKMPER_WGS72 * c + alt) * cos(lat);
-    m_date = date;
+   m_date = date;
-    m_pos.m_x = achcp * cos(theta);               // km
+   m_pos.m_x = achcp * cos(theta);               // km
-    m_pos.m_y = achcp * sin(theta);               // km
+   m_pos.m_y = achcp * sin(theta);               // km
-    m_pos.m_z = (XKMPER_WGS72 * s + alt) * sin(lat);   // km
+   m_pos.m_z = (XKMPER_WGS72 * s + alt) * sin(lat);   // km
-    m_pos.m_w = sqrt(sqr(m_pos.m_x) +
+   m_pos.m_w = sqrt(sqr(m_pos.m_x) +
-                     sqr(m_pos.m_y) +
+                    sqr(m_pos.m_y) +
-                     sqr(m_pos.m_z));            // range, km
+                    sqr(m_pos.m_z));            // range, km
-    m_vel.m_x = -mfactor * m_pos.m_y;            // km / sec
+   m_vel.m_x = -mfactor * m_pos.m_y;            // km / sec
-    m_vel.m_y =  mfactor * m_pos.m_x;
+   m_vel.m_y =  mfactor * m_pos.m_x;
-    m_vel.m_z = 0.0;
+   m_vel.m_z = 0.0;
-    m_vel.m_w = sqrt(sqr(m_vel.m_x) +            // range rate km/sec^2
+   m_vel.m_w = sqrt(sqr(m_vel.m_x) +            // range rate km/sec^2
-                     sqr(m_vel.m_y));
+                    sqr(m_vel.m_y));
  }
  //////////////////////////////////////////////////////////////////////////////
 Line 834 
 cEci::cEci(const cCoordGeo &geo, const c
  // Reference: www.celestrak.com (Dr. TS Kelso)
  cCoordGeo cEci::toGeo()
  {
-    ae2km(); // Vectors must be in kilometer-based units
+   ae2km(); // Vectors must be in kilometer-based units
-    double theta = AcTan(m_pos.m_y, m_pos.m_x);
+   double theta = AcTan(m_pos.m_y, m_pos.m_x);
-    double lon   = fmod(theta - m_date.toGMST(), TWOPI);
+   double lon   = fmod(theta - m_date.toGMST(), TWOPI);
-    if (lon < 0.0)
+   if (lon < 0.0)
-       lon += TWOPI;  // "wrap" negative modulo
+     lon += TWOPI;  // "wrap" negative modulo
-    double r   = sqrt(sqr(m_pos.m_x) + sqr(m_pos.m_y));
+   double r   = sqrt(sqr(m_pos.m_x) + sqr(m_pos.m_y));
-    double e2  = F * (2.0 - F);
+   double e2  = F * (2.0 - F);
-    double lat = AcTan(m_pos.m_z, r);
+   double lat = AcTan(m_pos.m_z, r);
-    const double delta = 1.0e-07;
+   const double delta = 1.0e-07;
-    double phi;
+   double phi;
-    double c;
+   double c;
-    do
+   do
-    {
+     {
        phi = lat;
        c   = 1.0 / sqrt(1.0 - e2 * sqr(sin(phi)));
        lat = AcTan(m_pos.m_z + XKMPER_WGS72 * c * e2 * sin(phi), r);
-    }
+     }
-    while (fabs(lat - phi) > delta);
+   while (fabs(lat - phi) > delta);
-    double alt = r / cos(lat) - XKMPER_WGS72 * c;
+   double alt = r / cos(lat) - XKMPER_WGS72 * c;
-    return cCoordGeo(lat, lon, alt); // radians, radians, kilometers
+   return cCoordGeo(lat, lon, alt); // radians, radians, kilometers
  }
  //////////////////////////////////////////////////////////////////////////////
 Line 869 
 cCoordGeo cEci::toGeo()
  // to kilometer based units.
  void cEci::ae2km()
  {
-    if (UnitsAreAe())
+   if (UnitsAreAe())
-    {
+     {
        MulPos(XKMPER_WGS72 / AE);                       // km
        MulVel((XKMPER_WGS72 / AE) * (MIN_PER_DAY / 86400));  // km/sec
        m_VecUnits = UNITS_KM;
+     }
+ }
+ //
+ // cNoradBase.cpp
+ //
+ // Historical Note:
+ // The equations used here (and in derived classes) to determine satellite
+ // ECI coordinates/velocity come from the December, 1980 NORAD document
+ // "Space Track Report No. 3". The report details 6 orbital models and
+ // provides FORTRAN IV implementations of each. The classes here
+ // implement only two of the orbital models: SGP4 and SDP4. These two models,
+ // one for "near-earth" objects and one for "deep space" objects, are widely
+ // used in satellite tracking software and can produce very accurate results
+ // when used with current NORAD two-line element datum.
+ //
+ // The NORAD FORTRAN IV SGP4/SDP4 implementations were converted to Pascal by
+ // Dr. TS Kelso in 1995. In 1996 these routines were ported in a straight-
+ // forward manner to C++ by Varol Okan. The SGP4/SDP4 classes here were
+ // written by Michael F. Henry in 2002-03 and are a modern C++ re-write of
+ // the work done by Okan. In addition to introducing an object-oriented
+ // architecture, the last residues of the original FORTRAN code (such as
+ // labels and gotos) were eradicated.
+ //
+ // For excellent information on the underlying physics of orbits, visible
+ // satellite observations, current NORAD TLE data, and other related material,
+ // see http://www.celestrak.com which is maintained by Dr. TS Kelso.
+ //
+ // Copyright (c) 2003 Michael F. Henry
+ //
+ // mfh 12/07/2003
+ //
+ //////////////////////////////////////////////////////////////////////////////
+ cNoradBase::cNoradBase(const cOrbit &orbit) :
+   m_Orbit(orbit)
+ {
+   Initialize();
+ }
+ cNoradBase& cNoradBase::operator=(const cNoradBase &b)
+ {
+   // m_Orbit is a "const" member var, so cast away its
+   // "const-ness" in order to complete the assigment.
+   *(const_cast<cOrbit*>(&m_Orbit)) = b.m_Orbit;
+   return *this;
+ }
+ //////////////////////////////////////////////////////////////////////////////
+ // Initialize()
+ // Perform the initialization of member variables, specifically the variables
+ // used by derived-class objects to calculate ECI coordinates.
+ void cNoradBase::Initialize()
+ {
+   // Initialize any variables which are time-independent when
+   // calculating the ECI coordinates of the satellite.
+   m_satInc = m_Orbit.Inclination();
+   m_satEcc = m_Orbit.Eccentricity();
+   m_cosio  = cos(m_satInc);
+   m_theta2 = m_cosio * m_cosio;
+   m_x3thm1 = 3.0 * m_theta2 - 1.0;
+   m_eosq   = m_satEcc * m_satEcc;
+   m_betao2 = 1.0 - m_eosq;
+   m_betao  = sqrt(m_betao2);
+   // The "recovered" semi-minor axis and mean motion.
+   m_aodp  = m_Orbit.SemiMinor();
+   m_xnodp = m_Orbit.mnMotionRec();
+   // For perigee below 156 km, the values of S and QOMS2T are altered.
+   m_perigee = XKMPER_WGS72 * (m_aodp * (1.0 - m_satEcc) - AE);
+   m_s4      = S;
+   m_qoms24  = QOMS2T;
+   if (m_perigee < 156.0)
+     {
+       m_s4 = m_perigee - 78.0;
+       if (m_perigee <= 98.0)
+         {
+           m_s4 = 20.0;
+         }
+       m_qoms24 = pow((120.0 - m_s4) * AE / XKMPER_WGS72, 4.0);
+       m_s4 = m_s4 / XKMPER_WGS72 + AE;
+     }
+   const double pinvsq = 1.0 / (m_aodp * m_aodp * m_betao2 * m_betao2);
+   m_tsi   = 1.0 / (m_aodp - m_s4);
+   m_eta   = m_aodp * m_satEcc * m_tsi;
+   m_etasq = m_eta * m_eta;
+   m_eeta  = m_satEcc * m_eta;
+   const double psisq  = fabs(1.0 - m_etasq);
+   m_coef  = m_qoms24 * pow(m_tsi,4.0);
+   m_coef1 = m_coef / pow(psisq,3.5);
+   const double c2 = m_coef1 * m_xnodp *
+     (m_aodp * (1.0 + 1.5 * m_etasq + m_eeta * (4.0 + m_etasq)) +
+.75 * CK2 * m_tsi / psisq * m_x3thm1 *
+      (8.0 + 3.0 * m_etasq * (8.0 + m_etasq)));
+   m_c1    = m_Orbit.BStar() * c2;
+   m_sinio = sin(m_satInc);
+   const double a3ovk2 = -XJ3 / CK2 * pow(AE,3.0);
+   m_c3     = m_coef * m_tsi * a3ovk2 * m_xnodp * AE * m_sinio / m_satEcc;
+   m_x1mth2 = 1.0 - m_theta2;
+   m_c4     = 2.0 * m_xnodp * m_coef1 * m_aodp * m_betao2 *
+     (m_eta * (2.0 + 0.5 * m_etasq) +
+      m_satEcc * (0.5 + 2.0 * m_etasq) -
+.0 * CK2 * m_tsi / (m_aodp * psisq) *
+      (-3.0 * m_x3thm1 * (1.0 - 2.0 * m_eeta + m_etasq * (1.5 - 0.5 * m_eeta)) +
+.75 * m_x1mth2 *
+       (2.0 * m_etasq - m_eeta * (1.0 + m_etasq)) *
+       cos(2.0 * m_Orbit.ArgPerigee())));
+   const double theta4 = m_theta2 * m_theta2;
+   const double temp1  = 3.0 * CK2 * pinvsq * m_xnodp;
+   const double temp2  = temp1 * CK2 * pinvsq;
+   const double temp3  = 1.25 * CK4 * pinvsq * pinvsq * m_xnodp;
+   m_xmdot = m_xnodp + 0.5 * temp1 * m_betao * m_x3thm1 +
+.0625 * temp2 * m_betao *
+     (13.0 - 78.0 * m_theta2 + 137.0 * theta4);
+   const double x1m5th = 1.0 - 5.0 * m_theta2;
+   m_omgdot = -0.5 * temp1 * x1m5th + 0.0625 * temp2 *
+     (7.0 - 114.0 * m_theta2 + 395.0 * theta4) +
+     temp3 * (3.0 - 36.0 * m_theta2 + 49.0 * theta4);
+   const double xhdot1 = -temp1 * m_cosio;
+   m_xnodot = xhdot1 + (0.5 * temp2 * (4.0 - 19.0 * m_theta2) +
+.0 * temp3 * (3.0 - 7.0 * m_theta2)) * m_cosio;
+   m_xnodcf = 3.5 * m_betao2 * xhdot1 * m_c1;
+   m_t2cof  = 1.5 * m_c1;
+   m_xlcof  = 0.125 * a3ovk2 * m_sinio *
+     (3.0 + 5.0 * m_cosio) / (1.0 + m_cosio);
+   m_aycof  = 0.25 * a3ovk2 * m_sinio;
+   m_x7thm1 = 7.0 * m_theta2 - 1.0;
+ }
+ //////////////////////////////////////////////////////////////////////////////
+ bool cNoradBase::FinalPosition(double incl, double  omega,
+                                double    e, double      a,
+                                double   xl, double  xnode,
+                                double   xn, double tsince,
+                                cEci &eci)
+ {
+   if ((e * e) > 1.0)
+     {
+       // error in satellite data
+       return false;
+     }
+   double beta = sqrt(1.0 - e * e);
+   // Long period periodics
+   double axn  = e * cos(omega);
+   double temp = 1.0 / (a * beta * beta);
+   double xll  = temp * m_xlcof * axn;
+   double aynl = temp * m_aycof;
+   double xlt  = xl + xll;
+   double ayn  = e * sin(omega) + aynl;
+   // Solve Kepler's Equation
+   double capu   = Fmod2p(xlt - xnode);
+   double temp2  = capu;
+   double temp3  = 0.0;
+   double temp4  = 0.0;
+   double temp5  = 0.0;
+   double temp6  = 0.0;
+   double sinepw = 0.0;
+   double cosepw = 0.0;
+   bool   fDone  = false;
+   for (int i = 1; (i <= 10) && !fDone; i++)
+     {
+       sinepw = sin(temp2);
+       cosepw = cos(temp2);
+       temp3 = axn * sinepw;
+       temp4 = ayn * cosepw;
+       temp5 = axn * cosepw;
+       temp6 = ayn * sinepw;
+       double epw = (capu - temp4 + temp3 - temp2) /
+         (1.0 - temp5 - temp6) + temp2;
+       if (fabs(epw - temp2) <= E6A)
+         fDone = true;
+       else
+         temp2 = epw;
+     }
+   // Short period preliminary quantities
+   double ecose = temp5 + temp6;
+   double esine = temp3 - temp4;
+   double elsq  = axn * axn + ayn * ayn;
+   temp  = 1.0 - elsq;
+   double pl = a * temp;
+   double r  = a * (1.0 - ecose);
+   double temp1 = 1.0 / r;
+   double rdot  = XKE * sqrt(a) * esine * temp1;
+   double rfdot = XKE * sqrt(pl) * temp1;
+   temp2 = a * temp1;
+   double betal = sqrt(temp);
+   temp3 = 1.0 / (1.0 + betal);
+   double cosu  = temp2 * (cosepw - axn + ayn * esine * temp3);
+   double sinu  = temp2 * (sinepw - ayn - axn * esine * temp3);
+   double u     = AcTan(sinu, cosu);
+   double sin2u = 2.0 * sinu * cosu;
+   double cos2u = 2.0 * cosu * cosu - 1.0;
+   temp  = 1.0 / pl;
+   temp1 = CK2 * temp;
+   temp2 = temp1 * temp;
+   // Update for short periodics
+   double rk = r * (1.0 - 1.5 * temp2 * betal * m_x3thm1) +
+.5 * temp1 * m_x1mth2 * cos2u;
+   double uk = u - 0.25 * temp2 * m_x7thm1 * sin2u;
+   double xnodek = xnode + 1.5 * temp2 * m_cosio * sin2u;
+   double xinck  = incl + 1.5 * temp2 * m_cosio * m_sinio * cos2u;
+   double rdotk  = rdot - xn * temp1 * m_x1mth2 * sin2u;
+   double rfdotk = rfdot + xn * temp1 * (m_x1mth2 * cos2u + 1.5 * m_x3thm1);
+   // Orientation vectors
+   double sinuk  = sin(uk);
+   double cosuk  = cos(uk);
+   double sinik  = sin(xinck);
+   double cosik  = cos(xinck);
+   double sinnok = sin(xnodek);
+   double cosnok = cos(xnodek);
+   double xmx = -sinnok * cosik;
+   double xmy = cosnok * cosik;
+   double ux  = xmx * sinuk + cosnok * cosuk;
+   double uy  = xmy * sinuk + sinnok * cosuk;
+   double uz  = sinik * sinuk;
+   double vx  = xmx * cosuk - cosnok * sinuk;
+   double vy  = xmy * cosuk - sinnok * sinuk;
+   double vz  = sinik * cosuk;
+   // Position
+   double x = rk * ux;
+   double y = rk * uy;
+   double z = rk * uz;
+   cVector vecPos(x, y, z);
+   // Validate on altitude
+   double altKm = (vecPos.Magnitude() * (XKMPER_WGS72 / AE));
+   if ((altKm < XKMPER_WGS72) || (altKm > (2 * GEOSYNC_ALT)))
+     return false;
+   // Velocity
+   double xdot = rdotk * ux + rfdotk * vx;
+   double ydot = rdotk * uy + rfdotk * vy;
+   double zdot = rdotk * uz + rfdotk * vz;
+   cVector vecVel(xdot, ydot, zdot);
+   cJulian gmt = m_Orbit.Epoch();
+   gmt.addMin(tsince);
+   eci = cEci(vecPos, vecVel, gmt);
+   return true;
+ }
+ //
+ // cNoradSGP4.cpp
+ //
+ // NORAD SGP4 implementation. See historical note in cNoradBase.cpp
+ // Copyright (c) 2003 Michael F. Henry
+ //
+ // mfh 12/07/2003
+ //
+ //////////////////////////////////////////////////////////////////////////////
+ cNoradSGP4::cNoradSGP4(const cOrbit &orbit) :
+   cNoradBase(orbit)
+ {
+   m_c5     = 2.0 * m_coef1 * m_aodp * m_betao2 *
+     (1.0 + 2.75 * (m_etasq + m_eeta) + m_eeta * m_etasq);
+   m_omgcof = m_Orbit.BStar() * m_c3 * cos(m_Orbit.ArgPerigee());
+   m_xmcof  = -TWOTHRD * m_coef * m_Orbit.BStar() * AE / m_eeta;
+   m_delmo  = pow(1.0 + m_eta * cos(m_Orbit.mnAnomaly()), 3.0);
+   m_sinmo  = sin(m_Orbit.mnAnomaly());
+ }
+ //////////////////////////////////////////////////////////////////////////////
+ // getPosition()
+ // This procedure returns the ECI position and velocity for the satellite
+ // in the orbit at the given number of minutes since the TLE epoch time
+ // using the NORAD Simplified General Perturbation 4, near earth orbit
+ // model.
+ //
+ // tsince - Time in minutes since the TLE epoch (GMT).
+ // eci    - ECI object to hold position information.
+ //           To convert the returned ECI position vector to km,
+ //           multiply each component by:
+ //              (XKMPER_WGS72 / AE).
+ //           To convert the returned ECI velocity vector to km/sec,
+ //           multiply each component by:
+ //              (XKMPER_WGS72 / AE) * (MIN_PER_DAY / 86400).
+ bool cNoradSGP4::getPosition(double tsince, cEci &eci)
+ {
+   // For m_perigee less than 220 kilometers, the isimp flag is set and
+   // the equations are truncated to linear variation in sqrt a and
+   // quadratic variation in mean anomaly.  Also, the m_c3 term, the
+   // delta omega term, and the delta m term are dropped.
+   bool isimp = false;
+   if ((m_aodp * (1.0 - m_satEcc) / AE) < (220.0 / XKMPER_WGS72 + AE))
+     {
+       isimp = true;
+     }
+   double d2 = 0.0;
+   double d3 = 0.0;
+   double d4 = 0.0;
+   double t3cof = 0.0;
+   double t4cof = 0.0;
+   double t5cof = 0.0;
+   if (!isimp)
+     {
+       double c1sq = m_c1 * m_c1;
+       d2 = 4.0 * m_aodp * m_tsi * c1sq;
+       double temp = d2 * m_tsi * m_c1 / 3.0;
+       d3 = (17.0 * m_aodp + m_s4) * temp;
+       d4 = 0.5 * temp * m_aodp * m_tsi *
+         (221.0 * m_aodp + 31.0 * m_s4) * m_c1;
+       t3cof = d2 + 2.0 * c1sq;
+       t4cof = 0.25 * (3.0 * d3 + m_c1 * (12.0 * d2 + 10.0 * c1sq));
+       t5cof = 0.2 * (3.0 * d4 + 12.0 * m_c1 * d3 + 6.0 *
+                      d2 * d2 + 15.0 * c1sq * (2.0 * d2 + c1sq));
+     }
+   // Update for secular gravity and atmospheric drag.
+   double xmdf   = m_Orbit.mnAnomaly() + m_xmdot * tsince;
+   double omgadf = m_Orbit.ArgPerigee() + m_omgdot * tsince;
+   double xnoddf = m_Orbit.RAAN() + m_xnodot * tsince;
+   double omega  = omgadf;
+   double xmp    = xmdf;
+   double tsq    = tsince * tsince;
+   double xnode  = xnoddf + m_xnodcf * tsq;
+   double tempa  = 1.0 - m_c1 * tsince;
+   double tempe  = m_Orbit.BStar() * m_c4 * tsince;
+   double templ  = m_t2cof * tsq;
+   if (!isimp)
+     {
+       double delomg = m_omgcof * tsince;
+       double delm = m_xmcof * (pow(1.0 + m_eta * cos(xmdf), 3.0) - m_delmo);
+       double temp = delomg + delm;
+       xmp = xmdf + temp;
+       omega = omgadf - temp;
+       double tcube = tsq * tsince;
+       double tfour = tsince * tcube;
+       tempa = tempa - d2 * tsq - d3 * tcube - d4 * tfour;
+       tempe = tempe + m_Orbit.BStar() * m_c5 * (sin(xmp) - m_sinmo);
+       templ = templ + t3cof * tcube + tfour * (t4cof + tsince * t5cof);
+     }
+   double a  = m_aodp * sqr(tempa);
+   double e  = m_satEcc - tempe;
+   double xl = xmp + omega + xnode + m_xnodp * templ;
+   double xn = XKE / pow(a, 1.5);
+   return FinalPosition(m_satInc, omgadf, e, a, xl, xnode, xn, tsince, eci);
+ }
+ //
+ // cNoradSDP4.cpp
+ //
+ // NORAD SDP4 implementation. See historical note in cNoradBase.cpp
+ // Copyright (c) 2003 Michael F. Henry
+ //
+ // mfh 12/07/2003
+ //
+ const double zns    =  1.19459E-5;     const double c1ss   =  2.9864797E-6;
+ const double zes    =  0.01675;        const double znl    =  1.5835218E-4;
+ const double c1l    =  4.7968065E-7;   const double zel    =  0.05490;
+ const double zcosis =  0.91744867;     const double zsinis =  0.39785416;
+ const double zsings = -0.98088458;     const double zcosgs =  0.1945905;
+ const double q22    =  1.7891679E-6;   const double q31    =  2.1460748E-6;
+ const double q33    =  2.2123015E-7;   const double g22    =  5.7686396;
+ const double g32    =  0.95240898;     const double g44    =  1.8014998;
+ const double g52    =  1.0508330;      const double g54    =  4.4108898;
+ const double root22 =  1.7891679E-6;   const double root32 =  3.7393792E-7;
+ const double root44 =  7.3636953E-9;   const double root52 =  1.1428639E-7;
+ const double root54 =  2.1765803E-9;   const double thdt   =  4.3752691E-3;
+ //////////////////////////////////////////////////////////////////////////////
+ cNoradSDP4::cNoradSDP4(const cOrbit &orbit) :
+    cNoradBase(orbit)
+ {
+    m_sing = sin(m_Orbit.ArgPerigee());
+    m_cosg = cos(m_Orbit.ArgPerigee());
+    dp_savtsn = 0.0;
+    dp_zmos = 0.0;
+    dp_se2 = 0.0;
+    dp_se3 = 0.0;
+    dp_si2 = 0.0;
+    dp_si3 = 0.0;
+    dp_sl2 = 0.0;
+    dp_sl3 = 0.0;
+    dp_sl4 = 0.0;
+    dp_sghs = 0.0;
+    dp_sgh2 = 0.0;
+    dp_sgh3 = 0.0;
+    dp_sgh4 = 0.0;
+    dp_sh2 = 0.0;
+    dp_sh3 = 0.0;
+    dp_zmol = 0.0;
+    dp_ee2 = 0.0;
+    dp_e3 = 0.0;
+    dp_xi2 = 0.0;
+    dp_xi3 = 0.0;
+    dp_xl2 = 0.0;
+    dp_xl3 = 0.0;
+    dp_xl4 = 0.0;
+    dp_xgh2 = 0.0;
+    dp_xgh3 = 0.0;
+    dp_xgh4 = 0.0;
+    dp_xh2 = 0.0;
+    dp_xh3 = 0.0;
+    dp_xqncl = 0.0;
+    dp_thgr = 0.0;
+    dp_omegaq = 0.0;
+    dp_sse = 0.0;
+    dp_ssi = 0.0;
+    dp_ssl = 0.0;
+    dp_ssh = 0.0;
+    dp_ssg = 0.0;
+    dp_d2201 = 0.0;
+    dp_d2211 = 0.0;
+    dp_d3210 = 0.0;
+    dp_d3222 = 0.0;
+    dp_d4410 = 0.0;
+    dp_d4422 = 0.0;
+    dp_d5220 = 0.0;
+    dp_d5232 = 0.0;
+    dp_d5421 = 0.0;
+    dp_d5433 = 0.0;
+    dp_xlamo = 0.0;
+    dp_del1 = 0.0;
+    dp_del2 = 0.0;
+    dp_del3 = 0.0;
+    dp_fasx2 = 0.0;
+    dp_fasx4 = 0.0;
+    dp_fasx6 = 0.0;
+    dp_xfact = 0.0;
+    dp_xli = 0.0;
+    dp_xni = 0.0;
+    dp_atime = 0.0;
+    dp_stepp = 0.0;
+    dp_stepn = 0.0;
+    dp_step2 = 0.0;
+    dp_iresfl = false;
+    dp_isynfl = false;
+ }
+ /////////////////////////////////////////////////////////////////////////////
+ bool cNoradSDP4::DeepInit(double *eosq,  double *sinio,  double *cosio,
+                           double *betao, double *aodp,   double *theta2,
+                           double *sing,  double *cosg,   double *betao2,
+                           double *xmdot, double *omgdot, double *xnodott)
+ {
+    eqsq   = *eosq;
+    siniq  = *sinio;
+    cosiq  = *cosio;
+    rteqsq = *betao;
+    ao     = *aodp;
+    cosq2  = *theta2;
+    sinomo = *sing;
+    cosomo = *cosg;
+    bsq    = *betao2;
+    xlldot = *xmdot;
+    omgdt  = *omgdot;
+    xnodot = *xnodott;
+    // Deep space initialization
+    cJulian jd = m_Orbit.Epoch();
+    dp_thgr = jd.toGMST();
+    double eq   = m_Orbit.Eccentricity();
+    double aqnv = 1.0 / ao;
+    dp_xqncl = m_Orbit.Inclination();
+    double xmao   = m_Orbit.mnAnomaly();
+    double xpidot = omgdt + xnodot;
+    double sinq   = sin(m_Orbit.RAAN());
+    double cosq   = cos(m_Orbit.RAAN());
+    dp_omegaq = m_Orbit.ArgPerigee();
+    // Initialize lunar solar terms
+    double day = jd.FromJan1_12h_1900();
+    if (day != dpi_day)
+    {
+       dpi_day    = day;
+       dpi_xnodce = 4.5236020 - 9.2422029E-4 * day;
+       dpi_stem   = sin(dpi_xnodce);
+       dpi_ctem   = cos(dpi_xnodce);
+       dpi_zcosil = 0.91375164 - 0.03568096 * dpi_ctem;
+       dpi_zsinil = sqrt(1.0 - dpi_zcosil * dpi_zcosil);
+       dpi_zsinhl = 0.089683511 *dpi_stem / dpi_zsinil;
+       dpi_zcoshl = sqrt(1.0 - dpi_zsinhl * dpi_zsinhl);
+       dpi_c      = 4.7199672 + 0.22997150 * day;
+       dpi_gam    = 5.8351514 + 0.0019443680 * day;
+       dp_zmol    = Fmod2p(dpi_c - dpi_gam);
+       dpi_zx     = 0.39785416 * dpi_stem / dpi_zsinil;
+       dpi_zy     = dpi_zcoshl * dpi_ctem + 0.91744867 * dpi_zsinhl * dpi_stem;
+       dpi_zx     = AcTan(dpi_zx,dpi_zy) + dpi_gam - dpi_xnodce;
+       dpi_zcosgl = cos(dpi_zx);
+       dpi_zsingl = sin(dpi_zx);
+       dp_zmos    = 6.2565837 + 0.017201977 * day;
+       dp_zmos    = Fmod2p(dp_zmos);
+    }
+    dp_savtsn = 1.0e20;
+    double zcosg = zcosgs;
+    double zsing = zsings;
+    double zcosi = zcosis;
+    double zsini = zsinis;
+    double zcosh = cosq;
+    double zsinh = sinq;
+    double cc  = c1ss;
+    double zn  = zns;
+    double ze  = zes;
+    double zmo = dp_zmos;
+    double xnoi = 1.0 / m_xnodp;
+    double a1;  double a3;  double a7;  double a8;  double a9;  double a10;
+    double a2;  double a4;  double a5;  double a6;  double x1;  double x2;
+    double x3;  double x4;  double x5;  double x6;  double x7;  double x8;
+    double z31; double z32; double z33; double z1;  double z2;  double z3;
+    double z11; double z12; double z13; double z21; double z22; double z23;
+    double s3;  double s2;  double s4;  double s1;  double s5;  double s6;
+    double s7;
+    double se  = 0.0;  double si = 0.0;  double sl = 0.0;
+    double sgh = 0.0;  double sh = 0.0;
+    // Apply the solar and lunar terms on the first pass, then re-apply the
+    // solar terms again on the second pass.
+    for (int pass = 1; pass <= 2; pass++)
+    {
+       // Do solar terms
+       a1 =  zcosg * zcosh + zsing * zcosi * zsinh;
+       a3 = -zsing * zcosh + zcosg * zcosi * zsinh;
+       a7 = -zcosg * zsinh + zsing * zcosi * zcosh;
+       a8 = zsing * zsini;
+       a9 = zsing * zsinh + zcosg * zcosi * zcosh;
+       a10 = zcosg * zsini;
+       a2 = cosiq * a7 +  siniq * a8;
+       a4 = cosiq * a9 +  siniq * a10;
+       a5 = -siniq * a7 +  cosiq * a8;
+       a6 = -siniq * a9 +  cosiq * a10;
+       x1 = a1 * cosomo + a2 * sinomo;
+       x2 = a3 * cosomo + a4 * sinomo;
+       x3 = -a1 * sinomo + a2 * cosomo;
+       x4 = -a3 * sinomo + a4 * cosomo;
+       x5 = a5 * sinomo;
+       x6 = a6 * sinomo;
+       x7 = a5 * cosomo;
+       x8 = a6 * cosomo;
+       z31 = 12.0 * x1 * x1 - 3.0 * x3 * x3;
+       z32 = 24.0 * x1 * x2 - 6.0 * x3 * x4;
+       z33 = 12.0 * x2 * x2 - 3.0 * x4 * x4;
+       z1 = 3.0 * (a1 * a1 + a2 * a2) + z31 * eqsq;
+       z2 = 6.0 * (a1 * a3 + a2 * a4) + z32 * eqsq;
+       z3 = 3.0 * (a3 * a3 + a4 * a4) + z33 * eqsq;
+       z11 = -6.0 * a1 * a5 + eqsq*(-24.0 * x1 * x7 - 6.0 * x3 * x5);
+       z12 = -6.0 * (a1 * a6 + a3 * a5) +
+             eqsq * (-24.0 * (x2 * x7 + x1 * x8) - 6.0 * (x3 * x6 + x4 * x5));
+       z13 = -6.0 * a3 * a6 + eqsq * (-24.0 * x2 * x8 - 6.0 * x4 * x6);
+       z21 = 6.0 * a2 * a5 + eqsq * (24.0 * x1 * x5 - 6.0 * x3 * x7);
+       z22 = 6.0*(a4 * a5 + a2 * a6) +
+             eqsq * (24.0 * (x2 * x5 + x1 * x6) - 6.0 * (x4 * x7 + x3 * x8));
+       z23 = 6.0 * a4 * a6 + eqsq*(24.0 * x2 * x6 - 6.0 * x4 * x8);
+       z1 = z1 + z1 + bsq * z31;
+       z2 = z2 + z2 + bsq * z32;
+       z3 = z3 + z3 + bsq * z33;
+       s3  = cc * xnoi;
+       s2  = -0.5 * s3/rteqsq;
+       s4  = s3 * rteqsq;
+       s1  = -15.0 * eq * s4;
+       s5  = x1 * x3 + x2 * x4;
+       s6  = x2 * x3 + x1 * x4;
+       s7  = x2 * x4 - x1 * x3;
+       se  = s1 * zn * s5;
+       si  = s2 * zn * (z11 + z13);
+       sl  = -zn * s3 * (z1 + z3 - 14.0 - 6.0 * eqsq);
+       sgh = s4 * zn * (z31 + z33 - 6.0);
+       sh  = -zn * s2 * (z21 + z23);
+       if (dp_xqncl < 5.2359877E-2)
+          sh = 0.0;
+       dp_ee2 = 2.0 * s1 * s6;
+       dp_e3  = 2.0 * s1 * s7;
+       dp_xi2 = 2.0 * s2 * z12;
+       dp_xi3 = 2.0 * s2 * (z13 - z11);
+       dp_xl2 = -2.0 * s3 * z2;
+       dp_xl3 = -2.0 * s3 * (z3 - z1);
+       dp_xl4 = -2.0 * s3 * (-21.0 - 9.0 * eqsq) * ze;
+       dp_xgh2 = 2.0 * s4 * z32;
+       dp_xgh3 = 2.0 * s4 * (z33 - z31);
+       dp_xgh4 = -18.0 * s4 * ze;
+       dp_xh2 = -2.0 * s2 * z22;
+       dp_xh3 = -2.0 * s2 * (z23 - z21);
+       if (pass == 1)
+       {
+          // Do lunar terms
+          dp_sse = se;
+          dp_ssi = si;
+          dp_ssl = sl;
+          dp_ssh = sh / siniq;
+          dp_ssg = sgh - cosiq * dp_ssh;
+          dp_se2 = dp_ee2;
+          dp_si2 = dp_xi2;
+          dp_sl2 = dp_xl2;
+          dp_sgh2 = dp_xgh2;
+          dp_sh2 = dp_xh2;
+          dp_se3 = dp_e3;
+          dp_si3 = dp_xi3;
+          dp_sl3 = dp_xl3;
+          dp_sgh3 = dp_xgh3;
+          dp_sh3 = dp_xh3;
+          dp_sl4 = dp_xl4;
+          dp_sgh4 = dp_xgh4;
+          zcosg = dpi_zcosgl;
+          zsing = dpi_zsingl;
+          zcosi = dpi_zcosil;
+          zsini = dpi_zsinil;
+          zcosh = dpi_zcoshl * cosq + dpi_zsinhl * sinq;
+          zsinh = sinq * dpi_zcoshl - cosq * dpi_zsinhl;
+          zn = znl;
+          cc = c1l;
+          ze = zel;
+          zmo = dp_zmol;
+       }
+    }
+    dp_sse = dp_sse + se;
+    dp_ssi = dp_ssi + si;
+    dp_ssl = dp_ssl + sl;
+    dp_ssg = dp_ssg + sgh - cosiq / siniq * sh;
+    dp_ssh = dp_ssh + sh / siniq;
+    // Geopotential resonance initialization for 12 hour orbits
+    dp_iresfl = false;
+    dp_isynfl = false;
+    bool   bInitOnExit = true;
+    double g310;
+    double f220;
+    double bfact = 0.0;
+    if ((m_xnodp >= 0.0052359877) || (m_xnodp <= 0.0034906585))
+    {
+       if ((m_xnodp < 8.26E-3) || (m_xnodp > 9.24E-3) || (eq < 0.5))
+       {
+          bInitOnExit = false;
+       }
+       else
+       {
+          dp_iresfl = true;
+          double eoc  = eq * eqsq;
+          double g201 = -0.306 - (eq - 0.64) * 0.440;
+          double g211;   double g322;
+          double g410;   double g422;
+          double g520;
+          if (eq <= 0.65)
+          {
+             g211 = 3.616 - 13.247 * eq + 16.290 * eqsq;
+             g310 = -19.302 + 117.390 * eq - 228.419 * eqsq + 156.591 * eoc;
+             g322 = -18.9068 + 109.7927 * eq - 214.6334 * eqsq + 146.5816 * eoc;
+             g410 = -41.122 + 242.694 * eq - 471.094 * eqsq + 313.953 * eoc;
+             g422 = -146.407 + 841.880 * eq - 1629.014 * eqsq + 1083.435 * eoc;
+             g520 = -532.114 + 3017.977 * eq - 5740.0 * eqsq + 3708.276 * eoc;
+          }
+          else
+          {
+             g211 = -72.099 + 331.819 * eq - 508.738 * eqsq + 266.724 * eoc;
+             g310 = -346.844 + 1582.851 * eq - 2415.925 * eqsq + 1246.113 * eoc;
+             g322 = -342.585 + 1554.908 * eq - 2366.899 * eqsq + 1215.972 * eoc;
+             g410 = -1052.797 + 4758.686 * eq - 7193.992 * eqsq + 3651.957 * eoc;
+             g422 = -3581.69 + 16178.11 * eq - 24462.77 * eqsq + 12422.52 * eoc;
+             if (eq <= 0.715)
+                g520 = 1464.74 - 4664.75 * eq + 3763.64 * eqsq;
+             else
+                g520 = -5149.66 + 29936.92 * eq - 54087.36 * eqsq + 31324.56 * eoc;
+          }
+          double g533;
+          double g521;
+          double g532;
+          if (eq < 0.7)
+          {
+             g533 = -919.2277 + 4988.61 * eq - 9064.77 * eqsq + 5542.21 * eoc;
+             g521 = -822.71072 + 4568.6173 * eq - 8491.4146 * eqsq + 5337.524 * eoc;
+             g532 = -853.666 + 4690.25 * eq - 8624.77 * eqsq + 5341.4 * eoc;
+          }
+          else
+          {
+             g533 = -37995.78 + 161616.52 * eq - 229838.2 * eqsq + 109377.94 * eoc;
+             g521 = -51752.104 + 218913.95 * eq - 309468.16 * eqsq + 146349.42 * eoc;
+             g532 = -40023.88 + 170470.89 * eq - 242699.48 * eqsq + 115605.82 * eoc;
+          }
+          double sini2 = siniq * siniq;
+          f220 = 0.75*(1.0 + 2.0 * cosiq + cosq2);
+          double f221 = 1.5 * sini2;
+          double f321 = 1.875 * siniq*(1.0 - 2.0 * cosiq - 3.0 * cosq2);
+          double f322 = -1.875 * siniq*(1.0 + 2.0 * cosiq - 3.0 * cosq2);
+          double f441 = 35.0 * sini2 * f220;
+          double f442 = 39.3750 * sini2 * sini2;
+          double f522 = 9.84375 * siniq*(sini2*(1.0 - 2.0 * cosiq - 5.0 * cosq2) +
+.33333333*(-2.0 + 4.0 * cosiq + 6.0 * cosq2));
+          double f523 = siniq*(4.92187512 * sini2*(-2.0 - 4.0 * cosiq + 10.0 * cosq2) +
+.56250012 * (1.0 + 2.0 * cosiq - 3.0 * cosq2));
+          double f542 = 29.53125 * siniq*(2.0 - 8.0 * cosiq + cosq2 * (-12.0 + 8.0 * cosiq + 10.0 * cosq2));
+          double f543 = 29.53125 * siniq*(-2.0 - 8.0 * cosiq + cosq2 * (12.0 + 8.0 * cosiq - 10.0 * cosq2));
+          double xno2 = m_xnodp * m_xnodp;
+          double ainv2 = aqnv * aqnv;
+          double temp1 = 3.0 * xno2 * ainv2;
+          double temp = temp1 * root22;
+          dp_d2201 = temp * f220 * g201;
+          dp_d2211 = temp * f221 * g211;
+          temp1 = temp1 * aqnv;
+          temp = temp1 * root32;
+          dp_d3210 = temp * f321 * g310;
+          dp_d3222 = temp * f322 * g322;
+          temp1 = temp1 * aqnv;
+          temp = 2.0 * temp1 * root44;
+          dp_d4410 = temp * f441 * g410;
+          dp_d4422 = temp * f442 * g422;
+          temp1 = temp1 * aqnv;
+          temp  = temp1 * root52;
+          dp_d5220 = temp * f522 * g520;
+          dp_d5232 = temp * f523 * g532;
+          temp = 2.0 * temp1 * root54;
+          dp_d5421 = temp * f542 * g521;
+          dp_d5433 = temp * f543 * g533;
+          dp_xlamo = xmao + m_Orbit.RAAN() + m_Orbit.RAAN() - dp_thgr - dp_thgr;
+          bfact = xlldot + xnodot + xnodot - thdt - thdt;
+          bfact = bfact + dp_ssl + dp_ssh + dp_ssh;
+       }
+    }
+    else
+    {
+       // Synchronous resonance terms initialization
+       dp_iresfl = true;
+       dp_isynfl = true;
+       double g200 = 1.0 + eqsq * (-2.5 + 0.8125 * eqsq);
+       g310 = 1.0 + 2.0 * eqsq;
+       double g300 = 1.0 + eqsq * (-6.0 + 6.60937 * eqsq);
+       f220 = 0.75 * (1.0 + cosiq) * (1.0 + cosiq);
+       double f311 = 0.9375 * siniq * siniq * (1.0 + 3 * cosiq) - 0.75 * (1.0 + cosiq);
+       double f330 = 1.0 + cosiq;
+       f330 = 1.875 * f330 * f330 * f330;
+       dp_del1 = 3.0 * m_xnodp * m_xnodp * aqnv * aqnv;
+       dp_del2 = 2.0 * dp_del1 * f220 * g200 * q22;
+       dp_del3 = 3.0 * dp_del1 * f330 * g300 * q33 * aqnv;
+       dp_del1 = dp_del1 * f311 * g310 * q31 * aqnv;
+       dp_fasx2 = 0.13130908;
+       dp_fasx4 = 2.8843198;
+       dp_fasx6 = 0.37448087;
+       dp_xlamo = xmao + m_Orbit.RAAN() + m_Orbit.ArgPerigee() - dp_thgr;
+       bfact = xlldot + xpidot - thdt;
+       bfact = bfact + dp_ssl + dp_ssg + dp_ssh;
+    }
+    if (bInitOnExit)
+    {
+       dp_xfact = bfact - m_xnodp;
+       // Initialize integrator
+       dp_xli = dp_xlamo;
+       dp_xni = m_xnodp;
+       dp_atime = 0.0;
+       dp_stepp = 720.0;
+       dp_stepn = -720.0;
+       dp_step2 = 259200.0;
+    }
+    *eosq   = eqsq;
+    *sinio  = siniq;
+    *cosio  = cosiq;
+    *betao  = rteqsq;
+    *aodp   = ao;
+    *theta2 = cosq2;
+    *sing   = sinomo;
+    *cosg   = cosomo;
+    *betao2 = bsq;
+    *xmdot  = xlldot;
+    *omgdot = omgdt;
+    *xnodott = xnodot;
+    return true;
+ }
+ //////////////////////////////////////////////////////////////////////////////
+ bool cNoradSDP4::DeepCalcDotTerms(double *pxndot, double *pxnddt, double *pxldot)
+ {
+     // Dot terms calculated
+    if (dp_isynfl)
+    {
+       *pxndot = dp_del1 * sin(dp_xli - dp_fasx2) +
+                 dp_del2 * sin(2.0 * (dp_xli - dp_fasx4)) +
+                 dp_del3 * sin(3.0 * (dp_xli - dp_fasx6));
+       *pxnddt = dp_del1 * cos(dp_xli - dp_fasx2) +
+.0 * dp_del2 * cos(2.0 * (dp_xli - dp_fasx4)) +
+.0 * dp_del3 * cos(3.0 * (dp_xli - dp_fasx6));
+    }
+    else
+    {
+       double xomi  = dp_omegaq + omgdt * dp_atime;
+       double x2omi = xomi + xomi;
+       double x2li  = dp_xli + dp_xli;
+       *pxndot = dp_d2201 * sin(x2omi + dp_xli - g22) +
+                 dp_d2211 * sin(dp_xli - g22)         +
+                 dp_d3210 * sin(xomi + dp_xli - g32)  +
+                 dp_d3222 * sin(-xomi + dp_xli - g32) +
+                 dp_d4410 * sin(x2omi + x2li - g44)   +
+                 dp_d4422 * sin(x2li - g44)           +
+                 dp_d5220 * sin(xomi + dp_xli - g52)  +
+                 dp_d5232 * sin(-xomi + dp_xli - g52) +
+                 dp_d5421 * sin(xomi + x2li - g54)    +
+                 dp_d5433 * sin(-xomi + x2li - g54);
+       *pxnddt = dp_d2201 * cos(x2omi + dp_xli - g22) +
+                 dp_d2211 * cos(dp_xli - g22)         +
+                 dp_d3210 * cos(xomi + dp_xli - g32)  +
+                 dp_d3222 * cos(-xomi + dp_xli - g32) +
+                 dp_d5220 * cos(xomi + dp_xli - g52)  +
+                 dp_d5232 * cos(-xomi + dp_xli - g52) +
+.0 * (dp_d4410 * cos(x2omi + x2li - g44) +
+                 dp_d4422 * cos(x2li - g44)         +
+                 dp_d5421 * cos(xomi + x2li - g54)  +
+                 dp_d5433 * cos(-xomi + x2li - g54));
+    }
+    *pxldot = dp_xni + dp_xfact;
+    *pxnddt = (*pxnddt) * (*pxldot);
+    return true;
+ }
+ //////////////////////////////////////////////////////////////////////////////
+ void cNoradSDP4::DeepCalcIntegrator(double *pxndot, double *pxnddt,
+                                     double *pxldot, const double &delt)
+ {
+    DeepCalcDotTerms(pxndot, pxnddt, pxldot);
+    dp_xli = dp_xli + (*pxldot) * delt + (*pxndot) * dp_step2;
+    dp_xni = dp_xni + (*pxndot) * delt + (*pxnddt) * dp_step2;
+    dp_atime = dp_atime + delt;
+ }
+ //////////////////////////////////////////////////////////////////////////////
+ bool cNoradSDP4::DeepSecular(double *xmdf, double *omgadf, double *xnode,
+                              double *emm,  double *xincc,  double *xnn,
+                              double *tsince)
+ {
+    xll    = *xmdf;
+    omgasm = *omgadf;
+    xnodes = *xnode;
+    xn     = *xnn;
+    t      = *tsince;
+    // Deep space secular effects
+    xll    = xll + dp_ssl * t;
+    omgasm = omgasm + dp_ssg * t;
+    xnodes = xnodes + dp_ssh * t;
+    _em    = m_Orbit.Eccentricity() + dp_sse * t;
+    xinc   = m_Orbit.Inclination()  + dp_ssi * t;
+    if (xinc < 0.0)
+    {
+       xinc   = -xinc;
+       xnodes = xnodes + PI;
+       omgasm = omgasm - PI;
+    }
+    double xnddt = 0.0;
+    double xndot = 0.0;
+    double xldot = 0.0;
+    double ft    = 0.0;
+    double delt  = 0.0;
+    bool fDone = false;
+    if (dp_iresfl)
+    {
+       while (!fDone)
+       {
+          if ((dp_atime == 0.0)                ||
+             ((t >= 0.0) && (dp_atime <  0.0)) ||
+             ((t <  0.0) && (dp_atime >= 0.0)))
+          {
+             if (t < 0)
+                delt = dp_stepn;
+             else
+                delt = dp_stepp;
+             // Epoch restart
+             dp_atime = 0.0;
+             dp_xni = m_xnodp;
+             dp_xli = dp_xlamo;
+             fDone = true;
+          }
+          else
+          {
+             if (fabs(t) < fabs(dp_atime))
+             {
+                delt = dp_stepp;
+                if (t >= 0.0)
+                   delt = dp_stepn;
+                DeepCalcIntegrator(&xndot, &xnddt, &xldot, delt);
+             }
+             else
+             {
+                delt = dp_stepn;
+                            delt = dp_stepp;
+                fDone = true;
+             }
+          }
+       }
+       while (fabs(t - dp_atime) >= dp_stepp)
+       {
+          DeepCalcIntegrator(&xndot, &xnddt, &xldot, delt);
+       }
+       ft = t - dp_atime;
+       DeepCalcDotTerms(&xndot, &xnddt, &xldot);
+       xn = dp_xni + xndot * ft + xnddt * ft * ft * 0.5;
+       double xl   = dp_xli + xldot * ft + xndot * ft * ft * 0.5;
+       double temp = -xnodes + dp_thgr + t * thdt;
+       xll = xl - omgasm + temp;
+       if (!dp_isynfl)
+          xll = xl + temp + temp;
+    }
+    *xmdf   = xll;
+    *omgadf = omgasm;
+    *xnode  = xnodes;
+    *emm    = _em;
+    *xincc  = xinc;
+    *xnn    = xn;
+    *tsince = t;
+    return true;
+ }
+ //////////////////////////////////////////////////////////////////////////////
+ bool cNoradSDP4::DeepPeriodics(double *e,      double *xincc,
+                                double *omgadf, double *xnode,
+                                double *xmam)
+ {
+    _em    = *e;
+    xinc   = *xincc;
+    omgasm = *omgadf;
+    xnodes = *xnode;
+    xll    = *xmam;
+    // Lunar-solar periodics
+    double sinis = sin(xinc);
+    double cosis = cos(xinc);
+    double sghs = 0.0;
+    double shs  = 0.0;
+    double sh1  = 0.0;
+    double pe   = 0.0;
+    double pinc = 0.0;
+    double pl   = 0.0;
+    double sghl = 0.0;
+    if (fabs(dp_savtsn - t) >= 30.0)
+    {
+       dp_savtsn = t;
+       double zm = dp_zmos + zns * t;
+       double zf = zm + 2.0 * zes * sin(zm);
+       double sinzf = sin(zf);
+       double f2  = 0.5 * sinzf * sinzf - 0.25;
+       double f3  = -0.5 * sinzf * cos(zf);
+       double ses = dp_se2 * f2 + dp_se3 * f3;
+       double sis = dp_si2 * f2 + dp_si3 * f3;
+       double sls = dp_sl2 * f2 + dp_sl3 * f3 + dp_sl4 * sinzf;
+       sghs = dp_sgh2 * f2 + dp_sgh3 * f3 + dp_sgh4 * sinzf;
+       shs = dp_sh2 * f2 + dp_sh3 * f3;
+       zm = dp_zmol + znl * t;
+       zf = zm + 2.0 * zel * sin(zm);
+       sinzf = sin(zf);
+       f2 = 0.5 * sinzf * sinzf - 0.25;
+       f3 = -0.5 * sinzf * cos(zf);
+       double sel  = dp_ee2 * f2 + dp_e3 * f3;
+       double sil  = dp_xi2 * f2 + dp_xi3 * f3;
+       double sll  = dp_xl2 * f2 + dp_xl3 * f3 + dp_xl4 * sinzf;
+       sghl = dp_xgh2 * f2 + dp_xgh3 * f3 + dp_xgh4 * sinzf;
+       sh1  = dp_xh2 * f2 + dp_xh3 * f3;
+       pe   = ses + sel;
+       pinc = sis + sil;
+       pl   = sls + sll;
+    }
+    double pgh  = sghs + sghl;
+    double ph   = shs + sh1;
+    xinc = xinc + pinc;
+    _em  = _em + pe;
+    if (dp_xqncl >= 0.2)
+    {
+       // Apply periodics directly
+       ph  = ph / siniq;
+       pgh = pgh - cosiq * ph;
+       omgasm = omgasm + pgh;
+       xnodes = xnodes + ph;
+       xll = xll + pl;
     }
+    else
+    {
+       // Apply periodics with Lyddane modification
+       double sinok = sin(xnodes);
+       double cosok = cos(xnodes);
+       double alfdp = sinis * sinok;
+       double betdp = sinis * cosok;
+       double dalf  =  ph * cosok + pinc * cosis * sinok;
+       double dbet  = -ph * sinok + pinc * cosis * cosok;
+       alfdp = alfdp + dalf;
+       betdp = betdp + dbet;
+       double xls = xll + omgasm + cosis * xnodes;
+       double dls = pl + pgh - pinc * xnodes * sinis;
+       xls    = xls + dls;
+       xnodes = AcTan(alfdp, betdp);
+       xll    = xll + pl;
+       omgasm = xls - xll - cos(xinc) * xnodes;
+    }
+    *e      = _em;
+    *xincc  = xinc;
+    *omgadf = omgasm;
+    *xnode  = xnodes;
+    *xmam   = xll;
+    return true;
+ }
+ //////////////////////////////////////////////////////////////////////////////
+ // getPosition()
+ // This procedure returns the ECI position and velocity for the satellite
+ // in the orbit at the given number of minutes since the TLE epoch time
+ // using the NORAD Simplified General Perturbation 4, "deep space" orbit
+ // model.
+ //
+ // tsince  - Time in minutes since the TLE epoch (GMT).
+ // pECI    - pointer to location to store the ECI data.
+ //           To convert the returned ECI position vector to km,
+ //           multiply each component by:
+ //              (XKMPER_WGS72 / AE).
+ //           To convert the returned ECI velocity vector to km/sec,
+ //           multiply each component by:
+ //              (XKMPER_WGS72 / AE) * (MIN_PER_DAY / 86400).
+ bool cNoradSDP4::getPosition(double tsince, cEci &eci)
+ {
+    DeepInit(&m_eosq, &m_sinio, &m_cosio,  &m_betao, &m_aodp,   &m_theta2,
+             &m_sing, &m_cosg,  &m_betao2, &m_xmdot, &m_omgdot, &m_xnodot);
+    // Update for secular gravity and atmospheric drag
+    double xmdf   = m_Orbit.mnAnomaly() + m_xmdot * tsince;
+    double omgadf = m_Orbit.ArgPerigee() + m_omgdot * tsince;
+    double xnoddf = m_Orbit.RAAN() + m_xnodot * tsince;
+    double tsq    = tsince * tsince;
+    double xnode  = xnoddf + m_xnodcf * tsq;
+    double tempa  = 1.0 - m_c1 * tsince;
+    double tempe  = m_Orbit.BStar() * m_c4 * tsince;
+    double templ  = m_t2cof * tsq;
+    double xn     = m_xnodp;
+    double em;
+    double xinc;
+    DeepSecular(&xmdf, &omgadf, &xnode, &em, &xinc, &xn, &tsince);
+    double a    = pow(XKE / xn, TWOTHRD) * sqr(tempa);
+    double e    = em - tempe;
+    double xmam = xmdf + m_xnodp * templ;
+    DeepPeriodics(&e, &xinc, &omgadf, &xnode, &xmam);
+    double xl = xmam + omgadf + xnode;
+    xn = XKE / pow(a, 1.5);
+    return FinalPosition(xinc, omgadf, e, a, xl, xnode, xn, tsince, eci);
+ }
+ // cOrbit.cpp
+ //
+ // Copyright (c) 2002-2003 Michael F. Henry
+ //
+ // mfh 11/15/2003
+ //
+ //////////////////////////////////////////////////////////////////////
+ cOrbit::cOrbit(const cTle &tle) :
+   m_tle(tle),
+   m_pNoradModel(NULL)
+ {
+   m_tle.Initialize();
+   int    epochYear = (int)m_tle.getField(cTle::FLD_EPOCHYEAR);
+   double epochDay  =      m_tle.getField(cTle::FLD_EPOCHDAY );
+   if (epochYear < 57)
+     epochYear += 2000;
+   else
+     epochYear += 1900;
+   m_jdEpoch = cJulian(epochYear, epochDay);
+   m_secPeriod = -1.0;
+   // Recover the original mean motion and semimajor axis from the
+   // input elements.
+   double mm     = mnMotion();
+   double rpmin  = mm * 2 * PI / MIN_PER_DAY;   // rads per minute
+   double a1     = pow(XKE / rpmin, TWOTHRD);
+   double e      = Eccentricity();
+   double i      = Inclination();
+   double temp   = (1.5 * CK2 * (3.0 * sqr(cos(i)) - 1.0) /
+                    pow(1.0 - e * e, 1.5));
+   double delta1 = temp / (a1 * a1);
+   double a0     = a1 *
+     (1.0 - delta1 *
+      ((1.0 / 3.0) + delta1 *
+       (1.0 + 134.0 / 81.0 * delta1)));
+   double delta0 = temp / (a0 * a0);
+   m_mnMotionRec        = rpmin / (1.0 + delta0);
+   m_aeAxisSemiMinorRec = a0 / (1.0 - delta0);
+   m_aeAxisSemiMajorRec = m_aeAxisSemiMinorRec / sqrt(1.0 - (e * e));
+   m_kmPerigeeRec       = XKMPER_WGS72 * (m_aeAxisSemiMajorRec * (1.0 - e) - AE);
+   m_kmApogeeRec        = XKMPER_WGS72 * (m_aeAxisSemiMajorRec * (1.0 + e) - AE);
+   if (2.0 * PI / m_mnMotionRec >= 225.0)
+     {
+       // SDP4 - period >= 225 minutes.
+       m_pNoradModel = new cNoradSDP4(*this);
+     }
+   else
+     {
+       // SGP4 - period < 225 minutes
+       m_pNoradModel = new cNoradSGP4(*this);
+     }
+ }
+ /////////////////////////////////////////////////////////////////////////////
+ cOrbit::~cOrbit()
+ {
+   delete m_pNoradModel;
+ }
+ //////////////////////////////////////////////////////////////////////////////
+ // Return the period in seconds
+ double cOrbit::Period() const
+ {
+   if (m_secPeriod < 0.0)
+     {
+       // Calculate the period using the recovered mean motion.
+       if (m_mnMotionRec == 0)
+         m_secPeriod = 0.0;
+       else
+         m_secPeriod = (2 * PI) / m_mnMotionRec * 60.0;
+     }
+   return m_secPeriod;
+ }
+ //////////////////////////////////////////////////////////////////////////////
+ // Returns elapsed number of seconds from epoch to given time.
+ // Note: "Predicted" TLEs can have epochs in the future.
+ double cOrbit::TPlusEpoch(const cJulian &gmt) const
+ {
+   return gmt.spanSec(Epoch());
+ }
+ //////////////////////////////////////////////////////////////////////////////
+ // Returns the mean anomaly in radians at given GMT.
+ // At epoch, the mean anomaly is given by the elements data.
+ double cOrbit::mnAnomaly(cJulian gmt) const
+ {
+   double span = TPlusEpoch(gmt);
+   double P    = Period();
+   assert(P != 0.0);
+   return fmod(mnAnomaly() + (TWOPI * (span / P)), TWOPI);
  }
+ //////////////////////////////////////////////////////////////////////////////
+ // getPosition()
+ // This procedure returns the ECI position and velocity for the satellite
+ // at "tsince" minutes from the (GMT) TLE epoch. The vectors returned in
+ // the ECI object are kilometer-based.
+ // tsince  - Time in minutes since the TLE epoch (GMT).
+ bool cOrbit::getPosition(double tsince, cEci *pEci) const
+ {
+   bool rc;
+   rc = m_pNoradModel->getPosition(tsince, *pEci);
+   pEci->ae2km();
+   return rc;
+ }
+ //////////////////////////////////////////////////////////////////////////////
+ // SatName()
+ // Return the name of the satellite. If requested, the NORAD number is
+ // appended to the end of the name, i.e., "ISS (ZARYA) #25544".
+ // The name of the satellite with the NORAD number appended is important
+ // because many satellites, especially debris, have the same name and
+ // would otherwise appear to be the same satellite in ouput data.
+ string cOrbit::SatName(bool fAppendId /* = false */) const
+ {
+   string str = m_tle.getName();
+   if (fAppendId)
+     {
+       string strId;
+       m_tle.getField(cTle::FLD_NORADNUM, cTle::U_NATIVE, &strId);
+       str = str + " #" + strId;
+     }
+   return str;
+ }

 Legend:



Removed from v.1.1
 


changed lines


 
Added in v.1.2
 Legend:



Removed from v.1.1
 


changed lines


 
Added in v.1.2
-Removed from v.1.1
+Added in v.1.2

	ViewVC Help
Powered by ViewVC 1.1.23