c   .    1    .    2    .    3    .    4    .    5    .    6    .    7..
c
c    program SINCAL1
c
c    Determination of the three parameters of an electromagnetic
c    seismometer: free period, damping and generator constant (gain)
c    from a record of 4 to 9 test frequencies.
c
c    The exact frequencies of the test signals must be specified
c    in the parameter file sincal.par. If they are not known, use
c    the program sincal2.
c
c    *******************************************************************
c
c    This program  assumes that you have recorded the response of a
c    simple seismometer (with a geophone-type transfer function) to
c    at least four sinusoidal test signals of equal amplitude. When
c    f0 is the eigenfrequency of the sensor, then the signal frequencies
c    0.2*f0, 0.5*f0, f0, 2*f0, and 5*f0 are a good choice. The amplitudes
c    must be equal in terms of calibration coil current or ground
c    acceleration. If your calibration-signal generator puts out
c    signals corresponding ground velocity, the amplitudes must be
c    chosen inversely proportional to the frequency. Exact values
c    of the test frequencies must be communicated to the program as the
c    first parameter t0 of the 'ban' commands in the 'sincal.par' file.
c
c    The records of the test signals must be contained in a single
c    data file. Each test signal should contain about 25 full periods,
c    depending on the half-width of the gaussian bandpass filters used
c    in the analysis. This half-width, measured in periods, is entered
c    as the second parameter, m, of the 'ban' command. m=3 is generally
c    sufficient; higher values result in more selective filters, thus
c    a better immunity to noise. Each test frequency must be present
c    for at least 5*m periods at constant amplitude (more if tapered).
c    amplitude (more if tapered). Between the test signals, quiet
c    intervals corresponding to m periods of the slower one of the
c    adjacent signals should be present.
c
c    The program is controlled by a parameter file 'sincal,par'
c    unless a different file name is specified in the runstring.
c    This file contains file names and processing commands. (see below).
c    Several data files can be processed in one run if the test
c    frequencies are the same.
c
c
c    *******************************************************************
c
c
c   Derived from program SEIFE written by e. wielandt, eth zurich, 1984
c   sinewave amplitudes inverted for seismometer parameters Oct. 2002
c   Updated July 2004.
c
c   this program has limited time-domain signal processing capabilities
c   like SEIFE, but is not a general data-processing tool, and produces
c   no output signal.
c
c   CAUTION! Since this program is a modification of SEIFE, some comments
c   may apply to the original program rather than the modified one. I have
c   however tried to update the comments as well.
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   Other supported formats are:
c   - single-column ASCII with header lines ignored (read only)
c   - MARS88 - ASCII with header line (read only)
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 files (up to 24; all
c   processes apply to all input files). each line of the command file
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 columns 
c   5 to 50.
c
c   the following codes are accepted:
c
c   blank first character, or rem: line is ignored     
c
c   commands for input and output:
c
c   fil  nam1: specify input files 
c        fil will read 'seife' files.
c   mar  nam1  like fil, but for MARS88-ASCII files from MARSDUMP
c   asc  nam1 dt lin    where dt is the sampling interval and lin the
c        number of header lines to be skipped. Will read any ascii data.
c        The last character of the last line must be a <cr>, or the last
c        sample will be skipped.        
c
c   lim  n: limit the length of the input series to n samples
c
c   end: end of input file, start the inversion for system parameters
c
c   Commands for signal processing:
c
c   ban  t0 m: gaussian FIR bandpass-filter with center period t0
c        and a total width of 2*m periods. The maximum amplitude of
c        the output signal is stored but the output signal itself
c        is discarded.
c   avg  n: remove average (determined from first n samples)
c   dec  n: decimate to every n'th sample
c   dif  differentiate the signal
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   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   rev  reverse signal in time (useful for backward filtering)
c   rms  determine rms amplitude (signal is unchanged)
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   tap  n1  n2: signal is cosine- tapered between samples n1 and n2
c        (signal is unchanged between samples 1 and n1)
c   tre  n: remove linear trend (determined from first n samples)
c   twi  tmin1  tsec1  tmin2  tsec2: window defined by time after midnight
c   swi  sec1 sec2: window defined by seconds after first sample
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
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]
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 Additional type: sez, inverse second-order high-pass filter (restitution)
c   sez t0 h: remove seismo response with free period t0 and damping h
c
c end of sincal1 info.


      program sincal1
c  the following statement is required with NAS-FortranPlus
      use nas_system
      parameter(ndim=120000,ndf=ndim/12)
      dimension x(ndim),f1(ndf),f2(ndf)
      dimension ix(ndim),fdt(32),lin(32),t0(9),resp(9)
      equivalence (x,ix)
      character typ*3,text*72,ttext*72,zeile*60,par*56,info*78
      character*24 nam1(32),nam3
      logical nil
      data fdt/32*0./, lin/32*0/
        write(*,*)
        write(*,*) 'HELP: >> sincal1 info <<  or  >> sincal1 -help <<'
        write(*,*)
      if(iargc().eq.0) then
        nam3='sincal.par'
      else
        call getarg(1,nam3)
      endif
      if(nam3(1:4).eq.'info'.or.nam3(1:5).eq.'-help') then
        open(7,file='sincal1.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)') nam3
      open(7,file=nam3,status='old')
      nmax=ndim
      nf=ndf
      call files(nam1,fdt,lin,nseis,zeile,nmax)
      close(7)
      do 1 ns=1,nseis
      n1=1
      dt=fdt(ns)
      n=ndim+1-n1
       n=min(n,nmax)
       call inpasc(nam1(ns),lin(ns),ttext,n,dt,ttmin,ttsec,x(n1))
      if(n1.eq.1) then
        tmin=ttmin
        tsec=ttsec
        text=ttext
      endif
      if(n.eq.0.or.dt.eq.0.) then
      write(6,'("n or dt unspecified - check data header")')
      stop
      endif
      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)
      par=zeile(5:60)
      if(typ.eq.'fil'.or.typ.eq.'asc'.or.typ.eq.'lim'
     & .or.typ.eq.'mar') goto 2
      if(typ.eq.'prt') then
        write(6,*) par(1:56)
        goto 2
      endif
      if(typ.eq.'rem') goto 2
      nil=.true.
      n=min(n,nmax)
      if(typ.eq.'sei') then
        ic=0
        read(par,*) ts,hs
        nil=.false.
      endif
      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.'ban') call band(nil,par,x,n,dt,f1,f2,nf,t0,resp,ic)
      if(typ.eq.'dif') call deriv(nil,par,x,n,dt)
      if(typ.eq.'spl') call interp(nil,par,x,n)
      if(typ.eq.'win'.or.typ.eq.'skp'.or.typ.eq.'tap'.or.typ.eq.'swi') 
     & call window(nil,typ,par,x,n,dt,tmin,tsec)
      if(typ.eq.'twi') call timwin(nil,par,x,n,dt,tmin,tsec)
      if(typ.eq.'nul') call zero(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.'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.'ext') call extrem(nil,x,n)
      if(typ.eq.'rms') call rms(nil,x,n)
      if(typ.eq.'pad') call pad(nil,par,x,n)
      if(nil) call filter(nil,zeile,x,n,dt,typ,par)
c  all done - store result
      if(typ.eq.'end') then
        call calib(par,ts,hs,t0,resp,ic)
        goto 1
      endif
c  produce error message if code unrecognized
      if(nil) then
        write(6,3) typ
    3   format('prozedur ',a3,' ist nicht definiert --- stop')
        stop
      else
        goto 2
      endif
    1 continue
      stop
    4 write(6,'("end-zeile fehlt im befehlsfile --- stop")')
      stop
      end

c
      subroutine calib(par,ts,hs,t00,resp0,ic0)
      dimension t00(9),t0(9),resp0(9),resp(9),tresp(9)
      dimension rest0(9),rest1(9),rest2(9),tb(9),hb(9),gb(9)
      character par*56
      logical v

c invert data (ib=0) and then perform bootstrap test
    1 format("Checkpoint ",i2)
      do ib=0,ic0
        if(ib.eq.0) then
          v=.true.
          ic=ic0
          do i=1,ic
            t0(i)=t00(i)
            resp(i)=resp0(i)
          enddo
        else
          v=.false.
          if(ic0.lt.4) then
            write(*,*)
            write(*,*) '*** insufficient number of test frequencies ***'
            stop
          endif
          if(ib.eq.1) then
            write(6,*)
            write(6,'(a)') ' Bootstrap test:'
          endif
          ic=ic0-1
          iadd=0
          do i=1,ic
            if(i.eq.ib) iadd=1
            t0(i)=t00(i+iadd)
            resp(i)=resp0(i+iadd)
          enddo
        endif

      if(v) write(6,*)
      if(v) write(6,'(a)') ' Fitting a resonance curve to the data:'
      do iter=1,3
      gain=0.
      do i=1,ic
        tau=t0(i)/ts
        tresp(i)=sqrt((tau-1./tau)**2+4.*hs**2)
        gain=gain+resp(i)*tresp(i)
      enddo
      gain=gain/ic
      err=0.
      do i=1,ic
        tresp(i)=gain/tresp(i)
        rest0(i)=resp(i)-tresp(i)
        err=err+rest0(i)**2
      enddo
      err=sqrt(err/ic)
      if(v.and.iter.eq.1) write(6,111) " ts=",ts,"   hs=",hs,
     &  "   gain=",gain,"   rms=",err
  111 format(a,f8.3,a,f7.4,a,f10.3,a,f8.3)
      gain=0.
      do i=1,ic
        tau=t0(i)/ts/1.001
        tresp(i)=sqrt((tau-1./tau)**2+4.*hs**2)
        gain=gain+resp(i)*tresp(i)
      enddo
      gain=gain/ic
      do i=1,ic
        tresp(i)=gain/tresp(i)
        rest1(i)=resp(i)-tresp(i)-rest0(i)
      enddo
      gain=0.
      do i=1,ic
        tau=t0(i)/ts
        tresp(i)=sqrt((tau-1./tau)**2+4.008*hs**2)
        gain=gain+resp(i)*tresp(i)
      enddo
      gain=gain/ic
      err=0.
      do i=1,ic
        tresp(i)=gain/tresp(i)
        rest2(i)=resp(i)-tresp(i)-rest0(i)
      enddo
      c11=0.
      c12=0.
      c22=0.
      b1=0.
      b2=0.
      do i=1,ic
        c11=c11+rest1(i)*rest1(i)
        c12=c12+rest1(i)*rest2(i)
        c22=c22+rest2(i)*rest2(i)
        b1=b1+rest1(i)*rest0(i)
        b2=b2+rest2(i)*rest0(i)
      enddo
      c21=c12
      det=c11*c22-c21*c12
      x1=(b1*c22-b2*c12)/det
      x2=(c11*b2-c21*b1)/det
      ts=ts*(1-0.001*x1)
      hs=hs*(1-0.001*x2)
      gain=0.
      do i=1,ic
        tau=t0(i)/ts
        tresp(i)=sqrt((tau-1./tau)**2+4.*hs**2)
        gain=gain+resp(i)*tresp(i)
      enddo
      gain=gain/ic
      err=0.
      do i=1,ic
        tresp(i)=gain/tresp(i)
        rest0(i)=resp(i)-tresp(i)
        err=err+rest0(i)**2
      enddo
      err=sqrt(err/ic)
      if(v) write(6,111) " ts=",ts,"   hs=",hs,"   gain=",gain,
     &  "   rms=",err
      enddo
      if(v) then
        ts0=ts
        hs0=hs
        gs0=gain
        write(6,*)
        write(6,'(a)') ' Final fitting errors:'
        do i=1,ic
          write(6,112) " t0=",t0(i),"   resp=",resp(i),"   theor=",
     &    tresp(i),"   err=",resp(i)-tresp(i)
  112     format(a,f10.3,a,f10.3,a,f10.3,a,f8.3)
        enddo
      endif
      if(.not.v) then
      write(6,113) " w/o t0=",t00(ib),": ts=",ts,
     &  "  hs=",hs,"  gain=",gain,"  rms=",err
  113 format(a,f8.3,a,f9.3,a,f7.4,a,f10.3,a,f8.3)
      tb(ib)=ts
      hb(ib)=hs
      gb(ib)=gain
      endif
      enddo

c evaluate bootstrap test
      tba=0.
      hba=0.
      gba=0.
      tbaq=0.
      hbaq=0.
      gbaq=0.
      do i=1,ic0
        tba=tba+tb(i)
        hba=hba+hb(i)
        gba=gba+gb(i)
      enddo
      ts1=tba/ic0
      hs1=hba/ic0
      gs1=gba/ic0
      trms=0.
      hrms=0.
      grms=0.
      do i=1,ic0
        trms=trms+(tb(i)-ts1)**2
        hrms=hrms+(hb(i)-hs1)**2
        grms=grms+(gb(i)-gs1)**2
      enddo
      ic1=ic0-1
      trms=sqrt(trms/ic1)
      hrms=sqrt(hrms/ic1)
      grms=sqrt(grms/ic1)
      write(*,*)
      write(*,*)   'RESULT:   best fit   boots.avg.    std.dev.'
      write(*,114) ' Period  ',ts0,ts1,'  +-',trms
      write(*,115) ' Damping ',hs0,hs1,'  +-',hrms
      write(*,116) ' Gain    ',gs0,gs1,'  +-',grms
  114 format(a,f10.3,f12.3,a,f8.3)
  115 format(a,f10.4,f12.4,a,f8.4)
  116 format(a,f10.2,f12.2,a,f8.2,'   normalized with nox?')
      return
      end

c
      subroutine inpasc(name,lin,text,n,dt,tmin,tsec,x)
      dimension x(n)
      character iform*20,name*24,text*72,pro*80
      write(6,*) '================================================'
      write(6,*) 'Opening file ',name
      open(7,file=name,status='old')
c  ASCII files with header
      namlen=index(name,' ')
      if(dt.gt.0.) then
        write(6,22) name(1:namlen),0,n,dt
        do i=1,lin
          read(7,'(a)') text
        enddo
        text=name

c NAS FortranPlus
        read(7,*,err=30,end=30) (x(j),j=1,n)
c Other compilers
c       read(7,*,err=31,end=31) (x(j),j=j,n)

c FortranPlus does not recognize the end of data and fills
c the series with zeros. Therefore, we remove trailing zeros.
   30   continue
        do j=n,1,-1
          if(x(j).ne.0.) exit
        enddo
        if(j.eq.0) then
          write(6,*) 'All samples are zero - stop'
          stop
       endif
       n=j
       write(6,32) n
       close(7)
       return

   31 continue
      n=j
      write(6,32) n
   32 format(' end of data after',i8,' samples')
      close(7)
      return

c  SEIFE format
      else if(dt.eq.0.) then
        read(7,'(a)') text
        write(6,'(1x,a)') text
   21   read(7,'(a)') pro
        if(pro(1:1).ne.'%'.and.pro(1:1).ne.'}') goto 20
        goto 21
   20   read(pro,1) nn,iform,dt,tmin,tsec
    1   format(i10,a20,3f10.3)
        write(6,22) name(1:namlen),nn,n,dt
   22   format(' reading file ',a,' ',i7,'/',i7,' 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') then
          read(7,*,err=25,end=23) (x(j),j=1,n)
        else
          read(7,iform,end=23) (x(j),j=1,n)
        endif
c  MARS88-ASCII
      else if(dt.lt.0.) 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) name,n,n,dt
        do j=1,n
          read(7,*,err=23,end=23) x(j)
        enddo
      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) name,j
   26 format(' Input error reading file ',a,' at sample # ',i7)
      stop
      end
c



      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,fdt,lin,nseis,zeile,npts)
      character nam1(32)*24,zeile*60,z3*3
      dimension fdt(32),lin(32)
      nseis=0
  101 read(7,'(a)',end=2) zeile
      write(6,'(a)') zeile(1:60)
      z3=zeile(1:3)
      if(z3.eq.'fil'.or.z3.eq.'asc'.or.z3.eq.'mar') then
        nseis=nseis+1
        if(z3.eq.'mar') then
            fdt(nseis)=-1.
            lin(nseis)=1
        endif
        do 30 i=5,60
   30     if (zeile(i:i).ne.' '.and.zeile(i:i).ne.',') goto 31  
   31   ii=i
        do 11 i=ii,60
   11     if (zeile(i:i).eq.' '.or.zeile(i:i).eq.',') goto 34
   34   if(i.gt.ii+24) goto 24
        nam1(nseis)=zeile(ii:i-1)
        if(z3.eq.'asc') then
        ii=i
          do 14 i=ii,60
   14       if (zeile(i:i).ne.' '.and.zeile(i:i).ne.',') goto 4
    4     if(i.gt.58) goto 26
          if(z3.eq.'asc') then
            read(zeile(i:60),*,err=26) fdt(nseis),lin(nseis)
          else
            read(zeile(i:60),*,err=26) lin(nseis)
            lin(nseis)=-lin(nseis)
          endif
        endif
      else if(z3.eq.'lim') then
        read(zeile,'(4x,i10)') npts
      else if(z3.eq.'end') then
        if(nseis.gt.0) then
          return
        else
          write(6,3)
    3     format(' Kein Filename angegeben - stop')
          stop
        endif
      endif
      goto 101
    2 write(6,'("end-zeile fehlt im befehlsfile --- stop")')
      stop
   24 write(6,25) zeile
   25 format(' Filename zu lang: ',a)
      stop
   26 write(6,27)
   27 format(' Anweisung asc nicht lesbar - stop')
      stop
      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  ",f10.3)') f
      nil=.false.
      return
   99 write(6,*) ' insufficient or wrong parameters for fac'
      stop
      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,a,a,e12.6)') typ,': normalization factor ',f
      nil=.false.
      return
   99 write(6,*) ' insufficient or wrong parameters for ',typ
      stop
      end

      subroutine pad(nil,par,x,n)
c  Nullen ergaenzen
      dimension x(n)
      character*56 par
      logical nil
      read(par,*,err=99,end=99) nul
      do 1 j=n+1,nul
    1 x(j)=0.
      n=nul
      write(6,2) nul
    2 format(" pad: Es wurden Nullen aufgefuellt bis Punkt",i8)
      nil=.false.
      return
   99 write(6,*) ' insufficient or wrong parameters for pad'
      stop
      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
      do 1 j=n1,n2
    1 x(j)=0.
      write(6,'(a,2i10)') ' null ',n1,n2
      nil=.false.
      return
   99 write(6,*) ' insufficient or wrong parameters for null'
      stop
      end

      subroutine deriv(nil,par,x,n,dt)
c  derivative (non-recursive, non-causal, symmetric-difference)
      dimension x(n)
      character*56 par
      logical nil
      read(par,*,err=99,end=99) tau
      if(tau.eq.0.) tau=1.
      twodt=2.*dt/tau
      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
   99 write(6,*) ' insufficient or wrong parameters for dif'
      stop
      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 - an",i4," stellen wurden",i5," werte >",f9.0,
     & " interpoliert. 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
      write(6,'(" *** spl ***  zu wenig stuetzwerte - stop")')
      stop
      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 write(6,*) ' insufficient or wrong parameters for spl'
      stop
      end

      subroutine window(nil,typ,par,x,n,dt,tmin,tsec)
