function earth_rad_geodetic,geodetic_lat

;-----------------------------------------------------------------------
; Computes the radius (distance from Earth-center) of the 
; ellipsoidal Earth at a specified geodetic latitude.
;
; Input        geodetic_lat = geodetic latitude (degrees)
;
; On return    earth_rad = radius of the ellipsoidal Earth (km)
;                          (distance from earth-center)
;
; If |lat|>90, then earth_rad = -1 is returned              
;
; Note: Geocentric and geodetic latitude differ.  The geocentric
; latitude is the angle between the  equatorial plane and the Earth-
; centered radius vector.  The geodetic latitude is the "map" 
; latitude, defined as the angle between the equatorial plane and
; the local vertical.  The difference is always less than 0.2 degrees.
;
; Marty McHugh,   GATS, Inc.  Dec 1999 (IDL conversion by M. Hervig, June 2007)
;-----------------------------------------------------------------------

      a   = 6378.137d0     ; equatorial radius (km)
      b   = 6356.752d0     ; polar radius (km)
      ba2 = (b*b)/(a*a) 

;-----------------------------------------------------------------------

; test the input
      
      if (abs(geodetic_lat) gt 90.0d0) then begin
      
        print,'error: |geodetic latitude| > 90'
        print,'subroutine: earth_rad_geodetic.pro'
        stop
	    
      endif


; calculate radius to point on ellipse

      cos_lat = cos(geodetic_lat)
      sin_lat = sin(geodetic_lat)

      c = a/sqrt((cos_lat*cos_lat) + ba2*(sin_lat*sin_lat))
      s = ba2*c
      x = c*cos_lat
      y = s*sin_lat

      earth_rad_geodetic = sqrt(x*x+y*y)
 
      return,earth_rad_geodetic
      end