00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include <math.h>
00022 #include <stdio.h>
00023 #include "infft.h"
00024 #include "wigner.h"
00025 #include "nfft3util.h"
00026
00027 double SO3_alpha(const int m1, const int m2, const int j)
00028 {
00029 const int M = MAX(ABS(m1),ABS(m2)), mini = MIN(ABS(m1),ABS(m2));
00030
00031 if (j < 0)
00032 return K(0.0);
00033 else if (j == 0)
00034 {
00035 if (m1 == 0 && m2 == 0)
00036 return K(1.0);
00037 if (m1 == m2)
00038 return K(0.5);
00039 else
00040 return IF((m1+m2)%2,K(0.0),K(-0.5));
00041 }
00042 else if (j < M - mini)
00043 return IF(j%2,K(0.5),K(-0.5));
00044 else if (j < M)
00045 return K(0.5) * SIGNF((R)m1)*SIGNF((R)m2);
00046 else
00047 return
00048 SQRT(((R)(j+1))/((R)(j+1-m1)))
00049 * SQRT(((R)(2*j+1))/((R)(j+1+m1)))
00050 * SQRT(((R)(j+1))/((R)(j+1-m2)))
00051 * SQRT(((R)(2*j+1))/((R)(j+1+m2)));
00052 }
00053
00054 double SO3_beta(const int m1, const int m2, const int j)
00055 {
00056 if (j < 0)
00057 return K(0.0);
00058 else if (j < MAX(ABS(m1),ABS(m2)))
00059 return K(0.5);
00060 else if (m1 == 0 || m2 == 0)
00061 return K(0.0);
00062 else
00063 {
00064 const R m1a = FABS((R)m1), m2a = FABS((R)m2);
00065 return -COPYSIGN(
00066 ((SQRT(m1a)*SQRT(m2a))/((R)j))
00067 * SQRT(m1a/((R)(j+1-m1)))
00068 * SQRT(((R)(2*j+1))/((R)(j+1+m1)))
00069 * SQRT(m2a/((R)(j+1-m2)))
00070 * SQRT(((R)(2*j+1))/((R)(j+1+m2))),
00071 SIGNF((R)m1)*SIGNF((R)m2));
00072 }
00073 }
00074
00075 double SO3_gamma(const int m1, const int m2, const int j)
00076 {
00077 if (MAX(ABS(m1),ABS(m2)) < j)
00078 return -(((R)(j+1))/((R)j)) * SQRT((((R)(j-m1))/((R)(j+1-m1)))
00079 *(((R)(j+m1))/((R)(j+1+m1)))*(((R)(j-m2))/((R)(j+1-m2)))
00080 *(((R)(j+m2))/((R)(j+1+m2))));
00081 else if (j == -1)
00082 return IF(m1 > m2 || !((m1 + m2) % 2), K(1.0), K(-1.0))
00083 * nfft_lambda2((R)ABS(m2 - m1),(R)ABS(m2 + m1));
00084 else
00085 return K(0.0);
00086 }
00087
00088
00089
00090 inline void SO3_alpha_row(double *alpha, int N, int k, int m)
00091 {
00092 int j;
00093 double *alpha_act = alpha;
00094 for (j = -1; j <= N; j++)
00095 *alpha_act++ = SO3_alpha(k, m, j);
00096 }
00097
00098 inline void SO3_beta_row(double *beta, int N, int k, int m)
00099 {
00100 int j;
00101 double *beta_act = beta;
00102 for (j = -1; j <= N; j++)
00103 *beta_act++ = SO3_beta(k, m, j);
00104 }
00105
00106 inline void SO3_gamma_row(double *gamma, int N, int k, int m)
00107 {
00108 int j;
00109 double *gamma_act = gamma;
00110 for (j = -1; j <= N; j++)
00111 *gamma_act++ = SO3_gamma(k, m, j);
00112 }
00113
00114
00115
00116 inline void SO3_alpha_matrix(double *alpha, int N, int m)
00117 {
00118 int i, j;
00119 double *alpha_act = alpha;
00120 for (i = -N; i <= N; i++)
00121 {
00122 for (j = -1; j <= N; j++)
00123 {
00124 *alpha_act = SO3_alpha(i, m, j);
00125 alpha_act++;
00126 }
00127 }
00128 }
00129
00130 inline void SO3_beta_matrix(double *alpha, int N, int m)
00131 {
00132 int i, j;
00133 double *alpha_act = alpha;
00134 for (i = -N; i <= N; i++)
00135 {
00136 for (j = -1; j <= N; j++)
00137 {
00138 *alpha_act = SO3_beta(i, m, j);
00139 alpha_act++;
00140 }
00141 }
00142 }
00143
00144 inline void SO3_gamma_matrix(double *alpha, int N, int m)
00145 {
00146 int i, j;
00147 double *alpha_act = alpha;
00148 for (i = -N; i <= N; i++)
00149 {
00150 for (j = -1; j <= N; j++)
00151 {
00152 *alpha_act = SO3_gamma(i, m, j);
00153 alpha_act++;
00154 }
00155 }
00156 }
00157
00158
00159
00160 inline void SO3_alpha_all(double *alpha, int N)
00161 {
00162 int q;
00163 int i, j, m;
00164 double *alpha_act = alpha;
00165 q = 0;
00166 for (m = -N; m <= N; m++)
00167 {
00168 for (i = -N; i <= N; i++)
00169 {
00170 for (j = -1; j <= N; j++)
00171 {
00172 *alpha_act = SO3_alpha(i, m, j);
00173 fprintf(stdout, "alpha_all_%d^[%d,%d]=%f\n", j, i, m,
00174 SO3_alpha(i, m, j));
00175 alpha_act++;
00176 q = q + 1;
00177
00178 }
00179 }
00180 }
00181 }
00182
00183 inline void SO3_beta_all(double *alpha, int N)
00184 {
00185 int i, j, m;
00186 double *alpha_act = alpha;
00187 for (m = -N; m <= N; m++)
00188 {
00189 for (i = -N; i <= N; i++)
00190 {
00191 for (j = -1; j <= N; j++)
00192 {
00193 *alpha_act = SO3_beta(i, m, j);
00194 alpha_act++;
00195 }
00196 }
00197 }
00198 }
00199
00200 inline void SO3_gamma_all(double *alpha, int N)
00201 {
00202 int i, j, m;
00203 double *alpha_act = alpha;
00204 for (m = -N; m <= N; m++)
00205 {
00206 for (i = -N; i <= N; i++)
00207 {
00208 for (j = -1; j <= N; j++)
00209 {
00210 *alpha_act = SO3_gamma(i, m, j);
00211 alpha_act++;
00212 }
00213 }
00214 }
00215 }
00216
00217 inline void eval_wigner(double *x, double *y, int size, int k, double *alpha,
00218 double *beta, double *gamma)
00219 {
00220
00221
00222
00223 int i, j;
00224 double a, b, x_val_act, a_old;
00225 double *x_act, *y_act;
00226 double *alpha_act, *beta_act, *gamma_act;
00227
00228
00229 x_act = x;
00230 y_act = y;
00231 for (i = 0; i < size; i++)
00232 {
00233 a = 1.0;
00234 b = 0.0;
00235 x_val_act = *x_act;
00236
00237 if (k == 0)
00238 {
00239 *y_act = 1.0;
00240 }
00241 else
00242 {
00243 alpha_act = &(alpha[k]);
00244 beta_act = &(beta[k]);
00245 gamma_act = &(gamma[k]);
00246 for (j = k; j > 1; j--)
00247 {
00248 a_old = a;
00249 a = b + a_old * ((*alpha_act) * x_val_act + (*beta_act));
00250 b = a_old * (*gamma_act);
00251 alpha_act--;
00252 beta_act--;
00253 gamma_act--;
00254 }
00255 *y_act = (a * ((*alpha_act) * x_val_act + (*beta_act)) + b);
00256 }
00257 x_act++;
00258 y_act++;
00259 }
00260 }
00261
00262 inline int eval_wigner_thresh(double *x, double *y, int size, int k,
00263 double *alpha, double *beta, double *gamma, double threshold)
00264 {
00265
00266 int i, j;
00267 double a, b, x_val_act, a_old;
00268 double *x_act, *y_act;
00269 double *alpha_act, *beta_act, *gamma_act;
00270
00271
00272 x_act = x;
00273 y_act = y;
00274 for (i = 0; i < size; i++)
00275 {
00276 a = 1.0;
00277 b = 0.0;
00278 x_val_act = *x_act;
00279
00280 if (k == 0)
00281 {
00282 *y_act = 1.0;
00283 }
00284 else
00285 {
00286 alpha_act = &(alpha[k]);
00287 beta_act = &(beta[k]);
00288 gamma_act = &(gamma[k]);
00289 for (j = k; j > 1; j--)
00290 {
00291 a_old = a;
00292 a = b + a_old * ((*alpha_act) * x_val_act + (*beta_act));
00293 b = a_old * (*gamma_act);
00294 alpha_act--;
00295 beta_act--;
00296 gamma_act--;
00297 }
00298 *y_act = (a * ((*alpha_act) * x_val_act + (*beta_act)) + b);
00299 if (fabs(*y_act) > threshold)
00300 {
00301 return 1;
00302 }
00303 }
00304 x_act++;
00305 y_act++;
00306 }
00307 return 0;
00308 }
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319 double wigner_start(int m1, int m2, double theta)
00320 {
00321
00322 int i, l, delta;
00323 int cosPower, sinPower;
00324 int absM1, absM2;
00325 double dl, dm1, dm2, normFactor, sinSign;
00326 double dCP, dSP;
00327 double max;
00328 double min;
00329
00330 max = (double) (ABS(m1) > ABS(m2) ? ABS(m1) : ABS(m2));
00331 min = (double) (ABS(m1) < ABS(m2) ? ABS(m1) : ABS(m2));
00332
00333 l = max;
00334 delta = l - min;
00335
00336 absM1 = ABS(m1);
00337 absM2 = ABS(m2);
00338 dl = (double) l;
00339 dm1 = (double) m1;
00340 dm2 = (double) m2;
00341 sinSign = 1.;
00342 normFactor = 1.;
00343
00344 for (i = 0; i < delta; i++)
00345 normFactor *= SQRT((2. * dl - ((double) i)) / (((double) i) + 1.));
00346
00347
00348
00349 normFactor *= SQRT((2. * dl + 1.) / 2.);
00350
00351 if (l == absM1)
00352 if (m1 >= 0)
00353 {
00354 cosPower = l + m2;
00355 sinPower = l - m2;
00356 if ((l - m2) % 2)
00357 sinSign = -1.;
00358 }
00359 else
00360 {
00361 cosPower = l - m2;
00362 sinPower = l + m2;
00363 }
00364 else if (m2 >= 0)
00365 {
00366 cosPower = l + m1;
00367 sinPower = l - m1;
00368 }
00369 else
00370 {
00371 cosPower = l - m1;
00372 sinPower = l + m1;
00373 if ((l + m1) % 2)
00374 sinSign = -1.;
00375 }
00376
00377 dCP = (double) cosPower;
00378 dSP = (double) sinPower;
00379
00380 return normFactor * sinSign * POW(sin(theta / 2), dSP) * POW(cos(theta / 2),
00381 dCP);
00382 }