c     program SEIFE
c
c  REVISIONS (there were many; this protocol starts November 2012)
c  Nov. 2012: the global 'frq' command permits specification of corner  
c             frequencies in place of corner periods
c             prewarping activated for butterworth filters
c  July 2013: "nul" parameters limited to 1..n
c  Sept. 2013: intermediate results can be stored with 'sav'
c  Oct. 2014: output format 'asc' implemented (MATLAB compatible)
c
c   time-domain signal processing (filtering, windowing, decimation, etc.)
c   Author: E. Wielandt
c   The original version was written in 1984 at ETH Zurich
c
c   convention for the 'seife' type file format:
c   - one header line, up to 72 characters
c   - up to 48 lines of signal description and processing log
c     these lines have '%' as the first character
c     the processing log is updated by the present program
c   - one line specifying the number of samples, the format, the sampling
c     interval, and the start time in min and sec after midnight. this
c     line may be repeated after the data when its contents is not known
c     at the time when the file header is written; however the format
c     must always be specified.
c   - data in any FORTRAN format or in free ASCII format

c   The following INPUT formats are also supported:    
c   - ASCII with a variable number of samples per line, header lines ignored (asc)
c   - single-column ASCII, header lines ignored (sca)
c   - MARS88 - ASCII with header line
c   - ASL (Peterson) ASCII format such as produced by Quanterra's CIMARRON
c   - BDF ASCII as used in the Albuquerque Seismological Lab
c   - PSN textfiles from WINQUAKE (Larry Cochrane)
c
c   OUTPUT can be either in SEIFE format or in one-column or two-column headerless
c   ASCII format compatible with GNUPLOT, MATLAB and LABVIEW
c
c   The desired time-domain processes are defined in a command file whose
c   name is entered either interactively or as a runstring. this file also
c   specifies the names of the input and output files (up to 24 pairs; all
c   processes apply to all input files). each line of the command file
c   defines one processing step (a time window, filter, etc.) and 
c   consists of a thee-letter code followed by one or more blanks and some
c   numerical parameters (such as corner period, damping, window limits).
c   parameters must be separated by blanks or commas and fit into
c   columns 5 to 60. The following codes are accepted:
c
c   blank first character, or rem: line is ignored     
c
c   commands for input and output:
c
c   Reading data files:
c   fil  nam1 nam2: read files in SEIFE format. nam1=input file, nam2=output file
c        nam2='_join_' will append the next input file to the present one
c   mar  nam1 nam2:  read MARS88-ASCII files from MARSDUMP
c   asl  nam1 nam2:  read ASL ascii format (from Cimarron)
c   bdf  nam1 nam2:  read BDF ascii (headers interspersed with data - caution, 
c                    anything that is not a number will be ignored)
c   psn  nam1 nam2:  read PSN text files from WINQUAKE (Larry Cochrane)
c   asc  nam1 nam2 dt lin:  read ASCII files where dt is the sampling interval and 
c        lin the number of header lines to be skipped. (Header lines are echoed
c        but not interpreted.) Will read almost any ASCII data.
c        The last character of the last line must be a <cr>, or the last
c        sample will be skipped. asc is slow; use sca if possible.
c   sca  (single-column ascii) like asc but reads only one sample per line (much 
c        faster than asc)
c
c  global commands:
c
c   frq  specify filter corners as frequencies (default is periods)
c
c   chp  nch: chop signal into sections of nch samples for output.
c        Append .00, .01 etc. to name of output file (maximum 24 files).
c
c   com  insert comment into file header
c   rem  insert comment into file header
c
c   lim  n: limit the length of the input series to n samples
c
c   ncl  no comment lines
c
c   npl  don't write a plot-parameter file
c
c   pfi  set a path for the input files. Don't use quotation marks.
c        The path may also be entered as part of the file name.
c
c   pfo  set a path for the output files. Don't use quotation marks.
c        The path may also be entered as part of the file name.
c
c   pli  set up winplot.par to plot input and output files. Works for
c        SEIFE format only, and for up to six input signals. Otherwise ignored.
c        Do not apply a time window together with pli (traces will not be
c        synchronous if you do).
c
c   plt w h:  set nonstandard plot width w and height h in winplot.par
c
c   plw n1 n2: set window for plot, rms, and ext. Set 'remove average' for plot.
c
c   cfi  specify name for the Labview cfg file. Default is labview.cfg
c
c   sav  save intermediate result; attach -1, -2 etc. to name of output file
c
c   end  [format]: write output file. 
c        The output file may be identical to the input file.
c     Legal format specifications are: 
c       (any fortran format) - specified format, in parentheses but without
c          quotation marks, must fit into 20 characters. Example:(5x,10f8.3)
c       fre or <blank< - fortran format (5f13.x) with x selected auto-
c          matically.
c     SEIFE can re-read the above formats automatically.
c       asc - pure data f13.x without header, with x selected automatically.
c       lab - two-column output (time and data) for labview or gnuplot
c     To make these two formats re-readable by SEIFE, header lines according
c       to the SEIFE format specifications must be added.
c
c   This version produces a plot-parameter file 'winplot.par' for the new
c   version of winplot (with from - to indices) if the output is in SEIFE
c   format, or a configuration file for plotting with LABVIEW if the output
c   is in LABVIEW format.
c
c   Commands for signal processing:
c
c   abs  absolute value
c   add  a: add a to each sample
c   avg  n: remove average (determined from first n samples)
c   chp  nch: chop signal into sections of nch samples for output.
c        Append .00, .01 etc. to name of output file (maximum 24 files).
c   clp  clp: clip signal at level clp
c   cos  n: cosine-taper over n points at both ends
c   csi  dt: change sampling interval to dt by linear interpolation
c   dec  n: decimate to every n'th sample
c   del  n: delay by n samples (opposite to skp)
c   dif  differentiate the signal
c   dis  n a: add nonlineqar distortion: x <== x + a*xmax*(x/xmax)**n
c   exp  n: remove exponential trend determined from first n samples
c   ext  determine extremal samples (signal is unchanged)
c   fac  f: multiply each sample by f
c   int  integrate the signal
c   lin  remove linear trend through first and last sample
c   nor  a: normalize amplitude to a
c   nop  a: normalize peak amplitude to a
c   nox  a: normalize extrema to +-a
c   npw  no prewarp for recursive filters
c   nul  n1  n2: nulls the signal between samples n1 and n1 (inclusive)
c   pad  n; pad with zeroes to length n
c   pol  n: remove polynomial trend of degree n
c   qnt  a: quantization to multiples of the quantization interval a; 
c           rounding to integers if a=1
c   rev  reverse signal in time (useful for backward filtering)
c   rms  determine rms amplitude (signal is unchanged)
c   set  n1 n2 a: set samples n1 through n2 to value a
c   sez  t0 h: remove seismo response with free period t0 and damping h
c   shl  n: shift signal left in time by n samples (pad with zeros at end)
c   shr  n: shift signal right in time by n samples (pad with zeros at start)
c   sig  signum function
c   sin  n1  n2: window from n1 to n2 with sine weight
c   sis  n1  n2: the same, with a squared-sine weight
c   skp  n: skip n samples
c   spl  clp,n1,n2,apex:  interpolate clipped samples (value >= clp) with
c        cubic spline from sample n1 to n2. clp=0 is replaced by a default
c        value of slightly less than 2**18. if apex is specified other than
c        zero, interplation is done with a parabola with specified height of
c        apex.
c   sqr  square each sample (used for the detection of quadratic distortions)
c   srt  heapsort
c   swi  sec1 sec2: window defined by seconds after first sample
c   zsp  n1 n2 n3 n4  detrend with a z-spline with corners n1 .. n4
c   tap  n1: signal is tapered from the first sample to sample n1 and 
c        from sample (last-n1) to last. The weight is a cosine half-cycle.
c        The amplitude is normalized so that the energy of a stationary 
c        signal is unchanged.
c   tid  [n]: remove tides. tides are interpolated over n samples.
c        default (n=0 or blank) is interpolation over 5 minutes.
c   tr3  n1 t0 h: fit and remove remove transient with eigenperiod t0 
c                 and damping h, after sample n1
c   tr6  n1 t0 h f1 f2 f3: remove transient using amplitudes determined 
c        by tr3 (possibly from a different signal)
c   tre  n: remove linear trend (determined from first n samples)
c   twi  tmin1  tsec1  tmin2  tsec2: window defined by time after midnight
c   tim  tmin  tsec  set time stamp
c   win  n1  n2: time window from sample n1 to n2 (n1 becomes first sample)
c
c   recursive (IIR) filters:
c
c   filters are specified by a three-character code for the type
c                            and 1 to 4 parameters.
c   type: one of lp1, hp1, le1, he1, lp2, hp2, bp2, le2, he2, lpb, hpb, sez
c       lp = low pass, hp= high pass, bp = band pass
c             parameters: corner period [,damping if second-order]
c       le = inverse low pass, he = inverse high pass (equalizers)
c             parameters: old corner period, [old damping], new corner
c             period, [new damping]. These are APPROXIMATE inverses.
c         1 or 2: number of poles, i.e. order of denominator polynomial
c         b:      Butterworth, specify period and order
c     Example: lpb 30. 6 is a 30 sec, sixth-order Butterworth low-pass filter.
c       sez t0 h is the exact inverse to hp2 t0 h
c end of seife info.


      program seife
