/[PAMELA software]/yodaUtility/sgp4/cNoradSGP4.cpp
ViewVC logotype

Annotation of /yodaUtility/sgp4/cNoradSGP4.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (hide annotations) (download)
Sun Apr 30 11:08:15 2006 UTC (18 years, 7 months ago) by kusanagi
Branch point for: MAIN
Initial revision

1 kusanagi 1.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    

  ViewVC Help
Powered by ViewVC 1.1.23