NFFT 3.5.3alpha
legendre.c
1/*
2 * Copyright (c) 2002, 2017 Jens Keiner, Stefan Kunis, Daniel Potts
3 *
4 * This program is free software; you can redistribute it and/or modify it under
5 * the terms of the GNU General Public License as published by the Free Software
6 * Foundation; either version 2 of the License, or (at your option) any later
7 * version.
8 *
9 * This program is distributed in the hope that it will be useful, but WITHOUT
10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12 * details.
13 *
14 * You should have received a copy of the GNU General Public License along with
15 * this program; if not, write to the Free Software Foundation, Inc., 51
16 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17 */
18
19#include <math.h>
20#include <stdio.h>
21#include "infft.h"
22#include "legendre.h"
23#include "infft.h"
24
25/* One over sqrt(pi) */
26DK(KSQRTPII,0.56418958354775628694807945156077258584405062932900);
27
28static inline R alpha_al(const int k, const int n)
29{
30 if (k > 0)
31 {
32 if (k < n)
33 return IF(k%2,K(1.0),K(-1.0));
34 else
35 return SQRT(((R)(2*k+1))/((R)(k-n+1)))*SQRT((((R)(2*k+1))/((R)(k+n+1))));
36 }
37 else if (k == 0)
38 {
39 if (n == 0)
40 return K(1.0);
41 else
42 return IF(n%2,K(0.0),K(-1.0));
43 }
44 return K(0.0);
45}
46
47static inline R beta_al(const int k, const int n)
48{
49 if (0 <= k && k < n)
50 return K(1.0);
51 else
52 return K(0.0);
53}
54
55static inline R gamma_al(const int k, const int n)
56{
57 if (k == -1)
58 return SQRT(KSQRTPII*nfft_lambda((R)(n),K(0.5)));
59 else if (k <= n)
60 return K(0.0);
61 else
62 return -SQRT(((R)(k-n))/((R)(k-n+1))*((R)(k+n))/((R)(k+n+1)));
63}
64
65void alpha_al_row(R *alpha, const int N, const int n)
66{
67 int j;
68 R *p = alpha;
69 for (j = -1; j <= N; j++)
70 *p++ = alpha_al(j,n);
71}
72
73void beta_al_row(R *beta, const int N, const int n)
74{
75 int j;
76 R *p = beta;
77 for (j = -1; j <= N; j++)
78 *p++ = beta_al(j,n);
79}
80
81void gamma_al_row(R *gamma, const int N, const int n)
82{
83 int j;
84 R *p = gamma;
85 for (j = -1; j <= N; j++)
86 *p++ = gamma_al(j,n);
87}
88
89inline void alpha_al_all(R *alpha, const int N)
90{
91 int i,j;
92 R *p = alpha;
93 for (i = 0; i <= N; i++)
94 for (j = -1; j <= N; j++)
95 *p++ = alpha_al(j,i);
96}
97
98inline void beta_al_all(R *alpha, const int N)
99{
100 int i,j;
101 R *p = alpha;
102 for (i = 0; i <= N; i++)
103 for (j = -1; j <= N; j++)
104 *p++ = beta_al(j,i);
105}
106
107inline void gamma_al_all(R *alpha, const int N)
108{
109 int i,j;
110 R *p = alpha;
111 for (i = 0; i <= N; i++)
112 for (j = -1; j <= N; j++)
113 *p++ = gamma_al(j,i);
114}
115
116void eval_al(R *x, R *y, const int size, const int k, R *alpha,
117 R *beta, R *gamma)
118{
119 /* Evaluate the associated Legendre polynomial P_{k,nleg} (l,x) for the vector
120 * of knots x[0], ..., x[size-1] by the Clenshaw algorithm
121 */
122 int i,j;
123 R a,b,x_val_act,a_old;
124 R *x_act, *y_act;
125 R *alpha_act, *beta_act, *gamma_act;
126
127 /* Traverse all nodes. */
128 x_act = x;
129 y_act = y;
130 for (i = 0; i < size; i++)
131 {
132 a = 1.0;
133 b = 0.0;
134 x_val_act = *x_act;
135
136 if (k == 0)
137 {
138 *y_act = 1.0;
139 }
140 else
141 {
142 alpha_act = &(alpha[k]);
143 beta_act = &(beta[k]);
144 gamma_act = &(gamma[k]);
145 for (j = k; j > 1; j--)
146 {
147 a_old = a;
148 a = b + a_old*((*alpha_act)*x_val_act+(*beta_act));
149 b = a_old*(*gamma_act);
150 alpha_act--;
151 beta_act--;
152 gamma_act--;
153 }
154 *y_act = (a*((*alpha_act)*x_val_act+(*beta_act))+b);
155 }
156 x_act++;
157 y_act++;
158 }
159}
160
161int eval_al_thresh(R *x, R *y, const int size, const int k, R *alpha,
162 R *beta, R *gamma, R threshold)
163{
164 /* Evaluate the associated Legendre polynomial P_{k,nleg} (l,x) for the vector
165 * of knots x[0], ..., x[size-1] by the Clenshaw algorithm
166 */
167 int i,j;
168 R a,b,x_val_act,a_old;
169 R *x_act, *y_act;
170 R *alpha_act, *beta_act, *gamma_act;
171
172 /* Traverse all nodes. */
173 x_act = x;
174 y_act = y;
175 for (i = 0; i < size; i++)
176 {
177 a = 1.0;
178 b = 0.0;
179 x_val_act = *x_act;
180
181 if (k == 0)
182 {
183 *y_act = 1.0;
184 }
185 else
186 {
187 alpha_act = &(alpha[k]);
188 beta_act = &(beta[k]);
189 gamma_act = &(gamma[k]);
190 for (j = k; j > 1; j--)
191 {
192 a_old = a;
193 a = b + a_old*((*alpha_act)*x_val_act+(*beta_act));
194 b = a_old*(*gamma_act);
195 alpha_act--;
196 beta_act--;
197 gamma_act--;
198 }
199 *y_act = (a*((*alpha_act)*x_val_act+(*beta_act))+b);
200 if (fabs(*y_act) > threshold)
201 {
202 return 1;
203 }
204 }
205 x_act++;
206 y_act++;
207 }
208 return 0;
209}
int eval_al_thresh(R *x, R *y, const int size, const int k, R *alpha, R *beta, R *gamma, R threshold)
Evaluates an associated Legendre polynomials using the Clenshaw-algorithm if it no exceeds a given t...
Definition legendre.c:161
void alpha_al_all(R *alpha, const int N)
Compute three-term-recurrence coefficients of associated Legendre functions for .
Definition legendre.c:89
void gamma_al_all(R *alpha, const int N)
Compute three-term-recurrence coefficients of associated Legendre functions for .
Definition legendre.c:107
void eval_al(R *x, R *y, const int size, const int k, R *alpha, R *beta, R *gamma)
Evaluates an associated Legendre polynomials using the Clenshaw-algorithm.
Definition legendre.c:116
void beta_al_all(R *alpha, const int N)
Compute three-term-recurrence coefficients of associated Legendre functions for .
Definition legendre.c:98
Internal header file for auxiliary definitions and functions.