c combine 4k*4k pic.
 	parameter (nfile=89)	
	real white(nfile),scale(nfile),d(nfile),s(nfile)
	integer*2 sx(nfile),sy(nfile)
	integer*2 bx(30000),by(30000)

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

	real starx(400),stary(400),starv(400)
	real p0(100)
	integer*2 px(100),py(100),ac
	real pz0(100),pz1(100)

	real b(4096*64,nfile)
	logical logi

        ac=iargc()
	if(ac.lt.2)then
	  write(*,*)
	  write(*,*)'         ******** Combine 4k*4k Pic. ***2006.3'
	  write(*,*)
	  write(*,*)'    cob batch_file out_file'
	  write(*,*)
	  write(*,*)'    if cob4k.par is not existed in default directory'
          write(*,*)'    default value will be taken as following'
	  write(*,*)'       2.0    ! gain' 
	  write(*,*)'       6.0    ! sigam for star'
	  write(*,*)'       6.0    ! sigma for background'
	  write(*,*)'       500000 ! if saturating'
	  write(*,*)'    if default value is not suitable, you can edit'
          write(*,*)'    a file named cob4k.par inclued 4 line in your'
          write(*,*)'    default directory, '
	  write(*,*)'                                      1997,11'
	  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.4)stop 'too less files, must>=3'
	close(51)
	m=i-1
	if(m.gt.nfile)m=nfile
        inquire(file='cob4k.par',exist=logi)
        if(logi)then
          open(51,file='cob4k.par',status='old')
          read(51,*)gain
	  read(51,*)s61
	  read(51,*)s62
	  read(51,*)s50
	  close (51)
	  write(*,*)'parameters come from your home_directory'
	  goto 1234
        endif
        inquire(file='/vega2/rhbin/cob4k.par',exist=logi)
        if(logi)then
          open(51,file='/vega2/rhbin/cob4k.par',status='old')
          read(51,*)gain
	  read(51,*)s61
	  read(51,*)s62
	  read(51,*)s50
	  close (51)
	  write(*,*)'parameters come from /vega2/rhbin'
	  goto 1234
        endif
	gain=2.0
	s61=6.
	s62=6.
        s50=5000000.
	write(*,*)'default parameters:'

1234	n49=s50-1000.+0.1
	write(*,*)'gain=',gain,'    n1(star)=',s61,'    n2(background)=',s62 
        write(*,*)'saturated data=',s50,'     high good data=',n49
	open(51,file=f0,status='old')
	write(*,2)
2	format(35x,'shift(x,y)           UT.     exp.    sky'/78('-'))
	nt=0
	do 30 i=1,m
	if(i.eq.1)then
	  read(51,'(a)')f1              	! mother file
	  call gethead(f1,head)
	  read(head(4)(20:),*)n1
	  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)-2048.5
        y2=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)-2048.5
        x2=-x2                                              ! just +-
	y2=-y2
	sx(i)=nint(x2)
	sy(i)=nint(y2)
	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(i)
	write(*,1)i,f2(1:22),x2,y2,sx(i),sy(i),
     *  head(ipos1)(16:27),head(ipos2)(24:30),white(i)
1	format(i2,':  ',a,2f6.1,2x,2i4,2a,f8.2)
30	continue
	close(51)

        write(*,*)" no bad_pixel corrected"

******* get each shift end ******
****** get mother pic star begin *****	
12345	call readata(a,f1,n1)
	write(*,*)
        call shrink4(a,a,n1)
	nstar=400
	call getstar(a,a(n1*n1/2+1),starx,stary,starv,nstar,white(1),n1)
	write(*,*)'nstar: ',nstar
	k1=nstar/3.
	k2=k1+99
	if(k2.gt.nstar)k2=nstar
c	  write(*,*)k1,k2
	np=0
	do 110 i=k1,k2
	np=np+1
	px(np)=starx(i)+0.5
110	py(np)=stary(i)+0.5	
	call readata(a,f1,n1)	
	write(*,*)f1(1:22),'platen file'
	call star100(a,sx(1),sy(1),white(1),px,py,pz0,np,n1)
	call badfile(a,nb,bx,by,0,0,1,n1)
	scale(1)=1.
	open(51,file=f0,status='old')
	read(51,'(a)')f2