c the following statement is compiler-specific, it makes runstring
c parameters available
      use nas_system
      parameter(ndim=200000000)
      dimension x(ndim)
      dimension fdt(32),lin(32)
      real*8 ddt,tcom1,tcom2,tcom1m,tcom1s,tcom2m,tcom2s
      character typ*3,text*72,ttext*72,zeile*99,par*56,info*78
      character*60 prolog(48)
      character*48 nam1(32),nam2(32),nam3,nam4,cfgname,inpath,outpath
      character*256 cfgline(32)
      character timestring*12,einsbisneun*2
      character*3 datform(32)
      logical nil,join,bin,bina,ncl,npl,npw,pli,frq
      data fdt/32*0./, lin/32*0/, breit/24./, hoch/16./
      data cfgname,inpath,outpath /3*'#'/
      data ncl,npl,npw,pli,frq /5*.false./

      write(*,*)
      write(*,*) 'Latest revision: Feb. 2013. For a list of options,'
      write(*,*) 'copy the source code into the working directory and'
      write(*,*) 'type "seife info". See also the program description',
     &           ' seife.doc'
      write(*,*) 

      if(iargc().eq.0) then
        nam3='seife.par'
      else
        call getarg(1,nam3)
      endif
      if(nam3(1:4).eq.'info') then
        open(7,file='seife.f',status='old')
        do i=1,199
          read(7,'(a)',end=178) info
          write(6,'(a)') info
          if (info(1:8).eq.'c end of') stop
        enddo
  178   stop
      endif
      write(6,'(" reading commands from file ",a)') trim(nam3)
      open(7,file=nam3,status='old')
      nchin=ndim
      nmax=ndim
      ntrend=0
      call files(nam1,nam2,fdt,lin,nseis,zeile,nmax,nchin,ncl,npl,npw,
     & pli,frq,cfgname,inpath,outpath,datform)
      close(7)
      tcom1=-1.d12
      tcom2=1.d12
      lines=0
      join=.false.
      npl1=1
      npl2=nmax
      
      
      do 1 ns=1,nseis
      nsav=0
      nch=nchin
      npro=0
      if(join) then
        n1=n+1
      else
        n1=1
      endif
      bina=bin
      bin=.false.
      dt=fdt(ns)
      if(join.and.(bina.ne.bin)) then
        stop ' Cannot join binary and ascii files'
      endif
      n=ndim+1-n1
      if(bin) then
       stop ' binary input not implemented'
      else if (lin(ns).ge.0) then
       n=min(n,nmax)
       call inpasc(nam1(ns),lin(ns),ttext,prolog,npro,n,dt,ttmin,ttsec,
     &   x(n1),ndim,datform(ns))
      else
       stop ' unknown format specification'
      endif
      npl2=min(npl2,n)
      tcom1=max(tcom1,60.d0*ttmin+ttsec)
      tcom2=min(tcom2,60.d0*ttmin+ttsec+(n-1)*dt)
      if(n1.eq.1) then
        tmin=ttmin
        tsec=ttsec
        text=ttext
      endif
      if(n.eq.0.or.dt.eq.0.) then
      stop ' n or dt unspecified - check data header'
      endif
      if(join) then
      n=n1+n-1
      write(6,11) n
   11 format(' data appended to previous data, total ',i8,' samples')
      endif
      join=nam2(ns).eq.'_join_'
      if(join) goto 1

      npro=npro+1
      write(prolog(npro),5) trim(nam1(ns)),trim(nam2(ns))
    5 format('from ',a,'  to ',a)
      open(7,file=nam3,status='old')
    2 continue
      read(7,'(a)',end=4) zeile
      if(zeile(1:1).eq.' ') goto 2
      typ=zeile(1:3)
      if(zeile(4:4).ne.' ') write(6,*) 'WARNING: parameter file has non-
     &blank entry in column 4'
      par=zeile(5:60)
      if(typ.eq.'fil'.or.typ.eq.'asc'.or.typ.eq.'lim'.or.typ.eq.'chp'
     & .or.typ.eq.'asl'.or.typ.eq.'bdf'.or.typ.eq.'ncl'.or.typ.eq.'mar'
     & .or.typ.eq.'pfi'.or.typ.eq.'pfo'.or.typ.eq.'cfi'.or.typ.eq.'npl'
     & .or.typ.eq.'npw'.or.typ.eq.'pli'.or.typ.eq.'sca'.or.typ.eq.'psn'
     & .or.typ.eq.'frq')
     &  goto 2
      npro=npro+1
      if(npro.gt.48) then
      stop ' prolog has more than 48 lines'
      endif
      write(prolog(npro),'(a)') zeile
      if(typ.eq.'ext'.or.typ.eq.'rms') npro=npro-1
      if(typ.eq.'rem'.or.typ.eq.'com') goto 2
      nil=.true.
      n=min(n,nmax)
      npl2=min(npl2,n)
      ddt=dt
      if(typ.eq.'sqr') call square (nil,x,n)
      if(typ.eq.'fac') call factor(nil,par,x,n)
      if(typ.eq.'nor'.or.typ.eq.'nop'.or.typ.eq.'nox') 
     &                                    call norm(nil,typ,par,x,n)
      if(typ.eq.'add') call add(nil,par,x,n)
      if(typ.eq.'dif') call deriv(nil,par,x,n,dt)
      if(typ.eq.'dis') call distort(nil,par,x,n)
      if(typ.eq.'spl') call interp(nil,par,x,n)
      if(typ.eq.'shl'.or.typ.eq.'shr') call shift(nil,typ,par,x,n)
      if(typ.eq.'csi') call intpo (nil,par,dt,x,n)
      if(typ.eq.'win'.or.typ.eq.'sin'.or.typ.eq.'skp'.or.typ.eq.'tap'
     & .or.typ.eq.'sis'.or.typ.eq.'cos'.or.typ.eq.'swi') 
     & call window(nil,pli,typ,par,x,n,dt,tmin,tsec)
      if(typ.eq.'twi') call timwin(nil,par,x,n,dt,tmin,tsec)
      if(typ.eq.'tim') call timeset(nil,par,tmin,tsec)
      if(typ.eq.'nul') call zero(nil,par,x,n)
      if(typ.eq.'set') call set(nil,par,x,n)
      if(typ.eq.'clp') call clip(nil,par,x,n)
      if(typ.eq.'dec') call decim(nil,par,x,n,dt)
      if(typ.eq.'avg') call mean(nil,par,x,n)
      if(typ.eq.'tre'.or.typ.eq.'lin') call trend(nil,typ,par,x,n)
      if(typ.eq.'zsp') call zspline(nil,typ,par,x,n)
      if(typ.eq.'pol') call polytrend(nil,par,x,n)
      if(typ.eq.'exp') call expon(nil,par,x,n)
      if(typ.eq.'rev') call reverse(nil,par,x,n)
      if(typ.eq.'tid') call tides(nil,par,x,n,ddt)
      if(typ.eq.'tr3'.or.typ.eq.'tr6') call trans(nil,typ,par,x,n,ddt)
      if(typ.eq.'ext') call extrem(nil,x,npl1,npl2,n)
      if(typ.eq.'rms') call rms(nil,x,npl1,npl2,n)
      if(typ.eq.'pad') call pad(nil,par,x,n,ndim)
      if(typ.eq.'del') call delay(nil,par,x,n)
      if(typ.eq.'srt') call heapsort(nil,x,n)
      if(typ.eq.'sig'.or.typ.eq.'abs') call sig(nil,zeile,x,n)
      if(typ.eq.'plt') call plotpar(nil,par,breit,hoch)
      if(typ.eq.'qnt') call quant(nil,par,x,n)
      if(typ.eq.'plw') call plotwindow(nil,par,npl1,npl2,n,ntrend)
      if(nil) call filter(nil,zeile,x,n,ddt,npw,frq)
c  all done - store result
      if(typ.eq.'sav'.or.typ.eq.'end') then
        nch=min(max(nch,int((n+9)/24)),n)
        nfile=int(n-1)/nch+1
        if(typ.eq.'sav'.and.nfile.gt.1) then
          write(6,*) 'sav ignored, it is incompatible with chp'
          goto 2
        endif
        if(typ.eq.'sav'.and.nseis.gt.1) then
          write(6,*) 'sav ignored, it accepts only one input file'
          goto 2
        endif
        do nfi=1,nfile
          n1=1+(nfi-1)*nch
          n2=min(nch,n-n1+1)
          nam4=nam2(ns)
          if(nfile.gt.1) then
            i=index(nam4,' ')
            write(nam4(i:i+2),6) nfi
    6       format('.',i2.2)
            tsec=60.*tmin+tsec+(n1-1)*dt
            tmin=int(tsec/60.)
            tsec=tsec-60.*tmin
          endif
          if(typ.eq.'sav') then 
            nsav=nsav+1
            if(nsav.gt.9) goto 2
            write(einsbisneun,'("-",i1)') nsav
            nam4=trim(nam2(ns))//einsbisneun
          endif
          call outform(nam4,par,text,prolog,npro,ncl,n2,dt,tmin,tsec,
     &        x(n1),cfgline,lines)
        enddo
        if(typ.eq.'sav') goto 2
        goto 1
      endif
c  produce error message if code unrecognized
      if(nil) then
        write(*,*) ' undefined procedure ',typ
        stop
      else
        goto 2
      endif
    1 continue
    
    
c print common start and end time
      tcom1m=int(tcom1/60.d0)
      tcom1s=tcom1-60.d0*tcom1m
      tcom2m=int(tcom2/60.d0)
      tcom2s=tcom2-60.d0*tcom2m
      write(6,20) tcom1m,tcom1s,tcom2m,tcom2s
   20 format(" common time interval: from ",f6.0," min ",f7.3," sec"/
     &       "                         to ",f6.0," min ",f7.3," sec")

      if(lines.eq.0) then
c   produce plot-parameter file
        open(9,file='winplot.par')        
        if(.not.npl.and.nsav.eq.0) then
          jcount=0
          do ns=1,ns
            if(trim(nam2(ns)).eq.'_join_') jcount=jcount+1
          enddo
          if(pli.and.jcount.eq.0) then
            write(9,7) ntrend,2*nseis,breit,hoch,npl1,npl2
            do ns=1,nseis
              write(9,'(a)') trim(nam1(ns))
              write(9,'(a)') trim(nam2(ns))
            enddo
          else
            write(9,7) ntrend,nseis-jcount,breit,hoch,npl1,npl2
    7       format(i3,',',i3,",",f5.1,",",f5.1,",",i8,",",i8,", 0.8")
            do ns=1,nseis
            if(trim(nam2(ns)).ne.'_join_') write(9,'(a)') trim(nam2(ns))
            enddo
          endif
        else if(.not.npl.and.nsav.gt.0) then
          nsav=min(nsav,9)
          
	  if(pli) then
	    write(9,7) ntrend,nsav+2,breit,hoch,npl1,npl2
	    write(9,'(a)') trim(nam1(1))
	    do i=1,nsav
	    write(einsbisneun,'("-",i1)') i
	      write(9,'(a)') trim(nam2(1))//einsbisneun
	    enddo
	    write(9,'(a)') trim(nam2(1))
	  else
	    write(9,7) ntrend,nsav+1,breit,hoch,npl1,npl2
	    do i=1,nsav
	      write(einsbisneun,'("-",i1)') i
	      write(9,'(a)') trim(nam2(1))//einsbisneun
	    enddo
	    write(9,'(a)') trim(nam2(1))
	  endif  
	  
	  close(9)
	  write(6,*) ' Plot-parameter file winplot.par was generated'
        else
          write(6,*) ' this should not happen - check program'
          stop
        endif
      else
c   produce cfg file for Labview
        if(index(cfgname,'#').ne.0) then
          cfgname='labview.cfg'
        endif
        if(index(outpath,'#').eq.0) then
          cfgname=trim(outpath)//'\'//cfgname
        endif
        cfgline(lines+1)='Type: Single'
        ithr1=int((tcom1m*60.d0+tcom1s)/3600.d0)
        itmin1=int((tcom1m*60.d0+tcom1s-3600.d0*ithr1)/60.d0)
        tcom1s=int((tcom1m*60.d0+tcom1s-3600.d0*ithr1-60.d0*itmin1))
        write(timestring,77)ithr1,itmin1,tcom1s
   77 format(i2,':',i2,':',f6.3)
        cfgline(lines+2)='Startzeit: 00/00/00 '//timestring
        ithr2=int((tcom2m*60.d0+tcom2s)/3600.d0)
        itmin2=int((tcom2m*60.d0+tcom2s-3600.d0*ithr2)/60.d0)
        itsec2=int((tcom2m*60.d0+tcom2s-3600.d0*ithr2-60.d0*itmin2))
        write(timestring,77)ithr2,itmin2,tcom2s
        cfgline(lines+3)='Endzeit:   00/00/00 '//timestring
        open(9,file=cfgname)
        do i=1,lines+3
          write(9,'(a)') trim(cfgline(i))
        enddo
        close(9)
        write(6,*) ' configuration file ',trim(cfgname),' was generated'
      endif
      stop

    4 write(*,*) ' end line missing in file ',trim(nam3)
      stop
      end
c
      subroutine inpasc(name,lin,text,prolog,npro,n,dt,tmin,tsec,x,ndim,
     &datform)
      dimension x(ndim),y(60),ix(60)
      character iform*20,name*48,text*72,pro*80,code*4,rest*40
      character prolog(48)*60,code1*6,code2*6,datform*3
      character zeile*120,b,zarr(120),zeil*40
      common /asl/ code1,code2,year,day,thr,tm,ts,t1000,srate,cpv
      real*8 srate8
      logical first
c
      b=' '
      cpv=1.000
      do j=1,n
        x(j)=0.
      enddo

      write(6,*) '==='
      write(6,*) ' Opening file ',name
      open(7,file=name,status='old')
c  read ASCII files
      if(datform.eq.'asc'.or.datform.eq.'sca') then
        write(6,22) trim(name),0,ndim,dt
        do i=1,lin
          read(7,'(a)') text
          write(6,*) ' header line: ',trim(text)
          prolog(i)=text
        enddo
        npro=lin
        text=name
        if(datform.eq.'sca') then
c single-column ascii
          do j=1,n
           read(7,*,err=232,end=231) x(j)
          enddo
  231     n=j-1
          write(6,*) n,' samples read from file ',trim(name)
          return
  232 write(6,*) 'error (maybe a non-numeric character) in data line '
     &  ,n,' of file ',trim(name)
          stop
        
        else
c multi- and variable-column ascii
c note: zeile is padded with blanks, therefore there is always a blank after the last digit
          nz=lin
          nval=0
          call cpu_time(cputime1)
          do
            read(7,'(a)',end=131) 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
            nz=nz+1
            read(zeile,*,err=132) (y(j),j=1,nn)
            do j=1,nn
              x(nval+j)=y(j)
            enddo
            nval=nval+nn
          enddo
  131     write(6,*) nval,' samples read from file ',trim(name)
          n=min(n,nval)
c          call cpu_time(cputime2)
c          write(*,'(a,f7.3,a)') ' in ',cputime2-cputime1,' seconds'
          close(7)
          return
  132   write(6,*) 'error (maybe a non-numeric character) in data line '
     &    ,nz,' of file ',trim(name)
          stop
        endif  

c  read data in SEIFE format
      else if(datform.eq.'fil') then
        read(7,'(a)') text
        write(6,*) ' header line: ',trim(text)
        npro=0
   21   read(7,'(a)') pro
        if(pro(1:1).ne.'%'.and.pro(1:1).ne.'}') goto 20
        npro=npro+1
        if(npro.gt.48) stop ' cannot handle more than 48 header lines'
        prolog(npro)=pro(3:62)
        goto 21
   20   read(pro,1) nn,iform,dt,tmin,tsec
    1   format(i10,a20,3f10.3)
        write(6,22) trim(name),nn,ndim,dt
   22   format('  reading file ',a,' ',i9,'/',i9,' samples','  dt=',
     & f10.6)
        if(n.lt.nn.or.nn.eq.0) then
c         write(6,'(" number of samples limited to",i7)') n
          nn=n
        endif
        n=nn
        if(iform(1:3).eq.'fre'.or.index(iform,'i').gt.0.or.
     &                            index(iform,'*').gt.0) then
          write(6,*) ' data will be read in free-format'
          read(7,*,err=25,end=23) (x(j),j=1,n)
        else
          read(7,iform,end=23) (x(j),j=1,n)
        endif

c  read data in MARS88 format
      else if(datform.eq.'mar') then
        read(7,'(a)') text
        text(61:72)='            '
        read(text(15:17),*) idt
        dt=idt/1000.
        read(text(53:54),*) thr
        read(text(56:57),*) tmin
        read(text(59:60),*) tsec
        tmin=tmin+60.*thr
        write(6,22) trim(name),n,ndim,dt
        do j=1,n
          read(7,*,err=23,end=23) x(j)
        enddo

      else if(datform.eq.'asl') then
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(1)),'(a)') code1
        read(zeile(ix(1)+2:ix(2)),'(a)') code2
        read(zeile(ix(2)+2:ix(7)),*) iyear,iday,ithr,itmin,itsec
        read(zeile(ix(7)+2:ix(11)),*) t1000,srate8,n,cpv 
        srate=srate8
        thr=ithr
        tmin=itmin
        tsec=itsec
        dt=nint(1.d6/srate8)/1.d6
        tmin=tmin+60.*ithr
        tsec=tsec+t1000/1000.
        nz=1
        nval=0
        do
          read(7,'(a)',end=31) zeile
c          read(zeile,'(120a)') zarr
          nn=0
          call cpu_time(cputime1)
          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)
        n=min(n,nval)
c        call cpu_time(cputime2)
c        write(*,'(a,f7.3,a)') ' in ',cputime2-cputime1,' seconds'
        close(7)
        return
   32   write(6,*) 'error (maybe a non-numeric character) in data line '
     &  ,nz,' of file ',trim(name)
        stop
 
      else if(datform.eq.'psn') then
