c overscan corrected 
c only for *.fit (16bit) 
c make mark at BZERO line end, it is "(overscan corrected)" 
c if command line "overscan filename", then display only
c 1994,9,15 
c /* 94,10,5 modify BZERO --> PZERO, BSCALE --> PSCALE */
c 94,11,17 fitshead
c 94,11,28 fit length -4 bytes fitsfile
	parameter (mbyte=2048*2048*2+5764)
	parameter (nbyte=2048*2048*2)
	character*80 head(72),f1,cmark*2
	byte buf(mbyte),a(nbyte)
	logical logi
	equivalence (buf,head),(buf(5761),a)

	if(iargc().eq.0)stop ' Usage:  overscan ccdfile [!]'
	call getarg(1,f1)	! get filename
	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!'
c check file : 16/32,  done/?
	open(3,file=f1,status='old',access='direct',recl=1440*4)
	read(3,rec=1)head
	close(3)
	read(head(2)(23:),*)k
	if(k.ne.16)stop ' not raw data'
	read(head(4)(23:),*)n1
	read(head(5)(23:),*)n2
	ipos=indexpos(head,'OBSERVER')
	if(head(ipos)(52:59).eq.'OVERSCAN')
     *  stop 'already done overscan corrected'
	call filesize(f1,mm)
c readdata
	call dumpread(f1,buf,mm)
	mm=n1*n2*2+5764                    !add 4 byte at end
	call swap2(a,n1*n2*2)
	call getarg(2,cmark)
	mark=1
	if(cmark.eq.'!')mark=0
	call overscan_correct(a,n1,n2,mark)
	head(ipos)(51:70)='(OVERSCAN CORRECTED)'
        ipos=indexpos(head,'BZERO   ')
	if(ipos.lt.72)head(ipos)(1:1)='P'
        ipos=indexpos(head,'BSCALE  ')
	if(ipos.lt.72)head(ipos)(1:1)='P'
