c   .    1    .    2    .    3    .    4    .    5    .    6    .    7..
c
c     TILTCAL8
c     absolute calibration of seismometers by stepwise alternating tilt
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     Apr. 2008:  TILTCAL revised according to the new version 8 of DISPCAL.
c                 Low-pass filtering is now 2nd order forward and backward,
c                  with damping selected for shortest ringing. The filter-
c                  order parameter in the parameter file is ignored and may 
c                  be omitted. 
c                  The determination of the straightness limit was improved. 
c                 Double precision is used throughout.
c                 Minimum safety distance is 1.5 times filter corner period
c                 Free-format reading (ASL format) has been reprogrammed
c                    because my FORTRAN compiler does not handle
c                    free-format reading properly.
c     Aug. 2011:  zlim redefined: averge of median and 97%ile
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 tiltcal9
      use nas_system
      implicit double precision (a-h,o-z)
      parameter(ndi=8,ndim=480000)
      character text*72,nam1*50,zeile*40,nam2*50,nam3*50,nam4*50,par*36
      logical trac,pulselist,onegap,acc
      dimension x(ndim),y(ndim),xx(ndim),z(ndim)
      dimension sstep(199),step(199),ja(199),na(199),ne(199)
      common /win/ a(ndi),fnh
      common /avstep/ avs,estep
      common /tra/ trac

      write(6,*)
      write(6,*) ' This is program TILTCAL, 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 6 quiet seconds before'
      write(6,*) ' 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='tiltcal.par'
      endif
      open(11,file=nam3,status='old')
      read(11,*) text
      write(6,1) ' ',text
    1 format(a,a,a,a,a)
      read(11,*) nam1
      ipunkt=index(nam3,'.')
      nam4=nam3(1:ipunkt)//'out'
      open(3,file=nam4)
      write(3,1) ' This is program TILTCAL, version 8, April 2008'
      write(3,1)
      write(3,1) ' ',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,*) tilt
      perc=0.
      read(11,*,err=2019) gap,perc
 
 2019 read(11,*) tfilt
      read(11,*) pulselist
      read(11,*) corr
      close(11)
 
      fac=volts/tilt*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,
     &'   tilt=',tilt,'   gap=',gap
      write(3,2020) 'ts0=',ts0,'   hs0=',hs0,'   volts=',volts,
     &'   tilt=',tilt,'   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.gt.0.) then
        onegap=.true.
        ngap=1
      else
        gap=0.25
        onegap=.false.
        ngap=5
      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
      do i=1,n
        z(i)=1.
      enddo

      call average(x,nint(6./dt),n)
c  store a copy of the input signal for plotting
      nam2='tiltcal.dat'
      call outform(nam2,text,n,dt,tmin,tsec,x)
      
      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

c  deconvolution
      if(.not.acc) then
        write(6,*) 'conversion to apparent ground velocity'
        write(3,*) 'conversion to apparent ground velocity'
        write(zeile,2002) 'sez',ts0,hs0
 2002   format(a3,2x,2(f9.3,1x),15x)
        call filter(zeile,x,n,dt)
        write(6,*) 'removing polynomial trend'
        write(3,*) 'removing polynomial trend'
        call polytrend(npol,x,n,z)
      endif
      nam2='tiltcal.vel'
      call outform(nam2,text,n,dt,tmin,tsec,x)
      if(.not.acc) then
        write(6,*) 'conversion to apparent ground acceleration'
        write(3,*) 'conversion to apparent ground acceleration'
        call deriv(x,n,dt)
      endif
c  deconvolved signal will be stored later, after final trend removal

c  x is the acceleration signal, nearly a square wave
c  save x as xx so that different values of gap can be tried
  
      do i=1,n
        xx(i)=x(i)
      enddo
 
      write(6,*) 'measuring curvature of the acceleration trace'
      write(3,*) 'measuring curvature of the acceleration trace'

