c 4k*4k ---> 2k*2k

	character*80 f1,f2,f3,head1(72),head2(72)
	real*8 a8(8),b8(8),adcoef(2,3)
	real a(2048,2048),b(4096,4096)
	real c(4000,4000)
	real*8 xi,xn,alp,del

	if(iargc().lt.3)then
	  write(*,*)
	  write(*,*)'         ******** 4k*4k ccd -->2k*2k ccd ********'
	  write(*,*)
          write(*,*)
          write(*,*)'      4kto2k refence(2k) input(4k) output(2k)'
          write(*,*)
          write(*,*)'            jiang zhaoji    2006,9,13'
          write(*,*)
          write(*,*)
	  stop
	endif

	call getarg(1,f1)
	open(1,file=f1,status='old',iostat=i,access='direct',recl=5760)
	if(i.ne.0)stop 'refence file not found!'
	read(1,rec=1)head1
	ipos=indexpos(head1,'NAXIS1   ')
	read(head1(ipos)(16:),*)n1
	if(n1.ne.2048)stop 'ref_file not a 2k*2k file'
	close(1)
	call geta8(head1,a8)

	call getarg(2,f2)
	open(1,file=f2,status='old',iostat=i,access='direct',recl=5760)
	if(i.ne.0)stop 'input file not found!'
	read(1,rec=1)head2
	ipos=indexpos(head2,'NAXIS1  ')
	read(head2(ipos)(16:),*)n1
	if(n1.ne.4096.and.n1.ne.4000)stop 'inp_file not a 4k*4k file'
	close(1)
	write(*,*)'size: ',n1
	if(n1.eq.4096)call readfits(f2,b,n1,n1)
	if(n1.eq.4000)call readfits(f2,c,n1,n1)
	call geta8(head2,b8)
	call xy2ad(b8,adcoef)

	y=a8(7)-b8(7)
	if(y.lt.0)y=-y
	if(y.gt.6.)y=y-2*3.141592654
	x=y**2+(a8(8)-b8(8))**2
	if(x.gt.0.0001)stop 'wrong coord frames !'

	do 10 i=1,2048
	if(i/100*100.eq.i)call dispdot()
	do 10 j=1,2048
	call xy_rad(float(i),float(j),a8,alpha,delta)
	alp=alpha
	del=delta
	call d_standc(b8(7),b8(8),alp,del,xi,xn)
        x=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
        y=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)
        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
        if(n1.eq.4096)a(2049-i,2049-j)=b(ix0,iy0)*x1+b(ix1,iy0)*x2+
     c                   b(ix0,iy1)*x3+b(ix1,iy1)*x4
        if(n1.eq.4000)a(2049-i,2049-j)=c(ix0,iy0)*x1+c(ix1,iy0)*x2+
     c                   c(ix0,iy1)*x3+c(ix1,iy1)*x4
10	continue	

	call getarg(3,f3)
	n1=2048
        call swap4(a,n1*n1*4)
	i=indexpos(head1,'IMNAME   ')
	j=indexpos(head2,'IMNAME   ')
	do 20 k=0,10
20	head1(k+i)(1:80)=head2(k+j)(1:80)
	i=indexpos(head1,'SHUTTER  ')
	j=indexpos(head2,'GAIN     ')
        if(i.ne.73)then
	do 21 k=0,2
21	head1(k+i)(1:80)=head2(k+j)(1:80)
	endif
	i=indexpos(head1,'TTIME   ')
	j=indexpos(head2,'TTIME   ')
	head1(i)(1:80)=head2(j)(1:80)
        i=1
	j=2048
	k=indexpos(head1,'NAXIS1  ')
        write(head1(k)(60:74),"('range',2i5)")i,j
        write(head1(k+1)(60:74),"('range',2i5)")i,j
        k=n1*n1*4+5760
        open(1,file=f3,status='unknown',access='direct',recl=k)
	write(1,rec=1)head1,a
	close(1)

	write(*,*)'re coordinating!'
        f1='coord '//f3(1:25)//' !'
        call system(f1)

	end
 
	subroutine exp4096(a,b)
	real a(4000*4000),b(4096*4096)
	k=4096*4096
	j=4000*4000
	do 20 i1=1,4000
	k=k-96
	do 10 i2=1,4000
	b(k)=a(j)
	k=k-1
	j=j-1
10	continue
20	continue
	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

	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)                            ! from fits_head "A81"--"A88"
        pi=4.*atan(1.)
        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=cos(dc)
        td=tan(dc)
        alpha=atan(xi/cd/(1.-xn*td))
        delta=atan((xn+td)*cos(alpha)/(1.-xn*td))
        alpha=alpha+rc
        if(abs(dc-delta).gt.2.)then
          delta=-delta
          alpha=alpha+pi
        endif
c	alpha=alpha*12./pi
c	delta=delta*180./pi
        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

        subroutine xy2ad(x,a)
c 2 set 6 plate coef convert
        real*8 x(2,3),a(2,3),z
        z=x(1,1)*x(2,2)-x(1,2)*x(2,1)
        a(1,1)=x(2,2)/z
        a(1,2)=-x(1,2)/z
        a(1,3)=(x(1,2)*x(2,3)-x(1,3)*x(2,2))/z
        a(2,1)=-x(2,1)/z
        a(2,2)=x(1,1)/z
        a(2,3)=(x(1,3)*x(2,1)-x(1,1)*x(2,3))/z
        end

        subroutine d_standc(rc,dc,ra,de,xi,xn)
c alpha delta --------> X Y
c rc dc : center of plate
c rc dc ra de in arc
        real*8 rc,dc,ra,de,xi,xn,tmp
        sd=dsin(dc)
        cd=dcos(dc)
        td=dtan(de)
        co=dcos(ra-rc)
        tmp=sd*td+cd*co
        xi=dsin(ra-rc)/tmp
        xn=(cd*td-sd*co)/tmp
        end

