c   .    1    .    2    .    3    .    4    .    5    .    6    .    7..
c
c     DISPCAL9
c     absolute calibration of seismometers by stepwise displacement
c
c     November 1998: ew@geophys.uni-stuttgart.de
c     Revisions:
c     July 2002: threshold for straight sections determined
c                 from heap-sorted straightness samples
c     Sep. 2002: five trial values for the gap parameter
c     Oct. 2005: straightness limit (zlim) set at the steepest slope
c                 of the sorted series
c                slope of velocity trace included in the straightness
c                 criterion so that a trace with a constant acceler-
c                 ation is not considered as representing rest
c                option "accelerometer" (int in place of sez
c                 when free period is specified as zero)
c     June 2006: read ASL-format files; name of parameter file may be
c                 entered as runstring; outputfile named after parameter
c                 file; other minor improvements
c     Jan. 2008: the plot-parameter file is named winplot.par
c     April 2008: the previously diverging versions "Edi" and "ASL" were
c                  recombined (as were the revision histories)
c                 The definition of the straighness threshold was improved
c                 Low-pass filtering is now 2nd order forward and backward
c                  (the order parameter in the parameter file is ignored).
c                  The determination of the straightness limit was improved. 
c                 Minimum safety distance is 1.5 times filter corner period
c                 Double precision is used throughout.
c                 Free-format reading (ASL format) has been reprogrammed
c                  because my FORTRAN compiler does not handle
c                  free-format reading properly.
c     Nov. 2008:  deconvolution is skipped when ts0 < 0. This is for
c                 calibrating bare velocity transducers, or analyzing
c                 data that have already been deconvolved.
c
c     March 2009: median statistics added
c     Sep. 2010:  'perc'entage of "quiet" time as optional parameter
c     Aug. 2011:  zlim redefined as (logarithmic) average of median and 97%ile. 
c                 'perc' parameter is obsolete and deleted deleted.           
c
c   There MUST BE AT LEAST 6 quiet seconds before the first pulse!
c
c   convention for the 'seife' type file format:
c   - one header line, up to 40 characters
c   - up to 48 lines of signal description and processing log
c     these lines have '%' as the first character
c   - one line specifying the number of samples, the format, the
c     sampling interval, and the start time in min and sec after
c     midnight. This line may be repeated after the data when its
c     contents is not known at the time when the file header is
c     written; however the format must always be specified.
c

      program dispcal8 
      use nas_system
      implicit double precision (a-h,o-z)
      parameter(ndim=480000)
      dimension x(ndim),y(ndim),xx(ndim),z(ndim)
      dimension sstep(199),step(199),step1(199),ja(199),na(199),ne(199)
      character text*72,zeile*40,par*36,unit*7
      character nam1*50,nam2*50,nam3*50,nam4*50
      common /avstep/ avs,estep
      logical trac,pulselist,onegap,acc,decon
      common /tra/ trac
      common /tre/ a,b

      write(6,*)
      write(6,*) '  This is program DISPCAL, version 9, August 2011'     
      write(6,*) '  reading data in SEIFE and Cimarron-ASL format'
      write(6,*) '  output files are in SEIFE format and can' 
      write(6,*) '  automatically be plotted by calling WINPLOT'
      write(6,*) '  Name of parameter file may be entered as runstring'
      write(6,*) '  Protocol file is named after parameter file'
      write(6,*)
      write(6,*) '=== There MUST BE AT LEAST 5 quiet seconds before the
     & first pulse! ==='
      write(6,*)

      npol=3
      trac=.false.
      n=ndim
      if(iargc().gt.0) then
        call getarg(1,nam3)
      if(index(nam3,'.par').eq.0) then
        write(6,1) ' stop - ',
     &' the name of the parameter file must have the extension .par'
        stop
      endif
      else
        nam3='dispcal.par'
      endif
      open(11,file=nam3,status='old')
      read(11,*) text
      write(6,'(1x,a)') text
      read(11,*) nam1
      ipunkt=index(nam3,'.')
      nam4=nam3(1:ipunkt)//'out'
      open(3,file=nam4)
      write(3,1) ' This is program DISPCAL, version 9, August 2011'
      write(3,'(1x,a)') text
      write(6,'(a,a)') ' Data File:   ',nam1
      write(3,'(a,a)') ' Data File:   ',nam1
      read(11,*) text
      write(6,'(a,a)') ' Seismometer: ',text
      write(3,'(a,a)') ' Seismometer: ',text
      read(11,*) ts0
      acc=ts0.eq.0.
      decon=ts0.ge.0.
      read(11,*) hs0
      read(11,*) volts
      read(11,*) displ
      read(11,'(a40)') zeile     
      read(zeile,*) gap      
      read(11,*) tfilt
      read(11,*) pulselist
      read(11,*) corr
      close(11)
 
      fac=volts/displ*0.001*corr
 
      write(6,'(a,a,a)') ' Parameters read from file ',trim(nam3),':'
      write(3,'(a,a,a)') ' Parameters read from file ',trim(nam3),':'
      write(6,2020) 'ts0=',ts0,'   hs0=',hs0,'   volts=',volts,
     &'   displ=',displ,'   gap=',gap
      write(3,2020) 'ts0=',ts0,'   hs0=',hs0,'   volts=',volts,
     &'   displ=',displ,'   gap=',gap
 2020 format(1x,a,f6.2,a,f5.3,a,f11.3,a,f6.3,a,f5.3,a,f3.0)
      write(6,2021) ' Butterworth low-pass filter To=',tfilt,' s'
      write(3,2021) ' Butterworth low-pass filter To=',tfilt,' s'
 2021 format(a,f7.2,a,i1)
 
       if(gap.gt.0.) then
        onegap=.true.
        ngap=1
      else
        gap=0.125
        onegap=.false.
        ngap=7
      endif
      if(gap.lt.1.5*tfilt) then
      gap=1.5*tfilt
      write(6,2020) 'safety distance ''gap'' was increased to ',gap,' s'
      write(3,2020) 'safety distance ''gap'' was increased to ',gap,' s'
      endif
      write(6,'(a,f6.3)') ' Correction factor (source+load)/load =',corr
      write(3,'(a,f6.3)') ' Correction factor (source+load)/load =',corr
      write(6,1) ' Parameters set in the source code:'
      write(3,1) ' Parameters set in the source code:'
      write(6,2022) ' ndim=',ndim,'   npol=',npol
      write(3,2022) ' ndim=',ndim,'   npol=',npol
 2022 format(a,i7,a,i1,a,f5.1,a,i2,a,i2,a)
      call input(nam1,text,n,dt,tmin,tsec,x)
      if(n.eq.0.or.dt.eq.0.) then
      write(6,1) 'n or dt unspecified - check data header'
      write(3,1) 'n or dt unspecified - check data header'
      stop
      endif