c Loop over five trial values of gap if gap = 0

      if(.not.onegap) then
      write(6,1)"======================================================"
      write(3,1)"======================================================"
      endif

      do 3000 ntry=1,ngap
      if(.not.onegap) trac=.false.
      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)
      ymin=1.e12
      do i=1,n
        ymin=min(ymin,y(i))
        y(i)=dlog10(max(y(i),1.d-6))
        z(i)=y(i)
      enddo
      if(ymin.le.1.d-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='tiltcal.str'
      if(onegap) call outform(nam2,text,n,dt,tmin,tsec,y)
c      nam2='tiltcal.srt'
c      if(onegap) call outform(nam2,text,n,dt,tmin,tsec,z)
 
      write(6,1010) 10.d0**zlim
      write(3,1010) 10.d0**zlim
 1010 format(' Threshold for quiet intervals set at',f12.3,' counts')
      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.
c      nam2='tiltcal.bol'
c      if(onegap) call outform(nam2,text,n,dt,tmin,tsec,z)


      if(.false.) then
      write(6,*)
      write(6,*) 'final trend removal on the acceleration signal:'
      mpol=npol-2
      if(mpol.lt.1) stop 'npol too small - stop'
      do i=1,n
        y(i)=x(i)
      enddo
      call deriv(y,n,dt)
      do i=1,n
        y(i)=y(i)*z(i)
      enddo
      call polytrend(mpol,y,n,z)
      write(6,1) '# re-integrating'
      xpolint=0.d0
      xavg=0.
      do i=1,n
        xpol=a(mpol+1)
        do j=mpol,1,-1
          xpol=xpol*(i/fnh-1.d0)+a(j)
        enddo
        xpolint=xpolint+xpol*dt
        x(i)=x(i)-xpolint
        xavg=xavg+x(i)
      enddo
      xavg=xavg/n
      do i=1,n
        x(i)=x(i)-xavg
      enddo      
      endif
      
      nam2='tiltcal.acc'
      call outform(nam2,text,n,dt,tmin,tsec,x)

      if(onegap) then
        write(6,*) 'locating and measuring steps of acceleration'
        write(3,*) 'locating and measuring steps of acceleration'
      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'
      do i=1,n
        y(i)=0.
      enddo
      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).gt.0.5)
         i=i+1
      enddo
      n2=i-1
c  this is the last sample of first quiet interval
      if(i.eq.n) stop 'no pulses found - check tiltcal.??? signals'
      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).gt.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

      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
      n1a=max(n1,n2-iavg)
      n4a=min(n4,n3+iavg)
      navg=min(n2-n1a,n4a-n3)
      n1a=n2-navg
      n4a=n3+navg
      do i=n1,n4
        y(i)=x(i)
      enddo 
      
      write(par,2003) n1a,n2,n3,n4a
 2003 format(4i9)
      call mspline(par,y,n)
c      write(par,2003) n1a,n2,n3,n4a
c      call steps(par,y,n) ! alte Methode, unnötig kompliziert
      n1=n3
      n2=n4
c      step(ns)=dabs(avs) ! alte Methode
      sstep(ns)=estep
      step(ns)=dabs(estep)
      
      if(onegap.and.pulselist) then
      write(6,'(a,i3,a,f16.6,a)') "Step # ",ns," has ",estep," counts"
      write(3,'(a,i3,a,f16.6,a)') "Step # ",ns," has ",estep," counts"
      endif
      
      if(i.lt.n-iavg.and.ns.lt.199) goto 1050
c  loop over steps ends here
 1080 continue
 
      write(6,*) '============ ',ns,' steps found.',' ============'
      write(3,*) '============ ',ns,' steps found.',' ============'
  
c  prepare residual for plotting
      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='tiltcal.res'
      call outform(nam2,text,n,dt,tmin,tsec,x)

