c combine 2048*2048 pic.
 	parameter (nfile=45)	
	parameter (nb=5259)
	real white(nfile),scale(nfile),d(nfile),s(nfile)
	integer*2 sx(nfile),sy(nfile)
	integer*2 bx(nb),by(nb)

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

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

	real b(2048*32,nfile)
	logical logi

	if(iargc().lt.2)then
	  write(*,*)
	  write(*,*)'         ******** Combine 2048*2048 Pic. ****2000.7'
	  write(*,*)'                      no bad_pixel delete'
	  write(*,*)'    cob batch_file out_file'
	  write(*,*)
	  write(*,*)'    if cob.par is not existed in default directory'
          write(*,*)'    default value will be taken as following'
	  write(*,*)'       3.3'
	  write(*,*)'       6.0'
	  write(*,*)'       6.0'
	  write(*,*)'       500000.'
	  write(*,*)'    if default value is not suitable, you can edit'
          write(*,*)'    a file named cob.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='/vega2/rhbin/cob.par',exist=logi)
        if(.not.logi)then
          open(51,file='cob.par',status='old')
          read(51,*)gain
	  read(51,*)s61
	  read(51,*)s62
	  read(51,*)s50
	  close (51)
	  write(*,*)'parameters come from file cob.par'
	else
	  gain=3.3
	  s61=6.
	  s62=6.
          s50=500000.
	  write(*,*)'default parameters:'
        endif
	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)
	  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)
	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)
	open(52,file='/vega2/rhbin/badpixel.bin',status='old',
     *       access='direct',recl=nb*4+4)
	read(52,rec=1)i,bx,by
	call swap2(bx,nb*2)
	call swap2(by,nb*2)
	close (52)
******* get each shift end ******
****** get mother pic star begin *****	
	call readata(a,f1)
        call shrink4(a,a)
	nstar=400
	call getstar(a,a(2048*1024+1),
     *      starx,stary,starv,nstar,white(1),sx,sy,m)
	write(*,*)'nstar: ',nstar
	k1=nstar/3.
	k2=k1+99
	if(k2.gt.nstar)k2=nstar
	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)	
	write(*,*)f1(1:22),'platen file'
	call star100(a,sx(1),sy(1),white(1),px,py,pz0,np)
	call badfile(a,nb,bx,by,0,0,1)
	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)	
	call star100(a,sx(i),sy(i),white(i),px,py,pz1,np)
	call badfile(a,nb,bx,by,sx(i),sy(i),i)
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
******** find scale end *******
	do 215 i=1,m
215	open(50+i,file=char(i/10+64)//char(mod(i,10)+64)//'.tmp',
     *  status='old',access='direct',recl=2048*32*4)
	do 216 i=1,m
	d(i)=sqrt(white(i)/gain)
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(50+i,l,b(1,i))	
	do 240 i=1,2048*32
	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(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)
	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)
	call badfile2(a,nb,bx,by) 
	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='rm '//char(i/10+64)//char(mod(i,10)+64)//'.tmp'
200	call system(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)
	real a(2048,2048)
	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=2048-ix2
        iy1=-iy1
        iy2=2048-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.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 badfile2(a,nb,x,y)
	real a(2048,2048)
	integer*2 x(1),y(1)
	bad=-999.
        do 10 k=1,nb
        i=x(k)
        j=y(k)
c	if(a(i,j).ne.bad)goto 10
	if(a(i,j).gt.0.)goto 10
        do 40 i1=i-1,1,-1
        if(a(i1,j).ne.bad)goto 41
40      continue
41      do 50 i2=i+1,2048
c       if(a(i2,j).ne.bad)goto 30
        if(a(i2,j).gt.0.)goto 30
50      continue
30      a(i,j)=(a(i1,j)+a(i2,j))*.5
10	continue
	end

	function allc(d,white,scale,m,sigma,s,s61,s50)
	real d(1),white(1),scale(1),s(1)
	real d1(45)			! normal value
	real d2(45)			! bad working unit
	bad=-999.			! s(i)=6*sigma(i)/scale(i)
c bad pixel
	allc=bad
	n=0
	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 do white hole

200	  do 202 i=1,m
	  if(d(i).lt.white(i)*0.8) d1(i)=x
202	  continue
	
c combine
201	x=0.
	do 210 i=1,m
        y= d1(i)*scale(i)+white(i)
        if(y.ge.29000.0)then
          x=s50
          goto 220
        endif
210     x=x+y
220     allc=x

	end	

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

	subroutine badfile(a,nb,bx,by,ix,iy,k)
	real a(2048,2048)
	integer*2 bx(1),by(1),ix,iy
c	do 10 i=1,nb
c10	a(bx(i),by(i))=-999.
c for no bad_pixel picture.
	i=iy*2048+ix
	call shift(a,i)
	open(52,file=char(k/10+64)//char(mod(k,10)+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

        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)
	integer*2 ix,iy,px(1),py(1)
	real a(2048,2048),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)
        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 shrink4(a,b)
        real a(1)
        integer*2 b(1)
        do 10 i=1,2048*2048
        x=a(i)
        if(x.gt.32767.)x=32767.
        if(x.lt.0.)x=0.
10      b(i)=x
	j=2048*2048
	do 20 i=1,2048*2048
20	b(j+i)=b(i)
        end

	subroutine getstar(a,b,starx,stary,starv,nstar,white,sx,sy,m)
	integer*2 a(2048,2048),b(1)
	real starx(1),stary(1),starv(1)
	integer*2 sx(1),sy(1)
	ix1=1
	ix2=-1
	iy1=1
	iy2=-1
	do 5 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)
5	if(sy(i).gt.iy2)iy2=sy(i)	
        ix1=-ix1
        ix2=2048-ix2
        iy1=-iy1
        iy2=2048-iy2
	ix1=ix1+20
	ix2=ix2-20
	iy1=iy1+20
	iy2=iy2-20

	kt=2000
	kk=0
	ktw=kt+white
	do 10 i=ix1,ix2
	do 10 j=iy1,iy2
	if(a(i,j).gt.ktw)then
	  call center9(i,j,a,ktw,ix,iy,iv,n9)
	  if(iv.lt.ktw*2)goto 10
	  x=ix
	  y=iy
	  call star(b,2048,2048,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)goto 11
	  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
	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 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

        subroutine star(map,n1,n2,xx,yy,ir,ierr,sky)
c when ir=4, 15*15 for center
        integer*2 map(n1,n2),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.n2)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.n2)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
 

