      complex function butter(f,fl,fh,npole)
c
c ...... bandpass filter  (nPBP Butterworth Filter)
c
c        convolve with npole Butterworth Bandpass filter
c
c        stolen from Seisgram - 10/28/92
c        - 6 pole output compared with subroutine bpbw
c          Pretty good, although slight differences exist
c
c        where -
c          f     - frequency of evaluation in rad/s
c          fl    - low frequency corner in rad/s
c          fh    - high frequency corner in rad/s
c          npole - number of poles in filter at each corner
c                  (not more than 20)
c
      integer*4 npole
      complex*8 sph(20),spl(20)
      complex*8 c0, c1,cjw
      complex*8 cph,cpl
c
c     include some useful stuff
c
      pi = 4.*atan(1.0)
c
      nop = npole-2*(npole/2)
      nepp = npole/2
      c0 = cmplx(0.,0.)
      c1 = cmplx(1.,0.)
      w = f
      wch = fh
      wcl = fl
      np = 0
      if (w .eq. 0.0) then
        butter = c0
        return
      end if
      if (nop.gt.0) then
        np = np+1
        sph(np) = c1
      endif
      if (nepp.gt.0) then
        do i = 1, nepp
          ak = 2.*sin((2.*float(i)-1.)*pi/(2.*float(npole)))
          ar = 0.5*ak*wch
          ai = 0.5*wch*sqrt(4.-ak*ak)
          np = np+1
          sph(np) = cmplx(-ar,-ai)
          np = np+1
          sph(np) = cmplx(-ar,+ai)
        end do
      endif
      np = 0
      if (nop.gt.0) then
        np = np+1
        spl(np) = c1
      endif
      if (nepp.gt.0) then
        do i = 1, npole/2
          ak = 2.*sin((2.*float(i)-1.)*pi/(2.*float(npole)))
          ar = 0.5*ak*wcl
          ai = 0.5*wcl*sqrt(4.-ak*ak)
          np = np+1
          spl(np) = cmplx(-ar,-ai)
          np = np+1
          spl(np) = cmplx(-ar,+ai)
        end do
      endif
      cjw = cmplx(0.,-w)
      cph = c1
      cpl = c1
      do j = 1, npole
        cph = cph*sph(j)/(sph(j)+cjw)
        cpl = cpl*cjw/(spl(j)+cjw)
      end do
      butter = cph*cpl
c
      return
      end
