c combine 2048*2048 pic.
 	parameter (nfile=45)	
	parameter (nb=5259)
	integer*2 sx(nfile),sy(nfile)
	integer*2 xx1(nfile),xx2(nfile),yy1(nfile),yy2(nfile)

	character*80 f0,f1,f2,head(72)
	real*8 a8(8),adcoef(2,3)
	
	real a(2048*2048)
	real b(2048*32,nfile)

	if(iargc().lt.2)then
	  write(*,*)
	  write(*,*)'         ******** Combine 2048*2048 Pic. ***2000.7'
	  write(*,*)
          write(*,*)
          write(*,*)'    sh-- mean fragment operate'
          write(*,*)'    cob      origion'
          write(*,*)'    cob50    if > 29000, put 500000'
          write(*,*)'    cob51    if >29000, process as bad point'
          write(*,*)'    cobad    add frames direct'
          write(*,*)
          write(*,*)'    cob batch_file out_file'
          write(*,*)
          write(*,*)'                1996,12,19'
	  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

	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)=2048
	yy1(i)=1
	yy2(i)=2048
	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)
	  f2=f1
	else
	  read(51,'(a)')f2
	  call gethead(f2,head)
	endif
  
	ipos=indexpos(head,'A87     ')
	read(head(ipos)(20:),*)alpha
	read(head(ipos+1)(20:),*)delta
        call standc(x,y,alpha,delta,xi,xn)
        x2=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)-1024.5
        y2=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)-1024.5
	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)	
	call shfile(a,sx(i),sy(i),i)
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=2048*32*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,2048*32
	xad=0.
	do 235 j=1,m
235	xad=xad+b(i,j)
	k=k+1
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   ')
	if(ipos.ne.73)write(head(ipos)(25:30),"('>',i5)")nt
          x1=statt(a,b,256,256)
          x2=statt(a,b,256,1280)
          x3=statt(a,b,768,768)
          x4=statt(a,b,1280,256)
          x5=statt(a,b,1280,1280)
          sky=amin1(x1,x2,x3,x4,x5)
          ipos=indexpos(head,'VOLT2   ')
          write(head(ipos)(22:),"(f9.3,' / 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=2048*2048*4+5760
	call f_blank(a,ix1,ix2,iy1,iy2)
	call swap4(a,2048*2048*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(2048,2048)
        if(ix1.gt.1)then
          do 10 j=1,2048
          do 10 i=1,ix1
10        a(i,j)=0.
        endif
        if(ix2.lt.2048)then
          do 20 j=1,2048
          do 20 i=ix2,2048
20        a(i,j)=0.
        endif
        if(iy1.gt.1)then
          do 30 i=1,2048
          do 30 j=1,iy1
30        a(i,j)=0.
        endif
        if(iy2.lt.2048)then
          do 40 i=1,2048
          do 40 j=iy2,2048
40        a(i,j)=0.
        endif
        end


	subroutine edge(a,sx,sy,m,ix1,ix2,iy1,iy2,xx1,xx2,yy1,yy2)
	real a(2048,2048)
	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=2048
	iy1=1
	iy2=2048
	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.2048)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.2048)a(i,iy2)=a(i,iy2)+5000.
        if(ix1.eq.0)ix1=1
	if(iy1.eq.0)iy1=1
	end

        function statt(a,b,m1,m2)
	real a(2048,2048)
        integer*2 b(512,512)
	y=amin1(a(m1,m2),a(m1+512,m2+512),a(m1,m2+512),
     *          a(m1+512,m2),a(m1+256,m2+256))
	y=amin1(y,a(m1+256,m2),a(m1+512,m2+256),
     *          a(m1+256,m2+512),a(m1,m2+256))
	y=y*0.8
        j1=m2
        do 10 j=1,512
        i1=m1
        j1=j1+1
        do 10 i=1,512
        i1=i1+1
        x=a(i1,j1)-y
	if(x.lt.-32767.)x=-32767.
	if(x.gt. 32767.)x= 32767.
10	b(i,j)=x
        call whitexblack(b,512,512,x,s)
        stat=x+y
        end

	subroutine readrec(k,l,b)
	real b(2048*32)
	read(k,rec=l)b
	end

        subroutine shfile(a,ix,iy,k)
        real a(2048,2048)
        integer*2 ix,iy
        i=iy*2048+ix
        call shift(a,i)
        open(52,file=char(k+64)//'.tmp',status='unknown',
     *       access='direct',recl=2048*2048*4)
        write(52,rec=1)a
        close(52)
        end
 
	subroutine shift(a,k)
	real a(2048*2048)
	mm=2048*2048
	if(k.gt.0)then
	  j=k
	  do 10 i=1,mm-k
	  j=j+1
10	  a(i)=a(j)
	endif
	if(k.lt.0)then
	  j=mm+k+1
	  do 20 i=mm,1-k,-1
	  j=j-1
20	  a(i)=a(j)
	endif
	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=2048   
        n2=2048
        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 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
	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

