KJ135
/
sgp4
mirror of http://172.16.0.12:5000/LAMPSARToolSoftware/sgp4.git
4
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Removed #define

feature/19
Daniel Warner 2011-05-20 22:09:23 +01:00
parent 32457f4120
commit f422fce314
10 changed files with 194 additions and 246 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
Eci.cpp
View File

@ -3,7 +3,7 @@
Eci::Eci(const Julian &date, const CoordGeodetic &geo) Eci::Eci(const Julian &date, const CoordGeodetic &geo)
: date_(date) { : date_(date) {
static const double mfactor = Globals::TWOPI() * (Globals::OMEGA_E() / Globals::SEC_PER_DAY()); static const double mfactor = kTWOPI * (kOMEGA_E / kSECONDS_PER_DAY);
const double latitude = geo.GetLatitude(); const double latitude = geo.GetLatitude();
const double longitude = geo.GetLongitude(); const double longitude = geo.GetLongitude();
const double altitude = geo.GetAltitude(); const double altitude = geo.GetAltitude();
@ -16,9 +16,9 @@ Eci::Eci(const Julian &date, const CoordGeodetic &geo)
/* /*
* take into account earth flattening * take into account earth flattening
*/ */
const double c = 1.0 / sqrt(1.0 + Globals::F() * (Globals::F() - 2.0) * pow(sin(latitude), 2.0)); const double c = 1.0 / sqrt(1.0 + kF * (kF - 2.0) * pow(sin(latitude), 2.0));
const double s = pow(1.0 - Globals::F(), 2.0) * c; const double s = pow(1.0 - kF, 2.0) * c;
const double achcp = (Globals::XKMPER() * c + altitude) * cos(latitude); const double achcp = (kXKMPER * c + altitude) * cos(latitude);
/* /*
* X position in km * X position in km
@ -28,7 +28,7 @@ Eci::Eci(const Julian &date, const CoordGeodetic &geo)
*/ */
position_.SetX(achcp * cos(theta)); position_.SetX(achcp * cos(theta));
position_.SetY(achcp * sin(theta)); position_.SetY(achcp * sin(theta));
position_.SetZ((Globals::XKMPER() * s + altitude) * sin(latitude)); position_.SetZ((kXKMPER * s + altitude) * sin(latitude));
position_.SetW(position_.GetMagnitude()); position_.SetW(position_.GetMagnitude());
/* /*
@ -58,16 +58,16 @@ Eci::~Eci(void) {
CoordGeodetic Eci::ToGeodetic() const { CoordGeodetic Eci::ToGeodetic() const {
const double theta = Globals::AcTan(position_.GetY(), position_.GetX()); const double theta = AcTan(position_.GetY(), position_.GetX());
/* /*
* changes lon to 0>= and <360 * changes lon to 0>= and <360
* const double lon = Globals::Fmod2p(theta - date_.ToGreenwichSiderealTime()); * const double lon = Globals::Fmod2p(theta - date_.ToGreenwichSiderealTime());
*/ */
const double lon = fmod(theta - date_.ToGreenwichSiderealTime(), Globals::TWOPI()); const double lon = fmod(theta - date_.ToGreenwichSiderealTime(), kTWOPI);
const double r = sqrt((position_.GetX() * position_.GetX()) + (position_.GetY() * position_.GetY())); const double r = sqrt((position_.GetX() * position_.GetX()) + (position_.GetY() * position_.GetY()));
static const double e2 = Globals::F() * (2.0 - Globals::F()); static const double e2 = kF * (2.0 - kF);
double lat = Globals::AcTan(position_.GetZ(), r); double lat = AcTan(position_.GetZ(), r);
double phi = 0.0; double phi = 0.0;
double c = 0.0; double c = 0.0;
int cnt = 0; int cnt = 0;
@ -76,11 +76,11 @@ CoordGeodetic Eci::ToGeodetic() const {
phi = lat; phi = lat;
const double sinphi = sin(phi); const double sinphi = sin(phi);
c = 1.0 / sqrt(1.0 - e2 * sinphi * sinphi); c = 1.0 / sqrt(1.0 - e2 * sinphi * sinphi);
lat = Globals::AcTan(position_.GetZ() + Globals::XKMPER() * c * e2 * sinphi, r); lat = AcTan(position_.GetZ() + kXKMPER * c * e2 * sinphi, r);
cnt++; cnt++;
} while (fabs(lat - phi) >= 1e-10 && cnt < 10); } while (fabs(lat - phi) >= 1e-10 && cnt < 10);
const double alt = r / cos(lat) - Globals::XKMPER() * c; const double alt = r / cos(lat) - kXKMPER * c;
return CoordGeodetic(lat, lon, alt); return CoordGeodetic(lat, lon, alt);
} }

39
Globals.cpp
View File

@ -1,41 +1,2 @@
#include "Globals.h" #include "Globals.h"
Globals::Globals(void) {
}
Globals::~Globals(void) {
}
double Globals::Fmod2p(const double arg) {
double modu = fmod(arg, TWOPI());
if (modu < 0.0)
modu += TWOPI();
return modu;
}
double Globals::Deg2Rad(const double deg) {
return deg * PI() / 180.0;
}
double Globals::Rad2Deg(const double rad) {
return rad * 180.0 / PI();
}
double Globals::AcTan(const double sinx, const double cosx) {
if (cosx == 0.0) {
if (sinx > 0.0)
return PI() / 2.0;
else
return 3.0 * PI() / 2.0;
} else {
if (cosx > 0.0)
return atan(sinx / cosx);
else
return PI() + atan(sinx / cosx);
}
}

159
Globals.h
View File

