	subroutine coinst2(npts,dt,ind,ind2,c1,cw,np,fl,fh,ierr)
c
c     subroutine for Fourier Domain operations in process_waveform
c
c     input:  
c             npts   - number of points in timeseries
c             dt     - sample spacing in sec
c             ind    - type of instrument operation to perform
c                 1 - Benioff 100 kg
c                 2 - 17s-23s 10 pole Butterworth
c                 3 - Wood Anderson
c                 4 - ULP instrument
c                 5 - SP WWSSN instrument
c                 6 - LP WWSSN instrument
c                 7 - LP HGLP instrument
c                 8 - Nominal STS1 response
c                 9 - Benioff 14kg 
c                10 - ASRO LP response
c                11 - STS2 LP response
c                20 - deconvolve to displacement
c                21 - deconvolve to velocity
c                22 - deconvolve to acceleration
c             ind2  - type of filter operation to perform
c                 1 - apply 6 pole Butterworth bandpass filter between fl and fh
c                 2 - apply 6 pole Butterworth phaseless bandpass filter between fl and fh
c                 3 - apply 6 pole Butterworth lowpass  filter with cutoff fl
c                 4 - apply 6 pole Butterworth phaseless lowpass  filter with cutoff fl
c                 5 - apply 6 pole Butterworth hipass   filter with cutoff fh
c                 6 - apply 6 pole Butterworth phaseless hipass   filter with cutoff fh
c                 7 - apply np pole Butterworth bandpass filter between fl and fh
c             c1    - input and output timeseries
c             cw    - input and output dummy array for storing spectra
c             np    - input number of poles of option 7
c             fl    - input low  freq corner of Butterworth filters (input in Hz)
c             fh    - input high freq corner of Butterworth filters (input in Hz)
c             ierr  - error return code
c
      real*4 c1(*)
c
	real*4 raw
      complex*8 cw(*), aw, rw, response
      complex*8 ben100, lpbptf_b, wood
      complex*8 ulp, wwsp, wwlp
      complex*8 hglp, sts1, ben14, asro, sts2
      complex*8 lobwth, hibwth, bpbwth, butter
c
      integer*4 ipow
c
c     compute useful numbers
c
      include 'numerical.h'
c
      io = 6
c
c     initialize error flag
c
      ierr = 0
c
c     detrend/deman
c
c     call detrend(c1, npts, dt, c1)
c
      call demean(c1, npts)
c
c     taper
c
c     call taper(c1, npts, 10.0)
c
      call taper(c1, npts, 5.0)
c
c     compute the number of points
c
c      Seistool option
c        n2 = ipow(npts)
c
       n2 = ipow(npts*2)
c
c     frequency increment & number of frequency samples
c
      dw = tpi / (float(n2) * dt)
      nt2h = (n2 + 2) / 2
c
c     zero excess points
c
      do ii = npts + 1,  n2
        c1(ii) = 0.0
      end do
c
c     forward transform
c
      ierr = 0
      call fftl (c1, n2, -1, dt, ierr)
      if (ierr .ne. 0) then
        ierr = 2
        write(io,*) 'Error in forward transform'
        return
      endif
c
c     form complex array
c
      do ii = 1, nt2h
        jj = 2*(ii - 1) + 1
        cw(ii) = cmplx(c1(jj),c1(jj+1))
      end do
c
c     Perform instrument operations
c
      itype = 1
      if (ind .eq. 0) then
c       write(io,*) 'No Instrument operations'
      elseif (ind .eq. 1) then
c
c      convolve with the Benioff 100kg response
c
        cw(1) = (0.0,0.0)
        do ii = 2, nt2h
          f = float(ii - 1)*dw
          rw = response(f, itype)
          aw = ben100(f)    
          cw(ii) = cw(ii) * (aw/rw)
        end do
      elseif (ind .eq. 2) then
c 
c      convolve with a 20 s filter
c 
        cw(1) = (0.0,0.0)
        do ii = 2, nt2h 
          f = float(ii - 1)*dw
          rw = response(f, itype)
          aw = lpbptf_b(f)
          cw(ii) = cw(ii)*(aw/rw)
        end do
      elseif (ind .eq. 3) then
c 
c      convolve with the Wood Anderson response
c 
        cw(1) = (0.0,0.0)
        do ii = 2, nt2h 
          f = float(ii - 1)*dw
          rw = response(f, itype)
          aw = wood(f)    
          cw(ii) = cw(ii)*(aw/rw)
        end do
      elseif (ind .eq. 4) then
c
c      convolve with the ULP response
c
        cw(1) = (0.0,0.0)
        do ii = 2, nt2h
          f = float(ii - 1)*dw
          rw = response(f, itype)
          aw = ulp(f)
          cw(ii) = cw(ii)*(aw/rw)
        end do
      elseif (ind .eq. 5) then
c
c      convolve with the SP WWSSN response
c
        cw(1) = (0.0,0.0)
        do ii = 2, nt2h
          f = float(ii - 1)*dw
          rw = response(f, itype)
          aw = wwsp(f)
          cw(ii) = cw(ii)*(aw/rw)
        end do
      elseif (ind .eq. 6) then
c
c      convolve with the LP WWSSN response
c
        cw(1) = (0.0,0.0)
        do ii = 2, nt2h
          f = float(ii - 1)*dw
          rw = response(f, itype)
          aw = wwlp(f)
          cw(ii) = cw(ii)*(aw/rw)
        end do
      elseif (ind .eq. 7) then
