	Program Eq2Galactic
	implicit real*8 (a-h,o-z)
	real*8 l(0:35,-9:9),b(0:35,-9:9)
	pi=3.141592653589d0
	open(1,'tabDif.txt',status='unknown')
	write(1,'(a/a,36(i4,2x))') 
     *'    Konwersja RA/DEC na wsp. galaktyczne (l,b) [w stopniach]'
     *,'Dec\RA',(i*10,i=0,35)

	do 1 i=0,35
	do 1 j=-9,9
	ra2000=i*10*pi/180
	d2000=j*10*pi/180
c	call Eq2Gal2k(ra2000,d2000,al,ab)
c	call Eq2Gal0(ra2000,d2000,2000d0,al1,ab1)
	call Eq2Gal(ra2000,d2000,2000d0,al1,ab1)
	call Gal2Eq(al1,ab1,RA2k,D2k,2000d0)
	dl=i*10-ra2k*180/pi
	if(dabs(dl).gt.180) dl=dl-dsign(360,dl)
	l(i,j)=dl*1d9
	b(i,j)=(d2000-d2k)*1d9*180/pi

1	continue
	print*,(l(i,0),b(i,0),i=0,35)
	write(1,'(i3,a3,36f6.1/3x,a3,36f6.1)') (j*10,'l',(l(i,j),i=0,35)
     *, ' b',(b(i,j),i=0,35),j=9,-9,-1)
c	write(1,'(i3,a3,36f6.1)') (j*10,'b',(b(i,j),i=0,35),j=9,-9,-1)
	end


	subroutine Gal2Eq(gl,b,RAepo,Depo,Epoch)

c Converts galactic coordinates (gl,b) given in radians 
c into the equatorial coordinates (RAepo,Depo) [radians]
c at the Julian epoch "Epoch".

	implicit real*8 (a-h,o-z)
	data B1950/2433282.42345905d0/,pi/3.141592653589d0/,ifirst/0/
c Standard galactic pole coordinates and center position at B1950 [degrees]: 
     * , RAp,Dp,gle/192.25d0,27.40d0,33d0/

	if(ifirst.ne.0) go to 1
	ifirst=ifirst+1
	pi_180=pi/180
	RApr=RAp*pi_180
	Dpr=Dp*pi_180
	gl0=gle*pi_180
	sinDp=dsin(Dpr)
	cosDp=dcos(Dpr)

1	sind=sinDp*dsin(b)+cosDp*dcos(b)*dsin(gl-gl0)
	D1950=dasin(sind)
	CdSdRA=dcos(b)*dcos(gl-gl0)*cosDp
	CdCdRA=dsin(b)-sind*sinDp
	RA1950=datan2(CdSdRA,CdCdRA)+RApr

	DJepo=2451545+365.25d0*(Epoch-2000d0)
	call PREnew(B1950,DJepo,RA1950,D1950,RAepo,Depo)

	end


	subroutine Eq2Gal2k(ra2000,d2000,gl,b)

c Procedure to convert the equatorial coordinates (ra2000,d2000) [radians]
c given for the J2000 epoch to the galactic coordinates (gl,b) [radians]

	implicit real*8 (a-h,o-z)
	data pi/3.141592653589d0/,ifirst/0/
     * ,RAp,Dp/192.8593357295d0,27.1282510272d0/,gl0/32.9319193453d0/
c Here above are galactic pole coordinates and location of the galactic
c center precessed to J2000 from standard values at B1950 (in degrees):
c RAp=192.25, Dp=27.40 and gl0=33.00. The last parameter, gl0, is 
c the galactic longitude of the equator.

	if(ifirst.ne.0) go to 1
	ifirst=ifirst+1
	pi_180=pi/180
	RApr=RAp*pi_180
	Dpr=Dp*pi_180
	sinDp=dsin(Dpr)
	cosDp=dcos(Dpr)
	gl0r=gl0*pi/180

1	sinb=sinDp*dsin(d2000)+cosDp*dcos(d2000)*dcos(ra2000-RApr)
	b=dasin(sinb)
	CbSdl=dsin(d2000)*cosDp-dcos(d2000)*sinDp*dcos(ra2000-RApr)
	CbCdl=dcos(d2000)*dsin(ra2000-RApr)
	gl=datan2(CbSdl,CbCdl)+gl0r
	gl=dmod(gl+pi+pi,pi+pi)

	end


	SUBROUTINE PREnew(DJE1,DJE2,ra1,d1,ra2,d2)

