subroutine scans(fsamps,fsteps,nsamps,nsteps)
      implicit none
      character*32 fsamps,fsteps
      integer nsamps,nsteps
      vector tsamp,ch1samp,ch2samp
      vector tstep,ch1,ch2,cur1,cur2
      double precision tsampd(99999)
      double precision tstepd(99999)
      integer step,samp
      real ch1ave,ch2ave
      integer nave
      call dpread(fsamps,nsamps,tsampd)
      do samp=1,nsamps
        tsamp(samp)=tsampd(samp)-tsampd(1)
      enddo
      call dpread(fsteps,nsteps,tstepd)
      do step=1,nsteps
        tstep(step)=tstepd(step)-tsampd(1)
      enddo
      samp=1
      do step=1,nsteps
        nave=0
        ch1ave=0
        ch2ave=0
        do while (tsamp(samp).lt.tstep(step))
          ch1ave=ch1ave+ch1samp(samp)
          ch2ave=ch2ave+ch2samp(samp)
          nave=nave+1
          if (cur1(step)+cur2(step).gt.10.0) then
            ch1(step)=ch1ave/nave*1e2*18/(cur1(step)+cur2(step))
            ch2(step)=ch2ave/nave*1e2*18/(cur1(step)+cur2(step))
          endif
          if (samp.lt.nsamps) then
            samp=samp+1
          else
            return
          endif
        enddo
      enddo
      end

      subroutine dpread(sfile,nrecs,td)
      implicit none
      character*32 sfile
      double precision td(*)
      integer nrecs
      integer rec
      open(unit=89,file=sfile)
      do rec=1,nrecs 
        read(89,*) td(rec)
      enddo
      close(unit=89)
      end

      real function blockf(x)
      implicit none
      real x
      real cent,r0(11),h0(11)
      parameter (cent=18.88)
* Best fit to scan f, right side
*     data r0/0.28,0.50,0.60,0.95,2.00,2.50,2.90,3.30,3.60,4.20,6.00/
*     data h0/0.00,1.00,1.40,1.00,0.25,0.33,0.10,0.32,0.10,0.02,0.005/
* Best fit to scan f, left side
*     data r0/0.28,0.50,0.60,0.95,1.80,2.50,2.90,3.30,3.60,4.20,6.00/
*     data h0/0.00,1.00,1.40,1.00,0.28,0.55,0.10,0.32,0.10,0.02,0.005/
* Best fit to scan g, right side
*     data r0/0.28,0.50,0.60,0.95,2.00,2.50,2.90,3.30,3.60,4.20,6.00/
*     data h0/0.00,1.00,1.10,1.00,0.25,0.33,0.10,0.32,0.10,0.02,0.005/
* Best fit to scan g, left side
      data r0/0.28,0.50,0.60,0.95,1.80,2.50,2.90,3.30,3.60,4.20,6.00/
      data h0/0.00,1.00,1.10,1.00,0.28,0.55,0.10,0.32,0.10,0.02,0.005/
*
      real r
      integer i
      r=abs(x-cent)
      do i=1,11
        if (r.lt.r0(i)) then
          blockf=h0(i)
          return
        endif
      enddo
      blockf=0
      end

      real function fuzzf(x)
      implicit none
      real x
      real blockf
      external blockf
      real dxlim,sigma
      integer nsteps
      parameter (dxlim=1.0,sigma=0.11)
      parameter (nsteps=250)
      integer step
      real y
      fuzzf=0
      do step=1,nsteps
        y=-dxlim+2*dxlim*step/nsteps
c       fuzzf=fuzzf+blockf(x-y)*exp(-0.5*(y/sigma)**2)
c    +             *(2*dxlim/nsteps)/(sqrt(2*3.1416)*sigma)
        fuzzf=fuzzf+blockf(x-y)/(1 + (y/sigma)**2)
     +             *(2*dxlim/nsteps)/(2*atan(2*dxlim/sigma))/sigma
      enddo
      end