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

Dos to unix

feature/19
Daniel Warner 2011-12-24 00:26:41 +00:00
parent 3ca8edbc96
commit 84bbd5fa54
8 changed files with 2203 additions and 2203 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

164
Eci.cpp
View File

@ -1,82 +1,82 @@
#include "Eci.h" #include "Eci.h"
#include "Util.h" #include "Util.h"
void Eci::ToEci(const Julian& date, const CoordGeodetic &g) void Eci::ToEci(const Julian& date, const CoordGeodetic &g)
{ {
/* /*
* set date * set date
*/ */
date_ = date; date_ = date;
static const double mfactor = kTWOPI * (kOMEGA_E / kSECONDS_PER_DAY); static const double mfactor = kTWOPI * (kOMEGA_E / kSECONDS_PER_DAY);
/* /*
* Calculate Local Mean Sidereal Time for observers longitude * Calculate Local Mean Sidereal Time for observers longitude
*/ */
const double theta = date_.ToLocalMeanSiderealTime(g.longitude); const double theta = date_.ToLocalMeanSiderealTime(g.longitude);
/* /*
* take into account earth flattening * take into account earth flattening
*/ */
const double c = 1.0 const double c = 1.0
/ sqrt(1.0 + kF * (kF - 2.0) * pow(sin(g.latitude), 2.0)); / sqrt(1.0 + kF * (kF - 2.0) * pow(sin(g.latitude), 2.0));
const double s = pow(1.0 - kF, 2.0) * c; const double s = pow(1.0 - kF, 2.0) * c;
const double achcp = (kXKMPER * c + g.altitude) * cos(g.latitude); const double achcp = (kXKMPER * c + g.altitude) * cos(g.latitude);
/* /*
* X position in km * X position in km
* Y position in km * Y position in km
* Z position in km * Z position in km
* W magnitude in km * W magnitude in km
*/ */
position_.x = achcp * cos(theta); position_.x = achcp * cos(theta);
position_.y = achcp * sin(theta); position_.y = achcp * sin(theta);
position_.z = (kXKMPER * s + g.altitude) * sin(g.latitude); position_.z = (kXKMPER * s + g.altitude) * sin(g.latitude);
position_.w = position_.GetMagnitude(); position_.w = position_.GetMagnitude();
/* /*
* X velocity in km/s * X velocity in km/s
* Y velocity in km/s * Y velocity in km/s
* Z velocity in km/s * Z velocity in km/s
* W magnitude in km/s * W magnitude in km/s
*/ */
velocity_.x = -mfactor * position_.y; velocity_.x = -mfactor * position_.y;
velocity_.y = mfactor * position_.x; velocity_.y = mfactor * position_.x;
velocity_.z = 0.0; velocity_.z = 0.0;
velocity_.w = velocity_.GetMagnitude(); velocity_.w = velocity_.GetMagnitude();
} }
CoordGeodetic Eci::ToGeodetic() const CoordGeodetic Eci::ToGeodetic() const
{ {
const double theta = Util::AcTan(position_.y, position_.x); const double theta = Util::AcTan(position_.y, position_.x);
// 0 >= lon < 360 // 0 >= lon < 360
// const double lon = Fmod2p(theta - date_.ToGreenwichSiderealTime()); // const double lon = Fmod2p(theta - date_.ToGreenwichSiderealTime());
// 180 >= lon < 180 // 180 >= lon < 180
const double lon = Util::WrapNegPosPI(theta - date_.ToGreenwichSiderealTime()); const double lon = Util::WrapNegPosPI(theta - date_.ToGreenwichSiderealTime());
const double r = sqrt((position_.x * position_.x) const double r = sqrt((position_.x * position_.x)
+ (position_.y * position_.y)); + (position_.y * position_.y));
static const double e2 = kF * (2.0 - kF); static const double e2 = kF * (2.0 - kF);
double lat = Util::AcTan(position_.z, r); double lat = Util::AcTan(position_.z, r);
double phi = 0.0; double phi = 0.0;
double c = 0.0; double c = 0.0;
int cnt = 0; int cnt = 0;
do do
{ {
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 = Util::AcTan(position_.z + kXKMPER * c * e2 * sinphi, r); lat = Util::AcTan(position_.z + 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) - kXKMPER * c; const double alt = r / cos(lat) - kXKMPER * c;
return CoordGeodetic(lat, lon, alt, true); return CoordGeodetic(lat, lon, alt, true);
} }

612
Julian.cpp
View File

@ -1,306 +1,306 @@
#include "Globals.h" #include "Globals.h"
#include "Julian.h" #include "Julian.h"
#include "Util.h" #include "Util.h"
#include <cmath> #include <cmath>
#include <ctime> #include <ctime>
#include <cassert> #include <cassert>
#ifdef WIN32 #ifdef WIN32
#include <windows.h> #include <windows.h>
#else #else
#include <sys/time.h> #include <sys/time.h>
#endif #endif
Julian::Julian() Julian::Julian()
{ {
#ifdef WIN32 #ifdef WIN32
SYSTEMTIME st; SYSTEMTIME st;
GetSystemTime(&st); GetSystemTime(&st);
Initialize(st.wYear, Initialize(st.wYear,
st.wMonth, st.wMonth,
st.wDay, st.wDay,
st.wHour, st.wHour,
st.wMinute, st.wMinute,
(double) st.wSecond + (double) st.wMilliseconds / 1000.0); (double) st.wSecond + (double) st.wMilliseconds / 1000.0);
#else #else
struct timeval tv; struct timeval tv;
gettimeofday(&tv, NULL); gettimeofday(&tv, NULL);
struct tm gmt; struct tm gmt;
gmtime_r(&tv.tv_sec, &gmt); gmtime_r(&tv.tv_sec, &gmt);
Initialize(gmt.tm_year + 1900, Initialize(gmt.tm_year + 1900,
gmt.tm_mon + 1, gmt.tm_mon + 1,
gmt.tm_mday, gmt.tm_mday,
gmt.tm_hour, gmt.tm_hour,
gmt.tm_min, gmt.tm_min,
(double) gmt.tm_sec + (double) tv.tv_usec / 1000000.0); (double) gmt.tm_sec + (double) tv.tv_usec / 1000000.0);
#endif #endif
} }
/* /*
* create julian date given time_t value * create julian date given time_t value
*/ */
Julian::Julian(const time_t t) Julian::Julian(const time_t t)
{ {
struct tm ptm; struct tm ptm;
#if WIN32 #if WIN32
if (gmtime_s(&ptm, &t)) if (gmtime_s(&ptm, &t))
{ {
assert(1); assert(1);
} }
#else #else
if (gmtime_r(&t, &ptm) == NULL) if (gmtime_r(&t, &ptm) == NULL)
{ {
assert(1); assert(1);
} }
#endif #endif
int year = ptm.tm_year + 1900; int year = ptm.tm_year + 1900;
double day = ptm.tm_yday + 1 + double day = ptm.tm_yday + 1 +
(ptm.tm_hour + (ptm.tm_hour +
((ptm.tm_min + ((ptm.tm_min +
(ptm.tm_sec / 60.0)) / 60.0)) / 24.0; (ptm.tm_sec / 60.0)) / 60.0)) / 24.0;
Initialize(year, day); Initialize(year, day);
} }
/* /*
* comparison * comparison
*/ */
bool Julian::operator==(const Julian &date) const bool Julian::operator==(const Julian &date) const
{ {
return date_ == date.date_ ? true : false; return date_ == date.date_ ? true : false;
} }
bool Julian::operator!=(const Julian &date) const bool Julian::operator!=(const Julian &date) const
{ {
return !(*this == date); return !(*this == date);
} }
bool Julian::operator>(const Julian &date) const bool Julian::operator>(const Julian &date) const
{ {
return date_ > date.date_ ? true : false; return date_ > date.date_ ? true : false;
} }
bool Julian::operator<(const Julian &date) const { bool Julian::operator<(const Julian &date) const {
return date_ < date.date_ ? true : false; return date_ < date.date_ ? true : false;
} }
bool Julian::operator>=(const Julian &date) const bool Julian::operator>=(const Julian &date) const
{ {
return date_ >= date.date_ ? true : false; return date_ >= date.date_ ? true : false;
} }
bool Julian::operator<=(const Julian &date) const bool Julian::operator<=(const Julian &date) const
{ {
return date_ <= date.date_ ? true : false; return date_ <= date.date_ ? true : false;
} }
/* /*
* assignment * assignment
*/ */
Julian& Julian::operator=(const Julian& b) Julian& Julian::operator=(const Julian& b)
{ {
if (this != &b) { if (this != &b) {
date_ = b.date_; date_ = b.date_;
} }
return (*this); return (*this);
} }
Julian& Julian::operator=(const double b) Julian& Julian::operator=(const double b)
{ {
date_ = b; date_ = b;
return (*this); return (*this);
} }
/* /*
* arithmetic * arithmetic
*/ */
Julian Julian::operator +(const Timespan& b) const Julian Julian::operator +(const Timespan& b) const
{ {
return Julian(*this) += b; return Julian(*this) += b;
} }
Julian Julian::operator-(const Timespan& b) const Julian Julian::operator-(const Timespan& b) const
{ {
return Julian(*this) -= b; return Julian(*this) -= b;
} }
Timespan Julian::operator-(const Julian& b) const Timespan Julian::operator-(const Julian& b) const
{ {
return Timespan(date_ - b.date_); return Timespan(date_ - b.date_);
} }
/* /*
* compound assignment * compound assignment
*/ */
Julian & Julian::operator +=(const Timespan& b) Julian & Julian::operator +=(const Timespan& b)
{ {
date_ += b; date_ += b;
return (*this); return (*this);
} }
Julian & Julian::operator -=(const Timespan& b) Julian & Julian::operator -=(const Timespan& b)
{ {
date_ -= b; date_ -= b;
return (*this); return (*this);
} }
/* /*
* create julian date from year and day of year * create julian date from year and day of year
*/ */
void Julian::Initialize(int year, double day) void Julian::Initialize(int year, double day)
{ {
year--; year--;
int A = (year / 100); int A = (year / 100);
int B = 2 - A + (A / 4); int B = 2 - A + (A / 4);
double new_years = static_cast<int> (365.25 * year) + double new_years = static_cast<int> (365.25 * year) +
static_cast<int> (30.6001 * 14) + static_cast<int> (30.6001 * 14) +
1720994.5 + B; 1720994.5 + B;
date_ = new_years + day; date_ = new_years + day;
} }
/* /*
* create julian date from individual components * create julian date from individual components
* year: 2004 * year: 2004
* mon: 1-12 * mon: 1-12
* day: 1-31 * day: 1-31
* hour: 0-23 * hour: 0-23
* min: 0-59 * min: 0-59
* sec: 0-59.99 * sec: 0-59.99
*/ */
void Julian::Initialize(int year, int mon, int day, void Julian::Initialize(int year, int mon, int day,
int hour, int min, double sec) int hour, int min, double sec)
{ {
// Calculate N, the day of the year (1..366) // Calculate N, the day of the year (1..366)
int N; int N;
int F1 = (int) ((275.0 * mon) / 9.0); int F1 = (int) ((275.0 * mon) / 9.0);
int F2 = (int) ((mon + 9.0) / 12.0); int F2 = (int) ((mon + 9.0) / 12.0);
if (IsLeapYear(year)) if (IsLeapYear(year))
{ {
// Leap year // Leap year
N = F1 - F2 + day - 30; N = F1 - F2 + day - 30;
} }
else else
{ {
// Common year // Common year
N = F1 - (2 * F2) + day - 30; N = F1 - (2 * F2) + day - 30;
} }
double dblDay = N + (hour + (min + (sec / 60.0)) / 60.0) / 24.0; double dblDay = N + (hour + (min + (sec / 60.0)) / 60.0) / 24.0;
Initialize(year, dblDay); Initialize(year, dblDay);
} }
/* /*
* converts time to time_t * converts time to time_t
* note: resolution to seconds only * note: resolution to seconds only
*/ */
time_t Julian::ToTime() const time_t Julian::ToTime() const
{ {
return static_cast<time_t> ((date_ - 2440587.5) * 86400.0); return static_cast<time_t> ((date_ - 2440587.5) * 86400.0);
} }
/* /*
* Greenwich Mean Sidereal Time * Greenwich Mean Sidereal Time
*/ */
double Julian::ToGreenwichSiderealTime() const double Julian::ToGreenwichSiderealTime() const
{ {
#if 0 #if 0
const double UT = fmod(jul + 0.5, 1.0); const double UT = fmod(jul + 0.5, 1.0);
const double TU = (jul - 2451545.0 - UT) / 36525.0; const double TU = (jul - 2451545.0 - UT) / 36525.0;
double GMST = 24110.54841 + TU * double GMST = 24110.54841 + TU *
(8640184.812866 + TU * (0.093104 - TU * 6.2e-06)); (8640184.812866 + TU * (0.093104 - TU * 6.2e-06));
GMST = fmod(GMST + SEC_PER_DAY * OMEGA_E * UT, SEC_PER_DAY); GMST = fmod(GMST + SEC_PER_DAY * OMEGA_E * UT, SEC_PER_DAY);
if (GMST < 0.0) if (GMST < 0.0)
{ {
GMST += SEC_PER_DAY; // "wrap" negative modulo value GMST += SEC_PER_DAY; // "wrap" negative modulo value
} }
return (TWOPI * (GMST / SEC_PER_DAY)); return (TWOPI * (GMST / SEC_PER_DAY));
#endif #endif
// tut1 = Julian centuries from 2000 Jan. 1 12h UT1 // tut1 = Julian centuries from 2000 Jan. 1 12h UT1
// (since J2000 which is 2451545.0) // (since J2000 which is 2451545.0)
// a Julian century is 36525 days // a Julian century is 36525 days
const double tut1 = (date_ - 2451545.0) / 36525.0; const double tut1 = (date_ - 2451545.0) / 36525.0;
// Rotation angle in arcseconds // Rotation angle in arcseconds
double theta = 67310.54841 + (876600.0 * 3600.0 + 8640184.812866) * tut1 double theta = 67310.54841 + (876600.0 * 3600.0 + 8640184.812866) * tut1
+ 0.093104 * pow(tut1, 2.0) - 0.0000062 * pow(tut1, 3.0); + 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 = Util::WrapTwoPI(Util::DegreesToRadians(theta / 240.0)); theta = Util::WrapTwoPI(Util::DegreesToRadians(theta / 240.0));
return theta; return theta;
#if 0 #if 0
static const double C1 = 1.72027916940703639e-2; static const double C1 = 1.72027916940703639e-2;
static const double THGR70 = 1.7321343856509374; static const double THGR70 = 1.7321343856509374;
static const double FK5R = 5.07551419432269442e-15; static const double FK5R = 5.07551419432269442e-15;
/* /*
* get integer number of days from 0 jan 1970 * get integer number of days from 0 jan 1970
*/ */
const double ts70 = date_ - 2433281.5 - 7305.0; const double ts70 = date_ - 2433281.5 - 7305.0;
const double ds70 = floor(ts70 + 1.0e-8); const double ds70 = floor(ts70 + 1.0e-8);
const double tfrac = ts70 - ds70; const double tfrac = ts70 - ds70;
/* /*
* find greenwich location at epoch * find greenwich location at epoch
*/ */
const double c1p2p = C1 + kTWOPI; const double c1p2p = C1 + kTWOPI;
double gsto = Util::WrapTwoPI(THGR70 + C1 * ds70 + c1p2p * tfrac double gsto = Util::WrapTwoPI(THGR70 + C1 * ds70 + c1p2p * tfrac
+ ts70 * ts70 * FK5R); + ts70 * ts70 * FK5R);
return gsto; return gsto;
#endif #endif
} }
/* /*
* 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, kTWOPI); return fmod(ToGreenwichSiderealTime() + lon, kTWOPI);
} }
void Julian::ToGregorian(struct DateTimeComponents* datetime) const void Julian::ToGregorian(struct DateTimeComponents* datetime) const
{ {
double jdAdj = GetDate() + 0.5; double jdAdj = GetDate() + 0.5;
int Z = (int) jdAdj; int Z = (int) jdAdj;
double F = jdAdj - Z; double F = jdAdj - Z;
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;
} }

126
Observer.cpp
View File

@ -1,63 +1,63 @@
#include "Observer.h" #include "Observer.h"
/* /*
* calculate lookangle between the observer and the passed in Eci object * calculate lookangle between the observer and the passed in Eci object
*/ */
CoordTopographic Observer::GetLookAngle(const Eci &eci) CoordTopographic Observer::GetLookAngle(const Eci &eci)
{ {
/* /*
* update the observers Eci to match the time of the Eci passed in * update the observers Eci to match the time of the Eci passed in
* if necessary * if necessary
*/ */
UpdateObserversEci(eci.GetDate()); UpdateObserversEci(eci.GetDate());
/* /*
* calculate differences * calculate differences
*/ */
Vector range_rate = eci.GetVelocity().Subtract(observers_eci_.GetVelocity()); Vector range_rate = eci.GetVelocity().Subtract(observers_eci_.GetVelocity());
Vector range = eci.GetPosition().Subtract(observers_eci_.GetPosition()); Vector range = eci.GetPosition().Subtract(observers_eci_.GetPosition());
range.w = range.GetMagnitude(); range.w = range.GetMagnitude();
/* /*
* Calculate Local Mean Sidereal Time for observers longitude * Calculate Local Mean Sidereal Time for observers longitude
*/ */
double theta = eci.GetDate().ToLocalMeanSiderealTime(geo_.longitude); double theta = eci.GetDate().ToLocalMeanSiderealTime(geo_.longitude);
double sin_lat = sin(geo_.latitude); double sin_lat = sin(geo_.latitude);
double cos_lat = cos(geo_.latitude); double cos_lat = cos(geo_.latitude);
double sin_theta = sin(theta); double sin_theta = sin(theta);
double cos_theta = cos(theta); double cos_theta = cos(theta);
double top_s = sin_lat * cos_theta * range.x double top_s = sin_lat * cos_theta * range.x
+ sin_lat * sin_theta * range.y - cos_lat * range.z; + sin_lat * sin_theta * range.y - cos_lat * range.z;
double top_e = -sin_theta * range.x double top_e = -sin_theta * range.x
+ cos_theta * range.y; + cos_theta * range.y;
double top_z = cos_lat * cos_theta * range.x double top_z = cos_lat * cos_theta * range.x
+ cos_lat * sin_theta * range.y + sin_lat * range.z; + cos_lat * sin_theta * range.y + sin_lat * range.z;
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 += kPI; az += kPI;
} }
if (az < 0.0) if (az < 0.0)
{ {
az += 2.0 * kPI; az += 2.0 * kPI;
} }
double el = asin(top_z / range.w); double el = asin(top_z / range.w);
double rate = range.Dot(range_rate) / range.w; double rate = range.Dot(range_rate) / range.w;
/* /*
* azimuth in radians * azimuth in radians
* elevation in radians * elevation in radians
* range in km * range in km
* range rate in km/s * range rate in km/s
*/ */
return CoordTopographic(az, return CoordTopographic(az,
el, el,
range.w, range.w,
rate); rate);
} }

