Misc formatting

Program now loops around and prints az/el for a set of satellites

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);

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
  */

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() {
     };
 

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;
+}

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
      */

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;
+    Julian julian = Epoch();
-    std::cout.width(17);
+    julian.AddMin(tsince);
-    std::cout << tsince << " ";
+    eci = Eci(julian, position, velocity);
-    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;
 }
 
 /*

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 FindPositionSDP4(Eci& eci, double tsince);
-    void FindPositionSGP4(double tsince);
+    void FindPositionSGP4(Eci& eci, double tsince);
-    void CalculateFinalPositionVelocity(const double& tsince, const double& e,
+    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,

SatelliteException.h

@@ -4,20 +4,22 @@
 #include <exception>
 #include <iostream>
 
-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()
+    SatelliteException(const char* message)
-	{
+    : message_(message) {
-		return message_.c_str();
+    }
-	}
+
+    virtual ~SatelliteException(void) throw () {
+    }
+
+    virtual const char* what() const throw () {
+        return message_.c_str();
+    }
 
 private:
-	std::string message_;
+    std::string message_;
 };
 
 #endif

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);
 }
 

main.cpp

@@ -2,6 +2,8 @@
 #include "Tle.h"
 #include "SGDP4.h"
 #include "Globals.h"
+#include "Observer.h"
+#include "Coord.h"
 
 #include <list>
 #include <string>
@@ -12,7 +14,6 @@ void RunTest();
 int main() {
 
     std::list<Tle> 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",
+    //tles.push_back(Tle("SDP4 Test",
-            "1 11801U          80230.29629788  .01431103  00000-0  14311-1      13",
+    //        "1 11801U          80230.29629788  .01431103  00000-0  14311-1      13",
-            "2 11801  46.7916 230.4354 7318036  47.4722  10.4117  2.28537848    13"));
+    //        "2 11801  46.7916 230.4354 7318036  47.4722  10.4117  2.28537848    13"));
-    std::list<Tle>::iterator itr;
+
-    Julian jul;
+    Observer obs(51.360242, 0.101473, 0.07);
+
+    std::list<Tle>::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;