LegendrePoly.f90 2.75 KB
Newer Older
Matthias Redies's avatar
Matthias Redies committed
1 2 3 4
module m_LegendrePoly
   implicit none

contains
5
   pure recursive function LegendrePoly(l,x) result(p)
Matthias Redies's avatar
Matthias Redies committed
6 7 8 9 10
      implicit none
      
      integer, intent(in) :: l
      real, intent(in)    :: x(:)
      real                :: p(size(x))
11
      real, parameter     :: one_128 = 1.0/128.0, one_256 = 1.0/256.0
Matthias Redies's avatar
Matthias Redies committed
12 13 14

      select case(l)
         case(0)
15
            p =             1.0
Matthias Redies's avatar
Matthias Redies committed
16
         case(1)
17
            p =             x
Matthias Redies's avatar
Matthias Redies committed
18
         case(2)
19
            p = 0.5     * ( 3   * x**2  - 1)
Matthias Redies's avatar
Matthias Redies committed
20
         case(3)
21
            p = 0.5     * ( 5   * x**3  -  3  * x)
Matthias Redies's avatar
Matthias Redies committed
22
         case(4)
23
            p = 0.125   * (35   * x**4  - 30  * x**2   + 3)
Matthias Redies's avatar
Matthias Redies committed
24
         case(5)
25
            p = 0.125   * (63   * x**5  - 70  * x**3   + 15  * x)
Matthias Redies's avatar
Matthias Redies committed
26
         case(6) 
Matthias Redies's avatar
Matthias Redies committed
27
            p = 0.0625  * (231  * x**6  -315  * x**4   + 105 * x**2   - 5)
Matthias Redies's avatar
Matthias Redies committed
28
         case(7)
29
            p = 0.0625  * (429  * x**7  -693  * x**5   + 315 * x**3   - 35    * x)
Matthias Redies's avatar
Matthias Redies committed
30
         case(8) 
31
            p = one_128 * (6435 * x**8  -12012 *x**6   + 6930* x**4   - 1260  * x**2   + 35)
Matthias Redies's avatar
Matthias Redies committed
32
         case(9)
33
            p = one_128 * (12155* x**9  -25740 *x**7   +18018* x**5   - 4620  * x**3   + 315 * x)
Matthias Redies's avatar
Matthias Redies committed
34
         case(10)
35
            p = one_256 * (46189* x**10 -109395*x**8   +90090* x**6   -30030  * x**4   +3465 * x**2 - 63)
36 37
         case default
            p = ( (2*l-1)*x*LegendrePoly(l-1,x)  - (l-1)*LegendrePoly(l-2,x) ) / l
Matthias Redies's avatar
Matthias Redies committed
38 39 40
      end select
   end function LegendrePoly
   
41
   pure recursive function LegendrePoly_scalar(l,x) result(p)
Matthias Redies's avatar
Matthias Redies committed
42 43 44 45 46
      implicit none
      
      integer, intent(in) :: l
      real, intent(in)    :: x
      real                :: p
47
      real, parameter     :: one_128 = 1.0/128.0, one_256 = 1.0/256.0
Matthias Redies's avatar
Matthias Redies committed
48 49 50

      select case(l)
         case(0)
51
            p =             1.0
Matthias Redies's avatar
Matthias Redies committed
52
         case(1)
53
            p =             x
Matthias Redies's avatar
Matthias Redies committed
54
         case(2)
55
            p = 0.5     * ( 3   * x**2  - 1)
Matthias Redies's avatar
Matthias Redies committed
56
         case(3)
57
            p = 0.5     * ( 5   * x**3  -  3  * x)
Matthias Redies's avatar
Matthias Redies committed
58
         case(4)
59
            p = 0.125   * (35   * x**4  - 30  * x**2   + 3)
Matthias Redies's avatar
Matthias Redies committed
60
         case(5)
61
            p = 0.125   * (63   * x**5  - 70  * x**3   + 15  * x)
Matthias Redies's avatar
Matthias Redies committed
62
         case(6) 
Matthias Redies's avatar
Matthias Redies committed
63
            p = 0.0625  * (231  * x**6  -315  * x**4   + 105 * x**2   - 5)
Matthias Redies's avatar
Matthias Redies committed
64
         case(7)
65
            p = 0.0625  * (429  * x**7  -693  * x**5   + 315 * x**3   - 35    * x)
Matthias Redies's avatar
Matthias Redies committed
66
         case(8) 
67
            p = one_128 * (6435 * x**8  -12012 *x**6   + 6930* x**4   - 1260  * x**2   + 35)
Matthias Redies's avatar
Matthias Redies committed
68
         case(9)
69
            p = one_128 * (12155* x**9  -25740 *x**7   +18018* x**5   - 4620  * x**3   + 315 * x)
Matthias Redies's avatar
Matthias Redies committed
70
         case(10)
71
            p = one_256 * (46189* x**10 -109395*x**8   +90090* x**6   -30030  * x**4   +3465 * x**2 - 63)
72 73
         case default
            p = ( (2*l-1)*x*LegendrePoly_scalar(l-1,x)  - (l-1)*LegendrePoly_scalar(l-2,x) ) / l
Matthias Redies's avatar
Matthias Redies committed
74 75 76 77
      end select
   end function LegendrePoly_scalar

end module m_LegendrePoly