[Proj] Geocentric vs. Geodetic latitude

[Frank Warmerdam]

Glynn Clements wrote:
> 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,

Well done!  I have tried (in python):

import math

a = 6378137.0
b = 6356752.3142

gc_lat = 32.98333333333333
gc_lat_rad = gc_lat * math.pi / 180.0

gd_lat_rad = math.atan( ((a*a)/(b*b)) * math.tan(gc_lat_rad) )

gd_lat = gd_lat_rad * 180.0 / math.pi

print gc_lat
print gd_lat

The results is 33.1593159397 as opposed to the  33.1593160 (rounded)
I was getting.  There would appear to be some difference in the 8th
decimal place (the 10th digit) but I don't know why that is.

In any event, your solutions is much briefer than mine, and I will replace
my solution with yours.

Thank you,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent