Comments only

Author: Martin D. Weinberg

expui/BiorthBasis.cc

@@ -601,7 +601,7 @@ namespace BasisClasses
     }
 
     double densfac = 1.0/(scale*scale*scale) * 0.25/M_PI;
-    double potlfac = G/scale;
+    double potlfac = G/scale;	// Grav constant G is 1 by default
 
     return
       {den0 * densfac,		// 0
@@ -707,7 +707,7 @@ namespace BasisClasses
       }
     }
 
-    double potlfac = G/scale;
+    double potlfac = G/scale;	// Grav constant G is 1 by default
 
     potr *= (-potlfac)/scale;
     pott *= (-potlfac);
@@ -1494,8 +1494,8 @@ namespace BasisClasses
 
     sl->accumulated_eval(R, z, phi, tpotl0, tpotl, tpotR, tpotz, tpotp);
 
-    tpotl0 *= G;
-    tpotl  *= G;
+    tpotl0 *= G;		// Apply gravitational constant to
+    tpotl  *= G;		// potential and forces
     tpotR  *= G;
     tpotz  *= G;
     tpotp  *= G;
@@ -1520,7 +1520,7 @@ namespace BasisClasses
 
     sl->accumulated_eval(R, z, phi, tpotl0, tpotl, tpotR, tpotz, tpotp);
 
-    tpotl0 *= G;
+    tpotl0 *= G;		// Apply G to potential and forces
     tpotl  *= G;
     tpotR  *= G;
     tpotz  *= G;
@@ -1551,6 +1551,7 @@ namespace BasisClasses
     double tpotx = tpotR*x/R - tpotp*y/R ;
     double tpoty = tpotR*y/R + tpotp*x/R ;
 
+    // Apply G to forces on return
     return {tpotx*G, tpoty*G, tpotz*G};
   }
 
@@ -1562,7 +1563,7 @@ namespace BasisClasses
     sl->accumulated_eval(R, z, phi, tpotl0, tpotl, tpotR, tpotz, tpotp);
     tdens = sl->accumulated_dens_eval(R, z, phi, tdens0);
 
-    tpotl0 *= G;
+    tpotl0 *= G;		// Apply G to potential and forces
     tpotl  *= G;
     tpotR  *= G;
     tpotz  *= G;
@@ -2053,6 +2054,7 @@ namespace BasisClasses
     if (R>ortho->getRtable() or fabs(z)>ortho->getRtable()) {
       double r2 = R*R + z*z;
       double r  = sqrt(r2);
+      // Apply G to potential and forces
       pot0 = -G*totalMass/r;
       rpot = -G*totalMass*R/(r*r2 + 10.0*std::numeric_limits<double>::min());
       zpot = -G*totalMass*z/(r*r2 + 10.0*std::numeric_limits<double>::min());
@@ -2128,7 +2130,7 @@ namespace BasisClasses
 
     den0 *= -1.0;
     den1 *= -1.0;
-    pot0 *= -G;
+    pot0 *= -G;			// Apply G to potential and forces
     pot1 *= -G;
     rpot *= -G;
     zpot *= -G;
@@ -2158,6 +2160,7 @@ namespace BasisClasses
       double r2 = R*R + z*z;
       double r  = sqrt(r2);
 
+      // Apply G to forces
       rpot = -G*totalMass*R/(r*r2 + 10.0*std::numeric_limits<double>::min());
       zpot = -G*totalMass*z/(r*r2 + 10.0*std::numeric_limits<double>::min());
 
@@ -2218,7 +2221,7 @@ namespace BasisClasses
       }
     }
 
-    rpot *= -G;
+    rpot *= -G;			// Apply G to forces
     zpot *= -G;
     ppot *= -G;
 
@@ -2910,7 +2913,7 @@ namespace BasisClasses
 
     den0 *= -1.0;
     den1 *= -1.0;
-    pot0 *= -G;
+    pot0 *= -G;			// Apply G to potential and forces
     pot1 *= -G;
     rpot *= -G;
     ppot *= -G;
