c combine 2048*2048 pic.
 	parameter (nfile=64)              ! cannot change it, "b" use it	
	real white(nfile),scale(nfile),d(nfile),s(nfile),g(nfile)
	integer*2 sx(nfile),sy(nfile),badtype(nfile)
        integer*2 sx3(nfile),sy3(nfile)
        real      sz(nfile)                        ! turn angle
	integer*2 bx(30000),by(30000)

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

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

	real b(2048*3*96,nfile)                   ! ==b(6144,6144)
	logical logi

        ac=iargc()
	if(ac.lt.2)then
	  write(*,*)
	  write(*,*)'      ******** Combine 2k Pic -> 6k pic ***2004.4'
	  write(*,*)'    cob batch_file out_file [!]'
	  write(*,*)'    if cob.par is not existed in default directory'
          write(*,*)'    default value will be taken as following'
	  write(*,*)'       3.4    ! gain, before 2003,2,5'
	  write(*,*)'       6.0    ! sigam for star'
	  write(*,*)'       6.0    ! sigma for background'
	  write(*,*)'       500000 ! default value, if saturating'
          write(*,*)'       8.26   ! gain, after 2003,2,5'
	  write(*,*)'    if default value is not suitable, you can edit'
          write(*,*)'    a file named cob.par inclued 5 line in your'
          write(*,*)'    default directory, '
	  write(*,*)'                                      1997,11'
          write(*,*)'    after 2003,2,5 fits, at end, data*=8.26/3.4'
          write(*,*)'                  2003,5; 2003,8 modify'
          write(*,*)'                  auto selected 3 bad_pixel tables'
          write(*,*)'                  zjzx.bad   ljl.bad   ccd5.bad'
          write(*,*)'    2k-->6k'
          write(*,*)'    notice: must left huge disk_space'
          write(*,*)'    ex: if combine 20 frames, disk_space must > 3G'
	  write(*,*)'                                      2004,4'
	  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='cob.par',exist=logi)
        if(logi)then
          open(51,file='cob.par',status='old')
          read(51,*)gain
	  read(51,*)s61
	  read(51,*)s62
	  read(51,*)s50
          read(51,*,end=501)gain2003
502	  close (51)
	  write(*,*)'parameters come from default directory'
	  goto 1234
501	  gain2003=8.0
          goto 502
        endif
        inquire(file='/vega2/rhbin/cob.par',exist=logi)
        if(logi)then
          open(51,file='/vega2/rhbin/cob.par',status='old')
          read(51,*)gain
	  read(51,*)s61
	  read(51,*)s62
	  read(51,*)s50
          read(51,*)gain2003
	  close (51)
	  write(*,*)'parameters come from /vega2/rhbin'
	  goto 1234
        endif
	gain=3.4
	s61=6.
	s62=6.
        s50=500000.
        gain2003=8.0
	write(*,*)'default parameters:'

1234	n49=s50-1000.+0.1
	write(*,*)'gain=',gain,'    n1(star)=',s61,'    n2(background)=',s62 
        write(*,*)'gain_2003=',gain2003
        write(*,*)'saturated data=',s50,'     high good data=',n49
	open(51,file=f0,status='old')
	write(*,2)
