//////////////////////////////////////////////////////////////////////
//
// Convert a cartexian XYZ coordinate to a geodetic lat/lon/alt. This
// is a "recursion relation". In essence, it iterates over the 2D
// part of sgGeodToCart refining its approximation at each step. The
// MAX_LAT_ERROR threshold is picked carefully to allow us to reach
// the full precision of an IEEE double. While this algorithm might
// look slow, it's not. It actually converges very fast indeed --
// I've never seen it take more than six iterations under normal
// conditions. Three or four is more typical. (It gets slower as the
// altitude/error gets larger; at 50000m altitude, it starts to need
// seven loops.) One caveat is that at *very* large altitudes, it
// starts making very poor guesses as to latitude. As altitude
// approaches infinity, it should be guessing with geocentric
// coordinates, not "local geodetic up" ones.
//
//////////////////////////////////////////////////////////////////////
void
sgCartToGeod ( const Point3D& CartPoint , Point3D& GeodPoint )
{
static const double MAX_LAT_ERROR = 8.881784197001252e-16;
double x = CartPoint[X];
double y = CartPoint[Y];
double z = CartPoint[Z];
double r = sqrt(x*x + y*y);
double lon, lat, alt;
double lat1, lat2 = localLat(r, z);
double r2, z2, dot;
lon = atan2(y, x);
do
{
lat1 = lat2;
double upr = cos(lat1);
double upz = sin(lat1);
surfRZ(upr, upz, &r2, &z2);
r2 = r - r2;
z2 = z - z2;
dot = r2*upr + z2*upz;
r2 = r - dot * upr;
z2 = z - dot * upz;
lat2 = localLat(r2, z2);
} while (fabs(lat2 - lat1) > MAX_LAT_ERROR);
lat = lat1;
double dr = r - r2;
double dz = z - z2;
double altsign = (dot > 0) ? 1 : -1;
alt = altsign * sqrt(dr*dr + dz*dz);
GeodPoint[Lat] = lon * SG_RADIANS_TO_DEGREES;
GeodPoint[Lon] = lat * SG_RADIANS_TO_DEGREES;
GeodPoint[Alt] = alt * SG_METER_TO_FEET;
} // sgCartToGeod()