        character*30 f1,line,c*1
        character*12 ac,dc
        real*8 a8(8)
        character*80 head(72)

        common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
        pi=4.*atan(1.0)
        cx=pi/12.
        cy=pi/180.
        phi=(40.+23/60.+36/3600.)*pi/180.       ! XL latitude
        sphi=sin(phi)
        cphi=cos(phi)

        k=iargc()
        if(k.eq.0)then
          write(*,*)'               usage: xy2az filename'
          write(*,*)'                                      97,6'
          write(*,*)
          write(*,*)'       context in file:'
          write(*,*)'        fits_filename'
          write(*,*)'        1 or 2049 (iraf data, put 2049)'
          write(*,*)'        x1 y1'

          write(*,*)'        x2 y2'
          write(*,*)'       ......'
          stop
        endif
        call getarg(1,f1)
        open(1,file=f1,status='old',iostat=i)
        if(i.ne.0)stop ' file not found!'
        read(1,'(a)')f1
        open(2,file=f1,status='old',access='direct',recl=5760)
        read(2,rec=1)head
        close(2)
        call getpar(head)                               ! get s0,iy,im,id,ut
        ipos=indexpos(head,'A81     ')
        if(ipos.eq.73)stop 'do position calibration first!'
        read(1,'(a)')f1
        read(f1(1:),*,end=11)ip
11      do 10 i=1,8
10      read(head(ipos+i-1)(11:),*)a8(i)

	ipos=indexpos(head,'RA      ')
	line=head(ipos)(12:31)
	line(20:20)=char(0)
	do 15 i=1,19
	if(line(i:i).eq.':')line(i:i)=','
15	if(ichar(line(i:i)).eq.39)line(i:i)=' '
	read(line(1:),*,err=155)i,j,x
155	alpha_c=i+j/60.+x/3600.
	a8(7)=alpha_c*cx

	ipos=indexpos(head,'DEC     ')
	line=head(ipos)(12:31)
	line(20:20)=char(0)
	c=' '
	do 16 i=1,19
	if(line(i:i).eq.'-')c='-'
	if(line(i:i).eq.'-')line(i:i)=' '
	if(line(i:i).eq.':')line(i:i)=','
16	if(ichar(line(i:i)).eq.39)line(i:i)=' '
	read(line(1:),*,err=156)i,j,x
156	delta_c=i+j/60.+x/3600.
	if(c.eq.'-')delta_c=-delta_c
	a8(8)=delta_c*cy

c	write(*,*)s0,iy,im,id,ut

20      read(1,*,err=100,end=100)x,y
        x1=x
        y1=y
        if(ip.eq.2049)then
          x1=2049.-x
          y1=2049.-y
        endif
        call xy_rad(x1,y1,a8,alpha,delta)
        call toms(alpha,ac,1)
        call toms(delta,dc,0)
c        write(*,2)x,y,ac,dc
2       format(2f12.3,'---> ',2a)
        call ad_azl(alpha,delta,x1,y1)
        write(*,1)x,y,x1,y1
1       format(2f12.3,'---> ',2f12.3)
        goto 20
100     end

        function indexpos(head,f1)
        character*80 head(72),f1*8
        byte bell
        data bell/7/
        do 10 indexpos=1,72
10      if(head(indexpos)(1:8).eq.f1)return
        if(indexpos.eq.73)write(*,'(1x,a)')bell
        end

        subroutine xy_rad(x,y,a8,alpha,delta)
c 2048*2048 ccd's X Y (in pixel) -------------> alpha delta (in h,d)
        real x,y,alpha,delta
        real*8 a8(8),pi,xi,xn,rc,dc,cd,td,alph,delt
        pi=4.*datan(1d0)
        xi=a8(1)*x+a8(3)*y+a8(5)
        xn=a8(2)*x+a8(4)*y+a8(6)
        rc=a8(7)                                ! center of ccd (in rad.)
        dc=a8(8)
        cd=dcos(dc)
        td=dtan(dc)
        alph=datan(xi/cd/(1.-xn*td))
        delt=datan((xn+td)*dcos(alph)/(1.-xn*td))
        alph=alph+rc
        if(abs(dc-delt).gt.2.)then
          delt=-delt
          alph=alph+pi
        endif
        alpha=alph*12./pi
        delta=delt*180./pi
        end



        subroutine ad_azl(alpha,delta,a,h)     ! hd,dd
        common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
        hour=(s0-alpha)*cx
        delt=delta*cy
        s_delta=sin(delt)
        c_delta=cos(delt)
        s_hour=sin(hour)
        c_hour=cos(hour)
        sinh=sphi*s_delta+cphi*c_delta*c_hour
        cosh=sqrt(1.-sinh*sinh)
        sinA=c_delta*s_hour/cosh
        cosA=(sphi*c_delta*c_hour-cphi*s_delta)/cosh
        h=atan(sinh/cosh)/cy
        a=atan(sinA/cosA)
        if(cosA.lt.0.)a=a+pi
        if(a.lt.0.)a=a+pi+pi
        if(a.gt.pi+pi)a=a-pi-pi
        a=a/cy
        end



        subroutine getpar(head)
c get s0 iy im id ut
        character*80 head(72)
        real*8 d1,d
        common /xinglong/pi,cx,cy,sphi,cphi,s0,iy,im,id,ut
        ipos=indexpos(head,'DATE-OBS')
        read(head(ipos)(22:),'(i2,1x,i2,1x,i2)')id,im,iy        ! read obs date
        iy=iy+1900
        if(iy.lt.1950)iy=iy+100
        ipos=indexpos(head,'TIME    ')
        read(head(ipos)(20:),'(i2,1x,i2,1x,f4.1)')i,j,x
        ipos=indexpos(head,'EXPOSURE')
        read(head(ipos)(20:),*)k                                ! exp time
        ut=i+j/60.+x/3600.+k/7200.                      ! center UT time, in h
        call fjd(d1,1992,12,31)                 ! 1992,12,31
        d2=(6+38/60.+40.1954/3600.)/24.         ! 6 38 40.1954
        x1=7+50/60.+18.344/3600.                ! 7:50:18.344 lamda of XL
        call fjd(d,iy,im,id)
        d=(d-d1)*1.0027379093d0+d2
        d=dmod(d,1d0)
        if(d.lt.0d0)d=d+1d0
        s0=d*24d0
        s0=s0+x1+ut*1.002738
        end

        subroutine fjd(d,iy,im,id)
        real*8 d
        i=(im+9)*0.09
        i=(i+iy-1900)*1.75+0.01
        j=30.56*im
        j=(iy-1950)*367+id+j-i
        d=2433338.5d0+j
        end