c  time window by sample number
      dimension x(n)
      character*56 par
      character*3 typ
      logical nil
      if(typ.eq.'skp') 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
      if(typ.ne.'tap') goto 16
      n=np2
      pih=2.*atan(1.)
      fak=pih/(np2-np1)
      do 14 j=np1,np2
   14 x(j)=x(j)*cos(fak*(j-np1))**2
      goto 15
   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
   13 if(typ.eq.'win'.or.typ.eq.'swi') then
      do 11 j=1,n
   11 x(j)=x(np1+j)
      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(6,*) ' 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 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
      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 write(6,*) ' insufficient or wrong parameters for twi'
      stop
      end

      subroutine decim(nil,par,x,n,dt)
c  decimate time series by a factor up to 100
      dimension x(n),g(100)
      character*56 par
      logical nil
      read(par,*,err=99,end=99) idez
      if(idez.gt.100) stop 100
      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,'(" dec  ",i10)') idez
      nil=.false.
      return
   99 write(6,*) ' insufficient or wrong parameters for dec'
      stop
      end

      subroutine mean(nil,par,x,n)
c  remove average
      character*56 par
      dimension x(n)
      logical nil
      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  ",e11.3)') sum
      nil=.false.
      return
   99 write(6,*) ' insufficient or wrong parameters for avg'
      stop
      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
        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  ",2f10.3)') a,b
      nil=.false.
      return
   99 write(6,*) ' insufficient or wrong parameters for ',typ
      stop
      end

      subroutine polytrend(nil,par,x,n)
c  remove polynomial trend
      parameter(ndi=6)
      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(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 write(6,*) ' insufficient or wrong parameters for pol'
      stop
      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 write(6,*) ' insufficient or wrong parameters for exp'
      stop
      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
c
      subroutine filter(nil,zeile,x,n,dts,typ,par)
      implicit real*8 (a-h,o-z)
      character typ*3,typ1*3,typ2*3,par*56,zeile*60
      real*4 x(n),dts
      logical nil
      dt=dts
      typ=zeile(1:3)
      par=zeile(5:60)
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
      it=11
      else if(typ.eq.'hp1') then
      read(par,*,err=99,end=99) t0
      it=12
      else if(typ.eq.'lp2') then
      read(par,*,err=99,end=99) t0,h
      it=21
      else if(typ.eq.'hp2') then
      read(par,*,err=99,end=99) t0,h
      it=22
      else if(typ.eq.'bp2') then
      read(par,*,err=99,end=99) 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,*,err=99,end=99) t0s,t0
      it=13
      else if(typ.eq.'he2'.or.typ.eq.'le2') then
        read(par,*,err=99,end=99) t0s,hs,t0,h
      it=23
      else if(typ.eq.'sez') then
        read(par,*,err=99,end=99) t0s,hs
        t0=1.e12
        h=1.
        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
      call rekfl(x,x,n,f0,f1,f2,g1,g2)
    8 format(1x,a3,2x,5f10.3)
      nil=.false.
