	character*30 f1
	character*12 ac,dc
	real*8 a8(8)
	character*80 head(72)
	k=iargc()
	if(k.eq.0)then
	  write(*,*)'               usage: xy2ad filename'
	  write(*,*)'                                      06,03'
	  write(*,*)'                  epoch = 2000.0'
	  write(*,*)
	  write(*,*)'       context in file:'
	  write(*,*)'        fits_filename'
	  write(*,*)'        1 or 2049 or 513,[epoch] (iraf data, put >1)'
	  write(*,*)'        x1 y1                    (4k*4k data, put 1)'
	  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)
	ipos=indexpos(head,'A81     ')
	if(ipos.eq.73)stop 'do position calibration first!'
	epoch=2000.0
	read(1,'(a)')f1
	ix=0	
	read(f1(1:),*,end=11)ip,x
	if(x.ne.0.)epoch=x
	ix=1
11      do 10 i=1,8
10      read(head(ipos+i-1)(11:),*)a8(i)
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
	if(ip.eq.513)then
	x1=513.-x
	y1=513.-y
	endif
        call xy_rad(x1,y1,a8,alpha,delta)
	if(ix.eq.0)then
	call toms(alpha,ac,1)
	call toms(delta,dc,0)
	write(*,1)x,y,ac,dc
1	format(2f12.3,'---> ',2a)
	else
	call astprs(alpha,delta,2000.,x1,y1,epoch)
	call toms(x1,ac,1)
	call toms(y1,dc,0)
	write(*,2)x,y,ac,dc,epoch
2	format(2f12.3,'---> ',2a,' (',f6.1,')')
	endif
	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