c	call pgbegin(0,'/su',2,2)
	do 120 i=2,m
	read(51,'(a)')f2
	call readata(a,f2,n1)	
	call star100(a,sx(i),sy(i),white(i),px,py,pz1,np,n1)
	call badfile(a,nb,bx,by,sx(i),sy(i),i,n1)
c	call pgenv(0.,120.,-10.,10.1,0,0)
	do 130 j=1,np
	p0(j)=j
130	pz1(j)=pz1(j)/pz0(j)
c	call pgline(np,p0,pz1)
	scale(i)=xmedian(pz1,np)
120	write(*,*)f2(1:22),'scale:',scale(i)
c	call pgend
	close(51)
******** find scale end *******
	do 215 i=1,m
215	open(10+i,file=char(i/10+64)//char(mod(i,10)+64)//'.tmp',
     *  status='old',access='direct',recl=n1*n1/16)
	do 216 i=1,m
c 20030808
	d(i)=sqrt(white(i))
216	s(i)=s62*d(i)/scale(i)
	sigma=xmedian(d,m)

	k=0
	do 220 l=1,64
        call dispdot()
	do 230 i=1,m
230	call readrec(10+i,l,b(1,i),n1)	
	do 240 i=1,n1*n1/64
	do 235 j=1,m
235	d(j)=b(i,j)
	k=k+1
240	a(k)=allc(d,white,scale,m,sigma,s,s61,s50)		
220	continue

	do 250 i=1,m
250	close(10+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   ')
        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,n1)
	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:77),"('frame',i5,i8)")m,n49
	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=n1*n1*4+5760
	call f_blank(a,ix1,ix2,iy1,iy2,n1)
	call swap4(a,n1*n1*4)
	open(50,file=f2,status='unknown',access='direct',recl=k)
	call wfile(head,a,n1)
	do 200 i=1,m
	  f0='rm '//char(i/10+64)//char(mod(i,10)+64)//'.tmp'
200	call system(f0)
	write(*,*)'re coordinating!'
        f0='coord9 '//f2(1:25)//' !'
	call system(f0)
	stop
	end

	subroutine wfile(head,a,n1)
	character*80 head(72)
	real a(n1*n1)
	write(50,rec=1)head,a
	close(50)
	end

	subroutine readrec(k,l,b,n1)
	real b(n1*n1/64)
	read(k,rec=l)b
	end

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

	subroutine edge(a,sx,sy,m,ix1,ix2,iy1,iy2,n1)
	real a(n1,n1)
	integer*2 sx(1),sy(1)
	ix1=1
	ix2=-1
	iy1=1
	iy2=-1
	do 10 i=1,m
	if(sx(i).lt.ix1)ix1=sx(i)
	if(sx(i).gt.ix2)ix2=sx(i)
	if(sy(i).lt.iy1)iy1=sy(i)
10	if(sy(i).gt.iy2)iy2=sy(i)	
        ix1=-ix1
        ix2=n1-ix2
        iy1=-iy1
        iy2=n1-iy2
	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.n1)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.n1)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

	function allc(d,white,scale,m,sigma,s,s61,s50)
	real d(1),white(1),scale(1),s(1)
	real d1(99)			! normal value
	real d2(99)			! bad working unit
	bad=-999.			! s(i)=6*sigma(i)/scale(i)
c bad pixel
	allc=bad
	n=0
        do 5 i=1,m
5	if(d(i).lt.white(i)*.5)d(i)=bad
	do 10 i=1,m
10	if(d(i).le.bad)n=n+1			! eq
	if(n.eq.m)return
	do 20 i=1,m
20	if(d(i).gt.bad)d1(i)=(d(i)-white(i))/scale(i)  !ne
	if(n.eq.0)goto 100
	n=0
	do 30 i=1,m
	if(d(i).le.bad)goto 30                    ! eq
	n=n+1
	d2(n)=d1(i)
30	continue
	x=xmedian(d2,n)
	do 40 i=1,m
40	if(d(i).le.bad)d1(i)=x                    ! eq
c cosmos ray
100	do 110 i=1,m
110	d2(i)=d1(i)
	x=xmedian(d2,m)
	if(x.gt.sigma)then		! star
	  if(d2(1).lt.s61*x)goto 200
	  do 120 i=1,m
120	  if(d1(i).eq.d2(1))d1(i)=x
	else				! background
	  do 130 i=1,m
