NFFT 3.5.3alpha
wigner.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 "wigner.h"
23#include "infft.h"
24
25double SO3_alpha(const int m1, const int m2, const int j)
26{
27 const int M = MAX(ABS(m1),ABS(m2)), mini = MIN(ABS(m1),ABS(m2));
28
29 if (j < 0)
30 return K(0.0);
31 else if (j == 0)
32 {
33 if (m1 == 0 && m2 == 0)
34 return K(1.0);
35 if (m1 == m2)
36 return K(0.5);
37 else
38 return IF((m1+m2)%2,K(0.0),K(-0.5));
39 }
40 else if (j < M - mini)
41 return IF(j%2,K(0.5),K(-0.5));
42 else if (j < M)
43 return K(0.5) * SIGNF((R)m1)*SIGNF((R)m2);
44 else
45 return
46 SQRT(((R)(j+1))/((R)(j+1-m1)))
47 * SQRT(((R)(2*j+1))/((R)(j+1+m1)))
48 * SQRT(((R)(j+1))/((R)(j+1-m2)))
49 * SQRT(((R)(2*j+1))/((R)(j+1+m2)));
50}
51
52double SO3_beta(const int m1, const int m2, const int j)
53{
54 if (j < 0)
55 return K(0.0);
56 else if (j < MAX(ABS(m1),ABS(m2)))
57 return K(0.5);
58 else if (m1 == 0 || m2 == 0)
59 return K(0.0);
60 else
61 {
62 const R m1a = FABS((R)m1), m2a = FABS((R)m2);
63 return -COPYSIGN(
64 ((SQRT(m1a)*SQRT(m2a))/((R)j))
65 * SQRT(m1a/((R)(j+1-m1)))
66 * SQRT(((R)(2*j+1))/((R)(j+1+m1)))
67 * SQRT(m2a/((R)(j+1-m2)))
68 * SQRT(((R)(2*j+1))/((R)(j+1+m2))),
69 SIGNF((R)m1)*SIGNF((R)m2));
70 }
71}
72
73double SO3_gamma(const int m1, const int m2, const int j)
74{
75 if (MAX(ABS(m1),ABS(m2)) < j)
76 return -(((R)(j+1))/((R)j)) * SQRT((((R)(j-m1))/((R)(j+1-m1)))
77 *(((R)(j+m1))/((R)(j+1+m1)))*(((R)(j-m2))/((R)(j+1-m2)))
78 *(((R)(j+m2))/((R)(j+1+m2))));
79 else if (j == -1)
80 return IF(m1 > m2 || !((m1 + m2) % 2), K(1.0), K(-1.0))
81 * nfft_lambda2((R)ABS(m2 - m1),(R)ABS(m2 + m1));
82 else
83 return K(0.0);
84}
85
86/*compute the coefficients for all degrees*/
87
88inline void SO3_alpha_row(double *alpha, int N, int k, int m)
89{
90 int j;
91 double *alpha_act = alpha;
92 for (j = -1; j <= N; j++)
93 *alpha_act++ = SO3_alpha(k, m, j);
94}
95
96inline void SO3_beta_row(double *beta, int N, int k, int m)
97{
98 int j;
99 double *beta_act = beta;
100 for (j = -1; j <= N; j++)
101 *beta_act++ = SO3_beta(k, m, j);
102}
103
104inline void SO3_gamma_row(double *gamma, int N, int k, int m)
105{
106 int j;
107 double *gamma_act = gamma;
108 for (j = -1; j <= N; j++)
109 *gamma_act++ = SO3_gamma(k, m, j);
110}
111
112/*compute for all degrees l and orders k*/
113
114inline void SO3_alpha_matrix(double *alpha, int N, int m)
115{
116 int i, j;
117 double *alpha_act = alpha;
118 for (i = -N; i <= N; i++)
119 {
120 for (j = -1; j <= N; j++)
121 {
122 *alpha_act = SO3_alpha(i, m, j);
123 alpha_act++;
124 }
125 }
126}
127
128inline void SO3_beta_matrix(double *alpha, int N, int m)
129{
130 int i, j;
131 double *alpha_act = alpha;
132 for (i = -N; i <= N; i++)
133 {
134 for (j = -1; j <= N; j++)
135 {
136 *alpha_act = SO3_beta(i, m, j);
137 alpha_act++;
138 }
139 }
140}
141
142inline void SO3_gamma_matrix(double *alpha, int N, int m)
143{
144 int i, j;
145 double *alpha_act = alpha;
146 for (i = -N; i <= N; i++)
147 {
148 for (j = -1; j <= N; j++)
149 {
150 *alpha_act = SO3_gamma(i, m, j);
151 alpha_act++;
152 }
153 }
154}
155
156/*compute all 3termrecurrence coeffs*/
157
158inline void SO3_alpha_all(double *alpha, int N)
159{
160 int q;
161 int i, j, m;
162 double *alpha_act = alpha;
163 q = 0;
164 for (m = -N; m <= N; m++)
165 {
166 for (i = -N; i <= N; i++)
167 {
168 for (j = -1; j <= N; j++)
169 {
170 *alpha_act = SO3_alpha(i, m, j);
171 fprintf(stdout, "alpha_all_%d^[%d,%d]=%f\n", j, i, m,
172 SO3_alpha(i, m, j));
173 alpha_act++;
174 q = q + 1;
175
176 }
177 }
178 }
179}
180
181inline void SO3_beta_all(double *alpha, int N)
182{
183 int i, j, m;
184 double *alpha_act = alpha;
185 for (m = -N; m <= N; m++)
186 {
187 for (i = -N; i <= N; i++)
188 {
189 for (j = -1; j <= N; j++)
190 {
191 *alpha_act = SO3_beta(i, m, j);
192 alpha_act++;
193 }
194 }
195 }
196}
197
198inline void SO3_gamma_all(double *alpha, int N)
199{
200 int i, j, m;
201 double *alpha_act = alpha;
202 for (m = -N; m <= N; m++)
203 {
204 for (i = -N; i <= N; i++)
205 {
206 for (j = -1; j <= N; j++)
207 {
208 *alpha_act = SO3_gamma(i, m, j);
209 alpha_act++;
210 }
211 }
212 }
213}
214
215inline void eval_wigner(double *x, double *y, int size, int k, double *alpha,
216 double *beta, double *gamma)
217{
218 /* Evaluate the wigner function d_{k,nleg} (l,x) for the vector
219 * of knots x[0], ..., x[size-1] by the Clenshaw algorithm
220 */
221 int i, j;
222 double a, b, x_val_act, a_old;
223 double *x_act, *y_act;
224 double *alpha_act, *beta_act, *gamma_act;
225
226 /* Traverse all nodes. */
227 x_act = x;
228 y_act = y;
229 for (i = 0; i < size; i++)
230 {
231 a = 1.0;
232 b = 0.0;
233 x_val_act = *x_act;
234
235 if (k == 0)
236 {
237 *y_act = 1.0;
238 }
239 else
240 {
241 alpha_act = &(alpha[k]);
242 beta_act = &(beta[k]);
243 gamma_act = &(gamma[k]);
244 for (j = k; j > 1; j--)
245 {
246 a_old = a;
247 a = b + a_old * ((*alpha_act) * x_val_act + (*beta_act));
248 b = a_old * (*gamma_act);
249 alpha_act--;
250 beta_act--;
251 gamma_act--;
252 }
253 *y_act = (a * ((*alpha_act) * x_val_act + (*beta_act)) + b);
254 }
255 x_act++;
256 y_act++;
257 }
258}
259
260inline int eval_wigner_thresh(double *x, double *y, int size, int k,
261 double *alpha, double *beta, double *gamma, double threshold)
262{
263
264 int i, j;
265 double a, b, x_val_act, a_old;
266 double *x_act, *y_act;
267 double *alpha_act, *beta_act, *gamma_act;
268
269 /* Traverse all nodes. */
270 x_act = x;
271 y_act = y;
272 for (i = 0; i < size; i++)
273 {
274 a = 1.0;
275 b = 0.0;
276 x_val_act = *x_act;
277
278 if (k == 0)
279 {
280 *y_act = 1.0;
281 }
282 else
283 {
284 alpha_act = &(alpha[k]);
285 beta_act = &(beta[k]);
286 gamma_act = &(gamma[k]);
287 for (j = k; j > 1; j--)
288 {
289 a_old = a;
290 a = b + a_old * ((*alpha_act) * x_val_act + (*beta_act));
291 b = a_old * (*gamma_act);
292 alpha_act--;
293 beta_act--;
294 gamma_act--;
295 }
296 *y_act = (a * ((*alpha_act) * x_val_act + (*beta_act)) + b);
297 if (fabs(*y_act) > threshold)
298 {
299 return 1;
300 }
301 }
302 x_act++;
303 y_act++;
304 }
305 return 0;
306}
307
308/************************************************************************/
309/* L2 normed wigner little d, WHERE THE DEGREE OF THE FUNCTION IS EQUAL
310 TO ONE OF ITS ORDERS. This is the function to use when starting the
311 three-term recurrence at orders (m1,m2)
312
313 Note that, by definition, since I am starting the recurrence with this
314 function, that the degree j of the function is equal to max(abs(m1), abs(m2) ).
315 */
316
317double wigner_start(int m1, int m2, double theta)
318{
319
320 int i, l, delta;
321 int cosPower, sinPower;
322 int absM1, absM2;
323 double dl, dm1, dm2, normFactor, sinSign;
324 double dCP, dSP;
325 double max;
326 double min;
327
328 max = (double) (ABS(m1) > ABS(m2) ? ABS(m1) : ABS(m2));
329 min = (double) (ABS(m1) < ABS(m2) ? ABS(m1) : ABS(m2));
330
331 l = max;
332 delta = l - min;
333
334 absM1 = ABS(m1);
335 absM2 = ABS(m2);
336 dl = (double) l;
337 dm1 = (double) m1;
338 dm2 = (double) m2;
339 sinSign = 1.;
340 normFactor = 1.;
341
342 for (i = 0; i < delta; i++)
343 normFactor *= SQRT((2. * dl - ((double) i)) / (((double) i) + 1.));
344
345 /* need to adjust to make the L2-norm equal to 1 */
346
347 normFactor *= SQRT((2. * dl + 1.) / 2.);
348
349 if (l == absM1)
350 if (m1 >= 0)
351 {
352 cosPower = l + m2;
353 sinPower = l - m2;
354 if ((l - m2) % 2)
355 sinSign = -1.;
356 }
357 else
358 {
359 cosPower = l - m2;
360 sinPower = l + m2;
361 }
362 else if (m2 >= 0)
363 {
364 cosPower = l + m1;
365 sinPower = l - m1;
366 }
367 else
368 {
369 cosPower = l - m1;
370 sinPower = l + m1;
371 if ((l + m1) % 2)
372 sinSign = -1.;
373 }
374
375 dCP = (double) cosPower;
376 dSP = (double) sinPower;
377
378 return normFactor * sinSign * POW(sin(theta / 2), dSP) * POW(cos(theta / 2),
379 dCP);
380}
Internal header file for auxiliary definitions and functions.
Header file for functions related to Wigner-d/D functions.