[Proj] Geocentric vs. Geodetic latitude

[Glynn Clements]

Frank Warmerdam wrote:

> I was recently suprised to encounter locations identified as being
> geocentric that were give as a latitude and longitude.  I had expected them
> to be x/y/z values in meters.
> 
> A client of mine, better versed in such things explained that there is in
> fact a distinction between geocentric and geodetic lat/long.  He wrote:
> 
>  >"The angle L' is called "geocentric latitude" and is defined as the
>  > angle between the equatorial plane and the radius from the geocenter.
>  >
>  > The angle L is called "geodetic latitude" and is defined as the angle
>  > between the equatorial plane and the normal to the surface of the
>  > ellipsoid.  The word "latitude" usually means geodetic latitude.  This
>  > is the basis for most of the maps and charts we use.  The normal to the
>  > surface is the direction that a plumb bob would hang were it not for
>  > local anomalies in the earth's gravitational field."
> 
> I need to implement some code to convert geocentric latitude to geodetic
> latitude.  I think I can do from the above description, but I am wondering
> if anyone can confirm the above description, provide a forumla or most
> importantly provide some sample latitudes in both systems I can check my
> work against.

Given the parametric equation for an ellipse:

        x = a.cos(t)
        y = b.sin(t)

the tangent vector is:

        dx/dt = -a.sin(t)
        dy/dt =  b.cos(t)

and thus the outward normal is:

        nx = b.cos(t)
        ny = a.sin(t)

By the above definitions, the geocentric latitude L' is given by:

        tan(L') = y / x
                = b.sin(t) / a.cos(t)
                = (b/a).tan(t)

while the geodetic latitude L is given by:

        tan(L)  = ny / nx
                = a.sin(t) / b.cos(t)
                = (a/b).tan(t)

Thus the ratio of the two is:

        tan(L) / tan(L')
                = (a/b).tan(t) / (b/a).tan(t)
                = (a/b) / (b/a)
                = a²/b²

So:
        L = atan((a²/b²).tan(L'))
and:
        L' = atan((b²/a²).tan(L))

where a is the equatorial axis and b the polar axis.

-- 
Glynn Clements <glynn.clements at virgin.net>