From bbed197a1981d7db628031fe10c5499d29d2d0d9 Mon Sep 17 00:00:00 2001
From: Daniel Warner <contact@danrw.com>
Date: Sat, 9 Apr 2011 21:40:45 +0100
Subject: [PATCH] Added ToGeodetic, to convert satellite position to
 latitude/longitude/altitude

---
 Eci.cpp | 28 ++++++++++++++++++++++++++++
 Eci.h   |  2 ++
 2 files changed, 30 insertions(+)

diff --git a/Eci.cpp b/Eci.cpp
index 6ddffb2..f95d773 100644
--- a/Eci.cpp
+++ b/Eci.cpp
@@ -50,3 +50,31 @@ Eci::Eci(const Julian &date, const Vector &position, const Vector &velocity)
 
 Eci::~Eci(void) {
 }
+
+CoordGeodetic Eci::ToGeodetic() const {
+
+    const double theta = Globals::AcTan(position_.GetY(), position_.GetX());
+    const double lon = Globals::Fmod2p(theta - date_.ToGreenwichSiderealTime());
+    const double r = sqrt((position_.GetX() * position_.GetX()) + (position_.GetY() * position_.GetY()));
+    const double e2 = Globals::F() * (2.0 - Globals::F());
+
+    double lat = Globals::AcTan(position_.GetZ(), r);
+    double phi = 0.0;
+    double c = 0.0;
+
+    do {
+        phi = lat;
+        const double sinphi = sin(phi);
+        c = 1 / sqrt(1 - e2 * sinphi * sinphi);
+        lat = Globals::AcTan(position_.GetZ() + Globals::XKMPER() * c * e2 * sinphi, r);
+    } while (fabs(lat - phi) >= 1e-10);
+
+    const double alt = r / cos(lat) - Globals::XKMPER() * c;
+
+    /*
+     * todo: check
+     */
+    //if (lat > pio2) geodetic->lat -= twopi;
+
+    return CoordGeodetic(lat, lon, alt);
+}
diff --git a/Eci.h b/Eci.h
index 12b8215..3af2ea4 100644
--- a/Eci.h
+++ b/Eci.h
@@ -26,6 +26,8 @@ public:
         return date_;
     }
 
+    CoordGeodetic ToGeodetic() const;
+
 private:
     Julian date_;
     Vector position_;