c  Data preparation: deconvolution to ground velocity
c  remove average determined from first 6 seconds
      do i=1,n
        z(i)=1.
      enddo
      
      call average(x,nint(4./dt),n)     
c  store a copy of the input signal for plotting
      nam2='dispcal.dat'
      call outform(nam2,text,n,dt,tmin,tsec,x)
 
 2001 format(a3,2x,i5,30x)


      if(tfilt.gt.0.) then
        write(zeile,2004) 'lp2',tfilt,0.9
 2004   format(a3,2x,f8.3,2x,f5.2,20x)
        write(6,*) 'forward filtration'
        write(3,*) 'forward filtration'
        call filter(zeile,x,n,dt)
        call reverse(x,n)
        write(6,*) 'backward filtration'
        write(3,*) 'backward filtration'
        call filter(zeile,x,n,dt)
        call reverse(x,n)
      endif
 
      if(decon) then
c       deconvolution
        if(.not.acc) then
          write(6,*) 'conversion to ground velocity'
          write(3,*) 'conversion to ground velocity'
          if(hs0.gt.0.) then
            write(zeile,2002) 'sez',ts0,hs0
 2002       format(a3,2x,2(f9.3,1x),15x)
          else
            write(zeile,2002) 'he1',ts0,1e6
          endif
        else
          write(zeile,2005) 'int'
 2005     format(a3,17x)
        endif  
        call filter(zeile,x,n,dt)
        write(6,*) 'removing polynomial trend'
        write(3,*) 'removing polynomial trend'
        call polytrend(npol,x,n,z)
      else
        write(6,*) 'no deconvolution'
        write(3,*) 'no deconvolution'
      endif

c  store deconvolved signal for visual check
      nam2='dispcal.v_1'
c     call outform(nam2,text,n,dt,tmin,tsec,x)
      do i=1,n
        xx(i)=x(i)
      enddo
 
c   x is now the ground-velocity signal, will be detrended and integrated
c   to give a displacement signal representing the steps
 
      write(6,*) 'measuring curvature of the velocity trace'
      write(3,*) 'measuring curvature of the velocity trace'
 
c Loop over five trial values of gap if gap = 0
 
      if(.not.onegap) then
      write(6,1)"======================================================"
      write(3,1)"======================================================"
    1 format(a,a,a,a,a)
      endif

      do 3000 ntry=1,ngap
      if(.not.onegap) then
      trac=.false.
      write(6,12) gap
      write(3,12) gap
   12 format(" distance parameter 'gap' set to ",f5.2," sec")
      endif
      avglen=5.*gap
      iab=nint(gap/dt)
      iavg=nint(avglen/dt)
      do i=1,n
        x(i)=xx(i)
        z(i)=1.
        y(i)=0.
      enddo
      do i=1,199
        ja(i)=1
      enddo
 
c Search for quiet (linear) sections of minimum length 2*iab+1 samples
 
      call krum(x,n,2*iab+1,y)
      pmin=1.e12
      do i=1,n
        pmin=min(pmin,y(i))
        y(i)=dlog10(max(y(i),1.e-6))
        z(i)=y(i)
      enddo
      if(pmin.le.1.e-6) then
        write(6,*) 'WARNING *************************************'
        write(6,*) 'This program may not work with noiseless data'
      endif
      
      call heapsort(n,z)
      zlim=(z(nint(0.5*n))+z(nint(0.97*n)))/2.
      
      do i=1,n
        y(i)=y(i)-zlim
        z(i)=z(i)-zlim
      enddo
 
c     store quietness signals for visual check
      nam2='dispcal.str'
      if(onegap) call outform(nam2,text,n,dt,tmin,tsec,y)
c      nam2='dispcal.srt'
c      if(onegap) call outform(nam2,text,n,dt,tmin,tsec,z)
 
      write(6,1190) iab*dt,iavg*dt
      write(3,1190) iab*dt,iavg*dt
 1190 format(' Quiet intervals: min. length ',f5.2,
     & ' s, max.',f6.2,' s evaluated')

c  z=1 if quiet, z=0 if pulse
      do 1005 i=1,n
        if(y(i).lt.0.) then
          z(i)=1.
        else
          z(i)=0.
        endif
 1005 continue
      z(n-1)=0.
      nam2='dispcal.bol'
c      if(onegap) call outform(nam2,text,n,dt,tmin,tsec,z)


c     remove polynomial trend from quiet segments (correcting
c     the previous trend removal that included the pulses)
c     this is done by using the z series as a weight
      if(onegap) write(6,*) 
     & 'removing polynomial vel. trend determined from quiet intervals'
      if(onegap) write(3,*) 
     & 'removing polynomial vel. trend determined from quiet intervals'
      call polytrend(npol,x,n,z)
c  store restored velocity signal and residual noise for visual check
       nam2='dispcal.vel'
       if(onegap) call outform(nam2,text,n,dt,tmin,tsec,x)
       
c  blank out velocity pulses and store series for a plot       
      do 1120 i=1,n
 1120   y(i)=x(i)*z(i)
      nam2='dispcal.w_2'
c      if(onegap) call outform(nam2,text,n,dt,tmin,tsec,y) 

      if(onegap) then
        write(6,*) 'locating pulses of velocity (steps of displacement)'
        write(3,*) 'locating pulses of velocity (steps of displacement)'
      endif
 
c  finding the steps. Before looping, determine first quiet interval.
      if(onegap.and.pulselist) write(6,'(a)') 'Pulses: times in seconds'
      if(onegap.and.pulselist) write(3,'(a)') 'Pulses: times in seconds'
      ns=0
      trac=.false.
      i=1
 1060 do while (i.lt.n.and.z(i).lt.0.5)
          i=i+1
      enddo
      n1=i
c  this is the first sample of the first quiet interval
      do while (i.lt.n.and.z(i).ge.0.5)
         i=i+1
      enddo
      n2=i-1
c  this is the last sample of first quiet interval
      if(i.eq.n) then
        ns=0
        goto 9000
      endif
      if(n2.lt.n1+iab) goto 1060
      
c  loop over steps begins here
 1050 continue
       do while (i.lt.n.and.z(i).lt.0.5)
        i=i+1
      enddo
      if(i.eq.n) goto 1080
      n3=i