2	format(35x,'shift(x,y)          turn     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)
	sx3(i)=nint(x2*3.)
	sy3(i)=nint(y2*3.)
c	ipos1=indexpos(head,'TIME    ')
          call geta8(head,b8)
          sz(i)=(a8(2)-a8(3)-b8(2)+b8(3))/(a8(1)+a8(4))
	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*3,y2*3,sx3(i),sy3(i),sz(i),
     *  head(ipos2)(24:30),white(i)
1	format(i2,':  ',a,2f6.1,2x,2i4,f11.5,a,f8.2)

        ipos=indexpos(head,'INSTRUME')
        if(ipos.gt.70)stop 'not our fits file**'
        badtype(i)=0                             !                 zjzx.bad
        if(head(ipos)(13:13).eq.'1')badtype(i)=1 ! new controller, ljl.bad
        if(indexpos(head,'MPP     ').gt.70)badtype(i)=0
        if(head(ipos)(13:13).eq.'5')badtype(i)=5 ! after 2003
        if(head(ipos)(17:17).eq.'4')badtype(i)=4 ! thin ccd, no bad pixel
        if(ac.ge.3)badtype(i)=-1

c 20030808
	ipos=indexpos(head,'DATE-OBS')
        g(i)=gain
        if(ipos.lt.70)then
	  read(head(ipos)(22:29),'(i2,1x,i2,1x,i2)')id,im,iy
          if(iy.lt.90)then
            j=iy*10000+im*100+id
            if(j.gt.30205)g(i)=gain2003
          endif
        endif
30	continue
	close(51)

        if(ac.ge.3)then
          write(*,*)" no bad_pixel corrected"
          goto 12345
	endif
          

        inquire(file='zjzx.bad',exist=logi)
        if(logi)then
          write(*,*)'read zjzx.bad from default directory'
	  open(52,file='zjzx.bad',status='old')
        else
          write(*,*)'read zjzx.bad from /vega2/rhbin'
	  open(52,file='/vega2/rhbin/zjzx.bad',status='old')
	endif
        nb=1
301     read(52,*,end=302)bx(nb),by(nb)
        if(bx(nb).lt.1.or.bx(nb).gt.2048)goto 301
        if(by(nb).lt.1.or.by(nb).gt.2048)goto 301
        if(nb.lt.10000)nb=nb+1
	goto 301	
302	nb=nb-1
        write(*,*)"readin badpixel table from file zjzx.bad:",nb 
	close (52)

        inquire(file='ljl.bad',exist=logi)
        if(logi)then
          write(*,*)'read ljl.bad from default directory'
	  open(52,file='ljl.bad',status='old')
        else
          write(*,*)'read ljl.bad from /vega2/rhbin'
	  open(52,file='/vega2/rhbin/ljl.bad',status='old')
	endif
        lb=10001
3011     read(52,*,end=3021)bx(lb),by(lb)
        if(bx(lb).lt.1.or.bx(lb).gt.2048)goto 3011
        if(by(lb).lt.1.or.by(lb).gt.2048)goto 3011
        if(lb.lt.20000)lb=lb+1
	goto 3011	
3021	lb=lb-1
        write(*,*)"readin badpixel table from file  ljl.bad:",lb-10000
	close (52)

        inquire(file='ccd5.bad',exist=logi)
        if(logi)then
          write(*,*)'read ccd5.bad from default directory'
	  open(52,file='ccd5.bad',status='old')
        else
          write(*,*)'read ccd5.bad from /vega2/rhbin'
	  open(52,file='/vega2/rhbin/ccd5.bad',status='old')
	endif
        mb=20001
4011     read(52,*,end=4021)bx(mb),by(mb)
        if(bx(mb).lt.1.or.bx(mb).gt.2048)goto 4011
        if(by(mb).lt.1.or.by(mb).gt.2048)goto 4011
        if(lb.lt.30000)mb=mb+1
	goto 4011	
4021	mb=mb-1
        write(*,*)"readin badpixel table from file ccd5.bad:",mb-20000
	close (52)

******* get each shift end ******
****** get mother pic star begin *****	
12345	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,g(1),gain)
	write(*,*)'star: ',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,b,sx(1),sy(1),white(1),px,py,pz0,np,0.)
	call badfile(a,b,nb,lb,mb,bx,by,badtype(1),0,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,b,sx(i),sy(i),white(i),px,py,pz1,np,sz(i))
	call badfile(a,b,nb,lb,mb,bx,by,badtype(i),sx3(i),sy3(i),sz(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
	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=6144*96*4)
	do 216 i=1,m
c 20030808
	d(i)=sqrt(white(i)/g(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))	
	do 240 i=1,6144*96
	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,g,gain)		
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   ')
	if(ipos.ne.73)write(head(ipos)(24:30),"(' >',i5)")nt
        call shrink9(a,6144,b)
        x1=statt(b,256,256,m)
        x2=statt(b,256,1280,m)
        x3=statt(b,768,768,m)
        x4=statt(b,1280,256,m)
        x5=statt(b,1280,1280,m)
        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(sx3,sy3,m,ix1,ix2,iy1,iy2)
	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
	if(ipos.ne.73)write(head(ipos)(27:30),"('6144')")
	ipos=indexpos(head,'NAXIS2  ')
	if(ipos.ne.73)write(head(ipos)(60:74),"('range',2i5)")iy1,iy2
	if(ipos.ne.73)write(head(ipos)(27:30),"('6144')")
	
	call getarg(2,f2)
	write(*,*)' combined file: '//f2(1:25)
	k=6144*6144*4+5760
	call swap4(a,6144*6144*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 shrink9(a,n,b)
        real a(n,n),b(2048,2048)
          ii=0
        do 10 i=2,n,3
          ii=ii+1
          jj=0
        do 10 j=2,n,3
          jj=jj+1
10      b(jj,ii)=a(j,i)
        end

	subroutine edge(sx,sy,m,ix1,ix2,iy1,iy2)
	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=6144-ix2
        iy1=-iy1
        iy2=6144-iy2
        if(ix1.eq.0)ix1=1
	if(iy1.eq.0)iy1=1
	end

        function statt(a,m1,m2,m)
	real a(2048,2048)
        integer*2 b(512,512)
        y=4./(m+1)
        j1=m2
        do 10 j=1,512
        i1=m1
        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)
        statt=x/y
        end

	function allc(d,white,scale,m,sigma,s,s61,s50,g,gain)
	real d(1),white(1),scale(1),s(1),g(1)
	real d1(99)			! normal value
	real d2(99)			! bad working unit
	bad=-999.			! s(i)=6*sigma(i)/scale(i)
        frame=-998.
