function geodetic2geocentric,geodetic
;geodetic=geodetic latitude in radians

f=1/298.257222101d
esquared=2*f-f^2
;GRS80 ellipsoid, use 1/298.257223563 for WGS84

geocentric=geodetic
geocentric[0]=atan((1-esquared)*tan(geodetic[0]))
return,geocentric

end