lpmns Subroutine

public subroutine lpmns(m, n, x, Pm, Pd)

Compute associated Legendre functions Pmn(x) and Pmn'(x) for a given order

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: m

Order of Pmn(x), m = 0,1,2,...,n
integer, intent(in) :: n

Degree of Pmn(x), n = 0,1,2,...,N
real(kind=wp), intent(in) :: x

Argument of Pmn(x)
real(kind=wp), intent(out) :: Pm(0:n)

Pmn(x)
real(kind=wp), intent(out) :: Pd(0:n)

Pmn'(x)

Called by

proc~~lpmns~~CalledByGraph proc~lpmns specfun_module::lpmns proc~aswfb specfun_module::aswfb proc~aswfb->proc~lpmns proc~rmn2sp specfun_module::rmn2sp proc~rmn2sp->proc~lpmns proc~rswfp specfun_module::rswfp proc~rswfp->proc~rmn2sp

Source Code

   subroutine lpmns(m,n,x,Pm,Pd)

      real(wp),intent(in) :: x !! Argument of `Pmn(x)`
      integer,intent(in) :: m !! Order of `Pmn(x),  m = 0,1,2,...,n`
      integer,intent(in) :: n !! Degree of `Pmn(x), n = 0,1,2,...,N`
      real(wp),intent(out) :: Pm(0:n) !! `Pmn(x)`
      real(wp),intent(out) :: Pd(0:n) !! `Pmn'(x)`

      integer :: k
      real(wp) :: pm0 , pm1 , pm2 , pmk , x0

      do k = 0 , n
         Pm(k) = 0.0_wp
         Pd(k) = 0.0_wp
      enddo
      if ( abs(x)==1.0_wp ) then
         do k = 0 , n
            if ( m==0 ) then
               Pm(k) = 1.0_wp
               Pd(k) = 0.5_wp*k*(k+1.0_wp)
               if ( x<0.0_wp ) then
                  Pm(k) = (-1)**k*Pm(k)
                  Pd(k) = (-1)**(k+1)*Pd(k)
               endif
            elseif ( m==1 ) then
               Pd(k) = 1.0e+300_wp
            elseif ( m==2 ) then
               Pd(k) = -0.25_wp*(k+2.0_wp)*(k+1.0_wp)*k*(k-1.0_wp)
               if ( x<0.0_wp ) Pd(k) = (-1)**(k+1)*Pd(k)
            endif
         enddo
         return
      endif
      x0 = abs(1.0_wp-x*x)
      pm0 = 1.0_wp
      pmk = pm0
      do k = 1 , m
         pmk = (2.0_wp*k-1.0_wp)*sqrt(x0)*pm0
         pm0 = pmk
      enddo
      pm1 = (2.0_wp*m+1.0_wp)*x*pm0
      Pm(m) = pmk
      Pm(m+1) = pm1
      do k = m + 2 , n
         pm2 = ((2.0_wp*k-1.0_wp)*x*pm1-(k+m-1.0_wp)*pmk)/(k-m)
         Pm(k) = pm2
         pmk = pm1
         pm1 = pm2
      enddo
      Pd(0) = ((1.0_wp-m)*Pm(1)-x*Pm(0))/(x*x-1.0_wp)
      do k = 1 , n
         Pd(k) = (k*x*Pm(k)-(k+m)*Pm(k-1))/(x*x-1.0_wp)
      enddo
      do k = 1 , n
         Pm(k) = (-1)**m*Pm(k)
         Pd(k) = (-1)**m*Pd(k)
      enddo

   end subroutine lpmns