Added comments on normalization

Author: Martin D. Weinberg

exputil/OrthoFunction.cc

@@ -54,6 +54,7 @@ double OrthoFunction::scalar_prod(int n, int m)
   return ans;
 }
 
+// Compute the unnormalized polynomial
 Eigen::VectorXd OrthoFunction::poly_eval(double t, int n)
 {
   Eigen::VectorXd ret(n+1);
@@ -68,6 +69,8 @@ Eigen::VectorXd OrthoFunction::poly_eval(double t, int n)
   return ret;
 }
 
+// Compute the orthogonality matrix.  A perfect result is the unit
+// matrix.
 Eigen::MatrixXd OrthoFunction::testOrtho()
 {
   Eigen::MatrixXd ret(nmax+1, nmax+1);
@@ -80,10 +83,10 @@ Eigen::MatrixXd OrthoFunction::testOrtho()
     double f = dx*lq->weight(i)*d_r_to_x(x)*pow(r, dof-1);
     double y = 2.0*r/scale; if (segment) y = x;
 
-    // Evaluate polynomial
+    // Evaluate the unnormalized polynomial
     Eigen::VectorXd p = poly_eval(y, nmax) * W(r);
 
-    // Compute scalar products
+    // Compute scalar products with the normalizations
     for (int n1=0; n1<=nmax; n1++) {
       for (int n2=0; n2<=nmax; n2++) {
 	ret(n1, n2) += f * p(n1) * p(n2) / sqrt(norm[n1]*norm[n2]);
@@ -110,7 +113,7 @@ void OrthoFunction::dumpOrtho(const std::string& filename)
       double r = x_to_r(x);
       double y = 2.0*r/scale; if (segment) y = x;
 
-      // Evaluate polynomial
+      // Evaluate polynomial and apply the normalization
       //
       Eigen::VectorXd p = poly_eval(y, nmax) * W(r);
       fout << std::setw(16) << r;
@@ -132,9 +135,11 @@ Eigen::VectorXd OrthoFunction::operator()(double r)
   double y = 2.0*r/scale;
   if (segment) y = r_to_x(r);
 
-  // Generate normalized orthogonal functions
+  // Generate normalized orthogonal functions with weighting
   auto p = poly_eval(y, nmax);
   for (int i=0; i<=nmax; i++) p(i) *= w/sqrt(norm[i]);
 
+  // The inner product of p(i), p(j) is Kronecker delta(i, j)
+
   return p;
 }