NFFT  3.3.1
wigner.c
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 }