25 double 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)));
52 double 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));
73 double 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));
91 double *alpha_act = alpha;
92 for (j = -1; j <= N; j++)
99 double *beta_act = beta;
100 for (j = -1; j <= N; j++)
107 double *gamma_act = gamma;
108 for (j = -1; j <= N; j++)
117 double *alpha_act = alpha;
118 for (i = -N; i <= N; i++)
120 for (j = -1; j <= N; j++)
131 double *alpha_act = alpha;
132 for (i = -N; i <= N; i++)
134 for (j = -1; j <= N; j++)
145 double *alpha_act = alpha;
146 for (i = -N; i <= N; i++)
148 for (j = -1; j <= N; j++)
162 double *alpha_act = alpha;
164 for (m = -N; m <= N; m++)
166 for (i = -N; i <= N; i++)
168 for (j = -1; j <= N; j++)
171 fprintf(stdout,
"alpha_all_%d^[%d,%d]=%f\n", j, i, m,
184 double *alpha_act = alpha;
185 for (m = -N; m <= N; m++)
187 for (i = -N; i <= N; i++)
189 for (j = -1; j <= N; j++)
201 double *alpha_act = alpha;
202 for (m = -N; m <= N; m++)
204 for (i = -N; i <= N; i++)
206 for (j = -1; j <= N; j++)
215 inline 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);
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)
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),
double wigner_start(int n1, int n2, double theta)
A method used for debugging, gives the values to start the "old" three-term recurrence generates WHE...
void SO3_gamma_row(double *gamma, int N, int m, int n)
Compute three-term-recurrence coefficients of Wigner-d functions for all degrees ...
void SO3_gamma_matrix(double *gamma, int N, int n)
Compute three-term-recurrence coefficients of Wigner-d functions for all order and degrees ...
void SO3_alpha_all(double *alpha, int N)
Compute three-term-recurrence coefficients of Wigner-d functions for all and .
void SO3_beta_row(double *beta, int N, int m, int n)
Compute three-term-recurrence coefficients of Wigner-d functions for all degrees ...
void SO3_alpha_matrix(double *alpha, int N, int n)
Compute three-term-recurrence coefficients of Wigner-d functions for all order and degrees ...
void SO3_gamma_all(double *gamma, int N)
Compute three-term-recurrence coefficients of Wigner-d functions for all and .
double SO3_alpha(int k, int m, int l)
Computes three-term recurrence coefficients of Wigner-d functions.
double SO3_beta(int k, int m, int l)
Computes three-term recurrence coefficients of Wigner-d functions.
double SO3_gamma(int k, int m, int l)
Computes three-term recurrence coefficients of Wigner-d functions.
void SO3_alpha_row(double *alpha, int N, int m, int n)
Compute three-term-recurrence coefficients of Wigner-d functions for all degrees ...
void SO3_beta_matrix(double *beta, int N, int n)
Compute three-term-recurrence coefficients of Wigner-d functions for all order and degrees ...
void SO3_beta_all(double *beta, int N)
Compute three-term-recurrence coefficients of Wigner-d functions for all and .
int eval_wigner_thresh(double *x, double *y, int size, int l, double *alpha, double *beta, double *gamma, double threshold)
Evaluates Wigner-d functions using the Clenshaw-algorithm if it not exceeds a given threshold...
Header file for functions related to Wigner-d/D functions.
void eval_wigner(double *x, double *y, int size, int l, double *alpha, double *beta, double *gamma)
Evaluates Wigner-d functions using the Clenshaw-algorithm.