--- trunk/SHAPES/SphereHarm.cpp	2004/06/23 20:48:08	1289
+++ trunk/SHAPES/SphereHarm.cpp	2004/07/20 20:02:15	1360
@@ -109,15 +109,30 @@ void SphereHarm::printShapesFileStart(char name[200], 
                 weights );
 }
 
-void SphereHarm::printShapesFileStart(char name[200], string particle,
+void SphereHarm::printShapesFileStart(char name[200], char particle[80],
                                       double mass, double momInert[3][3]){
-  ofstream shapes(name);
+  ofstream myShapes(name);
+  printShapesStreamStart(myShapes, particle, mass, momInert);
+}
+
+void SphereHarm::printShapesStreamStart(ostream& shapes, char particle[80],
+                                      double mass, double momInert[3][3]){
+  shapes << "begin ShapeInfo\n";
+  shapes << "#name\t\tmass\tI_xx\tI_yy\tI_zz\n";
   shapes << particle << "\t" << mass << "\t" << momInert[0][0] << "\t" 
-         << momInert[1][1] << "\t" << momInert[2][2] << "\n\n";
+         << momInert[1][1] << "\t" << momInert[2][2] << "\n";
+  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);
+
+  printToShapesStream(shapes, index, tolVal);
+}
+
+void SphereHarm::printToShapesStream(ostream& shapes, int index, double tolVal) {
   
   biggest = 0.0;
   nfuncs = 0;
@@ -127,11 +142,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);
@@ -144,15 +161,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++;
     }
   }
     
@@ -177,29 +196,54 @@ 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) 
-        shapes << l << "\t" << m << "\tcos\t" << cm << "\n";
-      if (fabs(sm) > 0.01 * biggest) 
-        shapes << l << "\t" << m << "\tsin\t" << sm << "\n";
+      if (fabs(cm) > tolVal * biggest) 
+        shapes << l << "\t" << m << "\tcos\t\t" << cm << "\n";
+      if (fabs(sm) > tolVal * biggest) 
+        shapes << l << "\t" << m << "\tsin\t\t" << sm << "\n";
     }
   }
   switch(index){
     case 0:{
-      shapes << "\nend ContactFunctions\n";
+      shapes << "end ContactFunctions\n";
     }; break;
     case 1:{
-      shapes << "\nend RangeFunctions\n";
+      shapes << "end RangeFunctions\n";
     }; break;
     case 2:{
-      shapes << "\nend StrengthFunctions\n";
+      shapes << "end StrengthFunctions\n";
     }; break;
   }
 }
 
+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);
+  }
+}