130	  if( d1(i).eq.d2(1) .and. d1(i).gt.s(i) )d1(i)=x
	endif
c combine
200	x=0.
	do 210 i=1,m
        y= d1(i)*scale(i)+white(i)
        if(y.ge.60000.0)then
          x=s50
          goto 220
        endif
210     x=x+y
220     allc=x
	end	

	subroutine badfile(a,nb,bx,by,ix,iy,k,n1)
	real a(n1,n1)
	integer*2 bx(1),by(1),ix,iy
	i=iy*n1+ix
	call shift(a,i,n1)
	open(52,file=char(k/10+64)//char(mod(k,10)+64)//'.tmp',
     *  status='unknown',access='direct',recl=n1*n1*4)
	write(52,rec=1)a
	close(52)
	end	

	subroutine shift(a,k,n1)
	real a(n1*n1)
	mm=n1*n1
	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

        function xmedian(a,n)
        real a(n)
        do 10 i=1,n-1
        do 10 j=i+1,n
        if(a(i).gt.a(j))goto 10
        x=a(i)   
        a(i)=a(j)
        a(j)=x  
10      continue
        m=n/2
        xmedian=a(m+1)
        if(m+m.ne.n)return
	xmedian=(a(m)+a(m+1))/2
        end

	subroutine star100(a,ix,iy,white,px,py,pz,n,n1)
	integer*2 ix,iy,px(1),py(1)
	real a(n1,n1),pz(1)
	do 10 i=1,n
	z=0.
	do 20 ii=px(i)+ix-7,px(i)+ix+7
	do 20 jj=py(i)+iy-7,py(i)+iy+7
20	z=z+a(ii,jj)-white
10	pz(i)=z
	end

        subroutine readata(a,f1,n1)
        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
        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*n1-1)/nn
        do 60 k=k1,k2
        if(mod(k,800).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*n1*4)
        end
 
        subroutine shrink4(a,b,n1)
        real a(1)
        integer*2 b(1)
        do 10 i=1,n1*n1
        x=a(i)
        if(x.gt.32767.)x=32767.
        if(x.lt.0.)x=0.
10      b(i)=x
	j=n1*n1
	do 20 i=1,n1*n1
20	b(j+i)=b(i)
        end

	subroutine getstar(a,b,starx,stary,starv,nstar,white,n1)
	integer*2 a(n1,n1),b(1)
	real starx(1),stary(1),starv(1)
	integer*2 sx(1),sy(1)
	ix1=600
	ix2=n1-600
	iy1=600
	iy2=n1-600
	kt=3000
5	kk=0
	if(white.lt.0.)white=0.
	ktw=kt+white
c	  write(*,*)ktw
	call putdata2(b,a,n1)
	do 10 i=ix1,ix2
	do 10 j=iy1,iy2
	if(a(i,j).gt.ktw)then
	  call center4k(i,j,a,ktw,ix,iy,iv,n9,n1)
	  if(iv.lt.ktw*2)goto 10
	  x=ix
	  y=iy
	  call star(b,n1,x,y,4,ierr,white)
	  if(ierr.ne.0)x=-99.
	  z=(ix-x)**2+(iy-y)**2
	  if(z.gt.5.)goto 10
	  if(kk.eq.nstar)then
	    if(kt.lt.3000)goto 11
            kt=kt+2000
	    if(kt.lt.15000)goto 5
	    goto 11
          endif
	  kk=kk+1
	  starx(kk)=x
	  stary(kk)=y
	  starv(kk)=iv
	  if(mod(kk,10).eq.1)call dispdot()
	endif
10	continue
11	nstar=kk
	if(kk.lt.40.and.kt.eq.3000)then
	  kt=1000
          if(kk.lt.20)kt=600
	  if(kk.lt.10)kt=300
	  goto 5
	endif
	do 20 i=1,nstar-1
	do 20 j=i+1,nstar
	if(starv(i).ge.starv(j))goto 20
	x=starx(i)
	starx(i)=starx(j)
	starx(j)=x
	x=stary(i)
	stary(i)=stary(j)
	stary(j)=x
	x=starv(i)
	starv(i)=starv(j)
	starv(j)=x
20	continue
	if(nstar.gt.300)nstar=300
	end

	subroutine putdata2(a,b,n1)
	integer*2 a(1),b(1)
	do 10 i=1,n1*n1