c writedata
	call swap2(a,n1*n2*2)
	call dumpwrite('1.fm',buf,mm)
	call unlink(f1)
	call system('mv 1.fm '//f1)
	write(*,*)' overscan corrected OK!'
	end

	subroutine overscan_correct(a,n1,n2,mark)
	integer*2 a(n1,n2)
	real x(2048),y(2048)
	integer*2 z(2048),u(2048),v(2048),k,iz1,iz2,kk
	real xx(2),yy(2)
	xn1=n1
	xn2=n2
	if(mark.eq.1)then
c raw image display
	  call pgbegin(0,'/xw',1,1)
	  call pgenv(1.,xn1,1.,xn2,1,1)
85	  write(*,*)' white,black: (-1,-1 for exit)'
	  read(*,*,err=85)w,b
	  if(w.lt.0.)goto 90
	  call pgxgray(a,n1,n2,1,n1,1,n2,w,b,2)
	  goto 85
90	  call pgend
	endif

c get overscan & mean value kk
	iz1=32000
	iz2=0
	w=0.
	do 10 i=1,n2
	x(i)=i	
	k=a(n1,i)
	z(i)=k
	w=w+k
	if(k.lt.iz1)iz1=k
	if(k.gt.iz2)iz2=k
10	y(i)=k
	z1=iz1
	z2=iz2
	kk=nint(w/n2)
c median51
	call median2048(z,v,51,n2) 

	if(mark.eq.1)then
c overscan curver
	  call pgbegin(0,'/xw',1,3)
	  call pgenv(0.,xn2,z1,z2,0,0)
	  call pgline(n2,x,y)
	  call pgenv(0.,xn2,z1,z2,0,0)
	  do 21 i=1,n2
21	  y(i)=v(i)
	  call pgline(n2,x,y)
	  call pgsci(2)
	  xx(1)=1
	  xx(2)=n2
	  yy(1)=kk
	  yy(2)=kk
	  call pgline(2,xx,yy)
	  call pgsci(1)
	endif

c edges modify & filter in 2nd step 
	do 40 i=1,n2
40	u(i)=kk
c to find 2 ADU dip
	n3=n2-25
	j=26
41	k=j
	do 50 i=k,n3
	if(abs(v(i)-kk).gt.1)then
	  do 51 j=i,k,-1
	  if(v(j).eq.kk)goto 52
	  if(j.eq.26)then
	    do 58 m=1,26
58	    u(m)=v(j)
	  endif
	  if(j.eq.n3)then
	    do 59 m=n3,n2
59	    u(m)=v(j)
	  endif
51	  u(j)=v(j)
52	  do 53 j=i,n3
	  if(v(j).eq.kk)goto 41
53	  u(j)=v(j)
	endif
50	continue
c statics number of pixel over average line 1 ADU
	k=0
	w=0.01
	do 54 i=26,n3
	if(v(i).lt.kk)w=w+1.
54	if(v(i).gt.kk)k=k+1
	if(k/w.lt.0.1)goto 70
	if(k.gt.300)then		! keep under average line pixel 
           if(v(26).lt.kk)then
	     do 55 i=1,26
55	     u(i)=v(26)
	   endif 	
	   if(v(n3).lt.kk)then
	     do 56 i=n3,n2
56	     u(i)=v(n3)
	   endif
	   do 57 i=26,n3
57	   if(v(i).lt.kk)u(i)=v(i)
	endif
c statics under line continue pixel
	v(n3+1)=kk
	k=0
	do 60 i=26,n3+1
	if(v(i).lt.kk)then
	  k=k+1
	  goto 60
	endif
	if(k.gt.200)then
	  do 61 j=i-1,i-k,-1	
61	  u(j)=v(j)
          if(i.eq.n3+1)then
	    do 62 j=n3,n2
62	    u(j)=v(n3)
	  endif
	  if(i-k.le.26)then
	    do 63 j=1,26
63	    u(j)=v(26)
	  endif
	endif
  	k=0
60	continue 
	  	  
70	call median2048(u,v,51,n2)

c add delta to image	
	do 20 j=1,n2
	k=z(j)-v(j)
	do 20 i=1,n1-1
20	a(i,j)=a(i,j)+k

c display
	if(mark.eq.1)then
	  do 22 i=1,n2
22	  y(i)=v(i)
	  call pgenv(0.,xn2,z1,z2,0,0)
	  call pgline(n2,x,y)
c corrected image display
	  call pgend
	  call pgbegin(0,'/xw',1,1)
	  call pgenv(1.,xn1,1.,xn2,1,1)
25	  write(*,*)' white,black: (-1,-1 for exit)'
	  read(*,*,err=25)w,b
	  if(w.lt.0.)goto 30
	  call pgxgray(a,n1,n2,1,n1,1,n2,w,b,2)
	  goto 25
30	  call pgend
	  stop 'I donot change anything!'
	endif

	end

	subroutine median2048(y,z,n,m)
	integer*2 y(2048),z(2048)
	do 10 i=1,m
10	z(i)=y(i)
	k=n/2
	n1=k+1
	n2=m-k
        int=2
11      int=int+int
        if(int.lt.n)goto 11
        int=min0(n,(3*int)/4-1)
	do 20 i=n1,n2
	ii=int
20	z(i)=median(y(i-k),n,ii)
	end

	function median(y,n,int)
        integer*2 y(n),a(199),t
	do 5 i=1,n
5	a(i)=y(i)
20      int=int/2
        ifin=n-int
        do 70 ii=1,ifin
        i=ii
        j=i+int
        if(a(i).le.a(j))goto 70
        t=a(j)
40      a(j)=a(i)
        j=i
        i=i-int
        if(i.le.0)goto 60
        if(a(i).gt.t)goto 40
60      a(j)=t
70      continue
        if(int.gt.1)goto 20
	median=a((n+1)/2)
	end
	
	subroutine dumpread(f1,buf,mbyte)
	character*(*) f1
	byte buf(mbyte)
	open(41,file=f1,status='old',access='direct',recl=mbyte)
	read(41,rec=1)buf
	close (41)
	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

	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
		
