1subroutine geocentric(alat,elev,hlt,erad)
2
3  implicit real*8 (a-h,o-z)
4
5! IAU 1976 flattening f, equatorial radius a
6  f = 1.d0/298.257d0
7  a = 6378140.d0
8  c = 1.d0/sqrt(1.d0 + (-2.d0 + f)*f*sin(alat)*sin(alat))
9  arcf = (a*c + elev)*cos(alat)
10  arsf = (a*(1.d0 - f)*(1.d0 - f)*c + elev)*sin(alat)
11  hlt = datan2(arsf,arcf)
12  erad = sqrt(arcf*arcf + arsf*arsf)
13  erad = 0.001d0*erad
14
15  return
16end subroutine geocentric
17
18