subroutine pe2xy(mjd,dpsi,deps,X,Y)

C-------------------------------------------------------------------------C
C									  C
C       Conversion of total dpsi, deps to X, Y                            C
C									  C
C-------------------------------------------------------------------------C
C	 								  C
C	References :                                                      C
C               IERS Conventions 2000, chap. 5                            C
C	 								  C
C-------------------------------------------------------------------------C
C	 								  C
C	Input :                                                           C
C		mjd : time in mjd                                         C
C		dpsi, deps : total nutation angles in arcsec              C
C	 		     = model [IAU1980] + dp, de [IAU1980]         C
C	 		    or model [IAU2000] + dp, de [IAU2000]         C
C	Output :                                                          C
C               X, Y : CIP coordinates in the ICRF in arcsec	          C
C	 								  C
C-------------------------------------------------------------------------C
C	 								  C
C	Subroutine interface : none                                       C
C                                                                         C
C-------------------------------------------------------------------------C
C	 								  C
C	Author : Sebastien Lambert (Sebastien.Lambert@obspm.fr)		  C
C	Last modified : 23/05/03	 				  C
C	 								  C
C-------------------------------------------------------------------------C

	implicit none

C	Input variables

	double precision               :: mjd,dpsi,deps

C 	Output variables

	double precision               :: X,Y
	
C	Other variables
	
	double precision               :: deps1,dpsi1
	double precision               :: psiA,omeA,epsA,chiA
	double precision               :: ksi0,eta0,dalpha0
	double precision               :: psi,omega
	double precision               :: Xb,Yb
	double precision               :: jc,eps0
        double precision, parameter    :: a2r=4.84813681109535993D-6
        double precision, parameter    :: r2a=2.06264806247096355D+5

C	Julian centuries from J2000.0

	jc=(mjd-51544.5D0)/36525.D0

C	Lieske et al. (1977) expressions in arcsec and converted in radians

	eps0=84381.448D0

	psiA=5038.47875D0*jc
     .      -1.07259D0*jc*jc
     .      -0.001147D0*jc*jc*jc

        omeA=eps0
     .      -0.02524D0*jc
     .      +0.05127D0*jc*jc
     .      -0.007726D0*jc*jc*jc

        epsA=eps0
     .      -46.8402D0*jc
     .      -0.00059D0*jc*jc
     .      +0.001813D0*jc*jc*jc

        chiA=10.5526D0*jc
     .      -2.38064D0*jc*jc
     .      -0.001125D0*jc*jc*jc

     	eps0=eps0*a2r
	psiA=psiA*a2r
	omeA=omeA*a2r
	epsA=epsA*a2r
	chiA=chiA*a2r
	
C	Right ascension of the mean equinox of J2000.0 in the ICRS in radians

	dalpha0=-0.01460D0*a2r
	
C	Celestial pole offsets at J2000.0 in radians

	ksi0=-0.0166170D0*a2r
	eta0=-0.0068192D0*a2r

C	Total nutation angles dpsi, deps in radians

	dpsi=dpsi*a2r
	deps=deps*a2r

C	Chapter 5, equation (22)

	dpsi1=(dpsi*dsin(epsA)*dcos(chiA)-deps*dsin(chiA))/dsin(omeA)
	deps1=dpsi*dsin(epsA)*dsin(chiA)+deps*dcos(chiA)

C	Chapter 5, equation (21)

	omega=omeA+deps1
	psi=psiA+dpsi1

C	Chapter 5, equation (20)

	Xb=dsin(omega)*dsin(psi)
	Yb=-dsin(eps0)*dcos(omega)+dcos(eps0)*dsin(omega)*dcos(psi)

C	Chapter 5, equation (18)

	X=Xb+ksi0-dalpha0*Yb
	Y=Yb+eta0+dalpha0*Xb

C	Conversion radians in arcsec

	X=X*r2a
	Y=Y*r2a

	end