c read PSN single-column ascii with multi-line header
          read(7,'(a)') zeil
          if(zeil(1:11).ne.'! PSN ASCII') stop 'not a PSN ASCII file'
          do
            read(7,'(a)') zeil
            if(zeil(1:11).eq.'Start Time:') then
              read(zeil(24:25),*) thr
              read(zeil(27:28),*) tmin
              read(zeil(30:35),*) tsec
              tmin=tmin+60.*thr
            else if(zeil(1:4).eq.'SPS:') then
              read(zeil(6:9),*) srate8
              srate=srate8
              dt=nint(1.d6/srate8)/1.d6
            else if(zeil(1:18).eq.'Number of Samples:') then
              read(zeil(20:30),*) n
            else if(zeil(1:5).eq.'Data:') then
              exit
            endif
          enddo
          do j=1,n
           read(7,*,err=233,end=234) x(j)
          enddo
          write(6,*) n,' samples read from file ',trim(name)
          return
  234     n=j-1
          write(6,*) n,' samples read from file ',trim(name)
          return
  233 write(6,*) 'error (maybe a non-numeric character) in data line '
     &  ,n,' of file ',trim(name)
          stop

      else if(datform.eq.'bdf') then
c read data in BDF format
        first=.true.
        nnn=0
   40   read(7,'(a)',end=23) text
        code=text(1:4)
        rest=text(5:44)
c        write(*,*) code,rest
        if(code.eq.'DATA') then
          goto 41
        else if(code.eq.'RATE') then
          read(rest,*) srate8
          srate=srate8
        dt=nint(1.d6/srate8)/1.d6
        else if(first.and.code.eq.'TIME') then
          read(rest(1:4),*) year
          read(rest(6:8),*) day
          read(rest(10:11),*) thr
          read(rest(13:14),*) tmin
          read(rest(16:22),*) tsec
          tmin=tmin+60.*thr
          first=.false.
        else if(code.eq.'NSAM') then
          read(rest,*) nsam
        endif
        goto 40
   41   read(7,*,err=26,end=23) (x(j),j=nnn+1,nnn+nsam)
        nnn=nnn+nsam
        goto  40
      endif
      close(7)
      return
   23 n=j-1
      write(6,24) n
   24 format(' end of file after ',i8,' samples')
      close(7)
      return
   25 write(6,26) trim(name),j
   26 format(' Input error reading file ',a,' at sample # ',i7)
      stop
      end

      subroutine outform(name,par,text,prol,npro,ncl,n,dt,tmin,tsec,x,
     & cfgline,lines)
c  formatted output
      dimension x(n)
      character name*48,text*72,par*56,iform*20,fform*22
      character*6 code1*6,code2*6,prol(48)*60,cfgline(32)*256     
      common /asl/ code1,code2,year,day,thr,tm,ts,t100,srate,cpv
      logical ncl
      if(npro.gt.48) stop ' cannot handle more than 48 header lines'
      iform='                    '
      xmax=0.
      do 2 j=1,n
    2 xmax=max(xmax,abs(x(j)))
      if(index(par,'lab').gt.0) then
        nvor=log10(max(xmax/cpv,2.))+1.
        ndec=max(0,10-nvor)
      else
        nvor=log10(max(xmax,2.))+1.
        ndec=max(0,10-nvor)
      endif
      i1=index(par,'(')
      i2=index(par,')')
      if(trim(par).eq.'') par(1:3)='fre'
      if(index(par,'fre').gt.0.) then
        write(iform,3) ndec
    3   format('(5f13.',i1,')')
      else if(index(par,'asc').gt.0.) then
        write(iform,7) ndec
    7   format('(f13.',i1,')')  
      else if(index(par,'lab').gt.0) then
        write(iform,4) ndec
    4   format('(1x,f10.3,1x,f13.',i1,') ')
      else if(index(par,'nil').gt.0) then
        return
      else if(i1.gt.0.and.i2.gt.i1+4) then
        read(par(i1:i2),'(a)') iform
      else
        write(*,*) ' illegal format ',par
        stop
      endif
      write(6,5) trim(name),n,iform
    5 format(' writing file ',a/i9,' samples in format ',a)
      if(index(par,'lab').gt.0) then
        open(8,file=name,status='new')
        write(8,iform) ((j-1)*dt,x(j)/cpv,j=1,n)
        call cfgfile(name,n,cfgline,lines,xmax/cpv,dt)
      else if(index(par,'asc').gt.0.) then
        open(8,file=name)
        write(8,iform) (x(j),j=1,n)
      else
        ndf=min(8,8-floor(alog10(dt)))
        write(fform,'(a13,i1,a8)') '(i10,a20,f10.',ndf,',2f10.3)'
        open(8,file=name)
        write(8,'(a)') text
        if(.not.ncl) write(8,'("% ",a)') (prol(j),j=1,npro)
        write(8,fform) n,iform,dt,tmin,tsec
    1   format(i10,a20,f10.5,2f10.3)
        write(8,iform) (x(j),j=1,n)
      endif
      close(8)
      write(6,'(" finished!")')
      return
   99 stop ' insufficient or wrong parameters for end command'
      end


      subroutine plotwindow(nil,par,npl1,npl2,n,ntrend)
      dimension x(n)
      character*56 par
      logical nil
      read(par,*,err=99,end=99) npl1,npl2
      npl1=max(npl1,1)
      npl2=min(npl2,n)
      if(npl1.gt.npl2) then
        stop 'plw: impossible plot window'
      endif
      write(6,*) 'plotwindow set: ',npl1,npl2
      nil=.false.
      ntrend=1
      return
   99 stop 'plw parameters missing or unreadable'
      end


      subroutine cfgfile(name,n,cfgline,lines,xmax,dt)
c     Erzeugt eine Zeile des cfg-Files
      character name*48,code1*6,code2*6
      character*256 cfgline(32)
      common /asl/ code1,code2,year,day,thr,tm,ts,t100,srate,cpv
      lines=lines+1
      g=1.
      ifilt=0
      tlp=2.5*dt
      thp=0.
      write(cfgline(lines),10) trim(name),2.*xmax,g,srate,ifilt,
     & tlp,thp,1./dt,n 
   10 format(a,'; Range:',e10.3,'; Gain:',f6.1,'; A/DFreq:',f8.3,
     & '; A/DFiltCode:',i2,'; BWTPCornerPer:',f8.3,'; BWHPCornerPer:',
     & f8.3,'; PrintFreq:',f8.3,'; Samples:',i7)
      return
      end

      subroutine gauss(aik,m,n,rs,f)
c  solve linear equations
      implicit real*8 (a-h,o-z)
      dimension aik(n,n),rs(n),f(n),h(14),imax(13)
      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)
      return
      end

      subroutine rekfl(x,y,n,f0,f1,f2,g1,g2)
