***************************************
* ベッセル関数を計算するサブ関数      *
* 参考文献：                          *
* http://www.math.meiji.ac.jp/~mk/labo/2004/computing-bessel-function/node34.html *
***************************************


      real function bessj(m,t)
      parameter( mmax=100 )    ! 数値積分用の配列
      integer n, m
      real x(mmax)
      real t
      real xmin, xmax, dx, pi
      pi=3.14159
      xmin=0.0
      xmax=2.0*pi
      dx=(xmax-xmin)/(mmax-1)
      do 10 i=1,mmax
         x(i)=xmin+dx*(i-1)
  10  continue

*-- 負の次数であった場合を分ける ---
      if(m.lt.0)then
         n=-m
      else
         n=m
      end if

*-- ベッセル関数の積分計算 ---
      bessj=0.5*dx*cos(t*sin(x(1))-real(n)*x(1))
      do 20 i=2,mmax-1
         bessj=bessj+0.5*dx*(cos(t*sin(x(i))-real(n)*x(i))
     &          +cos(t*sin(x(i+1))-real(n)*x(i+1)))
  20  continue
      bessj=bessj+0.5*dx*cos(t*sin(x(mmax))-real(n)*x(mmax))
      bessj=bessj/(2.0*pi)

*-- 負の次数であった場合を分ける ---
      if(m.lt.0)then
         bessj=((-1.0)**n)*bessj
      end if

      return
      end