diff --git a/expui/CMakeLists.txt b/expui/CMakeLists.txt
index 8690817f7..d506bea8d 100644
--- a/expui/CMakeLists.txt
+++ b/expui/CMakeLists.txt
@@ -34,7 +34,7 @@ endif()
 set(expui_SOURCES BasisFactory.cc BiorthBasis.cc FieldBasis.cc
   CoefContainer.cc CoefStruct.cc FieldGenerator.cc expMSSA.cc
   Coefficients.cc KMeans.cc Centering.cc ParticleIterator.cc
-  Koopman.cc BiorthBess.cc)
+  Koopman.cc BiorthBess.cc SGSmooth.cc)
 add_library(expui ${expui_SOURCES})
 set_target_properties(expui PROPERTIES OUTPUT_NAME expui)
 target_include_directories(expui PUBLIC ${common_INCLUDE})
diff --git a/expui/CoefContainer.H b/expui/CoefContainer.H
index 66a3c9951..0c8bf39d0 100644
--- a/expui/CoefContainer.H
+++ b/expui/CoefContainer.H
@@ -91,6 +91,7 @@ namespace MSSA
     void pack_sphere();
     void unpack_sphere();
     void restore_background_sphere();
+    //@}
 
     //@{
     //! Cylindrical coefficients
@@ -99,6 +100,18 @@ namespace MSSA
     void restore_background_cylinder();
     //@}
 
+    //@{
+    //! Spherical field coefficients
+    void pack_sphfld();
+    void unpack_sphfld();
+    //@}
+
+    //@{
+    //! Cylindrical coefficients
+    void pack_cylfld();
+    void unpack_cylfld();
+    //@}
+
     //@{
     //! Slab coefficients
     void pack_slab();
diff --git a/expui/CoefContainer.cc b/expui/CoefContainer.cc
index 03af71d92..d5deb13da 100644
--- a/expui/CoefContainer.cc
+++ b/expui/CoefContainer.cc
@@ -3,6 +3,7 @@
 #include <iomanip>
 #include <fstream>
 #include <cctype>
+#include <typeinfo>
 
 #include <yaml-cpp/yaml.h>
 
@@ -64,8 +65,14 @@ namespace MSSA
     else if (dynamic_cast<CoefClasses::TableData*>(coefs.get())) {
       unpack_table();
     }
+    else if (dynamic_cast<CoefClasses::SphFldCoefs*>(coefs.get())) {
+      unpack_sphfld();
+    }
+    else if (dynamic_cast<CoefClasses::CylFldCoefs*>(coefs.get())) {
+      unpack_cylfld();
+    }
     else {
-      throw std::runtime_error("CoefDB::unpack_channels(): can not reflect coefficient type");
+      throw std::runtime_error(std::string("CoefDB::unpack_channels(): can not reflect coefficient type=") + typeid(*coefs.get()).name());
     }
 
     return coefs;
@@ -79,8 +86,12 @@ namespace MSSA
       restore_background_cylinder();
     else if (dynamic_cast<CoefClasses::TableData*>(coefs.get()))
       { } // Do nothing
+    else if (dynamic_cast<CoefClasses::SphFldCoefs*>(coefs.get()))
+      { } // Do nothing
+    else if (dynamic_cast<CoefClasses::CylFldCoefs*>(coefs.get()))
+      { } // Do nothing
     else {
-      throw std::runtime_error("CoefDB::background(): can not reflect coefficient type");
+      throw std::runtime_error(std::string("CoefDB::background(): can not reflect coefficient type=") + typeid(*coefs.get()).name());
     }
   }
 
@@ -96,8 +107,12 @@ namespace MSSA
       pack_cube();
     else if (dynamic_cast<CoefClasses::TableData*>(coefs.get()))
       pack_table();
+    else if (dynamic_cast<CoefClasses::SphFldCoefs*>(coefs.get()))
+      pack_sphfld();
+    else if (dynamic_cast<CoefClasses::CylFldCoefs*>(coefs.get()))
+      pack_cylfld();
     else {
-      throw std::runtime_error("CoefDB::pack_channels(): can not reflect coefficient type");
+      throw std::runtime_error(std::string("CoefDB::pack_channels(): can not reflect coefficient type=") + typeid(*coefs.get()).name());
     }
   }
 
@@ -108,7 +123,8 @@ namespace MSSA
     times = cur->Times();
     complexKey = true;
 
-    auto cf = dynamic_cast<CoefClasses::CylStruct*>( cur->getCoefStruct(times[0]).get() );
+    auto cf = dynamic_cast<CoefClasses::CylStruct*>
+      ( cur->getCoefStruct(times[0]).get() );
 
     int mmax = cf->mmax;
     int nmax = cf->nmax;
