c combine 4096*4096 pic.
 	parameter (nfile=64)	
	integer*2 sx(nfile),sy(nfile)
	real      tx2(nfile),ty2(nfile),turn(nfile)
	integer*2 xx1(nfile),xx2(nfile),yy1(nfile),yy2(nfile)

	character*80 f0,f1,f2,f3,head(72)
	real*8 a8(8),b8(8),adcoef(2,3)
	
	real a(4096*4096)
	real b(4096*64,nfile)  ! need nother 4096*4096 array

	if(iargc().lt.2)then
	  write(*,*)
	  write(*,*)'         ******** Combine 4096*4096 Pic. ***2007,3'
	  write(*,*)
	  write(*,*)'  1. shift   2. turn   3. add '
          write(*,*)
          write(*,*)'    cobsta4k batch_file out_file [1]'
          write(*,*)'  if combined_sky>100000, 3rd parameter set 2'
          write(*,*)'  else pip25 does not work well'
          write(*,*)
	  stop
	endif

	call getarg(1,f0)
	open(51,file=f0,status='old')
	i=0
20	i=i+1
	read(51,'(a)',end=10)f2
	goto 20
10	if(i.lt.3)stop 'too less files, must>=2'
	close(51)
	m=i-1
	if(m.gt.nfile)m=nfile

	xf3=1.
	if(iargc().gt.2)then
	  call getarg(3,f3)
	  read(f3(1:),*)xf3
	endif
	write(*,*)'div: ',xf3
  

	open(51,file=f0,status='old')
	write(*,2)
2	format(35x,'shift(x,y)           UT.     exp.    sky'/78('-'))
	nt=0
	np=0			!  n frame
	do 30 i=1,m
	xx1(i)=1
	xx2(i)=4096
	yy1(i)=1
	yy2(i)=4096
	if(i.eq.1)then
	  read(51,'(a)')f1              	! mother file
	  call gethead(f1,head)
	  call geta8(head,a8)
	  call xytoad(a8,adcoef)
	  x=a8(7)
	  y=a8(8)
	  alpha=a8(7)
	  delta=a8(8)
	  f2=f1
	else
	  read(51,'(a)')f2
	  call gethead(f2,head)
	  call geta8(head,b8)
	  alpha=b8(7)
	  delta=b8(8)
	endif

        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
	tx2(i)=x2
	ty2(i)=y2
	turn(i)=(a8(2)-a8(3)-b8(2)+b8(3))/(a8(1)+a8(4))
	sx(i)=nint(x2)
	sy(i)=nint(y2)

        ipos=indexpos(head,'NAXIS   ')
        read(head(ipos)(65:69),'(i5)')n12
        if(n12.ne.0)then
        ipos=indexpos(head,'NAXIS1  ')
        read(head(ipos)(65:74),'(2i5)')xx1(i),xx2(i)
        ipos=indexpos(head,'NAXIS2  ')
        read(head(ipos)(65:74),'(2i5)')yy1(i),yy2(i)
	else
	n12=1
	endif
	np=np+n12

	ipos1=indexpos(head,'TIME    ')
	ipos2=indexpos(head,'EXPOSURE')
	read(head(ipos2)(26:),*)j
	nt=nt+j
	ipos3=indexpos(head,'VOLT2   ')
	if(head(ipos3)(42:46).ne."PIXEL")then
          write(*,*)"do statsky first !"
          stop
        endif
	read(head(ipos3)(22:),*)white
	write(*,1)i,f2(1:22),x2,y2,sx(i),sy(i),
     *  head(ipos1)(16:27),head(ipos2)(24:30),white
1	format(i2,':  ',a,2f6.1,2x,2i4,2a,f8.2)
30	continue
	close(51)
******* get each shift end ******
	open(51,file=f0,status='old')
	do 120 i=1,m
	read(51,'(a)')f2
	call readata(a,f2)	
        if(i.ne.1)call t_s(a,b,turn(i),tx2(i),ty2(i),4096)
        open(52,file=char(i+64)//'.tmp',status='unknown',
     *       access='direct',recl=4096*4096*4)
	if(i.eq.1)write(52,rec=1)a
        if(i.ne.1)write(52,rec=1)b
        close(52)
120	write(*,*)f2(1:22)
	close(51)
**********************************
	do 215 i=1,m
215	open(50+i,file=char(i+64)//'.tmp',status='old',
     *  access='direct',recl=4096*64*4)
	k=0
	do 220 l=1,64
        call dispdot()
	do 230 i=1,m
230	call readrec(50+i,l,b(1,i))	
	do 240 i=1,4096*4096/64
	xad=0.
	do 235 j=1,m
235	xad=xad+b(i,j)
	k=k+1
        if(xad.ge.5000000.)then
	  xad=5000000.
        else
	  xad=xad/xf3
	endif
240	a(k)=xad
220	continue
	do 250 i=1,m
250	close(50+i)
	write(*,*)
	open(51,file=f1,status='old',access='direct',recl=5760)
	read(51,rec=1)head
	close(51)
	ipos=indexpos(head,'END     ')
	if(ipos.lt.72)then
	  f2=head(ipos)(1:80)
	  head(ipos)(1:80)=head(72)(1:80)
	  head(72)(1:80)=f2
	endif	
	if(ipos.eq.73)head(72)(1:8)='END     '
	write(head(ipos2)(26:30),'(i5)')nt
	head(ipos2)(54:63)='(combined)'
	ipos=indexpos(head,'TTIME   ')
        n1=4096
        n0=n1/2
	if(ipos.ne.73)write(head(ipos)(25:30),"('>',i5)")nt
          x1=stat1(a,n1,n0-600,n0-600)
          x2=stat1(a,n1,n0-600,n0+600)
          x3=stat1(a,n1,n0,n0)
          x4=stat1(a,n1,n0+600,n0-600)
          x5=stat1(a,n1,n0+600,n0+600)                                          
          sky=amin1(x1,x2,x3,x4,x5)
          ipos=indexpos(head,'VOLT2   ')
          write(head(ipos)(22:),"(f9.1,' / SKY ADU/PIXEL')")sky