c  first sample of next quiet interval
      do while (i.lt.n.and.z(i).ge.0.5)
         i=i+1
      enddo
      if(i.eq.n) goto 1080
      n4=i-1
c  last sample of next quiet interval
      if(n4.lt.n3+iab) goto 1050
      ns=ns+1

c  store indices where steps begin and end as na and ne
      if(onegap.and.pulselist) then
      write(6,2006) 'pulse #',ns,'  from ',(n2-1)*dt,'   to ',
     & (n3-1)*dt,' quiet to ',(n4-1)*dt
      write(3,2006) 'pulse #',ns,'  from ',(n2-1)*dt,'   to ',
     & (n3-1)*dt,' quiet to ',(n4-1)*dt
 2006 format(a,i3,3(a,f7.2))
      endif
      na(ns)=n2
      ne(ns)=n3
      n1=n3
      n2=n4
      if(i.lt.n-iavg.and.ns.lt.199) goto 1050
c  loop over steps ends here
 1080 continue

      write(6,*) '============ ',ns,' pulses found.',' ============'
      write(3,*) '============ ',ns,' pulses found.',' ============'
      if(ns.lt.3) then
        write(6,1081) 
        write(3,1081) 
 1081 format(' At least 3 pulses are required. Check the perc parameter.
     &')
        stop
      endif

c  Construct a hypothetical baseline y for the velocity signal x,
c  assuming that the tilt-related acceleration is constant in the 
c  quiet intervals, and the resulting linear trend of the velocity 
c  in these intervals must be removed. The baseline is linearly 
c  interpolated in between.

      do 1090 j=1,n
 1090   y(j)=0.d0

c  in the next section: 
c  n1= start of pulse # k-1
c  n2 = end of pulse # k-1, =1 for k=1
c  n3 = start of pulse # k
c  xa = y(n1) = y(previous n3)
c  xe = y(n2)
c  note that a and b come from the last call of trend
      n2=1
      n3=na(1)
      call trend('nix',x,n2,n3)
      do j=n2,n3
        y(j)=a+b*(j-n2+1)
      enddo
      xa=a+b*(n3-n2+1)
      
      do k=2,ns
        n1=na(k-1)
        n2=ne(k-1)
        n3=na(k)
        call trend('nix',x,n2,n3)
        do j=n2,n3
	   y(j)=a+b*(j-n2+1)
	enddo
	xe=a+b
c   interpolate baseline of pulse # k-1
      dn=n2-n1
      dx=xe-xa
        do j=n1,n2
          y(j)=xa+dble(j-n1)/dn*dx
        enddo
        xa=a+b*(n3-n2+1)
      enddo
      
c   detrending the last quiet interval
      n2=ne(ns)
      call trend('nix',x,n2,n)
      do j=n2,n
        y(j)=a+b*(j-n2+1)      
      enddo
      xe=a+b
c   baseline of pulse # ns
      n1=na(ns)
      n2=ne(ns)
      dn=n2-n1
      dx=xe-xa
      do j=n1,n2
        y(j)=xa+dble(j-n1)/dn*dx
      enddo

      nam2='dispcal.vtr'
c      call outform(nam2,text,n,dt,tmin,tsec,y)

      if(onegap) then
        write(6,*) 'removing hypothetical baseline from velocity trace'
        write(3,*) 'removing hypothetical baseline from velocity trace'
        endif
      do j=1,n
        x(j)=x(j)-y(j)
      enddo

      nam2='dispcal.v_3'
c      call outform(nam2,text,n,dt,tmin,tsec,x)
      
      if(onegap) then
        write(6,*) 'integrating "cleaned" velocity to displacement'
        write(3,*) 'integrating "cleaned" velocity to displacement'
      endif
      write(zeile,2000) 'int'
 2000 format(a3,37x)
      call filter(zeile,x,n,dt)
      call average(x,n,n)
      
      nam2='dispcal.dis'
      if(onegap) call outform(nam2,text,n,dt,tmin,tsec,x)           
            
c  now, finally, we measure the steps. n1 .. n4 are redefined symmetric
c  to pulse k.
          
        if(onegap) write(6,*) 'measuring the steps of displacement'
        if(onegap) write(3,*) 'measuring the steps of displacement'

        do k=1,ns           
        if(k.eq.1) then
          n1=1
        else
          n1=ne(k-1)
        endif
        n2=na(k)
        n3=ne(k)
        if(k.eq.ns) then
          n4=n
        else
          n4=na(k+1)
        endif       
        n1a=max(n1,n2-iavg)
        n4a=min(n4,n3+iavg)
        navg=min(n2-n1a,n4a-n3)
        n1a=n2-navg
        n4a=n3+navg
  
        write(par,2003) n1a,n2,n3,n4a
 2003   format(4i9) 
        call mspline(par,x,n)
c       write(par,2003) n1a,n2,n3,n4a
c       call steps(par,x,n) ! alte Methode, unnoetig kompliziert
c       step(k)=dabs(avs) ! alte Methode
        sstep(k)=estep
        step(k)=dabs(estep)
        
        if(onegap.and.pulselist) then
        write(6,'(a,i3,a,f16.6,a)') "Step # ",k," has ",estep," counts"
        write(3,'(a,i3,a,f16.6,a)') "Step # ",k," has ",estep," counts"
        endif
      enddo
      
      if(onegap) then
        write(6,*) 'preparing the residual displacement for display'
        write(3,*) 'preparing the residual displacement for display'
      endif
      do i=1,ns-1
        call trend('tre',x,ne(i),na(i+1))
      enddo
      do 1085 i=1,na(1)
 1085   x(i)=0.
      do 1086 i=ne(ns),n
 1086    x(i)=0.
      do i=1,n
        x(i)=x(i)*z(i)
      enddo
      nam2='dispcal.res'
      if(onegap) call outform(nam2,text,n,dt,tmin,tsec,x)
    
c  Statistics
      total=0.
      do k=1,ns
        total=total+ja(k)*step(k)
        step1(k)=step(k)
c        write(*,*) k, step(k)*fac
      enddo
      avgs=total/ns
      total=0.
      do k=1,ns
        total=total+ja(k)*(step(k)-avgs)**2
      enddo
      sigma=sqrt(total/ns)
      if(onegap) then
      write(6,*)
      write(3,*)
      write(6,1163) ' Raw average abs. step:  ',avgs,' +-',sigma,
     &' counts rms'
      write(3,1163) ' Raw average abs. step:  ',avgs,' +-',sigma,
     &' counts rms'
 1163 format(a,f15.3,a,f14.3,a)
      endif
      	        
      call heapsort(ns,step1)
      nsh=int(ns/2)
      
      if(ns.eq.2*nsh) then
        amed=0.5*(step1(nsh)+step1(nsh+1))*fac
      else
        amed=step1(nsh+1)*fac
      endif

      if(acc) then
        write(6,1161) avgs*fac,sigma*fac
        write(3,1161) avgs*fac,sigma*fac
 1161 format(' Raw avg. generator constant:    ',f10.4,' +-',f9.4,
     & ' Vs^2/m')
        write(6,1164) amed
 	write(3,1164) amed
 1164   format(' Raw median generator constant: ',f11.4,' Vs^2/m')        
      else
        write(6,1159) avgs*fac,sigma*fac
        write(3,1159) avgs*fac,sigma*fac
 1159 format(' Raw avg. generator constant:    ',f10.2,' +-',f9.2,
     & ' Vs/m')
         write(6,1165) amed
 	write(3,1165) amed
 1165 format(' Raw median generator constant: ',f11.2,' Vs/m') 
      endif

      if (ns.gt.2) then
        if(onegap) write(6,*)
        if(onegap) write(3,*)
        do 1150 nrest=ns-1,int((ns+1)/2),-1
          siga=sigma
          offmax=0.
            do 1160 k=1,ns
              off=ja(k)*abs(step(k)-avgs)
              if(off.gt.offmax) then
                offmax=off
                kmax=k
              endif
 1160       continue
          ja(kmax)=0
          total=0.
          do k=1,ns
            total=total+ja(k)*step(k)
          enddo
          avgs=total/nrest
          total=0.
          do k=1,ns
            total=total+ja(k)*(step(k)-avgs)**2
          enddo
          sigma=sqrt(total/nrest)
        if(acc) then
          if(onegap) write(6,1172) kmax,nrest,avgs*fac,sigma*fac
          if(onegap) write(3,1172) kmax,nrest,avgs*fac,sigma*fac
        else
          if(onegap) write(6,1170) kmax,nrest,avgs*fac,sigma*fac
          if(onegap) write(3,1170) kmax,nrest,avgs*fac,sigma*fac
        endif
 1172   format(' Step #',i4,' elim.,',i3,' remain.  ',
     &  'Gen.Const.=',f10.4,'  +-',f9.4,' Vs^2/m')
 1170   format(' Step #',i4,' elim.,',i3,' remain.  ',
     &  'Gen.Const.=',f10.2,'  +-',f9.2,' Vs/m')
          if(sigma.lt.0.05*avgs.and.siga-sigma.lt.siga/nrest) exit
 1150 continue
      nrest=max(nrest,int((ns+1)/2))
      if(onegap) write(6,*)
      if(onegap) write(3,*)
      if(acc) then
        write(6,1173) nrest,avgs*fac,sigma*fac
        write(3,1173) nrest,avgs*fac,sigma*fac
 1173 format('                          Best ',i3,' Steps: ',f10.4,
     & '  +-',f9.4,' Vs^2/m'/'=================',
     & '=======================================================')
      else
        write(6,1171) nrest,avgs*fac,sigma*fac
        write(3,1171) nrest,avgs*fac,sigma*fac
 1171 format('                          Best ',i3,' Steps: ',f10.2,
     & '  +-',f9.2,' Vs/m'/'===============',
     & '=======================================================')
      endif
      endif
      
 9000 continue      
           
      if(onegap) then
      write(6,*)
      write(3,*)
      
c  step table special: is there a systematic difference between up-steps and down-steps?
      avgup=0.
      nup=0
      avgdn=0.
      ndn=0
      do k=1,ns
        if(sstep(k).gt.0.) then
          avgup=avgup+ja(k)*sstep(k)
          nup=nup+ja(k)
        else
          avgdn=avgdn+ja(k)*sstep(k)
          ndn=ndn+ja(k)
        endif
      enddo
      
      if(nup.gt.0.and.ndn.gt.0) then
        if(acc) then
          unit=' Vs/m^2'
        else
          unit=' Vs/m  '
        endif
        avgup=avgup/nup
        avgdn=avgdn/ndn
        write(6,1200) avgup*fac,unit
        write(3,1200) avgup*fac,unit
        write(6,1201) avgdn*fac,unit
        write(3,1201) avgdn*fac,unit
 1200   format(' Average of good   up-steps is ', f10.2,a)
 1201   format(' Average of good down-steps is ', f10.2,a)
 
      sigold=sigma*fac
      sigma=0.
      devgut=0.
      devall=0.
      do k=1,ns
        devall=max(devall,step(k)-avgs)
        devgut=max(devgut,ja(k)*(step(k)-avgs))
        if(sstep(k).gt.0.) then
          sigma=sigma+ja(k)*(sstep(k)-avgup)**2
        else
         sigma=sigma+ja(k)*(sstep(k)-avgdn)**2
        endif
      enddo
      sigma=sqrt(sigma/(nup+ndn))*fac

        if(sigma.lt.0.82*sigold) then
          write(6,1202) sigma,unit
          write(3,1202) sigma,unit
 1202   format(' When the asymmetry is substantial, it may indicate sens
     &or nonlinearity.'/' Without it, the scatter of the generator const
     &ant is +- ',f8.2,a)
      endif
        write(6,*)
        write(3,*)
      endif
      
c  produce time series of step variance
      nee=(na(1)+ne(1))/2
      do i=1,nee-1
        xx(i)=0.d0
      enddo
c  set clip level in percent for plot of step variance
      clip=1.d0
      do k=1,ns-1
        naa=nee
        nee=(na(k+1)+ne(k+1))/2
        xxx=(step(k)/avgs-1.d0)*1.d2
        xxx=max(-clip,min(clip,xxx))
        do i=naa,nee
          xx(i)=xxx
        enddo
      enddo
      xxx=(step(ns)/avgs-1.d0)*1.d2
      xxx=max(-clip,min(clip,xxx))
      do i=nee,n
        xx(i)=xxx
      enddo 
      nam2='dispcal.stp'
      call outform(nam2,text,n,dt,tmin,tsec,xx)
      
c  write a table of step sizes as file disptable.lst
      open(9,file='disptable.lst')
      pos=0.d0
      do k=1,ns        
        steptime=(na(k)+ne(k))/2.*dt
        stepmm=sstep(k)/avgs*displ
        pos=pos+stepmm
        write(9,1180) k,steptime,stepmm,pos
      enddo
        write(9,1180) 0,0.,0.,0.
      write(9,1181) displ
      close(9)
 1180 format(i10,f10.1,2f10.4)
 1181 format('    step #    time s   step mm    pos mm'/
     &'Note: The average step size was adjusted to',f8.4/
     &'      This table is good for step variance only.')      
      write(*,*) 'Steps are listed in file disptable.lst'
      write(*,*) 'time series dispcal.stp is in percent of avg, clipped
     & at 1%'
      
c  generate plot-parameter file
      open(4,file='winplot.par')
      if(ns.eq.0) then
        write(4,*) '0, 4, 24, 15, 1,',n, ', 0.8'
        write(4,4) 'dispcal.dat'
        write(4,4) 'dispcal.v_1'
        write(4,4) 'dispcal.str'
        write(4,4) 'dispcal.bol'
       else
        write(4,*) '0, 6, 24, 15, 1,',n, ', 0.8'
        write(4,4) 'dispcal.dat'
        write(4,4) 'dispcal.vel'
        write(4,4) 'dispcal.str'
        write(4,4) 'dispcal.dis'
        write(4,4) 'dispcal.res'
        write(4,4) 'dispcal.stp'
      endif
      write(4,4) 'Program DISPCAL v.8    Data File: '//trim(nam1)
    4 format(a)
      close(4)
      write(6,*) 'plot-parameter file winplot.par was generated'
      write(6,*) 'type ''winplot'' to automatically plot the signals'
c      write(6,*) 'You may also plot dispcal.srt and dispcal.bol'
      write(6,*) 'results were saved to file ',trim(nam4)
      endif
      write(6,*)"======================================================"
      write(3,*)"======================================================"
      gap=gap*sqrt(2.)
 3000 continue
      stop
      end
 
      subroutine input(name,text,n,dt,tmin,tsec,x)
      implicit double precision(a-h,o-z)
      parameter(ndim=480000)
      dimension x(ndim),ix(60)
      character iform*20,name*50,text*72,code*12
      character zeile*120,zarr(120),b
      logical seif
      b=' '
      write(6,*) 'Opening file ',trim(name)
      open(7,file=name,status='old')
      read(7,'(a)') zeile
      read(7,'(a)') zeile
      seif=index(zeile,'%').gt.0.or.index(zeile,'(').gt.0
      if(seif) then
        write(6,*) 'file ',trim(name),' assumed to be in SEIFE format'
      else
        write(6,*) 'file ',trim(name),' assumed to be in ASL format'
      endif
      rewind 7
      if(seif) then
c read data in SEIFE format 
        read(7,'(a)') text
        write(6,*) 'header: ',trim(text)
        write(3,*) 'header: ',trim(text)
   21   read(7,'(a)') zeile
        if(zeile(1:1).eq.'%') goto 21
   20   read(zeile,1) n,iform,dt,tmin,tsec
    1   format(i10,a20,3f10.3)
        write(6,22) n,trim(name)
        write(3,22) n,trim(name)
   22   format(' reading ',i7,' samples from file ',a)
   
        read(7,iform,err=25,end=23) (x(j),j=1,n)
        
        close(7)
	return
   23   n=j-1
        write(6,24) n
        write(3,24) n
   24   format('end of file after ',i7,' samples')
        close(7)
        return
   25   write(6,26) j
        write(3,26) j
   26   format(' Input error (non-numeric character?) at sample # ',i7)
      else
c read data in ASL format (such as written by Quanterra's Cimarron)
        read(7,'(a)') zeile
        text=zeile(1:72)
        write(6,*) 'header: ',trim(text)        
        read(zeile,'(120a)') zarr
	nn=0
	do i=1,119
	  if(zarr(i).ne.b.and.zarr(i+1).eq.b) then
	    nn=nn+1
	    ix(nn)=i
	  endif
	enddo
        read(zeile(1:ix(2)),'(a)') code
        read(zeile(ix(2)+2:ix(7)),*) iyear,iday,ithr,itmin,itsec
        read(zeile(ix(7)+2:ix(10)),*) t1000,srate,n        
c        read(text,30) code,iyear,iday,ithr,itmin,itsec,t1000,srate,n
c   30   format(a12,1x,i2,1x,i3,3(1x,i2),1x,f7.3,f7.2,i6)
        thr=ithr
        tmin=itmin
        tsec=itsec
        dt=1./srate
        tmin=tmin+60.*thr
        tsec=tsec+t1000/1000.
c  ???  read(7,*) (x(j),j=1,n)
c   schoen waer's, geht aber nicht zuverlaessig - Bug im Compiler?
        nz=1
        nval=0
        do
          read(7,'(a)',end=31) zeile
c         read(zeile,'(120a)') zarr
          nn=0
          do i=1,119
c           if(zarr(i).ne.b.and.zarr(i+1).eq.b) nn=nn+1
            if(zeile(i:i).ne.b.and.zeile(i+1:i+1).eq.b) nn=nn+1
          enddo
          read(zeile,*,err=32) (ix(j),j=1,nn)
          do j=1,nn
            x(nval+j)=ix(j)
          enddo
          nval=nval+nn
        enddo
   31   write(6,*) nval,' samples read from file ',trim(name)
        write(3,*) nval,' samples read from file ',trim(name)
        n=min(n,nval)
        close(7)
        return
   32   write(6,*) 'error (maybe a non-numeric character) in line ',nz,
     &  ' of file ',trim(name)
        stop
      endif
      end
      

      subroutine outform(name,text,n,dt,tmin,tsec,x)
c  formatted output
      implicit double precision (a-h,o-z)
      parameter(ndim=480000)
      dimension x(ndim)
      character name*50,text*72,iform*20
      xmax=0.
      do 2 j=1,n
    2 xmax=max(xmax,abs(x(j)))
      nvor=int(log10(max(xmax,2.))+1.)
      ndec=max(0,10-nvor)
      write(iform,3) ndec
    3 format('(5f13.',i1,')')
      write(6,5) n,trim(name)
      write(3,5) n,trim(name)
    5 format(' writing ',i7,' samples into file ',a)
      open(8,file=trim(name))
      write(8,'(a)') text
      write(8,1) n,iform,dt,tmin,tsec
    1 format(i10,a20,f10.3,2f10.3)
      write(8,iform) (x(j),j=1,n)
      close(8)
      return
      end
 

      subroutine filter(zeile,x,n,dt)
      implicit double precision (a-h,o-z)
      parameter(ndim=480000)
      dimension x(ndim)
      character typ*3,typ1*3,typ2*3,par*36,zeile*40
      logical trac
      common /tra/ trac
      t0s=0.d0
      hs=0.d0
      typ=zeile(1:3)
      par=zeile(5:40)
c  Subtract value of first sample before high-pass filtration
      if(typ(1:2).eq.'hp') then
        x0=x(1)
        do i=1,n
          x(i)=x(i)-x0
        enddo
      endif
c  decode filter type and read parameters
      if(typ.eq.'lp1') then
        read(par,*) t0
        it=11
      else if(typ.eq.'hp1') then
        read(par,*) t0
        it=12
      else if(typ.eq.'lp2') then
        read(par,*) t0,h
        it=21
      else if(typ.eq.'hp2') then
        read(par,*) t0,h
        it=22
      else if(typ.eq.'bp2') then
        read(par,*) t0,h
        it=24
      else if(typ.eq.'int') then
        t0=1.d0
        it= 1
      else if(typ.eq.'he1'.or.typ.eq.'le1') then
        read(par,*) t0s,t0
        it=13
      else if(typ.eq.'he2'.or.typ.eq.'le2') then
        read(par,*) t0s,hs,t0,h
        it=23
      else if(typ.eq.'sez') then
        read(par,*) t0s,hs
        t0=1.d12
        h=1.d0
        it=23
      else
      goto 2
      endif
c  calculate filter weights
      call rfk(it,t0/dt,h,t0s/dt,hs,f0,f1,f2,g1,g2)
      if(typ.eq.'le1') then
        fac=t0s/t0
      f0=f0*fac
      f1=f1*fac
      endif
      if(typ.eq.'le2') then
        fac=(t0s/t0)**2
      f0=f0*fac
      f1=f1*fac
      f2=f2*fac
      endif
c  perform recursive filtration
      if(typ(1:1).eq.'l') then
      call lrekfl(x,x,n,f0,f1,f2,g1,g2)
      else
      call rekfl(x,x,n,f0,f1,f2,g1,g2)
      endif
    8 format('# ',a3,2x,5f10.3)
c  confirm execution
      if(.not.trac) return
      if(it.eq.1) write(6,8) typ,t0
      if(it.eq.1) write(3,8) typ,t0
      if(it.eq.11.or.it.eq.12) write(6,8) typ,t0
      if(it.eq.11.or.it.eq.12) write(3,8) typ,t0
      if(it.eq.13) write(6,8) typ,t0s,t0
      if(it.eq.13) write(3,8) typ,t0s,t0
      if(it.eq.21.or.it.eq.22.or.it.eq.24) write(6,8) typ,t0,h
      if(it.eq.21.or.it.eq.22.or.it.eq.24) write(3,8) typ,t0,h
      if(it.eq.23) then
        if(typ.eq.'sez') write(6,8) typ,t0s,hs
        if(typ.eq.'sez') write(3,8) typ,t0s,hs
        if(typ.ne.'sez') write(6,8) typ,t0s,hs,t0,h
        if(typ.ne.'sez') write(3,8) typ,t0s,hs,t0,h
      endif
      return

c  Butterworth filters
    2 if(typ.eq.'lpb'.or.typ.eq.'hpb') then
        read(par,*) t0,m
        mm=m/2
        if(typ.eq.'lpb') then
          it1=11
          it2=21
          typ1='lp1'
          typ2='lp2'
        else if(typ.eq.'hpb') then
          it1=12
          it2=22
          typ1='hp1'
          typ2='hp2'
        else
          write(6,*) 'unknown filter type ',typ
          write(3,*) 'unknown filter type ',typ
          stop
        endif
        if(m.gt.2*mm) then
          call rfk(it1,t0/dt,h,t0s/dt,hs,f0,f1,f2,g1,g2)
      if(typ(1:1).eq.'l') then
      call lrekfl(x,x,n,f0,f1,f2,g1,g2)
      else
      call rekfl(x,x,n,f0,f1,f2,g1,g2)
      endif
          write(6,8) typ1,t0
          write(3,8) typ1,t0
        endif
        pih=2.d0*datan(1.d0)
        do 3 j=1,mm
          h=dsin(pih*(2*j-1)/m)
          call rfk(it2,t0/dt,h,t0s/dt,hs,f0,f1,f2,g1,g2)
      if(typ(1:1).eq.'l') then
      call lrekfl(x,x,n,f0,f1,f2,g1,g2)
      else
      call rekfl(x,x,n,f0,f1,f2,g1,g2)
      endif
          write(6,8) typ2,t0,h
          write(3,8) typ2,t0,h
    3   continue
      endif
      return
      end

      subroutine lrekfl(x,y,n,f0,f1,f2,g1,g2)
c  perform recursive low-pass filtering with transient suppression
      implicit double precision (a-h,o-z)
      parameter(ndim=480000)
      dimension x(ndim),y(ndim)
      xa=x(1)
      xaa=xa
      ya=xa
      yaa=xa
      do 1 j=1,n
      xn=dble(x(j))
      yn=f0*xn+f1*xa+f2*xaa+g1*ya+g2*yaa
      y(j)=yn
      xaa=xa
      xa=xn
      yaa=ya
    1 ya=yn
      return
      end

 
      subroutine rekfl(x,y,n,f0,f1,f2,g1,g2)
c  perform standard recursive filtering
      implicit double precision (a-h,o-z)
      parameter(ndim=480000)
      dimension x(ndim),y(ndim)
      xa=0.d0
      xaa=0.d0
      ya=0.d0
      yaa=0.d0
      do 1 j=1,n
      xn=dble(x(j))
      yn=f0*xn+f1*xa+f2*xaa+g1*ya+g2*yaa
      y(j)=yn
      xaa=xa
      xa=xn
      yaa=ya
    1 ya=yn
      return
      end
 
      subroutine rfk(it,t0,h,t0s,hs,f0,f1,f2,g1,g2)
c  determine coefficients for recursive filter
      implicit double precision (a-h,o-z)
      data zero,one,two,four,eight/0.d0,1.d0,2.d0,4.d0,8.d0/
      if(it.ne.1) goto 10
      f0=one/two/t0
      f1=f0
      f2=zero
      g1=one
      g2=zero
      return
   10 zpi=eight*datan(one)
      eps=zpi/t0
      f2=zero
      g2=zero
      if(it.gt.20) goto 20
      g1=(two-eps)/(two+eps)
      if(it.gt.11) goto 12
      f0=eps/(two+eps)
      f1=f0
      goto 14
   12 if(it.gt.12) goto 13
      f0=two/(two+eps)
      f1=-f0
      goto 14
   13 if(it.gt.13) return
      epss=zpi/t0s
      f0=(epss+two)/(eps+two)
      f1=(epss-two)/(eps+two)
   14 return
   20 epsq=eps*eps
      a=one-eps*h+epsq/four
      b=-two+epsq/two
      c=one+eps*h+epsq/four
      g1=-b/c
      g2=-a/c
      if(it.gt.21) goto 22
      f0=epsq/four/c
      f1=f0+f0
      f2=f0
      goto 25
   22 if(it.gt.22) goto 23
      f0=one/c
      f1=-f0-f0
      f2=f0
      goto 25
   23 if(it.gt.23) goto 24
      epss=zpi/t0s
      epssq=epss*epss
      as=one-epss*hs+epssq/four
      bs=-two+epssq/two
      cs=one+epss*hs+epssq/four
      f0=cs/c
      f1=bs/c
      f2=as/c
      goto 25
   24 if(it.gt.24) return
      f0=eps/two/c
      f1=zero
      f2=-f0
   25 return
      end
      
      
      subroutine average(x,n1,n)
c  remove average determined from first n1 samples
      implicit double precision(a-h,o-z)
      parameter(ndim=480000)
      dimension x(ndim)
      if(n1.lt.1.or.n1.gt.n) n1=n
      sum=0.
      do 1 j=1,n1
    1 sum=sum+x(j)
      sum=sum/n1
      do 2 j=1,n
    2 x(j)=x(j)-sum
c      write(6,'(a,i7,5x,f10.3)') '# avg     ',n1,sum
c      write(3,'(a,i7,5x,f10.3)') '# avg     ',n1,sum
      return
      end
 
      subroutine trend(typ,x,na,ne)
c  determine trend from samples na ..ne, 
c  remove it only when typ is not 'nix' Put coefficients into common.
c  parameters a and be refer to n1 as first point.
      implicit double precision(a-h,o-z)
      parameter(ndim=480000)
      dimension x(ndim)
      character typ*3
      logical trac
      common /tra/ trac
      common /tre/ a,b
      n2=ne-na+1
      gn = n2
      alpha = 0.5d0*gn*(gn+1.d0)
      beta = (2.d0*gn+1.d0)*(gn+1.d0)*gn/6.d0
      det = gn*beta-alpha*alpha
      sx = 0.d0
      sjx = 0.d0
      do 1001 j=1,n2
        sx = sx+x(j+na-1)
 1001   sjx = sjx+j*x(j+na-1)
      a = (sx*beta-sjx*alpha)/det
      b = (sjx*gn-sx*alpha)/det
      if(typ.eq.'nix') return
      do 1002 j=1,n2
 1002   x(j+na-1) = x(j+na-1)-a-b*j
      if(trac) write(6,'("# tre  ",2f10.3)') a,b
      if(trac) write(3,'("# tre  ",2f10.3)') a,b
      return
      end
  
      subroutine zspline(par,x,n)
c  remove z-trend with a z-spline
      implicit double precision (a-h,o-z)
      parameter(ndim=480000)
      dimension x(ndim)
      logical trac
      character par*36
      common /tre/ a,b
      common /tra/ trac      
      read(par,*) n1,n2,n3,n4
      fn32=n3-n2
      x2=x(n2)
      x3=x(n3)
      if(n1.lt.1.or.n2.lt.n1+1.or.n3.lt.n2+1.or.n4.lt.n3+1.or.n4.gt.n)
     & then
        write(6,*) n1,n2,n3,n4
        write(3,*) n1,n2,n3,n4
        write(6,'(a)') 'wrong parameters for zspline - stop'
        write(3,'(a)') 'wrong parameters for zspline - stop'
        stop
      endif
      call trend('nix',x,n1,n2)
      a1=a
      b1=b
      call trend('nix',x,n3,n4)
      a3=a
      b3=b
      do 2 j=n1,n2
    2 x(j) = x(j)-(a1+b1*(j-n1+1))
      do 3 j=n2+1,n3-1
    3 x(j) = x(j)-(x2+(j-n2)/fn32*(x3-x2))
      do 4 j=n3,n4
    4 x(j) = x(j)-(a3+b3*(j-n3+1))
      if(trac) write(6,'("# zspline eck ",4i10)') n1,n2,n3,n4
      if(trac) write(3,'("# zspline eck ",4i10)') n1,n2,n3,n4
      if(trac) write(6,'("# zspline coef",4f10.3)') a1,b1,a3,b3
      if(trac) write(3,'("# zspline coef",4f10.3)') a1,b1,a3,b3
      return
      end
 
 
      subroutine polytrend(npol,x,n,z)
c  remove polynomial trend
      implicit double precision (a-h,o-z)
      parameter(ndi=8,ndim=480000)
      dimension x(ndim),z(ndim),b(ndi),c(ndi,ndi)
      common /win/ a(ndi),fnh
      logical trac
      common /tra/ trac
      m=min(npol,ndi-1)
      fnh=n/2.d0
      one=1.d0
      if(trac) write(6,'("# pol  ",i5)') m
      if(trac) write(3,'("# pol  ",i5)') m
      do j=1, m+1
        do k=1, m+1
          c(j,k)=0
          do i=1, n
            if(z(i).gt.0.5) c(j,k)=c(j,k)+(dble(i)/fnh-one)**(j+k-2)
          enddo
        enddo
        b(j)=0
        do i=1,n
          if(z(i).gt.0.5) b(j)=b(j)+(dble(i)/fnh-one)**(j-1)*dble(x(i))
        enddo
      enddo
      call gauss(c,m+1,ndi,b,a)
      if(.false.) write(6,100) (j-1,a(j),j=1,m+1)
      if(.false.) write(3,100) (j-1,a(j),j=1,m+1)
  100 format(i5,e15.6)
      do i=1,n
        xpol=a(m+1)
        do j=m,1,-1
          xpol=xpol*(i/fnh-one)+a(j)
        enddo
        x(i)=x(i)-xpol
      enddo 
      return
      end
 
      subroutine gauss(aik,m,n,rs,f)
c  solve linear equations
      implicit double precision (a-h,o-z)
      common /tra/ trac
      logical trac
      dimension aik(n,n),rs(n),f(n),h(14),imax(13)
      if(.false.) then
        do i=1,m
          write(6,*) (aik(i,k),k=1,m), rs(i)
          write(3,*) (aik(i,k),k=1,m), rs(i)
        enddo
      endif
      do 1401 j=1,m
      aikmax=0.d0
      do 1402 k=1,m
      h(k)=aik(j,k)
      if(abs(h(k)).le.aikmax) go to 1402
      aikmax=abs(h(k))
      index=k
 1402 continue
      h(m+1)=rs(j)
      do 1403 k=1,m
      q=aik(k,index)/h(index)
      do 1404 l=1,m
 1404 aik(k,l)=aik(k,l)-q*h(l)
 1403 rs(k)=rs(k)-q*h(m+1)
      do 1405 k=1,m
 1405 aik(j,k)=h(k)
      rs(j)=h(m+1)
 1401 imax(j)=index
      do 1406 j=1,m
      index=imax(j)
 1406 f(index)=rs(j)/aik(j,index)
c     if(trac) write(6,*) 'gauss: ',m,' equations solved'
c     if(trac) write(3,*) 'gauss: ',m,' equations solved'
      return
      end
 
      subroutine heapsort(n,ra)
      implicit double precision (a-h,o-z)
      dimension ra(n)
	l=n/2+1
	ir=n
 10	continue
	if(l.gt.1) then
	  l=l-1
	  rra=ra(l)
	else
	  rra=ra(ir)
	  ra(ir)=ra(1)
	  ir=ir-1
	  if(ir.eq.1) then
	    ra(1)=rra
	    return
	  endif	
	endif
	i=l
	j=l+l
 20	if(j.le.ir) then
	  if(j.lt.ir) then
	    if(ra(j).lt.ra(j+1)) j=j+1
	  endif
	  if(rra.lt.ra(j)) then
	    ra(i)=ra(j)
	    i=j
	    j=j+j
	  else
       	  j=ir+1
	  endif
	  goto 20
	endif
	ra(i)=rra
	goto 10
 	end
 
 
      subroutine krum(x,n,ng,xk)
      implicit double precision (a-h,o-z)
      parameter(ndim=480000)
      dimension x(ndim),xk(ndim)
      double precision sj,sjq,sx,sjx,det,a,b,var
      nh = (ng - 1) / 2
      If (ng .lt. 3 .Or. ng .ne. 2 * nh + 1) Then
  	   Stop "ng ist keine ungerade Zahl >=3"
      EndIf
c Konstanten
      sj = ng * (ng + 1) / 2.
      sjq = ng * (ng + 1) * (2 * ng + 1) / 6.
c Erster Punkt des Krummheitsmasses
      sx = 0.
      sjx = 0.
      do j = 1, ng
        sx = sx + x(j)
        sjx = sjx + j * x(j)
      enddo
      det = ng * sjq - sj ** 2
      a = (sx * sjq - sjx * sj) / det
      b = (ng * sjx - sj * sx) / det
      var = 0.
      Do j = 1, ng
        var = var + (x(j) - a - j * b) ** 2
      enddo
      xk(1) = Sqrt(var / ng)
      Do j = 2, nh + 1
        xk(j) = xk(1)
      enddo
c Weitere Punkte rekursiv
      Do k = nh + 2, n - nh
        sjx = sjx - sx + ng * x(nh + k)
        sx = sx - x(k - nh - 1) + x(k + nh)
        a = (sx * sjq - sjx * sj) / det
        b = (ng * sjx - sj * sx) / det
        var = 0.
        Do j = k - nh, k + nh
          var = var + (x(j) - a - (j - k + nh + 1) * b) ** 2
        enddo
        xk(k) = Sqrt(var / ng + 0.12*(b*nh)**2)
c        xk(k) = Sqrt(var / ng)
      enddo
      Do k = n - nh + 1, n
        xk(k) = xk(n - nh)
      enddo
      return
      end
 
 
      subroutine reverse(x,n)
      implicit double precision (a-h,o-z)
      parameter(ndim=480000)
      dimension x(ndim)
      nh=n/2
      do i=1,nh
        xx=x(i)
        x(i)=x(n+1-i)
        x(n+1-i)=xx
      enddo
      return
      end

      subroutine mspline(par,x,n)
c  determine trends from samples n1 .. n2 and n3 .. n4,
c  and the step between the two trends. n1 .. n4 imported through par.
      implicit double precision(a-h,o-z)
      parameter(ndim=480000)
      dimension x(ndim)
      logical trac
      character typ*3,par*36
      common /tre/ a,b
      common /tra/ trac
      common /avstep/ avs,estep
      typ='nix'
      read(par,*) n1,n2,n3,n4
      fn32=n3-n2
      if(n1.lt.1.or.n2.lt.n1+1.or.n3.lt.n2+1.or.n4.lt.n3+1.or.n4.gt.n)
     & then
        write(6,*) n1,n2,n3,n4
        write(3,*) n1,n2,n3,n4
        write(6,'(a)') 'wrong parameters for mspline - stop'
        write(3,'(a)') 'wrong parameters for mspline - stop'
        stop
      endif
      call trend(typ,x,n1,n2)
      a1=a
      b1=b
      call trend(typ,x,n3,n4)
      a3=a
      b3=b
      if(trac) write(6,'("# mspline eck ",4i10)') n1,n2,n3,n4
      if(trac) write(3,'("# mspline eck ",4i10)') n1,n2,n3,n4
      if(trac) write(6,'("# mspline coef",4f10.3)') a1,b1,a3,b3
      if(trac) write(3,'("# mspline coef",4f10.3)') a1,b1,a3,b3
        fjj=(n2+n3)/2
        estep=(a3+b3*(fjj-n3+1))-(a1+b1*(fjj-n1+1))
      if(trac) write(6,*) 
     &   ' estimated step at sample ',int(fjj),' is ',estep
        return
      end

      subroutine steps(par,x,n)
c  determine average step
      implicit double precision(a-h,o-z)
      parameter(ndim=480000)
      dimension x(ndim)
      character par*36
      common /avstep/ avs,estep
      read(par,*) n1,n2,n3,n4
      av1=0.
      do i=n1,n2-1
      av1=av1+x(i)
      enddo
      av1=av1/(n2-n1)
      av2=0.
      do i=n3,n4-1
      av2=av2+x(i)
      enddo
      av2=av2/(n4-n3)
      avs=av2-av1
      jj=(n2+n3)/2
      write(6,*) '  measured step at sample ',jj,' is ',avs
c      write(6,*) '              step=',avs
c      write(3,*) '              step=',avs
      return
      end