diff --git a/Eci.cpp b/Eci.cpp index b1b34ad..4618f7b 100644 --- a/Eci.cpp +++ b/Eci.cpp @@ -3,7 +3,7 @@ Eci::Eci(const Julian &date, const CoordGeodetic &geo) : 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 longitude = geo.GetLongitude(); const double altitude = geo.GetAltitude(); @@ -16,9 +16,9 @@ Eci::Eci(const Julian &date, const CoordGeodetic &geo) /* * 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 s = pow(1.0 - Globals::F(), 2.0) * c; - const double achcp = (Globals::XKMPER() * c + altitude) * cos(latitude); + const double c = 1.0 / sqrt(1.0 + kF * (kF - 2.0) * pow(sin(latitude), 2.0)); + const double s = pow(1.0 - kF, 2.0) * c; + const double achcp = (kXKMPER * c + altitude) * cos(latitude); /* * X position in km @@ -28,7 +28,7 @@ Eci::Eci(const Julian &date, const CoordGeodetic &geo) */ position_.SetX(achcp * cos(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()); /* @@ -58,16 +58,16 @@ Eci::~Eci(void) { 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 * 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())); - 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 c = 0.0; int cnt = 0; @@ -76,11 +76,11 @@ CoordGeodetic Eci::ToGeodetic() const { phi = lat; const double sinphi = sin(phi); 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++; } 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); } diff --git a/Globals.cpp b/Globals.cpp index fd798ec..7da6d90 100644 --- a/Globals.cpp +++ b/Globals.cpp @@ -1,41 +1,2 @@ #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); - } -} \ No newline at end of file diff --git a/Globals.h b/Globals.h index 1dde008..3d27c78 100644 --- a/Globals.h +++ b/Globals.h @@ -3,77 +3,100 @@ #include -class Globals { -public: - Globals(void); - ~Globals(void); +const double kAE = 1.0; +const double kQ0 = 120.0; +const double kS0 = 78.0; +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 diff --git a/Julian.cpp b/Julian.cpp index 3985d76..f77e5e8 100644 --- a/Julian.cpp +++ b/Julian.cpp @@ -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); // 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 */ if (theta < 0.0) - theta += Globals::TWOPI(); + theta += kTWOPI; return theta; #endif @@ -299,10 +299,10 @@ double Julian::ToGreenwichSiderealTime() const { /* * find greenwich location at epoch */ - const double c1p2p = C1 + Globals::TWOPI(); - double gsto = fmod(THGR70 + C1 * ds70 + c1p2p * tfrac + ts70 * ts70 * FK5R, Globals::TWOPI()); + const double c1p2p = C1 + kTWOPI; + double gsto = fmod(THGR70 + C1 * ds70 + c1p2p * tfrac + ts70 * ts70 * FK5R, kTWOPI); if (gsto < 0.0) - gsto = gsto + Globals::TWOPI(); + gsto = gsto + kTWOPI; return gsto; } @@ -311,11 +311,11 @@ double Julian::ToGreenwichSiderealTime() const { * Local Mean Sideral Time */ 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 { - + double jdAdj = GetDate() + 0.5; int Z = (int) jdAdj; double F = jdAdj - Z; @@ -323,24 +323,24 @@ void Julian::GetDateTime(struct DateTimeComponents* datetime) const { int A = 0; if (Z < 2299161) { - A = static_cast(Z); + A = static_cast (Z); } else { - int a = static_cast((Z - 1867216.25) / 36524.25); - A = static_cast(Z + 1 + a - static_cast(a / 4)); + int a = static_cast ((Z - 1867216.25) / 36524.25); + A = static_cast (Z + 1 + a - static_cast (a / 4)); } int B = A + 1524; - int C = static_cast((B - 122.1) / 365.25); - int D = static_cast(365.25 * C); - int E = static_cast((B - D) / 30.6001); + int C = static_cast ((B - 122.1) / 365.25); + int D = static_cast (365.25 * C); + int E = static_cast ((B - D) / 30.6001); - datetime->hours = static_cast(F * 24.0); + datetime->hours = static_cast (F * 24.0); F -= datetime->hours / 24.0; - datetime->minutes = static_cast(F * 1440.0); + datetime->minutes = static_cast (F * 1440.0); F -= datetime->minutes / 1440.0; datetime->seconds = F * 86400.0; - datetime->days = B - D - static_cast(30.6001 * E); + datetime->days = B - D - static_cast (30.6001 * E); datetime->months = E < 14 ? E - 1 : E - 13; datetime->years = datetime->months > 2 ? C - 4716 : C - 4715; - } +} diff --git a/Julian.h b/Julian.h index 9c947b2..d81e4d9 100644 --- a/Julian.h +++ b/Julian.h @@ -48,15 +48,15 @@ public: double ToLocalMeanSiderealTime(const double& lon) const; double FromJan1_00h_1900() const { - return date_ - Globals::EPOCH_JAN1_00H_1900(); + return date_ - kEPOCH_JAN1_00H_1900; } double FromJan1_12h_1900() const { - return date_ - Globals::EPOCH_JAN1_12H_1900(); + return date_ - kEPOCH_JAN1_12H_1900; } 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; @@ -70,15 +70,15 @@ public: } void AddHour(double hr) { - date_ += (hr / Globals::HR_PER_DAY()); + date_ += (hr / kHOURS_PER_DAY); } void AddMin(double min) { - date_ += (min / Globals::MIN_PER_DAY()); + date_ += (min / kMINUTES_PER_DAY); } void AddSec(double sec) { - date_ += (sec / Globals::SEC_PER_DAY()); + date_ += (sec / kSECONDS_PER_DAY); } double SpanDay(const Julian& b) const { @@ -86,15 +86,15 @@ public: } 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 { - return SpanDay(b) * Globals::MIN_PER_DAY(); + return SpanDay(b) * kMINUTES_PER_DAY; } double SpanSec(const Julian& b) const { - return SpanDay(b) * Globals::SEC_PER_DAY(); + return SpanDay(b) * kSECONDS_PER_DAY; } static bool IsLeapYear(int y) { diff --git a/Observer.cpp b/Observer.cpp index 72e1a49..743e393 100644 --- a/Observer.cpp +++ b/Observer.cpp @@ -7,8 +7,8 @@ */ Observer::Observer(const double latitude, const double longitude, const double altitude) { - geo_.SetLatitude(Globals::Deg2Rad(latitude)); - geo_.SetLongitude(Globals::Deg2Rad(longitude)); + geo_.SetLatitude(DegreesToRadians(latitude)); + geo_.SetLongitude(DegreesToRadians(longitude)); geo_.SetAltitude(altitude); UpdateObserversEci(Julian()); @@ -69,10 +69,10 @@ CoordTopographic Observer::GetLookAngle(const Eci &eci) { double az = atan(-top_e / top_s); if (top_s > 0.0) - az += Globals::PI(); + az += kPI; if (az < 0.0) - az += 2.0 * Globals::PI(); + az += 2.0 * kPI; double el = asin(top_z / range.GetW()); double rate = range.Dot(range_rate) / range.GetW(); diff --git a/SGP4.cpp b/SGP4.cpp index 3aa6bdf..0e34478 100644 --- a/SGP4.cpp +++ b/SGP4.cpp @@ -6,42 +6,6 @@ #include #include -#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) { Reset(); @@ -65,7 +29,7 @@ void SGP4::SetTle(const Tle& tle) { argument_perigee_ = tle.ArgumentPerigee(false); eccentricity_ = tle.Eccentricity(); inclination_ = tle.Inclination(false); - mean_motion_ = tle.MeanMotion() * TWOPI / Globals::MIN_PER_DAY(); + mean_motion_ = tle.MeanMotion() * kTWOPI / kMINUTES_PER_DAY; bstar_ = tle.BStar(); epoch_ = tle.Epoch(); @@ -76,7 +40,7 @@ void SGP4::SetTle(const Tle& tle) { 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"); } @@ -94,7 +58,7 @@ void SGP4::Initialize() { * recover original mean motion (xnodp) and semimajor axis (aodp) * from input elements */ - const double a1 = pow(XKE / MeanMotion(), TWOTHIRD); + const double a1 = pow(kXKE / MeanMotion(), kTWOTHIRD); common_consts_.cosio = cos(Inclination()); common_consts_.sinio = sin(Inclination()); const double theta2 = common_consts_.cosio * common_consts_.cosio; @@ -102,7 +66,7 @@ void SGP4::Initialize() { const double eosq = Eccentricity() * Eccentricity(); const double betao2 = 1.0 - eosq; 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 a0 = a1 * (1.0 - del1 * (1.0 / 3.0 + del1 * (1.0 + del1 * 134.0 / 81.0))); const double del0 = temp / (a0 * a0); @@ -118,8 +82,8 @@ void SGP4::Initialize() { /* * find perigee and period */ - perigee_ = (RecoveredSemiMajorAxis() * (1.0 - Eccentricity()) - AE) * XKMPER; - period_ = TWOPI / RecoveredMeanMotion(); + perigee_ = (RecoveredSemiMajorAxis() * (1.0 - Eccentricity()) - kAE) * kXKMPER; + period_ = kTWOPI / RecoveredMeanMotion(); @@ -149,15 +113,15 @@ void SGP4::Initialize() { * for perigee below 156km, the values of * s4 and qoms2t are altered */ - double s4 = S; - double qoms24 = QOMS2T; + double s4 = kS; + double qoms24 = kQOMS2T; if (Perigee() < 156.0) { s4 = Perigee() - 78.0; if (Perigee() < 98.0) { s4 = 20.0; } - qoms24 = pow((120.0 - s4) * AE / XKMPER, 4.0); - s4 = s4 / XKMPER + AE; + qoms24 = pow((120.0 - s4) * kAE / kXKMPER, 4.0); + s4 = s4 / kXKMPER + kAE; } /* @@ -173,22 +137,22 @@ void SGP4::Initialize() { const double coef1 = coef / pow(psisq, 3.5); const double c2 = coef1 * RecoveredMeanMotion() * (RecoveredSemiMajorAxis() * (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 * (8.0 + etasq))); 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_.c4 = 2.0 * RecoveredMeanMotion() * coef1 * RecoveredSemiMajorAxis() * betao2 * (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 * (1.5 - 0.5 * eeta)) + 0.75 * common_consts_.x1mth2 * (2.0 * etasq - eeta * (1.0 + etasq)) * cos(2.0 * ArgumentPerigee()))); const double theta4 = theta2 * theta2; - const double temp1 = 3.0 * CK2 * pinvsq * RecoveredMeanMotion(); - const double temp2 = temp1 * CK2 * pinvsq; - const double temp3 = 1.25 * CK4 * pinvsq * pinvsq * RecoveredMeanMotion(); + const double temp1 = 3.0 * kCK2 * pinvsq * RecoveredMeanMotion(); + const double temp2 = temp1 * kCK2 * pinvsq; + const double temp3 = 1.25 * kCK4 * pinvsq * pinvsq * RecoveredMeanMotion(); common_consts_.xmdot = RecoveredMeanMotion() + 0.5 * temp1 * betao * common_consts_.x3thm1 + 0.0625 * temp2 * betao * (13.0 - 78.0 * theta2 + 137.0 * theta4); @@ -222,7 +186,7 @@ void SGP4::Initialize() { double c3 = 0.0; 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(); } @@ -232,7 +196,7 @@ void SGP4::Initialize() { nearspace_consts_.xmcof = 0.0; 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_.sinmo = sin(MeanAnomoly()); @@ -306,7 +270,7 @@ void SGP4::FindPositionSDP4(Eci* eci, double tsince) const { 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; /* @@ -340,8 +304,8 @@ void SGP4::FindPositionSDP4(Eci* eci, double tsince) const { */ if (xincl < 0.0) { xincl = -xincl; - xnode += PI; - omgadf = omgadf - PI; + xnode += kPI; + omgadf = omgadf - kPI; } xl = xmam + omgadf + xnode; @@ -471,7 +435,7 @@ void SGP4::CalculateFinalPositionVelocity(Eci* eci, const double& tsince, const double temp3; 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 */ @@ -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 * convergence problems. */ - const double capu = fmod(xlt - xnode, TWOPI); + const double capu = fmod(xlt - xnode, kTWOPI); double epw = capu; double sinepw = 0.0; @@ -554,8 +518,8 @@ void SGP4::CalculateFinalPositionVelocity(Eci* eci, const double& tsince, const const double r = a * (1.0 - ecose); temp1 = 1.0 / r; - const double rdot = XKE * sqrt(a) * esine * temp1; - const double rfdot = XKE * sqrt(pl) * temp1; + const double rdot = kXKE * sqrt(a) * esine * temp1; + const double rfdot = kXKE * sqrt(pl) * temp1; temp2 = a * temp1; const double betal = sqrt(temp); 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 cos2u = 2.0 * cosu * cosu - 1.0; temp = 1.0 / pl; - temp1 = CK2 * temp; + temp1 = kCK2 * temp; temp2 = temp1 * temp; /* @@ -602,13 +566,13 @@ void SGP4::CalculateFinalPositionVelocity(Eci* eci, const double& tsince, const /* * position and velocity */ - const double x = rk * ux * XKMPER; - const double y = rk * uy * XKMPER; - const double z = rk * uz * XKMPER; + const double x = rk * ux * kXKMPER; + const double y = rk * uy * kXKMPER; + const double z = rk * uz * kXKMPER; Vector position(x, y, z); - const double xdot = (rdotk * ux + rfdotk * vx) * XKMPER / 60.0; - const double ydot = (rdotk * uy + rfdotk * vy) * XKMPER / 60.0; - const double zdot = (rdotk * uz + rfdotk * vz) * XKMPER / 60.0; + const double xdot = (rdotk * ux + rfdotk * vx) * kXKMPER / 60.0; + const double ydot = (rdotk * uy + rfdotk * vy) * kXKMPER / 60.0; + const double zdot = (rdotk * uz + rfdotk * vz) * kXKMPER / 60.0; Vector velocity(xdot, ydot, zdot); 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 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 ctem = cos(xnodce_temp); 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 c = 4.7199672 + 0.22997150 * 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 zy = zcoshl * ctem + 0.91744867 * zsinhl * stem; zx = atan2(zx, zy); - zx = fmod(gam + zx - xnodce, TWOPI); + zx = fmod(gam + zx - xnodce, kTWOPI); const double zcosgl = cos(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 @@ -751,7 +715,7 @@ void SGP4::DeepSpaceInitialize(const double& eosq, const double& sinio, const do * with * 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; } else { 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; 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; } 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; 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; } @@ -1060,9 +1024,9 @@ void SGP4::DeepSpacePeriodics(const double& t, double* em, alfdp += dalf; betdp += dbet; - (*xnodes) = fmod((*xnodes), TWOPI); + (*xnodes) = fmod((*xnodes), kTWOPI); if ((*xnodes) < 0.0) - (*xnodes) += TWOPI; + (*xnodes) += kTWOPI; double xls = (*xll) + (*omgasm) + cosis * (*xnodes); double dls = pl + pgh - pinc * (*xnodes) * sinis; @@ -1075,18 +1039,18 @@ void SGP4::DeepSpacePeriodics(const double& t, double* em, (*xnodes) = atan2(alfdp, betdp); if ((*xnodes) < 0.0) - (*xnodes) += TWOPI; + (*xnodes) += kTWOPI; /* * Get perturbed xnodes in to same quadrant as original. * RAAN is in the range of 0 to 360 degrees * atan2 is in the range of -180 to 180 degrees */ - if (fabs(oldxnodes - (*xnodes)) > PI) { + if (fabs(oldxnodes - (*xnodes)) > kPI) { if ((*xnodes) < oldxnodes) - (*xnodes) += TWOPI; + (*xnodes) += kTWOPI; else - (*xnodes) = (*xnodes) - TWOPI; + (*xnodes) = (*xnodes) - kTWOPI; } (*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; 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) (*xll) = xl + temp - (*omgasm); diff --git a/SolarPosition.cpp b/SolarPosition.cpp index 447784e..c38e1aa 100644 --- a/SolarPosition.cpp +++ b/SolarPosition.cpp @@ -14,21 +14,21 @@ void SolarPosition::FindPosition(const Julian& j, Eci& eci) { const double mjd = j.FromJan1_12h_1900(); const double year = 1900 + mjd / 365.25; - const double T = (mjd + Delta_ET(year) / Globals::SEC_PER_DAY()) / 36525.0; - const double M = Globals::Deg2Rad(Modulus(358.47583 + Modulus(35999.04975 * T, 360.0) + const double T = (mjd + Delta_ET(year) / kSECONDS_PER_DAY) / 36525.0; + const double M = DegreesToRadians(Modulus(358.47583 + Modulus(35999.04975 * 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)); 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)); - const double O = Globals::Deg2Rad(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 nu = Modulus(M + C, Globals::TWOPI()); + const double O = DegreesToRadians(Modulus(259.18 - 1934.142 * T, 360.0)); + const double Lsa = Modulus(L + C - DegreesToRadians(0.00569 - 0.00479 * sin(O)), kTWOPI); + const double nu = Modulus(M + C, kTWOPI); 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)); - R = R * Globals::AU(); + R = R * kAU; Vector solar_position = Vector(R * cos(Lsa), R * sin(Lsa) * cos(eps), @@ -48,5 +48,5 @@ double SolarPosition::Modulus(double arg1, double arg2) 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); } diff --git a/Timespan.cpp b/Timespan.cpp index 210fdf3..a551665 100644 --- a/Timespan.cpp +++ b/Timespan.cpp @@ -25,11 +25,11 @@ double Timespan::GetTotalDays() const { } double Timespan::GetTotalHours() const { - return time_span_ * Globals::HR_PER_DAY(); + return time_span_ * kHOURS_PER_DAY; } double Timespan::GetTotalMinutes() const { - return time_span_ * Globals::MIN_PER_DAY(); + return time_span_ * kMINUTES_PER_DAY; } Timespan& Timespan::operator =(const Timespan& b) { diff --git a/Tle.h b/Tle.h index 2ba7847..279e445 100644 --- a/Tle.h +++ b/Tle.h @@ -59,14 +59,14 @@ public: if (in_degrees) return inclination_; else - return Globals::Deg2Rad(inclination_); + return DegreesToRadians(inclination_); } double RightAscendingNode(const bool in_degrees)const { if (in_degrees) return right_ascending_node_; else - return Globals::Deg2Rad(right_ascending_node_); + return DegreesToRadians(right_ascending_node_); } double Eccentricity()const { @@ -77,14 +77,14 @@ public: if (in_degrees) return argument_perigee_; else - return Globals::Deg2Rad(argument_perigee_); + return DegreesToRadians(argument_perigee_); } double MeanAnomaly(const bool in_degrees)const { if (in_degrees) return mean_anomaly_; else - return Globals::Deg2Rad(mean_anomaly_); + return DegreesToRadians(mean_anomaly_); } double MeanMotion()const {