![]() |
NFFT
3.3.1
|
00001 /* 00002 * Copyright (c) 2002, 2016 Jens Keiner, Stefan Kunis, Daniel Potts 00003 * 00004 * This program is free software; you can redistribute it and/or modify it under 00005 * the terms of the GNU General Public License as published by the Free Software 00006 * Foundation; either version 2 of the License, or (at your option) any later 00007 * version. 00008 * 00009 * This program is distributed in the hope that it will be useful, but WITHOUT 00010 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 00011 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 00012 * details. 00013 * 00014 * You should have received a copy of the GNU General Public License along with 00015 * this program; if not, write to the Free Software Foundation, Inc., 51 00016 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 00017 */ 00018 00019 #include <math.h> 00020 #include <stdio.h> 00021 #include "infft.h" 00022 #include "wigner.h" 00023 #include "infft.h" 00024 00025 double SO3_alpha(const int m1, const int m2, const int j) 00026 { 00027 const int M = MAX(ABS(m1),ABS(m2)), mini = MIN(ABS(m1),ABS(m2)); 00028 00029 if (j < 0) 00030 return K(0.0); 00031 else if (j == 0) 00032 { 00033 if (m1 == 0 && m2 == 0) 00034 return K(1.0); 00035 if (m1 == m2) 00036 return K(0.5); 00037 else 00038 return IF((m1+m2)%2,K(0.0),K(-0.5)); 00039 } 00040 else if (j < M - mini) 00041 return IF(j%2,K(0.5),K(-0.5)); 00042 else if (j < M) 00043 return K(0.5) * SIGNF((R)m1)*SIGNF((R)m2); 00044 else 00045 return 00046 SQRT(((R)(j+1))/((R)(j+1-m1))) 00047 * SQRT(((R)(2*j+1))/((R)(j+1+m1))) 00048 * SQRT(((R)(j+1))/((R)(j+1-m2))) 00049 * SQRT(((R)(2*j+1))/((R)(j+1+m2))); 00050 } 00051 00052 double SO3_beta(const int m1, const int m2, const int j) 00053 { 00054 if (j < 0) 00055 return K(0.0); 00056 else if (j < MAX(ABS(m1),ABS(m2))) 00057 return K(0.5); 00058 else if (m1 == 0 || m2 == 0) 00059 return K(0.0); 00060 else 00061 { 00062 const R m1a = FABS((R)m1), m2a = FABS((R)m2); 00063 return -COPYSIGN( 00064 ((SQRT(m1a)*SQRT(m2a))/((R)j)) 00065 * SQRT(m1a/((R)(j+1-m1))) 00066 * SQRT(((R)(2*j+1))/((R)(j+1+m1))) 00067 * SQRT(m2a/((R)(j+1-m2))) 00068 * SQRT(((R)(2*j+1))/((R)(j+1+m2))), 00069 SIGNF((R)m1)*SIGNF((R)m2)); 00070 } 00071 } 00072 00073 double SO3_gamma(const int m1, const int m2, const int j) 00074 { 00075 if (MAX(ABS(m1),ABS(m2)) < j) 00076 return -(((R)(j+1))/((R)j)) * SQRT((((R)(j-m1))/((R)(j+1-m1))) 00077 *(((R)(j+m1))/((R)(j+1+m1)))*(((R)(j-m2))/((R)(j+1-m2))) 00078 *(((R)(j+m2))/((R)(j+1+m2)))); 00079 else if (j == -1) 00080 return IF(m1 > m2 || !((m1 + m2) % 2), K(1.0), K(-1.0)) 00081 * nfft_lambda2((R)ABS(m2 - m1),(R)ABS(m2 + m1)); 00082 else 00083 return K(0.0); 00084 } 00085 00086 /*compute the coefficients for all degrees*/ 00087 00088 inline void SO3_alpha_row(double *alpha, int N, int k, int m) 00089 { 00090 int j; 00091 double *alpha_act = alpha; 00092 for (j = -1; j <= N; j++) 00093 *alpha_act++ = SO3_alpha(k, m, j); 00094 } 00095 00096 inline void SO3_beta_row(double *beta, int N, int k, int m) 00097 { 00098 int j; 00099 double *beta_act = beta; 00100 for (j = -1; j <= N; j++) 00101 *beta_act++ = SO3_beta(k, m, j); 00102 } 00103 00104 inline void SO3_gamma_row(double *gamma, int N, int k, int m) 00105 { 00106 int j; 00107 double *gamma_act = gamma; 00108 for (j = -1; j <= N; j++) 00109 *gamma_act++ = SO3_gamma(k, m, j); 00110 } 00111 00112 /*compute for all degrees l and orders k*/ 00113 00114 inline void SO3_alpha_matrix(double *alpha, int N, int m) 00115 { 00116 int i, j; 00117 double *alpha_act = alpha; 00118 for (i = -N; i <= N; i++) 00119 { 00120 for (j = -1; j <= N; j++) 00121 { 00122 *alpha_act = SO3_alpha(i, m, j); 00123 alpha_act++; 00124 } 00125 } 00126 } 00127 00128 inline void SO3_beta_matrix(double *alpha, int N, int m) 00129 { 00130 int i, j; 00131 double *alpha_act = alpha; 00132 for (i = -N; i <= N; i++) 00133 { 00134 for (j = -1; j <= N; j++) 00135 { 00136 *alpha_act = SO3_beta(i, m, j); 00137 alpha_act++; 00138 } 00139 } 00140 } 00141 00142 inline void SO3_gamma_matrix(double *alpha, int N, int m) 00143 { 00144 int i, j; 00145 double *alpha_act = alpha; 00146 for (i = -N; i <= N; i++) 00147 { 00148 for (j = -1; j <= N; j++) 00149 { 00150 *alpha_act = SO3_gamma(i, m, j); 00151 alpha_act++; 00152 } 00153 } 00154 } 00155 00156 /*compute all 3termrecurrence coeffs*/ 00157 00158 inline void SO3_alpha_all(double *alpha, int N) 00159 { 00160 int q; 00161 int i, j, m; 00162 double *alpha_act = alpha; 00163 q = 0; 00164 for (m = -N; m <= N; m++) 00165 { 00166 for (i = -N; i <= N; i++) 00167 { 00168 for (j = -1; j <= N; j++) 00169 { 00170 *alpha_act = SO3_alpha(i, m, j); 00171 fprintf(stdout, "alpha_all_%d^[%d,%d]=%f\n", j, i, m, 00172 SO3_alpha(i, m, j)); 00173 alpha_act++; 00174 q = q + 1; 00175 00176 } 00177 } 00178 } 00179 } 00180 00181 inline void SO3_beta_all(double *alpha, int N) 00182 { 00183 int i, j, m; 00184 double *alpha_act = alpha; 00185 for (m = -N; m <= N; m++) 00186 { 00187 for (i = -N; i <= N; i++) 00188 { 00189 for (j = -1; j <= N; j++) 00190 { 00191 *alpha_act = SO3_beta(i, m, j); 00192 alpha_act++; 00193 } 00194 } 00195 } 00196 } 00197 00198 inline void SO3_gamma_all(double *alpha, int N) 00199 { 00200 int i, j, m; 00201 double *alpha_act = alpha; 00202 for (m = -N; m <= N; m++) 00203 { 00204 for (i = -N; i <= N; i++) 00205 { 00206 for (j = -1; j <= N; j++) 00207 { 00208 *alpha_act = SO3_gamma(i, m, j); 00209 alpha_act++; 00210 } 00211 } 00212 } 00213 } 00214 00215 inline void eval_wigner(double *x, double *y, int size, int k, double *alpha, 00216 double *beta, double *gamma) 00217 { 00218 /* Evaluate the wigner function d_{k,nleg} (l,x) for the vector 00219 * of knots x[0], ..., x[size-1] by the Clenshaw algorithm 00220 */ 00221 int i, j; 00222 double a, b, x_val_act, a_old; 00223 double *x_act, *y_act; 00224 double *alpha_act, *beta_act, *gamma_act; 00225 00226 /* Traverse all nodes. */ 00227 x_act = x; 00228 y_act = y; 00229 for (i = 0; i < size; i++) 00230 { 00231 a = 1.0; 00232 b = 0.0; 00233 x_val_act = *x_act; 00234 00235 if (k == 0) 00236 { 00237 *y_act = 1.0; 00238 } 00239 else 00240 { 00241 alpha_act = &(alpha[k]); 00242 beta_act = &(beta[k]); 00243 gamma_act = &(gamma[k]); 00244 for (j = k; j > 1; j--) 00245 { 00246 a_old = a; 00247 a = b + a_old * ((*alpha_act) * x_val_act + (*beta_act)); 00248 b = a_old * (*gamma_act); 00249 alpha_act--; 00250 beta_act--; 00251 gamma_act--; 00252 } 00253 *y_act = (a * ((*alpha_act) * x_val_act + (*beta_act)) + b); 00254 } 00255 x_act++; 00256 y_act++; 00257 } 00258 } 00259 00260 inline int eval_wigner_thresh(double *x, double *y, int size, int k, 00261 double *alpha, double *beta, double *gamma, double threshold) 00262 { 00263 00264 int i, j; 00265 double a, b, x_val_act, a_old; 00266 double *x_act, *y_act; 00267 double *alpha_act, *beta_act, *gamma_act; 00268 00269 /* Traverse all nodes. */ 00270 x_act = x; 00271 y_act = y; 00272 for (i = 0; i < size; i++) 00273 { 00274 a = 1.0; 00275 b = 0.0; 00276 x_val_act = *x_act; 00277 00278 if (k == 0) 00279 { 00280 *y_act = 1.0; 00281 } 00282 else 00283 { 00284 alpha_act = &(alpha[k]); 00285 beta_act = &(beta[k]); 00286 gamma_act = &(gamma[k]); 00287 for (j = k; j > 1; j--) 00288 { 00289 a_old = a; 00290 a = b + a_old * ((*alpha_act) * x_val_act + (*beta_act)); 00291 b = a_old * (*gamma_act); 00292 alpha_act--; 00293 beta_act--; 00294 gamma_act--; 00295 } 00296 *y_act = (a * ((*alpha_act) * x_val_act + (*beta_act)) + b); 00297 if (fabs(*y_act) > threshold) 00298 { 00299 return 1; 00300 } 00301 } 00302 x_act++; 00303 y_act++; 00304 } 00305 return 0; 00306 } 00307 00308 /************************************************************************/ 00309 /* L2 normed wigner little d, WHERE THE DEGREE OF THE FUNCTION IS EQUAL 00310 TO ONE OF ITS ORDERS. This is the function to use when starting the 00311 three-term recurrence at orders (m1,m2) 00312 00313 Note that, by definition, since I am starting the recurrence with this 00314 function, that the degree j of the function is equal to max(abs(m1), abs(m2) ). 00315 */ 00316 00317 double wigner_start(int m1, int m2, double theta) 00318 { 00319 00320 int i, l, delta; 00321 int cosPower, sinPower; 00322 int absM1, absM2; 00323 double dl, dm1, dm2, normFactor, sinSign; 00324 double dCP, dSP; 00325 double max; 00326 double min; 00327 00328 max = (double) (ABS(m1) > ABS(m2) ? ABS(m1) : ABS(m2)); 00329 min = (double) (ABS(m1) < ABS(m2) ? ABS(m1) : ABS(m2)); 00330 00331 l = max; 00332 delta = l - min; 00333 00334 absM1 = ABS(m1); 00335 absM2 = ABS(m2); 00336 dl = (double) l; 00337 dm1 = (double) m1; 00338 dm2 = (double) m2; 00339 sinSign = 1.; 00340 normFactor = 1.; 00341 00342 for (i = 0; i < delta; i++) 00343 normFactor *= SQRT((2. * dl - ((double) i)) / (((double) i) + 1.)); 00344 00345 /* need to adjust to make the L2-norm equal to 1 */ 00346 00347 normFactor *= SQRT((2. * dl + 1.) / 2.); 00348 00349 if (l == absM1) 00350 if (m1 >= 0) 00351 { 00352 cosPower = l + m2; 00353 sinPower = l - m2; 00354 if ((l - m2) % 2) 00355 sinSign = -1.; 00356 } 00357 else 00358 { 00359 cosPower = l - m2; 00360 sinPower = l + m2; 00361 } 00362 else if (m2 >= 0) 00363 { 00364 cosPower = l + m1; 00365 sinPower = l - m1; 00366 } 00367 else 00368 { 00369 cosPower = l - m1; 00370 sinPower = l + m1; 00371 if ((l + m1) % 2) 00372 sinSign = -1.; 00373 } 00374 00375 dCP = (double) cosPower; 00376 dSP = (double) sinPower; 00377 00378 return normFactor * sinSign * POW(sin(theta / 2), dSP) * POW(cos(theta / 2), 00379 dCP); 00380 }