jyndd Subroutine

subroutine jyndd ( n, x, bjn, djn, fjn, byn, dyn, fyn )

  !*****************************************************************************80
  !
  !! JYNDD: Bessel functions Jn(x) and Yn(x), first and second derivatives.
  !
  !  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:
  !
  !    02 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 ) N, the order.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) BJN, DJN, FJN, BYN, DYN, FYN, the values of
  !    Jn(x), Jn'(x), Jn"(x), Yn(x), Yn'(x), Yn"(x).
  !
  implicit none

  real ( kind = 8 ) bj(102)
  real ( kind = 8 ) bjn
  real ( kind = 8 ) byn
  real ( kind = 8 ) bs
  real ( kind = 8 ) by(102)
  real ( kind = 8 ) djn
  real ( kind = 8 ) dyn
  real ( kind = 8 ) e0
  real ( kind = 8 ) ec
  real ( kind = 8 ) f
  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  real ( kind = 8 ) fjn
  real ( kind = 8 ) fyn
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  integer ( kind = 4 ) mt
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nt
  real ( kind = 8 ) s1
  real ( kind = 8 ) su
  real ( kind = 8 ) x

  do nt = 1, 900
     mt = int ( 0.5D+00 * log10 ( 6.28D+00 * nt ) &
          - nt * log10 ( 1.36D+00 * abs ( x ) / nt ) )
     if ( 20 < mt ) then
        exit
     end if
  end do

  m = nt
  bs = 0.0D+00
  f0 = 0.0D+00
  f1 = 1.0D-35
  su = 0.0D+00
  do k = m, 0, -1
     f = 2.0D+00 * ( k + 1.0D+00 ) * f1 / x - f0
     if ( k <= n + 1 ) then
        bj(k+1) = f
     end if
     if ( k == 2 * int ( k / 2 ) ) then
        bs = bs + 2.0D+00 * f
        if ( k /= 0 ) then
           su = su + ( -1.0D+00 ) ** ( k / 2 ) * f / k
        end if
     end if
     f0 = f1
     f1 = f
  end do

  do k = 0, n + 1
     bj(k+1) = bj(k+1) / ( bs - f )
  end do

  bjn = bj(n+1)
  ec = 0.5772156649015329D+00
  e0 = 0.3183098861837907D+00
  s1 = 2.0D+00 * e0 * ( log ( x / 2.0D+00 ) + ec ) * bj(1)
  f0 = s1 - 8.0D+00 * e0 * su / ( bs - f )
  f1 = ( bj(2) * f0 - 2.0D+00 * e0 / x ) / bj(1)

  by(1) = f0
  by(2) = f1
  do k = 2, n + 1 
     f = 2.0D+00 * ( k - 1.0D+00 ) * f1 / x - f0
     by(k+1) = f
     f0 = f1
     f1 = f
  end do

  byn = by(n+1)
  djn = - bj(n+2) + n * bj(n+1) / x
  dyn = - by(n+2) + n * by(n+1) / x
  fjn = ( n * n / ( x * x ) - 1.0D+00 ) * bjn - djn / x
  fyn = ( n * n / ( x * x ) - 1.0D+00 ) * byn - dyn / x

  return
end subroutine jyndd