c-------------------------------------------------
c     Author:  Brian Martin
c     modified last:  Dec. 22, 2004.
c-------------------------------------------------

c-------------------------------------------------
c     a program to calculate the sum of 
c     Bessel Functions.
c     The equation is:
c     V(x) = 2 *sum(n=-inf..inf)U(sqrt(x^2+(nw)^2))
c     Where U(r) is defined by:
c     U(r) = const. * Ko(r)
c     Ko is the 0th order bessel function
c-------------------------------------------------

      subroutine calcv(v,x,w,e)

c-------------------------------------------------
c     function declaration
      external u
c-------------------------------------------------

c-------------------------------------------------
c     Variable declarations
c-------------------------------------------------

c     this is the x in the equation
      real*8 x
c     v is the value that is returned
      real*8 v
c     w is the width of the slab
      real*8 w
c     e is the error tolerance
      real*8 e

c     index variable
      integer n
c-------------------------------------------------

c      write(0,*) "x: ",x
c      write(0,*) "w: ",w
c      write(0,*) "e: ",e

c-------------------------------------------------
c     we initialize v to be the 0th element in
c     the sum, then, because the function is
c     symmetric, we double the sum from
c     1 to infinity
c-------------------------------------------------
      dv    = 2.0d0 * u(x,0,w)
      v     = dv
      n     = 1
      
      do while (abs(dv) > e)
         dv = 4.0d0 * u(x,n,w)
         n  = n+1
         v  = v + dv
      end do

      return

      end


c-------------------------------------------------
c     this is the function u, as defined above.
c-------------------------------------------------
      real*8 function  u(x,n,w)
        implicit none
c-------------------------------------------------
c     variables
c-------------------------------------------------
c     reals for the arbitrary constants and
c     the bessel function.
c-------------------------------------------------
        real*8 x, y, w
        
        real*8 k, lambda, xi
        parameter (k = 1.0d0)
        parameter (lambda = 1.0d0)
        parameter (xi = lambda/2.0d1)
        
        real*8 dbesk0
c     same n as in the main routine
        integer n

c        write(0,*) 'w = ', w

        if (x .eq. 0.0d0 .and. n .eq. 0.0d0) then
c           write(0,*) "yo"
           u = log(lambda/xi)
c           write(0,*) "yo2"
           return
c           write(0,*) "yo3"
        end if
        
        y = sqrt(x**2 + (w**2)*(n**2))
c        write(0,*) n, y
        u = k * dbesk0(y/lambda)
        
c     for y values greater than this, dbesk0 may underflow
c     but this approximation is valid since the function
c     is approaching 0.
        if (y < 32d0) then
           u = u - k * dbesk0(y/xi)
        end if

        return
      end