--- trunk/SHAPES/SHFunc.cpp	2004/06/23 20:53:17	1290
+++ trunk/SHAPES/SHFunc.cpp	2004/06/23 22:43:54	1291
@@ -1,3 +1,5 @@
+#include <stdio.h>
+#include <cmath>
 #include "SHFunc.hpp"
 
 SHFunc::SHFunc() {
@@ -5,7 +7,124 @@ double SHFunc::getValueAt(double costheta, double phi)
 
 double SHFunc::getValueAt(double costheta, double phi) {
 
+  double f, p, phase;
+
+  // incredibly inefficient way to get the normalization, but
+  // we use a lookup table in the factorial code:
+
+  // normalization factor:
+  f = sqrt( (2*L+1)/(4.0*M_PI) * Fac(L-M) / Fac(L+M) );
+  // associated Legendre polynomial
+  p = LegendreP(L,M,costheta);
   
+  if (funcType == SH_SIN) {
+    phase = sin((double)M * phi);
+  } else {
+    phase = cos((double)M * phi);
+  }
 
+  
+  return coefficient*f*p*phase;
 
 }
+//-----------------------------------------------------------------------------//
+//
+// double LegendreP (int l, int m, double x);
+//
+// Computes the value of the associated Legendre polynomial P_lm (x)
+// of order l at a given point.
+//
+// Input:
+//   l  = degree of the polynomial  >= 0
+//   m  = parameter satisfying 0 <= m <= l,
+//   x  = point in which the computation is performed, range -1 <= x <= 1.
+// Returns:
+//   value of the polynomial in x
+//
+//-----------------------------------------------------------------------------//
+
+double SHFunc::LegendreP (int l, int m, double x) {
+  // check parameters
+  if (m < 0 || m > l || fabs(x) > 1.0) {
+    printf("LegendreP got a bad argument\n");
+    return NAN;
+  }
+  
+  double pmm = 1.0;
+  if (m > 0) {
+    double h = sqrt((1.0-x)*(1.0+x)),
+      f = 1.0;
+    for (int i = 1; i <= m; i++) {
+      pmm *= -f * h;
+      f += 2.0;
+    }
+  }
+  if (l == m)
+    return pmm;
+  else {
+    double pmmp1 = x * (2 * m + 1) * pmm;
+    if (l == (m+1))
+      return pmmp1;
+    else {
+      double pll = 0.0;
+      for (int ll = m+2; ll <= l; ll++) {
+        pll = (x * (2 * ll - 1) * pmmp1 - (ll + m - 1) * pmm) / (ll - m);
+        pmm = pmmp1;
+        pmmp1 = pll;
+      }
+      return pll;
+    }
+  }
+}
+
+double SHFunc::Fac (int n) {
+
+  static double facn[31] = {
+    1.0,
+    1.0,
+    2.0,
+    6.0,
+    24.0,
+    120.0,
+    720.0,
+    5040.0,
+    40320.0,
+    362880.0,
+    3628800.0,
+    39916800.0,
+    479001600.0,
+    6227020800.0,
+    87178291200.0,
+    1.307674368e12,
+    2.0922789888e13,
+    3.55687428096e14,
+    6.402373705728e15,
+    1.21645100408832e17,
+    2.43290200817664e18,
+    5.109094217170944e19,
+    1.12400072777760768e21,
+    2.585201673888497664e22,
+    6.2044840173323943936e23,
+    1.5511210043330985984e25,
+    4.03291461126605635584e26,
+    1.0888869450418352160768e28,
+    3.04888344611713860501504e29,
+    8.841761993739701954543616e30,
+    2.6525285981219105863630848e32
+  };
+  
+  
+  static int nmax = 0;
+  static double xmin, xmax;
+    
+  if (n < 0) {
+    printf("factorial of negative integer undefined\n");
+    return NAN;
+  }
+  
+  if (n <= 30) return facn[n];
+  else {
+    printf("n is so large that Fac(n) will overflow\n");
+    return NAN;
+  }  
+}