bess.f

Fortran的数学物理方程数值算法源程序。这是"Numerical Methods for Physics"第二版的Fortran源程序

      subroutine bess( m_max, x, jj )
      integer*4 m_max
      real*8 x, jj(*)
! Bessel function
! Inputs
!    m_max  Largest desired order
!    x = Value at which Bessel function J(x) is evaluated
! Output
!    jj = Vector of J(x) for order m = 0, 1, ..., m_max

      integer*4 MAXm_top
      parameter( MAXm_top = 1000 )
      integer*4 maxmx, m_top, m
      real*8 tinyNumber, j(MAXm_top+1), norm

      !* Perform downward recursion from initial guess
      if( m_max .gt. x ) then
        maxmx = m_max
      else             ! maxmx = Max(m_max,x)
        maxmx = int(x)
      endif
      ! Recursion is downward from m_top (which is even)
      m_top = 2*(int( (maxmx+15)/2 + 1 ))
      j(m_top+1) = 0.0
      j(m_top) = 1.0
      tinyNumber = 1e-16
      do m=(m_top-2), 0, -1       ! Downward recursion
        j(m+1) = 2*(m+1)/(x+tinyNumber)*j(m+2) - j(m+3)
      enddo

      !* Normalize using identity and return requested values
      norm = j(1)                ! NOTE: Be careful, m=0,1,... but
      do m=2, m_top, 2           ! vector goes j(1),j(2),...
        norm = norm + 2*j(m+1)
      enddo
      do m=0,m_max               ! Send back only the values for
        jj(m+1) = j(m+1)/norm    ! m=0,...,m_max and discard values
      enddo                      ! for m=m_max+1,...,m_top
      return
      end