c  perform recursive filtering
      implicit real*8 (a-h,o-z)
      real*4 x(n),y(n)
      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)=sngl(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 real*8 (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 files(nam1,nam2,fdt,lin,nseis,zeile,npts,nch,ncl,npl,
     & npw,pli,frq,cfgname,inpath,outpath,datform)
      character nam1(32)*48,nam2(32)*48,zeile*99,z3*3
      character*48 inpath,outpath,cfgname
      character*3 datform(32)
      dimension fdt(32),lin(32)
      logical ncl,npl,npw,pli,nopli,frq
      data nopli /.false./
      nseis=0
  101 read(7,'(a)',end=2) zeile
      if(zeile(1:1).ne.' ') write(6,'(1x,a)') trim(zeile)
      z3=zeile(1:3)
      if(z3.eq.'fil'.or.z3.eq.'asc'.or.z3.eq.'asl'.or.z3.eq.'mar'
     & .or.z3.eq.'bdf'.or.z3.eq.'sca'.or.z3.eq.'psn') then
        nseis=nseis+1
        datform(nseis)=z3
        if(z3.ne.'fil') nopli=.true.
       if(z3.eq.'mar'.or.z3.eq.'asl'.or.z3.eq.'bdf'.or.z3.eq.'psn') then
            lin(nseis)=1
        endif
        do 30 i=5,99
   30     if (zeile(i:i).ne.' '.and.zeile(i:i).ne.',') goto 31  
   31   ii=i
        do 11 i=ii,99
   11     if (zeile(i:i).eq.' '.or.zeile(i:i).eq.',') goto 34
   34   if(i.gt.ii+48) goto 24
        nam1(nseis)=zeile(ii:i-1)
        ii=i
        do 12 i=ii,99
   12     if (zeile(i:i).ne.' '.and.zeile(i:i).ne.',') goto 32  
   32   ii=i
        do 13 i=ii,99
   13     if (zeile(i:i).eq.' '.or.zeile(i:i).eq.',') goto 33
   33   if(i.gt.ii+48) goto 24
        nam2(nseis)=zeile(ii:i-1)
        if(z3.eq.'asc'.or.z3.eq.'sca') then
          ii=i
          do 14 i=ii,99
   14       if (zeile(i:i).ne.' '.and.zeile(i:i).ne.',') goto 4
    4     if(i.gt.97) goto 26
c          if(z3.eq.'asc') then
            read(zeile(i:99),*,err=26) fdt(nseis),lin(nseis)
c          else
c            read(zeile(i:99),*,err=26) lin(nseis)
c            lin(nseis)=-lin(nseis)
c          endif
        endif
      else if(z3.eq.'lim') then
        read(zeile,'(4x,i10)') npts
      else if(z3.eq.'chp') then
        read(zeile,'(4x,i10)') nch
      else if(z3.eq.'ncl') then
        ncl=.true.
      else if(z3.eq.'npl') then
        npl=.true.
      else if(z3.eq.'npw') then
	npw=.true.
      else if(z3.eq.'pli') then
	pli=.true.
      else if(z3.eq.'pfi') then
        inpath=zeile(5:28)
      else if(z3.eq.'pfo') then
        outpath=zeile(5:28)
      else if(z3.eq.'cfi') then
        cfgname=zeile(5:28)
      else if(z3.eq.'frq') then
         frq=.true.
      else if(z3.eq.'end') then
        if(pli.and.nopli) write(6,*) ' pli command will be ignored'
        pli=pli.and..not.nopli
        if(nseis.gt.0) then
c Pfade vor die Filenamen setzen
          if(inpath(1:1).ne.'#') then
            do i=1,nseis
              nam1(i)=trim(inpath)//'\'//nam1(i)
            enddo
          endif
          if(outpath(1:1).ne.'#') then
            do i=1,nseis
              nam2(i)=trim(outpath)//'\'//nam2(i)
            enddo
          endif
          return
        else
          stop ' Filename missing'
        endif
      endif
      goto 101
    2 stop 'missing end line or missing <cr> after end line'
   24 write(*,*) ' Filename too long: ',zeile
      stop
   26 stop ' asc or sca statement incorrect or incomplete'
      end

      subroutine factor(nil,par,x,n)
c  multiply by a constant factor
      dimension x(n)
      character*56 par
      logical nil
      read(par,*,err=99,end=99) f
      do 1 j=1,n
    1 x(j)=f*x(j)
      write(6,'(" fac  ",e13.6)') f
      nil=.false.
      return
   99 stop ' parameter for fac missing or unreadable'
      end

      subroutine norm(nil,typ,par,x,n)
c  normalize rms, peak, or peak-to-peak amplitude
      dimension x(n)
      character typ*3,par*56
      logical nil
      read(par,*,err=99,end=99) fn
      xrms=0.
      xmin=1e24
      xmax=-xmin
      do 1 j=1,n
        xrms=xrms+x(j)**2
        xmin=min(xmin,x(j))
    1   xmax=max(xmax,x(j))
      xrms=sqrt(xrms/n)
      xp=max(abs(xmin),abs(xmax))
      if(typ.eq.'nox') then
      f=2.*fn/(xmax-xmin)
      off=(xmax+xmin)/2.
        do 2 j=1,n
    2     x(j)=f*(x(j)-off)
      else
        if(typ.eq.'nor') then
          f=fn/xrms
        else
          f=fn/xp
        endif
        do 3 j=1,n
    3     x(j)=f*x(j)
      endif
      write(6,'(1x,a3,a9,e13.6)') typ,' Faktor ',f
      nil=.false.
      return
   99 write(*,*) ' illegal parameter for ',typ 
      stop
      end

      subroutine pad(nil,par,x,n,nd)
c  Nullen ergaenzen
      dimension x(nd)
      character*56 par
      logical nil
      read(par,*,err=99,end=99) nul
      if(nul.le.n.or.nul.gt.nd) goto 99
      do 1 j=n+1,nul
    1 x(j)=0.
      n=nul
      write(6,2) nul
    2 format(" pad: signal was padded with zeros to length",i8)
      nil=.false.
      return
   99 stop ' parameter for pad missing, illegal, or unreadable'
      end

      subroutine delay(nil,par,x,n)
c  vorne Nullen einfuegen
      dimension x(n)
      character*56 par
      logical nil
      read(par,*,err=99,end=99) nul
      if(nul.le.0) goto 99
      do 1 j=n,nul+1,-1
    1 x(j)=x(j-nul)
      do 3 j=nul,1,-1
    3 x(j)=0.
      write(6,2) nul
    2 format(" del: the signal was delayed by",i8," samples")
      nil=.false.
      return
   99 stop ' parameter for del missing, wrong, or unreadable'
      end

      subroutine plotpar(nil,par,breit,hoch)
c  Plotgroesse aendern
      character*56 par
      logical nil
      read(par,*,err=99,end=99) breit,hoch
      write(6,2) breit,hoch
    2 format(" plt: plotsize ",f5.1," * ",f5.1)
      nil=.false.
      return
   99 stop ' parameter for plt missing, wrong, or unreadable'
      end

       subroutine add(nil,par,x,n)
c  add a constant to each sample
      dimension x(n)
      character*56 par
      logical nil
      read(par,*,err=99,end=99) f
      do 1 j=1,n
    1 x(j)=f+x(j)
      write(6,'(" add  ",e13.6)') f
      nil=.false.
      return
   99 stop ' insufficient or wrong parameters for add'
      end

      subroutine distort(nil,par,x,n)
c  add nonlinear distortion
      dimension x(n)
      character*56 par
      logical nil
      read(par,*,err=99,end=99) k,a
      xmax=0.
      do 2 j=1,n
    2 xmax=max(xmax,abs(x(j)))
      do 1 j=1,n
    1 x(j)=x(j)+a*xmax*(x(j)/xmax)**k
      write(6,'(" distortion added",f12.6," * x**",i1)') a,k
      nil=.false.
      return
   99 stop ' insufficient or wrong parameters for dis'
      end

      subroutine zero(nil,par,x,n)
c  null time series from n1 to n2
      dimension x(n)
      character*56 par
      logical nil
      read(par,*,err=99,end=99) n1,n2
      n1=max(n1,1)
      n2=min(n2,n)
      if(n1.gt.n2) goto 99
      do 1 j=n1,n2
    1 x(j)=0.
      write(6,'(" null ",2i10)') n1,n2
      nil=.false.
      return
   99 stop ' insufficient or wrong parameters for null'
      end
      
      subroutine set(nil,par,x,n)
c  set time series from n1 to n2 to a
      dimension x(n)
      character*56 par
      logical nil
      read(par,*,err=99,end=99) n1,n2,a
      do 1 j=n1,n2
    1 x(j)=a
      write(6,'(" set ",2i10)') n1,n2
      nil=.false.
      return
   99 stop ' insufficient or wrong parameters for set'
      end

      subroutine clip(nil,par,x,n)
c  clip signal at level clp
      dimension x(n)
      character*56 par
      logical nil
      read(par,*,err=99,end=99) clp
      do 1 j=1,n
      x(j)=min(x(j),clp)
    1 x(j)=max(x(j),-clp)
      write(6,'(" clip ",e13.6)') clp
      nil=.false.
      return
   99 stop ' insufficient or wrong parameters for clip'
      end

      subroutine deriv(nil,par,x,n,dt)
c  derivative (non-recursive, non-causal, symmetric-difference)
      dimension x(n)
      character*56 par
      logical nil
      twodt=2.*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)
      write(6,'(" dif - derivative by symmetric differences")')
      nil=.false.
      return
      end

      subroutine interp(nil,par,x,n)
c  interpolate clipped samples
      real*8 aik(4,4),xx(4),f(4),t(4)
      dimension x(n)
      character*56 par
      logical nil
      data zero,one,two/0.d0,1.d0,2.d0/
c  define illegal samples
      read(par,*,err=99,end=99) clip,n1,n2,apex
      if(clip.eq.0.) clip=262000.
      n1=max(n1,1)
      if(n2.eq.0.or.n2.gt.n) n2=n
c  search for illegal samples
      naus=0
      nabt=0
   10 do 1 j=n1,n2
      if(abs(x(j)).gt.clip) goto 2
    1 continue
c  if no more, done
      nil=.false.
      write(6,7) naus,nabt,clip,apex
    7 format(" spl - in",i4," places, ",i5," values >",e13.6,
     & " were interpolated. apex =",f9.0)
      return
    2 j1=j
c  determine the last bad sample (two good ones must follow)
      do 3 j2=j1,n-2
      if(abs(x(j2+1)).lt.clip.and.abs(x(j2+2)).lt.clip) goto 4
    3 continue
      if(j1.lt.n1+2.or.j2.gt.n2-2) then
      stop ' *** spl ***  insufficient number of pivots'
      endif
    4 if(apex.eq.zero) then
c  compute spline
      xx(1)=x(j1-2)
      t(1)=-two
      xx(2)=x(j1-1)
      t(2)=-one
      xx(3)=x(j2+1)
      t(3)=j2-j1+1
      xx(4)=x(j2+2)
      t(4)=t(3)+one
      do 6 j=1,4
      aik(j,1)=one
      do 6 k=2,4
    6 aik(j,k)=t(j)*aik(j,k-1)
      else
c  compute parabola
      xx(1)=x(j1-1)
      t(1)=-one
      xx(3)=x(j2+1)
      t(3)=j2-j1+1
      xx(2)=apex
      t(2)=(t(1)+t(3))/two
      do 16 j=1,3
      aik(j,1)=one
      aik(j,4)=zero
      aik(4,j)=0.d0
      do 16 k=2,3
   16 aik(j,k)=t(j)*aik(j,k-1)
      aik(4,4)=one
      xx(4)=zero
      endif
      call gauss(aik,4,4,xx,f)
c  interpolate
      do 5 j=j1,j2
      tt=j-j1
    5 x(j)=f(1)+tt*(f(2)+tt*(f(3)+tt*f(4)))
      n1=j2+1
      naus=naus+1
      nabt=nabt+j2-j1+1
      goto 10
   99 stop ' insufficient or wrong parameters for spl'
      end

      subroutine window(nil,pli,typ,par,x,n,dt,tmin,tsec)
c  time window by sample number
      dimension x(n)
      character*56 par
      character*3 typ
      logical nil,pli
      if(pli.and.typ.ne.'tap') write(*,*) 
     &   " WARNING   when using PLI, don't apply any time window"
      if(typ.eq.'skp'.or.typ.eq.'cos'.or.typ.eq.'tap') then
        read(par,*,err=99,end=99) np1
        write(6,*) typ,np1
      else if(typ.eq.'swi') then
        read(par,*,err=99,end=99) sec1,sec2
        np1=max(nint(sec1/dt)+1,1)
        np2=min(nint(sec2/dt),n)
        write(6,*) typ,(np1-0.5)*dt,(np2-0.5)*dt
        np1=np1-1
        n=np2-np1
      else
        read(par,*,err=99,end=99) np1,np2
        np1=max(np1,1)
        np2=min(np2,n)
        write(6,*) typ,np1,np2
        np1=np1-1
        n=np2-np1
      endif
c  *** tap
      if(typ.ne.'tap') goto 16
      pi=4.*atan(1.)
      if(np1.gt.n/2) write(6,*) 'np1 set to n/2'
      np1=min(np1,n/2)
      fak=pi/np1
      do 14 j=1,np1
      tap=0.5-0.5*cos(fak*(j-1))
      x(j)=x(j)*tap
   14 x(n-j+1)=x(n-j+1)*tap
      acor=1./sqrt(1.-np1*1.25/n)
      write(6,19) acor
   19 format(' amplitude correction factor ',f6.3,' applied')
      do 18 j=1,n
   18 x(j)=x(j)*acor
      goto 15
c  *** cos
   16 if(typ.ne.'cos') goto 13
      n1=n+1
      pih=2.*atan(1.)
      fak=pih/np1
      do 17 j=1,np1
      taper=sin(fak*j)**2
      x(j)=x(j)*taper
   17 x(n1-j)=x(n1-j)*taper
      goto 15
c  *** win, swi
   13 if(typ.eq.'win'.or.typ.eq.'swi') then
      do 11 j=1,n
   11 x(j)=x(np1+j)
      endif
c  *** sin, sis
      if(typ.eq.'sin'.or.typ.eq.'sis') then
      pi=4.*atan(1.)
      f=pi/(n+1)
      w2=sqrt(2.)
      nex=1
      if(typ.eq.'sis') nex=2
      do 12 j=1,n
   12 x(j)=w2*x(np1+j)*sin(f*j)**nex
      endif
      if(typ.eq.'skp') then
      n=n-np1
      do j=1,n
      x(j)=x(np1+j)
      enddo
      endif
      tsec=60.*tmin+tsec+np1*dt
      tmin=int(tsec/60.)
      tsec=tsec-60.*tmin
   15 nil=.false.
      return
   99 write(*,*) ' insufficient or wrong parameters for ',typ
      stop
      end

      subroutine timwin(nil,par,x,n,dt,tmin,tsec)
c  time window by time after midnight
      dimension x(n)
      character*56 par
      logical nil
      read(par,*,err=99,end=99) tmin1,tsec1,tmin2,tsec2
      t0=60.*tmin+tsec
      t1=60.*tmin1+tsec1
      t2=60.*tmin2+tsec2
      np1=1.5+(t1-t0)/dt
      np2=0.5+(t2-t0)/dt
      if(np1.lt.1.or.np2.gt.n.or.np2.lt.np1) then
      write(6,1) tmin1,tsec1,tmin2,tsec2
    1 format(" twi  ",4f10.1/"      *** impossible time window ***")
      if(np1.lt.1) then
      tsec1=tsec1+(1-np1)*dt
      np1=1
      endif
      if(np2.gt.n) then
      tsec2=tsec2-(np2-n)*dt
      np2=n
      endif
      if(np2.lt.np1) stop ' impossible time window'
      write(6,'(" new window from sample ",i7," to",i7)') np1,np2
      endif
      write(6,2) tmin1,tsec1,tmin2,tsec2
    2 format(' twi  ',4f10.1)
      np1=np1-1
      n=np2-np1
      do 3 j=1,n
    3 x(j)=x(np1+j)
      tsec=t0+np1*dt
      tmin=int(tsec/60.)
      tsec=tsec-60.*tmin
      nil=.false.
      return
   99 stop ' insufficient or wrong parameters for twi'
      end


      subroutine timeset(nil,par,tmin,tsec)
      character*56 par
      logical nil
      read(par,*,err=99,end=99) tmin,tsec
      write(6,*) ' time stamp set to tmin=',tmin,', tsec=',tsec
      nil=.false.
      return
   99 stop ' insufficient or wrong parameters for tim'
      end


      subroutine decim(nil,par,x,n,dt)
c  decimate time series by a factor up to 100
      dimension x(n),g(200)
      character*56 par
      logical nil
      read(par,*,err=99,end=99) idez
      if(idez.gt.200) stop ' cannot decimate by more than 200'
      pih=2.*atan(1.)
      fidez=idez
      do 3 j=1,idez
    3 g(j)=(cos((j-idez)/fidez*pih))**2/fidez
      n=n/idez
      dt=dt*idez
      do 5 k=1,n-1
      jk=k*idez+1
      x(k)=g(idez)*x(jk)
      do 6 l=1,idez-1
    6 x(k)=x(k)+g(l)*(x(jk-idez+l)+x(jk+idez-l))
    5 continue
      do 7 j=n,2,-1
    7 x(j)=x(j-1)
      write(6,'(" time series was decimated by ",i10)') idez
      nil=.false.
      return
   99 stop ' insufficient or wrong parameters for dec'
      end
      

      subroutine mean(nil,par,x,n)
c  remove average
      character*56 par
      dimension x(n)
      logical nil
      if(trim(par).eq.'') write(par,*) 0
      read(par,*,err=99,end=99) n2
      if(n2.lt.1.or.n2.gt.n) n2=n
      sum=0.
      do 1 j=1,n2
    1 sum=sum+x(j)
      sum=sum/n2
      do 2 j=1,n
    2 x(j)=x(j)-sum
      write(6,'(" avg  ",e13.6)') sum
      nil=.false.
      return
   99 stop ' insufficient or wrong parameters for avg'
      end
      
      
      subroutine quant(nil,par,x,n)
c  quantization
      character*56 par
      dimension x(n)
      logical nil
      if(trim(par).eq.'') write(par,*) 1.
      read(par,*,err=99,end=99) a
      if(a.le.0.) goto 99
      do 1 j=1,n
    1 x(j)=a*anint(x(j)/a)
      write(6,'(" qnt: signal was quantized to integer *",f10.6)') a
      nil=.false.
      return
   99 stop ' quantization interval must be positive'
      end


      subroutine trend(nil,typ,par,x,n)
c  remove trend
      implicit real*8 (a-h,o-z)
      real*4 x(n)
      logical nil
      character typ*3,par*56
      common /tre/ a,b
      if(typ.eq.'tre') then
        if(trim(par).eq.'') write(par,*) 0
        read(par,*,err=99,end=99) n2
        if(n2.le.1.or.n2.gt.n) n2=n
        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)
 1001     sjx = sjx+x(j)*j
        a = (sx*beta-sjx*alpha)/det
        b = (sjx*gn-sx*alpha)/det
      else if(typ.eq.'zsp') then
        n2=n
        gn = n
        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 1003 j=1,n2
          sx = sx+x(j)
 1003     sjx = sjx+x(j)*j
        a = (sx*beta-sjx*alpha)/det
        b = (sjx*gn-sx*alpha)/det
        return
      else
        b=(x(n)-x(1))/float(n-1)
        a=x(1)-b
      endif
      do 1002 j=1,n
 1002   x(j) = x(j)-sngl(a+b*j)
      write(6,'(" tre  ",2f13.6)') a,b
      nil=.false.
      return
   99 write(*,*) ' insufficient or wrong parameters for ',typ
      stop
      end


      subroutine shift(nil,typ,par,x,n)