C Calculates general precession from DJE1 to DJE2 using new IAU theory
c Constants according to J.H.Lieske 1979, Astron. Astrophys. 73, 282.
	implicit real*8 (d,r)
	data DCSAR/4.848136812D-6/,J2000/2451545/,DC1/2306.2181D0/,
     * DC2/1.39656D0/,DC3/.30188D0/,DC4/.017998D0/,DC5/1.09468D0/
     * ,DC6/2004.3109D0/,DC7/-.8533D0/,DC8/-.42665D0/,DC9/-.041833D0/
	DT0=(Dje1-j2000)/36525
	DT=(Dje2-Dje1)/36525
	DTS=DT*DT
	DTC=DTS*DT
	DZETA=((DC1+(DC2-1.39d-4*dt0)*DT0)*DT+(DC3-3.44d-4*dt0)*DTS+
     * DC4*DTC)*DCSAR
	dzet =DZETA+((DC5-dc3+4.1d-4*dt0)*DTS+.000205*dtc)*DCSAR
	dth=((DC6+(DC7-2.17d-4*dt0)*DT0)*DT+(DC8-2.17d-4*dt0)*DTS+
     * DC9*DTC)*DCSAR

c convert Dj1 position to Dj2 using new IAU precession angles
	dc=cos(ra1+dzeta)
	ra2=dzet+datan2(sin(ra1+dzeta),-tan(d1)*sin(dth)+cos(dth)*dc)
	d2=asin(sin(d1)*cos(dth)+cos(d1)*sin(dth)*dc)
	ra2=dmod(ra2+2*3.141592653589d0,2*3.141592653589d0)
	end

	subroutine Eq2Gal0(RAepo,Depo,Epoch,gl,b)

c Converts equatorial coordinates (RAepo,Depo) given in radians
c for any Julian epoch "Epoch" (e.g. 2001.526) into the galactic
c coordinates (gl,b) in radians.

	implicit real*8 (a-h,l,o-z)
	data B1950/2433282.42345905d0/,pi/3.141592653589d0/,ifirst/0/
c Standard galactic center and pole coordinates at B1950 [degrees]: 
     * , RAp,Dp/192.25d0,27.40d0/,gl0/33d0/

	DJepo=2451545+365.25d0*(Epoch-2000d0)
	call PREnew(DJepo,B1950,RAepo,Depo,RA1950,D1950)

	if(ifirst.ne.0) go to 1
	ifirst=ifirst+1
	pi_180=pi/180
	RApr=RAp*pi_180
	sinDp=dsin(Dp*pi_180)
	cosDp=dcos(Dp*pi_180)
	gl0r=gl0*pi_180

1	sinb=sinDp*dsin(D1950)+cosDp*dcos(D1950)*dcos(RA1950-RApr)
	b=dasin(sinb)
	CbSdl=dsin(D1950)*cosDp-dcos(D1950)*sinDp*dcos(RA1950-RApr)
	CbCdl=dcos(D1950)*dsin(RA1950-RApr)
	gl=(datan2(CbSdl,CbCdl)+gl0r)
	gl=dmod(gl+pi+pi,pi+pi)

	end


	subroutine Eq2Gal(RAepo,Depo,Epoch,gl,b)

c Converts equatorial coordinates (RAepo,Depo) given in radians 
c from any Julian epoch "Epoch" (e.g. 2001.526) into the galactic
c coordinates (gl,b) [radians].

	implicit real*8 (a-h,o-z)
	data B1950/2433282.42345905d0/,pi/3.141592653589d0/,ifirst/0/
c Standard galactic center and pole coordinates at B1950 [degrees]: 
     * , RAp,Dp,gle/192.25d0,27.40d0,33d0/

	DJepo=2451545+365.25d0*(Epoch-2000d0)
	call PREnew(DJepo,B1950,RAepo,Depo,RA1950,D1950)

	if(ifirst.ne.0) go to 1
	ifirst=ifirst+1
	pi_180=pi/180
	RApr=RAp*pi_180
	Dpr=Dp*pi_180
c The angular distance: galactic center to equator
	gl0=gle*pi_180
	sinDp=dsin(Dpr)
	cosDp=dcos(Dpr)

1	sinb=sinDp*dsin(D1950)+cosDp*dcos(D1950)*dcos(RA1950-RApr)
	b=dasin(sinb)
	CbSdl=dsin(D1950)-sinb*sinDp
	CbCdl=dcos(D1950)*dsin(RA1950-RApr)*cosDp
	gl=datan2(CbSdl,CbCdl)+gl0
	gl=dmod(gl+pi+pi,pi+pi)

	end

