//////////////////////////////////////////////////////////////////////
//
// 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()