@ -3,77 +3,100 @@
#include <cmath> #include <cmath>
class Globals { const double kAE = 1.0;
public: const double kQ0 = 120.0;
Globals(void); const double kS0 = 78.0;
~Globals(void); const double kMU = 398600.8;
const double kXKMPER = 6378.135;
const double kXJ2 = 1.082616e-3;
const double kXJ3 = -2.53881e-6;
const double kXJ4 = -1.65597e-6;
static const double AE() { /*
return 1.0; * alternative XKE
* affects final results
* aiaa-2006-6573
* const double kXKE = 60.0 / sqrt(kXKMPER * kXKMPER * kXKMPER / kMU);
* dundee
* const double kXKE = 7.43669161331734132e-2;
*/
const double kXKE = 60.0 / sqrt(kXKMPER * kXKMPER * kXKMPER / kMU);
const double kCK2 = 0.5 * kXJ2 * kAE * kAE;
const double kCK4 = -0.375 * kXJ4 * kAE * kAE * kAE * kAE;
/*
* alternative QOMS2T
* affects final results
* aiaa-2006-6573
* #define QOMS2T (pow(((Q0 - S0) / XKMPER), 4.0))
* dundee
* #define QOMS2T (1.880279159015270643865e-9)
*/
const double kQOMS2T = pow(((kQ0 - kS0) / kXKMPER), 4.0);
const double kS = kAE * (1.0 + kS0 / kXKMPER);
const double kPI = 3.14159265358979323846264338327950288419716939937510582;
const double kTWOPI = 2.0 * kPI;
const double kTWOTHIRD = 2.0 / 3.0;
const double kTHDT = 4.37526908801129966e-3;
/*
* earth flattening
*/
const double kF = 1.0 / 298.26;
/*
* earth rotation per sideral day
*/
const double kOMEGA_E = 1.00273790934;
const double kAU = 1.49597870691e8;
const double kSECONDS_PER_DAY = 86400.0;
const double kMINUTES_PER_DAY = 1440.0;
const double kHOURS_PER_DAY = 24.0;
// Jan 1.0 1900 = Jan 1 1900 00h UTC
const double kEPOCH_JAN1_00H_1900 = 2415019.5;
// Jan 1.5 1900 = Jan 1 1900 12h UTC
const double kEPOCH_JAN1_12H_1900 = 2415020.0;
// Jan 1.5 2000 = Jan 1 2000 12h UTC
const double kEPOCH_JAN1_12H_2000 = 2451545.0;
inline double Fmod2p(const double arg) {
double modu = fmod(arg, kTWOPI);
if (modu < 0.0)
modu += kTWOPI;
return modu;
}
inline double DegreesToRadians(const double degrees) {
return degrees * kPI / 180.0;
}
inline double RadiansToDegrees(const double radians) {
return radians * 180.0 / kPI;
}
inline double AcTan(const double sinx, const double cosx) {
if (cosx == 0.0) {
if (sinx > 0.0)
return kPI / 2.0;
else
return 3.0 * kPI / 2.0;
} else {
if (cosx > 0.0)
return atan(sinx / cosx);
else
return kPI + atan(sinx / cosx);
} }
}
static const double F() {
/*
* earth flattening
*/
return 1.0 / 298.26;
}
static const double OMEGA_E() {
/*
* earth rotation per sideral day
*/
return 1.00273790934;
}
static const double AU() {
return 1.49597870691e8;
}
static const double XKMPER() {
return 6378.135;
}
static const double PI() {
return 3.14159265358979323846264338327950288419716939937510582;
}
static const double TWOPI() {
return 2.0 * 3.14159265358979323846264338327950288419716939937510582;
}
static const double SEC_PER_DAY() {
return 86400.0;
}
static const double MIN_PER_DAY() {
return 1440.0;
}
static const double HR_PER_DAY() {
return 24.0;
}
static const double EPOCH_JAN1_00H_1900() {
// Jan 1.0 1900 = Jan 1 1900 00h UTC
return 2415019.5;
}
static const double EPOCH_JAN1_12H_1900() {
// Jan 1.5 1900 = Jan 1 1900 12h UTC
return 2415020.0;
}
static const double EPOCH_JAN1_12H_2000() {
// Jan 1.5 2000 = Jan 1 2000 12h UTC
return 2451545.0;
}
static double Fmod2p(const double arg);
static double Deg2Rad(const double deg);
static double Rad2Deg(const double rad);
static double AcTan(const double sinx, const double cosx);
};
#endif #endif

32
Julian.cpp
View File

