NFFT Logo 3.2.2
legendre.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: legendre.c 3896 2012-10-10 12:19:26Z tovo $ */
00020 
00021 #include <math.h>
00022 #include <stdio.h>
00023 #include "infft.h"
00024 #include "nfft3util.h"
00025 #include "legendre.h"
00026 
00027 /* One over sqrt(pi) */
00028 DK(KSQRTPII,0.56418958354775628694807945156077258584405062932900);
00029 
00030 static inline R alpha_al(const int k, const int n)
00031 {
00032   if (k > 0)
00033   {
00034     if (k < n)
00035       return IF(k%2,K(1.0),K(-1.0));
00036     else
00037       return SQRT(((R)(2*k+1))/((R)(k-n+1)))*SQRT((((R)(2*k+1))/((R)(k+n+1))));
00038   }
00039   else if (k == 0)
00040   {
00041     if (n == 0)
00042       return K(1.0);
00043     else
00044       return IF(n%2,K(0.0),K(-1.0));
00045   }
00046   return K(0.0);
00047 }
00048 
00049 static inline R beta_al(const int k, const int n)
00050 {
00051   if (0 <= k && k < n)
00052     return K(1.0);
00053   else
00054     return K(0.0);
00055 }
00056 
00057 static inline R gamma_al(const int k, const int n)
00058 {
00059   if (k == -1)
00060     return SQRT(KSQRTPII*nfft_lambda((R)(n),K(0.5)));
00061   else if (k <= n)
00062     return K(0.0);
00063   else
00064     return -SQRT(((R)(k-n))/((R)(k-n+1))*((R)(k+n))/((R)(k+n+1)));
00065 }
00066 
00067 void alpha_al_row(R *alpha, const int N, const int n)
00068 {
00069   int j;
00070   R *p = alpha;
00071   for (j = -1; j <= N; j++)
00072     *p++ = alpha_al(j,n);
00073 }
00074 
00075 void beta_al_row(R *beta, const int N, const int n)
00076 {
00077   int j;
00078   R *p = beta;
00079   for (j = -1; j <= N; j++)
00080     *p++ = beta_al(j,n);
00081 }
00082 
00083 void gamma_al_row(R *gamma, const int N, const int n)
00084 {
00085   int j;
00086   R *p = gamma;
00087   for (j = -1; j <= N; j++)
00088     *p++ = gamma_al(j,n);
00089 }
00090 
00091 inline void alpha_al_all(R *alpha, const int N)
00092 {
00093   int i,j;
00094   R *p = alpha;
00095   for (i = 0; i <= N; i++)
00096     for (j = -1; j <= N; j++)
00097       *p++ = alpha_al(j,i);
00098 }
00099 
00100 inline void beta_al_all(R *alpha, const int N)
00101 {
00102   int i,j;
00103   R *p = alpha;
00104   for (i = 0; i <= N; i++)
00105     for (j = -1; j <= N; j++)
00106       *p++ = beta_al(j,i);
00107 }
00108 
00109 inline void gamma_al_all(R *alpha, const int N)
00110 {
00111   int i,j;
00112   R *p = alpha;
00113   for (i = 0; i <= N; i++)
00114     for (j = -1; j <= N; j++)
00115       *p++ = gamma_al(j,i);
00116 }
00117 
00118 void eval_al(R *x, R *y, const int size, const int k, R *alpha,
00119   R *beta, R *gamma)
00120 {
00121   /* Evaluate the associated Legendre polynomial P_{k,nleg} (l,x) for the vector
00122    * of knots  x[0], ..., x[size-1] by the Clenshaw algorithm
00123    */
00124   int i,j;
00125   R a,b,x_val_act,a_old;
00126   R *x_act, *y_act;
00127   R *alpha_act, *beta_act, *gamma_act;
00128 
00129   /* Traverse all nodes. */
00130   x_act = x;
00131   y_act = y;
00132   for (i = 0; i < size; i++)
00133   {
00134     a = 1.0;
00135     b = 0.0;
00136     x_val_act = *x_act;
00137 
00138     if (k == 0)
00139     {
00140       *y_act = 1.0;
00141     }
00142     else
00143     {
00144       alpha_act = &(alpha[k]);
00145       beta_act = &(beta[k]);
00146       gamma_act = &(gamma[k]);
00147       for (j = k; j > 1; j--)
00148       {
00149         a_old = a;
00150         a = b + a_old*((*alpha_act)*x_val_act+(*beta_act));
00151          b = a_old*(*gamma_act);
00152         alpha_act--;
00153         beta_act--;
00154         gamma_act--;
00155       }
00156       *y_act = (a*((*alpha_act)*x_val_act+(*beta_act))+b);
00157     }
00158     x_act++;
00159     y_act++;
00160   }
00161 }
00162 
00163 int eval_al_thresh(R *x, R *y, const int size, const int k, R *alpha,
00164   R *beta, R *gamma, R threshold)
00165 {
00166   /* Evaluate the associated Legendre polynomial P_{k,nleg} (l,x) for the vector
00167    * of knots  x[0], ..., x[size-1] by the Clenshaw algorithm
00168    */
00169   int i,j;
00170   R a,b,x_val_act,a_old;
00171   R *x_act, *y_act;
00172   R *alpha_act, *beta_act, *gamma_act;
00173 
00174   /* Traverse all nodes. */
00175   x_act = x;
00176   y_act = y;
00177   for (i = 0; i < size; i++)
00178   {
00179     a = 1.0;
00180     b = 0.0;
00181     x_val_act = *x_act;
00182 
00183     if (k == 0)
00184     {
00185      *y_act = 1.0;
00186     }
00187     else
00188     {
00189       alpha_act = &(alpha[k]);
00190       beta_act = &(beta[k]);
00191       gamma_act = &(gamma[k]);
00192       for (j = k; j > 1; j--)
00193       {
00194         a_old = a;
00195         a = b + a_old*((*alpha_act)*x_val_act+(*beta_act));
00196          b = a_old*(*gamma_act);
00197         alpha_act--;
00198         beta_act--;
00199         gamma_act--;
00200       }
00201       *y_act = (a*((*alpha_act)*x_val_act+(*beta_act))+b);
00202       if (fabs(*y_act) > threshold)
00203       {
00204         return 1;
00205       }
00206     }
00207     x_act++;
00208     y_act++;
00209   }
00210   return 0;
00211 }

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409