@@ -165,7 +181,9 @@ namespace MSSA
     }
 
     for (int t=0; t<ntimes; t++) {
-      cf = dynamic_cast<CoefClasses::CylStruct*>( cur->getCoefStruct(times[t]).get() );
+      cf = dynamic_cast<CoefClasses::CylStruct*>
+	( cur->getCoefStruct(times[t]).get() );
+
       for (auto k : keys) {
 	if (k[2]==0)
 	  data[k][t] = (*cf->coefs)(k[0], k[1]).real();
@@ -185,7 +203,8 @@ namespace MSSA
   void CoefDB::unpack_cylinder()
   {
     for (int i=0; i<times.size(); i++) {
-      auto cf = dynamic_cast<CoefClasses::CylStruct*>( coefs->getCoefStruct(times[i]).get() );
+      auto cf = dynamic_cast<CoefClasses::CylStruct*>
+	( coefs->getCoefStruct(times[i]).get() );
       
       for (auto k : keys0) {
 	auto c = k, s = k;
@@ -201,6 +220,89 @@ namespace MSSA
     // END time loop
   }
 
+  void CoefDB::pack_cylfld()
+  {
+    auto cur = dynamic_cast<CoefClasses::CylFldCoefs*>(coefs.get());
+
+    times = cur->Times();
+    complexKey = true;
+
+    auto cf = dynamic_cast<CoefClasses::CylFldStruct*>
+      ( cur->getCoefStruct(times[0]).get() );
+
+    int nfld = cf->nfld;
+    int mmax = cf->mmax;
+    int nmax = cf->nmax;
+    int ntimes = times.size();
+    
+    // Promote desired keys into c/s pairs
+    //
+    keys.clear();
+    for (auto v : keys0) {
+      // Sanity check rank
+      //
+      if (v.size() != 3) {
+	std::ostringstream sout;
+	sout << "CoefDB::pack_cylfld: key vector should have rank 3; "
+	     << "found rank " << v.size() << " instead";
+	throw std::runtime_error(sout.str());
+      }
+      // Sanity check values
+      //
+      else if (v[0] >= 0 and v[0] <  nfld and
+	       v[1] >= 0 and v[1] <= mmax and
+	       v[2] >= 0 and v[2] <= nmax ) {
+	auto c = v, s = v;
+	c.push_back(0);
+	s.push_back(1);
+	keys.push_back(c);
+	if (v[0]) keys.push_back(s);
+      } else {
+	throw std::runtime_error("CoefDB::pack_cylfld: key is out of bounds");
+      }
+    }
+
+    bkeys.clear();		// No background fields
+
+    // Only pack the keys in the list
+    //
+    for (auto k : keys) {
+      data[k].resize(ntimes);
+    }
+
+    for (int t=0; t<ntimes; t++) {
+      cf = dynamic_cast<CoefClasses::CylFldStruct*>
+	( cur->getCoefStruct(times[t]).get() );
+
+      for (auto k : keys) {
+	if (k[3]==0)
+	  data[k][t] = (*cf->coefs)(k[0], k[1], k[2]).real();
+	else
+	  data[k][t] = (*cf->coefs)(k[0], k[1], k[3]).imag();
+      }
+    }
+  }
+
+  void CoefDB::unpack_cylfld()
+  {
+    for (int i=0; i<times.size(); i++) {
+      auto cf = dynamic_cast<CoefClasses::CylFldStruct*>
+	( coefs->getCoefStruct(times[i]).get() );
+      
+      for (auto k : keys0) {
+	auto c = k, s = k;
+	c.push_back(0); s.push_back(1);
+
+	int f = k[1], m = k[1], n = k[2];
+
+	if (m==0) (*cf->coefs)(f, m, n) = {data[c][i], 0.0};
+	else      (*cf->coefs)(f, m, n) = {data[c][i], data[s][i]};
+      }
+      // END key loop
+    }
+    // END time loop
+  }
+
   void CoefDB::pack_sphere()
   {
     auto cur = dynamic_cast<CoefClasses::SphCoefs*>(coefs.get());
@@ -208,7 +310,8 @@ namespace MSSA
     times = cur->Times();
     complexKey = true;
 
-    auto cf = dynamic_cast<CoefClasses::SphStruct*>( cur->getCoefStruct(times[0]).get() );
+    auto cf = dynamic_cast<CoefClasses::SphStruct*>
+      ( cur->getCoefStruct(times[0]).get() );
 
     int lmax   = cf->lmax;
     int nmax   = cf->nmax;
@@ -272,7 +375,9 @@ namespace MSSA
     auto I = [](const Key& k) { return k[0]*(k[0]+1)/2 + k[1]; };
 
     for (int t=0; t<ntimes; t++) {
-      cf = dynamic_cast<CoefClasses::SphStruct*>( cur->getCoefStruct(times[t]).get() );
+      cf = dynamic_cast<CoefClasses::SphStruct*>
+	( cur->getCoefStruct(times[t]).get() );
+
       for (auto k : keys)  {
 	auto c = (*cf->coefs)(I(k), k[2]);
 	data[k][t] = c.real();
@@ -293,7 +398,8 @@ namespace MSSA
 
     for (int i=0; i<times.size(); i++) {
 
-      auto cf = dynamic_cast<CoefClasses::SphStruct*>( coefs->getCoefStruct(times[i]).get() );
+      auto cf = dynamic_cast<CoefClasses::SphStruct*>
+	( coefs->getCoefStruct(times[i]).get() );
       
       for (auto k : keys0) {
 	auto c = k, s = k;
@@ -309,6 +415,99 @@ namespace MSSA
     }
     // END time loop
   }
+
+  void CoefDB::pack_sphfld()
+  {
+    auto cur = dynamic_cast<CoefClasses::SphFldCoefs*>(coefs.get());
+
+    times = cur->Times();
+    complexKey = true;
+
+    auto cf = dynamic_cast<CoefClasses::SphFldStruct*>
+      ( cur->getCoefStruct(times[0]).get() );
+
+    int nfld   = cf->nfld;
+    int lmax   = cf->lmax;
+    int nmax   = cf->nmax;
+    int ntimes = times.size();
+    
+    // Make extended key list
+    //
+    keys.clear();
+    for (auto k : keys0) {
+      // Sanity check rank
+      //
+      if (k.size() != 4) {
+	std::ostringstream sout;
+	sout << "CoefDB::pack_sphfld: key vector should have rank 4; "
+	     << "found rank " << k.size() << " instead";
+	throw std::runtime_error(sout.str());
+      }
+      // Sanity check values
+      //
+      else if (k[0] <  nfld and k[0] >= 0 and
+	       k[1] <= lmax and k[1] >= 0 and
+	       k[2] <= k[1] and k[2] >= 0 and
+	       k[3] <  nmax and k[3] >= 0 ) {
+	
+	auto v = k;
+	v.push_back(0);
+	keys.push_back(v);
+	data[v].resize(ntimes);
+
+	if (k[2]>0) {
+	  v[4] = 1;
+	  keys.push_back(v);
+	  data[v].resize(ntimes);
+	}
+      }
+      else {
+	throw std::runtime_error("CoefDB::pack_sphfld: key is out of bounds");
+      }
+    }
+
+    bkeys.clear();		// No background for field quantities
+
+    auto I = [](const Key& k) { return k[1]*(k[1]+1)/2 + k[2]; };
+
+    for (int t=0; t<ntimes; t++) {
+
+      cf = dynamic_cast<CoefClasses::SphFldStruct*>
+	( cur->getCoefStruct(times[t]).get() );
+
+      for (auto k : keys)  {
+	auto c = (*cf->coefs)(k[0], I(k), k[3]);
+	data[k][t] = c.real();
+	if (k[4]) data[k][t] = c.imag();
+      }
+    }
+  }
+
+  void CoefDB::unpack_sphfld()
+  {
+    auto I = [](const Key& k) { return k[1]*(k[1]+1)/2 + k[2]; };
+
+    for (int i=0; i<times.size(); i++) {
+
+      auto cf = dynamic_cast<CoefClasses::SphFldStruct*>
+	( coefs->getCoefStruct(times[i]).get() );
+      
+      for (auto k : keys0) {
+
+	auto c = k, s = k;
+
+	c.push_back(0);
+	s.push_back(1);
+
+	int f = k[0], l = k[1], m = k[2], n = k[3];
+
+	if (m==0) (*cf->coefs)(f, I(k), n) = {data[c][i], 0.0       };
+	else      (*cf->coefs)(f, I(k), n) = {data[c][i], data[s][i]};
+      }
+      // END key loop
+    }
+    // END time loop
+  }
   
   void CoefDB::pack_slab()
   {
@@ -317,7 +516,8 @@ namespace MSSA
     times = cur->Times();
     complexKey = true;
 
-    auto cf = dynamic_cast<CoefClasses::SlabStruct*>( cur->getCoefStruct(times[0]).get() );
+    auto cf = dynamic_cast<CoefClasses::SlabStruct*>
+      ( cur->getCoefStruct(times[0]).get() );
 
     int nmaxx   = cf->nmaxx;
     int nmaxy   = cf->nmaxy;
@@ -376,7 +576,9 @@ namespace MSSA
     }
 
     for (int t=0; t<ntimes; t++) {
-      cf = dynamic_cast<CoefClasses::SlabStruct*>( cur->getCoefStruct(times[t]).get() );
+      cf = dynamic_cast<CoefClasses::SlabStruct*>
+	( cur->getCoefStruct(times[t]).get() );
+
       for (auto k : keys)  {
 	auto c = (*cf->coefs)(k[0], k[1], k[2]);
 	if (k[3]) data[k][t] = c.imag();
@@ -398,7 +600,8 @@ namespace MSSA
     times = cur->Times();
     complexKey = true;
 
-    auto cf = dynamic_cast<CoefClasses::CubeStruct*>( cur->getCoefStruct(times[0]).get() );
+    auto cf = dynamic_cast<CoefClasses::CubeStruct*>
+      ( cur->getCoefStruct(times[0]).get() );
 
     int nmaxx   = cf->nmaxx;
     int nmaxy   = cf->nmaxy;
@@ -457,7 +660,9 @@ namespace MSSA
     }
 
     for (int t=0; t<ntimes; t++) {
-      cf = dynamic_cast<CoefClasses::CubeStruct*>( cur->getCoefStruct(times[t]).get() );
+      cf = dynamic_cast<CoefClasses::CubeStruct*>
+	( cur->getCoefStruct(times[t]).get() );
+
       for (auto k : keys)  {
 	auto c = (*cf->coefs)(k[0], k[1], k[2]);
 	if (k[3]) data[k][t] = c.imag();
@@ -476,7 +681,8 @@ namespace MSSA
   {
     for (int i=0; i<times.size(); i++) {
 
-      auto cf = dynamic_cast<CoefClasses::SlabStruct*>( coefs->getCoefStruct(times[i]).get() );
+      auto cf = dynamic_cast<CoefClasses::SlabStruct*>
+	( coefs->getCoefStruct(times[i]).get() );
       
       for (auto k : keys0) {
 	auto c = k, s = k;
@@ -494,7 +700,8 @@ namespace MSSA
   {
     for (int i=0; i<times.size(); i++) {
 
-      auto cf = dynamic_cast<CoefClasses::CubeStruct*>( coefs->getCoefStruct(times[i]).get() );
+      auto cf = dynamic_cast<CoefClasses::CubeStruct*>
+	( coefs->getCoefStruct(times[i]).get() );
       
       for (auto k : keys0) {
 	auto c = k, s = k;
@@ -515,7 +722,8 @@ namespace MSSA
     times = cur->Times();
     complexKey = false;
 
-    auto cf = dynamic_cast<CoefClasses::TblStruct*>( cur->getCoefStruct(times[0]).get() );
+    auto cf = dynamic_cast<CoefClasses::TblStruct*>
+      ( cur->getCoefStruct(times[0]).get() );
 
     int cols    = cf->cols;
     int ntimes  = times.size();
@@ -552,7 +760,9 @@ namespace MSSA
       for (unsigned c=0; c<cols; c++) {
 	Key key = {c};
 
-	cf = dynamic_cast<CoefClasses::TblStruct*>( cur->getCoefStruct(times[t]).get() );
+	cf = dynamic_cast<CoefClasses::TblStruct*>
+	  ( cur->getCoefStruct(times[t]).get() );
+
 	data[key][t] = (*cf->coefs)(c).real();
       }
     }
@@ -562,7 +772,8 @@ namespace MSSA
   {
     for (int i=0; i<times.size(); i++) {
 
-      auto cf = dynamic_cast<CoefClasses::TblStruct*>( coefs->getCoefStruct(times[i]).get() );
+      auto cf = dynamic_cast<CoefClasses::TblStruct*>
+	( coefs->getCoefStruct(times[i]).get() );
 
       int cols = cf->cols;
 
@@ -680,6 +891,7 @@ namespace MSSA
     auto I = [](const Key& k) { return k[0]*(k[0]+1)/2 + k[1]; };
 
     for (int t=0; t<times.size(); t++) {
+
       // Cast the coefficient class to spherical
       auto cf = dynamic_cast<CoefClasses::SphStruct*>
 	(cur->getCoefStruct(times[t]).get());
@@ -700,6 +912,7 @@ namespace MSSA
     auto cur = dynamic_cast<CoefClasses::CylCoefs*>(coefs.get());
 
     for (int t=0; t<times.size(); t++) {
+
       // Cast the coefficient class to cylindrical
       auto cf = dynamic_cast<CoefClasses::CylStruct*>
 	(cur->getCoefStruct(times[t]).get());
diff --git a/expui/FieldBasis.H b/expui/FieldBasis.H
index d405c3179..25ef5dbe3 100644
--- a/expui/FieldBasis.H
+++ b/expui/FieldBasis.H
@@ -53,9 +53,9 @@ namespace BasisClasses
   
     //@{
     //! Parameters
-    std::string model, filename;
+    std::string model, modelname;
     int lmax, nmax;
-    double rmin, rmax, ascl, scale, delta;
+    double rmin, rmax, ascl, rmapping, delta;
     //@}
 
     //@{
@@ -116,8 +116,8 @@ namespace BasisClasses
     //! Register phase-space functionoid
     void addPSFunction(PSFunction func, std::vector<std::string>& labels);
 
-    //! Prescaling factor
-    void set_scale(const double scl) { scale = scl; }
+    //! Coordinate mapping factor
+    void set_scale(const double scl) { rmapping = scl; }
     
     //! Zero out coefficients to prepare for a new expansion
     void reset_coefs(void);
diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc
index eecdcc520..32533400d 100644
--- a/expui/FieldBasis.cc
+++ b/expui/FieldBasis.cc
@@ -11,16 +11,23 @@
 #include <cfenv>
 #endif
 
-extern double Ylm01(int ll, int mm);
 extern double plgndr(int l, int m, double x);
 
+static double Ylm(int l, int m)
+{
+  m = abs(m);
+  return sqrt( (2.0*l+1)/(4.0*M_PI) ) *
+    exp(0.5*(lgamma(1.0+l-m) - lgamma(1.0+l+m)));
+}
+
+
 namespace BasisClasses
 {
   const std::set<std::string>
   FieldBasis::valid_keys = {
-    "filename",
+    "modelname",
     "dof",
-    "scale",
+    "rmapping",
     "rmin",
     "rmax",
     "ascl",
@@ -60,18 +67,18 @@ namespace BasisClasses
   //
   void FieldBasis::configure()
   {
-    nfld     = 2;		// Weight and density fields by default
-    lmax     = 4;
-    nmax     = 10;
-    rmin     = 1.0e-4;
-    rmax     = 2.0;
-    ascl     = 0.01;
-    delta    = 0.005;
-    scale    = 0.05;
-    dof      = 3;
-    model    = "file";
-    name     = "field";
-    filename = "SLGridSph.model";
+    nfld      = 2;		// Weight and density fields by default
+    lmax      = 4;
+    nmax      = 10;
+    rmin      = 1.0e-4;
+    rmax      = 2.0;
+    ascl      = 0.01;
+    delta     = 0.005;
+    rmapping  = 0.05;
+    dof       = 3;
+    model     = "file";
+    name      = "field";
+    modelname = "SLGridSph.model";
 
     initialize();
 
@@ -111,8 +118,8 @@ namespace BasisClasses
     //
     if (model == "file") {
       std::vector<double> r, d;
-      std::ifstream in(filename);
-      if (not in) throw std::runtime_error("Error opening file: " + filename);
+      std::ifstream in(modelname);
+      if (not in) throw std::runtime_error("Error opening file: " + modelname);
     
       std::string line;
       while (std::getline(in, line)) {
@@ -156,13 +163,30 @@ namespace BasisClasses
     
     // Generate the orthogonal function instance
     //
-    ortho = std::make_shared<OrthoFunction>(nmax, densfunc, rmin, rmax, scale, dof);
+    ortho = std::make_shared<OrthoFunction>(nmax, densfunc, rmin, rmax, rmapping, dof);
 
     // Initialize fieldlabels
     //
     fieldLabels.clear();
     fieldLabels.push_back("weight");
     fieldLabels.push_back("density");
+
+    // Debug
+    //
+    if (true) {
+      auto tst = ortho->testOrtho();
+      double worst = 0.0;
+      for (int i=0; i<tst.rows(); i++) {
+	for (int j=0; j<tst.rows(); j++) {
+	  if (i==j) worst = std::max<double>(worst, fabs(1.0 - tst(i, j)));
+	  else      worst = std::max<double>(worst, fabs(tst(i, j)));
+	}
+      }
+      if (myid==0) {
+	std::cout << "FieldBasis::orthoTest: worst=" << worst << std::endl;
+	ortho->dumpOrtho("fieldbasis_ortho.dump");
+      }
+    }
   }
 
   void FieldBasis::allocateStore() 
@@ -187,17 +211,17 @@ namespace BasisClasses
     // Assign values from YAML
     //
     try {
-      if (conf["filename"])	filename = conf["filename"].as<std::string>();
-      if (conf["nfld"    ])     nfld     = conf["nfld"    ].as<int>();
-      if (conf["lmax"    ])     lmax     = conf["lmax"    ].as<int>();
-      if (conf["nmax"    ])     nmax     = conf["nmax"    ].as<int>();
-      if (conf["dof"     ])     dof      = conf["dof"     ].as<int>();
-      if (conf["rmin"    ])     rmin     = conf["rmin"    ].as<double>();
-      if (conf["rmax"    ])     rmax     = conf["rmax"    ].as<double>();
-      if (conf["ascl"    ])     ascl     = conf["ascl"    ].as<double>();
-      if (conf["delta"   ])     delta    = conf["delta"   ].as<double>();
-      if (conf["scale"   ])     scale    = conf["scale"   ].as<double>();
-      if (conf["model"   ])     model    = conf["model"   ].as<std::string>();
+      if (conf["modelname"])	modelname = conf["modelname"].as<std::string>();
+      if (conf["model"    ])    model     = conf["model"    ].as<std::string>();
+      if (conf["nfld"     ])    nfld      = conf["nfld"     ].as<int>();
+      if (conf["lmax"     ])    lmax      = conf["lmax"     ].as<int>();
+      if (conf["nmax"     ])    nmax      = conf["nmax"     ].as<int>();
+      if (conf["dof"      ])    dof       = conf["dof"      ].as<int>();
+      if (conf["rmin"     ])    rmin      = conf["rmin"     ].as<double>();
+      if (conf["rmax"     ])    rmax      = conf["rmax"     ].as<double>();
+      if (conf["ascl"     ])    ascl      = conf["ascl"     ].as<double>();
+      if (conf["delta"    ])    delta     = conf["delta"    ].as<double>();
+      if (conf["rmapping" ])    rmapping  = conf["rmapping" ].as<double>();
     }
     catch (YAML::Exception & error) {
       if (myid==0) std::cout << "Error parsing parameters in FieldBasis: "
@@ -254,9 +278,33 @@ namespace BasisClasses
     if (dof==2) {
       auto p = dynamic_cast<CoefClasses::CylFldStruct*>(c.get());
       coefs[0] = p->coefs;
+      store[0] = p->store;
+
+      // Sanity test dimensions
+      if (nfld!=p->nfld || lmax!=p->mmax || nmax!=p->nmax) {
+	std::ostringstream serr;
+	serr << "FieldBasis::set_coefs: dimension error! "
+	     << " nfld [" << nfld << "!= " << p->nfld << "]"
+	     << " mmax [" << lmax << "!= " << p->mmax << "]"
+	     << " nmax [" << nmax << "!= " << p->nmax << "]";
+	throw std::runtime_error(serr.str());
+      }
+
     } else {
       auto p = dynamic_cast<CoefClasses::SphFldStruct*>(c.get());
       coefs[0] = p->coefs;
+      store[0] = p->store;
+
+      // Sanity test dimensions
+      if (nfld!=p->nfld || lmax!=p->lmax || nmax!=p->nmax) {
+	std::ostringstream serr;
+	serr << "FieldBasis::set_coefs: dimension error! "
+	     << " nfld [" << nfld << "!= " << p->nfld << "]"
+	     << " lmax [" << lmax << "!= " << p->lmax << "]"
+	     << " nmax [" << nmax << "!= " << p->nmax << "]";
+	throw std::runtime_error(serr.str());
+      }
+
     }
   }
 
@@ -265,6 +313,8 @@ namespace BasisClasses
 			      double u, double v, double w)
   {
     constexpr std::complex<double> I(0, 1);
+    constexpr double fac0 = 1.0/sqrt(4*M_PI);
+
     int tid = omp_get_thread_num();
     PS3 pos{x, y, z}, vel{u, v, w};
 
@@ -287,12 +337,16 @@ namespace BasisClasses
       
       auto p = (*ortho)(R);
     
-      (*coefs[tid])(0, 0, 0) += mass*p(0);
+      (*coefs[tid])(0, 0, 0) += mass*p(0)*fac0;
 
       for (int m=0; m<=lmax; m++) {
-	std::complex<double> P = std::exp(I*(phi*m));
+	
+	std::complex<double> P = std::exp(-I*(phi*m));
+	
 	for (int n=0; n<nmax; n++) {
+
 	  (*coefs[tid])(1, m, n) += mass*P*p(n);
+
 	  for (int k=0; k<vec.size(); k++)
 	    (*coefs[tid])(k+2, m, n) += mass*P*p(n)*vec[k];
 	}
@@ -304,15 +358,24 @@ namespace BasisClasses
       
       (*coefs[tid])(0, 0, 0) += mass*p(0);
 
-      for (int l=0, k=0; l<=lmax; l++) {
-	for (int m=0; m<=l; m++, k++) {
+      for (int l=0, lm=0; l<=lmax; l++) {
+
+	double s = 1.0;
+
+	for (int m=0; m<=l; m++, lm++) {
+	  
+	  // Spherical harmonic value
 	  std::complex<double> P =
-	    std::exp(I*(phi*m))*Ylm01(l, m)*plgndr(l, m, cth);
-	
+	    std::exp(-I*(phi*m))*Ylm(l, m)*plgndr(l, m, cth) * s;
+
+	  s *= -1.0;		// Flip sign for next evaluation
+
 	  for (int n=0; n<nmax; n++) {
-	    (*coefs[tid])(1, m, n) += mass*P*p(n);
+
+	    (*coefs[tid])(1, lm, n) += mass*P*p(n);
+
 	    for (int k=0; k<vec.size(); k++)
-	      (*coefs[tid])(k+2, m, n) += mass*P*p(n)*vec[k];
+	      (*coefs[tid])(k+2, lm, n) += mass*P*p(n)*vec[k];
 	  }
 	}
       }
@@ -387,21 +450,38 @@ namespace BasisClasses
       auto p = (*ortho)(r*sqrt(fabs(1.0 - costh*costh)));
       for (int i=0; i<nfld; i++) {
 	for (int m=0; m<=lmax; m++) {
+	  double fac = m>0 ? 2.0 : 1.0;
 	  for (int n=0; n<nmax; n++) {
-	    ret[i] += std::real((*coefs[0])(i, m, n)*exp(I*(phi*m)))*p[n];
+	    ret[i] += fac *
+	      std::real((*coefs[0])(i, m, n)*exp(I*(phi*m)))*p[n];
 	  }
 	}
       }
       return ret;
     } else {
+      static bool firstime = true;
+      if (firstime) {
+	int tid = omp_get_thread_num();
+	std::ostringstream file;
+	file << "field.bin." << tid;
+	std::ofstream fout(file.str());
+	const auto& d = coefs[0]->dimensions();
+	fout << "Dim size: " << d.size();
+	for (int i=0; i<d.size(); i++) fout << ", dim " << i << ": " << d[i];
+	fout << std::endl;
+	for (auto v : store[0]) fout << v << std::endl;
+	firstime = false;
+      }
+
       std::vector<double> ret(nfld, 0);
       auto p = (*ortho)(r);
       for (int i=0; i<nfld; i++) {
 	for (int l=0, L=0; l<=lmax; l++) {
 	  for (int m=0; m<=l; m++, L++) {
+	    double fac = m>0 ? 2.0 : 1.0;
 	    for (int n=0; n<nmax; n++) {
 	      ret[i] += std::real((*coefs[0])(i, L, n)*exp(I*(phi*m)))*p[n]*
-		Ylm01(l, m)*plgndr(l, m, costh);
+		Ylm(l, m)*plgndr(l, m, costh) * fac;
 	    }
 	  }
 	}
diff --git a/expui/Koopman.H b/expui/Koopman.H
index a84b31dc7..921ac4089 100644
--- a/expui/Koopman.H
+++ b/expui/Koopman.H
@@ -119,7 +119,7 @@ namespace MSSA
       return L;
     }
 
-    //! Return the EDMD modes
+    //! Return the EDMD modes, an approximation to the Koopman eigenfunctions
     Eigen::MatrixXcd getModes()
     {
       if (not computed) koopman_analysis();
diff --git a/expui/Koopman.cc b/expui/Koopman.cc
index b32f00b92..cc85cca1f 100644
--- a/expui/Koopman.cc
+++ b/expui/Koopman.cc
@@ -194,7 +194,7 @@ namespace MSSA {
       Eigen::VectorXcd B  = Phi.inverse() * xx;
       Eigen::MatrixXcd LL = I.asDiagonal();
 	
-      // Propate the solution with the operator
+      // Propagate the solution using the approximate Koopman operator
       //
       for (int i=0; i<numT; i++) {
 	Y.row(i) = (Phi*LL*B).real();
diff --git a/expui/SGSmooth.H b/expui/SGSmooth.H
new file mode 100644
index 000000000..6615bffa2
--- /dev/null
+++ b/expui/SGSmooth.H
@@ -0,0 +1,13 @@
+#ifndef __SGSMOOTH_HPP__
+#define __SGSMOOTH_HPP__
+
+#include <vector>
+
+//! savitzky golay smoothing.
+std::vector<double> sg_smooth(const std::vector<double> &v, const int w, const int deg);
+
+//! numerical derivative based on savitzky golay smoothing.
+std::vector<double> sg_derivative(const std::vector<double> &v, const int w,
+                                const int deg, const double h=1.0);
+
+#endif // __SGSMOOTH_HPP__
diff --git a/expui/SGSmooth.cc b/expui/SGSmooth.cc
new file mode 100644
index 000000000..1b5150b53
--- /dev/null
+++ b/expui/SGSmooth.cc
@@ -0,0 +1,549 @@
+//!
+// Sliding window signal processing (and linear algebra toolkit).
+//
+// supported operations:
+// <ul>
+// <li> Savitzky-Golay smoothing.
+// <li> computing a numerical derivative based of Savitzky-Golay smoothing.
+// <li> required linear algebra support for SG smoothing using STL based
+//      vector/matrix classes
+// </ul>
+//
+// \brief Linear Algebra "Toolkit".
+//
+// modified by Rob Patro, 2016
+// modified by MDW, 2024
+
+// system headers
+#include <iostream>
+#include <sstream>
+#include <cstddef>             // for size_t
+#include <cmath>               // for fabs
+#include <vector>
+
+//! default convergence
+static const double TINY_FLOAT = 1.0e-300;
+
+//! comfortable array of doubles
+using float_vect = std::vector<double>;
+//! comfortable array of ints;
+using int_vect = std::vector<int>;
+
+/*! matrix class.
+ *
+ * This is a matrix class derived from a vector of float_vects.  Note that
+ * the matrix elements indexed [row][column] with indices starting at 0 (c
+ * style). Also note that because of its design looping through rows should
+ * be faster than looping through columns.
+ *
+ * \brief two dimensional floating point array
+ */
+class float_mat : public std::vector<float_vect> {
+private:
+    //! disable the default constructor
+    explicit float_mat() {};
+    //! disable assignment operator until it is implemented.
+    float_mat &operator =(const float_mat &) { return *this; };
+public:
+    //! constructor with sizes
+    float_mat(const size_t rows, const size_t cols, const double def=0.0);
+    //! copy constructor for matrix
+    float_mat(const float_mat &m);
+    //! copy constructor for vector
+    float_mat(const float_vect &v);
+
+    //! use default destructor
+    // ~float_mat() {};
+
+    //! get size
+    size_t nr_rows(void) const { return size(); };
+    //! get size
+    size_t nr_cols(void) const { return front().size(); };
+};
+
+
+
+// constructor with sizes
+float_mat::float_mat(const size_t rows,const size_t cols,const double defval)
+        : std::vector<float_vect>(rows) {
+    int i;
+    for (i = 0; i < rows; ++i) {
+        (*this)[i].resize(cols, defval);
+    }
+    if ((rows < 1) || (cols < 1)) {
+        std::ostringstream msg;
+        msg << "cannot build matrix with " << rows << " rows and " << cols
+               << "columns";
+        throw std::runtime_error(msg.str());
+    }
+}
+
+// copy constructor for matrix
+float_mat::float_mat(const float_mat &m) : std::vector<float_vect>(m.size()) {
+
+    float_mat::iterator inew = begin();
+    float_mat::const_iterator iold = m.begin();
+    for (/* empty */; iold < m.end(); ++inew, ++iold) {
+        const size_t oldsz = iold->size();
+        inew->resize(oldsz);
+        const float_vect oldvec(*iold);
+        *inew = oldvec;
+    }
+}
+
+// copy constructor for vector
+float_mat::float_mat(const float_vect &v)
+        : std::vector<float_vect>(1) {
+
+    const size_t oldsz = v.size();
+    front().resize(oldsz);
+    front() = v;
+}
+
+//////////////////////
+// Helper functions //
+//////////////////////
+
+//! permute() orders the rows of A to match the integers in the index array.
+void permute(float_mat &A, int_vect &idx)
+{
+    int_vect i(idx.size());
+    int j,k;
+
+    for (j = 0; j < A.nr_rows(); ++j) {
+        i[j] = j;
+    }
+
+    // loop over permuted indices
+    for (j = 0; j < A.nr_rows(); ++j) {
+        if (i[j] != idx[j]) {
+
+            // search only the remaining indices
+            for (k = j+1; k < A.nr_rows(); ++k) {
+                if (i[k] ==idx[j]) {
+                    std::swap(A[j],A[k]); // swap the rows and
+                    i[k] = i[j];     // the elements of
+                    i[j] = idx[j];   // the ordered index.
+                    break; // next j
+                }
+            }
+        }
+    }
+}
+
+/*! \brief Implicit partial pivoting.
+ *
+ * The function looks for pivot element only in rows below the current
+ * element, A[idx[row]][column], then swaps that row with the current one in
+ * the index map. The algorithm is for implicit pivoting (i.e., the pivot is
+ * chosen as if the max coefficient in each row is set to 1) based on the
+ * scaling information in the vector scale. The map of swapped indices is
+ * recorded in swp. The return value is +1 or -1 depending on whether the
+ * number of row swaps was even or odd respectively. */
+static int partial_pivot(float_mat &A, const size_t row, const size_t col,
+                         float_vect &scale, int_vect &idx, double tol)
+{
+    if (tol <= 0.0)
+        tol = TINY_FLOAT;
+
+    int swapNum = 1;
+
+    // default pivot is the current position, [row,col]
+    int pivot = row;
+    double piv_elem = fabs(A[idx[row]][col]) * scale[idx[row]];
+
+    // loop over possible pivots below current
+    int j;
+    for (j = row + 1; j < A.nr_rows(); ++j) {
+
+        const double tmp = fabs(A[idx[j]][col]) * scale[idx[j]];
+
+        // if this elem is larger, then it becomes the pivot
+        if (tmp > piv_elem) {
+            pivot = j;
+            piv_elem = tmp;
+        }
+    }
+
+#if 0
+    if (piv_elem < tol) {
+        throw std::runtime_error("partial_pivot(): Zero pivot encountered.");
+#endif
+
+    if(pivot > row) {           // bring the pivot to the diagonal
+        j = idx[row];           // reorder swap array
+        idx[row] = idx[pivot];
+        idx[pivot] = j;
+        swapNum = -swapNum;     // keeping track of odd or even swap
+    }
+    return swapNum;
+}
+
+/*! \brief Perform backward substitution.
+ *
+ * Solves the system of equations A*b=a, ASSUMING that A is upper
+ * triangular. If diag==1, then the diagonal elements are additionally
+ * assumed to be 1.  Note that the lower triangular elements are never
+ * checked, so this function is valid to use after a LU-decomposition in
+ * place.  A is not modified, and the solution, b, is returned in a. */
+static void lu_backsubst(float_mat &A, float_mat &a, bool diag=false)
+{
+    int r,c,k;
+
+    for (r = (A.nr_rows() - 1); r >= 0; --r) {
+        for (c = (A.nr_cols() - 1); c > r; --c) {
+            for (k = 0; k < A.nr_cols(); ++k) {
+                a[r][k] -= A[r][c] * a[c][k];
+            }
+        }
+        if(!diag) {
+            for (k = 0; k < A.nr_cols(); ++k) {
+                a[r][k] /= A[r][r];
+            }
+        }
+    }
+}
+
+/*! \brief Perform forward substitution.
+ *
+ * Solves the system of equations A*b=a, ASSUMING that A is lower
+ * triangular. If diag==1, then the diagonal elements are additionally
+ * assumed to be 1.  Note that the upper triangular elements are never
+ * checked, so this function is valid to use after a LU-decomposition in
+ * place.  A is not modified, and the solution, b, is returned in a. */
+static void lu_forwsubst(float_mat &A, float_mat &a, bool diag=true)
+{
+    int r,k,c;
+    for (r = 0;r < A.nr_rows(); ++r) {
+        for(c = 0; c < r; ++c) {
+            for (k = 0; k < A.nr_cols(); ++k) {
+                a[r][k] -= A[r][c] * a[c][k];
+            }
+        }
+        if(!diag) {
+            for (k = 0; k < A.nr_cols(); ++k) {
+                a[r][k] /= A[r][r];
+            }
+        }
+    }
+}
+
+/*! \brief Performs LU factorization in place.
+ *
+ * This is Crout's algorithm (cf., Num. Rec. in C, Section 2.3).  The map of
+ * swapped indeces is recorded in idx. The return value is +1 or -1
+ * depending on whether the number of row swaps was even or odd
+ * respectively.  idx must be preinitialized to a valid set of indices
+ * (e.g., {1,2, ... ,A.nr_rows()}). */
+static int lu_factorize(float_mat &A, int_vect &idx, double tol=TINY_FLOAT)
+{
+    if ( tol <= 0.0)
+        tol = TINY_FLOAT;
+
+    if ((A.nr_rows() == 0) || (A.nr_rows() != A.nr_cols())) {
+        throw std::runtime_error("lu_factorize(): cannot handle empty "
+                                 "or nonsquare matrices.");
+    }
+
+    float_vect scale(A.nr_rows());  // implicit pivot scaling
+    int i,j;
+    for (i = 0; i < A.nr_rows(); ++i) {
+        double maxval = 0.0;
+        for (j = 0; j < A.nr_cols(); ++j) {
+            if (fabs(A[i][j]) > maxval)
+                maxval = fabs(A[i][j]);
+        }
+        if (maxval == 0.0) {
+            throw std::runtime_error("lu_factorize(): zero pivot found.");
+        }
+        scale[i] = 1.0 / maxval;
+    }
+
+    int swapNum = 1;
+    int c,r;
+    for (c = 0; c < A.nr_cols() ; ++c) {            // loop over columns
+        swapNum *= partial_pivot(A, c, c, scale, idx, tol); // bring pivot to diagonal
+        for(r = 0; r < A.nr_rows(); ++r) {      //  loop over rows
+            int lim = (r < c) ? r : c;
+            for (j = 0; j < lim; ++j) {
+                A[idx[r]][c] -= A[idx[r]][j] * A[idx[j]][c];
+            }
+            if (r > c)
+                A[idx[r]][c] /= A[idx[c]][c];
+        }
+    }
+    permute(A,idx);
+    return swapNum;
+}
+
+/*! \brief Solve a system of linear equations.
+ * Solves the inhomogeneous matrix problem with lu-decomposition. Note that
+ * inversion may be accomplished by setting a to the identity_matrix. */
+static float_mat lin_solve(const float_mat &A, const float_mat &a,
+                           double tol=TINY_FLOAT)
+{
+    float_mat B(A);
+    float_mat b(a);
+    int_vect idx(B.nr_rows());
+    int j;
+
+    for (j = 0; j < B.nr_rows(); ++j) {
+        idx[j] = j;  // init row swap label array
+    }
+    lu_factorize(B,idx,tol); // get the lu-decomp.
+    permute(b,idx);          // sort the inhomogeneity to match the lu-decomp
+    lu_forwsubst(B,b);       // solve the forward problem
+    lu_backsubst(B,b);       // solve the backward problem
+    return b;
+}
+
+///////////////////////
+// related functions //
+///////////////////////
+
+//! Returns the inverse of a matrix using LU-decomposition.
+static float_mat invert(const float_mat &A)
+{
+    const int n = A.size();
+    float_mat E(n, n, 0.0);
+    float_mat B(A);
+    int i;
+
+    for (i = 0; i < n; ++i) {
+        E[i][i] = 1.0;
+    }
+
+    return lin_solve(B, E);
+}
+
+//! returns the transposed matrix.
+static float_mat transpose(const float_mat &a)
+{
+    float_mat res(a.nr_cols(), a.nr_rows());
+    int i,j;
+
+    for (i = 0; i < a.nr_rows(); ++i) {
+        for (j = 0; j < a.nr_cols(); ++j) {
+            res[j][i] = a[i][j];
+        }
+    }
+    return res;
+}
+
+//! matrix multiplication.
+float_mat operator *(const float_mat &a, const float_mat &b)
+{
+    float_mat res(a.nr_rows(), b.nr_cols());
+    if (a.nr_cols() != b.nr_rows()) {
+        throw std::runtime_error("incompatible matrices in multiplication");
+    }
+
+    int i,j,k;
+
+    for (i = 0; i < a.nr_rows(); ++i) {
+        for (j = 0; j < b.nr_cols(); ++j) {
+            double sum(0.0);
+            for (k = 0; k < a.nr_cols(); ++k) {
+                sum += a[i][k] * b[k][j];
+            }
+            res[i][j] = sum;
+        }
+    }
+    return res;
+}
+
+
+//! calculate savitzky golay coefficients.
+static float_vect sg_coeff(const float_vect &b, const size_t deg)
+{
+    const size_t rows(b.size());
+    const size_t cols(deg + 1);
+    float_mat A(rows, cols);
+    float_vect res(rows);
+
+    // generate input matrix for least squares fit
+    int i,j;
+    for (i = 0; i < rows; ++i) {
+        for (j = 0; j < cols; ++j) {
+            A[i][j] = pow(double(i), double(j));
+        }
+    }
+
+    float_mat c(invert(transpose(A) * A) * (transpose(A) * transpose(b)));
+
+    for (i = 0; i < b.size(); ++i) {
+        res[i] = c[0][0];
+        for (j = 1; j <= deg; ++j) {
+            res[i] += c[j][0] * pow(double(i), double(j));
+        }
+    }
+    return res;
+}
+
+/*! \brief savitzky golay smoothing.
+ *
+ * This method means fitting a polynome of degree 'deg' to a sliding window
+ * of width 2w+1 throughout the data.  The needed coefficients are
+ * generated dynamically by doing a least squares fit on a "symmetric" unit
+ * vector of size 2w+1, e.g. for w=2 b=(0,0,1,0,0). evaluating the polynome
+ * yields the sg-coefficients.  at the border non symmectric vectors b are
+ * used. */
+float_vect sg_smooth(const float_vect &v, const int width, const int deg)
+{
+    float_vect res(v.size(), 0.0);
+    if ((width < 1) || (deg < 0) || (v.size() < (2 * width + 2))) {
+        throw std::runtime_error("sgsmooth: parameter error.");
+    }
+
+    const int window = 2 * width + 1;
+    const int endidx = v.size() - 1;
+
+    // do a regular sliding window average
+    int i,j;
+    if (deg == 0) {
+        // handle border cases first because we need different coefficients
+#if defined(_OPENMP)
+#pragma omp parallel for private(i,j) schedule(static)
+#endif
+        for (i = 0; i < width; ++i) {
+	    const double scale = 1.0/double(i+1);
+            const float_vect c1(width, scale);
+            for (j = 0; j <= i; ++j) {
+                res[i]          += c1[j] * v[j];
+                res[endidx - i] += c1[j] * v[endidx - j];
+            }
+        }
+
+        // now loop over rest of data. reusing the "symmetric" coefficients.
+	const double scale = 1.0/double(window);
+        const  float_vect c2(window, scale);
+#if defined(_OPENMP)
+#pragma omp parallel for private(i,j) schedule(static)
+#endif
+        for (i = 0; i <= (v.size() - window); ++i) {
+            for (j = 0; j < window; ++j) {
+                res[i + width] += c2[j] * v[i + j];
+            }
+        }
+        return res;
+    }
+
+    // handle border cases first because we need different coefficients
+#if defined(_OPENMP)
+#pragma omp parallel for private(i,j) schedule(static)
+#endif
+    for (i = 0; i < width; ++i) {
+        float_vect b1(window, 0.0);
+        b1[i] = 1.0;
+
+        const float_vect c1(sg_coeff(b1, deg));
+        for (j = 0; j < window; ++j) {
+            res[i]          += c1[j] * v[j];
+            res[endidx - i] += c1[j] * v[endidx - j];
+        }
+    }
+
+    // now loop over rest of data. reusing the "symmetric" coefficients.
+    float_vect b2(window, 0.0);
+    b2[width] = 1.0;
+    const float_vect c2(sg_coeff(b2, deg));
+
+#if defined(_OPENMP)
+#pragma omp parallel for private(i,j) schedule(static)
+#endif
+    for (i = 0; i <= (v.size() - window); ++i) {
+        for (j = 0; j < window; ++j) {
+            res[i + width] += c2[j] * v[i + j];
+        }
+    }
+    return res;
+}
+
+/*! least squares fit a polynome of degree 'deg' to data in 'b'.
+ *  then calculate the first derivative and return it. */
+static float_vect lsqr_fprime(const float_vect &b, const int deg)
+{
+    const int rows(b.size());
+    const int cols(deg + 1);
+    float_mat A(rows, cols);
+    float_vect res(rows);
+
+    // generate input matrix for least squares fit
+    int i,j;
+    for (i = 0; i < rows; ++i) {
+        for (j = 0; j < cols; ++j) {
+            A[i][j] = pow(double(i), double(j));
+        }
+    }
+
+    float_mat c(invert(transpose(A) * A) * (transpose(A) * transpose(b)));
+
+    for (i = 0; i < b.size(); ++i) {
+        res[i] = c[1][0];
+        for (j = 1; j < deg; ++j) {
+            res[i] += c[j + 1][0] * double(j+1)
+                * pow(double(i), double(j));
+        }
+    }
+    return res;
+}
+
+/*! \brief savitzky golay smoothed numerical derivative.
+ *
+ * This method means fitting a polynome of degree 'deg' to a sliding window
+ * of width 2w+1 throughout the data.
+ *
+ * In contrast to the sg_smooth function we do a brute force attempt by
+ * always fitting the data to a polynome of degree 'deg' and using the
+ * result. */
+float_vect sg_derivative(const float_vect &v, const int width,
+                         const int deg, const double h)
+{
+    float_vect res(v.size(), 0.0);
+    if ((width < 1) || (deg < 1) || (v.size() < (2 * width + 2))) {
+        throw std::runtime_error("sgsderiv: parameter error");
+    }
+
+    const int window = 2 * width + 1;
+
+    // handle border cases first because we do not repeat the fit
+    // lower part
+    float_vect b(window, 0.0);
+    int i,j;
+
+    for (i = 0; i < window; ++i) {
+        b[i] = v[i] / h;
+    }
+    const float_vect c(lsqr_fprime(b, deg));
+    for (j = 0; j <= width; ++j) {
+        res[j] = c[j];
+    }
+    // upper part. direction of fit is reversed
+    for (i = 0; i < window; ++i) {
+        b[i] = v[v.size() - 1 - i] / h;
+    }
+    const float_vect d(lsqr_fprime(b, deg));
+    for (i = 0; i <= width; ++i) {
+        res[v.size() - 1 - i] = -d[i];
+    }
+
+    // now loop over rest of data. wasting a lot of least squares calcs
+    // since we only use the middle value.
+#if defined(_OPENMP)
+#pragma omp parallel for private(i,j) schedule(static)
+#endif
+    for (i = 1; i < (v.size() - window); ++i) {
+        for (j = 0; j < window; ++j) {
+            b[j] = v[i + j] / h;
+        }
+        res[i + width] = lsqr_fprime(b, deg)[width];
+    }
+    return res;
+}
+
+// Local Variables:
+// mode: c++
+// c-basic-offset: 4
+// fill-column: 76
+// indent-tabs-mode: nil
+// End:
diff --git a/expui/expMSSA.H b/expui/expMSSA.H
index 4c2bfa9f7..a5cdb1606 100644
--- a/expui/expMSSA.H
+++ b/expui/expMSSA.H
@@ -50,6 +50,13 @@ namespace MSSA
     //! Right singular vectors
     Eigen::MatrixXd U;
 
+    //@{
+    //! Koopman modes
+    Eigen::VectorXcd L;
+    Eigen::MatrixXcd Phi;
+    int window;
+    //@}
+
     //! Parameters
     //@{
     bool flip, verbose, powerf;
@@ -279,6 +286,14 @@ namespace MSSA
       for (auto v : mean) ret.push_back(v.first);
       return ret;
     }
+
+    //! Estimate Koopman modes from the trajectory eigenvectors
+    std::tuple<Eigen::VectorXcd, Eigen::MatrixXcd>
+    getKoopmanModes(const double tol, int window, bool debug);
+
+    //! Return the reconstructed Koopman modes
+    std::map<std::string, CoefClasses::CoefsPtr> getReconstructedKoopman(int mode);
+
   };
 
 
diff --git a/expui/expMSSA.cc b/expui/expMSSA.cc
index 5adf38f57..7a19c6f6d 100644
--- a/expui/expMSSA.cc
+++ b/expui/expMSSA.cc
@@ -1536,13 +1536,283 @@ namespace MSSA {
 	reconstructed = true;
       }
 
-    } catch (HighFive::Exception& err) {
+    } catch (HighFive::Exception& err) {    // Number of channels
+    //
+    nkeys = mean.size();
+
+    // Make sure parameters are sane
+    //
+    if (numW<=0) numW = numT/2;
+    if (numW > numT/2) numW = numT/2;
+
+    numK = numT - numW + 1;
+
+
       std::cerr << "**** Error opening or reading H5 file ****" << std::endl;
       throw;
     }
 
   }
 
+  std::tuple<Eigen::VectorXcd, Eigen::MatrixXcd>
+  expMSSA::getKoopmanModes(double tol, int D, bool debug)
+  {
+    bool use_fullKh = true;	// Use the non-reduced computation of
+				// Koopman/eDMD
+    // Number of channels
+    //
+    nkeys = mean.size();
+
+    // Make sure parameters are sane
+    //
+    if (numW<=0) numW = numT/2;
+    if (numW > numT/2) numW = numT/2;
+
+    numK = numT - numW + 1;
+
+    Eigen::VectorXd S1;
+    Eigen::MatrixXd Y1;
+    Eigen::MatrixXd V1;
+    Eigen::MatrixXd VT1;
+    Eigen::MatrixXd VT2;
+
+    // Make a new trajetory matrix with smoothing
+    //
+    Y1.resize(numK, numW*nkeys + D*(nkeys-1));
+
+    size_t n=0, offset=0;
+    for (auto k : mean) {
+      for (int i=0; i<numK; i++) {
+	if (n > 0) {
+	  // Back blending
+	  for (int j=0; j<D; j++) {
+	    Y1(i, offset-j) += data[k.first][i + j] *
+	      static_cast<double>(D-j)/D;
+	  }
+	}
+	// Main series
+	for (int j=0; j<numW; j++) {
+	  if (i + j < numT) Y1(i, j + offset) = data[k.first][i + j];
+	}
+	// Forward blending
+	if (n<nkeys-1) {
+	  for (int j=0; j<D; j++) {
+	    Y1(i, j + numW + offset) = data[k.first][i + numW - j - 2] *
+	      static_cast<double>(D-j)/D;
+	  }
+	}
+      }
+      offset += numW + D;
+      n++;
+    }
+      
+    double Scale = Y1.norm();
+
+    // auto YY = Y1/Scale;
+    auto YY = Y1;
+
+    // Use one of the built-in Eigen3 algorithms
+    //
+    /*
+      if (params["Jacobi"]) {
+      // -->Using Jacobi
+      Eigen::JacobiSVD<Eigen::MatrixXd>
+      svd(YY, Eigen::ComputeThinU | Eigen::ComputeThinV);
+      S1 = svd.singularValues();
+      V1 = svd.matrixV();
+      } else if (params["BDCSVD"]) {
+    */
+    // -->Using BDC
+    Eigen::BDCSVD<Eigen::MatrixXd>
+      svd(YY, Eigen::ComputeFullU | Eigen::ComputeFullV);
+      // svd(YY, Eigen::ComputeThinU | Eigen::ComputeThinV);
+    S1 = svd.singularValues();
+    V1 = svd.matrixV();
+    /*
+      } else {
+      // -->Use Random approximation algorithm from Halko, Martinsson,
+      //    and Tropp
+      int srank = std::min<int>(YY.cols(), YY.rows());
+      RedSVD::RedSVD<Eigen::MatrixXd> svd(YY, srank);
+      S1 = svd.singularValues();
+      V1 = svd.matrixV();
+      }
+    */
+    
+    std::cout << "shape V1 = " << V1.rows() << " x "
+	      << V1.cols() << std::endl;
+    
+    std::cout << "shape Y1 = " << Y1.rows() << " x "
+	      << Y1.cols() << std::endl;
+      
+    int lags = V1.rows();
+    int rank = V1.cols();
+
+    std::ofstream out;
+    if (debug) out.open("debug.txt");
+
+    if (out) out << "rank=" << rank << " lags=" << lags << std::endl;
+
+    VT1.resize(rank, lags-1);
+    VT2.resize(rank, lags-1);
+
+    for (int j=0; j<rank; j++) {
+      std::vector<double> uu;
+      for (int i=0; i<lags-1; i++) {
+	VT1(j, i) = V1(i+0, j);
+	VT2(j, i) = V1(i+1, j);
+      }
+    }
+    
+    // Singular values
+    Eigen::VectorXd SS;
+
+    // Left & right singular vectors
+    Eigen::MatrixXd UU, VV;
+
+    if (out) {
+      out << std::endl << "Y shape: " << YY.rows()
+	  << " x " << YY.cols() << std::endl << std::endl
+	  << "VT1" << std::endl << VT1.transpose() << std::endl << std::endl
+	  << "VT2" << std::endl << VT2.transpose() << std::endl << std::endl;
+    }
+
+    // SVD
+    // Use one of the built-in Eigen3 algorithms
+    //
+    /*
+    if (params["Jacobi"]) {
+      // -->Using Jacobi
+      Eigen::JacobiSVD<Eigen::MatrixXd>
+	// svd(VT1, Eigen::ComputeThinU | Eigen::ComputeThinV);
+	svd(VT1, Eigen::ComputeFullU | Eigen::ComputeFullV);
+      SS = svd.singularValues();
+      UU = svd.matrixU();
+      VV = svd.matrixV();
+    } else if (params["BDCSVD"]) {
+    */
+    {
+      // -->Using BDC
+      Eigen::BDCSVD<Eigen::MatrixXd>
+	// svd(VT1, Eigen::ComputeThinU | Eigen::ComputeThinV);
+	svd(VT1, Eigen::ComputeFullU | Eigen::ComputeFullV);
+      SS = svd.singularValues();
+      UU = svd.matrixU();
+      VV = svd.matrixV();
+    }
+      /*
+    } else {
+      // -->Use Random approximation algorithm from Halko, Martinsson,
+      //    and Tropp
+      // RedSVD::RedSVD<Eigen::MatrixXd> svd(VT1, std::min<int>(rank, numK-1));
+      RedSVD::RedSVD<Eigen::MatrixXd> svd(VT1, std::max<int>(VT1.rows(), VT2.cols()));
+      SS = svd.singularValues();
+      UU = svd.matrixU();
+      VV = svd.matrixV();
+    }
+      */
+
+    if (out) out << "Singular values" << std::endl << SS << std::endl;
+
+    // Compute inverse
+    for (int i=0; i<SS.size(); i++) {
+      if (std::fabs(SS(i)) > tol) SS(i) = 1.0/SS(i);
+      else                        SS(i) = 0.0;
+    }
+
+    // Compute full Koopman operator
+    if (use_fullKh) {
+
+      Eigen::MatrixXd DD(VV.cols(), UU.cols());
+      DD.setZero();
+      for (int i=0; i<SS.size(); i++) DD(i, i) = SS(i);
+
+      Eigen::MatrixXd AT = VT2 * VV * DD * UU.transpose();
+
+      std::cout << "Size of AT: " << AT.rows() << " x " << AT.cols() << std::endl;
+
+      // Compute spectrum
+      Eigen::EigenSolver<Eigen::MatrixXd> es(AT);
+
+      L    = es.eigenvalues();
+      Phi  = es.eigenvectors();
+      
+      if (out) {
+	out << std::endl << "Eigenvalues"  << std::endl << L   << std::endl
+	    << std::endl << "Eigenvectors" << std::endl << Phi << std::endl;
+      }
+
+    } 
+    // Compute the reduced Koopman operator
+    else {
+
+      Eigen::MatrixXd AT = UU.transpose() * (VT2 * VV) * SS.asDiagonal();
+
+      // Compute spectrum
+      Eigen::EigenSolver<Eigen::MatrixXd> es(AT, true);
+
+      L      = es.eigenvalues();
+      auto W = es.eigenvectors();
+
+      // Compute the EDMD modes
+      //
+      Eigen::VectorXcd Linv(L);
+      for (int i=0; i<L.size(); i++) {
+	if (std::abs(Linv(i)) > tol) Linv(i) = 1.0/Linv(i);
+	else Linv(i) = 0.0;
+      }
+
+      Phi  = VT2 * VV * SS.asDiagonal() * W * Linv.asDiagonal();
+      
+      if (out) {
+	out << std::endl << "Eigenvalues"  << std::endl << L << std::endl
+	    << std::endl << "Eigenvectors" << std::endl << Phi << std::endl;
+      }
+    }
+
+    // Cache window size
+    //
+    window = D;
+
+    return {L, Phi};
+  }
+
+  std::map<std::string, CoefClasses::CoefsPtr>
+  expMSSA::getReconstructedKoopman(int mode)
+  {
+    // Copy the original map for return
+    //
+    auto newdata = data;
+
+    size_t n=0, offset=0;
+
+    for (auto u : mean) {
+      
+      double disp = totVar;
+      if (type == TrendType::totPow) disp = totPow;
+      if (disp==0.0) disp = var[u.first];
+      
+      std::complex phase = 1.0;
+      for (int i=0; i<numW; i++, phase *= L(mode))
+  	newdata[u.first][i] = std::real(Phi(mode, offset+i)*phase)*disp + u.second;
+
+      offset += numW + window;
+      n++;
+    }
+
+    // Copy to the working data
+    //
+    for (auto v : newdata) {
+      if (verbose) std::cout << "Updating for: " << v.first << std::endl;
+      coefDB.setData(v.first, v.second);
+    }
+
+    // Copies working data back to the coefficient structures
+    //
+    return coefDB.endUpdate();
+  }
+
+
   expMSSA::expMSSA(const mssaConfig& config, int nW, int nPC, const std::string flags) : numW(nW), npc(nPC), trajectory(true)
   {
     // Parse the YAML string
@@ -1564,8 +1834,8 @@ namespace MSSA {
     //
     if (params["Traj"]) trajectory = params["Traj"].as<bool>();
 
-    std::cout << "Trajectory is " << std::boolalpha << trajectory
-	      << std::endl;
+    // std::cout << "Trajectory is " << std::boolalpha << trajectory
+    // << std::endl;
 
     // Eigen OpenMP reporting
     //
diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc
index 23d579c0f..e69be02aa 100644
--- a/exputil/EmpCylSL.cc
+++ b/exputil/EmpCylSL.cc
@@ -216,6 +216,128 @@ EmpCylSL::EmpCylSL(int nmax, int lmax, int mmax, int nord,
   maxSNR       = 0.0;
 }
 
+template <typename U>
+U getH5(std::string name, HighFive::File& file)
+{
+  if (file.hasAttribute(name)) {
+    U v; 
+    HighFive::Attribute vv = file.getAttribute(name);
+    vv.read(v);
+    return v;
+  } else {
+    std::ostringstream sout;
+    sout << "EmpCylSL: could not find <" << name << ">";
+    throw std::runtime_error(sout.str());
+  }
+}
+
+EmpCylSL::EmpCylSL(int mlim, std::string cachename)
+{
+  // Use default name?
+  //
+  if (cachename.size()==0)
+    throw std::runtime_error("EmpCylSL: you must specify a cachename");
+
+  cachefile = cachename;
+
+  // Open and read the cache file to get the needed input parameters
+  //
+  
+  try {
+    // Silence the HDF5 error stack
+    //
+    HighFive::SilenceHDF5 quiet;
+    
+    // Open the hdf5 file
+    //
+    HighFive::File file(cachename, HighFive::File::ReadOnly);
+    
+    // For basis ID
+    std::string forceID("Cylinder"), geometry("cylinder");
+    std::string model = EmpModelLabs[mtype];
+      
+    // Version check
+    //
+    if (file.hasAttribute("Version")) {
+      auto ver = getH5<std::string>(std::string("Version"), file);
+      if (ver.compare(Version))
+	throw std::runtime_error("EmpCylSL: version mismatch");
+    } else {
+      throw std::runtime_error("EmpCylSL: outdated cache");
+    }
+
+    NMAX    = getH5<int>("nmaxfid", file);
+    MMAX    = getH5<int>("mmax",    file);
+    LMAX    = getH5<int>("lmaxfid", file);
+    NORDER  = getH5<int>("nmax",    file);
+    CMAPR   = getH5<int>("cmapr",   file);
+    CMAPZ   = getH5<int>("cmapz",   file);
+    MMIN    = 0;
+    MLIM    = std::min<int>(mlim, MMAX);
+    NMIN    = 0;
+    NLIM    = std::numeric_limits<int>::max();
+    Neven   = getH5<int>("neven",   file);
+    Nodd    = getH5<int>("nodd",    file);
+    ASCALE  = getH5<double>("ascl", file);
+    HSCALE  = getH5<double>("hscl", file);
+  }
+  catch (YAML::Exception& error) {
+    std::ostringstream sout;
+    sout << "EmpCylSL::getHeader: invalid cache file <" << cachefile << ">. ";
+    sout << "YAML error in getHeader: " << error.what() << std::endl;
+    throw GenericError(sout.str(), __FILE__, __LINE__, 1038, false);
+  }
+
+  // Set EvenOdd if values seem sane
+  //
+  EvenOdd = false;
+  if (Nodd>=0 and Nodd<=NORDER and Nodd+Neven==NORDER) {
+    EvenOdd  = true;
+  }
+    
+  pfac     = 1.0/sqrt(ASCALE);
+  ffac     = pfac/ASCALE;
+  dfac     = ffac/ASCALE;
+  
+  EVEN_M   = false;
+
+  // Check whether MPI is initialized
+  //
+  int flag;
+  MPI_Initialized(&flag);
+  if (flag) use_mpi = true;
+  else      use_mpi = false;
+  
+  // Enable MPI code in SLGridSph
+  //
+  if (use_mpi and numprocs>1) SLGridSph::mpi = 1;
+  
+  // Choose table dimension
+  //
+  MPItable = 4;
+    
+  // Initialize storage and values
+  //
+  coefs_made.resize(multistep+1);
+  std::fill(coefs_made.begin(), coefs_made.end(), false);
+  
+  eof_made   = false;
+  
+  sampT        = 1;
+  defSampT     = 1;
+  tk_type      = None;
+  
+  cylmass      = 0.0;
+  cylmass1     = std::vector<double>(nthrds);
+  cylmass_made = false;
+  
+  hallfile     = "";
+  minSNR       = std::numeric_limits<double>::max();
+  maxSNR       = 0.0;
+
+  read_cache();
+}
+
 
 void EmpCylSL::reset(int numr, int lmax, int mmax, int nord, 
 		     double ascale, double hscale, int nodd,
@@ -815,7 +937,7 @@ std::string compare_out(std::string str, U one, U two)
   return sout.str();
 }
 
-int EmpCylSL::cache_grid(int readwrite, string cachename)
+int EmpCylSL::cache_grid(int readwrite, std::string cachename)
 {
 
   // Option to reset cache file name
diff --git a/exputil/OrthoFunction.cc b/exputil/OrthoFunction.cc
index 7a1badc20..e0bf4f1f9 100644
--- a/exputil/OrthoFunction.cc
+++ b/exputil/OrthoFunction.cc
@@ -1,3 +1,4 @@
+#include <fstream>
 #include <OrthoFunction.H>
 
 OrthoFunction::OrthoFunction
@@ -93,6 +94,32 @@ Eigen::MatrixXd OrthoFunction::testOrtho()
   return ret;
 }
 
+void OrthoFunction::dumpOrtho(const std::string& filename)
+{
+  std::ofstream fout(filename);
+
+  if (fout) {
+    fout << "# OrthoFunction dump" << std::endl;
+
+    const int number = 1000;
+
+    for (int i=0; i<number; i++) {
+      // Coordinate transformation
+      //
+      double x = xmin + dx*i/(number-1);
+      double r = x_to_r(x);
+      double y = 2.0*r/scale; if (segment) y = x;
+
+      // Evaluate polynomial
+      //
+      Eigen::VectorXd p = poly_eval(y, nmax) * W(r);
+      fout << std::setw(16) << r;
+      for (int n=0; n<nmax; n++) fout << std::setw(16) << p(n);
+      fout << std::endl;
+    }
+  }
+}
+
 Eigen::VectorXd OrthoFunction::operator()(double r)
 {
   // Enforce bounds
diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc
index 0984fff00..7924551b7 100644
--- a/exputil/SLGridMP2.cc
+++ b/exputil/SLGridMP2.cc
@@ -168,6 +168,42 @@ SLGridSph::SLGridSph(std::shared_ptr<SphericalModelTable> mod,
 }
 
 
+SLGridSph::SLGridSph(std::string cachename)
+{
+  if (cachename.size()) sph_cache_name  = cachename;
+  else throw std::runtime_error("SLGridSph: you must specify a cachename");
+
+  tbdbg = false;
+
+  int LMAX, NMAX, NUMR, CMAP, DIVERGE=0;
+  double RMIN, RMAX, RMAP, DFAC=1.0;
+
+  try {
+    
+    auto node = getHeader(cachename);
+
+    LMAX  = node["lmax"].as<int>();
+    NMAX  = node["nmax"].as<int>();
+    NUMR  = node["numr"].as<int>();
+    CMAP  = node["cmap"].as<int>();
+    RMIN  = node["rmin"].as<double>();
+    RMAX  = node["rmax"].as<double>();
+    RMAP  = node["rmapping"].as<double>();
+
+    model_file_name = node["model"].as<std::string>();
+    model = SphModTblPtr(new SphericalModelTable(model_file_name, diverge, dfac));
+  }
+  catch (YAML::Exception& error) {
+    std::ostringstream sout;
+    sout << "SLGridMP2: error parsing parameters from getHeader: "
+	 << error.what();
+    throw GenericError(sout.str(), __FILE__, __LINE__, 1039, false);
+  }
+
+  initialize(LMAX, NMAX, NUMR, RMIN, RMAX, false, CMAP, RMAP);
+}
+
+
 std::map<std::string, std::string>
 SLGridSph::cacheInfo(const std::string& cachefile, bool verbose)
 {
diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H
index 30a4a6587..233039ef7 100644
--- a/include/EmpCylSL.H
+++ b/include/EmpCylSL.H
@@ -501,6 +501,9 @@ public:
 	   double ascale, double hscale, int Nodd,
 	   std::string cachename);
 
+  //! Construct from cache file
+  EmpCylSL(int mlim, const std::string cache);
+
   //! Destructor
   ~EmpCylSL(void);
 
diff --git a/include/OrthoFunction.H b/include/OrthoFunction.H
index b6f725ad6..6cc90ad5d 100644
--- a/include/OrthoFunction.H
+++ b/include/OrthoFunction.H
@@ -100,6 +100,9 @@ public:
 
   //! Test orthogonality
   Eigen::MatrixXd testOrtho();
+
+  //! Dump the orthogonal function table
+  void dumpOrtho(const std::string& filename);
 };
 
 
diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H
index 95a9bf6b7..f83d98fe9 100644
--- a/include/SLGridMP2.H
+++ b/include/SLGridMP2.H
@@ -116,6 +116,10 @@ public:
 	    std::string cachename=".slgrid_sph_cache",
 	    bool Verbose=false);
 
+  //! Constructor (from cache file)
+  SLGridSph(std::string cachename);
+
+
   //! Destructor
   virtual ~SLGridSph();
 
@@ -216,6 +220,11 @@ public:
   double getRmax() { return rmax; }
   //@}
 
+  //@{
+  //! Get expansion limits
+  int getLmax() { return lmax; }
+  int getNmax() { return nmax; }
+  //@}
 };
 
 
diff --git a/pyEXP/MSSAWrappers.cc b/pyEXP/MSSAWrappers.cc
index 9fdc5658f..aa5ef5c57 100644
--- a/pyEXP/MSSAWrappers.cc
+++ b/pyEXP/MSSAWrappers.cc
@@ -326,7 +326,6 @@ void MSSAtoolkitClasses(py::module &m) {
         dict({id: Coefs},...)
              reconstructed time series in the original coefficient form
 
-
         Notes
         -----
         The reconstructed data will overwrite the memory of the original coefficient 
@@ -543,6 +542,48 @@ void MSSAtoolkitClasses(py::module &m) {
 	)");
 
 
+  f.def("getKoopmanModes", &expMSSA::getKoopmanModes,
+	R"(
+        Compute the Koopman mode estimate from the right-singular vectors
+
+        Uses eDMD to estimate the modes
+
+        Parameters
+        ----------
+        tol : double
+            singular value truncation level
+        window: int
+            Smoothing between serialized channels (0 for no smoothing)
+        debug : bool
+            flag indicating whether to print debug information
+
+        Notes
+        -----
+        Use getReconstructedKoopman() to copy the reconstruction for a
+        particular mode back to the coefficient db
+
+        Returns
+        -------
+        tuple(numpy.ndarray, numpy.ndarray)
+            vector of eigenvalues and modes
+    )", py::arg("tol")=1.0e-12, py::arg("window")=0, py::arg("debug")=false);
+
+  f.def("getReconstructedKoopman", &expMSSA::getReconstructedKoopman,
+	R"(
+        Reconstruct the coefficients for a particular Koopman mode
+
+        Parameters
+        ----------
+        mode: int
+            The index of the mode to be reconstructed
+
+        Returns
+        -------
+        dict({id: Coefs},...)
+             reconstructed time series in the original coefficient form
+
+    )", py::arg("mode"));
+
   f.def("getRC", &expMSSA::getRC,
 	R"(
         Access the detrended reconstructed channel series by internal key
diff --git a/src/OutVel.H b/src/OutVel.H
index 4fc907543..f2704120d 100644
--- a/src/OutVel.H
+++ b/src/OutVel.H
@@ -21,7 +21,7 @@
 
 /** Dump velocity flow coefficients at each interval
     
-    @param filename of the coefficient file
+    @param modelname is the file specifying the density model
 
     @param name of the component to dump
   
@@ -33,7 +33,7 @@
 
     @param nmax is the maximum radial order
 
-    @param scale is the scale parameter for the density field
+    @param rmapping is the coordinate mapping parameter for the expansion
 
     @param rmin is the minimum radius for the density field
 
@@ -50,7 +50,7 @@ class OutVel : public Output
 
 private:
 
-  std::string filename, model, outfile;
+  std::string modelname, model, outfile;
   double prev = -std::numeric_limits<double>::max();
   Component *tcomp;
   CoefClasses::CoefsPtr coefs;
diff --git a/src/OutVel.cc b/src/OutVel.cc
index a10a732ce..cd877e4d3 100644
--- a/src/OutVel.cc
+++ b/src/OutVel.cc
@@ -9,12 +9,12 @@
 
 const std::set<std::string>
 OutVel::valid_keys = {
-  "filename",
+  "modelname",
   "nint",
   "nintsub",
   "name",
   "dof",
-  "scale",
+  "rmapping",
   "rmin",
   "rmax",
   "ascl",
@@ -97,7 +97,7 @@ void OutVel::initialize()
       model    = conf["model"].as<std::string>();
     else {
       std::string message = "OutVel: no model specified. Please specify "
-	"either 'file' with the model 'filename' or 'expon' for the\n"
+	"either 'file' with the model 'modelname' or 'expon' for the\n"
 	"exponential disk model (i.e. Laguerre polynomials)";
       throw std::runtime_error(message);
     }
@@ -119,7 +119,7 @@ void OutVel::initialize()
       throw std::runtime_error(message);
     }
 
-    if (conf["filename"]) filename = conf["filename"].as<std::string>();
+    if (conf["modelname"]) modelname = conf["modelname"].as<std::string>();
 
   }
   catch (YAML::Exception & error) {