10	b(i)=a(i)
	end

        subroutine center4k(i1,j1,a,kt,ix,iy,iv,n,n1)
        integer t1(200)
        integer*2 c1(200),a(n1,n1)
        iv=0
        n=0
        i=i1
        j=j1
10      k=0
        iu=j
        id=j
        if(a(i,j).lt.kt)goto 100        ! find end
        k=k+a(i,j)
        a(i,j)=-a(i,j)

20      if(a(i,iu-1).lt.kt)goto 30      ! up find end
        iu=iu-1
          if(iu.lt.1)goto 30
        k=k+a(i,iu)
        a(i,iu)=-a(i,iu)
        goto 20
30      if(a(i,id+1).lt.kt)goto 40      ! down find end
        id=id+1
          if(id.gt.n1)goto 40
        k=k+a(i,id)
        a(i,id)=-a(i,id)
        goto 30
40      iv=iv+k
        if(n.lt.200)n=n+1
        t1(n)=k
        c1(n)=iu+id
        i=i+1
        j=c1(n)*.5
        if(a(i,j).lt.a(i,j+1))j=j+1
        goto 10

100     iu=t1(1)
        id=1
        do 110 k=2,n
        if(t1(k).gt.iu)then
          iu=t1(k)
          id=k
        endif
110     continue
        ix=i1+id-1
        iu=0
        do 120 k=1,n
120     iu=iu+c1(k)
        iy=iu/2/n
        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 4k*4k 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

        subroutine star(map,n1,xx,yy,ir,ierr,sky)
c when ir=4, 15*15 for center
        integer*2 map(n1,n1),maps(-25:25,-25:25)
        real b(-25:25),xmap(-25:25,-25:25)
c sky first
        n25=25
        x=xx
        y=yy
        kx=nint(x)
        ky=nint(y)
        if(kx.le.n25)kx=n25+1
        if(kx.gt.n1-n25)kx=n1-n25
        if(ky.le.n25)ky=n25+1
        if(ky.gt.n2-n25)ky=n2-n25
        do 10 j=-n25,n25
        do 10 i=-n25,n25
10      maps(i,j)=map(kx+i,ky+j)
        n51=n25+n25+1
        ierr=1                                  ! error =1
        i5=ir*1.1+2
        do 1000 ll=1,7
         if(x.lt.1..or.x.gt.n1)return
         if(y.lt.1..or.y.gt.n1)return
        ix=nint(x)
        iy=nint(y)
cccccc x center    weight_center method
        do 30 i=ix-i5,ix+i5
        z=0.
        do 40 j=iy-i5,iy+i5
40      z=z+map(i,j)-sky          
30      b(i-ix)=z                       ! for weight_center cal.
        z=0.
        s=0.
        do 50 i=-i5,i5
        z=z+b(i)
50      s=s+b(i)*i
        s=s/z
        x=ix+s
cccccc y center
        do 130 j=iy-i5,iy+i5
        z=0.
        do 140 i=ix-i5,ix+i5
140     z=z+map(i,j)-sky
130     b(j-iy)=z
        z=0.
        s=0.
        do 150 i=-i5,i5
        z=z+b(i)
150     s=s+b(i)*i
        s=s/z
        y=iy+s
1000    continue
        ix=nint(x)
        iy=nint(y)
        do 101 j=-n25,n25
        do 101 i=-n25,n25
101     xmap(i,j)=map(ix+i,iy+j)
        x=ix-25+center_1(xmap,sky,51,1)
        y=iy-25+center_1(xmap,sky,51,0)
         if(x.lt.1..or.x.gt.n1)return
         if(y.lt.1..or.y.gt.n1)return
        ierr=0
          xx=x
          yy=y
        return
        end
 
        function center_1(a,sky,n,ip)
        real a(51,51)
        integer b(51)
        skyn=sky*n
        if(ip.eq.1)then
        do 10 i=1,n
        x=0.
        do 20 j=1,n
20      x=x+a(i,j)
10      b(i)=x-skyn
        else
        do 110 i=1,n
        x=0.
        do 120 j=1,n
120     x=x+a(j,i)
110     b(i)=x-skyn
        endif
        call histat(b,mode,maxh,xmean,xpeak,fsigma,gsigma,n)
        center_1=xpeak
        end