othpl Subroutine

subroutine othpl ( kf, n, x, pl, dpl )

  !*****************************************************************************80
  !
  !! OTHPL computes orthogonal polynomials Tn(x), Un(x), Ln(x) or Hn(x).
  !
  !  Discussion:
  !
  !    This procedure computes orthogonal polynomials: Tn(x) or Un(x),
  !    or Ln(x) or Hn(x), and their 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:
  !
  !    08 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, integer ( kind = 4 ) KT, the function code:
  !    1 for Chebyshev polynomial Tn(x)
  !    2 for Chebyshev polynomial Un(x)
  !    3 for Laguerre polynomial Ln(x)
  !    4 for Hermite polynomial Hn(x)
  !
  !    Input, integer ( kind = 4 ) N, the order.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) PL(0:N), DPL(0:N), the value and derivative of
  !    the polynomials of order 0 through N at X.
  !
  implicit none

  integer n

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  real ( kind = 8 ) c
  real ( kind = 8 ) dpl(0:n)
  real ( kind = 8 ) dy0
  real ( kind = 8 ) dy1
  real ( kind = 8 ) dyn
  integer ( kind = 4 ) k
  integer ( kind = 4 ) kf
  real ( kind = 8 ) pl(0:n)
  real ( kind = 8 ) x
  real ( kind = 8 ) y0
  real ( kind = 8 ) y1
  real ( kind = 8 ) yn

  a = 2.0D+00
  b = 0.0D+00
  c = 1.0D+00
  y0 = 1.0D+00
  y1 = 2.0D+00 * x
  dy0 = 0.0D+00
  dy1 = 2.0D+00
  pl(0) = 1.0D+00
  pl(1) = 2.0D+00 * x
  dpl(0) = 0.0D+00
  dpl(1) = 2.0D+00

  if ( kf == 1 ) then
     y1 = x
     dy1 = 1.0D+00
     pl(1) = x
     dpl(1) = 1.0D+00
  else if ( kf == 3 ) then
     y1 = 1.0D+00 - x
     dy1 = -1.0D+00
     pl(1) = 1.0D+00 - x
     dpl(1) = -1.0D+00
  end if

  do k = 2, n

     if ( kf == 3 ) then
        a = -1.0D+00 / k
        b = 2.0D+00 + a
        c = 1.0D+00 + a
     else if ( kf == 4 ) then
        c = 2.0D+00 * ( k - 1.0D+00 )
     end if

     yn = ( a * x + b ) * y1 - c * y0
     dyn = a * y1 + ( a * x + b ) * dy1 - c * dy0
     pl(k) = yn
     dpl(k) = dyn
     y0 = y1
     y1 = yn
     dy0 = dy1
     dy1 = dyn

  end do

  return
end subroutine othpl