c
c      convolve with the LP HGLP response
c
        cw(1) = (0.0,0.0)
        do ii = 2, nt2h
          f = float(ii - 1)*dw
          rw = response(f, itype)
          aw = hglp(f)
          cw(ii) = cw(ii)*(aw/rw)
        end do
      elseif (ind .eq. 8) then
c
c      convolve with the nominal STS1 response
c
        cw(1) = (0.0,0.0)
        do ii = 2, nt2h
          f = float(ii - 1)*dw
          rw = response(f, itype)
          aw = sts1(f)
          cw(ii) = cw(ii)*(aw/rw)
        end do
      elseif (ind .eq. 9) then
c
c      convolve with a Benioff 14 kg
c
        cw(1) = (0.0,0.0) 
        do ii = 2, nt2h 
          f = float(ii - 1)*dw 
          rw = response(f, itype) 
          aw = ben14(f)    
          cw(ii) = cw(ii)*(aw/rw) 
        end do 
      elseif (ind .eq. 10) then
c          
c      convolve with a ASRO LP response
c
        cw(1) = (0.0,0.0)
        do ii = 2, nt2h
          f = float(ii - 1)*dw
          rw = response(f, itype)
          aw = asro(f)
          cw(ii) = cw(ii)*(aw/rw)
        end do
      elseif (ind .eq. 11) then
c
c      convolve with a STS2 LP response
c
        cw(1) = (0.0,0.0)
        do ii = 2, nt2h
          f = float(ii - 1)*dw
          rw = response(f, itype)
          aw = sts2(f)
          cw(ii) = cw(ii)*(aw/rw)
        end do
      elseif (ind .gt. 11 .and. ind .lt. 20) then
        write(io,*) 'Non-existent response'
        ierr = 1
        return
c
c     remove the instrument only - no reconvolution
c
      elseif (ind .eq. 20) then
c
c      remove instrument to displacement
c
        itype = 1
        cw(1) = (0.0,0.0)
        do ii = 2, nt2h
          f = float(ii - 1)*dw
          rw = response(f, itype)
          cw(ii) = cw(ii) / rw
        end do
      elseif (ind .eq. 21) then 
c 
c      remove instrument to velocity
c 
        itype = 2 
        cw(1) = (0.0,0.0) 
        do ii = 2, nt2h 
          f = float(ii - 1)*dw 
          rw = response(f, itype) 
          cw(ii) = cw(ii) / rw 
        end do 
      elseif (ind .eq. 22) then 
c 
c      remove instrument to acceleration
c 
        itype = 3 
        cw(1) = (0.0,0.0) 
        do ii = 2, nt2h 
          f = float(ii - 1)*dw 
          rw = response(f, itype) 
          cw(ii) = cw(ii) / rw 
        end do 
      elseif (ind .ge. 23) then
        write(io,*) 'Unknown option'
        ierr = 1
        return
      end if
c
c     Check for filter operations
c
      wl = fl*tpi
      wh = fh*tpi
      if (ind2 .eq. 0) then
c       write(io,*) 'No filtering'
      elseif (ind2 .eq. 1) then
c
c     apply a 6 pole bandpass Butterworth filter
c
        do ii = 1, nt2h
          f = float(ii - 1)*dw
          cw(ii) = cw(ii)*bpbwth(f, wl, wh)
        end do
      elseif (ind2 .eq. 2) then
c
c     apply a 6 pole bandpass phaseless Butterworth filter
c
        do ii = 1, nt2h
          f = float(ii - 1)*dw
          aw = bpbwth (f, wl, wh)
          raw = aw * conjg(aw)
          cw(ii) = cw(ii)*sqrt(raw)
        end do
      elseif (ind2 .eq. 3) then
c
c     apply a 6 pole lowpass Butterworth filter
c
        do ii = 1, nt2h
          f = float(ii - 1)*dw
          cw(ii) = cw(ii)*lobwth(f, wl)
        end do
      elseif (ind2 .eq. 4) then
c
c     apply a 6 pole lowpass phaseless Butterworth filter
c
        do ii = 1, nt2h  
          f = float(ii - 1)*dw
          aw = lobwth (f, wl)
          raw = aw * conjg(aw)
          cw(ii) = cw(ii)*sqrt(raw)
        end do    
      elseif (ind2 .eq. 5) then
c
c     apply a 6 pole hipass Butterworth filter
c
        do ii = 1, nt2h
          f = float(ii - 1)*dw
          cw(ii) = cw(ii)*hibwth(f, wh)
        end do
      elseif (ind2 .eq. 6) then  
c
c     apply a 6 pole hipass phaseless Butterworth filter
c
        do ii = 1, nt2h
          f = float(ii - 1)*dw
          aw = hibwth (f, wh)
          raw = aw * conjg(aw)
          cw(ii) = cw(ii)*sqrt(raw)
        end do
      elseif (ind2 .eq. 7) then
c
c     apply a np pole bandpass Butterworth filter
c
        do ii = 1, nt2h
          f = float(ii - 1)*dw
          cw(ii) = cw(ii)*butter(f,wl,wh,np)
        end do
c
      elseif (ind2 .ge. 8) then
        write(io,*) 'Unknown option'
        ierr = 1
        return
      end if
c
c     reform the complex array into a real
c
      do ii = 1, nt2h 
        jj = 2*(ii - 1) + 1 
        c1(jj) = real(cw(ii))
        c1(jj+1) = aimag(cw(ii))
      end do 
c
c     inverse transform filtered function
c
      ierr = 0
      call fftl (c1, n2, -2, dt, ierr)
      if (ierr .ne. 0) then
        ierr = 2
        write(io,*) 'Error in forward transform'
        return
      endif
c
c     all done
c
      return
      end
c
