From d497544ca21b032a8f05ab4fe8cd6c43a8ef31b2 Mon Sep 17 00:00:00 2001 From: Daniel Warner Date: Wed, 30 Mar 2011 18:16:37 +0100 Subject: [PATCH] Misc formatting Program now loops around and prints az/el for a set of satellites --- Eci.h | 3 ++ Julian.cpp | 8 +++++ Julian.h | 3 ++ Observer.cpp | 69 ++++++++++++++++++++++++++++++++++++++++++-- Observer.h | 16 ++++++++++ SGDP4.cpp | 31 +++++++++----------- SGDP4.h | 9 +++--- SatelliteException.h | 22 +++++++------- Tle.cpp | 4 +-- main.cpp | 57 +++++++++++++++++++++++++++--------- 10 files changed, 172 insertions(+), 50 deletions(-) diff --git a/Eci.h b/Eci.h index a14940a..6b5cd85 100644 --- a/Eci.h +++ b/Eci.h @@ -7,6 +7,9 @@ class Eci { public: + + Eci() { + }; Eci(const Julian &date, const CoordGeodetic &geo); Eci(const Julian &date, const Vector &position, const Vector &velocity); virtual ~Eci(void); diff --git a/Julian.cpp b/Julian.cpp index 5692a2b..3985fcf 100644 --- a/Julian.cpp +++ b/Julian.cpp @@ -80,6 +80,14 @@ Julian::Julian(int year, int mon, int day, int hour, int min, double sec) { Initialize(year, mon, day, hour, min, sec); } +bool Julian::operator==(const Julian &date) const { + return date_ == date.date_ ? true : false; +} + +bool Julian::operator!=(const Julian &date) const { + return date_ == date.date_ ? false : true; +} + /* * create julian date from year and day of year */ diff --git a/Julian.h b/Julian.h index b101f9d..5a4bf96 100644 --- a/Julian.h +++ b/Julian.h @@ -13,6 +13,9 @@ public: Julian(int year, double day); Julian(int year, int mon, int day, int hour, int min, double sec); + bool operator==(const Julian &date) const; + bool operator!=(const Julian &date) const; + ~Julian() { }; diff --git a/Observer.cpp b/Observer.cpp index b6920f1..9a4d3a1 100644 --- a/Observer.cpp +++ b/Observer.cpp @@ -5,16 +5,81 @@ /* * in degrees! */ -Observer::Observer(const double latitude, const double longitude, const double altitude) { +Observer::Observer(const double latitude, const double longitude, const double altitude) +: observers_eci_(Julian(), CoordGeodetic()) { geo_.SetLatitude(Globals::Deg2Rad(latitude)); geo_.SetLongitude(Globals::Deg2Rad(longitude)); geo_.SetAltitude(altitude); } Observer::Observer(const CoordGeodetic &geo) -: geo_(geo) { +: geo_(geo), observers_eci_(Julian(), geo) { } Observer::~Observer(void) { } + +void Observer::UpdateObserversEci(const Julian &date) { + + if (observers_eci_.GetDate() != date) { + observers_eci_ = Eci(date, geo_); + } +} + +CoordTopographic Observer::GetLookAngle(const Eci &eci) { + + Julian date = eci.GetDate(); + + /* + * update observers eci value if the date is different to the eci passed in + */ + UpdateObserversEci(date); + + Vector range_rate(eci.GetVelocity().GetX() - observers_eci_.GetVelocity().GetX(), + eci.GetVelocity().GetY() - observers_eci_.GetVelocity().GetY(), + eci.GetVelocity().GetZ() - observers_eci_.GetVelocity().GetZ()); + + const double x = eci.GetPosition().GetX() - observers_eci_.GetPosition().GetX(); + const double y = eci.GetPosition().GetY() - observers_eci_.GetPosition().GetY(); + const double z = eci.GetPosition().GetZ() - observers_eci_.GetPosition().GetZ(); + const double w = sqrt(pow(x, 2.0) + pow(y, 2.0) + pow(z, 2.0)); + + Vector range(x, y, z, w); + + // The site's Local Mean Sidereal Time at the time of interest. + double theta = date.ToLMST(geo_.GetLongitude()); + + double sin_lat = sin(geo_.GetLatitude()); + double cos_lat = cos(geo_.GetLatitude()); + double sin_theta = sin(theta); + double cos_theta = cos(theta); + + double top_s = sin_lat * cos_theta * range.GetX() + + sin_lat * sin_theta * range.GetY() - + cos_lat * range.GetZ(); + double top_e = -sin_theta * range.GetX() + + cos_theta * range.GetY(); + double top_z = cos_lat * cos_theta * range.GetX() + + cos_lat * sin_theta * range.GetY() + + sin_lat * range.GetZ(); + double az = atan(-top_e / top_s); + + if (top_s > 0.0) + az += Globals::PI(); + + if (az < 0.0) + az += 2.0 * Globals::PI(); + + double el = asin(top_z / range.GetW()); + double rate = (range.GetX() * range_rate.GetX() + + range.GetY() * range_rate.GetY() + + range.GetZ() * range_rate.GetZ()) / range.GetW(); + + CoordTopographic topo(az, // azimuth, radians + el, // elevation, radians + range.GetW(), // range, km + rate); // rate, km / sec + + return topo; +} diff --git a/Observer.h b/Observer.h index fee247a..e6acc27 100644 --- a/Observer.h +++ b/Observer.h @@ -2,6 +2,8 @@ #define OBSERVER_H_ #include "Coord.h" +#include "Julian.h" +#include "Eci.h" class Observer { public: @@ -11,13 +13,27 @@ public: void SetLocation(const CoordGeodetic& geo) { geo_ = geo; + observers_eci_ = Eci(observers_eci_.GetDate(), geo_); } CoordGeodetic GetLocation() const { return geo_; } + Eci GetEciPosition(const Julian &date) const { + return Eci(date, geo_); + } + + CoordTopographic GetLookAngle(const Eci &eci); + private: + void UpdateObserversEci(const Julian &date); + + /* + * the observers eci for a particular time + */ + Eci observers_eci_; + /* * the observers position */ diff --git a/SGDP4.cpp b/SGDP4.cpp index a62b3cb..337206c 100644 --- a/SGDP4.cpp +++ b/SGDP4.cpp @@ -237,20 +237,21 @@ void SGDP4::Initialize(const double& theta2, const double& betao2, const double& } } - FindPosition(0.0); + Eci eci; + FindPosition(eci, 0.0); first_run_ = false; } -void SGDP4::FindPosition(double tsince) { +void SGDP4::FindPosition(Eci& eci, double tsince) { if (i_use_deep_space_) - FindPositionSDP4(tsince); + FindPositionSDP4(eci, tsince); else - FindPositionSGP4(tsince); + FindPositionSGP4(eci, tsince); } -void SGDP4::FindPositionSDP4(double tsince) { +void SGDP4::FindPositionSDP4(Eci& eci, double tsince) { /* * the final values @@ -320,7 +321,7 @@ void SGDP4::FindPositionSDP4(double tsince) { /* * using calculated values, find position and velocity */ - CalculateFinalPositionVelocity(tsince, e, + CalculateFinalPositionVelocity(eci, tsince, e, a, omega, xl, xnode, xincl, perturbed_xlcof, perturbed_aycof, perturbed_x3thm1, perturbed_x1mth2, perturbed_x7thm1, @@ -328,7 +329,7 @@ void SGDP4::FindPositionSDP4(double tsince) { } -void SGDP4::FindPositionSGP4(double tsince) { +void SGDP4::FindPositionSGP4(Eci& eci, double tsince) { /* * the final values @@ -381,7 +382,7 @@ void SGDP4::FindPositionSGP4(double tsince) { * using calculated values, find position and velocity * we can pass in constants from Initialize() as these dont change */ - CalculateFinalPositionVelocity(tsince, e, + CalculateFinalPositionVelocity(eci, tsince, e, a, omega, xl, xnode, xincl, i_xlcof_, i_aycof_, i_x3thm1_, i_x1mth2_, i_x7thm1_, @@ -389,7 +390,7 @@ void SGDP4::FindPositionSGP4(double tsince) { } -void SGDP4::CalculateFinalPositionVelocity(const double& tsince, const double& e, +void SGDP4::CalculateFinalPositionVelocity(Eci& eci, const double& tsince, const double& e, const double& a, const double& omega, const double& xl, const double& xnode, const double& xincl, const double& xlcof, const double& aycof, const double& x3thm1, const double& x1mth2, const double& x7thm1, @@ -548,15 +549,9 @@ void SGDP4::CalculateFinalPositionVelocity(const double& tsince, const double& e const double zdot = (rdotk * uz + rfdotk * vz) * constants_.XKMPER / 60.0; Vector velocity(xdot, ydot, zdot); - std::cout << std::setprecision(8) << std::fixed; - std::cout.width(17); - std::cout << tsince << " "; - std::cout.width(17); - std::cout << position.GetX() << " "; - std::cout.width(17); - std::cout << position.GetY() << " "; - std::cout.width(17); - std::cout << position.GetZ() << std::endl; + Julian julian = Epoch(); + julian.AddMin(tsince); + eci = Eci(julian, position, velocity); } /* diff --git a/SGDP4.h b/SGDP4.h index 9da3696..32b2d56 100644 --- a/SGDP4.h +++ b/SGDP4.h @@ -2,6 +2,7 @@ #define SGDP4_H_ #include "Tle.h" +#include "Eci.h" class SGDP4 { public: @@ -16,7 +17,7 @@ public: void SetConstants(EnumConstants constants); void SetTle(const Tle& tle); - void FindPosition(double tsince); + void FindPosition(Eci& eci, double tsince); private: void Initialize(const double& theta2, const double& betao2, const double& betao, const double& eosq); @@ -27,9 +28,9 @@ private: double& omgasm, double& xnodes, double& xll); void DeepSecular(const double& t, double& xll, double& omgasm, double& xnodes, double& em, double& xinc, double& xn); - void FindPositionSDP4(double tsince); - void FindPositionSGP4(double tsince); - void CalculateFinalPositionVelocity(const double& tsince, const double& e, + void FindPositionSDP4(Eci& eci, double tsince); + void FindPositionSGP4(Eci& eci, double tsince); + void CalculateFinalPositionVelocity(Eci& eci, const double& tsince, const double& e, const double& a, const double& omega, const double& xl, const double& xnode, const double& xincl, const double& xlcof, const double& aycof, const double& x3thm1, const double& x1mth2, const double& x7thm1, diff --git a/SatelliteException.h b/SatelliteException.h index 487d506..fc35b08 100644 --- a/SatelliteException.h +++ b/SatelliteException.h @@ -4,20 +4,22 @@ #include #include -class SatelliteException : public std::exception -{ +class SatelliteException : public std::exception { public: - SatelliteException(const char* message) - : message_(message) {} - virtual ~SatelliteException(void) throw() {} - virtual const char* what() const throw() - { - return message_.c_str(); - } + SatelliteException(const char* message) + : message_(message) { + } + + virtual ~SatelliteException(void) throw () { + } + + virtual const char* what() const throw () { + return message_.c_str(); + } private: - std::string message_; + std::string message_; }; #endif \ No newline at end of file diff --git a/Tle.cpp b/Tle.cpp index 686f339..c235628 100644 --- a/Tle.cpp +++ b/Tle.cpp @@ -125,7 +125,7 @@ void Tle::Initialize() { TrimRight(line_one_); TrimLeft(line_two_); TrimRight(line_two_); - + assert(!name_.empty()); assert(!line_one_.empty()); assert(!line_two_.empty()); @@ -198,7 +198,7 @@ void Tle::Initialize() { else year += 1900; double day = GetField(Tle::FLD_EPOCHDAY); - + date_ = Julian(year, day); } diff --git a/main.cpp b/main.cpp index e95e450..f45ce42 100644 --- a/main.cpp +++ b/main.cpp @@ -2,6 +2,8 @@ #include "Tle.h" #include "SGDP4.h" #include "Globals.h" +#include "Observer.h" +#include "Coord.h" #include #include @@ -12,7 +14,6 @@ void RunTest(); int main() { std::list tles; -#if 0 tles.push_back(Tle("ALSAT 1 ", "1 27559U 02054A 11089.12872679 .00000243 00000-0 48843-4 0 4131", @@ -47,25 +48,52 @@ int main() { tles.push_back(Tle("UK-DMC 2 ", "1 35683U 09041C 11089.11558659 .00000272 00000-0 54146-4 0 8712", "2 35683 98.0762 348.1067 0001434 99.8921 260.2456 14.69414094 89293")); -#endif + //tles.push_back(Tle("SGP4 Test", // "1 88888U 80275.98708465 .00073094 13844-3 66816-4 0 8", // "2 88888 72.8435 115.9689 0086731 52.6988 110.5714 16.05824518 105")); - tles.push_back(Tle("SDP4 Test", - "1 11801U 80230.29629788 .01431103 00000-0 14311-1 13", - "2 11801 46.7916 230.4354 7318036 47.4722 10.4117 2.28537848 13")); - std::list::iterator itr; - Julian jul; + //tles.push_back(Tle("SDP4 Test", + // "1 11801U 80230.29629788 .01431103 00000-0 14311-1 13", + // "2 11801 46.7916 230.4354 7318036 47.4722 10.4117 2.28537848 13")); + + Observer obs(51.360242, 0.101473, 0.07); + + std::list::iterator itr; + + while (true) { + Julian date_now; + + for (itr = tles.begin(); itr != tles.end(); itr++) { + SGDP4 model; + double tsince = date_now.SpanMin((*itr).GetEpoch()); + model.SetTle(*itr); + Eci eci; + model.FindPosition(eci, tsince); + CoordTopographic topo = obs.GetLookAngle(eci); + + std::cout << std::setprecision(8) << std::fixed; + std::cout.width(17); + std::cout << tsince << " "; + std::cout.width(17); + std::cout << Globals::Rad2Deg(topo.GetAzimuth()) << " "; + std::cout.width(17); + std::cout << Globals::Rad2Deg(topo.GetElevation()) << std::endl; + - for (itr = tles.begin(); itr != tles.end(); itr++) { - SGDP4 model; - model.SetTle(*itr); - for (int i = 0; i < 5; i++) { - model.FindPosition(i * 360.0); } } - + /* + std::cout << std::setprecision(8) << std::fixed; + std::cout.width(17); + std::cout << tsince << " "; + std::cout.width(17); + std::cout << position.GetX() << " "; + std::cout.width(17); + std::cout << position.GetY() << " "; + std::cout.width(17); + std::cout << position.GetZ() << std::endl; + */ //RunTest(); return 0; } @@ -78,7 +106,8 @@ void RunTle(Tle tle, double start, double end, double inc) { std::cout << "Satellite: " << tle.GetName() << std::endl << std::endl; while (running) { try { - model.FindPosition(current); + Eci eci; + model.FindPosition(eci, current); } catch (std::exception* ex) { std::cout << ex->what() << std::endl; break;