94
OrbitalElements.cpp
View File

@ -1,47 +1,47 @@
#include "OrbitalElements.h" #include "OrbitalElements.h"
OrbitalElements::OrbitalElements(const Tle& tle) OrbitalElements::OrbitalElements(const Tle& tle)
{ {
/* /*
* extract and format tle data * extract and format tle data
*/ */
mean_anomoly_ = tle.MeanAnomaly(false); mean_anomoly_ = tle.MeanAnomaly(false);
ascending_node_ = tle.RightAscendingNode(false); ascending_node_ = tle.RightAscendingNode(false);
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() * kTWOPI / kMINUTES_PER_DAY; mean_motion_ = tle.MeanMotion() * kTWOPI / kMINUTES_PER_DAY;
bstar_ = tle.BStar(); bstar_ = tle.BStar();
epoch_ = tle.Epoch(); epoch_ = tle.Epoch();
/* /*
* 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(kXKE / MeanMotion(), kTWOTHIRD); const double a1 = pow(kXKE / MeanMotion(), kTWOTHIRD);
const double cosio = cos(Inclination()); const double cosio = cos(Inclination());
const double theta2 = cosio * cosio; const double theta2 = cosio * cosio;
const double x3thm1 = 3.0 * theta2 - 1.0; const double x3thm1 = 3.0 * theta2 - 1.0;
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 * kCK2) * x3thm1 / (betao * betao2); const double temp = (1.5 * kCK2) * 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);
recovered_mean_motion_ = MeanMotion() / (1.0 + del0); recovered_mean_motion_ = MeanMotion() / (1.0 + del0);
/* /*
* alternative way to calculate * alternative way to calculate
* doesnt affect final results * doesnt affect final results
* recovered_semi_major_axis_ = pow(XKE / RecoveredMeanMotion(), TWOTHIRD); * recovered_semi_major_axis_ = pow(XKE / RecoveredMeanMotion(), TWOTHIRD);
*/ */
recovered_semi_major_axis_ = a0 / (1.0 - del0); recovered_semi_major_axis_ = a0 / (1.0 - del0);
/* /*
* find perigee and period * find perigee and period
*/ */
perigee_ = (RecoveredSemiMajorAxis() * (1.0 - Eccentricity()) - kAE) * kXKMPER; perigee_ = (RecoveredSemiMajorAxis() * (1.0 - Eccentricity()) - kAE) * kXKMPER;
period_ = kTWOPI / RecoveredMeanMotion(); period_ = kTWOPI / RecoveredMeanMotion();
} }

