c   .    1    .    2    .    3    .    4    .    5    .    6    .    7..
c
c     DISPCALx3
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     May 2008:  Analysing three data sets with the location of pulses
c                determined from the first one (for measuring crosstalk)
c                 
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 step(99),sstep(99),ja(99),na(99),ne(99)
      character text*72,zeile*40,par*36,ff*50
      character*50 nam(3),nam0,nam1,nam2,nam3,nam4,namdat,namvel,namdis
      common /avstep/ avs,estep
      logical trac,pulselist,acc
      common /tra/ trac
      common /tre/ a,b

      write(6,*)
      write(6,*) '  This is program DISPCALx3, version May 2008'     
      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,*) '  This version analyses three-component data for' 
      write(6,*) '  crosstalk. The parameter file must contain two '
      write(6,*) '  additional lines with the file names of the second'
      write(6,*) '  and third component. Pulses are identified only'
      write(6,*) '  in the first component.'
      write(6,*)
      write(6,*) '=== There MUST BE AT LEAST 6 quiet seconds before the
     & first pulse! ==='
      write(6,*)
    1 format(a,a,a,a,a)

      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='dispcalx3.par'
      endif
      open(11,file=nam3,status='old')
      read(11,*) text
      write(6,'(1x,a)') text
      read(11,*) nam(1)
      read(11,*) nam(2)
      read(11,*) nam(3)
      ipunkt=index(nam3,'.')
      nam4=nam3(1:ipunkt)//'out'
      open(3,file=nam4)
      write(3,1) ' This is program DISPCAL, version 8, April 2008'
      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.
      read(11,*) hs0
      read(11,*) volts
      read(11,*) displ
      read(11,*) 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)
      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.lt.1.5*tfilt) then
      gap=1.5*tfilt
      write(6,2020) 'safety distance ''gap'' was increased to ',gap
      write(3,2020) 'safety distance ''gap'' was increased to ',gap
      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)
 

c *************** Loop over three components ****************

      do 3333 iseis=1,3
      nam1=trim(nam(iseis))
      in=index(nam1,'.')
      if(in.eq.0) then
      write(*,*) ' sorry, the data file ',nam1
      write(*,*) ' must have a three-letter extension such as *.dat'
      stop
      endif
      nam0=nam1(1:in)
      namdat=trim(nam0)//'dat'
      namvel=trim(nam0)//'vel'
      namdis=trim(nam0)//'dis'
      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
      
      call average(x,nint(6./dt),n)     
c  store a copy of the input signal for plotting
      call outform(namdat,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)
        call filter(zeile,x,n,dt)
        call reverse(x,n)
        call filter(zeile,x,n,dt)
        call reverse(x,n)
      endif

c  deconvolution
      if(.not.acc) then
        write(zeile,2002) 'sez',ts0,hs0
 2002   format(a3,2x,2(f9.3,1x),15x)
      else
        write(zeile,2005) 'int'
 2005   format(a3,17x)
      endif  
      call filter(zeile,x,n,dt)
      do i=1,n
        xx(i)=1.
      enddo
      call polytrend(npol,x,n,xx)

c  store deconvolved signal for visual check
c     nam2='dispcal.v_1'
c     call outform(nam2,text,n,dt,tmin,tsec,x)
 
c   x is now the ground-velocity signal, will be detrended and integrated
c   to give a displacement signal representing the steps
 
 
      if(iseis.eq.1) then                                               ! if iseis=1
 
      avglen=5.*gap
      iab=nint(gap/dt)
      iavg=nint(avglen/dt)
      do i=1,n
        z(i)=1.
        y(i)=0.
      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)
      imax=0
      zsl=-1e12
      ndr=n/3
      nf=n/512
      nn=n-nf-3
      do i=ndr,nn
        if(z(i+nf)-z(i-nf).gt.zsl) then
          imax=i
          zsl=z(i+nf)-z(i-nf)
        endif
      enddo
      zlim=z(imax)  
      do i=1,n
        y(i)=y(i)-zlim
        z(i)=z(i)-zlim
      enddo
 
c     store quietness signals for visual check
c      nam2='dispcal.str'
c      call outform(nam2,text,n,dt,tmin,tsec,y)
c      nam2='dispcal.srt'
c      call outform(nam2,text,n,dt,tmin,tsec,z)
 
      if(iseis.eq.1) then
      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')
      write(6,1010) 10.**zlim
      write(3,1010) 10.**zlim
 1010 format(' Threshold for quiet intervals set at',f12.3,' counts')
      endif

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.
c      nam2='dispcal.bol'
c      call outform(nam2,text,n,dt,tmin,tsec,z)

      endif                                                             ! endif iseis=1

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
      call polytrend(npol,x,n,z)
c  store restored velocity signal and residual noise for visual check
c       call outform(namvel,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)
c      nam2='dispcal.w_2'
c      call outform(nam2,text,n,dt,tmin,tsec,y) 

c  finding the steps. Before looping, determine first quiet interval.
      if(pulselist) write(6,'(a)') 'Pulses: times in seconds'
      if(pulselist) write(3,'(a)') 'Pulses: times in seconds'

      if(iseis.eq.1) then                                               ! if iseis=1
      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(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,i2,3(a,f7.2))
      endif
      na(ns)=n2
      ne(ns)=n3
      n1=n3
      n2=n4
      if(i.lt.n-iavg.and.ns.lt.99) goto 1050
      endif                                                             ! endif iseis=1
