--- trunk/SHAPES/SHFunc.cpp	2004/06/23 20:53:17	1290
+++ trunk/SHAPES/SHFunc.cpp	2004/06/24 15:31:52	1295
@@ -1,3 +1,5 @@
+#include <stdio.h>
+#include <cmath>
 #include "SHFunc.hpp"
 
 SHFunc::SHFunc() {
@@ -5,7 +7,81 @@ 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
+
+  // normalization factor:
+  f = sqrt( (2*L+1)/(4.0*M_PI) * Fact(L-M) / Fact(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::Fact(double n) {
+
+  if (n < 0.0) return NAN;
+  else {
+    if (n < 2.0) return 1.0;
+    else
+      return n*Fact(n-1.0);
+  }
+  
+}