@ -275,13 +275,13 @@ double Julian::ToGreenwichSiderealTime() const {
theta = 67310.54841 + (876600.0 * 3600.0 + 8640184.812866) * tut1 + 0.093104 * pow(tut1, 2.0) - 0.0000062 * pow(tut1, 3.0); theta = 67310.54841 + (876600.0 * 3600.0 + 8640184.812866) * tut1 + 0.093104 * pow(tut1, 2.0) - 0.0000062 * pow(tut1, 3.0);
// 360.0 / 86400.0 = 1.0 / 240.0 // 360.0 / 86400.0 = 1.0 / 240.0
theta = fmod(Globals::Deg2Rad(theta / 240.0), Globals::TWOPI()); theta = fmod(DegreesToRadians(theta / 240.0), kTWOPI);
/* /*
* check quadrants * check quadrants
*/ */
if (theta < 0.0) if (theta < 0.0)
theta += Globals::TWOPI(); theta += kTWOPI;
return theta; return theta;
#endif #endif
@ -299,10 +299,10 @@ double Julian::ToGreenwichSiderealTime() const {
/* /*
* find greenwich location at epoch * find greenwich location at epoch
*/ */
const double c1p2p = C1 + Globals::TWOPI(); const double c1p2p = C1 + kTWOPI;
double gsto = fmod(THGR70 + C1 * ds70 + c1p2p * tfrac + ts70 * ts70 * FK5R, Globals::TWOPI()); double gsto = fmod(THGR70 + C1 * ds70 + c1p2p * tfrac + ts70 * ts70 * FK5R, kTWOPI);
if (gsto < 0.0) if (gsto < 0.0)
gsto = gsto + Globals::TWOPI(); gsto = gsto + kTWOPI;
return gsto; return gsto;
} }
@ -311,7 +311,7 @@ double Julian::ToGreenwichSiderealTime() const {
* Local Mean Sideral Time * Local Mean Sideral Time
*/ */
double Julian::ToLocalMeanSiderealTime(const double& lon) const { double Julian::ToLocalMeanSiderealTime(const double& lon) const {
return fmod(ToGreenwichSiderealTime() + lon, Globals::TWOPI()); return fmod(ToGreenwichSiderealTime() + lon, kTWOPI);
} }
void Julian::GetDateTime(struct DateTimeComponents* datetime) const { void Julian::GetDateTime(struct DateTimeComponents* datetime) const {
@ -323,24 +323,24 @@ void Julian::GetDateTime(struct DateTimeComponents* datetime) const {
int A = 0; int A = 0;
if (Z < 2299161) { if (Z < 2299161) {
A = static_cast<int>(Z); A = static_cast<int> (Z);
} else { } else {
int a = static_cast<int>((Z - 1867216.25) / 36524.25); int a = static_cast<int> ((Z - 1867216.25) / 36524.25);
A = static_cast<int>(Z + 1 + a - static_cast<int>(a / 4)); A = static_cast<int> (Z + 1 + a - static_cast<int> (a / 4));
} }
int B = A + 1524; int B = A + 1524;
int C = static_cast<int>((B - 122.1) / 365.25); int C = static_cast<int> ((B - 122.1) / 365.25);
int D = static_cast<int>(365.25 * C); int D = static_cast<int> (365.25 * C);
int E = static_cast<int>((B - D) / 30.6001); int E = static_cast<int> ((B - D) / 30.6001);
datetime->hours = static_cast<int>(F * 24.0); datetime->hours = static_cast<int> (F * 24.0);
F -= datetime->hours / 24.0; F -= datetime->hours / 24.0;
datetime->minutes = static_cast<int>(F * 1440.0); datetime->minutes = static_cast<int> (F * 1440.0);
F -= datetime->minutes / 1440.0; F -= datetime->minutes / 1440.0;
datetime->seconds = F * 86400.0; datetime->seconds = F * 86400.0;
datetime->days = B - D - static_cast<int>(30.6001 * E); datetime->days = B - D - static_cast<int> (30.6001 * E);
datetime->months = E < 14 ? E - 1 : E - 13; datetime->months = E < 14 ? E - 1 : E - 13;
datetime->years = datetime->months > 2 ? C - 4716 : C - 4715; datetime->years = datetime->months > 2 ? C - 4716 : C - 4715;
} }

18
Julian.h
View File

@ -48,15 +48,15 @@ public:
double ToLocalMeanSiderealTime(const double& lon) const; double ToLocalMeanSiderealTime(const double& lon) const;
double FromJan1_00h_1900() const { double FromJan1_00h_1900() const {
return date_ - Globals::EPOCH_JAN1_00H_1900(); return date_ - kEPOCH_JAN1_00H_1900;
} }
double FromJan1_12h_1900() const { double FromJan1_12h_1900() const {
return date_ - Globals::EPOCH_JAN1_12H_1900(); return date_ - kEPOCH_JAN1_12H_1900;
} }
double FromJan1_12h_2000() const { double FromJan1_12h_2000() const {
return date_ - Globals::EPOCH_JAN1_12H_2000(); return date_ - kEPOCH_JAN1_12H_2000;
} }
void GetComponent(int& year, int& month, double& dom) const; void GetComponent(int& year, int& month, double& dom) const;
@ -70,15 +70,15 @@ public:
} }
void AddHour(double hr) { void AddHour(double hr) {
date_ += (hr / Globals::HR_PER_DAY()); date_ += (hr / kHOURS_PER_DAY);
} }
void AddMin(double min) { void AddMin(double min) {
date_ += (min / Globals::MIN_PER_DAY()); date_ += (min / kMINUTES_PER_DAY);
} }
void AddSec(double sec) { void AddSec(double sec) {
date_ += (sec / Globals::SEC_PER_DAY()); date_ += (sec / kSECONDS_PER_DAY);
} }
double SpanDay(const Julian& b) const { double SpanDay(const Julian& b) const {
@ -86,15 +86,15 @@ public:
} }
double SpanHour(const Julian& b) const { double SpanHour(const Julian& b) const {
return SpanDay(b) * Globals::HR_PER_DAY(); return SpanDay(b) * kHOURS_PER_DAY;
} }
double SpanMin(const Julian& b) const { double SpanMin(const Julian& b) const {
return SpanDay(b) * Globals::MIN_PER_DAY(); return SpanDay(b) * kMINUTES_PER_DAY;
} }
double SpanSec(const Julian& b) const { double SpanSec(const Julian& b) const {
return SpanDay(b) * Globals::SEC_PER_DAY(); return SpanDay(b) * kSECONDS_PER_DAY;
} }
static bool IsLeapYear(int y) { static bool IsLeapYear(int y) {

8
Observer.cpp
View File

@ -7,8 +7,8 @@
*/ */
Observer::Observer(const double latitude, const double longitude, const double altitude) { Observer::Observer(const double latitude, const double longitude, const double altitude) {
geo_.SetLatitude(Globals::Deg2Rad(latitude)); geo_.SetLatitude(DegreesToRadians(latitude));
geo_.SetLongitude(Globals::Deg2Rad(longitude)); geo_.SetLongitude(DegreesToRadians(longitude));
geo_.SetAltitude(altitude); geo_.SetAltitude(altitude);
UpdateObserversEci(Julian()); UpdateObserversEci(Julian());
@ -69,10 +69,10 @@ CoordTopographic Observer::GetLookAngle(const Eci &eci) {
double az = atan(-top_e / top_s); double az = atan(-top_e / top_s);
if (top_s > 0.0) if (top_s > 0.0)
az += Globals::PI(); az += kPI;
if (az < 0.0) if (az < 0.0)
az += 2.0 * Globals::PI(); az += 2.0 * kPI;
double el = asin(top_z / range.GetW()); double el = asin(top_z / range.GetW());
double rate = range.Dot(range_rate) / range.GetW(); double rate = range.Dot(range_rate) / range.GetW();

128
SGP4.cpp
View File

@ -6,42 +6,6 @@
#include <cmath> #include <cmath>
#include <iomanip> #include <iomanip>
#define AE (1.0)
#define Q0 (120.0)
#define S0 (78.0)
#define MU (398600.8)
#define XKMPER (6378.135)
#define XJ2 (1.082616e-3)
#define XJ3 (-2.53881e-6)
#define XJ4 (-1.65597e-6)
/*
* alternative XKE
* affects final results
* aiaa-2006-6573
* #define XKE (60.0 / sqrt(XKMPER * XKMPER * XKMPER / MU))
* dundee
* #define XKE (7.43669161331734132e-2)
*/
#define XKE (60.0 / sqrt(XKMPER * XKMPER * XKMPER / MU))
#define CK2 (0.5 * XJ2 * AE * AE)
#define CK4 (-0.375 * XJ4 * AE * AE * AE * AE)
/*
* alternative QOMS2T
* affects final results
* aiaa-2006-6573
* #define QOMS2T (pow(((Q0 - S0) / XKMPER), 4.0))
* dundee
* #define QOMS2T (1.880279159015270643865e-9)
*/
#define QOMS2T (pow(((Q0 - S0) / XKMPER), 4.0))
#define S (AE * (1.0 + S0 / XKMPER))
#define PI (3.14159265358979323846264338327950288419716939937510582)
#define TWOPI (2.0 * PI)
#define TWOTHIRD (2.0 / 3.0)
#define THDT (4.37526908801129966e-3)
SGP4::SGP4(void) { SGP4::SGP4(void) {
Reset(); Reset();
@ -65,7 +29,7 @@ void SGP4::SetTle(const Tle& tle) {
argument_perigee_ = tle.ArgumentPerigee(false); argument_perigee_ = tle.ArgumentPerigee(false);
eccentricity_ = tle.Eccentricity(); eccentricity_ = tle.Eccentricity();
inclination_ = tle.Inclination(false); inclination_ = tle.Inclination(false);
mean_motion_ = tle.MeanMotion() * TWOPI / Globals::MIN_PER_DAY(); mean_motion_ = tle.MeanMotion() * kTWOPI / kMINUTES_PER_DAY;
bstar_ = tle.BStar(); bstar_ = tle.BStar();
epoch_ = tle.Epoch(); epoch_ = tle.Epoch();
@ -76,7 +40,7 @@ void SGP4::SetTle(const Tle& tle) {
throw new SatelliteException("Eccentricity out of range"); throw new SatelliteException("Eccentricity out of range");
} }
if (inclination_ < 0.0 || eccentricity_ > PI) { if (inclination_ < 0.0 || eccentricity_ > kPI) {
throw new SatelliteException("Inclination out of range"); throw new SatelliteException("Inclination out of range");
} }
@ -94,7 +58,7 @@ void SGP4::Initialize() {
* recover original mean motion (xnodp) and semimajor axis (aodp) * recover original mean motion (xnodp) and semimajor axis (aodp)
* from input elements * from input elements
*/ */
const double a1 = pow(XKE / MeanMotion(), TWOTHIRD); const double a1 = pow(kXKE / MeanMotion(), kTWOTHIRD);
common_consts_.cosio = cos(Inclination()); common_consts_.cosio = cos(Inclination());
common_consts_.sinio = sin(Inclination()); common_consts_.sinio = sin(Inclination());
const double theta2 = common_consts_.cosio * common_consts_.cosio; const double theta2 = common_consts_.cosio * common_consts_.cosio;
@ -102,7 +66,7 @@ void SGP4::Initialize() {
const double eosq = Eccentricity() * Eccentricity(); const double eosq = Eccentricity() * Eccentricity();
const double betao2 = 1.0 - eosq; const double betao2 = 1.0 - eosq;
const double betao = sqrt(betao2); const double betao = sqrt(betao2);
const double temp = (1.5 * CK2) * common_consts_.x3thm1 / (betao * betao2); const double temp = (1.5 * kCK2) * common_consts_.x3thm1 / (betao * betao2);
const double del1 = temp / (a1 * a1); const double del1 = temp / (a1 * a1);
const double a0 = a1 * (1.0 - del1 * (1.0 / 3.0 + del1 * (1.0 + del1 * 134.0 / 81.0))); const double a0 = a1 * (1.0 - del1 * (1.0 / 3.0 + del1 * (1.0 + del1 * 134.0 / 81.0)));
const double del0 = temp / (a0 * a0); const double del0 = temp / (a0 * a0);
@ -118,8 +82,8 @@ void SGP4::Initialize() {
/* /*
* find perigee and period * find perigee and period
*/ */
perigee_ = (RecoveredSemiMajorAxis() * (1.0 - Eccentricity()) - AE) * XKMPER; perigee_ = (RecoveredSemiMajorAxis() * (1.0 - Eccentricity()) - kAE) * kXKMPER;
period_ = TWOPI / RecoveredMeanMotion(); period_ = kTWOPI / RecoveredMeanMotion();
@ -149,15 +113,15 @@ void SGP4::Initialize() {
* for perigee below 156km, the values of * for perigee below 156km, the values of
* s4 and qoms2t are altered * s4 and qoms2t are altered
*/ */
double s4 = S; double s4 = kS;
double qoms24 = QOMS2T; double qoms24 = kQOMS2T;
if (Perigee() < 156.0) { if (Perigee() < 156.0) {
s4 = Perigee() - 78.0; s4 = Perigee() - 78.0;
if (Perigee() < 98.0) { if (Perigee() < 98.0) {
s4 = 20.0; s4 = 20.0;
} }
qoms24 = pow((120.0 - s4) * AE / XKMPER, 4.0); qoms24 = pow((120.0 - s4) * kAE / kXKMPER, 4.0);
s4 = s4 / XKMPER + AE; s4 = s4 / kXKMPER + kAE;
} }
/* /*
@ -173,22 +137,22 @@ void SGP4::Initialize() {
const double coef1 = coef / pow(psisq, 3.5); const double coef1 = coef / pow(psisq, 3.5);
const double c2 = coef1 * RecoveredMeanMotion() * (RecoveredSemiMajorAxis() * const double c2 = coef1 * RecoveredMeanMotion() * (RecoveredSemiMajorAxis() *
(1.0 + 1.5 * etasq + eeta * (4.0 + etasq)) + (1.0 + 1.5 * etasq + eeta * (4.0 + etasq)) +
0.75 * CK2 * tsi / psisq * 0.75 * kCK2 * tsi / psisq *
common_consts_.x3thm1 * (8.0 + 3.0 * etasq * common_consts_.x3thm1 * (8.0 + 3.0 * etasq *
(8.0 + etasq))); (8.0 + etasq)));
common_consts_.c1 = BStar() * c2; common_consts_.c1 = BStar() * c2;
common_consts_.a3ovk2 = -XJ3 / CK2 * pow(AE, 3.0); common_consts_.a3ovk2 = -kXJ3 / kCK2 * pow(kAE, 3.0);
common_consts_.x1mth2 = 1.0 - theta2; common_consts_.x1mth2 = 1.0 - theta2;
common_consts_.c4 = 2.0 * RecoveredMeanMotion() * coef1 * RecoveredSemiMajorAxis() * betao2 * common_consts_.c4 = 2.0 * RecoveredMeanMotion() * coef1 * RecoveredSemiMajorAxis() * betao2 *
(common_consts_.eta * (2.0 + 0.5 * etasq) + Eccentricity() * (0.5 + 2.0 * etasq) - (common_consts_.eta * (2.0 + 0.5 * etasq) + Eccentricity() * (0.5 + 2.0 * etasq) -
2.0 * CK2 * tsi / (RecoveredSemiMajorAxis() * psisq) * 2.0 * kCK2 * tsi / (RecoveredSemiMajorAxis() * psisq) *
(-3.0 * common_consts_.x3thm1 * (1.0 - 2.0 * eeta + etasq * (-3.0 * common_consts_.x3thm1 * (1.0 - 2.0 * eeta + etasq *
(1.5 - 0.5 * eeta)) + 0.75 * common_consts_.x1mth2 * (2.0 * etasq - eeta * (1.5 - 0.5 * eeta)) + 0.75 * common_consts_.x1mth2 * (2.0 * etasq - eeta *
(1.0 + etasq)) * cos(2.0 * ArgumentPerigee()))); (1.0 + etasq)) * cos(2.0 * ArgumentPerigee())));
const double theta4 = theta2 * theta2; const double theta4 = theta2 * theta2;
const double temp1 = 3.0 * CK2 * pinvsq * RecoveredMeanMotion(); const double temp1 = 3.0 * kCK2 * pinvsq * RecoveredMeanMotion();
const double temp2 = temp1 * CK2 * pinvsq; const double temp2 = temp1 * kCK2 * pinvsq;
const double temp3 = 1.25 * CK4 * pinvsq * pinvsq * RecoveredMeanMotion(); const double temp3 = 1.25 * kCK4 * pinvsq * pinvsq * RecoveredMeanMotion();
common_consts_.xmdot = RecoveredMeanMotion() + 0.5 * temp1 * betao * common_consts_.xmdot = RecoveredMeanMotion() + 0.5 * temp1 * betao *
common_consts_.x3thm1 + 0.0625 * temp2 * betao * common_consts_.x3thm1 + 0.0625 * temp2 * betao *
(13.0 - 78.0 * theta2 + 137.0 * theta4); (13.0 - 78.0 * theta2 + 137.0 * theta4);
@ -222,7 +186,7 @@ void SGP4::Initialize() {
double c3 = 0.0; double c3 = 0.0;
if (Eccentricity() > 1.0e-4) { if (Eccentricity() > 1.0e-4) {
c3 = coef * tsi * common_consts_.a3ovk2 * RecoveredMeanMotion() * AE * c3 = coef * tsi * common_consts_.a3ovk2 * RecoveredMeanMotion() * kAE *
common_consts_.sinio / Eccentricity(); common_consts_.sinio / Eccentricity();
} }
@ -232,7 +196,7 @@ void SGP4::Initialize() {
nearspace_consts_.xmcof = 0.0; nearspace_consts_.xmcof = 0.0;
if (Eccentricity() > 1.0e-4) if (Eccentricity() > 1.0e-4)
nearspace_consts_.xmcof = -TWOTHIRD * coef * BStar() * AE / eeta; nearspace_consts_.xmcof = -kTWOTHIRD * coef * BStar() * kAE / eeta;
nearspace_consts_.delmo = pow(1.0 + common_consts_.eta * (cos(MeanAnomoly())), 3.0); nearspace_consts_.delmo = pow(1.0 + common_consts_.eta * (cos(MeanAnomoly())), 3.0);
nearspace_consts_.sinmo = sin(MeanAnomoly()); nearspace_consts_.sinmo = sin(MeanAnomoly());
@ -306,7 +270,7 @@ void SGP4::FindPositionSDP4(Eci* eci, double tsince) const {
throw new SatelliteException("Error: #2 (xn <= 0.0)"); throw new SatelliteException("Error: #2 (xn <= 0.0)");
} }
a = pow(XKE / xn, TWOTHIRD) * pow(tempa, 2.0); a = pow(kXKE / xn, kTWOTHIRD) * pow(tempa, 2.0);
e = e - tempe; e = e - tempe;
/* /*
@ -340,8 +304,8 @@ void SGP4::FindPositionSDP4(Eci* eci, double tsince) const {
*/ */
if (xincl < 0.0) { if (xincl < 0.0) {
xincl = -xincl; xincl = -xincl;
xnode += PI; xnode += kPI;
omgadf = omgadf - PI; omgadf = omgadf - kPI;
} }
xl = xmam + omgadf + xnode; xl = xmam + omgadf + xnode;
@ -471,7 +435,7 @@ void SGP4::CalculateFinalPositionVelocity(Eci* eci, const double& tsince, const
double temp3; double temp3;
const double beta = sqrt(1.0 - e * e); const double beta = sqrt(1.0 - e * e);
const double xn = XKE / pow(a, 1.5); const double xn = kXKE / pow(a, 1.5);
/* /*
* long period periodics * long period periodics
*/ */
@ -491,7 +455,7 @@ void SGP4::CalculateFinalPositionVelocity(Eci* eci, const double& tsince, const
* - The fmod saves reduction of angle to +/-2pi in sin/cos() and prevents * - The fmod saves reduction of angle to +/-2pi in sin/cos() and prevents
* convergence problems. * convergence problems.
*/ */
const double capu = fmod(xlt - xnode, TWOPI); const double capu = fmod(xlt - xnode, kTWOPI);
double epw = capu; double epw = capu;
double sinepw = 0.0; double sinepw = 0.0;
@ -554,8 +518,8 @@ void SGP4::CalculateFinalPositionVelocity(Eci* eci, const double& tsince, const
const double r = a * (1.0 - ecose); const double r = a * (1.0 - ecose);
temp1 = 1.0 / r; temp1 = 1.0 / r;
const double rdot = XKE * sqrt(a) * esine * temp1; const double rdot = kXKE * sqrt(a) * esine * temp1;
const double rfdot = XKE * sqrt(pl) * temp1; const double rfdot = kXKE * sqrt(pl) * temp1;
temp2 = a * temp1; temp2 = a * temp1;
const double betal = sqrt(temp); const double betal = sqrt(temp);
temp3 = 1.0 / (1.0 + betal); temp3 = 1.0 / (1.0 + betal);
@ -565,7 +529,7 @@ void SGP4::CalculateFinalPositionVelocity(Eci* eci, const double& tsince, const
const double sin2u = 2.0 * sinu * cosu; const double sin2u = 2.0 * sinu * cosu;
const double cos2u = 2.0 * cosu * cosu - 1.0; const double cos2u = 2.0 * cosu * cosu - 1.0;
temp = 1.0 / pl; temp = 1.0 / pl;
temp1 = CK2 * temp; temp1 = kCK2 * temp;
temp2 = temp1 * temp; temp2 = temp1 * temp;
/* /*
@ -602,13 +566,13 @@ void SGP4::CalculateFinalPositionVelocity(Eci* eci, const double& tsince, const
/* /*
* position and velocity * position and velocity
*/ */
const double x = rk * ux * XKMPER; const double x = rk * ux * kXKMPER;
const double y = rk * uy * XKMPER; const double y = rk * uy * kXKMPER;
const double z = rk * uz * XKMPER; const double z = rk * uz * kXKMPER;
Vector position(x, y, z); Vector position(x, y, z);
const double xdot = (rdotk * ux + rfdotk * vx) * XKMPER / 60.0; const double xdot = (rdotk * ux + rfdotk * vx) * kXKMPER / 60.0;
const double ydot = (rdotk * uy + rfdotk * vy) * XKMPER / 60.0; const double ydot = (rdotk * uy + rfdotk * vy) * kXKMPER / 60.0;
const double zdot = (rdotk * uz + rfdotk * vz) * XKMPER / 60.0; const double zdot = (rdotk * uz + rfdotk * vz) * kXKMPER / 60.0;
Vector velocity(xdot, ydot, zdot); Vector velocity(xdot, ydot, zdot);
Julian julian = Epoch(); Julian julian = Epoch();
@ -663,7 +627,7 @@ void SGP4::DeepSpaceInitialize(const double& eosq, const double& sinio, const do
const double jday = Epoch().FromJan1_12h_1900(); const double jday = Epoch().FromJan1_12h_1900();
const double xnodce = 4.5236020 - 9.2422029e-4 * jday; const double xnodce = 4.5236020 - 9.2422029e-4 * jday;
const double xnodce_temp = fmod(xnodce, TWOPI); const double xnodce_temp = fmod(xnodce, kTWOPI);
const double stem = sin(xnodce_temp); const double stem = sin(xnodce_temp);
const double ctem = cos(xnodce_temp); const double ctem = cos(xnodce_temp);
const double zcosil = 0.91375164 - 0.03568096 * ctem; const double zcosil = 0.91375164 - 0.03568096 * ctem;
@ -672,15 +636,15 @@ void SGP4::DeepSpaceInitialize(const double& eosq, const double& sinio, const do
const double zcoshl = sqrt(1.0 - zsinhl * zsinhl); const double zcoshl = sqrt(1.0 - zsinhl * zsinhl);
const double c = 4.7199672 + 0.22997150 * jday; const double c = 4.7199672 + 0.22997150 * jday;
const double gam = 5.8351514 + 0.0019443680 * jday; const double gam = 5.8351514 + 0.0019443680 * jday;
deepspace_consts_.zmol = Globals::Fmod2p(c - gam); deepspace_consts_.zmol = Fmod2p(c - gam);
double zx = 0.39785416 * stem / zsinil; double zx = 0.39785416 * stem / zsinil;
double zy = zcoshl * ctem + 0.91744867 * zsinhl * stem; double zy = zcoshl * ctem + 0.91744867 * zsinhl * stem;
zx = atan2(zx, zy); zx = atan2(zx, zy);
zx = fmod(gam + zx - xnodce, TWOPI); zx = fmod(gam + zx - xnodce, kTWOPI);
const double zcosgl = cos(zx); const double zcosgl = cos(zx);
const double zsingl = sin(zx); const double zsingl = sin(zx);
deepspace_consts_.zmos = Globals::Fmod2p(6.2565837 + 0.017201977 * jday); deepspace_consts_.zmos = Fmod2p(6.2565837 + 0.017201977 * jday);
/* /*
* do solar terms * do solar terms
@ -751,7 +715,7 @@ void SGP4::DeepSpaceInitialize(const double& eosq, const double& sinio, const do
* with * with
* shdq = (-zn * s2 * (z21 + z23)) / sinio * shdq = (-zn * s2 * (z21 + z23)) / sinio
*/ */
if (Inclination() < 5.2359877e-2 || Inclination() > PI - 5.2359877e-2) { if (Inclination() < 5.2359877e-2 || Inclination() > kPI - 5.2359877e-2) {
shdq = 0.0; shdq = 0.0;
} else { } else {
shdq = (-zn * s2 * (z21 + z23)) / sinio; shdq = (-zn * s2 * (z21 + z23)) / sinio;
@ -834,7 +798,7 @@ void SGP4::DeepSpaceInitialize(const double& eosq, const double& sinio, const do
deepspace_consts_.del1 = deepspace_consts_.del1 * f311 * g310 * Q31 * aqnv; deepspace_consts_.del1 = deepspace_consts_.del1 * f311 * g310 * Q31 * aqnv;
integrator_consts_.xlamo = MeanAnomoly() + AscendingNode() + ArgumentPerigee() - deepspace_consts_.gsto; integrator_consts_.xlamo = MeanAnomoly() + AscendingNode() + ArgumentPerigee() - deepspace_consts_.gsto;
bfact = xmdot + xpidot - THDT; bfact = xmdot + xpidot - kTHDT;
bfact += deepspace_consts_.ssl + deepspace_consts_.ssg + deepspace_consts_.ssh; bfact += deepspace_consts_.ssl + deepspace_consts_.ssg + deepspace_consts_.ssh;
} else if (RecoveredMeanMotion() < 8.26e-3 || RecoveredMeanMotion() > 9.24e-3 || Eccentricity() < 0.5) { } else if (RecoveredMeanMotion() < 8.26e-3 || RecoveredMeanMotion() > 9.24e-3 || Eccentricity() < 0.5) {
@ -931,7 +895,7 @@ void SGP4::DeepSpaceInitialize(const double& eosq, const double& sinio, const do
deepspace_consts_.d5433 = temp * f543 * g533; deepspace_consts_.d5433 = temp * f543 * g533;
integrator_consts_.xlamo = MeanAnomoly() + AscendingNode() + AscendingNode() - deepspace_consts_.gsto - deepspace_consts_.gsto; integrator_consts_.xlamo = MeanAnomoly() + AscendingNode() + AscendingNode() - deepspace_consts_.gsto - deepspace_consts_.gsto;
bfact = xmdot + xnodot + xnodot - THDT - THDT; bfact = xmdot + xnodot + xnodot - kTHDT - kTHDT;
bfact = bfact + deepspace_consts_.ssl + deepspace_consts_.ssh + deepspace_consts_.ssh; bfact = bfact + deepspace_consts_.ssl + deepspace_consts_.ssh + deepspace_consts_.ssh;
} }
@ -1060,9 +1024,9 @@ void SGP4::DeepSpacePeriodics(const double& t, double* em,
alfdp += dalf; alfdp += dalf;
betdp += dbet; betdp += dbet;
(*xnodes) = fmod((*xnodes), TWOPI); (*xnodes) = fmod((*xnodes), kTWOPI);
if ((*xnodes) < 0.0) if ((*xnodes) < 0.0)
(*xnodes) += TWOPI; (*xnodes) += kTWOPI;
double xls = (*xll) + (*omgasm) + cosis * (*xnodes); double xls = (*xll) + (*omgasm) + cosis * (*xnodes);
double dls = pl + pgh - pinc * (*xnodes) * sinis; double dls = pl + pgh - pinc * (*xnodes) * sinis;
@ -1075,18 +1039,18 @@ void SGP4::DeepSpacePeriodics(const double& t, double* em,
(*xnodes) = atan2(alfdp, betdp); (*xnodes) = atan2(alfdp, betdp);
if ((*xnodes) < 0.0) if ((*xnodes) < 0.0)
(*xnodes) += TWOPI; (*xnodes) += kTWOPI;
/* /*
* Get perturbed xnodes in to same quadrant as original. * Get perturbed xnodes in to same quadrant as original.
* RAAN is in the range of 0 to 360 degrees * RAAN is in the range of 0 to 360 degrees
* atan2 is in the range of -180 to 180 degrees * atan2 is in the range of -180 to 180 degrees
*/ */
if (fabs(oldxnodes - (*xnodes)) > PI) { if (fabs(oldxnodes - (*xnodes)) > kPI) {
if ((*xnodes) < oldxnodes) if ((*xnodes) < oldxnodes)
(*xnodes) += TWOPI; (*xnodes) += kTWOPI;
else else
(*xnodes) = (*xnodes) - TWOPI; (*xnodes) = (*xnodes) - kTWOPI;
} }
(*xll) += pl; (*xll) += pl;
@ -1172,7 +1136,7 @@ void SGP4::DeepSpaceSecular(const double& t, double* xll, double* omgasm,
*/ */
(*xn) = integrator_params_.xni + integrator_params_.xndot_t * ft + integrator_params_.xnddt_t * ft * ft * 0.5; (*xn) = integrator_params_.xni + integrator_params_.xndot_t * ft + integrator_params_.xnddt_t * ft * ft * 0.5;
const double xl = integrator_params_.xli + integrator_params_.xldot_t * ft + integrator_params_.xndot_t * ft * ft * 0.5; const double xl = integrator_params_.xli + integrator_params_.xldot_t * ft + integrator_params_.xndot_t * ft * ft * 0.5;
const double temp = -(*xnodes) + deepspace_consts_.gsto + t * THDT; const double temp = -(*xnodes) + deepspace_consts_.gsto + t * kTHDT;
if (deepspace_consts_.synchronous_flag) if (deepspace_consts_.synchronous_flag)
(*xll) = xl + temp - (*omgasm); (*xll) = xl + temp - (*omgasm);

20
SolarPosition.cpp
View File

@ -14,21 +14,21 @@ void SolarPosition::FindPosition(const Julian& j, Eci& eci) {
const double mjd = j.FromJan1_12h_1900(); const double mjd = j.FromJan1_12h_1900();
const double year = 1900 + mjd / 365.25; const double year = 1900 + mjd / 365.25;
const double T = (mjd + Delta_ET(year) / Globals::SEC_PER_DAY()) / 36525.0; const double T = (mjd + Delta_ET(year) / kSECONDS_PER_DAY) / 36525.0;
const double M = Globals::Deg2Rad(Modulus(358.47583 + Modulus(35999.04975 * T, 360.0) const double M = DegreesToRadians(Modulus(358.47583 + Modulus(35999.04975 * T, 360.0)
- (0.000150 + 0.0000033 * T) * T * T, 360.0)); - (0.000150 + 0.0000033 * T) * T * T, 360.0));
const double L = Globals::Deg2Rad(Modulus(279.69668 + Modulus(36000.76892 * T, 360.0) const double L = DegreesToRadians(Modulus(279.69668 + Modulus(36000.76892 * T, 360.0)
+ 0.0003025 * T*T, 360.0)); + 0.0003025 * T*T, 360.0));
const double e = 0.01675104 - (0.0000418 + 0.000000126 * T) * T; const double e = 0.01675104 - (0.0000418 + 0.000000126 * T) * T;
const double C = Globals::Deg2Rad((1.919460 - (0.004789 + 0.000014 * T) * T) * sin(M) const double C = DegreesToRadians((1.919460 - (0.004789 + 0.000014 * T) * T) * sin(M)
+ (0.020094 - 0.000100 * T) * sin(2 * M) + 0.000293 * sin(3 * M)); + (0.020094 - 0.000100 * T) * sin(2 * M) + 0.000293 * sin(3 * M));
const double O = Globals::Deg2Rad(Modulus(259.18 - 1934.142 * T, 360.0)); const double O = DegreesToRadians(Modulus(259.18 - 1934.142 * T, 360.0));
const double Lsa = Modulus(L + C - Globals::Deg2Rad(0.00569 - 0.00479 * sin(O)), Globals::TWOPI()); const double Lsa = Modulus(L + C - DegreesToRadians(0.00569 - 0.00479 * sin(O)), kTWOPI);
const double nu = Modulus(M + C, Globals::TWOPI()); const double nu = Modulus(M + C, kTWOPI);
double R = 1.0000002 * (1 - e * e) / (1 + e * cos(nu)); double R = 1.0000002 * (1 - e * e) / (1 + e * cos(nu));
const double eps = Globals::Deg2Rad(23.452294 - (0.0130125 + (0.00000164 - const double eps = DegreesToRadians(23.452294 - (0.0130125 + (0.00000164 -
0.000000503 * T) * T) * T + 0.00256 * cos(O)); 0.000000503 * T) * T) * T + 0.00256 * cos(O));
R = R * Globals::AU(); R = R * kAU;
Vector solar_position = Vector(R * cos(Lsa), Vector solar_position = Vector(R * cos(Lsa),
R * sin(Lsa) * cos(eps), R * sin(Lsa) * cos(eps),
@ -48,5 +48,5 @@ double SolarPosition::Modulus(double arg1, double arg2) const {
double SolarPosition::Delta_ET(double year) const { double SolarPosition::Delta_ET(double year) const {
return 26.465 + 0.747622 * (year - 1950) + 1.886913 * sin(Globals::TWOPI() * (year - 1975) / 33); return 26.465 + 0.747622 * (year - 1950) + 1.886913 * sin(kTWOPI * (year - 1975) / 33);
} }

4
Timespan.cpp
View File

@ -25,11 +25,11 @@ double Timespan::GetTotalDays() const {
} }
double Timespan::GetTotalHours() const { double Timespan::GetTotalHours() const {
return time_span_ * Globals::HR_PER_DAY(); return time_span_ * kHOURS_PER_DAY;
} }
double Timespan::GetTotalMinutes() const { double Timespan::GetTotalMinutes() const {
return time_span_ * Globals::MIN_PER_DAY(); return time_span_ * kMINUTES_PER_DAY;
} }
Timespan& Timespan::operator =(const Timespan& b) { Timespan& Timespan::operator =(const Timespan& b) {

8
Tle.h
View File

@ -59,14 +59,14 @@ public:
if (in_degrees) if (in_degrees)
return inclination_; return inclination_;
else else
return Globals::Deg2Rad(inclination_); return DegreesToRadians(inclination_);
} }
double RightAscendingNode(const bool in_degrees)const { double RightAscendingNode(const bool in_degrees)const {
if (in_degrees) if (in_degrees)
return right_ascending_node_; return right_ascending_node_;
else else
return Globals::Deg2Rad(right_ascending_node_); return DegreesToRadians(right_ascending_node_);
} }
double Eccentricity()const { double Eccentricity()const {
@ -77,14 +77,14 @@ public:
if (in_degrees) if (in_degrees)
return argument_perigee_; return argument_perigee_;
else else
return Globals::Deg2Rad(argument_perigee_); return DegreesToRadians(argument_perigee_);
} }
double MeanAnomaly(const bool in_degrees)const { double MeanAnomaly(const bool in_degrees)const {
if (in_degrees) if (in_degrees)
return mean_anomaly_; return mean_anomaly_;
else else
return Globals::Deg2Rad(mean_anomaly_); return DegreesToRadians(mean_anomaly_);
} }
double MeanMotion()const { double MeanMotion()const {