real b(2048)
        real aa(2048,2048)
	real bb(2048*2048+2160)	
	character*60 f1,f2
	character*80 head(108)

        equivalence (bb,head),(bb(2161),aa)

	write(*,*)
	write(*,*)'    ******** delete bad points (median)***********'
	write(*,*)
	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(*,*)' n*sigma (ex. 8*0.005=0.04 ) n= '
	read(*,*)s
	call readfits(f1,aa,n1,n2)
	 open(1,file=f1,status='old',access='direct',recl=2880*3)
	 read(1,rec=1)head
	 close(1)
	 ipos=indexpos(head,'BITPIX  ',72)
	 read(head(ipos)(24:),*)i24
	 if(i24.eq.16)then
	   call change4(aa,aa,n1,n2)
	   head(ipos)(28:30)='-32'
	   ipos=indexpos(head,'END     ',72)
	   if(ipos.le.72)then
	     do 331 i=ipos,107
331	     head(i)='                                     '
	     head(108)(1:8)='END                            '
	   endif
	 endif

	n=n1
	iz=0	
	s=s*s

	xx=xcol(aa,n,16)

	call pgbegin(0,'/xw',1,2)
	call pgenv(0.,2048.,0.9*xx,1.1*xx,0,0)
	call pglabel(' ',' ','ROW MEDIAN VALUE')
	do 50 j=1,n
	do 40 i=1,n
40	b(i)=aa(i,j)
	x=xmedian(b,n)
	sm=sigma(b,n,x)**2*s
	do 41 i=1,n
	if((aa(i,j)-x)**2.gt.sm)then
	   iz=iz+1
	  aa(i,j)=x
	endif
41	continue
50	call pgpoint(1,float(j),x,-1)

	call pgenv(0.,2048.,0.9*xx,1.1*xx,0,0)
	call pglabel(' ',' ','COLUMN MEDIAN VALUE')
	do 150 j=1,n
	do 140 i=1,n
140	b(i)=aa(j,i)
	x=xmedian(b,n)
	s=sigma(b,n,x)**2*s
	do 141 i=1,n
	if((aa(j,i)-x)**2.gt.sm)then
	   iz=iz+1
	   aa(j,i)=x
	endif
141	continue
150	  call pgpoint(1,float(j),x,-1)
	call pgend
	write(*,*)iz,' points substuded!'
	write(*,*)

29	write(*,*)' output fits file: '
	read(*,'(a)')f2
	k=index(f2,'.')
	if(k.eq.0)f2=f2(1:lnblnk(f2))//'.fit'
	inquire(file=f2,exist=i)
	if(i.ne.0)then
	  write(*,*)' file exists !'
	  goto 29
	endif
        call swap4(aa,n*n*4)
	m=n*n*4+2880*3
        open(3,file=f2,status='new',access='direct',recl=m)
        write(3,rec=1)bb
        close(3)
	end

        function xcol(a,n,m)
        real a(n,n),b(2048)
        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(300)+b(800)+b(1100)+b(1400)+b(1700))*.2
        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 sigma(a,n,x)
	real a(n)
	xx=0
	do 10 i=11,n-10
10	xx=xx+(a(i)-x)**2
	sigma=sqrt(xx/(n-20))
	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