c  confirm execution
      if(it.eq.1) write(6,8) typ,t0
      if(it.eq.11.or.it.eq.12) write(6,8) typ,t0
      if(it.eq.13) write(6,9) typ,t0s,t0
    9 format(1x,a3,2x,f10.3,f16.3)
      if(it.eq.21.or.it.eq.22.or.it.eq.24) write(6,8) typ,t0,h
      if(it.eq.23) then
        if(typ.eq.'sez') then
          write(6,8) typ,t0s,hs
        else
          write(6,7) typ,t0s,hs,t0,h
    7 format(1x,a3,2x,2f10.3,f16.3,f10.3)
        endif
      endif
      return

c  Butterworth filters
    2 if(typ.eq.'lpb'.or.typ.eq.'hpb') then
        read(par,*,err=99,end=99) t0,m
        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
          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,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)
          call rekfl(x,x,n,f0,f1,f2,g1,g2)
          write(6,8) typ2,t0,h
    3   continue
      nil=.false.
      endif
      return
   99 write(6,*) ' insufficient or wrong parameters for ',typ
      stop
      end

      subroutine extrem(nil,x,n)
      dimension x(n)
      logical nil
      nil=.false.
      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
    2 format(" ext   xmin=",e11.3," (j=",i6,")   xmax=",e11.3," (j=",i6,
     &")")
      return
      end

      subroutine rms(nil,x,n)
      dimension x(n)
      logical nil
      nil=.false.
      xrms=0.
      do 1 j=1,n
    1 xrms=xrms+x(j)**2
      write(6,2) sqrt(xrms/n)
    2 format(" rms  =",e11.3)
      return
      end

      subroutine band(nil,par,x,n,dt,f1,f2,nf,t0,resp,ic)
      character*56 par
      logical nil
      dimension x(n),f1(nf),f2(nf),t0(9),resp(9)
      ic=ic+1
      if(ic.eq.1) write(6,'(/,a)')
     & ' Bandpass filtration and amplitude measurement:'
      if(ic.gt.9) then
        write(*,*) "This program accepts only 9 test frequencies"
        stop
      endif
      read(par,*,err=99,end=99) t0(ic),m
      lenh=int(m*t0(ic)/dt)
      if(lenh.gt.nf) then
        write(6,*) 'ban:  FIR filter is too long - stop'
        stop
      endif
      twopi=8.*atan(1.)
      om=twopi/t0(ic)
      fn1=1.
      fn2=0.
      do k=1,lenh
        co=cos(om*k*dt)
        si=sin(om*k*dt)
        f1(k)=co*exp(-(2.5*k/lenh)**2)
        f2(k)=si*exp(-(2.5*k/lenh)**2)
        fn1=fn1+2.*f1(k)*co
        fn2=fn2+2.*f2(k)*si
      enddo
      do k=1,lenh
        f1(k)=f1(k)/fn1
        f2(k)=f2(k)/fn2
      enddo
      xmax=0.
      do i=1,n
        yy=x(i)/fn1
        zz=0.
        do k=1,lenh
          if(i.gt.k) then
            yy=yy+f1(k)*x(i-k)
            zz=zz-f2(k)*x(i-k)
          endif
          if(i+k.le.n) then
            yy=yy+f1(k)*x(i+k)
            zz=zz+f2(k)*x(i+k)
          endif
        enddo
        xmax=max(xmax,yy**2+zz**2)
      enddo
      xmax=sqrt(xmax)
      resp(ic)=xmax
      nil=.false.
      write(6,90) " filter t0=",t0(ic),"   m=",m,"   resp=",resp(ic)
   90 format(a,f8.3,a,i2,a,f10.3)
      return
   99 write(6,*) ' insufficient or wrong parameters for ban'
      stop
      end