476
RunTest.cpp
View File

@ -1,238 +1,238 @@
#include "Julian.h" #include "Julian.h"
#include "Tle.h" #include "Tle.h"
#include "SGP4.h" #include "SGP4.h"
#include "Globals.h" #include "Globals.h"
#include "Util.h" #include "Util.h"
#include "Observer.h" #include "Observer.h"
#include "CoordGeodetic.h" #include "CoordGeodetic.h"
#include "CoordTopographic.h" #include "CoordTopographic.h"
#include <list> #include <list>
#include <string> #include <string>
#include <iomanip> #include <iomanip>
#include <iostream> #include <iostream>
#include <fstream> #include <fstream>
#include <vector> #include <vector>
#include <cstdlib> #include <cstdlib>
void RunTle(Tle tle, double start, double end, double inc) void RunTle(Tle tle, double start, double end, double inc)
{ {
double current = start; double current = start;
SGP4 model(tle); SGP4 model(tle);
bool running = true; bool running = true;
bool first_run = true; bool first_run = true;
std::cout << " " << std::setprecision(0) << tle.NoradNumber() << " xx" << std::endl; std::cout << " " << std::setprecision(0) << tle.NoradNumber() << " xx" << std::endl;
while (running) while (running)
{ {
try try
{ {
double val; double val;
if (first_run && current != 0.0) if (first_run && current != 0.0)
{ {
/* /*
* make sure first run is always as zero * make sure first run is always as zero
*/ */
val = 0.0; val = 0.0;
} }
else else
{ {
/* /*
* otherwise run as normal * otherwise run as normal
*/ */
val = current; val = current;
} }
Eci eci = model.FindPosition(val); Eci eci = model.FindPosition(val);
Vector position = eci.GetPosition(); Vector position = eci.GetPosition();
Vector velocity = eci.GetVelocity(); Vector velocity = eci.GetVelocity();
std::cout << std::setprecision(8) << std::fixed; std::cout << std::setprecision(8) << std::fixed;
std::cout.width(17); std::cout.width(17);
std::cout << val << " "; std::cout << val << " ";
std::cout.width(16); std::cout.width(16);
std::cout << position.x << " "; std::cout << position.x << " ";
std::cout.width(16); std::cout.width(16);
std::cout << position.y << " "; std::cout << position.y << " ";
std::cout.width(16); std::cout.width(16);
std::cout << position.z << " "; std::cout << position.z << " ";
std::cout << std::setprecision(9) << std::fixed; std::cout << std::setprecision(9) << std::fixed;
std::cout.width(14); std::cout.width(14);
std::cout << velocity.x << " "; std::cout << velocity.x << " ";
std::cout.width(14); std::cout.width(14);
std::cout << velocity.y << " "; std::cout << velocity.y << " ";
std::cout.width(14); std::cout.width(14);
std::cout << velocity.z << std::endl; std::cout << velocity.z << std::endl;
} }
catch (SatelliteException& e) catch (SatelliteException& e)
{ {
std::cout << e.what() << std::endl; std::cout << e.what() << std::endl;
running = false; running = false;
} }
if ((first_run && current == 0.0) || !first_run) if ((first_run && current == 0.0) || !first_run)
{ {
if (current == end) if (current == end)
{ {
running = false; running = false;
} }
else if (current + inc > end) else if (current + inc > end)
{ {
current = end; current = end;
} }
else else
{ {
current += inc; current += inc;
} }
} }
first_run = false; first_run = false;
} }
} }
void tokenize(const std::string& str, std::vector<std::string>& tokens) void tokenize(const std::string& str, std::vector<std::string>& tokens)
{ {
const std::string& delimiters = " "; const std::string& delimiters = " ";
/* /*
* skip delimiters at beginning * skip delimiters at beginning
*/ */
std::string::size_type last_pos = str.find_first_not_of(delimiters, 0); std::string::size_type last_pos = str.find_first_not_of(delimiters, 0);
/* /*
* find first non-delimiter * find first non-delimiter
*/ */
std::string::size_type pos = str.find_first_of(delimiters, last_pos); std::string::size_type pos = str.find_first_of(delimiters, last_pos);
while (std::string::npos != pos || std::string::npos != last_pos) while (std::string::npos != pos || std::string::npos != last_pos)
{ {
/* /*
* add found token to vector * add found token to vector
*/ */
tokens.push_back(str.substr(last_pos, pos - last_pos)); tokens.push_back(str.substr(last_pos, pos - last_pos));
/* /*
* skip delimiters * skip delimiters
*/ */
last_pos = str.find_first_not_of(delimiters, pos); last_pos = str.find_first_not_of(delimiters, pos);
/* /*
* find next non-delimiter * find next non-delimiter
*/ */
pos = str.find_first_of(delimiters, last_pos); pos = str.find_first_of(delimiters, last_pos);
} }
} }
void RunTest(const char* infile) void RunTest(const char* infile)
{ {
std::ifstream file; std::ifstream file;
file.open(infile); file.open(infile);
if (!file.is_open()) if (!file.is_open())
{ {
std::cerr << "Error opening file" << std::endl; std::cerr << "Error opening file" << std::endl;
return; return;
} }
bool got_first_line = false; bool got_first_line = false;
std::string line1; std::string line1;
std::string line2; std::string line2;
std::string parameters; std::string parameters;
while (!file.eof()) while (!file.eof())
{ {
std::string line; std::string line;
std::getline(file, line); std::getline(file, line);
Util::Trim(line); Util::Trim(line);
/* /*
* skip blank lines or lines starting with # * skip blank lines or lines starting with #
*/ */
if (line.length() == 0 || line[0] == '#') if (line.length() == 0 || line[0] == '#')
{ {
got_first_line = false; got_first_line = false;
continue; continue;
} }
/* /*
* find first line * find first line
*/ */
if (!got_first_line) if (!got_first_line)
{ {
try try
{ {
Tle::IsValidLine(line, 1); Tle::IsValidLine(line, 1);
/* /*
* store line and now read in second line * store line and now read in second line
*/ */
got_first_line = true; got_first_line = true;
line1 = line; line1 = line;
} }
catch (TleException& e) catch (TleException& e)
{ {
std::cerr << "Error: " << e.what() << std::endl; std::cerr << "Error: " << e.what() << std::endl;
std::cerr << line << std::endl; std::cerr << line << std::endl;
} }
} }
else else
{ {
/* /*
* no second chances, second line should follow the first * no second chances, second line should follow the first
*/ */
got_first_line = false; got_first_line = false;
/* /*
* split line, first 69 is the second line of the tle * split line, first 69 is the second line of the tle
* the rest is the test parameters, if there is any * the rest is the test parameters, if there is any
*/ */
line2 = line.substr(0, 69); line2 = line.substr(0, 69);
double start = 0.0; double start = 0.0;
double end = 1440.0; double end = 1440.0;
double inc = 120.0; double inc = 120.0;
if (line.length() > 69) if (line.length() > 69)
{ {
std::vector<std::string> tokens; std::vector<std::string> tokens;
parameters = line.substr(70, line.length() - 69); parameters = line.substr(70, line.length() - 69);
tokenize(parameters, tokens); tokenize(parameters, tokens);
if (tokens.size() >= 3) if (tokens.size() >= 3)
{ {
start = atof(tokens[0].c_str()); start = atof(tokens[0].c_str());
end = atof(tokens[1].c_str()); end = atof(tokens[1].c_str());
inc = atof(tokens[2].c_str()); inc = atof(tokens[2].c_str());
} }
} }
/* /*
* following line must be the second line * following line must be the second line
*/ */
try try
{ {
Tle::IsValidLine(line2, 2); Tle::IsValidLine(line2, 2);
Tle tle("Test", line1, line2); Tle tle("Test", line1, line2);
RunTle(tle, start, end, inc); RunTle(tle, start, end, inc);
} }
catch (TleException& e) catch (TleException& e)
{ {
std::cerr << "Error: " << e.what() << std::endl; std::cerr << "Error: " << e.what() << std::endl;
std::cerr << line << std::endl; std::cerr << line << std::endl;
} }
} }
} }
/* /*
* close file * close file
*/ */
file.close(); file.close();
return; return;
} }
int main() int main()
{ {
const char* file_name = "SGP4-VER.TLE"; const char* file_name = "SGP4-VER.TLE";
RunTest(file_name); RunTest(file_name);
return 1; return 1;
} }