c  shift signal left or right
      implicit real*8 (a-h,o-z)
      real*4 x(n)
      logical nil
      character typ*3,par*56
      read(par,*,err=99,end=99) ns
      if(ns.eq.0) then
        nil=.false.
        return
      else if(ns.lt.0) then
        goto 99
      endif
      if(typ.eq.'shl') then
        do i=1,n-ns
          x(i)=x(i+ns)
        enddo
        do i=n-ns+1,n
          x(i)=0.
        enddo
      else if(typ.eq.'shr') then
        do i=n,ns+1,-1
          x(i)=x(i-ns)
        enddo
        do i=ns,1,-1
          x(i)=0.
        enddo
      endif
      write(6,'(1x,a3,i7,"signal shifted in time")') typ,ns
      nil=.false.
      return
   99 write(*,*) ' insufficient or wrong parameters for ',typ
      stop
      end

      subroutine zspline(nil,typ,par,x,n)
c  remove z-trend with a z-spline
      implicit real*8 (a-h,o-z)
      real*4 x(n),fn32,x2,x3
      logical nil
      character typ*3,par*56
      common /tre/ a,b
      read(par,*,err=99,end=99) 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
        stop ' illegal parameters for zspline'
      endif
      call trend(nil,typ,par,x(n1),n2-n1+1)
      a1=a
      b1=b
      call trend(nil,typ,par,x(n3),n4-n3+1)
      a3=a
      b3=b
      do 2 j=1,n2
    2 x(j) = x(j)-sngl(a1+b1*j)
      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)-sngl(a3+b3*(j-n3+1))
      write(6,'(" zspline eck ",4i10,)') n1,n2,n3,n4
      write(6,'(" zspline coef",4f10.6)') a1,b1,a3,b3
      nil=.false.
      return
   99 stop ' insufficient or wrong parameters for zspline'
      end


      subroutine polytrend(nil,par,x,n)
