	program COORD
	implicit real*8 (a-h,p-z)
	character rec*15
	pi=3.141592653589d0
	print*,'   Coordinates conversion program (x,y,z) -> (r,fi,lon)'
c	print*,' Rectangular coordinates: '
c	read(*,*) x,y,z
c Wspolrz. RT32 ocenione z polozenia wzgl. RT15 (KB)
c	x=3638553.54d0
c	y=1221965.38d0
c	z=5077034.00d0
	x=3638558.d0
	y=1221967.d0
	z=5077041.d0
1	print'(a,f13.3,a)',' Enter the X coordnate in m [ENTER=',x,']: '
	read(*,'(a)') rec
	if(rec.eq.' ') go to 2
	read(rec,*) x
	go to 1
2	print'(a,f13.3,a)',' Enter the Y coordnate in m [ENTER=',y,']: '
	read(*,'(a)') rec
	if(rec.eq.' ') go to 3
	read(rec,*) y
	go to 2
3	print'(a,f13.3,a)',' Enter the Z coordnate in m [ENTER=',z,']: '
	read(*,'(a)') rec
	if(rec.eq.' ') go to 4
	read(rec,*) z
	go to 3

4	r=sqrt(x*x+y*y)
	rho=sqrt(x*x+y*y+z*z)
	along=datan2(y,x)*180/pi
	gclat=datan(z/r)*180/pi
	call geod(r,z,gdlat,h)
	gdlat=gdlat*180/pi

	print'(6(a7,f11.2/),3(a7,f11.5/))',
     * ' x =',x,' y =',y,' z =',z,' r =',r,' R =',rho,' h =',h,
     *' Long',along,' gdLat',gdlat,' gcLat',gclat

	Print*,' Write these results to a file [ENTER for none]: '
	read(*,'(a)') rec
	if(rec.eq.' ') stop
	open(1,file=rec)
	write(1,'(6(a7,f11.2/),3(a7,f11.5/))')
     * ' x =',x,' y =',y,' z =',z,' r =',r,' R =',rho,' h =',h,
     *' Long',along,' gdLat',gdlat,' gcLat',gclat

	end

      subroutine GEOD(r,z,fi,h)
c  Program to transform Cartesian to geodetic coordinates based 
c  on the exact solution (K.M. Borkowski,1989, Bull. Geod. 63, 50-56)
c  Input :  r, z = equatorial [m] and polar [m] components
c  Output:  fi, h = geodetic coord's (latitude [rad], height [m])
      implicit real*8(a-h,o-z)
c  GRS80 ellipsoid:  semimajor axis (a) and inverse flattening (fr)
      data a,fr /6378137.d0,298.257222101d0/
      b = dsign(a - a/fr,z)
      E = ((z + b)*b/a - a)/r
      F = ((z - b)*b/a + a)/r
c  Find solution to: t**4 + 2*E*t**3 + 2*F*t - 1 = 0
      P = (E*F + 1.)*4d0/3.d0
      Q = (E*E - F*F)*2.d0
      D = P*P*P + Q*Q
           if(D.ge.0d0) then
      s = dsqrt(D) + Q
      s = dsign(dexp(dlog(dabs(s))/3d0),s)
      v = P/s - s
c Improve the accuracy of numeric value of v
      v = -(Q + Q + v*v*v)/(3*P)
           else
      v = 2.*dsqrt(-P)*dcos(dacos(Q/P/dsqrt(-P))/3.)
           endif
      G = .5*(E + dsqrt(E*E + v))
      t = dsqrt(G*G + (F - v*G)/(G + G - E)) - G
      fi = datan((1.d0 - t*t)*a/(2*b*t))
      h = (r - a*t)*dcos(fi) + (z - b)*dsin(fi)
      return
      end
