--- trunk/SHAPES/SphereHarm.cpp	2004/06/24 15:47:52	1298
+++ trunk/SHAPES/SphereHarm.cpp	2004/06/26 15:40:49	1313
@@ -119,7 +119,7 @@ void SphereHarm::printToShapesFile(char name[200], int
   shapes << "end ShapeInfo\n";
 }
 
-void SphereHarm::printToShapesFile(char name[200], int index){
+void SphereHarm::printToShapesFile(char name[200], int index, double tolVal){
   ofstream shapes(name, ios::app);
   
   biggest = 0.0;
@@ -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,15 +149,17 @@ 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++;
-      if (fabs(sm) > 0.01 * biggest) nfuncs++;
+      if (fabs(cm) > tolVal * biggest) nfuncs++;
+      if (fabs(sm) > tolVal * biggest) nfuncs++;
     }
   }
     
@@ -180,16 +184,18 @@ 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) 
+      if (fabs(cm) > tolVal * biggest) 
         shapes << l << "\t" << m << "\tcos\t\t" << cm << "\n";
-      if (fabs(sm) > 0.01 * biggest) 
+      if (fabs(sm) > tolVal * biggest) 
         shapes << l << "\t" << m << "\tsin\t\t" << sm << "\n";
     }
   }
@@ -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);
+  }
+}