cjk Subroutine

subroutine cjk ( km, a )

  !*****************************************************************************80
  !
  !! CJK: asymptotic expansion coefficients for Bessel functions of large order.
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    01 August 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  !
  !  Parameters:
  !
  !    Input, integer ( kind = 4 ) KM, the maximum value of K.
  !
  !    Output, real ( kind = 8 ) A(L), the value of Cj(k) where j and k are 
  !    related to L by L = j+1+[k*(k+1)]/2; j,k = 0,1,...,Km.
  !
  implicit none

  real ( kind = 8 ) a(*)
  real ( kind = 8 ) f
  real ( kind = 8 ) f0
  real ( kind = 8 ) g
  real ( kind = 8 ) g0
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) km
  integer ( kind = 4 ) l1
  integer ( kind = 4 ) l2
  integer ( kind = 4 ) l3
  integer ( kind = 4 ) l4

  a(1) = 1.0D+00
  f0 = 1.0D+00
  g0 = 1.0D+00
  do k = 0, km - 1
     l1 = ( k + 1 ) * ( k + 2 ) / 2 + 1
     l2 = ( k + 1 ) * ( k + 2 ) / 2 + k + 2
     f = ( 0.5D+00 * k + 0.125D+00 / ( k + 1 ) ) * f0
     g = - ( 1.5D+00 * k + 0.625D+00 &
          / ( 3.0D+00 * ( k + 1.0D+00 ) ) ) * g0
     a(l1) = f
     a(l2) = g
     f0 = f
     g0 = g
  end do

  do k = 1, km - 1
     do j = 1, k
        l3 = k * ( k + 1 ) / 2 + j + 1
        l4 = ( k + 1 ) * ( k + 2 ) / 2 + j + 1
        a(l4) = ( j + 0.5D+00 * k + 0.125D+00 &
             / ( 2.0D+00 * j + k + 1.0D+00 ) ) * a(l3) &
             - ( j + 0.5D+00 * k - 1.0D+00 + 0.625D+00 &
             / ( 2.0D+00 * j + k + 1.0D+00 ) ) * a(l3-1)
     end do
  end do

  return
end subroutine cjk