	real a(2048,2048),b(1024*1024),c(1024*1024)
	
	character*60 f1

	character*1 c4*12
	character*80 head(36)
	write(*,*)
	write(*,*)'    ******** s98 ******'
	write(*,*)
	write(*,*)'    not included at edge 10 col & row'
	write(*,*)
	if(iargc().ne.0)then
	  call getarg(1,f1)
	  goto 10
	endif
20	write(*,'(a,$)')'Input *.fits_names: '
	read(*,'(a)')f1
10	k=index(f1,'.')
	if(k.eq.0)f1=f1(1:lnblnk(f1))//'.fit'
	inquire(file=f1,exist=i)
	if(i.eq.0)then
	  write(*,*)' file not found !'
	  goto 20
	endif
31	write(*,*)'binning: (ex. 2,4,8,16,32,64,128,256) '
	read(*,*)m
	if(m.lt.2.or.m.gt.256)goto 31
	call readfits(f1,a,n1,n2)
         open(1,file=f1,status='old',access='direct',recl=2880)
         read(1,rec=1)head
         close(1)
         ipos=indexpos(head,'BITPIX  ',36)
         read(head(ipos)(24:),*)i24
         if(i24.eq.16) call change4(a,a,n1,n2)
	n=n1
	k=0
	x=xcol(a,n,b,m)
	do 39 i=1,n
	do 39 j=1,n
39	a(i,j)=a(i,j)/x
	do 30 j=10,n-10,m
	do 30 i=10,n-10,m
	k=k+1
30	b(k)=bin(a,n,i,j,m)
	write(*,*)' total point: ',k,'   sigma=',sigma(b,k)
	write(c4(1:),"(i4,'*',i4)")m,m
	do 32 i=1,k
32	c(i)=i
	call pgbegin(0,'/xw',1,1)
	call pgsch(1.2)
	call pgenv(1.,float(k),0.95,1.05,0,0)
	call pgline(k,c,b)
	call pglabel(' ',' ','BINNING'//c4)
	call pgend
	end	

        subroutine change4(aa,bb,n1,n2)
        integer*2 aa(1)
        real bb(1)
        n=n1*n2
        do 10 i=n,1,-1
10      bb(i)=aa(i)
        end

        function xcol(a,n,b,m)
        real a(n,n),b(n)
        do 10 i=1,n
        x=0
        do 20 j=m+1,n-m
20      x=x+a(j,i)
10      b(i)=x/(n-m-m)
        xcol=(b(100)+b(200))*.5
        end


        function sigma(b,k)
        real b(k)
        x=0
        do 20 j=1,k
20      x=x+b(j)
        x=x/k
        xx=0
        do 30 j=1,k
30      xx=xx+(b(j)-x)**2
        sigma=sqrt(xx/k)
        end


	function bin(a,n,ii,jj,m)
	real a(n,n)
	real b(1024*1024)
	k=0
	do 10 i=ii,ii+m
	do 10 j=jj,jj+m
	k=k+1
10	b(k)=a(i,j)
	bin=xmedian(b,k)
	end


       function xmedian(ia,n)
       real ia(n),it
       int=2
   10  int=2*int
       if(int.lt.n)goto 10
       int=min0(n,(3*int)/4-1)
   20  int=int/2
       ifin=n-int
       do 70 ii=1,ifin
       i=ii
       j=i+int
       if(ia(i).le.ia(j))goto 70
       it=ia(j)
   40  ia(j)=ia(i)
       j=i
       i=i-int
       if(i.le.0)goto 60
       if(ia(i).gt.it)goto 40
   60  ia(j)=it
   70  continue
       if(int.gt.1)goto 20
        m=n/2
        xmedian=ia(m+1)
        if(m+m.ne.n)return
        xmedian=(ia(m)+ia(m+1))/2
       end

        function indexpos(head,f1,n)
        character*80 head(108),f1*8
        do 10 indexpos=1,n
10      if(head(indexpos)(1:8).eq.f1)return
        end