@@ -2977,10 +2980,10 @@ namespace BasisClasses
       }
     }
 
+    // Apply G to forces
     rpot *= -G;
     ppot *= -G;
 
-
     double potx = rpot*x/R - ppot*y/R;
     double poty = rpot*y/R + ppot*x/R;
 
@@ -3472,6 +3475,7 @@ namespace BasisClasses
       }
     }
 
+    // Apply G to forces on return
     return {G*accx.real(), G*accy.real(), G*accz.real()};
   }
 
@@ -3483,6 +3487,7 @@ namespace BasisClasses
 
     auto [pot, den, frcx, frcy, frcz] = eval(x, y, z);
 
+    // Apply G to potential and forces on return
     return {0, den, den, 0, pot*G, pot*G, frcx*G, frcy*G, frcz*G};
   }
 
@@ -3500,10 +3505,11 @@ namespace BasisClasses
     double potp = -frcx*sin(phi) + frcy*cos(phi);
     double potz =  frcz;
 
-    potR *= -G;
+    potR *= -G;			// Apply G to forces
     potp *= -G;
     potz *= -G;
 
+    // Apply Go to potential on return
     return {0, den, den, 0, pot*G, pot*G, potR, potz, potp};
   }
 
@@ -3522,10 +3528,11 @@ namespace BasisClasses
     double pott =  frcx*cos(phi)*costh + frcy*sin(phi)*costh - frcz*sinth;
     double potp = -frcx*sin(phi)       + frcy*cos(phi);
 
-    potr *= -G;
+    potr *= -G;			// Apply G to forces
     pott *= -G;
     potp *= -G;
 
+    // Apply G to potential on return
     return {0, den, den, 0, pot*G, pot*G, potr, pott, potp};
   }
 
@@ -3841,6 +3848,7 @@ namespace BasisClasses
     double frcy = -frc(1).real();
     double frcz = -frc(2).real();
 
+    // Apply G to potential and forces on return
     return {0, den1, den1, 0, pot1*G, pot1*G, frcx*G, frcy*G, frcz*G};
   }
 
@@ -3855,6 +3863,7 @@ namespace BasisClasses
     // Get the basis fields
     auto frc = ortho->get_force(expcoef, pos);
 
+    // Apply G to forces on return
     return {-G*frc(0).real(), -G*frc(1).real(), -G*frc(2).real()};
   }
 
@@ -3871,10 +3880,12 @@ namespace BasisClasses
 
     // Get the basis fields
     double den1 = ortho->get_dens(expcoef, pos).real();
-    double pot1 = ortho->get_pot (expcoef, pos).real();
+    double pot1 = ortho->get_pot (expcoef, pos).real() * G;
 
-    auto frc = ortho->get_force(expcoef, pos)*G;
+    auto frc = ortho->get_force(expcoef, pos) * G;
 
+    // Gravitational constant G applied to potenial and forces above
+
     double frcx = frc(0).real(), frcy = frc(1).real(), frcz = frc(2).real();
 
     double potR =  frcx*cos(phi) + frcy*sin(phi);
@@ -3904,7 +3915,9 @@ namespace BasisClasses
     double den1 = ortho->get_dens(expcoef, pos).real();
     double pot1 = ortho->get_pot (expcoef, pos).real() * G;
 
-    auto frc = ortho->get_force(expcoef, pos);
+    auto frc = ortho->get_force(expcoef, pos) * G;
+
+    // Gravitational constant G applied to potential and forces above
 
     double frcx = frc(0).real();
     double frcy = frc(1).real();
@@ -3914,9 +3927,9 @@ namespace BasisClasses
     double pott =  frcx*cos(phi)*costh + frcy*sin(phi)*costh - frcz*sinth;
     double potp = -frcx*sin(phi)       + frcy*cos(phi);
 
-    potr *= -G;
-    pott *= -G;
-    potp *= -G;
+    potr *= -1.0;
+    pott *= -1.0;
+    potp *= -1.0;
 
     return {0, den1, den1, 0, pot1, pot1, potr, pott, potp};
   }