c  Statistics
 
      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)
      if(onegap) then
      write(6,*)
      write(3,*)
      write(6,1163) ' Raw average step:  ',avgs,' +-',sigma,
     &' counts rms'
      write(3,1163) ' Raw average step:  ',avgs,' +-',sigma,
     &' counts rms'
 1163 format(a,f17.3,a,f14.3,a)
      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(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)*abs(step(k))
          enddo
          avgs=total/nrest
          total=0.
          do k=1,ns
            total=total+ja(k)*(abs(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.5.d-2*avgs.and.siga-sigma.lt.siga/nrest) exit
 1150 continue
      nrest=max(nrest,int((ns+1)/2))
      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')
      endif
      
 9000 continue
      
      if(onegap) then
      write(6,*)
      write(3,*)
            
c  produce time series of absolute steps
      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='tiltcal.stp'
      call outform(nam2,text,n,dt,tmin,tsec,xx) 
      write(3,*) 'Note: step variance in percent and limited to 1%'
      write(6,*) 'Note: step variance in percent and limited to 1%'
     
c  write a table of step sizes as file tilttable.lst
      open(9,file='tilttable.lst')
      pos=0.d0
      do k=1,ns        
        steptime=(na(k)+ne(k))/2.*dt
        stepmm=sstep(k)/avgs*tilt
        pos=pos+stepmm
        write(9,1180) k,steptime,stepmm,pos
      enddo
        write(9,1180) 0,0.,0.,0.
      write(9,1181)
      close(9)
 1180 format(i10,f10.1,2f10.4)
 1181 format('    step #    time s    mm/s^2   pos mm/s^2'/
     &'Note: The absolute step size may be inaccurate.'/
     &'Only the relative size should be evaluated.')      
      write(*,*) 'Steps are listed in file tilttable.lst'           
      
c  generate plot-parameter file
      if(acc) then
      open(4,file='winplot.par')
      write(4,*) '0, 5, 24, 16, 1,',n, ', 0.8'
      write(4,4) 'tiltcal.dat'
      write(4,4) 'tiltcal.acc'
      write(4,4) 'tiltcal.str'
c      write(4,4) 'tiltcal.srt'
c      write(4,4) 'tiltcal.bol'
      write(4,4) 'tiltcal.res'
      write(4,4) 'tiltcal.stp'
      write(4,4) 'Program TILTCAL v.8 (accel.)   Data File: '//nam1
    4 format(a)
      close(4)
      else
      open(4,file='winplot.par')
      write(4,*) '0, 6, 24, 16, 1,',n, ', 0.8'
      write(4,4) 'tiltcal.dat'
      write(4,4) 'tiltcal.vel'
      write(4,4) 'tiltcal.acc'
      write(4,4) 'tiltcal.str'
c      write(4,4) 'tiltcal.srt'
c      write(4,4) 'tiltcal.bol'
      write(4,4) 'tiltcal.res'
      write(4,4) 'tiltcal.stp'
      write(4,4) 'Program TILTCAL v.8    Data File: '//nam1
      close(4)
      endif
      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 tiltcal.srt and tiltcal.bol'
      write(6,*) 'results were saved to file ',trim(nam4)
      endif
      write(6,*)"======================================================"
      write(3,*)"======================================================"
      gap=gap*sqrt(2.)
 3000 continue
      stop
      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 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
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=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 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 deriv(x,n,dt)          
c  derivative (non-recursive, non-causal, symmetric-difference) 
      implicit double precision (a-h,o-z)
      parameter(ndim=480000)
      dimension x(ndim)            
      twodt=2.d0*dt      
      do 1 j=1,n-2    
    1 x(j)=(x(j+2)-x(j))/twodt          
      do 2 j=n-1,2,-1           
    2 x(j)=x(j-1)       
      x(n)=x(n-1)       
c      write(6,'("# dif - derivative by symmetric differences")')   
c      write(3,'("# dif - derivative by symmetric differences")')         
      return            
      end

      subroutine filter(zeile,x,n,dt)
      implicit double precision(a-h,o-z)
      parameter(ndim=480000)
      character typ*3,typ1*3,typ2*3,par*36,zeile*40
      logical trac
      dimension x(ndim)
      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=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=xa
      ya=xa
      yaa=xa
      do 1 j=1,n
      xn=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, remove it only when typ is 'tre'
c  parameters a and be refer to n1 as first point.
      implicit double precision(a-h,o-z)
      parameter(ndim=480000)
      dimension x(ndim)
      logical trac
      common /tra/ trac
      character typ*3
      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.'msp') 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 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

      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)+(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)+(i/fnh-one)**(j-1)*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 extrem(x,n)
      implicit double precision(a-h,o-z)
      parameter(ndim=480000)
      dimension x(ndim)
      xmin=1e24
      xmax=-xmin
      do 1 j=1,n
      if(x(j).lt.xmax) goto 3
      xmax=x(j)
      nmax=j
    3 if(x(j).gt.xmin) goto 1
      xmin=x(j)
      nmin=j
    1 continue
      write(6,2) xmin,nmin,xmax,nmax
      write(3,2) xmin,nmin,xmax,nmax
    2 format('# ext   xmin='f10.2' (j='i6')   xmax='f10.2' (j='i6')')
      return
      end

      subroutine rms(x,n)
      implicit double precision(a-h,o-z)
      parameter(ndim=480000)
      dimension x(ndim)
      xrms=0.
      do 1 j=1,n
    1 xrms=xrms+x(j)**2
      write(6,2) sqrt(xrms/n)
      write(3,2) sqrt(xrms/n)
    2 format('# rms  ='f10.3)
      return
      end

      subroutine gauss(aik,m,n,rs,f)
c  solve linear equations
      implicit double precision(a-h,o-z)
      logical trac
      common /tra/ 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)
      nh = (ng - 1) / 2
      If (ng .lt. 3 .Or. ng .ne. 2 * nh + 1) Then
  	   Stop " STOP: ng must be an odd number >=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)
      enddo
      Do k = n - nh + 1, n
        xk(k) = xk(n - nh)
      enddo
      return
      end