c	conversion cartesien-geographique et r‚ciproque
c	pour divers ellipsoides
c
	real*8 x(3),r(4),a,e,pi,f,u,v,w
	pi=3.1415926535898
	open(9,file='convert.out')

   1	write(*,'(a,$)')' WGS72(1),WGS84(2),WGS66(4),Clarke(3): '
	read(*,*) i
	if(i.eq.1) then
	a=6378135.
	f=1./298.26
	endif
	if(i.eq.3) then
	a=6378249.145
	f=1./293.465
	endif
	if(i.eq.4) then
	a=6378145.
	f=1./298.25
	endif
	if(i.eq.2) then
	a=6378137.
	f=1./298.257223563
	endif
	e=dsqrt(1-(1-f)*(1-f))

	write(*,*) a,e	
	write(*,'(a,$)') ' Cart>Geo (1), Geo>Cart (2) : '
	read(*,*) i
	if(i.eq.1) then
	write(*,'(a,$)') ' X = '
	read(*,*) x(1)
	write(*,'(a,$)') ' Y = '
	read(*,*) x(2)
	write(*,'(a,$)') ' Z = '
	read(*,*) x(3)
	write(9,'(3f15.4)') x(1),x(2),x(3)
	call cargeo(x,r,a,e)
	r(3)=r(3)*180.d0/pi
	r(4)=r(4)*180.d0/pi
	u=dint(r(3))
	v=dint(60.d0*(r(3)-u))
	w=3600.d0*(r(3)-u-v/60.d0)
	write(9,'(f6.0,f6.0,f12.4)') u,v,w
	write(*,'(a11,2f6.0,f12.3)') ' Longitude : ',u,v,w
	u1=dint(r(4))
	v1=dint(60.d0*(r(4)-u1))
	w1=3600.d0*(r(4)-u1-v1/60.d0)
	write(*,'(a12,2f6.0,f12.3)') ' Latitude : ',u1,v1,w1
	write(*,'(a10,f10.3)') ' Altitude ',r(2)
	write(9,'(2f6.0,f12.3,2f6.0,f12.3,f12.3)') u1,v1,w1,u,v,w,r(2)
	goto 1
	else
	write(*,'(a,$)') ' Latitude = '
	read(*,*) u,v,w
	r(4)=(u+(v+w/60.d0)/60.d0)*pi/180.d0
	write(*,'(a,$)') ' Longitude = '
	read(*,*) u1,v1,w1
	r(3)=(u1+(v1+w1/60.d0)/60.d0)*pi/180.d0
	write(*,'(a,$)') ' Altitude = '
	read(*,*) r(2)
	write(9,'(2f6.0,f12.3,2f6.0,f12.3,f12.3)') u,v,w,u1,v1,w1,r(2)
	call geocar(x,r,a,e)
	write(*,*) ' Coordonnees tridimensionnelles'
	write(*,'(3f15.3)') x(1),x(2),x(3)
	write(9,'(3f15.3)') x(1),x(2),x(3)
	goto 1
	endif

	end

c ----------------------------
	subroutine cargeo(x,r,a,e)
c ----------------------------
c
c	passage des coordonnees cartesiennes aux geographiques
c	x : tableau des coodonnees cartesiennes
c	r : R,h,longitude,latitude
c
	real*8 x(3),r(4),a,e,p,n,pi
	pi=3.1415926535898
c
	r(1)=dsqrt(x(1)**2+x(2)**2+x(3)**2)
	r(3)=datan(x(2)/x(1))
	if(x(1).le.0) r(3)=r(3)+pi
	if(r(3).ge.pi) r(3)=r(3)-2.*pi
	p=dsqrt(x(1)**2+x(2)**2)
	r(4)=datan(x(3)/p/(1.d+00-a*e*e/r(1)))
	do 1 i=1,5
	  n=a/dsqrt(1.d+00-(e*dsin(r(4)))**2)
	  r(2)=p/dcos(r(4))-n
   1	  r(4)=datan(x(3)/p/(1.d+00-n*e*e/(n+r(2))))
	end
c
c ----------------------------
	subroutine geocar(x,r,a,e)
c ----------------------------
c
c	passage des coordonnees geographiques aux cartesiennes
c	x : tableau des coordonnees cartesiennes
c	r : R, h, longitude, latitude
c
	real*8 x(3),r(4),a,e,n
	n=a/dsqrt(1.d+00-(e*dsin(r(4)))**2)
	x(1)=(n+r(2))*dcos(r(4))*dcos(r(3))
	x(2)=(n+r(2))*dcos(r(4))*dsin(r(3))
	x(3)=(n*(1.d+00-e*e)+r(2))*dsin(r(4))
	end