c  remove polynomial trend
      parameter(ndi=13)
      implicit real*8 (a-h,o-z)
      real*4 x(n)
      dimension b(ndi),c(ndi,ndi),a(ndi)
      logical nil
      character*56 par
      read(par,*,err=99,end=99) m
      m=min(m,ndi-1)
      fnh=n/2.
      one=1.d0

      do j=1, m+1
        do k=1, m+1
          c(j,k)=0
          do i=1, n
            c(j,k)=c(j,k)+(dble(i)/fnh-one)**(j+k-2)
          enddo
        enddo
        b(j)=0
        do i=1,n
          b(j)=b(j)+(dble(i)/fnh-one)**(j-1)*x(i)
        enddo
      enddo
      call gauss(c,m+1,ndi,b,a)
      write(6,100) (j-1,a(j),j=1,m+1)
  100 format(1x,i5,e15.6)
      do i=1,n
        xpol=a(m+1)
        do j=m,1,-1
          xpol=xpol*(dble(i)/fnh-one)+a(j)
        enddo
        x(i)=x(i)-sngl(xpol)
      enddo

      write(6,'(" pol  ",i5)') m
      nil=.false.
      return
   99 stop ' insufficient or wrong parameters for pol'
      end

      subroutine expon(nil,par,x,n)
c  remove exponential trend
      implicit real*8 (a-h,o-z)
      real*4 x(n)
      logical nil
      character*56 par
      dimension sum(3)
      read(par,*,err=99,end=99) nex
      if(nex.lt.3.or.nex.gt.n) nex=n
      n3=nex/3
      do 1 i=1,3
        sum(i)=0.d0
        ii=(i-1)*n3
        do 1 j=1,n3
  1       sum(i)=sum(i)+x(j+ii)
      q=(sum(2)-sum(3))/(sum(1)-sum(2))
      c=-dlog(q)/dble(n3)
      b=c*(sum(1)-sum(2))/(1.d0-q)**2
      a=(sum(1)+sum(2)+sum(3)-b/c*(1.d0-q**3))/dble(n3)/3.d0
      do 2 j=1,n
  2     x(j) = x(j)-sngl(a+b*dexp(-c*dble(j)))
      write(6,'(" exp  ",2f10.3,f10.6)') a,b,c
      nil=.false.
      return
   99 stop ' insufficient or wrong parameters for exp'
      end


      subroutine sig(nil,zeile,x,n)
c  signum function or absolute value
      real*4 x(n)
      logical nil
      character zeile*99
      if(zeile(1:3).eq.'sig') then
        do i=1,n
          x(i)=sign(1.,x(i))
        enddo
      write(6,'(a)') ' sig: signum function applied'
        else
        do i=1,n
          x(i)=abs(x(i))
        enddo
      write(6,'(a)') ' abs: absolute value taken'
      endif
      nil=.false.
      return
      end


      subroutine reverse(nil,par,x,n)
c  reverse signal in time
      logical nil
      dimension x(n)
      character*56 par
      do j=1,n/2
        jj=n+1-j
        xx=x(j)
        x(j)=x(jj)
        x(jj)=xx
      enddo
      write(6,'(" rev - signal was reversed")')
      nil=.false.
      return
      end

      subroutine tides(nil,par,x,n,dt)
c  remove tides. number of frequencies is automatically chosen according
c  to the total length of the record.
      implicit real*8 (a-h,o-z)
      real*4 x(n)
      logical nil
      character*56 par
      dimension a(13,13),rs(13),c(13),d(13),e(13),f(13),
     & omega(6),omeg(6),cor(6)
      data omega/1.93227d0,0.92954d0,2.00000d0,1.00274d0,1.89567d0,
     & 0.89293d0/
      data zero,one,two/0.d0,1.d0,2.d0/
      ndi=13
      dth=dt/two
      read(par,*,err=99,end=99) nstep
      if(nstep.eq.0) nstep=300./dts
      nstep=max(nstep,1)
      step=nstep
      tstep2=(step-one)*dth
      tint=tstep2+dth
c  determine the number of frequencies required for a good fit
      dur=n*dt/3600.d0
      nfreq=6
      if(dur.lt.35.d0) nfreq=5
      if(dur.lt.18.d0) nfreq=4
      if(dur.lt.14.d0) nfreq=3
      if(dur.lt. 5.d0) nfreq=2
      nco=2*nfreq+1
      omeg0=8.d0*datan(one)/86400.d0
      do 101 i=1,nfreq
  101 omeg(i)=omeg0*omega(i)
c  determine partial amplitudes by least-squares fit
      do 1 i=1,nco
      rs(i)=zero
      do 1 k=1,nco
    1 a(i,k)=zero
c  correction for averaging over nstep samples
      do 102 j=1,nfreq
  102 cor(j)=step*dsin(omeg(j)*dth)/dsin(step*omeg(j)*dth)
c  set up system of linear equations
      c(1)=one
      do 2 j=1,n-nstep+1,nstep
      sx=zero
      do 12 jj=j,j+nstep-1
   12 sx=sx+x(jj)
      sx=sx/step
      t=(j-1)*dt+tstep2
      do 103 k=1,nfreq
      c(2*k)  =dcos(omeg(k)*t)
  103 c(2*k+1)=dsin(omeg(k)*t)
      do 2 i=1,nco
      rs(i)=rs(i)+sx*c(i)
      do 2 k=1,nco
    2 a(i,k)=a(i,k)+c(i)*c(k)
c  solve for partial amplitudes
      call gauss(a,nco,ndi,rs,f)
      do 3 j=1,n-nstep+1,nstep
      t=(j-1)*dt+tstep2
c  remove average and tides
      do 104 k=1,nfreq
      k2=2*k
      k3=k2+1
      c(k2)=dcos(omeg(k)*t)
      c(k3)=dsin(omeg(k)*t)
      d(k2)=-omeg(k)*c(k3)
      d(k3)= omeg(k)*c(k2)
      e(k2)=-omeg(k)*d(k3)
  104 e(k3)= omeg(k)*d(k2)
      xgez1=f(1)
      xgez2=zero
      xgez3=zero
      do 105 k=1,nfreq
      k2=2*k
      k3=k2+1
      xgez1=xgez1+cor(k)*(f(k2)*c(k2)+f(k3)*c(k3))
      xgez2=xgez2+cor(k)*(f(k2)*d(k2)+f(k3)*d(k3))
  105 xgez3=xgez3+cor(k)*(f(k2)*e(k2)+f(k3)*e(k3))/two
      if(j.gt.n-2*nstep+1) nstep=n+1-j
      do 3 jj=j,j+nstep-1
      tdif=(jj-j)*dt-tstep2
    3 x(jj)=x(jj)-sngl(xgez1+xgez2*tdif+xgez3*tdif*tdif)
      write(6,9) f(1),(omega(j),f(2*j),f(2*j+1),j=1,nfreq)
    9 format(" tid   avg=",f10.0/(5x,"freq=",f10.5,"  cos=",
     & e11.3,"  sin=",e11.3))
      nil=.false.
      return
   99 stop ' insufficient or wrong parameters for tides'
      end
c
      subroutine trans(nil,typ,par,x,n,dt)
c  remove transient response
      implicit real*8 (a-h,o-z)
      real*4 x(n),dts
      logical nil
      character typ*3, par*56
      dimension a(3,3),rs(3),c(3),f(3)
      data zero,one/0.d0,1.d0/
      twopi=8.d0*datan(one)
      if(typ.eq.'tr3') then
        read(par,*,err=99,end=99) n1,t0,h
        if(h.ge.one) then
          stop ' tra: h must be less than 1'
        endif
c  determine partial amplitudes by least-squares fit
      wh=dsqrt(one-h**2)
      omeg=twopi/t0
      do 1 i=1,3
      rs(i)=zero
      do 1 k=1,3
    1 a(i,k)=zero
