25double SO3_alpha(
const int m1,
const int m2,
const int j)
27 const int M = MAX(ABS(m1),ABS(m2)), mini = MIN(ABS(m1),ABS(m2));
33 if (m1 == 0 && m2 == 0)
38 return IF((m1+m2)%2,K(0.0),K(-0.5));
40 else if (j < M - mini)
41 return IF(j%2,K(0.5),K(-0.5));
43 return K(0.5) * SIGNF((R)m1)*SIGNF((R)m2);
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)));
52double SO3_beta(
const int m1,
const int m2,
const int j)
56 else if (j < MAX(ABS(m1),ABS(m2)))
58 else if (m1 == 0 || m2 == 0)
62 const R m1a = FABS((R)m1), m2a = FABS((R)m2);
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));
73double SO3_gamma(
const int m1,
const int m2,
const int j)
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))));
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));
215inline void eval_wigner(
double *x,
double *y,
int size,
int k,
double *alpha,
216 double *beta,
double *gamma)
222 double a, b, x_val_act, a_old;
223 double *x_act, *y_act;
224 double *alpha_act, *beta_act, *gamma_act;
229 for (i = 0; i < size; i++)
241 alpha_act = &(alpha[k]);
242 beta_act = &(beta[k]);
243 gamma_act = &(gamma[k]);
244 for (j = k; j > 1; j--)
247 a = b + a_old * ((*alpha_act) * x_val_act + (*beta_act));
248 b = a_old * (*gamma_act);
253 *y_act = (a * ((*alpha_act) * x_val_act + (*beta_act)) + b);
260inline int eval_wigner_thresh(
double *x,
double *y,
int size,
int k,
261 double *alpha,
double *beta,
double *gamma,
double threshold)
265 double a, b, x_val_act, a_old;
266 double *x_act, *y_act;
267 double *alpha_act, *beta_act, *gamma_act;
272 for (i = 0; i < size; i++)
284 alpha_act = &(alpha[k]);
285 beta_act = &(beta[k]);
286 gamma_act = &(gamma[k]);
287 for (j = k; j > 1; j--)
290 a = b + a_old * ((*alpha_act) * x_val_act + (*beta_act));
291 b = a_old * (*gamma_act);
296 *y_act = (a * ((*alpha_act) * x_val_act + (*beta_act)) + b);
297 if (fabs(*y_act) > threshold)
317double wigner_start(
int m1,
int m2,
double theta)
321 int cosPower, sinPower;
323 double dl, dm1, dm2, normFactor, sinSign;
328 max = (double) (ABS(m1) > ABS(m2) ? ABS(m1) : ABS(m2));
329 min = (double) (ABS(m1) < ABS(m2) ? ABS(m1) : ABS(m2));
342 for (i = 0; i < delta; i++)
343 normFactor *= SQRT((2. * dl - ((
double) i)) / (((double) i) + 1.));
347 normFactor *= SQRT((2. * dl + 1.) / 2.);
375 dCP = (double) cosPower;
376 dSP = (double) sinPower;
378 return normFactor * sinSign * POW(sin(theta / 2), dSP) * POW(cos(theta / 2),