NFFT Logo 3.2.2
wigner.c
00001 /*
00002  * Copyright (c) 2002, 2012 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 /* $Id: wigner.c 3896 2012-10-10 12:19:26Z tovo $ */
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 /*compute the coefficients for all degrees*/
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 /*compute for all degrees l and orders k*/
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 /*compute all 3termrecurrence coeffs*/
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   /* Evaluate the wigner function d_{k,nleg} (l,x) for the vector
00221    * of knots  x[0], ..., x[size-1] by the Clenshaw algorithm
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   /* Traverse all nodes. */
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   /* Traverse all nodes. */
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 /* L2 normed wigner little d, WHERE THE DEGREE OF THE FUNCTION IS EQUAL
00312  TO ONE OF ITS ORDERS. This is the function to use when starting the
00313  three-term recurrence at orders (m1,m2)
00314 
00315  Note that, by definition, since I am starting the recurrence with this
00316  function, that the degree j of the function is equal to max(abs(m1), abs(m2) ).
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   /* need to adjust to make the L2-norm equal to 1 */
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 }

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409