c delete 4 edge
	call edge(a,sx,sy,m,ix1,ix2,iy1,iy2,xx1,xx2,yy1,yy2)
	write(*,*)' Edge_frame: ccd(i,j)=ccd(i,j)+5000.'
	write(*,*)' line_number: col~~',ix1,ix2,'   raw~~',iy1,iy2

c add frame numbers, ix1,ix2,  iy1,iy2  96,12
	ipos=indexpos(head,'NAXIS   ')
	if(ipos.ne.73)write(head(ipos)(60:69),"('frame',i5)")np
	ipos=indexpos(head,'NAXIS1  ')
	if(ipos.ne.73)write(head(ipos)(60:74),"('range',2i5)")ix1,ix2
	ipos=indexpos(head,'NAXIS2  ')
	if(ipos.ne.73)write(head(ipos)(60:74),"('range',2i5)")iy1,iy2
	
	call getarg(2,f2)
	write(*,*)' combined file: '//f2(1:25)
	k=4096*4096*4+5760
	call f_blank(a,ix1,ix2,iy1,iy2)
	call swap4(a,4096*4096*4)
	open(50,file=f2,status='unknown',access='direct',recl=k)
	write(50,rec=1)head,a
	close(50)
	do 200 i=1,m
	  f0=char(i+64)//'.tmp'
200	call unlink(f0)
	write(*,*)'re coordinating!'
        f0='coord '//f2(1:25)//' !'
	call system(f0)
	stop
	end

        subroutine f_blank(a,ix1,ix2,iy1,iy2)
        real a(4096,4096)
        if(ix1.gt.1)then
          do 10 j=1,4096
          do 10 i=1,ix1
10        a(i,j)=0.
        endif
        if(ix2.lt.4096)then
          do 20 j=1,4096
          do 20 i=ix2,4096
20        a(i,j)=0.
        endif
        if(iy1.gt.1)then
          do 30 i=1,4096
          do 30 j=1,iy1
30        a(i,j)=0.
        endif
        if(iy2.lt.4096)then
          do 40 i=1,4096
          do 40 j=iy2,4096
40        a(i,j)=0.
        endif
        end


	subroutine edge(a,sx,sy,m,ix1,ix2,iy1,iy2,xx1,xx2,yy1,yy2)
	real a(4096,4096)
	integer*2 sx(1),sy(1)
	integer*2 xx1(1),xx2(1),yy1(1),yy2(1)
	do 5 i=1,m
	xx1(i)=xx1(1)-sx(i)
	xx2(i)=xx2(i)-sx(i)
	yy1(i)=yy1(i)-sy(i)
5	yy2(i)=yy2(i)-sy(i)

	ix1=1
	ix2=4096
	iy1=1
	iy2=4096
	do 10 i=1,m
	if(xx1(i).gt.ix1)ix1=xx1(i)
	if(xx2(i).lt.ix2)ix2=xx2(i)
	if(yy1(i).gt.iy1)iy1=yy1(i)
10	if(yy2(i).lt.iy2)iy2=yy2(i)

	k=iy1
	if(k.lt.1)k=1
        do 20 j=k,iy2
        if(ix1.ne.0)a(ix1,j)=a(ix1,j)+5000.
20      if(ix2.ne.4096)a(ix2,j)=a(ix2,j)+5000.
	k=ix1
	if(k.lt.1)k=1
        do 30 i=k,ix2
        if(iy1.ne.0)a(i,iy1)=a(i,iy1)+5000.
30      if(iy2.ne.4096)a(i,iy2)=a(i,iy2)+5000.
        if(ix1.eq.0)ix1=1
	if(iy1.eq.0)iy1=1
	end

        function stat1(a,n1,m1,m2)
        real a(n1,n1)
        integer*2 b(512,512)
        y=1.
        if(a(m1,m2).gt.20000.)y=0.1
        j1=m2-256
        do 10 j=1,512
        i1=m1-256
        j1=j1+1
        do 10 i=1,512
        i1=i1+1
10      b(i,j)=a(i1,j1)*y
        call whitexblack(b,512,512,x,s)
        stat1=x/y
        end
     
	subroutine readrec(k,l,b)
	real b(4096*64)
	read(k,rec=l)b
	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 readata(a,f1)
        integer c(720),a(1)
        character*80 f1,ch*2880
        equivalence (c,ch)
        nn=720         
        open(1,file=f1,status='old',access='direct',recl=nn*4)
        read(1,rec=1)ch
        n1=4096   
        n2=4096
        n=1   
55      k=index(ch,'END          ')
        if(k.eq.0)then
          n=n+1    
          read(1,rec=n)ch
          goto 55  
        endif
        i=0
        k1=n+1
        k2=n+(n1*n2-1)/nn
        do 60 k=k1,k2
        if(mod(k,200).eq.0)call dispdot()
        read(1,rec=k)ch
        do 60 j=1,nn
        i=i+1
60      a(i)=c(j)
        close (1)
	call swap4(a,n1*n2*4)
        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

	subroutine xy_rad(x,y,a8,alpha,delta)
c 4096*4096 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
	alpha=alpha*12./pi
	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

