implicit real*8 (a-h,o-z)
      parameter(a=6378137.d0,fm1=298.257224d0,w=0.72921151467d-4)
      real*8 lat,lon,latdeg,latmin,latsec,londeg,lonmin,lonsec
      logical xyzgeo
      pi=4.d0*atan(1.d0)
      print*,'Input x, y and z'
      read(5,*)x,y,z
      if(.not.xyzgeo(x,y,z,lat,lon,h))then
         print*,'Conversion failed'
      else
         lat=180.d0*lat/pi
         lon=180.d0*lon/pi
         print*,'lat = ',lat,'  lon = ',lon,'  h = ',h
         call dms(lat,lon,latdeg,latmin,latsec,londeg,lonmin,lonsec)
         print*,'Lat ',latdeg,' deg, ',latmin,' min',latsec,' sec'
         print*,'Lon ',londeg,' deg, ',lonmin,' min',lonsec,' sec'
      endif
      stop
      end
c
c
      subroutine geoxyz(gla,glo,eht,x,y,z)
*** convert geodetic lat, long, ellip ht. to x,y,z
*** input units of radians, output in meters
******NOTE: e2 is from wgs84  ***********CHECK THIS
      implicit double precision(a-h,o-z)
      parameter(a=6378137.d0,e2=0.00669437999013d0)
      sla=dsin(gla)
      cla=dcos(gla)
      w2=1.d0-e2*sla*sla
      w=dsqrt(w2)
      en=a/w
      x=(en+eht)*cla*dcos(glo)
      y=(en+eht)*cla*dsin(glo)
      z=(en*(1.d0-e2)+eht)*sla
      return
      end
c
c
      logical function xyzgeo(x,y,z,glat,glon,eht)
*** convert x,y,z into geodetic lat, lon, and ellip. ht
*** input units of meters, output in radians
******NOTE: e2 is from wgs84  ***********CHECK THIS
*** ref: eq a.4b, p. 132, appendix a, osu #370
*** ref: geom geod notes gs 658, rapp
      implicit double precision(a-h,o-z)
      parameter(maxint=10,tol=5.d-15)
      parameter(a=6378137.d0,e2=0.00669437999013d0)
      ae2=a*e2
*** compute initial estimate of reduced latitude  (eht=0)
      p=dsqrt(x*x+y*y)
      icount=0
      tgla=z/p/(1.d0-e2)
*** iterate to convergence, or to max # iterations
    1 if(icount.le.maxint) then
        tglax=tgla
        tgla=z/(p-(ae2/dsqrt(1.d0+(1.d0-e2)*tgla*tgla)))
        icount=icount+1
        if(dabs((tgla-tglax)/tgla).gt.tol) go to 1
*** convergence achieved
        xyzgeo=.true.
        glat=datan(tgla)
        slat=dsin(glat)
        clat=dcos(glat)
        glon=datan2(y,x)
        w=dsqrt(1.d0-e2*slat*slat)
        en=a/w
        if(dabs(glat).le.0.7854d0) then
          eht=p/clat-en
        else
          eht=z/slat-en+e2*en
        endif
        glon=datan2(y,x)
*** too many iterations
      else
        xyzgeo=.false.
        glat=0.d0
        glon=0.d0
        eht=0.d0
      endif
      return
      end
c
c
      subroutine dms(lat,lon,latdeg,latmin,latsec,londeg,lonmin,lonsec)
      implicit real*8 (a-h,o-z)
      real*8 lat,lon,latdeg,latmin,latsec,londeg,lonmin,lonsec
      latdeg=dfloat(int(lat))
      rem=lat-latdeg
      latmin=dfloat(int(rem*60))
      rem=rem-latmin/60.d0
      latsec=rem*3600
      londeg=dfloat(int(lon))
      rem=lon-londeg
      lonmin=dfloat(int(rem*60))
      rem=rem-lonmin/60.d0
      lonsec=rem*3600
      return
      end