| 1 |
kusanagi |
1.1 |
//
|
| 2 |
|
|
// cOrbit.cpp
|
| 3 |
|
|
//
|
| 4 |
|
|
// Copyright (c) 2002-2003 Michael F. Henry
|
| 5 |
|
|
//
|
| 6 |
|
|
// mfh 11/15/2003
|
| 7 |
|
|
//
|
| 8 |
|
|
#include "stdafx.h"
|
| 9 |
|
|
|
| 10 |
|
|
#include "cOrbit.h"
|
| 11 |
|
|
#include "math.h"
|
| 12 |
|
|
#include "time.h"
|
| 13 |
|
|
#include "cVector.h"
|
| 14 |
|
|
#include "cEci.h"
|
| 15 |
|
|
#include "coord.h"
|
| 16 |
|
|
#include "cJulian.h"
|
| 17 |
|
|
#include "cNoradSGP4.h"
|
| 18 |
|
|
#include "cNoradSDP4.h"
|
| 19 |
|
|
|
| 20 |
|
|
//////////////////////////////////////////////////////////////////////
|
| 21 |
|
|
cOrbit::cOrbit(const cTle &tle) :
|
| 22 |
|
|
m_tle(tle),
|
| 23 |
|
|
m_pNoradModel(NULL)
|
| 24 |
|
|
{
|
| 25 |
|
|
m_tle.Initialize();
|
| 26 |
|
|
|
| 27 |
|
|
int epochYear = (int)m_tle.getField(cTle::FLD_EPOCHYEAR);
|
| 28 |
|
|
double epochDay = m_tle.getField(cTle::FLD_EPOCHDAY );
|
| 29 |
|
|
|
| 30 |
|
|
if (epochYear < 57)
|
| 31 |
|
|
epochYear += 2000;
|
| 32 |
|
|
else
|
| 33 |
|
|
epochYear += 1900;
|
| 34 |
|
|
|
| 35 |
|
|
m_jdEpoch = cJulian(epochYear, epochDay);
|
| 36 |
|
|
|
| 37 |
|
|
m_secPeriod = -1.0;
|
| 38 |
|
|
|
| 39 |
|
|
// Recover the original mean motion and semimajor axis from the
|
| 40 |
|
|
// input elements.
|
| 41 |
|
|
double mm = mnMotion();
|
| 42 |
|
|
double rpmin = mm * 2 * PI / MIN_PER_DAY; // rads per minute
|
| 43 |
|
|
|
| 44 |
|
|
double a1 = pow(XKE / rpmin, TWOTHRD);
|
| 45 |
|
|
double e = Eccentricity();
|
| 46 |
|
|
double i = Inclination();
|
| 47 |
|
|
double temp = (1.5 * CK2 * (3.0 * sqr(cos(i)) - 1.0) /
|
| 48 |
|
|
pow(1.0 - e * e, 1.5));
|
| 49 |
|
|
double delta1 = temp / (a1 * a1);
|
| 50 |
|
|
double a0 = a1 *
|
| 51 |
|
|
(1.0 - delta1 *
|
| 52 |
|
|
((1.0 / 3.0) + delta1 *
|
| 53 |
|
|
(1.0 + 134.0 / 81.0 * delta1)));
|
| 54 |
|
|
|
| 55 |
|
|
double delta0 = temp / (a0 * a0);
|
| 56 |
|
|
|
| 57 |
|
|
m_mnMotionRec = rpmin / (1.0 + delta0);
|
| 58 |
|
|
m_aeAxisSemiMinorRec = a0 / (1.0 - delta0);
|
| 59 |
|
|
m_aeAxisSemiMajorRec = m_aeAxisSemiMinorRec / sqrt(1.0 - (e * e));
|
| 60 |
|
|
m_kmPerigeeRec = XKMPER_WGS72 * (m_aeAxisSemiMajorRec * (1.0 - e) - AE);
|
| 61 |
|
|
m_kmApogeeRec = XKMPER_WGS72 * (m_aeAxisSemiMajorRec * (1.0 + e) - AE);
|
| 62 |
|
|
|
| 63 |
|
|
if (2.0 * PI / m_mnMotionRec >= 225.0)
|
| 64 |
|
|
{
|
| 65 |
|
|
// SDP4 - period >= 225 minutes.
|
| 66 |
|
|
m_pNoradModel = new cNoradSDP4(*this);
|
| 67 |
|
|
}
|
| 68 |
|
|
else
|
| 69 |
|
|
{
|
| 70 |
|
|
// SGP4 - period < 225 minutes
|
| 71 |
|
|
m_pNoradModel = new cNoradSGP4(*this);
|
| 72 |
|
|
}
|
| 73 |
|
|
}
|
| 74 |
|
|
|
| 75 |
|
|
/////////////////////////////////////////////////////////////////////////////
|
| 76 |
|
|
cOrbit::~cOrbit()
|
| 77 |
|
|
{
|
| 78 |
|
|
delete m_pNoradModel;
|
| 79 |
|
|
}
|
| 80 |
|
|
|
| 81 |
|
|
//////////////////////////////////////////////////////////////////////////////
|
| 82 |
|
|
// Return the period in seconds
|
| 83 |
|
|
double cOrbit::Period() const
|
| 84 |
|
|
{
|
| 85 |
|
|
if (m_secPeriod < 0.0)
|
| 86 |
|
|
{
|
| 87 |
|
|
// Calculate the period using the recovered mean motion.
|
| 88 |
|
|
if (m_mnMotionRec == 0)
|
| 89 |
|
|
m_secPeriod = 0.0;
|
| 90 |
|
|
else
|
| 91 |
|
|
m_secPeriod = (2 * PI) / m_mnMotionRec * 60.0;
|
| 92 |
|
|
}
|
| 93 |
|
|
|
| 94 |
|
|
return m_secPeriod;
|
| 95 |
|
|
}
|
| 96 |
|
|
|
| 97 |
|
|
//////////////////////////////////////////////////////////////////////////////
|
| 98 |
|
|
// Returns elapsed number of seconds from epoch to given time.
|
| 99 |
|
|
// Note: "Predicted" TLEs can have epochs in the future.
|
| 100 |
|
|
double cOrbit::TPlusEpoch(const cJulian &gmt) const
|
| 101 |
|
|
{
|
| 102 |
|
|
return gmt.spanSec(Epoch());
|
| 103 |
|
|
}
|
| 104 |
|
|
|
| 105 |
|
|
//////////////////////////////////////////////////////////////////////////////
|
| 106 |
|
|
// Returns the mean anomaly in radians at given GMT.
|
| 107 |
|
|
// At epoch, the mean anomaly is given by the elements data.
|
| 108 |
|
|
double cOrbit::mnAnomaly(cJulian gmt) const
|
| 109 |
|
|
{
|
| 110 |
|
|
double span = TPlusEpoch(gmt);
|
| 111 |
|
|
double P = Period();
|
| 112 |
|
|
|
| 113 |
|
|
assert(P != 0.0);
|
| 114 |
|
|
|
| 115 |
|
|
return fmod(mnAnomaly() + (TWOPI * (span / P)), TWOPI);
|
| 116 |
|
|
}
|
| 117 |
|
|
|
| 118 |
|
|
//////////////////////////////////////////////////////////////////////////////
|
| 119 |
|
|
// getPosition()
|
| 120 |
|
|
// This procedure returns the ECI position and velocity for the satellite
|
| 121 |
|
|
// at "tsince" minutes from the (GMT) TLE epoch. The vectors returned in
|
| 122 |
|
|
// the ECI object are kilometer-based.
|
| 123 |
|
|
// tsince - Time in minutes since the TLE epoch (GMT).
|
| 124 |
|
|
bool cOrbit::getPosition(double tsince, cEci *pEci) const
|
| 125 |
|
|
{
|
| 126 |
|
|
bool rc;
|
| 127 |
|
|
|
| 128 |
|
|
rc = m_pNoradModel->getPosition(tsince, *pEci);
|
| 129 |
|
|
|
| 130 |
|
|
pEci->ae2km();
|
| 131 |
|
|
|
| 132 |
|
|
return rc;
|
| 133 |
|
|
}
|
| 134 |
|
|
|
| 135 |
|
|
//////////////////////////////////////////////////////////////////////////////
|
| 136 |
|
|
// SatName()
|
| 137 |
|
|
// Return the name of the satellite. If requested, the NORAD number is
|
| 138 |
|
|
// appended to the end of the name, i.e., "ISS (ZARYA) #25544".
|
| 139 |
|
|
// The name of the satellite with the NORAD number appended is important
|
| 140 |
|
|
// because many satellites, especially debris, have the same name and
|
| 141 |
|
|
// would otherwise appear to be the same satellite in ouput data.
|
| 142 |
|
|
string cOrbit::SatName(bool fAppendId /* = false */) const
|
| 143 |
|
|
{
|
| 144 |
|
|
string str = m_tle.getName();
|
| 145 |
|
|
|
| 146 |
|
|
if (fAppendId)
|
| 147 |
|
|
{
|
| 148 |
|
|
string strId;
|
| 149 |
|
|
|
| 150 |
|
|
m_tle.getField(cTle::FLD_NORADNUM, cTle::U_NATIVE, &strId);
|
| 151 |
|
|
str = str + " #" + strId;
|
| 152 |
|
|
}
|
| 153 |
|
|
|
| 154 |
|
|
return str;
|
| 155 |
|
|
}
|
| 156 |
|
|
|
Parent Directory
| 
Revision Log