c  loop over steps ends here
 1080 continue

      write(6,*) '============ ',ns,' pulses found.',' ============'
      write(3,*) '============ ',ns,' pulses found.',' ============'
      if(ns.eq.0) goto 9000

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
c      nam2='dispcal.vtr'
c      call outform(nam2,text,n,dt,tmin,tsec,y)
      do j=1,n
        x(j)=x(j)-y(j)
      enddo

      call outform(namvel,text,n,dt,tmin,tsec,x)
      
      write(zeile,2000) 'int'
 2000 format(a3,37x)
      call filter(zeile,x,n,dt)
      call average(x,n,n)
      
      call outform(namdis,text,n,dt,tmin,tsec,x)           
            
c  now, finally, we measure the steps. n1 .. n4 are redefined symmetric
c  to pulse k.
 
        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
        if(iseis.eq.1) then
          step(k)=dabs(estep)
          sstep(k)=sign(1.d0,estep)
        else
          step(k)=estep*sstep(k)
        endif
        if(pulselist) then
        write(6,'(a,i2,a,f16.6,a)') "Step # ",k," has ",estep," counts"
        write(3,'(a,i2,a,f16.6,a)') "Step # ",k," has ",estep," counts"
        endif
      enddo
      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
c      nam2='dispcal.res'
c      call outform(nam2,text,n,dt,tmin,tsec,x)
 
  
c  Statistics
 
      do i=1,99
        ja(i)=1
      enddo
      total=0.
      do k=1,ns
        total=total+ja(k)*step(k)
      enddo
      avgs=total/ns
      total=0.
      do k=1,ns
        total=total+ja(k)*(step(k)-avgs)**2
      enddo
      sigma=sqrt(total/ns)
      write(6,*)
      write(3,*)
      if(iseis.eq.1) then
      write(6,1163) ' Raw average abs. step:  ',avgs,' +-',sigma,
     &' counts'
      write(3,1163) ' Raw average abs. step:  ',avgs,' +-',sigma,
     &' counts'
 1163 format(a,f12.3,a,f14.3,a)
      else
        write(6,1163) ' Raw average signed step:  ',avgs,' +-',sigma,
     &' counts'
        write(3,1163) ' Raw average signed step:  ',avgs,' +-',sigma,
     &' counts'
      endif
      if(acc) then
        write(6,1161) avgs*fac,sigma*fac
        write(3,1161) avgs*fac,sigma*fac
 1161 format(' Raw generator constant: ',f10.4,' +-',f9.4,' Vs^2/m')
      else
        write(6,1159) avgs*fac,sigma*fac
        write(3,1159) avgs*fac,sigma*fac
      endif
 1159 format( ' Raw generator constant: ',f10.2,' +-',f9.2,' Vs/m')

      if (ns.gt.2) then
        write(6,*)
        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
          write(6,1172) kmax,nrest,avgs*fac,sigma*fac
          write(3,1172) kmax,nrest,avgs*fac,sigma*fac
        else
          write(6,1170) kmax,nrest,avgs*fac,sigma*fac
          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(iseis.eq.1) then
        write(6,*)
        write(3,*)
        if(acc) then
          write(6,1173) nrest,avgs*fac,sigma*fac
          write(3,1173) nrest,avgs*fac,sigma*fac
        else
          write(6,1171) nrest,avgs*fac,sigma*fac
          write(3,1171) nrest,avgs*fac,sigma*fac
        endif
 1173   format('                          Best ',i3,' Steps: ',f10.4,
     &  '  +-',f9.4,' Vs^2/m')
 1171   format('                          Best ',i3,' Steps: ',f10.2,
     &  '  +-',f9.2,' Vs/m')
      else
        write(6,1190) 
        write(3,1190)
 1190 format(/' Note: In the absence of real outliers, the raw average m
     &ay be the best'/' estimate of the crosstalk amplitude. Use outlier
     & elimination only to'/' assess the stability of the result.')
      endif
      endif
 9000 continue      
           
      write(6,*)
      write(3,*)
c  generate plot-parameter file
      if(iseis.eq.1) then                                               ! if iseis=1
      open(4,file='winplot.par')
        write(4,*) '0, 9, 24, 16, 1,',n, ', 0.8'
        ff=namdat
        write(4,4) namdat
        write(4,4) namvel
        write(4,4) namdis
      else                                                              ! else
        write(4,4) namdat
        write(4,4) namvel
        write(4,4) namdis
      endif                                                             ! endif iseis=1
    4 format(a)
      write(6,*)"======================================================"
      write(3,*)"======================================================"
 3333 continue
      write(4,4) 'Program DISPCAL v.8    first data file: '//trim(ff)
      
      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)
       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   schön wär's, geht aber nicht zuverlässig - Bug im Compiler?
        nz=1
        nval=0
        do
          read(7,'(a)',end=31) zeile
          read(zeile,'(120a)') zarr
          nn=0
          do i=1,119
            if(zarr(i).ne.b.and.zarr(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)
      parameter(ndim=480000)
      dimension ra(ndim)
	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