diff --git a/exputil/EXPmath.cc b/exputil/EXPmath.cc
index 74189d649..feed74357 100644
--- a/exputil/EXPmath.cc
+++ b/exputil/EXPmath.cc
@@ -61,12 +61,12 @@ namespace AltMath
     return ans;
   }
 
-#define ACC 40.0
-#define BIGNO 1.0e10
-#define BIGNI 1.0e-10
-
   double bessj(int n,double x)
   {
+    const double ACC = 40.0;
+    const double BIGNO = 1.0e10;
+    const double BIGNI = 1.0e-10;
+
     int j,jsum,m;
     double ax,bj,bjm,bjp,sum,tox,ans;
     
@@ -114,10 +114,6 @@ namespace AltMath
     return  x < 0.0 && n%2 == 1 ? -ans : ans;
   }
   
-#undef ACC
-#undef BIGNO
-#undef BIGNI
-
   double bessk0(double x)
   {
     double y,ans;
@@ -174,12 +170,12 @@ namespace AltMath
     return ans;
   }
 
-#define ACC 40.0
-#define BIGNO 1.0e10
-#define BIGNI 1.0e-10
-
   double bessi(int n,double x)
   {
+    const double ACC = 40.0;
+    const double BIGNO = 1.0e10;
+    const double BIGNI = 1.0e-10;
+
     int j;
     double bi,bim,bip,tox,ans;
     double bessi0(double );
@@ -207,10 +203,6 @@ namespace AltMath
     }
   }
 
-#undef ACC
-#undef BIGNO
-#undef BIGNI
-
   double bessi0(double x)
   {
     double ax,ans;
@@ -253,12 +245,12 @@ namespace AltMath
     return x < 0.0 ? -ans : ans;
   }
 
-#define ACC 40.0
-#define BIGNO 1.0e10
-#define BIGNI 1.0e-10
-
   double jn_sph(int n, double x)
   {
+    const double ACC = 40.0;
+    const double BIGNO = 1.0e10;
+    const double BIGNI = 1.0e-10;
+
     int j,m;
     double ax,bj,bjm,bjp,tox,ans;
     double jn_sph0(double x),jn_sph1(double x),jn_sphm1(double x);
@@ -306,13 +298,10 @@ namespace AltMath
     return  x < 0.0 && n%2 == 1 ? -ans : ans;
   }
   
