Line data    Source code

       1             : module m_legendre_poly
       2             :    use m_types
       3             :    use m_constants
       4             : contains
       5           0 :    subroutine legendre_poly(x, P)
       6             :       implicit none
       7             :       real, intent(in)                 :: x(:)
       8             :       type(t_mat), intent(inout)       :: P
       9             : 
      10             :       integer :: n_max, n
      11             : 
      12           0 :       n_max =  P%matsize2 - 1
      13           0 :       if(P%matsize1 /= size(x))  then 
      14           0 :          write (*,*) "P:", P%matsize1, P%matsize2 
      15           0 :          write (*,*) "x:", size(x) 
      16           0 :          write (*,*) "n_max+1:", n_max+1 
      17             : 
      18           0 :          call judft_error("Leg dimensions don't agree")
      19             :       endif 
      20             :  
      21           0 :       if (n_max >= 0) then
      22           0 :          if(P%l_real) then
      23           0 :             P%data_r(:, 1) = 1.0
      24             :          else 
      25           0 :             P%data_c(:, 1) = cmplx_1
      26             :          endif
      27             :       endif
      28             : 
      29           0 :       if (n_max >= 1) then
      30           0 :          if(P%l_real) then
      31           0 :             P%data_r(:, 2) = x
      32             :          else 
      33           0 :             P%data_c(:, 2) = x
      34             :          endif
      35             :       endif
      36             : 
      37           0 :       do n = 1, n_max - 1
      38           0 :          if(P%l_real) then
      39           0 :             P%data_r(:, n + 2) = ((2*n + 1)*x*P%data_r(:, n+1) - n*P%data_r(:, n))/(n + 1.0)
      40             :          else
      41           0 :             P%data_c(:, n + 2) = ((2*n + 1)*x*P%data_c(:, n+1) - n*P%data_c(:, n))/(n + 1.0)
      42             :          endif
      43             :       enddo
      44           0 :    end subroutine legendre_poly
      45             : end module m_legendre_poly