2494
SGP4.cpp
View File

File diff suppressed because it is too large Load Diff

92
SolarPosition.cpp
View File

@ -1,46 +1,46 @@
#include "SolarPosition.h" #include "SolarPosition.h"
#include "Globals.h" #include "Globals.h"
#include "Util.h" #include "Util.h"
#include <cmath> #include <cmath>
Eci SolarPosition::FindPosition(const Julian& j) Eci SolarPosition::FindPosition(const Julian& j)
{ {
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) / kSECONDS_PER_DAY) / 36525.0; const double T = (mjd + Delta_ET(year) / kSECONDS_PER_DAY) / 36525.0;
const double M = Util::DegreesToRadians(Util::Wrap360(358.47583 const double M = Util::DegreesToRadians(Util::Wrap360(358.47583
+ Util::Wrap360(35999.04975 * T) + Util::Wrap360(35999.04975 * T)
- (0.000150 + 0.0000033 * T) * T * T)); - (0.000150 + 0.0000033 * T) * T * T));
const double L = Util::DegreesToRadians(Util::Wrap360(279.69668 const double L = Util::DegreesToRadians(Util::Wrap360(279.69668
+ Util::Wrap360(36000.76892 * T) + Util::Wrap360(36000.76892 * T)
+ 0.0003025 * T*T)); + 0.0003025 * T*T));
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 = Util::DegreesToRadians((1.919460 const double C = Util::DegreesToRadians((1.919460
- (0.004789 + 0.000014 * T) * T) * sin(M) - (0.004789 + 0.000014 * T) * T) * sin(M)
+ (0.020094 - 0.000100 * T) * sin(2 * M) + (0.020094 - 0.000100 * T) * sin(2 * M)
+ 0.000293 * sin(3 * M)); + 0.000293 * sin(3 * M));
const double O = Util::DegreesToRadians( const double O = Util::DegreesToRadians(
Util::Wrap360(259.18 - 1934.142 * T)); Util::Wrap360(259.18 - 1934.142 * T));
const double Lsa = Util::WrapTwoPI(L + C const double Lsa = Util::WrapTwoPI(L + C
- Util::DegreesToRadians(0.00569 - 0.00479 * sin(O))); - Util::DegreesToRadians(0.00569 - 0.00479 * sin(O)));
const double nu = Util::WrapTwoPI(M + C); const double nu = Util::WrapTwoPI(M + C);
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 = Util::DegreesToRadians(23.452294 - (0.0130125 const double eps = Util::DegreesToRadians(23.452294 - (0.0130125
+ (0.00000164 - 0.000000503 * T) * T) * T + 0.00256 * cos(O)); + (0.00000164 - 0.000000503 * T) * T) * T + 0.00256 * cos(O));
R = R * kAU; 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),
R * sin(Lsa) * sin(eps), R * sin(Lsa) * sin(eps),
R); R);
return Eci(j, solar_position); return Eci(j, solar_position);
} }
double SolarPosition::Delta_ET(double year) const double SolarPosition::Delta_ET(double year) const
{ {
return 26.465 + 0.747622 * (year - 1950) + 1.886913 return 26.465 + 0.747622 * (year - 1950) + 1.886913
* sin(kTWOPI * (year - 1975) / 33); * sin(kTWOPI * (year - 1975) / 33);
} }

