	parameter (k4=4096)
	parameter (maxbyte=2*2880+k4*k4*4)
	parameter (nbyte=k4*k4*4)
	character*80 head(72)
	byte buf(maxbyte),a(nbyte)
	integer*2 b(k4*k4)
	character*60 f1
	logical logi
	equivalence (buf,head),(buf(5761),a)

        ii=iargc()
	if(ii.ne.0)then
	  call getarg(1,f1)
	  goto 10
	endif
	write(*,*)'  **********  static sky+balance **********2006.3'
	write(*,*)'   after run flaten, to get a smooth sky'
	write(*,*)
20	write(*,*)'Input *.fits_names: '
	read(*,'(a)')f1
10	k=index(f1,'.')
	if(k.eq.0)f1=f1(1:lnblnk(f1))//'.fit'
	inquire(file=f1,exist=logi)
	if(.not.logi)stop ' file not found !'
	call readata(f1,a,n1,n2,ii,nc1)

        if(n1.lt.4000)write(*,"('**statsky: ',a)")f1
	if(n1.ge.4000)write(*,"('**statsky+balence: ',a)")f1 
	call mean(a,n1,n2,x,nc1)        ! x is right level, dummy
	call ch2(a,b,n1,n2)	
	if(n1.ge.4000.and.n2.ge.4000)then
  	  x1=stat1(b,n1,n2,1400,1400)
	  x2=stat1(b,n1,n2,1400,2600)
	  x3=stat1(b,n1,n2,2000,2000)
	  x4=stat1(b,n1,n2,2600,1400)
	  x5=stat1(b,n1,n2,2600,2600)
c	  write(*,*)x1,x2,x3,x4,x5
	  sky=amin1(x1,x2,x3,x4,x5)
	else
	  call whitexblack(b,n1,n2,sky,s)
	endif
c	write(*,*)sky
        open(3,file=f1,status='old',access='direct',recl=1440*4)
        read(3,rec=1)head
        close(3)
        ipos=indexpos(head,'EXPOSURE')
	read(head(ipos)(25:),*)i
c	sky=sky/i
        ipos=indexpos(head,'VOLT2   ')
	write(head(ipos)(20:),"(f11.3,' / SKY ADU/PIXEL')")sky
        ipos=indexpos(head,'END     ')
        if(ipos.eq.73)head(72)(1:8)='END     '
	call swap4(a,n1*n2*4)			
	mbyte=5760+n1*n2*4
	call dumpwrite('1.fm',buf,mbyte)
	call unlink(f1)
	call system('mv 1.fm '//f1)
	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

	function stat1(a,n1,n2,m1,m2)
	integer*2 a(n1,n2),b(512,512)
	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)
	call whitexblack(b,512,512,x,s)
	stat1=x
	end
	
	subroutine readata(f1,a,n1,n2,ii,nc1)
	byte c(2880),a(1)
	character*60 f1,ch*2880
	equivalence (c,ch)
	nn=2880
	is=0
	open(3,file=f1,status='old',access='direct',recl=nn)	
	read(3,rec=1)ch
	k=index(ch,'BITPIX')+26
	read(ch(k:),'(i4)')ibit
        if(ibit.eq.16)stop' file must be p*.fit'
	k=index(ch,'NAXIS1')+26
	read(ch(k:),'(i4)')n1
	k=index(ch,'NAXIS2')+26
	read(ch(k:),'(i4)')n2
        k=index(ch,'CRVAL1')+26
	read(ch(k:),'(i4)')nc1
	k=index(ch,'VOLT2')+41
	if(ch(k:k+2).eq.'PIX' .and. ii.lt.2)stop 'already done !'
	n=1
55	k=index(ch,'END          ')
	if(k.eq.0)then
	n=n+1
	read(3,rec=n)ch
	goto 55
	endif
	i=0
	do 60 k=n+1,n+(n1*n2*4-1)/nn+1
	read(3,rec=k,err=61)ch
	do 60 j=1,nn
	i=i+1
60	a(i)=c(j)
61	close (3)
	call swap4(a,n1*n2*4)
	end

	subroutine ch2(a,b,n1,n2)
	real a(n1,n2)
	integer*2 b(n1,n2)
	do 10 j=1,n2
	do 10 i=1,n1
	k=a(i,j)
	if(k.lt.0)k=0
	if(k.gt.32767)k=32767
10	b(i,j)=k
	end
		

	subroutine mean(a,n1,n2,x,nc1)
c take 2:2,2:3,3:2,3:3  1/4 of whole ccd frame do statistic, mode method
	real a(n1,n2)
	integer*4 ihist(30000)
	x=0
	if(nc1.ge.2048)return
	if(nc1+n1.lt.2048)return

        do 10 i=1,30000
10      ihist(i)=0
        n11=n1/2-512
	n22=n2/2+512
	if(n11.lt.1)n11=1
	if(n22.gt.n2)n22=n2
        n3=2048-nc1
	n4=n3-512
	if(n4.lt.1)n4=1
	n5=n3+512+1
	if(n5.gt.n1)n5=n1
c n4--n3, n3+1,n5    1--2048, 2049-4096
        i=n3-n4
	j=n5-(n3+1)
	if(i.gt.j)n4=n3-j
	if(j.gt.i)n5=n3+i+1
        n=0
	nn=0
	do 70 j=n11,n22
	do 70 i=n4,n3
        k=a(i,j)*.5+101
        if(k.lt.1.or.k.gt.30000)goto 70
        ihist(k)=ihist(k)+1
        n=max(k,n)
	nn=nn+1
70      continue
	if(nn.lt.1000)return
        call histat(ihist,mode,maxh,xmean,xpeak,sigma,gsigma,n)
        xleft=(xpeak-100.0)*2.

        do 110 i=1,30000
110     ihist(i)=0
        n=0
	nn=0
	do 170 j=n11,n22
	do 170 i=n3+1,n5
        k=a(i,j)*.5+101
        if(k.lt.1.or.k.gt.30000)goto 170
        ihist(k)=ihist(k)+1
        n=max(k,n)
	nn=nn+1
170     continue
	if(nn.lt.1000)return
        call histat(ihist,mode,maxh,xmean,xpeak,sigma,gsigma,n)
        x=(xpeak-100.0)*2.

	y=x/xleft
	do 270 j=1,n2
	do 270 i=1,n3
270     a(i,j)=a(i,j)*y

	write(*,*)xleft,x,y,n3

	end

	subroutine dumpwrite(f1,buf,mbyte)
	character*(*) f1
	byte buf(mbyte)
	call unlink(f1)
	open(41,file=f1,status='new',access='direct',recl=mbyte) 
	write(41,rec=1)buf
	close (41)
	end
