V3FIT
build_matrices_legendre.f
1  SUBROUTINE build_matrices_legendre(n, a, b, a_inv, b_inv)
2  USE stel_kinds
3  IMPLICIT NONE
4 !-----------------------------------------------
5 ! D u m m y A r g u m e n t s
6 !-----------------------------------------------
7  INTEGER, INTENT(IN) :: n
8  REAL(rprec), INTENT(OUT), DIMENSION(0:n,0:n) :: a, b,
9  1 a_inv, b_inv
10 !-----------------------------------------------
11  INTEGER :: idx, i, j
12  REAL(rprec) :: factorial
13  REAL(rprec), DIMENSION(13), PARAMETER :: norm_legend =
14  1 ( / 1, 1, 3, 5, 35, 63, 231, 429, 6435,
15  2 12155, 46189, 88179, 676039 / )
16  REAL(rprec), DIMENSION(13), PARAMETER :: norm_inv_legend =
17  1 ( / 1, 1, 2, 2, 8, 8, 16, 16, 128,
18  2 128, 256, 256, 1024 / )
19  REAL(rprec), DIMENSION(13*13), PARAMETER :: b_legend =
20  1 ( / 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
21  2 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
22  3 1, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
23  4 0, 3, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
24  5 7, 0, 20, 0, 8, 0, 0, 0, 0, 0, 0, 0, 0,
25  6 0, 27, 0, 28, 0, 8, 0, 0, 0, 0, 0, 0, 0,
26  7 33, 0, 110, 0, 72, 0, 16, 0, 0, 0, 0, 0, 0,
27  8 0, 143, 0, 182, 0, 88, 0, 16, 0, 0, 0, 0, 0,
28  9 715, 0, 2600, 0, 2160, 0, 832, 0, 128, 0, 0, 0, 0,
29  a 0, 3315, 0, 4760, 0, 2992, 0, 960, 0, 128, 0, 0, 0,
30  b 4199, 0, 16150, 0, 15504, 0, 7904, 0, 2176, 0,
31  c 256, 0, 0,
32  d 0, 20349, 0, 31654, 0, 23408, 0, 10080, 0, 2432,
33  e 0, 256, 0,
34  f 52003, 0, 208012, 0, 220248, 0, 133952, 0, 50048,
35  g 0, 10752, 0, 1024 / )
36  REAL(rprec), DIMENSION(13*13), PARAMETER :: b_legend_inv =
37  1 ( / 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
38  2 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
39  3 -1, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
40  4 0, -3, 0, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0,
41  5 3, 0, -30, 0, 35, 0, 0, 0, 0, 0, 0, 0, 0,
42  6 0, 15, 0, -70, 0, 63, 0, 0, 0, 0, 0, 0, 0,
43  7 -5, 0, 105, 0, -315, 0, 231, 0, 0, 0, 0, 0, 0,
44  8 0, -35, 0, 315, 0, -693, 0, 429, 0, 0, 0, 0, 0,
45  9 35, 0, -1260, 0, 6930, 0, -12012, 0, 6435, 0, 0,
46  a 0, 0,
47  b 0, 315, 0, -4620, 0, 18018, 0, -25740, 0, 12155,
48  c 0, 0, 0,
49  d -63, 0, 3465, 0, -30030, 0, 90090, 0, -109395, 0,
50  e 46189, 0, 0,
51  f 0, -693, 0, 15015, 0, -90090, 0, 218790, 0, -230945,
52  g 0, 88179, 0,
53  h 231, 0, -18018, 0, 225225, 0, -1021020, 0, 2078505,
54  i 0, -1939938, 0, 676039 / )
55 
56 !--------------------------------------------------------------------------------------------------
57 ! P_n(y):: Legendre polynomial in [-1,1]. Matrix coefficients from
58 ! "Handbook of Mathematical Functions", Abramowitz and Stegun, 1973, p.795.
59 ! If higher degree polyn. needed, then data set must be extended.
60 !
61 ! DESCRIPTION OF MATRICES:
62 !
63 ! A, A_inv: convert legendre Pol. to and from y-powers in [-1,1]
64 !
65 ! norm_legend(m)*y**m = sum_n=0^m (b_legend{mn} L_n)
66 ! norm_legend_inv(m)*L_m = sum_n=0^N (b_legend_inv{mn} y**n)
67 !
68 ! ==> A_{mn} = b_legend_{mn}/norm_legend(m) IF m<=n or 0. otherwise
69 ! ==> A_inv_{mn} = b_legend_inv_{mn}/norm_inv_legend(m) IF m<=n or 0. otherwise
70 !
71 ! B, B_inv: convert y-powers in [-1,1] to and from x-powers in [0,1]
72 !
73 ! x = (y + 1)/2; x**m = sum_n=0^m 2**(-m) m!/(n!(m-n)!) y**n
74 ! y = 2*x - 1; y**m = sum_n=0^m 2**n (-1)**(m-n) m!/(n!(m-n)!) x**n
75 !
76 ! ==> B_{mn} = 2**(-m) m!/(n!(m-n)!) IF m<=n or 0. otherwise
77 ! ==> B_inv_{mn} = 2**n (-1)**(m-n) m!/(n!(m-n)!) IF m<=n or 0. otherwise
78 !
79 !
80 ! NOTICE that: A_inv = (A)**(-1); B_inv = (B)**(-1)
81 !--------------------------------------------------------------------------------------------------
82 
83  IF (n > 12) stop 'N(legendre) CANNOT be larger than 12!!!'
84 
85  a = 0
86  a_inv = 0
87  b = 0
88  b_inv = 0
89 
90  DO i = 0, n
91  idx = 13*i
92  DO j = 0, i
93  a(i,j) = b_legend(idx+j+1)/norm_legend(i+1)
94  a_inv(i,j) = b_legend_inv(idx+j+1)/norm_inv_legend(i+1)
95  b(i,j) = 2._dp**(-i) * factorial(i)/
96  1 (factorial(j)*factorial(i-j))
97  b_inv(i,j) = 2._dp**j*(-1._dp)**(i-j)* factorial(i)/
98  1 (factorial(j)*factorial(i-j))
99  END DO
100  END DO
101 
102  END SUBROUTINE build_matrices_legendre