27 double SO3_alpha(
const int m1,
const int m2,
const int j)
29 const int M = MAX(ABS(m1),ABS(m2)), mini = MIN(ABS(m1),ABS(m2));
35 if (m1 == 0 && m2 == 0)
40 return IF((m1+m2)%2,K(0.0),K(-0.5));
42 else if (j < M - mini)
43 return IF(j%2,K(0.5),K(-0.5));
45 return K(0.5) * SIGNF((R)m1)*SIGNF((R)m2);
48 SQRT(((R)(j+1))/((R)(j+1-m1)))
49 * SQRT(((R)(2*j+1))/((R)(j+1+m1)))
50 * SQRT(((R)(j+1))/((R)(j+1-m2)))
51 * SQRT(((R)(2*j+1))/((R)(j+1+m2)));
54 double SO3_beta(
const int m1,
const int m2,
const int j)
58 else if (j < MAX(ABS(m1),ABS(m2)))
60 else if (m1 == 0 || m2 == 0)
64 const R m1a = FABS((R)m1), m2a = FABS((R)m2);
66 ((SQRT(m1a)*SQRT(m2a))/((R)j))
67 * SQRT(m1a/((R)(j+1-m1)))
68 * SQRT(((R)(2*j+1))/((R)(j+1+m1)))
69 * SQRT(m2a/((R)(j+1-m2)))
70 * SQRT(((R)(2*j+1))/((R)(j+1+m2))),
71 SIGNF((R)m1)*SIGNF((R)m2));
75 double SO3_gamma(
const int m1,
const int m2,
const int j)
77 if (MAX(ABS(m1),ABS(m2)) < j)
78 return -(((R)(j+1))/((R)j)) * SQRT((((R)(j-m1))/((R)(j+1-m1)))
79 *(((R)(j+m1))/((R)(j+1+m1)))*(((R)(j-m2))/((R)(j+1-m2)))
80 *(((R)(j+m2))/((R)(j+1+m2))));
82 return IF(m1 > m2 || !((m1 + m2) % 2), K(1.0), K(-1.0))
83 * nfft_lambda2((R)ABS(m2 - m1),(R)ABS(m2 + m1));
93 double *alpha_act = alpha;
94 for (j = -1; j <= N; j++)
101 double *beta_act = beta;
102 for (j = -1; j <= N; j++)
109 double *gamma_act = gamma;
110 for (j = -1; j <= N; j++)
119 double *alpha_act = alpha;
120 for (i = -N; i <= N; i++)
122 for (j = -1; j <= N; j++)
133 double *alpha_act = alpha;
134 for (i = -N; i <= N; i++)
136 for (j = -1; j <= N; j++)
147 double *alpha_act = alpha;
148 for (i = -N; i <= N; i++)
150 for (j = -1; j <= N; j++)
164 double *alpha_act = alpha;
166 for (m = -N; m <= N; m++)
168 for (i = -N; i <= N; i++)
170 for (j = -1; j <= N; j++)
173 fprintf(stdout,
"alpha_all_%d^[%d,%d]=%f\n", j, i, m,
186 double *alpha_act = alpha;
187 for (m = -N; m <= N; m++)
189 for (i = -N; i <= N; i++)
191 for (j = -1; j <= N; j++)
203 double *alpha_act = alpha;
204 for (m = -N; m <= N; m++)
206 for (i = -N; i <= N; i++)
208 for (j = -1; j <= N; j++)
217 inline void eval_wigner(
double *x,
double *y,
int size,
int k,
double *alpha,
218 double *beta,
double *gamma)
224 double a, b, x_val_act, a_old;
225 double *x_act, *y_act;
226 double *alpha_act, *beta_act, *gamma_act;
231 for (i = 0; i < size; i++)
243 alpha_act = &(alpha[k]);
244 beta_act = &(beta[k]);
245 gamma_act = &(gamma[k]);
246 for (j = k; j > 1; j--)
249 a = b + a_old * ((*alpha_act) * x_val_act + (*beta_act));
250 b = a_old * (*gamma_act);
255 *y_act = (a * ((*alpha_act) * x_val_act + (*beta_act)) + b);
263 double *alpha,
double *beta,
double *gamma,
double threshold)
267 double a, b, x_val_act, a_old;
268 double *x_act, *y_act;
269 double *alpha_act, *beta_act, *gamma_act;
274 for (i = 0; i < size; i++)
286 alpha_act = &(alpha[k]);
287 beta_act = &(beta[k]);
288 gamma_act = &(gamma[k]);
289 for (j = k; j > 1; j--)
292 a = b + a_old * ((*alpha_act) * x_val_act + (*beta_act));
293 b = a_old * (*gamma_act);
298 *y_act = (a * ((*alpha_act) * x_val_act + (*beta_act)) + b);
299 if (fabs(*y_act) > threshold)
323 int cosPower, sinPower;
325 double dl, dm1, dm2, normFactor, sinSign;
330 max = (double) (ABS(m1) > ABS(m2) ? ABS(m1) : ABS(m2));
331 min = (double) (ABS(m1) < ABS(m2) ? ABS(m1) : ABS(m2));
344 for (i = 0; i < delta; i++)
345 normFactor *= SQRT((2. * dl - ((
double) i)) / (((double) i) + 1.));
349 normFactor *= SQRT((2. * dl + 1.) / 2.);
377 dCP = (double) cosPower;
378 dSP = (double) sinPower;
380 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.