Correct date set in SGP4. Added original GMST calculation to Julian.cpp

feature/19
Daniel Warner 2011-12-11 22:33:57 +00:00
parent 9c936a279d
commit f659c173bf
3 changed files with 27 additions and 11 deletions

11
Eci.cpp
View File

@ -1,5 +1,7 @@
#include "Eci.h" #include "Eci.h"
#include "Globals.h"
Eci::Eci(const Julian &date, const CoordGeodetic &geo) Eci::Eci(const Julian &date, const CoordGeodetic &geo)
: date_(date) { : date_(date) {
@ -59,11 +61,12 @@ Eci::~Eci(void) {
CoordGeodetic Eci::ToGeodetic() const { CoordGeodetic Eci::ToGeodetic() const {
const double theta = AcTan(position_.y, position_.x); const double theta = AcTan(position_.y, position_.x);
/*
* changes lon to 0>= and <360 // 0 >= lon < 360
* const double lon = Globals::Fmod2p(theta - date_.ToGreenwichSiderealTime()); // const double lon = Fmod2p(theta - date_.ToGreenwichSiderealTime());
*/ // 180 >= lon < 180
const double lon = fmod(theta - date_.ToGreenwichSiderealTime(), kTWOPI); const double lon = fmod(theta - date_.ToGreenwichSiderealTime(), kTWOPI);
const double r = sqrt((position_.x * position_.x) + (position_.y * position_.y)); const double r = sqrt((position_.x * position_.x) + (position_.y * position_.y));
static const double e2 = kF * (2.0 - kF); static const double e2 = kF * (2.0 - kF);

View File

@ -247,15 +247,27 @@ time_t Julian::ToTime() const {
double Julian::ToGreenwichSiderealTime() const { double Julian::ToGreenwichSiderealTime() const {
#if 0 #if 0
double theta; const double UT = fmod(jul + 0.5, 1.0);
double tut1; const double TU = (jul - 2451545.0 - UT) / 36525.0;
double GMST = 24110.54841 + TU *
(8640184.812866 + TU * (0.093104 - TU * 6.2e-06));
GMST = fmod(GMST + SEC_PER_DAY * OMEGA_E * UT, SEC_PER_DAY);
if (GMST < 0.0)
GMST += SEC_PER_DAY; // "wrap" negative modulo value
return (TWOPI * (GMST / SEC_PER_DAY));
#endif
#if 0
// tut1 = Julian centuries from 2000 Jan. 1 12h UT1 (since J2000 which is 2451545.0) // tut1 = Julian centuries from 2000 Jan. 1 12h UT1 (since J2000 which is 2451545.0)
// a Julian century is 36525 days // a Julian century is 36525 days
tut1 = (date_ - 2451545.0) / 36525.0; const double tut1 = (date_ - 2451545.0) / 36525.0;
// Rotation angle in arcseconds // Rotation angle in arcseconds
theta = 67310.54841 + (876600.0 * 3600.0 + 8640184.812866) * tut1 + 0.093104 * pow(tut1, 2.0) - 0.0000062 * pow(tut1, 3.0); double 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(DegreesToRadians(theta / 240.0), kTWOPI); theta = fmod(DegreesToRadians(theta / 240.0), kTWOPI);
@ -285,7 +297,7 @@ double Julian::ToGreenwichSiderealTime() const {
const double c1p2p = C1 + kTWOPI; const double c1p2p = C1 + kTWOPI;
double gsto = fmod(THGR70 + C1 * ds70 + c1p2p * tfrac + ts70 * ts70 * FK5R, kTWOPI); double gsto = fmod(THGR70 + C1 * ds70 + c1p2p * tfrac + ts70 * ts70 * FK5R, kTWOPI);
if (gsto < 0.0) if (gsto < 0.0)
gsto = gsto + kTWOPI; gsto += kTWOPI;
return gsto; return gsto;
} }

View File

@ -519,8 +519,9 @@ void SGP4::CalculateFinalPositionVelocity(Eci* eci, const double& tsince, const
const double zdot = (rdotk * uz + rfdotk * vz) * kXKMPER / 60.0; const double zdot = (rdotk * uz + rfdotk * vz) * kXKMPER / 60.0;
Vector velocity(xdot, ydot, zdot); Vector velocity(xdot, ydot, zdot);
Timespan diff(tsince); Julian julian(elements_.Epoch());
Julian julian = elements_.Epoch() + diff; julian.AddMin(tsince);
(*eci) = Eci(julian, position, velocity); (*eci) = Eci(julian, position, velocity);
} }