From: 103706.1551@compuserve.com (Pamela  M. Foote)
Newsgroups: sci.math
Subject: Re: Spherical / Oblate Earth Coordinate Conversions
Date: 17 Feb 1996 00:53:00 GMT

A posting from DragonCS@aol.com asked
> I am looking for  algorithms which will convert spherical earth
> eci coordinates to oblate earth eci coordinates and the inverse.
> Also algorithms to convert oblate earth eci to goedetic latitude,
> longitude, and altitude.

This stuff is straightforward but fairly hard to find.  Look under
"geodesy" for references.  For an exhaustive writeup, find a copy
of The Defense Mapping Agency World Geodetic System 1984 (WGS-84)
report.

A synopsis without derivations:

The defining constants (and their WGS-84 values) are
   a = 6378137 meters = equatorial radius 
   b = polar radius
   f = (a-b)/a = 1/298.257223563 = flattening 
   e2 = f*(2-f) = eccentricity squared
   w = 0.72921151467e-4  rad s^-1 = Earth rate

 To go from geodetic Latitude, longitude, altitude = (L,g,h)
 to Earth centered rotating (x,y,z)
    N = a*/sqrt(1-e2*sin(L)^2)
    x = (N + h).*cos(L)*cos(g)
    y = (N + h).*cos(L)*sin(g)
    z = ((1-e2)*N + h)*sin(L)

The inverse is trickier, since it requires the root of a cubic
equation.  Although there _is_ a closed solution, in practice an
iterative computation is easier and more efficient, namely
   g = atan2(y,x))
Iterative solution for t=tan(L)
   k = (1-e2);
   r = sqrt(x^2 + y^2);

   t = zr;
   for i=1:niter
      t = z/r + t*a*e2/(r*sqrt(1+k*t^2));
   end;
For altitudes in the atmosphere 2 or 3 iterations are plenty.  Further
out, you
need more.
Convert to actual coordinates
   L = atan(t)
   h = sqrt(1+t^2)*(r - a*/sqrt(1+k*t^2))

To go back and forth between ECR and ECI (Earth centered interial) you
need to
pick a reference time and compute T=time since then, and rotate the
cartesian
coordinates.
   C = [ cos(w*T)  sin(w*T)  0
        -sin(w*T)  cos(w*T)  0
            0         0      1 ]

  R_ecr=C*R_eci
  R_eci=transpose(C)*R_ecr
where R is the x,y,z cartesian vector.

For a spherical Earth, set f=0.  (So N=a and tan(L)=z/r, etc.)
Use the cartesian coordinates as a bridge.

            hope this helps -- MFoote