Fix round off boundary and never ignore nout

Author: Martin D. Weinberg

expui/BiorthBasis.H

@@ -1189,7 +1189,7 @@ namespace BasisClasses
   std::tuple<Eigen::VectorXd, Eigen::Tensor<float, 3>>
   IntegrateOrbits (double tinit, double tfinal, double h,
 		   Eigen::MatrixXd ps, std::vector<BasisCoef> bfe,
-		   AccelFunctor F, int nout=std::numeric_limits<int>::max());
+		   AccelFunctor F, int nout=0);
 
   using BiorthBasisPtr = std::shared_ptr<BiorthBasis>;
 }

expui/BiorthBasis.cc

@@ -3914,23 +3914,33 @@ namespace BasisClasses
 
     // Number of steps
     //
-    int numT = std::ceil( (tfinal - tinit)/h );
+    int numT = std::ceil( (tfinal - tinit)/h + 0.5);
 
     // Want both end points in the output at minimum
     //
     numT = std::max(2, numT);
-    nout = std::max(2, nout);
 
-    // Number of output steps.  Compute a stride>1 if nout<numT.
+    // Number of output steps
     //
-    int stride = std::ceil(static_cast<double>(numT)/static_cast<double>(nout));
-    if (stride>1) numT = (nout-1) * stride + 1;
-    else          nout = numT;
+    int stride = 1;		// Default stride
+    if (nout>0) {		// User has specified output count...
+      nout = std::max(2, nout);
+      stride = std::ceil(static_cast<double>(numT)/static_cast<double>(nout));
+      numT = (nout-1) * stride + 1;
+    } else {			// Otherwise, use the default output number
+      nout = numT;		// with the default stride
+    }
 
     // Compute the interval-matching step
     //
     h = (tfinal - tinit)/(numT-1);
 
+    // DEBUG
+    if (false) 
+      std::cout << "BasisClasses::IntegrateOrbits: choosing nout=" << nout
+		<< " numT=" << numT << " h=" << h << " stride=" << stride
+		<< std::endl;
+
     // Return data
     //
     Eigen::Tensor<float, 3> ret;

pyEXP/BasisWrappers.cc

@@ -2194,11 +2194,13 @@ void BasisFactoryClasses(py::module &m)
 	R"(
         Compute particle orbits in gravitational field from the bases
 
-        Integrate a list of initial conditions from tinit to tfinal with a
-	step size of h using the list of basis and coefficient pairs. Every
-	step will be included in return unless you provide an explicit
-	value for 'nout', the number of desired output steps.  This will
-	choose the 'nout' points closest to the desired time.
+        Integrate a list of initial conditions from 'tinit' to 'tfinal' with
+	a step size of 'h' using the list of basis and coefficient pairs. The
+        step size will be adjusted to provide uniform sampling.  Every
+	step will be returned unless you provide an explicit value for 'nout',
+        the number of desired output steps.  In this case, the code will
+        choose new set step size equal or smaller to the supplied step size
+        with a stride to provide exactly 'nout' output times.
 
         Parameters
         ----------
@@ -2216,7 +2218,7 @@ void BasisFactoryClasses(py::module &m)
         func : AccelFunctor
             the force function
         nout : int 
-            the number of output intervals if (tfinal - tinit) / h > nout
+            the number of output points, if specified
 
         Returns
         -------
@@ -2225,5 +2227,5 @@ void BasisFactoryClasses(py::module &m)
         )",
 	py::arg("tinit"), py::arg("tfinal"), py::arg("h"),
 	py::arg("ps"), py::arg("basiscoef"), py::arg("func"),
-	py::arg("nout")=std::numeric_limits<int>::max());
+	py::arg("nout")=0);
 }