character*80 f0,f1,f2,head(72)
	real*8 a8(8),b8(8),adcoef(2,3)
        real a(4096,4096),b(4096*4096)

	if(iargc().lt.3)then
	  write(*,*)
	  write(*,*)'         ******** shift_turn 4k Pic. ********'
	  write(*,*)
	  write(*,*)'    sht4k reference_file file out_file'
	  write(*,*)
          write(*,*)'                                     2007.11 '
	  stop
	endif

	call getarg(1,f1)         ! refenence
	call getarg(2,f2)         ! orient
	call getarg(3,f0)         ! out_put file

	write(*,2)
2	format(24x,'shift(x,y)        turn     exp.'/57('-'))

	call gethead(f1,head)
	call geta8(head,a8)
	call xytoad(a8,adcoef)
	x=a8(7)
	y=a8(8)

	call gethead(f2,head)
	call geta8(head,b8)
	ipos=indexpos(head,'A87     ')
	alpha=b8(7)
	delta=b8(8)
        call standc(x,y,alpha,delta,xi,xn)
        x2=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)-2048.5
        y2=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)-2048.5
	x2=-x2
        y2=-y2	
        ipos2=indexpos(head,'EXPOSURE')
	turn=(a8(2)-a8(3)-b8(2)+b8(3))/(a8(1)+a8(4))
	write(*,1)f2(1:22),x2,y2,turn,head(ipos2)(24:30)
1	format(a,2f6.1,2x,f11.6,a)

	call readfits(f2,a,n1,n2)
        call t_s(a,b,turn,x2,y2,n1)
	call swap4(b,n1*n2*4)
	ipos=indexpos(head,'A81     ')
        do 11 i=1,6
11      write(head(ipos-1+i)(11:30),'(E20.12)')a8(i)
        write(head(ipos+6)(11:30),'(f20.7)')a8(7)
        write(head(ipos+7)(11:30),'(f20.7)')a8(8)

        i=lnblnk(f0)
	f0(i+1:i+1)=char(0)
	call writefits(f0,head,80*72,1)
	call writefits(f0,b,n1*n2*4,0)
c do shift end
        write(*,*)'re coord !'
        call system("coord9 "//f0(1:i)//" !")
	end

        subroutine t_s(a,b,turn,z2,y2,n)
        real a(n,n),b(n,n)
        cosa=cos(turn)
        sina=sin(turn)
        cx=n/2+.5-z2
        cy=n/2+.5-y2
	do 10 j=1,n
         yy=j-cy
	do 10 i=1,n
        b(i,j)=0.
         xx=i-cx
         y= xx*sina + yy*cosa + cy + y2
	if(y.lt.1.or.y.gt.n)goto 10
         x= xx*cosa - yy*sina + cx + z2
	if(x.lt.1.or.x.gt.n)goto 10
	iy0=y
	iy1=iy0+1
	wy=y-iy0
        zy=1.-wy
	ix0=x
	ix1=ix0+1
	wx=x-ix0
        zx=1.-wx
        x1=zx*zy
        x2=wx*zy
        x3=zx*wy
        x4=wx*wy
	b(i,j)=a(ix0,iy0)*x1+a(ix1,iy0)*x2+a(ix0,iy1)*x3+a(ix1,iy1)*x4
10	continue
        end

	subroutine gethead(f1,head)
	character*80 f1,head(72)
	open(52,file=f1,status='old',access='direct',recl=5760)
	read(52,rec=1)head
	close (52)
	end
	
	subroutine geta8(head,a8)
	real*8 a8(8)
	character*80 head(72)
	ipos=indexpos(head,'A81     ')-1
	do 10 i=1,8
10	read(head(i+ipos)(12:),*)a8(i)
	end

        function indexpos(head,f1)
        character*80 head(72),f1*8
        do 10 indexpos=1,72
10      if(head(indexpos)(1:8).eq.f1)return
        end