yodaUtility/sgp4/cNoradSGP4.cpp

//
// cNoradSGP4.cpp
//
// NORAD SGP4 implementation. See historical note in cNoradBase.cpp
// Copyright (c) 2003 Michael F. Henry
//
// mfh 12/07/2003
//
#include "stdafx.h"

#include "cNoradSGP4.h"
#include "math.h"
#include "cJulian.h"
#include "cOrbit.h"
#include "cVector.h"
#include "coord.h"

//////////////////////////////////////////////////////////////////////////////
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());
}

cNoradSGP4::~cNoradSGP4(void)
{
}

//////////////////////////////////////////////////////////////////////////////
// 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);
}

1	//
2	// cNoradSGP4.cpp
3	//
4	// NORAD SGP4 implementation. See historical note in cNoradBase.cpp
5	// Copyright (c) 2003 Michael F. Henry
6	//
7	// mfh 12/07/2003
8	//
9	#include "stdafx.h"
10
11	#include "cNoradSGP4.h"
12	#include "math.h"
13	#include "cJulian.h"
14	#include "cOrbit.h"
15	#include "cVector.h"
16	#include "coord.h"
17
18	//////////////////////////////////////////////////////////////////////////////
19	cNoradSGP4::cNoradSGP4(const cOrbit &orbit) :
20	cNoradBase(orbit)
21	{
22	m_c5 = 2.0 * m_coef1 * m_aodp * m_betao2 *
23	(1.0 + 2.75 * (m_etasq + m_eeta) + m_eeta * m_etasq);
24	m_omgcof = m_Orbit.BStar() * m_c3 * cos(m_Orbit.ArgPerigee());
25	m_xmcof = -TWOTHRD * m_coef * m_Orbit.BStar() * AE / m_eeta;
26	m_delmo = pow(1.0 + m_eta * cos(m_Orbit.mnAnomaly()), 3.0);
27	m_sinmo = sin(m_Orbit.mnAnomaly());
28	}
29
30	cNoradSGP4::~cNoradSGP4(void)
31	{
32	}
33
34	//////////////////////////////////////////////////////////////////////////////
35	// getPosition()
36	// This procedure returns the ECI position and velocity for the satellite
37	// in the orbit at the given number of minutes since the TLE epoch time
38	// using the NORAD Simplified General Perturbation 4, near earth orbit
39	// model.
40	//
41	// tsince - Time in minutes since the TLE epoch (GMT).
42	// eci - ECI object to hold position information.
43	// To convert the returned ECI position vector to km,
44	// multiply each component by:
45	// (XKMPER_WGS72 / AE).
46	// To convert the returned ECI velocity vector to km/sec,
47	// multiply each component by:
48	// (XKMPER_WGS72 / AE) * (MIN_PER_DAY / 86400).
49
50	bool cNoradSGP4::getPosition(double tsince, cEci &eci)
51	{
52	// For m_perigee less than 220 kilometers, the isimp flag is set and
53	// the equations are truncated to linear variation in sqrt a and
54	// quadratic variation in mean anomaly. Also, the m_c3 term, the
55	// delta omega term, and the delta m term are dropped.
56	bool isimp = false;
57	if ((m_aodp * (1.0 - m_satEcc) / AE) < (220.0 / XKMPER_WGS72 + AE))
58	{
59	isimp = true;
60	}
61
62	double d2 = 0.0;
63	double d3 = 0.0;
64	double d4 = 0.0;
65
66	double t3cof = 0.0;
67	double t4cof = 0.0;
68	double t5cof = 0.0;
69
70	if (!isimp)
71	{
72	double c1sq = m_c1 * m_c1;
73
74	d2 = 4.0 * m_aodp * m_tsi * c1sq;
75
76	double temp = d2 * m_tsi * m_c1 / 3.0;
77
78	d3 = (17.0 * m_aodp + m_s4) * temp;
79	d4 = 0.5 * temp * m_aodp * m_tsi *
80	(221.0 * m_aodp + 31.0 * m_s4) * m_c1;
81	t3cof = d2 + 2.0 * c1sq;
82	t4cof = 0.25 * (3.0 * d3 + m_c1 * (12.0 * d2 + 10.0 * c1sq));
83	t5cof = 0.2 * (3.0 * d4 + 12.0 * m_c1 * d3 + 6.0 *
84	d2 * d2 + 15.0 * c1sq * (2.0 * d2 + c1sq));
85	}
86
87	// Update for secular gravity and atmospheric drag.
88	double xmdf = m_Orbit.mnAnomaly() + m_xmdot * tsince;
89	double omgadf = m_Orbit.ArgPerigee() + m_omgdot * tsince;
90	double xnoddf = m_Orbit.RAAN() + m_xnodot * tsince;
91	double omega = omgadf;
92	double xmp = xmdf;
93	double tsq = tsince * tsince;
94	double xnode = xnoddf + m_xnodcf * tsq;
95	double tempa = 1.0 - m_c1 * tsince;
96	double tempe = m_Orbit.BStar() * m_c4 * tsince;
97	double templ = m_t2cof * tsq;
98
99	if (!isimp)
100	{
101	double delomg = m_omgcof * tsince;
102	double delm = m_xmcof * (pow(1.0 + m_eta * cos(xmdf), 3.0) - m_delmo);
103	double temp = delomg + delm;
104
105	xmp = xmdf + temp;
106	omega = omgadf - temp;
107
108	double tcube = tsq * tsince;
109	double tfour = tsince * tcube;
110
111	tempa = tempa - d2 * tsq - d3 * tcube - d4 * tfour;
112	tempe = tempe + m_Orbit.BStar() * m_c5 * (sin(xmp) - m_sinmo);
113	templ = templ + t3cof * tcube + tfour * (t4cof + tsince * t5cof);
114	}
115
116	double a = m_aodp * sqr(tempa);
117	double e = m_satEcc - tempe;
118
119
120	double xl = xmp + omega + xnode + m_xnodp * templ;
121	double xn = XKE / pow(a, 1.5);
122
123	return FinalPosition(m_satInc, omgadf, e, a, xl, xnode, xn, tsince, eci);
124	}
125