c bad pixel
	allc=bad
	n=0
	do 10 i=1,m
10	if(d(i).eq.bad)n=n+1			! eq
	if(n.eq.m)return                        ! if all bad, return -999.
c frame
        allc=0
	do 15 i=1,m
15	if(d(i).eq.frame)return                 ! if out of frame
c
	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.29000.0)then
          x=s50
          goto 220
        endif
c if after 20030205's data, y*=8.0/3.4
c 20030808
        y=y*g(i)/gain
210     x=x+y
220     allc=x

	end	

	subroutine readrec(k,l,b)
	real b(6144*96)
	read(k,rec=l)b
	end

	subroutine badfile(a,b,nb,lb,mb,bx,by,badtype,ix,iy,sz,k)
	real a(2048,2048),b(6144,6144)
	integer*2 bx(1),by(1),badtype,ix,iy
c check file, if <-990
        do 1 j=1,2048
        do 1 i=1,2048
1	if(a(i,j).lt.-990.)a(i,j)=-990.
c put bad point, -999.
        if(badtype.eq.0)then
	  do 10 i=1,nb
10 	  a(bx(i),by(i))=-999.
        endif
        if(badtype.eq.1)then
	  do 20 i=10001,lb
20 	  a(bx(i),by(i))=-999.
        endif
        if(badtype.eq.5)then
	  do 30 i=20001,mb
30 	  a(bx(i),by(i))=-999.
        endif
c expand pic, 2k-->6k
        do 40 i=1,2048
         i6=(i-1)*3
        do 40 j=1,2048
         j6=(j-1)*3
        x=a(i,j)
        do 40 ii=1,3
        do 40 jj=1,3