-#undef ACC
-#undef BIGNO
-#undef BIGNI
-
-#define TOL 1.0e-7
   double jn_sph0(double x)
   {
+    const double TOL = 1.0e-7;
+
     if (fabs(x)<TOL)
       return 1.0-x*x/6.0;
     else
@@ -321,6 +310,8 @@ namespace AltMath
 
   double jn_sph1(double x)
   {
+    const double TOL = 1.0e-7;
+
     if (fabs(x)<TOL)
       return x/3.0;
     else
@@ -331,7 +322,6 @@ namespace AltMath
   {
     return cos(x)/x;
   }
-#undef TOL
 
   double check_nu(double nu)
   {
diff --git a/exputil/EmpCyl2d.cc b/exputil/EmpCyl2d.cc
index efd29bcdf..37f36cfb0 100644
--- a/exputil/EmpCyl2d.cc
+++ b/exputil/EmpCyl2d.cc
@@ -780,19 +780,31 @@ bool EmpCyl2d::ReadH5Cache()
     auto checkInt = [&file](int value, std::string name)
     {
       int v; HighFive::Attribute vv = file.getAttribute(name); vv.read(v);
-      if (value == v) return true; return false;
+      if (value == v) return true;
+      if (myid==0) std::cout << "---- EmpCyl2d::ReadH5Cache: "
+			     << name << " expected " << value << " but found "
+			     << v << std::endl;
+      return false;
     };
 
     auto checkDbl = [&file](double value, std::string name)
     {
       double v; HighFive::Attribute vv = file.getAttribute(name); vv.read(v);
-      if (fabs(value - v) < 1.0e-16) return true; return false;
+      if (fabs(value - v) < 1.0e-16) return true;
+      if (myid==0) std::cout << "---- EmpCyl2d::ReadH5Cache: "
+			     << name << " expected " << value << " but found "
+			     << v << std::endl;
+      return false;
     };
 
     auto checkStr = [&file](std::string value, std::string name)
     {
       std::string v; HighFive::Attribute vv = file.getAttribute(name); vv.read(v);
-      if (value.compare(v)==0) return true; return false;
+      if (value.compare(v)==0) return true;
+      if (myid==0) std::cout << "---- EmpCyl2d::ReadH5Cache: "
+			     << name << " expected " << value << " but found "
+			     << v << std::endl;
+      return false;
     };
 
     //
@@ -803,8 +815,8 @@ bool EmpCyl2d::ReadH5Cache()
       if (not checkStr(Version, "Version"))  return false;
     } else {
       if (myid==0)
-	std::cout << "---- EmpCyl2d::ReadH5Cache: "
-		  << "recomputing cache for HighFive API change"
+	std::cout << "---- EmpCyl2d::ReadH5Cache: " << "<" << cache_name_2d
+		  << "> recomputing cache for HighFive API change"
 		  << std::endl;
       return false;
     }
diff --git a/exputil/massmodel_dist.cc b/exputil/massmodel_dist.cc
index 27e2d877f..618a4bd8e 100644
--- a/exputil/massmodel_dist.cc
+++ b/exputil/massmodel_dist.cc
@@ -49,14 +49,14 @@
 #include <massmodel.H>
 #include <interp.H>
 
-#define OFFSET 1.0e-3
-#define OFFTOL 1.2
+const double OFFSET=1.0e-3;
+const double OFFTOL=1.2;
 
 extern double gint_0(double a, double b, std::function<double(double)> f, int NGauss);
 extern double gint_2(double a, double b, std::function<double(double)> f, int NGauss);
 
-#define TSTEP 1.0e-8
-#define NGauss 96
+const double TSTEP=1.0e-8;
+const int    NGauss=96;
 
 static int DIVERGE=0;
 
@@ -109,6 +109,10 @@ void SphericalModelTable::setup_df(int NUM, double RA)
   rhoQy.resize(num);
   rhoQy2.resize(num);
 
+  rhoQx. setZero();
+  rhoQy. setZero();
+  rhoQy2.setZero();
+
   for (int i=0; i<num; i++) {
     x = density.x[i];
     rhoQx[i] = pot.y[i];
@@ -138,6 +142,11 @@ void SphericalModelTable::setup_df(int NUM, double RA)
     dfc.Q  .resize(NUM);
     dfc.fQ .resize(NUM);
     dfc.ffQ.resize(NUM);
+
+    dfc.Q  .setZero();
+    dfc.fQ .setZero();
+    dfc.ffQ.setZero();
+
     dfc.num = NUM;
     dfc.ra2 = ra2;
 
@@ -145,7 +154,8 @@ void SphericalModelTable::setup_df(int NUM, double RA)
     Qmin = get_pot(pot.x[0]);
     dQ = (Qmax - Qmin)/(double)(dfc.num-1);
   
-    foffset = -std::numeric_limits<double>::max();
+    // foffset = -std::numeric_limits<double>::max();
+    foffset = -1.0e42;
     dfc.Q[dfc.num-1] = Qmax;
     dfc.ffQ[dfc.num-1] = 0.0;
     fac = 1.0/(sqrt(8.0)*M_PI*M_PI);
@@ -184,6 +194,13 @@ void SphericalModelTable::setup_df(int NUM, double RA)
     df.ffQ .resize(NUM);
     df.fQ2 .resize(NUM);
     df.ffQ2.resize(NUM);
+
+    df.Q   .setZero();
+    df.fQ  .setZero();
+    df.ffQ .setZero();
+    df.fQ2 .setZero();
+    df.ffQ2.setZero();
+
     df.num = NUM;
     df.ra2 = ra2;
 
@@ -194,7 +211,8 @@ void SphericalModelTable::setup_df(int NUM, double RA)
     df.Q[df.num-1] = Qmax;
     df.ffQ[df.num-1] = 0.0;
     fac = 1.0/(sqrt(8.0)*M_PI*M_PI);
-    foffset = -std::numeric_limits<double>::max();
+    // foffset = -std::numeric_limits<double>::max();
+    foffset = -1.0e42;
     for (int i=df.num-2; i>=0; i--) {
       df.Q[i] = df.Q[i+1] - dQ;
       Q = df.Q[i];
@@ -233,7 +251,7 @@ void SphericalModelTable::setup_df(int NUM, double RA)
 
   dist_defined = true;
 
-  debug_fdist();
+  // debug_fdist();
 }
 
 
@@ -294,7 +312,7 @@ double SphericalModelTable::distf(double E, double L)
 
   if (!dist_defined) bomb("distribution function not defined");
 
-  double d, g;
+  double d=0.0, g=0.0;
 
   if (chebyN) {
 
@@ -341,7 +359,7 @@ double SphericalModelTable::dfde(double E, double L)
 
   if (!dist_defined) bomb("distribution function not defined");
 
-  double d, g, h, d1;
+  double d=0, g=0, h=0, d1=0;
 
   if (chebyN) {
     
@@ -399,7 +417,7 @@ double SphericalModelTable::dfdl(double E, double L)
 
   if (!dist_defined) bomb("distribution function not defined");
 
-  double d, g, h, d1;
+  double d=0, g=0, h=0, d1=0;
 
   if (chebyN) {
     
@@ -451,7 +469,7 @@ double SphericalModelTable::d2fde2(double E, double L)
 {
   if (!dist_defined) bomb("distribution function not defined");
 
-  double d, g, h, k, d2;
+  double d=0, g=0, h=0, k=0, d2=0;
 
   if (chebyN) {
 
diff --git a/exputil/realize_model.cc b/exputil/realize_model.cc
index 0adfc7b0a..14ffc230c 100644
--- a/exputil/realize_model.cc
+++ b/exputil/realize_model.cc
@@ -37,6 +37,7 @@
 #include <localmpi.H>
 #include <massmodel.H>
 #include <interp.H>
+#include <euler.H>
 
 #ifdef DEBUG
 #include <orbit.H>
@@ -69,12 +70,10 @@ int      AxiSymModel::gen_itmax = 20000;
 const bool verbose = true;
 const double ftol = 0.01;
 
-Eigen::VectorXd AxiSymModel::gen_point_2d(int& ierr)
+AxiSymModel::PSret AxiSymModel::gen_point_2d()
 {
   if (!dist_defined) {
-    std::cerr << "AxiSymModel: must define distribution before realizing!"
-	      << std::endl;
-    exit (-1);
+    throw std::runtime_error("AxiSymModel: must define distribution before realizing!");
   }
 
   int it;			// Iteration counter
@@ -87,7 +86,6 @@ Eigen::VectorXd AxiSymModel::gen_point_2d(int& ierr)
   double Emin = get_pot(rmin);
   double Emax = get_pot(get_max_radius());
 
-
   if (gen_EJ) {
 
     if (gen_firstime) {
@@ -97,8 +95,9 @@ Eigen::VectorXd AxiSymModel::gen_point_2d(int& ierr)
       double dx = (Emax - Emin - 2.0*tol)/(numr-1);
       double dy = (1.0 - gen_kmin - 2.0*tol)/(numj-1);
 
+#ifdef DEBUG
       std::cout << "gen_point_2d[" << ModelID << "]: " << get_max_radius() << std::endl;
-      
+#endif
       gen_fomax = 0.0;
 
       for (int j=0; j<numr; j++) {
@@ -109,9 +108,9 @@ Eigen::VectorXd AxiSymModel::gen_point_2d(int& ierr)
 
 	  gen_orb.new_orbit(xxx, yyy);
 	    
-	  double zzz = distf(xxx, gen_orb.AngMom())*gen_orb.Jmax()
-	    /gen_orb.get_freq(1);
-	  gen_fomax = zzz>gen_fomax ? zzz : gen_fomax;
+	  double zzz = distf(xxx, gen_orb.AngMom())*gen_orb.Jmax()/gen_orb.get_freq(1);
+
+	  if (zzz>gen_fomax) gen_fomax = zzz;
 	}
       }
 
@@ -119,8 +118,8 @@ Eigen::VectorXd AxiSymModel::gen_point_2d(int& ierr)
     }
 
     
-				// Trial velocity point
-    
+    // Trial velocity point
+    //
     for (it=0; it<gen_itmax; it++) {
 
       double xxx = Emin + tol + (Emax - Emin - 2.0*tol)*Unit(random_gen);
@@ -166,7 +165,9 @@ Eigen::VectorXd AxiSymModel::gen_point_2d(int& ierr)
       gen_rloc.resize(gen_N);
       gen_fmax.resize(gen_N);
 
+#ifdef DEBUG
       std::cout << "gen_point_2d[" << ModelID << "]: " << get_max_radius() << std::endl;
+#endif
 
       if (rmin <= 1.0e-16) gen_logr = 0;
 
@@ -203,7 +204,7 @@ Eigen::VectorXd AxiSymModel::gen_point_2d(int& ierr)
 	    eee = pot + 0.5*(vr*vr + vt*vt);
 
 	    double zzz = distf(eee, gen_rloc[i]*vt);
-	    fmax = zzz>fmax ? zzz : fmax;
+	    if (zzz > fmax) fmax = zzz;
 	  }
 	}
 	gen_fmax[i] = fmax*(1.0 + ftol);
@@ -219,8 +220,8 @@ Eigen::VectorXd AxiSymModel::gen_point_2d(int& ierr)
     pot = get_pot(r);
     vmax = sqrt(2.0*fabs(Emax - pot));
 
-                                // Trial velocity point
-    
+    // Trial velocity point
+    //
     for (it=0; it<gen_itmax; it++) {
 
       double xxx = sqrt(Unit(random_gen));
@@ -232,7 +233,7 @@ Eigen::VectorXd AxiSymModel::gen_point_2d(int& ierr)
 
       if (Unit(random_gen) > distf(eee, r*vt)/fmax ) continue;
 
-      if (Unit(random_gen)<0.5) vr *= -1.0;
+      if (Unit(random_gen) < 0.5) vr *= -1.0;
     
       phi = 2.0*M_PI*Unit(random_gen);
 
@@ -254,13 +255,10 @@ Eigen::VectorXd AxiSymModel::gen_point_2d(int& ierr)
 		<< std::setw(12) << fmax
 		<< " %=" << std::setw(12)
 		<< static_cast<double>(++toomany)/totcnt << std::endl;
-    ierr = 1;
     out.setZero();
-    return out;
+    return {out, 1};
   }
 
-  ierr = 0;
-  
   cosp = cos(phi);
   sinp = sin(phi);
 
@@ -272,16 +270,14 @@ Eigen::VectorXd AxiSymModel::gen_point_2d(int& ierr)
   out[5] = vr * sinp + vt * cosp;
   out[6] = 0.0;
     
-  return out;
+  return {out, 0};
 }
 
 
-Eigen::VectorXd AxiSymModel::gen_point_2d(double r, int& ierr)
+AxiSymModel::PSret AxiSymModel::gen_point_2d(double r)
 {
   if (!dist_defined) {
-    std::cerr << "AxiSymModel: must define distribution before realizing!"
-	      << std::endl;
-    exit (-1);
+    throw std::runtime_error("AxiSymModel: must define distribution before realizing!");
   }
 
   int it;			// Iteration counter
@@ -300,8 +296,10 @@ Eigen::VectorXd AxiSymModel::gen_point_2d(double r, int& ierr)
     gen_rloc.resize(gen_N);
     gen_fmax.resize(gen_N);
 
+#ifdef DEBUG
     std::cout << "gen_point_2d[" << ModelID << "]: " << rmin
 	 << ", " << get_max_radius() << std::endl;
+#endif
 
     double dr = (get_max_radius() - rmin)/gen_N;
 
@@ -324,7 +322,7 @@ Eigen::VectorXd AxiSymModel::gen_point_2d(double r, int& ierr)
 	  eee = pot + 0.5*(vr*vr + vt*vt);
 
 	  double zzz = distf(eee, gen_rloc[i]*vt);
-	  fmax = zzz>fmax ? zzz : fmax;
+	  if (zzz>fmax) fmax = zzz;
 	}
       }
       gen_fmax[i] = fmax*(1.0 + ftol);
@@ -337,8 +335,8 @@ Eigen::VectorXd AxiSymModel::gen_point_2d(double r, int& ierr)
   pot = get_pot(r);
   vmax = sqrt(2.0*fabs(Emax - pot));
   
-                                // Trial velocity point
-    
+  // Trial velocity point
+  //
   for (it=0; it<gen_itmax; it++) {
 
     double xxx = sqrt(Unit(random_gen));
@@ -350,7 +348,7 @@ Eigen::VectorXd AxiSymModel::gen_point_2d(double r, int& ierr)
 
     if (Unit(random_gen) > distf(eee, r*vt)/fmax ) continue;
     
-    if (Unit(random_gen)<0.5) vr *= -1.0;
+    if (Unit(random_gen) < 0.5) vr *= -1.0;
     
     phi = 2.0*M_PI*Unit(random_gen);
 
@@ -370,13 +368,10 @@ Eigen::VectorXd AxiSymModel::gen_point_2d(double r, int& ierr)
 		<< std::setw(12) << fmax
 		<< " %=" << std::setw(12)
 		<< static_cast<double>(++toomany)/totcnt << std::endl;
-    ierr = 1;
     out.setZero();
-    return out;
+    return {out, 1};
   }
 
-  ierr = 0;
-  
   cosp = cos(phi);
   sinp = sin(phi);
 
@@ -388,16 +383,14 @@ Eigen::VectorXd AxiSymModel::gen_point_2d(double r, int& ierr)
   out[5] = vr * sinp + vt * cosp;
   out[6] = 0.0;
     
-  return out;
+  return {out, 0};
 }
 
 
-Eigen::VectorXd AxiSymModel::gen_point_3d(int& ierr)
+AxiSymModel::PSret AxiSymModel::gen_point_3d()
 {
   if (!dist_defined) {
-    std::cerr << "AxiSymModel: must define distribution before realizing!"
-	      << std::endl;
-    exit (-1);
+    throw std::runtime_error("AxiSymModel: must define distribution before realizing!");
   }
 
 #ifdef DEBUG
@@ -405,7 +398,7 @@ Eigen::VectorXd AxiSymModel::gen_point_3d(int& ierr)
 #endif
 
   double r, pot, vmax, vr=0.0, vt, eee, vt1=0.0, vt2=0.0, fmax;
-  double phi, sint, cost, sinp, cosp, azi;
+  double phi, sint, cost, sinp, cosp;
 
   double rmin = max<double>(get_min_radius(), gen_rmin);
   double Emax = get_pot(get_max_radius());
@@ -461,7 +454,7 @@ Eigen::VectorXd AxiSymModel::gen_point_3d(int& ierr)
 	  eee = pot + 0.5*(vr*vr + vt*vt);
 
 	  double zzz = distf(eee, r*vt);
-	  fmax = zzz>fmax ? zzz : fmax;
+	  if (zzz>fmax) fmax = zzz;
 	}
       }
       gen_fmax[i] = fmax*(1.0 + ftol);
@@ -525,7 +518,9 @@ Eigen::VectorXd AxiSymModel::gen_point_3d(int& ierr)
 
     if (Unit(random_gen)<0.5) vr *= -1.0;
     
-    azi = 2.0*M_PI*Unit(random_gen);
+    // Orientation of tangential velocity vector
+    //
+    double azi = 2.0*M_PI*Unit(random_gen);
     vt1 = vt*cos(azi);
     vt2 = vt*sin(azi);
 
@@ -544,13 +539,10 @@ Eigen::VectorXd AxiSymModel::gen_point_3d(int& ierr)
 		<< std::setw(12) << fmax
 		<< " %=" << std::setw(12)
 		<< static_cast<double>(++toomany)/totcnt << std::endl;
-    ierr = 1;
     out.setZero();
-    return out;
+    return {out, 1};
   }
 
-  ierr = 0;
-  
   if (Unit(random_gen)>=0.5) vr *= -1.0;
 
   phi = 2.0*M_PI*Unit(random_gen);
@@ -583,17 +575,284 @@ Eigen::VectorXd AxiSymModel::gen_point_3d(int& ierr)
 #endif
 
 
-  return out;
+  return {out, 0};
 }
 
 
-Eigen::VectorXd AxiSymModel::gen_point_3d(double Emin, double Emax, 
-					  double Kmin, double Kmax, int& ierr)
+AxiSymModel::PSret AxiSymModel::gen_point_3d(double r, double theta, double phi)
 {
   if (!dist_defined) {
-    std::cerr << "AxiSymModel: must define distribution before realizing!"
-	      << std::endl;
-    exit (-1);
+    throw std::runtime_error("AxiSymModel: must define distribution before realizing!");
+  }
+
+  int it;			// Iteration counter
+  double pot, vmax, vr=0.0, vt=0.0, eee, fmax;
+  double rmin = max<double>(get_min_radius(), gen_rmin);
+  double Emax = get_pot(get_max_radius());
+
+  if (gen_firstime) {
+
+    double tol = 1.0e-5;
+    double dx = (1.0 - 2.0*tol)/(numr-1);
+    double dy = (1.0 - 2.0*tol)/(numj-1);
+
+    gen_mass.resize(gen_N);
+    gen_rloc.resize(gen_N);
+    gen_fmax.resize(gen_N);
+
+#ifdef DEBUG
+    std::cout << "gen_point_3d[" << ModelID << "]: " << rmin
+	 << ", " << get_max_radius() << std::endl;
+#endif
+
+    double dr = (get_max_radius() - rmin)/gen_N;
+
+    for (int i=0; i<gen_N; i++) {
+      gen_rloc[i] = rmin + dr*(0.5+i);
+      gen_mass[i] = get_mass(gen_rloc[i]);
+
+      pot = get_pot(gen_rloc[i]);
+      vmax = sqrt(2.0*fabs(Emax - pot));
+
+      fmax = 0.0;
+      for (int j=0; j<numr; j++) {
+	double xxx = sqrt(tol + dx*j);
+
+	for (int k=0; k<numj; k++) {
+	  double yyy = 0.5*M_PI*(tol + dy*k);
+
+	  vr = vmax*xxx*cos(yyy);
+	  vt = vmax*xxx*sin(yyy);
+	  eee = pot + 0.5*(vr*vr + vt*vt);
+
+	  double zzz = distf(eee, gen_rloc[i]*vt);
+	  if (zzz>fmax) fmax = zzz;
+	}
+      }
+      gen_fmax[i] = fmax*(1.0 + ftol);
+
+    }
+    gen_firstime = false;
+
+    std::ofstream out("gen_point_3d.dat");
+    if (out) {
+      out << "# gen_point_3d[" << ModelID << "]" << std::endl
+	  << "# " << std::setw(14) << "r" << std::setw(16) << "mass"
+	  << std::endl;
+      for (int i=0; i<gen_N; i++) {
+	out << std::setw(16) << gen_rloc[i]
+	    << std::setw(16) << gen_mass[i]
+	    << std::endl;
+      }
+    }
+  }
+
+  fmax = odd2(r, gen_rloc, gen_fmax, 1);
+  pot  = get_pot(r);
+  vmax = sqrt(2.0*fabs(Emax - pot));
+  
+  // Trial velocity point
+  //
+  for (it=0; it<gen_itmax; it++) {
+
+    double xxx = pow(Unit(random_gen), 0.3333333333333333);
+    double cosq = 2.0*Unit(random_gen)-1.0;
+    double sinq = sqrt(fabs(1.0 - cosq*cosq));
+
+    vr = vmax*xxx*cosq;
+    vt = vmax*xxx*sinq;
+    eee = pot + 0.5*(vr*vr + vt*vt);
+
+    if (Unit(random_gen) > distf(eee, r*vt)/fmax ) continue;
+    
+    if (Unit(random_gen)<0.5) vr *= -1.0;
+    
+    break;
+  }
+
+  Eigen::VectorXd out(7);
+
+  if (std::isnan(vr) or std::isnan(vt)) {
+    if (verbose) std::cout << "NaN found in AxiSymModel::gen_point_3d with r="
+			   << r << " theta=" << theta << " phi=" << phi
+			   << std::endl;
+    out.setZero();
+    return {out, 1};
+  }
+  
+  static unsigned totcnt = 0, toomany = 0;
+  totcnt++;
+
+  if (it==gen_itmax) {
+    if (verbose)
+      std::cerr << "Velocity selection failed [" << myid << "]: r="
+		<< std::setw(12) << r
+		<< std::setw(12) << fmax
+		<< " %=" << std::setw(12)
+		<< static_cast<double>(++toomany)/totcnt << std::endl;
+    out.setZero();
+    return {out, 1};
+  }
+
+  double tv  = 2.0*M_PI*Unit(random_gen);
+  double vth = vt*cos(tv);
+  double vph = vt*sin(tv);
+
+  double cost = cos(theta);
+  double sint = sin(theta);
+
+  double cosp = cos(phi);
+  double sinp = sin(phi);
+
+  out[0] = 1.0;
+  out[1] = r*sint*cosp;
+  out[2] = r*sint*sinp;
+  out[3] = r*cost;
+  out[4] = vr*sint*cosp - vph*sinp + vth*cost*cosp;
+  out[5] = vr*sint*sinp + vph*cosp + vth*cost*sinp;
+  out[6] = vr*cost - vth*sint;
+    
+  return {out, 0};
+}
+
+AxiSymModel::PSret AxiSymModel::gen_point_3d_iso
+(double r, double theta, double phi, int N, int nv, int na,
+ double PHI, double THETA, double PSI)
+{
+  if (!dist_defined) {
+    throw std::runtime_error("AxiSymModel: must define distribution before realizing!");
+  }
+
+  double Emax = get_pot(get_max_radius());
+  double vpot = get_pot(r);
+  double v3   = pow(2.0*fabs(Emax - vpot), 1.5);
+  int knots   = 20;
+
+  // Sanity check
+  assert((N < nv*na));
+
+  // Rotation matrix
+  static Eigen::Matrix3d rotate;
+
+  // Assume isotropic and generate new r slice
+  if (fabs(gen_lastr-r) > 1.0e-14) {
+    gen_mass.resize(gen_N+1);
+    gen_rloc.resize(gen_N+1);
+    gen_lastr = r;
+
+#ifdef DEBUG
+    std::cout << "gen_point_3d_iso[" << ModelID << "]: " << rmin
+	 << ", " << get_max_radius() << std::endl;
+#endif
+
+    LegeQuad lv(knots);
+
+    double dv3 = v3/gen_N;
+
+    gen_rloc[0] = 0.0;
+    gen_mass[0] = 0.0;
+    for (int i=0; i<gen_N; i++) {
+      double vvv = dv3*(0.5+i);
+      gen_rloc[i+1] = vvv;
+
+      double sum = 0.0;
+      for (int l=0; l<knots; l++) {
+	double v = pow(vvv + dv3*(lv.knot(l)-0.5), 0.3333333333333333);
+	double eee = vpot + 0.5*v*v;
+	sum += dv3*lv.weight(l)*distf(eee, r*v);
+      }
+      gen_mass[i+1] = sum;
+    }
+
+    // Cumulate mass
+    for (int i=1; i<=gen_N; i++) gen_mass[i] += gen_mass[i-1];
+
+    std::ofstream out("gen_point_3d_iso.dat");
+    if (out) {
+      out << "# gen_point_3d_iso[" << ModelID << "]" << std::endl
+	  << "# " << std::setw(14) << "r" << std::setw(16) << "mass"
+	  << std::endl;
+      for (int i=0; i<gen_N; i++) {
+	out << std::setw(16) << gen_rloc[i]
+	    << std::setw(16) << gen_mass[i]
+	    << std::endl;
+	if (gen_rloc[i] <0.0 or gen_mass[i] < 0.0) {
+	  std::cout << "ERROR: gen_point_3d_iso[" << ModelID << "]: "
+		    << "r=" << gen_rloc[i] << " mass=" << gen_mass[i]
+		    << std::endl;
+	}
+      }
+    }
+
+    rotate = return_euler(PHI, THETA, PSI, 1);
+  }
+
+  double dmas = gen_mass[gen_N]/nv;
+
+  int jv = N/na;
+  int ja = N - jv*na;
+
+  // Sanity checks
+  //
+  assert((jv>=0 and jv<nv));
+  assert((ja>=0 and ja<na));
+
+  // Fibonacci lattice for velocities
+  //
+  const double goldenRatio = 1.618033988749895;
+
+  double vphi  = 2.0*M_PI * ja / goldenRatio;
+  double vcost = 1.0 - 2.0*ja/na;
+  double vsint = sqrt(fabs(1.0 - vcost*vcost));
+  
+  double vvv = odd2(dmas*(jv+0.5), gen_mass, gen_rloc);
+  double v   = pow(vvv, 0.3333333333333333);
+
+  Eigen::VectorXd out(7);
+
+  if (std::isnan(v)) {
+    if (verbose) std::cout << "NaN found in AxiSymModel::gen_point_3d_iso with r="
+			   << r << " theta=" << theta << " phi=" << phi
+			   << " v3=" << vvv  << " jv=" << jv << " m="
+			   << dmas*(jv+0.5) << " [" << gen_mass[0]
+			   << ", " <<  gen_mass[gen_N] << "]" << std::endl;
+    out.setZero();
+    return {out, 1};
+  }
+  
+  double vr  = v*vcost;
+  double vph = v*vsint*cos(vphi);
+  double vth = v*vsint*sin(vphi);
+
+  Eigen::Vector3d V = rotate * Eigen::Vector3d(vr, vphi, vth);
+
+  vr  = V[0];
+  vph = V[1];
+  vth = V[2];
+
+  double cost = cos(theta);
+  double sint = sin(theta);
+
+  double cosp = cos(phi);
+  double sinp = sin(phi);
+
+  out[0] = 1.0;
+  out[1] = r*sint*cosp;
+  out[2] = r*sint*sinp;
+  out[3] = r*cost;
+  out[4] = vr*sint*cosp - vph*sinp + vth*cost*cosp;
+  out[5] = vr*sint*sinp + vph*cosp + vth*cost*sinp;
+  out[6] = vr*cost - vth*sint;
+    
+  return {out, 0};
+}
+
+
+AxiSymModel::PSret AxiSymModel::gen_point_3d(double Emin, double Emax, 
+				double Kmin, double Kmax)
+{
+  if (!dist_defined) {
+    throw std::runtime_error("AxiSymModel: must define distribution before realizing!");
   }
 
 #ifdef DEBUG
@@ -602,7 +861,7 @@ Eigen::VectorXd AxiSymModel::gen_point_3d(double Emin, double Emax,
 #endif
 
   double r, vr, vt, vt1, vt2, E, K, J, jmax, w1t, eee, pot;
-  double phi, sint, cost, sinp, cosp, azi;
+  double phi, sint, cost, sinp, cosp;
   double rmin = max<double>(get_min_radius(), gen_rmin);
   double rmax = get_max_radius();
 
@@ -688,7 +947,8 @@ Eigen::VectorXd AxiSymModel::gen_point_3d(double Emin, double Emax,
 #endif    
   }
 
-				// Enforce limits
+  // Enforce limits
+  //
   Emin = max<double>(Emin, Emin_grid);
   Emin = min<double>(Emin, Emax_grid);
   Emax = max<double>(Emax, Emin_grid);
@@ -733,9 +993,10 @@ Eigen::VectorXd AxiSymModel::gen_point_3d(double Emin, double Emax,
       
   pot = get_pot(r);
   vt = J/r;
-  ierr = 0;
-				// Interpolation check (should be rare)
-				// Error condition is set
+  int ierr = 0;
+
+  // Interpolation check (should be rare).  Error condition is set.
+  //
   if (2.0*(E - pot) - vt*vt < 0.0) {
     ierr = 1;
     if (E < pot) E = pot;
@@ -745,7 +1006,9 @@ Eigen::VectorXd AxiSymModel::gen_point_3d(double Emin, double Emax,
 
   if (Unit(random_gen)<0.5) vr *= -1.0;
     
-  azi = 2.0*M_PI*Unit(random_gen);
+  // Orientation of tangential velocity vector
+  //
+  double azi = 2.0*M_PI*Unit(random_gen);
   vt1 = vt*cos(azi);
   vt2 = vt*sin(azi);
 
@@ -767,7 +1030,7 @@ Eigen::VectorXd AxiSymModel::gen_point_3d(double Emin, double Emax,
     
 
 #ifdef DEBUG
-  if (ierr) return out;
+  if (ierr) return {out, ierr};
 
   eee = pot + 0.5*(out[4]*out[4]+out[5]*out[5]+out[6]*out[6]);
   orb.new_orbit(eee, 0.5);
@@ -787,15 +1050,15 @@ Eigen::VectorXd AxiSymModel::gen_point_3d(double Emin, double Emax,
        << std::endl;
 #endif
 
-  return out;
+  return {out, ierr};
 }
 
 
 
-Eigen::VectorXd AxiSymModel::gen_point_jeans_3d(int& ierr)
+AxiSymModel::PSret AxiSymModel::gen_point_jeans_3d()
 {
   double r, d, vr, vt, vt1, vt2, vv, vtot;
-  double phi, sint, cost, sinp, cosp, azi;
+  double phi, sint, cost, sinp, cosp;
   double rmin = max<double>(get_min_radius(), gen_rmin);
 
   if (gen_firstime_jeans) {
@@ -838,8 +1101,8 @@ Eigen::VectorXd AxiSymModel::gen_point_jeans_3d(int& ierr)
 
 
     // Debug
-    
-    ofstream test("test.grid");
+    //
+    std::ofstream test("test.grid");
     if (test) {
 
       test << "# [Jeans] Rmin=" << rmin
@@ -885,7 +1148,9 @@ Eigen::VectorXd AxiSymModel::gen_point_jeans_3d(int& ierr)
   vr = vtot*xxx;
   vt = vtot*sqrt(yyy);
 
-  azi = 2.0*M_PI*Unit(random_gen);
+  // Orientation of tangential velocity vector
+  //
+  double azi = 2.0*M_PI*Unit(random_gen);
   vt1 = vt*cos(azi);
   vt2 = vt*sin(azi);
 
@@ -907,24 +1172,20 @@ Eigen::VectorXd AxiSymModel::gen_point_jeans_3d(int& ierr)
   out[5] = vr * sint*sinp + vt1 * cost*sinp + vt2*cosp;
   out[6] = vr * cost      - vt1 * sint;
     
-  ierr = 0;
-
-  return out;
+  return {out, 0};
 }
 
-void AxiSymModel::gen_velocity(double* pos, double* vel, int& ierr)
+int AxiSymModel::gen_velocity(double* pos, double* vel)
 {
   if (dof()!=3)
       bomb( "AxiSymModel: gen_velocity only implemented for 3d model!" );
-
+  
   if (!dist_defined) {
-    std::cerr << "AxiSymModel: must define distribution before realizing!"
-	      << std::endl;
-    exit (-1);
+    throw std::runtime_error("AxiSymModel: must define distribution before realizing!");
   }
 
   double r, pot, vmax, vr=0.0, vt, eee, vt1=0.0, vt2=0.0, fmax;
-  double phi, sint, cost, sinp, cosp, azi;
+  double phi, sint, cost, sinp, cosp;
   double rmin = max<double>(get_min_radius(), gen_rmin);
 
   double Emax = get_pot(get_max_radius());
@@ -976,7 +1237,7 @@ void AxiSymModel::gen_velocity(double* pos, double* vel, int& ierr)
 	  eee = pot + 0.5*(vr*vr + vt*vt);
 
 	  double zzz = distf(eee, r*vt);
-	  fmax = zzz>fmax ? zzz : fmax;
+	  if (zzz>fmax) fmax = zzz;
 	}
       }
       gen_fmax[i] = fmax*(1.0 + ftol);
@@ -1021,7 +1282,9 @@ void AxiSymModel::gen_velocity(double* pos, double* vel, int& ierr)
 
     if (Unit(random_gen)<0.5) vr *= -1.0;
     
-    azi = 2.0*M_PI*Unit(random_gen);
+    // Orientation of tangential velocity vector
+    //
+    double azi = 2.0*M_PI*Unit(random_gen);
     vt1 = vt*cos(azi);
     vt2 = vt*sin(azi);
 
@@ -1038,12 +1301,9 @@ void AxiSymModel::gen_velocity(double* pos, double* vel, int& ierr)
 		<< std::setw(12) << fmax
 		<< " %=" << std::setw(12)
 		<< static_cast<double>(++toomany)/totcnt << std::endl;
-    ierr = 1;
-    return;
+    return 1;
   }
 
-  ierr = 0;
-  
   if (Unit(random_gen)>=0.5) vr *= -1.0;
 
   phi = atan2(pos[1], pos[0]);
@@ -1055,18 +1315,19 @@ void AxiSymModel::gen_velocity(double* pos, double* vel, int& ierr)
   vel[0] = vr * sint*cosp + vt1 * cost*cosp - vt2*sinp;
   vel[1] = vr * sint*sinp + vt1 * cost*sinp + vt2*cosp;
   vel[2] = vr * cost      - vt1 * sint;
+
+  return 0;
 }
 
-Eigen::VectorXd SphericalModelMulti::gen_point(int& ierr)
+AxiSymModel::PSret SphericalModelMulti::gen_point()
 {
   if (!real->dist_defined || !fake->dist_defined) {
-    std::cerr << "SphericalModelMulti: input distribution functions must be defined before realizing!" << std::endl;
-    exit (-1);
+    throw std::runtime_error("SphericalModelMulti: input distribution functions must be defined before realizing!");
   }
 
   double r, pot, vmax;
   double vr=0.0, vt=0.0, eee=0.0, vt1=0.0, vt2=0.0, fmax, emax;
-  double mass, phi, sint, cost, sinp, cosp, azi;
+  double mass, phi, sint, cost, sinp, cosp;
 
   double Emax = get_pot(get_max_radius());
 
@@ -1271,6 +1532,7 @@ Eigen::VectorXd SphericalModelMulti::gen_point(int& ierr)
   for (it=0; it<gen_itmax; it++) {
 #else
   // Trial phase-space point
+  //
   for (it=0; it<gen_itmax; it++) {
 
     // Generate a radius (inside the loop)
@@ -1308,7 +1570,9 @@ Eigen::VectorXd SphericalModelMulti::gen_point(int& ierr)
 
     if (Unit(random_gen) < 0.5) vr *= -1.0;
     
-    azi = 2.0*M_PI*Unit(random_gen);
+    // Orientation of tangential velocity vector
+    //
+    double azi = 2.0*M_PI*Unit(random_gen);
     vt1 = vt*cos(azi);
     vt2 = vt*sin(azi);
 
@@ -1320,7 +1584,6 @@ Eigen::VectorXd SphericalModelMulti::gen_point(int& ierr)
   static unsigned totcnt = 0, toomany = 0;
   totcnt++;
 
-
   if (it==gen_itmax) {
     if (verbose) {
       std::cerr << "Velocity selection failed [" << myid << "]: r="
@@ -1331,13 +1594,10 @@ Eigen::VectorXd SphericalModelMulti::gen_point(int& ierr)
 		<< " %=" << std::setw(12)
 		<< static_cast<double>(++toomany)/totcnt << std::endl;
     }
-    ierr = 1;
     out.setZero();
-    return out;
+    return {out, 1};
   }
 
-  ierr = 0;
-  
   if (Unit(random_gen)>=0.5) vr *= -1.0;
 
   phi  = 2.0*M_PI*Unit(random_gen);
@@ -1375,24 +1635,23 @@ Eigen::VectorXd SphericalModelMulti::gen_point(int& ierr)
   out[5] = vr * sint*sinp + vt1 * cost*sinp + vt2*cosp;
   out[6] = vr * cost      - vt1 * sint;
     
+  int ierr = 0;
   for (int i=0; i<7; i++) {
     if (std::isnan(out[i]) || std::isinf(out[i])) ierr = 1;
   }
 
-  return out;
+  return {out, ierr};
 }
 
 
-Eigen::VectorXd SphericalModelMulti::gen_point(double radius, int& ierr)
+AxiSymModel::PSret SphericalModelMulti::gen_point(double radius)
 {
   if (!real->dist_defined || !fake->dist_defined) {
-    std::cerr << "SphericalModelMulti: input distribution functions must be defined before realizing!" << std::endl;
-    exit (-1);
+    throw std::runtime_error("SphericalModelMulti: input distribution functions must be defined before realizing!");
   }
 
   double r, pot, vmax;
   double vr=0.0, vt=0.0, eee=0.0, vt1=0.0, vt2=0.0, fmax;
-  double phi, sint, cost, sinp, cosp, azi;
 
   double Emax = get_pot(get_max_radius());
 
@@ -1441,7 +1700,8 @@ Eigen::VectorXd SphericalModelMulti::gen_point(double radius, int& ierr)
 	  eee = pot + 0.5*(vr*vr + vt*vt);
 
 	  double zzz = fake->distf(eee, r*vt);
-	  fmax = zzz>fmax ? zzz : fmax;
+
+	  if (zzz>fmax) fmax = zzz;
 	}
       }
       gen_fmax[i] = fmax*(1.0 + ftol);
@@ -1500,7 +1760,9 @@ Eigen::VectorXd SphericalModelMulti::gen_point(double radius, int& ierr)
 
     if (Unit(random_gen)<0.5) vr *= -1.0;
     
-    azi = 2.0*M_PI*Unit(random_gen);
+    // Orientation of the tangential velocity vector
+    //
+    double azi = 2.0*M_PI*Unit(random_gen);
     vt1 = vt*cos(azi);
     vt2 = vt*sin(azi);
 
@@ -1519,20 +1781,17 @@ Eigen::VectorXd SphericalModelMulti::gen_point(double radius, int& ierr)
 		<< " reject="  << std::setw( 6) << reject
 		<< std::endl;
     }
-    ierr = 1;
     out.setZero();
-    return out;
+    return {out, 1};
   }
 
-  ierr = 0;
-  
   if (Unit(random_gen)>=0.5) vr *= -1.0;
 
-  phi  = 2.0*M_PI*Unit(random_gen);
-  cost = 2.0*(Unit(random_gen) - 0.5);
-  sint = sqrt(1.0 - cost*cost);
-  cosp = cos(phi);
-  sinp = sin(phi);
+  double phi  = 2.0*M_PI*Unit(random_gen);
+  double cost = 2.0*(Unit(random_gen) - 0.5);
+  double sint = sqrt(1.0 - cost*cost);
+  double cosp = cos(phi);
+  double sinp = sin(phi);
 
   double dfr = real->distf(eee, r*vt);
   double dfn = fake->distf(eee, r*vt);
@@ -1560,25 +1819,206 @@ Eigen::VectorXd SphericalModelMulti::gen_point(double radius, int& ierr)
   out[5] = vr * sint*sinp + vt1 * cost*sinp + vt2*cosp;
   out[6] = vr * cost      - vt1 * sint;
     
+  int ierr = 0;
   for (int i=0; i<7; i++) {
     if (std::isnan(out[i]) || std::isinf(out[i])) ierr = 1;
   }
   
-  return out;
+  return {out, ierr};
 }
 
+AxiSymModel::PSret
+SphericalModelMulti::gen_point(double radius, double theta, double phi)
+{
+  if (!real->dist_defined || !fake->dist_defined) {
+    throw std::runtime_error("SphericalModelMulti: input distribution functions must be defined before realizing!");
+  }
+
+  double r, pot, vmax, fmax;
+  double vr=0.0, vt=0.0, eee=0.0, vt1=0.0, vt2=0.0;
+
+  double Emax = get_pot(get_max_radius());
+
+  if (gen_firstime) {
+    double tol = 1.0e-5;
+    double dx = (1.0 - 2.0*tol)/(numr-1);
+    double dy = (1.0 - 2.0*tol)/(numj-1);
+    double dr;
+
+    gen_mass.resize(gen_N);
+    gen_rloc.resize(gen_N);
+    gen_fmax.resize(gen_N);
+
+    if (rmin_gen <= 1.0e-16) gen_logr = 0;
+    
+    if (gen_logr)
+      dr = (log(rmax_gen) - log(rmin_gen))/(gen_N-1);
+    else
+      dr = (rmax_gen - rmin_gen)/(gen_N-1);
+
+    for (int i=0; i<gen_N; i++) {
+
+      if (gen_logr) {
+	gen_rloc[i] = log(rmin_gen) + dr*i;
+	r = exp(gen_rloc[i]);
+      }
+      else {
+	gen_rloc[i] = rmin_gen + dr*i;
+	r = gen_rloc[i];
+      }
+
+      gen_mass[i] = fake->get_mass(r);
+
+      pot = get_pot(r);
+      vmax = sqrt(2.0*fabs(Emax - pot));
+
+      fmax = 0.0;
+      for (int j=0; j<numr; j++) {
+	double xxx = tol + dx*j;
+
+	for (int k=0; k<numj; k++) {
+	  double yyy = tol + dy*k;
+
+	  vr = vmax*xxx;
+	  vt = vmax*sqrt((1.0 - xxx*xxx)*yyy);
+	  eee = pot + 0.5*(vr*vr + vt*vt);
+
+	  double zzz = fake->distf(eee, r*vt);
+
+	  if (zzz>fmax) fmax = zzz;
+	}
+      }
+      gen_fmax[i] = fmax*(1.0 + ftol);
+
+    }
+
+    // Debug
+    //
+    ofstream test("test_multi.grid");
+    if (test) {
+
+      test << "# Rmin=" << rmin_gen
+	   << "  Rmax=" << rmax_gen
+	   << std::endl;
+
+      for (int i=0; i<gen_N; i++) {
+	test << std::setw(15) << gen_rloc[i]
+	     << std::setw(15) << gen_mass[i]
+	     << std::setw(15) << gen_fmax[i]
+	     << std::setw(15) << get_pot(exp(gen_rloc[i]))
+	     << std::setw(15) << fake->distf(get_pot(exp(gen_rloc[i])), 0.5)
+	     << std::endl;
+      }
+    }
 
-Eigen::VectorXd
-  SphericalModelMulti::gen_point(double Emin, double Emax, double Kmin,
-				 double Kmax, int& ierr)
+    gen_firstime = false;
+  }
+
+  r    = max<double>(rmin_gen, min<double>(rmax_gen, radius));
+  fmax = odd2(r, gen_rloc, gen_fmax, 1);
+  if (gen_logr) r = exp(r);
+  
+  pot  = get_pot(r);
+  vmax = sqrt(2.0*max<double>(Emax - pot, 0.0));
+
+  // Trial velocity point
+  //
+				// Diagnostic counters
+  int reject=0;
+  int it;			// Iteration counter
+  for (it=0; it<gen_itmax; it++) {
+
+    double xxx = 2.0*sin(asin(Unit(random_gen))/3.0);
+    double yyy = (1.0 - xxx*xxx)*Unit(random_gen);
+
+    vr = vmax*xxx;
+    vt = vmax*sqrt(yyy);
+    eee = pot + 0.5*(vr*vr + vt*vt);
+
+    if (fmax<=0.0) continue;
+
+    if (Unit(random_gen) > fake->distf(eee, r*vt)/fmax ) {
+      reject++;
+      continue;
+    }
+
+    if (Unit(random_gen)<0.5) vr *= -1.0;
+    
+    // Orientation of tangential velocity vector
+    //
+    double azi = 2.0*M_PI*Unit(random_gen);
+    vt1 = vt*cos(azi);
+    vt2 = vt*sin(azi);
+
+    break;
+  }
+                
+  Eigen::VectorXd out(7);
+
+  if (it==gen_itmax) {
+    if (verbose) {
+      static unsigned totcnt = 0;
+      std::cerr << "Velocity selection failed [" << std::setw(7) << ++totcnt
+		<< "," << std::setw(4) << myid << "]: r="
+		<< std::setw(12) << r
+		<< std::setw(12) << fmax
+		<< " reject="  << std::setw( 6) << reject
+		<< std::endl;
+    }
+    out.setZero();
+    return {out, 1};
+  }
+
+  if (Unit(random_gen)>=0.5) vr *= -1.0;
+
+  double dfr = real->distf(eee, r*vt);
+  double dfn = fake->distf(eee, r*vt);
+  double rat = dfr/dfn;
+
+  // Deep debug
+  //
+#ifdef DEBUG
+  if (rat <= 0.0) {
+    std::cout << "[" << std::setw(3) << myid << "] Bad mass: rat="
+	      << std::setw(16) << rat
+	      << " df(M)=" << std::setw(16) << dfr
+	      << " df(N)=" << std::setw(16) << dfn
+	      << " r="     << std::setw(16) << r
+	      << " E="     << std::setw(16) << eee
+	      << std::endl;
+  }
+#endif
+
+  double cost = cos(theta);
+  double sint = sin(theta);
+  double cosp = cos(phi);
+  double sinp = sin(phi);
+
+  out[0] = rat;
+  out[1] = r * sint*cosp;
+  out[2] = r * sint*sinp;
+  out[3] = r * cost;
+  out[4] = vr * sint*cosp + vt1 * cost*cosp - vt2*sinp;
+  out[5] = vr * sint*sinp + vt1 * cost*sinp + vt2*cosp;
+  out[6] = vr * cost      - vt1 * sint;
+    
+  int ierr = 0;
+  for (int i=0; i<7; i++) {
+    if (std::isnan(out[i]) || std::isinf(out[i])) ierr = 1;
+  }
+  
+  return {out, ierr};
+}
+
+AxiSymModel::PSret SphericalModelMulti::gen_point
+(double Emin, double Emax, double Kmin,  double Kmax)
 {
   if (!real->dist_defined || !fake->dist_defined) {
-    std::cerr << "SphericalModelMulti: input distribution functions must be defined before realizing!" << std::endl;
-    exit (-1);
+    throw std::runtime_error("SphericalModelMulti: input distribution functions must be defined before realizing!");
   }
   
   double r, vr, vt, vt1, vt2, E, K, J, jmax, w1t, pot;
-  double phi, sint, cost, sinp, cosp, azi;
+  double phi, sint, cost, sinp, cosp;
   double rmin = max<double>(get_min_radius(), gen_rmin);
   double rmax = get_max_radius();
 
@@ -1619,7 +2059,9 @@ Eigen::VectorXd
 
       vector<WRgrid> wrvec;
       for (int j=0; j<gen_K; j++) {
-				// Trapezoidal rule factor
+
+	// Trapezoidal rule factor
+	//
 	if (j==0 || j==gen_K) tfac = 0.5;
 	else tfac = 1.0;
 
@@ -1669,7 +2111,8 @@ Eigen::VectorXd
     gen_firstime_E = false;
   }
 
-				// Enforce limits
+  // Enforce limits
+  //
   Emin = max<double>(Emin, Emin_grid);
   Emin = min<double>(Emin, Emax_grid);
   Emax = max<double>(Emax, Emin_grid);
@@ -1714,9 +2157,10 @@ Eigen::VectorXd
       
   pot = get_pot(r);
   vt = J/r;
-  ierr = 0;
-				// Interpolation check (should be rare)
-				// Error condition is set
+
+  // Interpolation check (should be rare). Error condition is set.
+  //
+  int ierr = 0;
   if (2.0*(E - pot) - vt*vt < 0.0) {
     ierr = 1;
     if (E < pot) E = pot;
@@ -1726,7 +2170,9 @@ Eigen::VectorXd
 
   if (Unit(random_gen)<0.5) vr *= -1.0;
     
-  azi = 2.0*M_PI*Unit(random_gen);
+  // Orientation of the tangential velocity vector
+  //
+  double azi = 2.0*M_PI*Unit(random_gen);
   vt1 = vt*cos(azi);
   vt2 = vt*sin(azi);
 
@@ -1750,7 +2196,7 @@ Eigen::VectorXd
     if (std::isnan(out[i]) || std::isinf(out[i])) ierr = 1;
   }
 
-  return out;
+  return {out, ierr};
 }
 
 
diff --git a/include/massmodel.H b/include/massmodel.H
index 88939322d..233e2f257 100644
--- a/include/massmodel.H
+++ b/include/massmodel.H
@@ -4,6 +4,7 @@
 #include <vector>
 #include <memory>
 #include <random>
+#include <tuple>
 
 #include <Eigen/Eigen>
 
@@ -111,6 +112,10 @@ public:
 class AxiSymModel : public MassModel,
 		    public std::enable_shared_from_this<AxiSymModel>
 {
+public:
+
+  using PSret = std::tuple<Eigen::VectorXd, int>;
+
 protected:
   //@{
   //! Stuff for gen_point
@@ -123,14 +128,18 @@ protected:
   bool gen_firstime_jeans;
   Eigen::VectorXd gen_rloc, gen_mass, gen_fmax;
   SphericalOrbit gen_orb;
-  double gen_fomax, gen_ecut;
+  double gen_fomax, gen_ecut, gen_lastr=-1.0;
   //@}
 
-  Eigen::VectorXd gen_point_2d(int& ierr);
-  Eigen::VectorXd gen_point_2d(double r, int& ierr);
-  Eigen::VectorXd gen_point_3d(int& ierr);
-  Eigen::VectorXd gen_point_3d(double Emin, double Emax, double Kmin, double Kmax, int& ierr);
-  Eigen::VectorXd gen_point_jeans_3d(int& ierr);
+  PSret gen_point_2d();
+  PSret gen_point_2d(double r);
+  PSret gen_point_3d();
+  PSret gen_point_3d(double Emin, double Emax, double Kmin, double Kmax);
+  PSret gen_point_jeans_3d();
+  PSret gen_point_3d(double r, double theta, double phi);
+  PSret gen_point_3d_iso(double r, double theta, double phi,
+			 int N, int nv, int na,
+			 double PHI=0, double THETA=0, double PSI=0);
   
   double Emin_grid, Emax_grid, dEgrid, dKgrid;
   vector<double> Egrid, Kgrid, EgridMass;
@@ -212,55 +221,89 @@ public:
   void set_Ecut(double cut) { gen_ecut = cut; }
 
   //! Generate a phase-space point
-  virtual Eigen::VectorXd gen_point(int& ierr) {
+  virtual PSret gen_point()
+  {
     if (dof()==2)
-      return gen_point_2d(ierr);
+      return gen_point_2d();
     else if (dof()==3)
-      return gen_point_3d(ierr);
+      return gen_point_3d();
     else
       bomb( "dof must be 2 or 3" );
     
-    return Eigen::VectorXd();
+    return {Eigen::VectorXd(), 1};
   }
   
   //! Generate a phase-space point using Jeans' equations
-  virtual Eigen::VectorXd gen_point_jeans(int& ierr) {
+  virtual PSret gen_point_jeans()
+  {
     if (dof()==2)
       bomb( "AxiSymModel: gen_point_jeans_2d(ierr) not implemented!" );
     else if (dof()==3)
-      return gen_point_jeans_3d(ierr);
+      return gen_point_jeans_3d();
     else
       bomb( "dof must be 2 or 3" );
     
-    return Eigen::VectorXd();
+    return {Eigen::VectorXd(), 1};
   }
   
   //! Generate a phase-space point at a particular radius
-  virtual Eigen::VectorXd gen_point(double r, int& ierr) {
+  virtual PSret gen_point(double r) {
+    if (dof()==2)
+      return gen_point_2d(r);
+    else if (dof()==3)
+      bomb( "AxiSymModel: gen_point(r) not implemented!" );
+    else
+      bomb( "AxiSymModel: dof must be 2 or 3" );
+    
+    return {Eigen::VectorXd(), 1};
+  }
+  
+  //! Generate a phase-space point at a particular radius, theta, and phi
+  virtual PSret
+  gen_point(double r, double theta, double phi)
+  {
     if (dof()==2)
-      return gen_point_2d(r, ierr);
+      bomb( "AxiSymModel: gen_point(r, theta, phi) not implemented!" );
     else if (dof()==3)
-      bomb( "AxiSymModel: gen_point(r, ierr) not implemented!" );
+      return gen_point_3d(r, theta, phi);
     else
       bomb( "AxiSymModel: dof must be 2 or 3" );
     
-    return Eigen::VectorXd();
+    return {Eigen::VectorXd(), 1};
+  }
+  
+  //! Generate a phase-space point at a particular radius, theta, and phi
+  virtual PSret
+  gen_point_iso(double r, double theta, double phi, int N, int nv, int na,
+		double PHI, double THETA, double PSI)
+  {
+    if (dof()==2)
+      bomb( "AxiSymModel: gen_point(r, theta, phi) not implemented!" );
+    else if (dof()==3)
+      return gen_point_3d_iso(r, theta, phi, N, nv, na,
+			      PHI, THETA, PSI);
+    else
+      bomb( "AxiSymModel: dof must be 2 or 3" );
+    
+    return {Eigen::VectorXd(), 1};
   }
   
   //! Generate a phase-space point at a particular energy and angular momentum
-  virtual Eigen::VectorXd gen_point(double Emin, double Emax, double Kmin, double Kmax, int& ierr) {
+  virtual PSret
+  gen_point(double Emin, double Emax, double Kmin, double Kmax)
+  {
     if (dof()==2)
-      bomb( "AxiSymModel: gen_point(r, ierr) not implemented!" );
+      bomb( "AxiSymModel: gen_point(r) not implemented!" );
     else if (dof()==3)
-      return gen_point_3d(Emin, Emax, Kmin, Kmax, ierr);
+      return gen_point_3d(Emin, Emax, Kmin, Kmax);
     else
       bomb( "AxiSymModel: dof must be 2 or 3" );
     
-    return Eigen::VectorXd();
+    return {Eigen::VectorXd(), 1};
   }
   
   //! Generate a the velocity variate from a position
-  virtual void gen_velocity(double *pos, double *vel, int& ierr);
+  virtual int gen_velocity(double *pos, double *vel);
   
 };
 
@@ -501,9 +544,10 @@ public:
 
   //@{
   //! Overloaded to provide mass distribution from Real and Number distribution from Fake
-  Eigen::VectorXd gen_point(int& ierr);
-  Eigen::VectorXd gen_point(double r, int& ierr);
-  Eigen::VectorXd gen_point(double Emin, double Emax, double Kmin, double Kmax, int& ierr);
+  PSret gen_point();
+  PSret gen_point(double r);
+  PSret gen_point(double Emin, double Emax, double Kmin, double Kmax);
+  PSret gen_point(double r, double theta, double phi);
   //@}
 
 
diff --git a/src/Basis.cc b/src/Basis.cc
index e7f7a6c6a..1db4dc2da 100644
--- a/src/Basis.cc
+++ b/src/Basis.cc
@@ -4,7 +4,7 @@
 
 // Machine constant for Legendre
 //
-constexpr double MINEPS = 20.0*std::numeric_limits<double>::min();
+constexpr double MINEPS = 3.0*std::numeric_limits<double>::epsilon();
 
 Basis::Basis(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf)
 {
diff --git a/src/cudaSphericalBasis.cu b/src/cudaSphericalBasis.cu
index 38d8bf008..2acb86066 100644
--- a/src/cudaSphericalBasis.cu
+++ b/src/cudaSphericalBasis.cu
@@ -19,7 +19,7 @@
 // Global symbols for coordinate transformation in SphericalBasis
 //
 __device__ __constant__
-cuFP_t sphScale, sphRscale, sphHscale, sphXmin, sphXmax, sphDxi, sphCen[3];
+cuFP_t sphScale, sphRscale, sphHscale, sphXmin, sphXmax, sphDxi, sphCen[3], sphEPS;
 
 __device__ __constant__
 int    sphNumr, sphCmap;
@@ -80,7 +80,7 @@ void legendre_v(int lmax, cuFP_t x, cuFP_t* p)
 }
 
 
-__host__ __device__
+__device__
 void legendre_v2(int lmax, cuFP_t x, cuFP_t* p, cuFP_t* dp)
 {
   cuFP_t fact, somx2, pll, pl1, pl2;
@@ -106,9 +106,9 @@ void legendre_v2(int lmax, cuFP_t x, cuFP_t* p, cuFP_t* dp)
     }
   }
 
-  if (1.0-fabs(x) < MINEPS) {
-    if (x>0) x =   1.0 - MINEPS;
-    else     x = -(1.0 - MINEPS);
+  if (1.0-fabs(x) < sphEPS) {
+    if (x>0) x =   1.0 - sphEPS;
+    else     x = -(1.0 - sphEPS);
   }
 
   somx2 = 1.0/(x*x - 1.0);
@@ -134,6 +134,7 @@ void testConstantsSph()
   printf("   Dxi    = %f\n", sphDxi   );
   printf("   Numr   = %d\n", sphNumr  );
   printf("   Cmap   = %d\n", sphCmap  );
+  printf("   Feps   = %d\n", sphEPS   );
   printf("-------------------------\n");
 }
 
@@ -231,6 +232,15 @@ void SphericalBasis::initialize_mapping_constants()
 		 __FILE__, __LINE__, "Error copying sphEVEN_M");
   cuda_safe_call(cudaMemcpyToSymbol(sphM0only, &M0_only,  sizeof(bool),  size_t(0), cudaMemcpyHostToDevice),
 		 __FILE__, __LINE__, "Error copying sphM0only");
+
+#if cuREAL == 4
+  cuFP_t feps = 3.0*std::numeric_limits<float>::epsilon();
+#else
+  cuFP_t feps = 3.0*std::numeric_limits<double>::epsilon();
+#endif
+
+  cuda_safe_call(cudaMemcpyToSymbol(sphEPS, &feps,  sizeof(cuFP_t),  size_t(0), cudaMemcpyHostToDevice),
+		 __FILE__, __LINE__, "Error copying sphEPS");
 }
 
 
diff --git a/utils/ICs/CMakeLists.txt b/utils/ICs/CMakeLists.txt
index 0c89b9942..281695e52 100644
--- a/utils/ICs/CMakeLists.txt
+++ b/utils/ICs/CMakeLists.txt
@@ -2,7 +2,7 @@
 set(bin_PROGRAMS shrinkics gensph gendisk gendisk2d gsphere pstmod
   empinfo empdump eofbasis eofcomp testcoefs testcoefs2 testdeval
   forcetest hdf52accel forcetest2 cylcache modelfit test2d addsphmod
-  cubeics zangics slabics)
+  cubeics zangics slabics addring)
 
 set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX yaml-cpp exputil
   ${VTK_LIBRARIES} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES}
@@ -88,6 +88,8 @@ add_executable(zangics ZangICs.cc)
 
 add_executable(slabics genslab.cc massmodel1d.cc)
 
+add_executable(addring addring.cc)
+
 foreach(program ${bin_PROGRAMS})
   target_link_libraries(${program} ${common_LINKLIB})
   target_include_directories(${program} PUBLIC ${common_INCLUDE})
diff --git a/utils/ICs/Disk2dHalo.cc b/utils/ICs/Disk2dHalo.cc
index f923e4792..ad0e9a61b 100644
--- a/utils/ICs/Disk2dHalo.cc
+++ b/utils/ICs/Disk2dHalo.cc
@@ -401,7 +401,7 @@ void Disk2dHalo::set_halo(vector<Particle>& phalo, int nhalo, int npart)
   for (int i=0; i<npart; i++) {
 
     do {
-      ps = multi->gen_point(ierr);
+      std::tie(ps, ierr) = multi->gen_point();
       if (ierr) count1++;
     } while (ierr);
     
@@ -2329,7 +2329,7 @@ void Disk2dHalo::set_vel_halo(vector<Particle>& part)
 				// Use Eddington
     
     if (DF && 0.5*(1.0+erf((r-R_DF)/DR_DF)) > rndU(gen)) {
-      halo2->gen_velocity(&p.pos[0], &p.vel[0], nok);
+      nok = halo2->gen_velocity(&p.pos[0], &p.vel[0]);
       
       if (nok) {
 	std::cout << "gen_velocity failed: "
diff --git a/utils/ICs/DiskHalo.cc b/utils/ICs/DiskHalo.cc
index c64a75c6a..0ff393af6 100644
--- a/utils/ICs/DiskHalo.cc
+++ b/utils/ICs/DiskHalo.cc
@@ -411,7 +411,7 @@ void DiskHalo::set_halo(vector<Particle>& phalo, int nhalo, int npart)
   for (int i=0; i<npart; i++) {
 
     do {
-      ps = multi->gen_point(ierr);
+      std::tie(ps, ierr) = multi->gen_point();
       if (ierr) count1++;
     } while (ierr);
     
@@ -631,7 +631,7 @@ set_halo_coordinates(vector<Particle>& phalo, int nhalo, int npart)
     p.pos[2] = r*costh;
 
     do {
-      ps = halo2->gen_point(ierr);
+      std::tie(ps, ierr) = halo2->gen_point();
     } while (ierr);
 
     massp1 += p.mass;
@@ -2560,7 +2560,8 @@ void DiskHalo::set_vel_halo(vector<Particle>& part)
 				// Use Eddington
     
     if (DF && 0.5*(1.0+erf((r-R_DF)/DR_DF)) > rndU(gen)) {
-      halo2->gen_velocity(&p.pos[0], &p.vel[0], nok);
+
+      nok = halo2->gen_velocity(&p.pos[0], &p.vel[0]);
       
       if (nok) {
 	std::cout << "gen_velocity failed: "
diff --git a/utils/ICs/ZangICs.cc b/utils/ICs/ZangICs.cc
index 638c6c755..660a7633c 100644
--- a/utils/ICs/ZangICs.cc
+++ b/utils/ICs/ZangICs.cc
@@ -27,6 +27,7 @@ main(int ac, char **av)
   //=====================
 
   int          N;		// Number of particles
+  int          Nrepl;		// Number of particle replicates per orbit
   double       mu, nu, Ri, Ro;	// Taper paramters
   double       Rmin, Rmax;      // Radial range
   double       sigma;		// Velocity dispersion
@@ -59,6 +60,8 @@ main(int ac, char **av)
      cxxopts::value<double>(sigma)->default_value("1.0"))
     ("s,seed",    "Random number seed. Default: use /dev/random",
      cxxopts::value<unsigned>(seed))
+    ("q,Nrepl",   "Number of particle replicates per orbit",
+     cxxopts::value<int>(Nrepl)->default_value("1"))
     ("f,file",    "Output body file",
      cxxopts::value<std::string>(bodyfile)->default_value("zang.bods"))
     ;
@@ -84,6 +87,10 @@ main(int ac, char **av)
     seed = std::random_device{}();
   }
 
+  // Set particle number consitent with even replicants
+  if (Nrepl<1) Nrepl = 1;
+  if (Nrepl>1) N = (N/Nrepl)*Nrepl;
+
   // Make sure N>0
   if (N<=0) {
     std::cerr << av[0] << ": you must requiest at least one body"
@@ -107,6 +114,11 @@ main(int ac, char **av)
 						   1.0, Rmin, Rmax);
   model->setup_df(sigma);
 
+  // Replicant logic
+  //
+  int Number = N/Nrepl;
+  double dPhi = 2.0*M_PI/Nrepl;
+
   // Progress bar
   //
   std::shared_ptr<progress::progress_display> progress;
@@ -115,7 +127,7 @@ main(int ac, char **av)
   {
     nomp = omp_get_num_threads();
     if (omp_get_thread_num()==0) {
-      progress = std::make_shared<progress::progress_display>(N/nomp);
+      progress = std::make_shared<progress::progress_display>(Number/nomp);
     }
   }
 
@@ -190,7 +202,7 @@ main(int ac, char **av)
   // Generation loop with OpenMP
   //
 #pragma omp parallel for reduction(+:over)
-  for (int n=0; n<N; n++) {
+  for (int n=0; n<Number; n++) {
     // Thread id
     int tid = omp_get_thread_num();
 
@@ -224,23 +236,27 @@ main(int ac, char **av)
 
     if (w1 > M_PI) vr *= -1.0;	// Branch of radial motion
 
-    // Convert from polar to Cartesian
-    //
-    pos[n][0] = r*cos(phi);
-    pos[n][1] = r*sin(phi);
-    pos[n][2] = 0.0;
-
-    vel[n][0] = vr*cos(phi) - vt*sin(phi);
-    vel[n][1] = vr*sin(phi) + vt*cos(phi);
-    vel[n][2] = 0.0;
-
-    // Accumulate mean position and velocity
-    //
-    for (int k=0; k<3; k++) {
-      zeropos[tid][k] += pos[n][k];
-      zerovel[tid][k] += vel[n][k];
+    for (int nn=0; nn<Nrepl; nn++) {
+      int indx = n*Nrepl + nn;
+      double Phi = phi + dPhi*nn;
+      // Convert from polar to Cartesian
+      //
+      pos[indx][0] = r*cos(Phi);
+      pos[indx][1] = r*sin(Phi);
+      pos[indx][2] = 0.0;
+
+      vel[indx][0] = vr*cos(Phi) - vt*sin(Phi);
+      vel[indx][1] = vr*sin(Phi) + vt*cos(Phi);
+      vel[indx][2] = 0.0;
+
+      // Accumulate mean position and velocity
+      //
+      for (int k=0; k<3; k++) {
+	zeropos[tid][k] += pos[indx][k];
+	zerovel[tid][k] += vel[indx][k];
+      }
     }
-
+      
     // Print progress bar
     if (tid==0) ++(*progress);
   }
diff --git a/utils/ICs/addring.cc b/utils/ICs/addring.cc
new file mode 100644
index 000000000..6e4cedafe
--- /dev/null
+++ b/utils/ICs/addring.cc
@@ -0,0 +1,192 @@
+/*
+  Add a ring of particles to an existing realization
+*/
+
+#include <filesystem>
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <sstream>
+#include <cstdlib>
+#include <string>
+#include <random>
+#include <vector>
+#include <cmath>
+
+#include <Eigen/Eigen>
+
+#include <fenv.h>
+
+#include <config_exp.h>
+#ifdef HAVE_OMP_H
+#include <omp.h>
+#endif
+
+// EXP classes
+//
+#include <numerical.H>
+#include <gaussQ.H>
+#include <isothermal.H>
+#include <hernquist_model.H>
+#include <model3d.H>
+#include <biorth.H>
+#include <SphericalSL.H>
+#include <interp.H>
+#include <AddSpheres.H>
+#include <libvars.H>		// Library globals
+#include <cxxopts.H>		// Command-line parsing
+#include <EXPini.H>		// Ini-style config
+
+                                // Local headers
+#include <SphericalSL.H>
+#include <localmpi.H>
+
+int
+main(int ac, char **av)
+{
+  //====================
+  // Inialize MPI stuff
+  //====================
+
+  local_init_mpi(ac, av);
+
+  //====================
+  // Begin opt parsing
+  //====================
+
+  std::string  halofile;
+  int          DIVERGE, N;
+  double       DIVERGE_RFAC, mass, radius, width, vdisp;
+  std::string  outfile, config;
+
+  const std::string mesg("Make a model from the combination of two spherical models.  E.g. a DM halo and a bulge.\n");
+
+  cxxopts::Options options(av[0], mesg);
+
+  options.add_options()
+    ("h,help", "Print this help message")
+    ("s,spline", "Set interplation in SphericalModelTable to 'spline' (default: linear)")
+    ("T,template", "Write template options file with current and all default values")
+    ("c,config", "Parameter configuration file",
+     cxxopts::value<string>(config))
+    ("N,number", "Number of particles",
+     cxxopts::value<int>(N)->default_value("100000"))
+    ("DIVERGE", "Cusp power-law density extrapolation (0 means off)",
+     cxxopts::value<int>(DIVERGE)->default_value("0"))
+    ("DIVERGE_RFAC", "Exponent for inner cusp extrapolation",
+     cxxopts::value<double>(DIVERGE_RFAC)->default_value("1.0"))
+    ("i,halofile", "Halo spherical model table",
+     cxxopts::value<string>(halofile)->default_value("SLGridSph.one"))
+    ("o,outfile", "Output particle file",
+     cxxopts::value<string>(outfile)->default_value("addring.bods"))
+    ("M,mass", "Mass of ring",
+     cxxopts::value<double>(mass)->default_value("0.1"))
+    ("R,radius", "Radius of ring",
+     cxxopts::value<double>(radius)->default_value("1.0"))
+    ("W,width", "Width of ring",
+     cxxopts::value<double>(width)->default_value("0.1"))
+    ("V,disp", "Velocity dispersion",
+     cxxopts::value<double>(vdisp)->default_value("0.1"))
+    ;
+
+  cxxopts::ParseResult vm;
+
+  try {
+    vm = options.parse(ac, av);
+  } catch (cxxopts::OptionException& e) {
+    if (myid==0) std::cout << "Option error: " << e.what() << std::endl;
+    MPI_Finalize();
+    exit(-1);
+  }
+
+  // Write YAML template config file and exit
+  //
+  if (vm.count("template")) {
+    // Write template file
+    //
+    if (myid==0) SaveConfig(vm, options, "template.yaml");
+
+    MPI_Finalize();
+    return 0;
+  }
+
+  // Print help message and exit
+  //
+  if (vm.count("help")) {
+    if (myid == 0) {
+      std::cout << options.help() << std::endl << std::endl
+		<< "Examples: " << std::endl
+		<< "\t" << "Use parameters read from a YAML config file"  << std::endl
+		<< "\t" << av[0] << " --config=addsphmod.config"  << std::endl << std::endl
+		<< "\t" << "Generate a template YAML config file from current defaults"  << std::endl
+		<< "\t" << av[0] << " --template=template.yaml" << std::endl << std::endl;
+    }
+    MPI_Finalize();
+    return 0;
+  }
+
+  // Read parameters fron the YAML config file
+  //
+  if (vm.count("config")) {
+    try {
+      vm = LoadConfig(options, config);
+    } catch (cxxopts::OptionException& e) {
+      if (myid==0) std::cout << "Option error in configuration file: "
+			     << e.what() << std::endl;
+      MPI_Finalize();
+      return 0;
+    }
+  }
+
+  if (vm.count("spline")) {
+    SphericalModelTable::linear = 0;
+  } else {
+    SphericalModelTable::linear = 1;
+  }
+
+  // Open output file
+  //
+  std::ofstream out(outfile);
+  if (not out) throw std::runtime_error("Cannot open output file " + outfile);
+
+  // Model
+  //
+  SphericalModelTable model(halofile, DIVERGE, DIVERGE_RFAC);
+
+
+  // Random setup
+  //
+  std::random_device rd{};
+  std::mt19937 gen{rd()};
+ 
+  // Unit normal distribution
+  //
+  std::normal_distribution<> d{0.0, 1.0};
+  std::uniform_real_distribution<> u{0.0, 1.0};
+  
+ 
+  double pmass = mass/N;
+  Eigen::Vector3d pos, vel;
+  
+  for (int i=0; i<N; i++) {
+
+    double R  = radius + 0.5*width*d(gen);
+    double v0 = sqrt(model.get_mass(radius)/radius);
+    double phi = 2.0*M_PI*u(gen);
+
+    pos(0) = R*cos(phi);
+    pos(1) = R*sin(phi);
+    pos(2) = 0.5*width*d(gen);
+    
+    vel(0) = vdisp*d(gen) - v0*sin(phi);
+    vel(1) = vdisp*d(gen) + v0*cos(phi);
+    vel(2) = vdisp*d(gen);
+    
+    out << std::setw(18) << pmass;
+    for (int k=0; k<3; k++) out << std::setw(18) << pos(k);
+    for (int k=0; k<3; k++) out << std::setw(18) << vel(k);
+    out << std::endl;
+  }
+
+  return 0;
+}
diff --git a/utils/ICs/gensph.cc b/utils/ICs/gensph.cc
index b42b397bc..045d48f28 100644
--- a/utils/ICs/gensph.cc
+++ b/utils/ICs/gensph.cc
@@ -61,13 +61,16 @@ using namespace __EXP__;
 #include <fftw3.h>
 #endif
 
+#include <euler.H>
+
 // Global variables
 
 int
 main(int argc, char **argv)
 {
   int HMODEL, N, NUMDF, NUMR, NUMJ, NUME, NUMG, NREPORT, SEED, ITMAX;
-  int  NUMMODEL, RNUM, DIVERGE, DIVERGE2, LINEAR, NUMINT, NI, ND;
+  int NUMMODEL, RNUM, DIVERGE, DIVERGE2, LINEAR, NUMINT, NI, ND;
+  int Nrepl=1, Nfib=1;
   double DIVERGE_RFAC, DIVERGE_RFAC2, NN, MM, RA, RMODMIN, RMOD, EPS;
   double X0, Y0, Z0, U0, V0, W0, TOLE;
   double Emin0, Emax0, Kmin0, Kmax0, RBAR, MBAR, BRATIO, CRATIO, SMOOTH;
@@ -75,6 +78,8 @@ main(int argc, char **argv)
   bool VTEST;
   std::string INFILE, MMFILE, OUTFILE, OUTPS, config;
 
+  const double goldenRatio = 1.618033988749895;
+
 #ifdef DEBUG
   set_fpu_handler();
 #endif
@@ -179,6 +184,10 @@ main(int argc, char **argv)
      cxxopts::value<bool>(ELIMIT)->default_value("false"))
     ("VTEST", "Test gen_velocity() generation",
      cxxopts::value<bool>(VTEST)->default_value("false"))
+    ("Nrepl", "Number of replicates in orbital plane (1 skips the Sellwood 1997 algorithm)",
+     cxxopts::value<int>(Nrepl)->default_value("1"))
+    ("Nfib", "Number of points on the sphere for each orbit. Replicate the orbits by tiling the angular momentum direction on the sphere using a Fibonnaci sequence.  Default value of 1 implies one plane per orbit.",
+     cxxopts::value<int>(Nfib)->default_value("1"))
     ("Emin0", "Minimum energy (if ELIMIT=true)",
      cxxopts::value<double>(Emin0)->default_value("-3.0"))
     ("Emax0", "Maximum energy (if ELIMIT=true)",
@@ -254,6 +263,11 @@ main(int argc, char **argv)
   if (vm.count("verbose")) VERBOSE = true;
   else                     VERBOSE = false;
 
+  // Use Fibonnaci for space-filling grid points
+  //
+  Nfib  = std::max<int>(1, Nfib );
+  Nrepl = std::max<int>(1, Nrepl);
+
   // Prepare output streams and create new files
   //
   std::ostringstream sout;
@@ -520,6 +534,22 @@ main(int argc, char **argv)
 
   }
 
+  int nfib=0, nangle=0, nshell=0;
+
+  // For replication algorithm; Nrepl=1 is the standard algorithm.
+  //
+  Nrepl *= Nfib;
+  int nplane = N/Nrepl;
+  N = Nrepl * nplane;
+  NREPORT = std::max<int>(1, NREPORT/Nrepl);
+
+  if (Nrepl > 1 and myid==0) {
+    std::cout << std::setw(60) << std::setfill('-') << '-'
+	      << std::setfill(' ') << std::endl
+	      << "Replication: N=" << N << " Nrepl=" << Nrepl/Nfib
+	      << " Nfib=" << Nfib << " Norbits=" << nplane << std::endl;
+  }
+
   double mass = hmodel->get_mass(hmodel->get_max_radius())/N;
   AxiSymModPtr rmodel = hmodel;
 
@@ -691,73 +721,195 @@ main(int argc, char **argv)
   double TT=0.0, WW=0.0, VC=0.0;
 
   if (myid==0)
-    std::cout << "-----------" << std::endl
+    std::cout << std::setw(60) << std::setfill('-') << '-' << std::endl
+	      << std::setfill('-')
 	      << "Body count:" << std::endl;
 
-  int npernode = N/numprocs;
+  int npernode = nplane/numprocs;
   int beg = myid*npernode;
   int end = beg + npernode;
 
-  if (myid==numprocs-1) end = N;
+  if (myid==numprocs-1) end = nplane;
 
   std::vector<Eigen::VectorXd> PS;
   Eigen::VectorXd zz = Eigen::VectorXd::Zero(7);
 
+  std::vector<double> rmass;
+
   for (int n=beg; n<end; n++) {
 
     do {
       if (ELIMIT)
-	ps = rmodel->gen_point(Emin0, Emax0, Kmin0, Kmax0, ierr);
+	std::tie(ps, ierr) = rmodel->gen_point(Emin0, Emax0, Kmin0, Kmax0);
       else if (VTEST) {
-	ps = rmodel->gen_point(ierr);
-	rmodel->gen_velocity(&ps[1], &ps[4], ierr);
-	if (ierr) {
-	  std::cout << "gen_velocity failed: "
-		    << ps[0] << " "
-		    << ps[1] << " "
-		    << ps[2] << "\n";
+	std::tie(ps, ierr) = rmodel->gen_point();
+	if (ierr==0) {
+	  ierr = rmodel->gen_velocity(&ps[1], &ps[4]);
+	  if (ierr) {
+	    std::cout << "gen_velocity failed: "
+		      << ps[0] << " "
+		      << ps[1] << " "
+		      << ps[2] << "\n";
+	  }
 	}
       }
       else
-	ps = rmodel->gen_point(ierr);
+	std::tie(ps, ierr) = rmodel->gen_point();
       if (ierr) count++;
     } while (ierr);
 
     if (ps[0] <= 0.0) negms++;
 
-    double RR=0.0;
-    for (int i=1; i<=3; i++) {
-      RR += ps[i]*ps[i];
-      TT += 0.5*mass*ps[0]*ps[3+i]*ps[3+i];
-    }
-    RR  =  sqrt(RR);
-    if (RR>=rmin) {
-      VC += -mass*ps[0]*rmodel->get_mass(RR)/RR;
-      WW +=  0.5*mass*ps[0]*rmodel->get_pot(RR);
-    }
+    Eigen::Vector3d pos(ps[1], ps[2], ps[3]), pos1;
+    Eigen::Vector3d vel(ps[4], ps[5], ps[6]), vel1;
+    Eigen::Matrix3d proj, Iprj, rot = Eigen::Matrix3d::Identity();
+    Eigen::Matrix3d projM, IprjM;
+
+    double dq = 2.0*M_PI * Nfib / Nrepl, mfac = ps[0];
 
-    if (zeropos or zerovel) {
-      ps[0] *= mass;
-      PS.push_back(ps);
-      zz[0] += ps[0];
-      if (zeropos) for (int i=1; i<3; i++) zz[i] -= ps[0]*ps[i];
-      if (zerovel) for (int i=4; i<7; i++) zz[i] -= ps[0]*ps[i];
+    // Compute orbital plane basis
+    //
+    if (Nrepl>1) {
+      auto L = pos.cross(vel);	// The angular momentum vector
+      if (pos.norm() < 1.0e-10 or L.norm() < 1.0e-10) {
+	proj = Iprj = Eigen::Matrix3d::Identity();
+      } else {
+	auto X = pos/pos.norm(); // Radial unit vector, X.
+	auto Y = L.cross(X);	 // Choose a right-hand3d coordinate
+	Y = Y/Y.norm();		 // system perpendicular to X and L.
+	auto Z = L/L.norm();     // Unit vector in the ang. mom. direction.
+	
+	// Azimuth of vertical plane containing L
+	double Phi   = atan2(Z(1), Z(0)) + 0.5*M_PI;
+
+	// Tilt of plane w.r.t to polar axis
+	double Theta = acos(Z(2));
+
+	// Location of X in orbital plane w.r.t line of nodes
+	Eigen::Vector3d lon(cos(Phi), sin(Phi), 0.0);
+	double dotp = lon.dot(X);
+	double Psi = acos(dotp);
+	auto dir = lon.cross(X);
+
+	proj.row(0) = X;	// Project to orbital plane frame
+	proj.row(1) = Y;
+	proj.row(2) = Z;
+	
+	Iprj.col(0) = X;	// Project from orbital plane to original
+	Iprj.col(1) = Y;
+	Iprj.col(2) = Z;
+
+	if (Nfib>1) {
+				// Quadrant geometry
+	  int sgn1 = 1, sgn2 = 1;
+	  if (Theta > 0.5*M_PI) sgn1 = -1;
+	  if (dir(2) < 0.0)     sgn2 = -1;
+	  
+				// Sanity check for debugging
+	  if (true) {
+	    double norm = 0.0;
+	    if (sgn1*sgn2==1)
+	      norm = (proj - return_euler(Phi, Theta,  Psi, 0)).norm();
+	    else
+	      norm = (proj - return_euler(Phi, Theta,  -Psi, 0)).norm();
+	    
+	    if (fabs(norm) > 1.0e-10) {
+	      std::cout << "Theta=" << Theta
+			<< " dir(2)=" << dir(2) << std::endl;
+	    }
+	  }
+
+	  if (sgn1*sgn2==1) {
+	    projM = return_euler(Phi, -Theta,  Psi, 0);
+	    IprjM = return_euler(Phi, -Theta,  Psi, 1);
+	  } else {
+	    projM = return_euler(Phi, -Theta,  -Psi, 0);
+	    IprjM = return_euler(Phi, -Theta,  -Psi, 1);
+	  }	  
+	}
+      }
     }
-    else {
-      out << std::setw(20) << mass * ps[0];
-      for (int i=1; i<=6; i++) out << std::setw(20) << ps[i]+ps0[i];
 
-      if (NI) {
-	for (int n=0; n<NI; n++) out << std::setw(10) << 0;
+    // Orbital plane loop
+    //
+    for (int q=0; q<Nrepl; q++) {
+
+      double RR=0.0;
+      for (int i=1; i<=3; i++) {
+	RR += ps[i]*ps[i];
+	TT += 0.5*mass*ps[0]*ps[3+i]*ps[3+i];
       }
-      if (ND) {
-	for (int n=0; n<ND; n++) out << std::setw(20) << 0.0;
+      RR  =  sqrt(RR);
+      if (RR>=rmin) {
+	VC += -mass*ps[0]*rmodel->get_mass(RR)/RR;
+	WW +=  0.5*mass*ps[0]*rmodel->get_pot(RR);
       }
 
-      out << std::endl;
+      if (zeropos or zerovel) {
+	ps[0] *= mass;
+	PS.push_back(ps);
+	zz[0] += ps[0];
+	if (zeropos) for (int i=1; i<3; i++) zz[i] -= ps[0]*ps[i];
+	if (zerovel) for (int i=4; i<7; i++) zz[i] -= ps[0]*ps[i];
+      }
+      else {
+	out << std::setw(20) << mass * ps[0];
+	for (int i=1; i<=6; i++) out << std::setw(20) << ps[i]+ps0[i];
+	
+	if (NI) {
+	  for (int n=0; n<NI; n++) out << std::setw(10) << 0;
+	}
+	if (ND) {
+	  for (int n=0; n<ND; n++) out << std::setw(20) << 0.0;
+	}
+
+	out << std::endl;
+      }
+
+      // Rotate particle in orbital plane
+      if (Nrepl>1) {
+
+	Eigen::Matrix3d invt;
+	double Q = dq*(q/Nfib+1), phi=0.0, cost=0.0;
+	int nfib = q % Nfib;
+
+	if (Nfib>1) {
+	  phi  = 2.0*M_PI * nfib / goldenRatio;
+	  cost = 1.0 - 2.0 * nfib/Nfib;
+	  cost = (cost < -1.0 ? -1.0 : (cost > 1.0 ? 1.0 : cost));
+	  invt = return_euler(phi, acos(cost),  0, 1);
+	}
+
+	// Update rotation matrix
+	rot(0, 0) = rot(1, 1) = cos(Q);
+	rot(0, 1) = -sin(Q);
+	rot(1, 0) =  sin(Q);
+
+	// Full rotation matrix
+	Eigen::Matrix3d trans;
+	if (Nfib>1)  trans = invt * rot * proj;
+	else         trans = Iprj * rot * proj;
+
+	// Rotate position and velocity
+	pos1 =  trans * pos;
+	vel1 = -trans * vel;
+
+	// Reset the phase-space vector
+	//
+	ps[0] = mfac;
+	
+	ps[1] = pos1(0);
+	ps[2] = pos1(1);
+	ps[3] = pos1(2);
+
+	ps[4] = vel1(0);
+	ps[5] = vel1(1);
+	ps[6] = vel1(2);
+      }
     }
 
-    if (myid==0 and !((n+1)%NREPORT)) cout << '\r' << (n+1)*numprocs << flush;
+    if (myid==0 and !((n+1)%NREPORT)) cout << '\r' << (n+1)*numprocs*Nrepl
+					   << flush;
   }
 
   if (zeropos or zerovel) {
@@ -797,16 +949,20 @@ main(int argc, char **argv)
     MPI_Reduce(MPI_IN_PLACE, &WW,    1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
     MPI_Reduce(MPI_IN_PLACE, &VC,    1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
 
-    cout << std::endl << "-----------" << std::endl << std::endl;
-    cout << std::setw(20) << "States rejected: " << count      << std::endl;
-    cout << std::setw(20) << "Negative masses: " << negms      << std::endl;
-    cout << std::setw(60) << setfill('-') << '-' << std::endl       << setfill(' ');
-    cout << std::setw(20) << "KE="               << TT         << std::endl;
-    cout << std::setw(20) << "PE="               << WW         << std::endl;
-    cout << std::setw(20) << "VC="               << VC         << std::endl;
-    cout << std::setw(20) << "Ratio (-2T/W)="    << -2.0*TT/WW << std::endl;
-    cout << std::setw(20) << "Ratio (-2T/C)="    << -2.0*TT/VC << std::endl;
-    std::cout << std::setw(60) << std::setfill('-') << '-' << std::endl << std::setfill(' ');
+    std::cout << std::endl
+	      << std::setw(60) << std::setfill('-') << '-' << std::endl
+	      << std::setfill(' ')
+	      << std::setw(20) << "States rejected: " << count      << std::endl
+	      << std::setw(20) << "Negative masses: " << negms      << std::endl
+	      << std::setw(60) << setfill('-') << '-' << std::endl
+	      << setfill(' ')
+	      << std::setw(20) << "KE="               << TT         << std::endl
+	      << std::setw(20) << "PE="               << WW         << std::endl
+	      << std::setw(20) << "VC="               << VC         << std::endl
+	      << std::setw(20) << "Ratio (-2T/W)="    << -2.0*TT/WW << std::endl
+	      << std::setw(20) << "Ratio (-2T/C)="    << -2.0*TT/VC << std::endl
+	      << std::setw(60) << std::setfill('-') << '-' << std::endl
+	      << std::setfill(' ');
   }
 
   out.close();