subroutine geodetic(r, z, phi, h)
!
!     Given Cartesian coordinates r (= sqrt(x**2+y**2)) and z, where
!     x, y, z are rectangular ECI coordinates, return geodetic 
!     latitude phi (radians), and geodetic height h (km).
!
!     B. Knapp, 2001-01-10
!
!     Reference: T. Fukushima, J. Geod. 73 (1999), 603-610.
!
C
C     RCS DATA
C     
C     $Header$
C     
C     $Log$
C
C
      implicit none
!
!     Input
      real*8 r, z
!
!     Output
      real*8 phi, h
!
!     Earth ellipsoid parameters
      real*8 RE, F, E_SQ, EPRIME, C
      parameter (RE=6378.135d0, F=1.d0/298.257d0, E_SQ=F*(2.d0-F),
     &    C=RE*E_SQ)
!
!     Constants
      real*8 TOL
      parameter (TOL=1.0d-12)
!
!     Local variables
      real*8 zprime, u, v, tm, fm, t, dt, r4, u3, t2, t2m, t2p, ept2
!
!     Ellipsoid parameter e-prime (should be a constant, but not all
!     compilers allow function calls within a parameter statement)
      EPRIME=sqrt(1.d0-E_SQ)
!
!     Quartic function, f(t) = r*t^4 + u*t^3 + v*t - r = 0
      zprime = abs(z)*EPRIME
      u = 2.d0*(zprime-C)
      v = 2.d0*(zprime+C)
!
!     Find Newton-Raphson starting point
      tm = (C-zprime)/r
      if (0.d0 .lt. tm .and. tm .lt. 1.d0) then
         fm = ((r*tm+u)*tm**2+v)*tm-r
      else
         fm = 0.d0
      endif
      if (tm .le. 0.d0 .or. fm .lt. 0.d0) then
         t = (r-C+zprime)/(r-C+2.d0*zprime)
      else
         t = r/(C+zprime)
      endif
!
!     Iterate to convergence
      r4 = r*4.d0
      u3 = u*3.d0
      dt = 1.d0
      t2 = t**2
      do while (abs(dt) .gt. TOL)
         dt = (((r*t+u)*t2+v)*t-r)/((r4*t+u3)*t2+v)
         t = t-dt
         t2 = t**2
      enddo
!
!     Calculate results
      t2m = 1.d0-t2
      t2p = 1.d0+t2
      ept2 = 2.d0*EPRIME*t
      phi = atan2(t2m,ept2)
      if (z .lt. 0.d0) then
         phi = -phi
      endif
      h = (r*ept2+abs(z)*t2m-RE*EPRIME*t2p)/sqrt(t2p**2-4.d0*E_SQ*t2)
!
      return
      end