--- trunk/SHAPES/SphereHarm.cpp	2004/06/25 18:58:12	1308
+++ trunk/SHAPES/SphereHarm.cpp	2004/06/25 20:10:53	1309
@@ -130,11 +130,13 @@ void SphereHarm::printToShapesFile(char name[200], int
       dummy2 = seanindex(-m, l, bw);
         
       if (m == 0) {
-        cm = rcoeffs[dummy1];
-        sm = icoeffs[dummy1];
+        cm = normFactor(l,m)*rcoeffs[dummy1];
+        sm = normFactor(l,m)*icoeffs[dummy1];
       } else {
-        cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2];
-        sm = pow(-1.0,(double)m)*icoeffs[dummy1] - icoeffs[dummy2];
+        cm = normFactor(l,m)*(pow(-1.0,(double)m)*rcoeffs[dummy1] 
+                              + rcoeffs[dummy2]);
+        sm = normFactor(l,m)*(pow(-1.0,(double)m)*icoeffs[dummy1] 
+                              - icoeffs[dummy2]);
       }
         
       if (fabs(cm) > biggest) biggest = fabs(cm);
@@ -147,11 +149,13 @@ void SphereHarm::printToShapesFile(char name[200], int
       dummy2 = seanindex(-m, l, bw);
         
       if (m == 0) {
-        cm = rcoeffs[dummy1];
-        sm = icoeffs[dummy1];
+        cm = normFactor(l,m)*rcoeffs[dummy1];
+        sm = normFactor(l,m)*icoeffs[dummy1];
       } else {
-        cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2];
-        sm = pow(-1.0,(double)m)*icoeffs[dummy1] - icoeffs[dummy2];
+        cm = normFactor(l,m)*(pow(-1.0,(double)m)*rcoeffs[dummy1] 
+                              + rcoeffs[dummy2]);
+        sm = normFactor(l,m)*(pow(-1.0,(double)m)*icoeffs[dummy1] 
+                              - icoeffs[dummy2]);
       }
         
       if (fabs(cm) > 0.01 * biggest) nfuncs++;
@@ -180,11 +184,13 @@ void SphereHarm::printToShapesFile(char name[200], int
       dummy2 = seanindex(-m, l, bw);
         
       if (m == 0) {
-        cm = rcoeffs[dummy1];
-        sm = icoeffs[dummy1];
+        cm = normFactor(l,m)*rcoeffs[dummy1];
+        sm = normFactor(l,m)*icoeffs[dummy1];
       } else {
-        cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2];
-        sm = pow(-1.0,(double)m)*icoeffs[dummy1] - icoeffs[dummy2];
+        cm = normFactor(l,m)*(pow(-1.0,(double)m)*rcoeffs[dummy1] 
+                              + rcoeffs[dummy2]);
+        sm = normFactor(l,m)*(pow(-1.0,(double)m)*icoeffs[dummy1] 
+                              - icoeffs[dummy2]);
       }
         
       if (fabs(cm) > 0.01 * biggest) 
@@ -206,3 +212,26 @@ void SphereHarm::printToShapesFile(char name[200], int
   }
 }
 
+double SphereHarm::normFactor(int L, int M){
+  // normalization factor:
+  if (L+M > 170){
+    printf("Warning: A coefficient was omitted because l + m > 170.\n"
+           "\tThe double buffer overflows with factorial calculations\n"
+           "\tof 170 and higher.  You should consider using a smaller\n"
+           "\tbandwidth if you aren't okay with the loss of the %i, %i\n"
+           "\tspherical harmonic.\n", L, M);
+    return 0.0;
+  }
+  else
+    return sqrt( (2*L+1)/(4.0*M_PI)*factorialFunc((double)(L-M)) 
+                 / factorialFunc(double(L+M)) );
+}
+
+double SphereHarm::factorialFunc(double n) {
+  if (n < 0.0) return NAN;
+  else {
+    if (n < 2.0) return 1.0;
+    else
+      return n*factorialFunc(n-1.0);
+  }
+}