c  set up system of linear equations
      c(1)=one
      do 2 j=n1,n
      t=(j-1)*dt
      ex=dexp(-omeg*h*t)
      c(2)=dcos(omeg*wh*t)*ex
      c(3)=dsin(omeg*wh*t)*ex
      do 2 i=1,3
      rs(i)=rs(i)+x(j)*c(i)
      do 2 k=1,3
    2 a(i,k)=a(i,k)+c(i)*c(k)
c  solve for partial amplitudes
      call gauss(a,3,3,rs,f)
      else
c  use previously determined partial amplitudes
        read(par,*,err=99,end=99) n1,t0,h,f(1),f(2),f(3)
      wh=dsqrt(one-h**2)
      omeg=twopi/t0
      endif 
      do 5 j=1,n1-1
    5   x(j)=0.
      do 3 j=n1,n
      t=(j-1)*dt
      ex=dexp(-omeg*h*t)
      omt=omeg*wh*t
    3 x(j)=x(j)-sngl(f(1)+f(2)*dcos(omt)*ex+f(3)*dsin(omt)*ex)
      write(6,9) f
    9 format(" transient:  avg=",f13.3,"   f2,f3=",2f13.3)
      nil=.false.
      return
   99 stop ' insufficient or wrong parameters for trans'
      end


      subroutine filter(nil,zeile,x,n,dt,npw,frq)
      implicit real*8 (a-h,o-z)
      character typ*3,typ1*3,typ2*3,par*56,zeile*99,warp*13
      real*4 x(n)
      logical nil,npw,frq
      typ=zeile(1:3)
      par=zeile(5:60)
      pi=4.d0*datan(1.d0)
c  Subtract value of first sample before high-pass filtration
      if(typ(1:1).eq.'h') 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,*,err=99,end=99) t0
      if(frq) t0=1./t0
      it=11
      else if(typ.eq.'hp1') then
      read(par,*,err=99,end=99) t0
      if(frq) t0=1./t0
      it=12
      else if(typ.eq.'lp2') then
      read(par,*,err=99,end=99) t0,h
      if(frq) t0=1./t0
      it=21
      else if(typ.eq.'hp2') then
      read(par,*,err=99,end=99) t0,h
      if(frq) t0=1./t0
      it=22
      else if(typ.eq.'bp2') then
      read(par,*,err=99,end=99) t0,h
      if(frq) t0=1./t0
      it=24
      else if(typ.eq.'int') then
      it= 1
      else if(typ.eq.'he1'.or.typ.eq.'le1') then
      read(par,*,err=99,end=99) t0s,t0
      if(frq) t0s=1./t0s
      if(frq) t0=1./t0
      it=13
      else if(typ.eq.'he2'.or.typ.eq.'le2') then
      read(par,*,err=99,end=99) t0s,hs,t0,h
      if(frq) t0s=1./t0s
      if(frq) t0=1./t0
      it=23
      else if(typ.eq.'sez') then
        read(par,*,err=99,end=99) t0s,hs
      if(frq) t0s=1./t0s
        t0=1.e12
        h=1.
        it=23
      else
      goto 2
      endif
c  prewarp and calculate filter weights
        if(it.eq.1) then
          t0p=1.d0/dt
          t0sp=t0p
        else
        if(npw) then
           t0p=t0/dt
           t0sp=t0s/dt
           warp='  no prewarp '
        else
          t0p=pi/dtan(pi/t0*dt)
          t0sp=pi/dtan(pi/t0s*dt)
          warp='  prewarp    '
        endif
      endif
      call rfk(it,t0p,h,t0sp,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
      call rekfl(x,x,n,f0,f1,f2,g1,g2)
      nil=.false.
c  confirm execution
      if(it.eq.1) write(6,'(1x,a3)') typ
      if(it.eq.11.or.it.eq.12) write(6,8) typ,warp,t0
      if(it.eq.13) write(6,9) typ,warp,t0s,t0
    7 format(1x,a3,a13,2f10.3,f16.3,f10.3)
    8 format(1x,a3,a13,5f10.3)
    9 format(1x,a3,a13,f10.3,f16.3)
      if(it.eq.21.or.it.eq.22.or.it.eq.24) write(6,8) typ,warp,t0,h
      if(it.eq.23) then
        if(typ.eq.'sez') then
          write(6,8) typ,warp,t0s,hs
        else
          write(6,7) typ,warp,t0s,hs,t0,h
        endif
      endif
      write(*,11) " f0,f1,f2,g1,g2= ",f0,f1,f2,g1,g2
   11 format(a,5f12.6)
      write(*,*)
      return

c  Butterworth filters
    2 if(typ.eq.'lpb'.or.typ.eq.'hpb') then
        read(par,*,err=99,end=99) t0,m
        if(frq) t0=1./t0
        mm=m/2
        if(typ.eq.'lpb') then
          it1=11
          it2=21
          typ1='lp1'
          typ2='lp2'
        else
          it1=12
          it2=22
          typ1='hp1'
          typ2='hp2'
        endif
        if(m.gt.2*mm) then

c  prewarp and calculate filter weights
        if(it.eq.1) then
          t0p=1.d0/dt
          t0sp=t0p
        else
        if(npw) then
           t0p=t0/dt
           t0sp=t0s/dt
           warp='  no prewarp '
        else
          t0p=pi/dtan(pi/t0*dt)
          t0sp=pi/dtan(pi/t0s*dt)
          warp='  prewarp    '
        endif
      endif

          call rfk(it1,t0/dt,h,t0s/dt,hs,f0,f1,f2,g1,g2)
          call rekfl(x,x,n,f0,f1,f2,g1,g2)
          write(6,8) typ1,warp,t0
        endif
        pih=2.d0*datan(1.d0)
        do 3 j=1,mm
          h=dsin(pih*(2*j-1)/m)
          
c  prewarp and calculate filter weights
        if(it.eq.1) then
          t0p=1.d0/dt
          t0sp=t0p
        else
        if(npw) then
           t0p=t0/dt
           t0sp=t0s/dt
           warp='  no prewarp '
        else
          t0p=pi/dtan(pi/t0*dt)
          t0sp=pi/dtan(pi/t0s*dt)
          warp='  prewarp    '
        endif
      endif          
          
          call rfk(it2,t0/dt,h,t0s/dt,hs,f0,f1,f2,g1,g2)
          call rekfl(x,x,n,f0,f1,f2,g1,g2)
          write(6,8) typ2,warp,t0,h
    3   continue
      nil=.false.
      endif
      return
   99 write(*,*) ' insufficient or wrong parameters for ',typ
      stop
      end

      subroutine extrem(nil,x,npl1,npl2,n)
      dimension x(n)
      logical nil
      nil=.false.
      xmin=1e96
      xmax=-xmin
      amin=xmin
      do 1 j=npl1,npl2
      if(abs(x(j)).gt.amin) goto 4
      amin=abs(x(j))
      namin=j
    4 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
      if(xmax.ge.-xmin) then
      amax=xmax
      namax=nmax
      else
      amax=-xmin
      namax=nmin
      endif
      write(6,*) 'ext: extrema between samples ',npl1,' and ',npl2,':'
      write(6,2) xmin,nmin,xmax,nmax
    2 format(" xmin=",e11.3," (j=",i7,")   xmax=",e11.3," (j=",i7,")")
      write(6,5) amin,namin,amax,namax
    5 format(" amin=",e11.3," (j=",i7,")   amax=",e11.3," (j=",i7,")")
      write(6,6) xmax-xmin
    6 format(" peak-to-peak ",e11.3)
      return
      end

      subroutine rms(nil,x,npl1,npl2,n)
      dimension x(n)
      logical nil
      nil=.false. 
      avg=0.
      do 3 j=npl1,npl2
    3   avg=avg+x(j)
      avg=avg/(npl2-npl1+1)
      xrms=0.
      do 1 j=npl1,npl2
    1 xrms=xrms+(x(j)-avg)**2
      write(6,2) npl1,npl2,sqrt(xrms/(npl2-npl1+1))
    2 format(" rms over de-meaned samples ",i8," to ",i8,": ",e11.3)
      return
      end

      subroutine square(nil,x,n)
      dimension x(n)
      logical nil
      nil=.false.
      do 1 j=1,n
    1 x(j)=x(j)**2
      write(6,2)
    2 format(" square: the signal has been squared")
      return
      end

      subroutine intpo(nil,par,dt1,x,np1)
      dimension x(np1)
      logical nil
      character*56 par
      nil=.false.
      read(par,*,err=99,end=99) dt2
      if(dt2.le.dt1) then
      stop ' sampling inteval can be increased but not decreased'
      endif
      j2=0
    1 j2=j2+1
      t=(j2-1)*dt2
      j10=t/dt1
      frac=t/dt1-j10
      j11=j10+1
      j12=j11+1
      if(j12.gt.np1) goto 2
      x(j2)=(1.-frac)*x(j11)+frac*x(j12)
      goto 1
    2 np1=j2-1
      write(6,3) dt1,dt2
    3 format(" csi  ",2f10.3)
      dt1=dt2
      return
   99 stop ' insufficient or wrong parameters for intpo'
      end
c
      subroutine heapsort(nil,ra,n)
      dimension ra(n)
      logical nil
      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
          goto 30
        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
 30   nil=.false.
      write(*,*) ' series was heap-sorted, is no longer a signal'
      return
      end