c substract 2 file
	character*80 fa, head(72)
	real a(2048,2048)
        real ww(16,16),xx(16),yy(16),xk(16),xks(16)
	real x2(2),y2(2)
	real dx(2048),dy(2048)
        integer*2 ma(128*128)
	logical logi
	ik=iargc()
	if(ik.lt.2)then
	  write(*,*)
	  write(*,*)'         ******** adjust sky_background *********'
	  write(*,*)
	  write(*,*)'    slope file.fit out.fit [!]'
	  write(*,*)
	  write(*,*)
	  stop
	endif

	call getarg(1,fa)
        k=index(fa,'.')-1
	if(k.le.0)k=lnblnk(fa)
	fa=fa(1:k)//'.fit'
	inquire(file=fa,exist=logi)
        if(.not.logi)then
	  write(*,*)fa,' not found!'
	  stop
	endif
	call gethead(fa,head)
        ipos=indexpos(head,'VOLT2   ')
        if(head(ipos)(42:46).ne."PIXEL")then
          write(*,*)fa(1:25)," do statsky first !"
          stop
        endif
        read(head(ipos)(22:),*)awhite
	if(awhite.eq.0.)stop "do statsky first !"
	write(*,*)"origon_sky: ",awhite
	awhite=awhite*1.2
	dwhite=awhite*0.03

	call readata(a,fa)
	 ik=0
	do 10 i=1,2048,128
	 ik=ik+1
         ij=0
	do 10 j=1,2048,128
         ij=ij+1
        ip=0
	do 20 ii=1,128
	do 20 jj=1,128
	x=a(ii+i,jj+j)+0.5
	if(x.gt.1..and.x.lt.32767.)then
	  ip=ip+1
	  ma(ip)=x
	endif
20	continue
	call whitexblack(ma,1,ip,white,sigma)
10      ww(ik,ij)=white
	
        do 29 i=1,16	
29      xx(i)=i

	ik=iargc()
	if(ik.gt.2)call pgbegin(0,"/xw",2,1)
	if(ik.gt.2)call pgenv(0.,17.,0.,awhite,0,0)

        do 30 i=1,16
        do 40 j=1,16
40	yy(j)=ww(i,j)-i*dwhite
	call yabx(xx,yy,1,16,xk(i),xm,xks(i),xms)
	if(ik.gt.2)call pgline(16,xx,yy)
30	continue

	yslop=slop16(xk,xks,16,ik)
        x2(1)=1.
	x2(2)=16.
	y=yy(16)-dwhite*3
        y2(1)=yslop*x2(1)+y
	y2(2)=yslop*x2(2)+y
	if(ik.gt.2)call pgsci(2)
	if(ik.gt.2)call pgline(2,x2,y2)
	if(ik.gt.2)call pgsci(1)
	

	if(ik.gt.2)call pgenv(0.,17.,0.,awhite,0,0)
        do 31 i=1,16
        do 41 j=1,16
41	yy(j)=ww(j,i)-i*dwhite
	call yabx(xx,yy,1,16,xk(i),xm,xks(i),xms)
	if(ik.gt.2)call pgline(16,xx,yy)
31	continue
	xslop=slop16(xk,xks,16,ik)
        x2(1)=1.
	x2(2)=16.
	y=yy(16)-dwhite*3
        y2(1)=xslop*x2(1)+y
	y2(2)=xslop*x2(2)+y
	if(ik.gt.2)call pgsci(2)
	if(ik.gt.2)call pgline(2,x2,y2)
	if(ik.gt.2)call pgsci(1)

	xslop=xslop*15./2047.
	yslop=yslop*15./2047.
	do 50 i=1,2048
	dx(i)=(i-1024.5)*xslop
50	dy(i)=(i-1024.5)*yslop

	do 60 i=1,2048
	do 60 j=1,2048
60	a(i,j)=a(i,j)-dx(i)-dy(j)

	call getarg(2,fa)
        k=index(fa,'.')-1
	if(k.le.0)k=lnblnk(fa)
	fa=fa(1:k)//'.fit'
        k=2048*2048*4+5760
        call swap4(a,2048*2048*4)
        open(50,file=fa,status='unknown',access='direct',recl=k)
        write(50,rec=1)head,a
        close(50)
	call system("statsky "//fa//" !")
	

	if(ik.gt.2)call pgend
	end

	function slop16(a,b,n,ik)
	real a(1),b(1)
	do 10 i=1,n-1
	do 10 j=i+1,n
	if(b(i).gt.b(j))then
	  x=a(i)
	  a(i)=a(j)
	  a(j)=x
	  x=b(i)
	  b(i)=b(j)
	  b(j)=x
	endif
10	continue
	x=0
	do 15 j=2,n
	if(b(j).gt.b(1)*2)goto 16
15	continue
16	j=j-1
	do 20 i=1,j
20	x=x+a(i)
	slop16=x/j
	if(ik.le.2)return
	do 30 i=1,j
30	write(*,*)a(i),b(i)
	write(*,*)"slope: ",slop16
	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
c        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 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
	
        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 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 YABX(AX,AY,N1,N2,A,B,AS,BS)
        REAL AX(1),AY(1)
        X=0.
        Y=0.
        XX=0.
        XY=0.
        DO 10 I=N1,N2
        A=AX(I)
        B=AY(I)
        X=X+A
        Y=Y+B
        XX=XX+A*A
10      XY=XY+A*B
        AS=N2-N1+1
        BS=AS*XX-X*X
        A=(AS*XY-X*Y)/BS
        B=(Y-A*X)/AS
        X=0.
        DO 20 I=N1,N2
20      X=X+(AX(I)*A+B-AY(I))**2
        X=X/(AS-2.)
        AS=SQRT(AS/BS*X)
        BS=SQRT(XX/BS*X)
        RETURN
        END