      subroutine bessel_fcoeff (n,b)
c
c Given filter order N, computes the Bessel-N filter coefficients B(0..N)
c
      implicit real*8 (a-h,o-z)
      dimension b(0:*)
      logical qcall,qshow,qtrace
      common /ccbess/ qcall,qshow,qtrace
      
      b(n)=1.d0
      do k=n,1,-1
         b(k-1)=b(k)*(2*n-k+1.d0)*k/(dfloat(2*(n-k+1)))
      enddo
      qcall=.true.
      if (qshow) then
         print '(a,1p,10d12.4,:/(12x,1p,10d12.4))'
     &,' <BessFc>>> ', (b(k),k=0,n)
         qshow=.false.
      endif
      return
      end

      complex*16 function cbessel_f (omega,n,b)
c
c Given filter coefficients B(0..N) and the ratio Omega = f/f0
c where f0 is the "corner frequency", computes the complex spectral
c filter response at frequency f
c
c Source: http://en.wikipedia.org/wiki/Bessel_filter
c
c omega = 2pi f / 2pi f0 where f0 is the corner frequency
c (which is not a -3dB corner). E.g. n=2 => |H(1)| = 3/sqrt(13)
c
      implicit complex*16 (z), real*8 (a-h,o-y)
      dimension b(0:*)
      logical qcall,qshow,qtrace
      common /ccbess/ qcall,qshow,qtrace

      data zi/(0.0d0,1.0d0)/

      if (.not.qcall) then
         call bessel_fcoeff (n,b)
      endif
      
      zr=b(0)
      
      zp=b(n)
      do k=1,n
         zp=zp*zi*omega+b(n-k)
      enddo

      cbessel_f = zr/zp
      if (qtrace) print *,'<CbessF>t>', omega,zr/zp
      return
      
      end

      subroutine reset_bessel_f ()
      logical qcall,qshow,qtrace
      common /ccbess/ qcall,qshow,qtrace
      qcall=.false.
      return

      entry show_bessel_f ()
      qshow=.true.
      return

      entry trace_bessel_f ()
      qtrace=.true.
      return
      end
      
      block data bcbess
      logical qcall,qshow,qtrace
      common /ccbess/ qcall,qshow,qtrace
      data qcall,qshow,qtrace/3*.false./
      end

