pbwa Subroutine

subroutine pbwa ( a, x, w1f, w1d, w2f, w2d )

  !*****************************************************************************80
  !
  !! PBWA computes parabolic cylinder functions W(a,x) and 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:
  !
  !    29 July 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, real ( kind = 8 ) A, the parameter.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) W1F, W1D, W2F, W2D, the values of
  !    W(a,x), W'(a,x), W(a,-x), W'(a,-x).
  !
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) d(100)
  real ( kind = 8 ) d1
  real ( kind = 8 ) d2
  real ( kind = 8 ) dl
  real ( kind = 8 ) eps
  real ( kind = 8 ) f1
  real ( kind = 8 ) f2
  real ( kind = 8 ) g1
  real ( kind = 8 ) g2
  real ( kind = 8 ) h(100)
  real ( kind = 8 ) h0
  real ( kind = 8 ) h1
  real ( kind = 8 ) hl
  integer ( kind = 4 ) k
  integer ( kind = 4 ) l1
  integer ( kind = 4 ) l2
  integer ( kind = 4 ) m
  real ( kind = 8 ) p0
  real ( kind = 8 ) r
  real ( kind = 8 ) r1
  real ( kind = 8 ) ugi
  real ( kind = 8 ) ugr
  real ( kind = 8 ) vgi
  real ( kind = 8 ) vgr
  real ( kind = 8 ) w1d
  real ( kind = 8 ) w1f
  real ( kind = 8 ) w2d
  real ( kind = 8 ) w2f
  real ( kind = 8 ) x
  real ( kind = 8 ) x1
  real ( kind = 8 ) x2
  real ( kind = 8 ) y1
  real ( kind = 8 ) y1d
  real ( kind = 8 ) y1f
  real ( kind = 8 ) y2d
  real ( kind = 8 ) y2f

  eps = 1.0D-15
  p0 = 0.59460355750136D+00

  if ( a == 0.0D+00 ) then
     g1 = 3.625609908222D+00
     g2 = 1.225416702465D+00
  else
     x1 = 0.25D+00
     y1 = 0.5D+00 * a
     call cgama ( x1, y1, 1, ugr, ugi )
     g1 = sqrt ( ugr * ugr + ugi * ugi )
     x2 = 0.75D+00
     call cgama ( x2, y1, 1, vgr, vgi )
     g2 = sqrt ( vgr * vgr + vgi * vgi )
  end if

  f1 = sqrt ( g1 / g2 )
  f2 = sqrt ( 2.0D+00 * g2 / g1 )
  h0 = 1.0D+00
  h1 = a
  h(1) = a
  do l1 = 4, 200, 2
     m = l1 / 2
     hl = a * h1 - 0.25D+00 * ( l1 - 2.0D+00 ) * ( l1 - 3.0D+00 ) * h0
     h(m) = hl
     h0 = h1
     h1 = hl
  end do
  y1f = 1.0D+00
  r = 1.0D+00
  do k = 1, 100
     r = 0.5D+00 * r * x * x / ( k * ( 2.0D+00 * k - 1.0D+00 ) )
     r1 = h(k) * r
     y1f = y1f + r1
     if ( abs ( r1 / y1f ) <= eps .and. 30 < k ) then
        exit
     end if
  end do

  y1d = a
  r = 1.0D+00
  do k = 1, 100
     r = 0.5D+00 * r * x * x / ( k * ( 2.0D+00 * k + 1.0D+00 ) )
     r1 = h(k+1) * r
     y1d = y1d + r1
     if ( abs ( r1 / y1d ) <= eps .and. 30 < k ) then
        exit
     end if
  end do

  y1d = x * y1d
  d1 = 1.0D+00
  d2 = a
  d(1) = 1.0D+00
  d(2) = a
  do l2 = 5, 160, 2
     m = ( l2 + 1 ) / 2
     dl = a * d2 - 0.25D+00 * ( l2 - 2.0D+00 ) * ( l2 - 3.0D+00 ) * d1
     d(m) = dl
     d1 = d2
     d2 = dl
  end do

  y2f = 1.0D+00
  r = 1.0D+00
  do k = 1, 100
     r = 0.5D+00 * r * x * x / ( k * ( 2.0D+00 * k + 1.0D+00 ) )
     r1 = d(k+1) * r
     y2f = y2f + r1
     if ( abs ( r1 / y2f ) <= eps .and. 30 < k ) then
        exit
     end if
  end do

  y2f = x * y2f
  y2d = 1.0D+00
  r = 1.0D+00
  do k = 1, 100
     r = 0.5D+00 * r * x * x / ( k * ( 2.0D+00 * k - 1.0D+00 ) )
     r1 = d(k+1) * r
     y2d = y2d + r1
     if ( abs ( r1 / y2d ) <= eps .and. 30 < k ) then
        exit
     end if
  end do

  w1f = p0 * ( f1 * y1f - f2 * y2f )
  w2f = p0 * ( f1 * y1f + f2 * y2f )
  w1d = p0 * ( f1 * y1d - f2 * y2d )
  w2d = p0 * ( f1 * y1d + f2 * y2d )

  return
end subroutine pbwa