348
Timespan.cpp
View File

@ -1,175 +1,175 @@
#include "Timespan.h" #include "Timespan.h"
#include "Globals.h" #include "Globals.h"
Timespan::Timespan() Timespan::Timespan()
: time_span_(0.0) { : time_span_(0.0) {
} }
Timespan::Timespan(const unsigned int days, const unsigned int hours, Timespan::Timespan(const unsigned int days, const unsigned int hours,
const unsigned int minutes, const double seconds) { const unsigned int minutes, const double seconds) {
SetValue(days, hours, minutes, seconds); SetValue(days, hours, minutes, seconds);
} }
Timespan::Timespan(const double b) { Timespan::Timespan(const double b) {
time_span_ = b; time_span_ = b;
} }
Timespan::Timespan(const Timespan& b) { Timespan::Timespan(const Timespan& b) {
time_span_ = b.time_span_; time_span_ = b.time_span_;
} }
Timespan::~Timespan(void) { Timespan::~Timespan(void) {
} }
void Timespan::SetValue(const unsigned int days, const unsigned int hours, void Timespan::SetValue(const unsigned int days, const unsigned int hours,
const unsigned int minutes, const double seconds) { const unsigned int minutes, const double seconds) {
time_span_ = static_cast<double> (days); time_span_ = static_cast<double> (days);
AddHours(hours); AddHours(hours);
AddMinutes(minutes); AddMinutes(minutes);
AddSeconds(seconds); AddSeconds(seconds);
} }
void Timespan::AddDays(const unsigned int days) { void Timespan::AddDays(const unsigned int days) {
time_span_ += static_cast<double> (days); time_span_ += static_cast<double> (days);
} }
void Timespan::AddHours(const unsigned int hours) { void Timespan::AddHours(const unsigned int hours) {
time_span_ += (static_cast<double> (hours) / kHOURS_PER_DAY); time_span_ += (static_cast<double> (hours) / kHOURS_PER_DAY);
} }
void Timespan::AddMinutes(const unsigned int minutes) { void Timespan::AddMinutes(const unsigned int minutes) {
time_span_ += (static_cast<double> (minutes) / kMINUTES_PER_DAY); time_span_ += (static_cast<double> (minutes) / kMINUTES_PER_DAY);
} }
void Timespan::AddSeconds(const double seconds) { void Timespan::AddSeconds(const double seconds) {
time_span_ += (seconds / kSECONDS_PER_DAY); time_span_ += (seconds / kSECONDS_PER_DAY);
} }
double Timespan::GetTotalDays() const { double Timespan::GetTotalDays() const {
return time_span_; return time_span_;
} }
double Timespan::GetTotalHours() const { double Timespan::GetTotalHours() const {
return time_span_ * kHOURS_PER_DAY; return time_span_ * kHOURS_PER_DAY;
} }
double Timespan::GetTotalMinutes() const { double Timespan::GetTotalMinutes() const {
return time_span_ * kMINUTES_PER_DAY; return time_span_ * kMINUTES_PER_DAY;
} }
double Timespan::GetTotalSeconds() const { double Timespan::GetTotalSeconds() const {
return time_span_ * kSECONDS_PER_DAY; return time_span_ * kSECONDS_PER_DAY;
} }
Timespan& Timespan::operator =(const Timespan& b) { Timespan& Timespan::operator =(const Timespan& b) {
if (this != &b) { if (this != &b) {
time_span_ = b.time_span_; time_span_ = b.time_span_;
} }
return (*this); return (*this);
} }
Timespan Timespan::operator +(const Timespan& b) const { Timespan Timespan::operator +(const Timespan& b) const {
return Timespan(*this) += b; return Timespan(*this) += b;
} }
Timespan Timespan::operator -(const Timespan& b) const { Timespan Timespan::operator -(const Timespan& b) const {
return Timespan(*this) -= b; return Timespan(*this) -= b;
} }
Timespan Timespan::operator/(const double b) const { Timespan Timespan::operator/(const double b) const {
return Timespan(*this) /= b; return Timespan(*this) /= b;
} }
Timespan Timespan::operator*(const double b) const { Timespan Timespan::operator*(const double b) const {
return Timespan(*this) *= b; return Timespan(*this) *= b;
} }
Timespan & Timespan::operator+=(const Timespan& b) { Timespan & Timespan::operator+=(const Timespan& b) {
time_span_ += b.time_span_; time_span_ += b.time_span_;
return (*this); return (*this);
} }
Timespan & Timespan::operator-=(const Timespan& b) { Timespan & Timespan::operator-=(const Timespan& b) {
time_span_ -= b.time_span_; time_span_ -= b.time_span_;
return (*this); return (*this);
} }
Timespan & Timespan::operator/=(const double b) { Timespan & Timespan::operator/=(const double b) {
time_span_ /= b; time_span_ /= b;
return (*this); return (*this);
} }
Timespan & Timespan::operator*=(const double b) { Timespan & Timespan::operator*=(const double b) {
time_span_ *= b; time_span_ *= b;
return (*this); return (*this);
} }
bool Timespan::operator ==(const Timespan& b) const { bool Timespan::operator ==(const Timespan& b) const {
if (time_span_ == b.time_span_) if (time_span_ == b.time_span_)
return true; return true;
else else
return false; return false;
} }
bool Timespan::operator !=(const Timespan& b) const { bool Timespan::operator !=(const Timespan& b) const {
return !(*this == b); return !(*this == b);
} }
bool Timespan::operator>(const Timespan& b) const { bool Timespan::operator>(const Timespan& b) const {
if (time_span_ > b.time_span_) if (time_span_ > b.time_span_)
return true; return true;
else else
return false; return false;
} }
bool Timespan::operator<(const Timespan& b) const { bool Timespan::operator<(const Timespan& b) const {
if (time_span_ < b.time_span_) if (time_span_ < b.time_span_)
return true; return true;
else else
return false; return false;
} }
bool Timespan::operator >=(const Timespan& b) const { bool Timespan::operator >=(const Timespan& b) const {
if (time_span_ >= b.time_span_) if (time_span_ >= b.time_span_)
return true; return true;
else else
return false; return false;
} }
bool Timespan::operator <=(const Timespan & b) const { bool Timespan::operator <=(const Timespan & b) const {
if (time_span_ <= b.time_span_) if (time_span_ <= b.time_span_)
return true; return true;
else else
return false; return false;
} }
double& operator +=(double& a, const Timespan& b) { double& operator +=(double& a, const Timespan& b) {
a += b.time_span_; a += b.time_span_;
return a; return a;
} }
double& operator -=(double& a, const Timespan& b) { double& operator -=(double& a, const Timespan& b) {
a -= b.time_span_; a -= b.time_span_;
return a; return a;
} }