	Program Ecliptic_Equator
c this is a test program
	implicit real*8 (a-h,o-z)
	data pi/3.141592653589793d0/
	RA=17.1/12*pi
	dec=-37.9/180*pi
	obliq=eps(2451045.D0)
	call Eq2Ecl(obliq,RA,dec,alambda,beta)

	print*,alambda*180/pi,beta*180/pi

	call Ecl2Eq(obliq,alambda,beta,RA1,dec1)
	print'(a,2f15.7)',' Error in round conversion (RA,dec)["]:'
     *,(RA1-RA)/pi*180*3600,(dec1-dec)/pi*180*3600

	end


	subroutine Ecl2Eq(epsilon,elng,eltd,RA,dec)
c Converts ecliptic coordinates 
c  elng - longitute, and
c  eltd - latitude
c to equatorial coordinates 
c   RA  - right ascension
c   dec - declination
c (all angles in radians), for given
c epsilon - the obliquity of equator with respect to ecliptic.
	implicit real*8 (a-h,o-z)
	data pi/3.141592653589793d0/
	sep=dsin(epsilon)
	cep=dsqrt(1-sep*sep)
	dec=dasin(cep*dsin(eltd)+sep*dcos(eltd)*dsin(elng))
	RA=datan2(cep*dsin(elng)-sep*dtan(eltd),dcos(elng))
	if(RA.lt.0d0) RA=RA+pi+pi
	end

	subroutine Eq2Ecl(epsilon,RA,dec,elng,eltd)
	implicit real*8 (a-h,o-z)
	data pi/3.141592653589793d0/
	sep=dsin(epsilon)
	cep=dsqrt(1-sep*sep)
	eltd=dasin(cep*dsin(dec)-sep*dcos(dec)*dsin(RA))
	elng=datan2(cep*dsin(RA)+sep*dtan(dec),dcos(RA))
	if(elng.lt.0d0) elng=elng+pi+pi
	end

	double precision function eps(dj)
	implicit real*8 (a-h,o-z)
	data pi/3.141592653589793d0/
	t=(dj-2451545.D0)/36525d0
	eps=(84381.448d0+t*(-46.815d0+t*(-.00059d0+t*.001813d0)))*
     *pi/(3600*180)
	end

	subroutine EqR2Ecl(epsilon,xyz,elng,eltd,exyz)
c Converts rectangular equatorial coordinates to ecliptic coordinates 
c rectangular and angular (longitute and latitude, in radians)
	implicit real*8 (a-h,o-z)
	dimension xyz(3),exyz(3)
	data pi/3.141592653589793d0/
	sep=dsin(epsilon)
	cep=dsqrt(1-sep*sep)
c*	Rp=dsqrt(Xyz(1)**2 + xYz(2)**2 + xyZ(3)**2)
	ye=xYz(2)*CEp + xyZ(3)*SEp
	ze=xyZ(3)*CEp - xYz(2)*SEp
c*	eltd = dASIN(ze/Rp)
	eltd = dAtan(ze/dsqrt(Xyz(1)**2+ye**2))
	elng = dMOD(dATAN2(ye,Xyz(1)) + 2*pi,2*pi)
c*	exyz(1)=Rp*dcos(eltd)*dcos(elng)
	exyz(1)=xyz(1)
c*	exyz(2)=Rp*dcos(eltd)*dsin(elng)
	exyz(2)=ye
c*	exyz(3)=Rp*dsin(eltd)
	exyz(3)=ze
	end