40      b(i6+ii,j6+jj)=x
c shfit, and put -998. as frame edge.
	i=iy*6144+ix
	call shift(b,i)              ! shift & y_clear
        i=-ix
        call clearx(b,i)             !         x_clear
c turn
        call t_s(b,a,sz,ix,iy)             ! a is dummy

	open(52,file=char(k/10+64)//char(mod(k,10)+64)//'.tmp',
     *  status='unknown',access='direct',recl=6144*6144*4)
	write(52,rec=1)b
	close(52)
	end	

        subroutine t_s(a,b,turn,ix,iy)
        real a(6144,6144),b(6144,6144)   ! b is dummy
        integer*2 ix,iy
        if(turn.eq.0.)return
        bad=-999.
        frame=-998.
        do 2 j=1,6144
	do 2 i=1,6144
2	b(i,j)=frame

        cosa=cos(turn)
        sina=sin(turn)
        cx=3072.5-ix
	cy=3072.5-iy
        do 10 j=1,6144
         yy=j-cy
        do 10 i=1,6144
         xx=i-cx
c        x=i-yy*turn
c        y=j+xx*turn
         y= xx*sina + yy*cosa + cy
        if(y.lt.1.or.y.gt.6144.)goto 10
         x= xx*cosa - yy*sina + cx
        if(x.lt.1.or.x.gt.6144.)goto 10
        iy0=y
        iy1=iy0+1
        wy=y-iy0
        zy=1.-wy
        ix0=x
        ix1=ix0+1
        if(a(ix0,iy0).eq.frame)goto 10
        if(a(ix0,iy1).eq.frame)goto 10
        if(a(ix1,iy0).eq.frame)goto 10
        if(a(ix1,iy1).eq.frame)goto 10
        b(i,j)=bad
        if(a(ix0,iy0).eq.bad)goto 10
        if(a(ix0,iy1).eq.bad)goto 10
        if(a(ix1,iy0).eq.bad)goto 10
        if(a(ix1,iy1).eq.bad)goto 10
        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
        do 20 j=1,6144
	do 20 i=1,6144
20	a(i,j)=b(i,j)
        return
        end

	subroutine shift(a,k)
	real a(6144*6144)
	mm=6144*6144
	if(k.gt.0)then
	  j=k
	  do 10 i=1,mm-k
	  j=j+1
10	  a(i)=a(j)
          do 11 i=mm-k+1,mm
11	  a(i)=-998.
	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)
          do 21 i=1,-k
21	  a(i)=-998.
	endif
	end

	subroutine clearx(b,ix)
	real b(6144,6144)
        if(ix.gt.0)then
	do 10 j=1,6144
	do 10 i=1,ix
10	b(i,j)=-998.
	endif
	if(ix.lt.0)then
        do 20 j=1,6144
	do 20 i=6145+ix,6144
20	b(i,j)=-998.
	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,b,ix,iy,white,px,py,pz,n,turn)
	integer*2 ix,iy,px(1),py(1)
	real a(2048,2048),b(2048,2048),pz(1)
        x=turn
        if(x.lt.0.)x=-x
        if(x.lt.0.001)then
	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
        else
c simple turn
        do 30 j=1,2048
        do 30 i=1,2048
        jy=j+(i-1024.5)*turn+0.5
        jx=i-(j-1024.5)*turn+0.5
        if(jx.lt.1)jx=1
        if(jx.gt.2048)jx=2048
        if(jy.lt.1)jy=1
	if(jy.gt.2048)jy=2048        
30      b(i,j)=a(jx,jy)
	do 40 i=1,n
	z=0.
	do 50 ii=px(i)+ix-7,px(i)+ix+7
	do 50 jj=py(i)+iy-7,py(i)+iy+7
50	z=z+b(ii,jj)-white
40	pz(i)=z
	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 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,
     c                     white,sx,sy,m,g1,gain)
	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
c 20030808
	kt=2000*gain/g1
	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