--- a/src/include/bvp.hpp
+++ b/src/include/bvp.hpp
@@ -18,14 +18,13 @@
 #include "func.hpp"
 
 
-/// Boundary value problems namespace
-namespace bvp{
+namespace kwantxi{
   using std::map;
   using std::set;
   using std::string;
-  using linalg::point;
-  using linalg::vector;
-  using linalg::matrix;
+  using kwantxi::point;
+  using kwantxi::vector;
+  using kwantxi::matrix;
   using boost::shared_ptr;
 
   /*! \brief A domain is a bunch of interior points, a bunch of boundary points,

--- a/src/include/ddm.hpp
+++ b/src/include/ddm.hpp
@@ -16,7 +16,7 @@
 #include "diff_op.hpp"
 #include "func.hpp"
 
-namespace bvp{
+namespace kwantxi{
   using std::map;
   using std::set;
   using boost::shared_ptr;

--- a/src/include/diff_op.hpp
+++ b/src/include/diff_op.hpp
@@ -13,8 +13,7 @@
 #include "linalg.hpp"
 #include "func.hpp"
 
-namespace bvp{
-  using namespace linalg;
+namespace kwantxi{
   /*
    *Class graph: (other leaves may also be present)
    *

--- a/src/include/error.hpp
+++ b/src/include/error.hpp
@@ -18,8 +18,7 @@
 #include <gsl/gsl_errno.h>
 #include <string>
 
-/// Exception classes in this namespace
-namespace error_handling{
+namespace kwantxi{
   //Structs to be thrown as exceptions
   using std::string;
 
@@ -83,7 +82,7 @@
 
 
 //A few more details about the structs we're throwing as exceptions.
-namespace error_handling{
+namespace kwantxi{
   using std::string;
 
   ///GSL_CONTINUE = -2,  /* iteration has not converged */

--- a/src/include/func.hpp
+++ b/src/include/func.hpp
@@ -17,8 +17,7 @@
 #include "error.hpp"
 #include <gsl/gsl_deriv.h>
 
-namespace bvp{
-  using namespace linalg;
+namespace kwantxi{
 
   /*! \brief  A real-valued function from \f$\mathbb{R}^n\f$ to
    *  \f$\mathbb{R}\f$ .

--- a/src/include/interp_values.hpp
+++ b/src/include/interp_values.hpp
@@ -4,9 +4,8 @@
 #include "linalg.hpp"
 #include <map>
 
-namespace bvp{
+namespace kwantxi{
   using std::map;
-  using namespace linalg;
 
   class bdry_values;
 
@@ -31,7 +30,7 @@
   private:
     template<typename RBF> friend class interpolator;
     ///Only constructor available for interpolator class
-    normals(size_t rbfs_hash_in, const map<point, linalg::vector>&
+    normals(size_t rbfs_hash_in, const map<point, kwantxi::vector>&
             normals_in, size_t n_in);
     ///Does nothing
     normals(){};
@@ -47,7 +46,7 @@
     size_t m;
 
     ///Store normals in a matrix
-    linalg::matrix N;
+    kwantxi::matrix N;
   };
 
 
@@ -88,7 +87,7 @@
 
   protected:
     /// Interpolator creates interp_values with this
-    interp_values(const linalg::vector& v_in, size_t rbfs_hash_in,
+    interp_values(const kwantxi::vector& v_in, size_t rbfs_hash_in,
                   size_t n_in, size_t m_in) 
       : v(v_in), rbfs_hash(rbfs_hash_in), n(n_in), m(m_in) {};
 
@@ -100,7 +99,7 @@
     badArgument different_rbfs(int line, string file) const;
 
     ///Interpolated values go here
-    linalg::vector v;
+    kwantxi::vector v;
     ///Two interp_values can be added iff these two values matches.
     size_t rbfs_hash;
     ///Number of interior points
@@ -123,7 +122,7 @@
   private:
     friend class normals;
     ///Same as for parent class
-    bdry_values(const linalg::vector& v_in, size_t rbfs_hash_in,
+    bdry_values(const kwantxi::vector& v_in, size_t rbfs_hash_in,
                 size_t n_in, size_t m_in) 
       : interp_values(v_in, rbfs_hash_in, n_in, m_in) {};
   };
@@ -137,7 +136,7 @@
     intr_values(const interp_values& v_in);
   private:
     ///Same as for parent class
-    intr_values(const linalg::vector& v_in, size_t rbfs_hash_in,
+    intr_values(const kwantxi::vector& v_in, size_t rbfs_hash_in,
                 size_t n_in, size_t m_in) 
       : interp_values(v_in, rbfs_hash_in, n_in, m_in) {};
   };

--- a/src/include/interpolator.hpp
+++ b/src/include/interpolator.hpp
@@ -21,7 +21,7 @@
 #include "diff_op.hpp"
 #include "interp_values.hpp"
 
-namespace bvp{
+namespace kwantxi{
   using std::map;
   using boost::shared_ptr;
 
@@ -184,7 +184,7 @@
      */
     normals get_normals() const;
     
-    /** \brief This will put into an ostream a linalg::matrix where
+    /** \brief This will put into an ostream a kwantxi::matrix where
      *   the first \f$n-1\f$ columns will be all the points where the
      *   interpolator is defined and the \f$n^{th}\f$ column is the
      *   values that the interpolator takes on that point.
@@ -219,12 +219,12 @@
     /** @name Exception builders
      */
     //@{
-    error_handling::badArgument not_initted(int line, string file) const; 
-    error_handling::badArgument not_precomputed(int line, string file) const;
+    kwantxi::badArgument not_initted(int line, string file) const; 
+    kwantxi::badArgument not_precomputed(int line, string file) const;
     //@}
     
     ///Coefficients of the RBFs
-    linalg::vector coeffs;
+    kwantxi::vector coeffs;
     ///The RBFS
     std::vector<RBF> rbfs;
     ///Their hash
@@ -264,10 +264,10 @@
     mutable map<std::vector<size_t>, map<point, double> >
     precomp_values;
     /// Same, but storing the vectors.
-    mutable map<std::vector<size_t>, linalg::vector >
+    mutable map<std::vector<size_t>, kwantxi::vector >
     precomp_values_vec;
     /// Convert vector of values to a map
-    map<point, double> valsvec2map(const linalg::vector& values) const;
+    map<point, double> valsvec2map(const  kwantxi::vector& values) const;
   };
 
 

--- a/src/include/linalg.hpp
+++ b/src/include/linalg.hpp
@@ -20,9 +20,7 @@
 
 #include <boost/operators.hpp>
 
-/// Linear algebra namespace
-namespace linalg{
-  using namespace error_handling;
+namespace kwantxi{
 
   class slice; 
   class vector;
@@ -396,7 +394,7 @@
   typedef vector point; //Useful alias.
 }
 
-namespace linalg{ //Non-member functions.
+namespace kwantxi{ //Non-member functions.
   
   /*! @name I/O
    */
@@ -438,7 +436,7 @@
 }
 
 //Inlined functions
-namespace linalg{
+namespace kwantxi{
   inline double& matrix::operator()(const size_t i, const size_t j){
     try{
       if(LUfactored)

--- a/src/include/rbf.hpp
+++ b/src/include/rbf.hpp
@@ -13,7 +13,7 @@
 #include "func.hpp"
 
 
-namespace error_handling{
+namespace kwantxi{
   ///An exception struct for RBFs when attempting to use the wrong dimension.
   struct badDimension : public badArgument{
     badDimension() {};
@@ -23,12 +23,10 @@
   };
 }
 
-/// Radial basis function namespace
-namespace rbf{
-  using namespace linalg;
+namespace kwantxi{
 
   ///Base abstract class.
-  class radial_basis_function : public bvp::realfunc{
+  class radial_basis_function : public kwantxi::realfunc{
   public:
     /// Default constructor, centred at zero by default. 
     radial_basis_function();
@@ -84,7 +82,7 @@
 }
 
 //Two important subclasses
-namespace rbf{
+namespace kwantxi{
   /// Piecewise smooth RBFs
   class piecewise_smooth_rbf : public radial_basis_function{
   public:
@@ -110,7 +108,7 @@
 
 
 //Specific rbf's follow...
-namespace rbf{
+namespace kwantxi{
 
   /// \f$r^n\f$ with \f$n\f$ odd
   class piecewise_polynomial : public piecewise_smooth_rbf{
@@ -139,7 +137,7 @@
   typedef piecewise_polynomial conical; 
 }
 
-namespace rbf{
+namespace kwantxi{
   /// \f$r^n \log(r)\f$ with \f$n\f$ even
   class thin_plate_spline : public piecewise_smooth_rbf{
   public:
@@ -163,7 +161,7 @@
 }
 
 
-namespace rbf{
+namespace kwantxi{
   /// \f$\sqrt{1+(\epsilon r)^2}\f$ with \f$\epsilon > 0\f$
   class multiquadric : public c_infty_rbf{
   public:
@@ -181,7 +179,7 @@
   };
 }
 
-namespace rbf{
+namespace kwantxi{
   /// \f$1/\sqrt{1 + (\epsilon r)^2}\f$ with \f$\epsilon > 0\f$
   class inverse_multiquadric : public c_infty_rbf{
   public:
@@ -199,7 +197,7 @@
   };
 }
 
-namespace rbf{
+namespace kwantxi{
   /// \f$1/(1+(\epsilon r)^2)\f$ with \f$\epsilon > 0 \f$
   class inverse_quadratic : public c_infty_rbf{
   public:
@@ -217,7 +215,7 @@
   };
 }
 
-namespace rbf{
+namespace kwantxi{
   /// \f$ e^{- (\epsilon r)^2}\f$ with \f$\epsilon > 0\f$.
   class gaussian : public c_infty_rbf{
   public:

--- a/src/include/utils.hpp
+++ b/src/include/utils.hpp
@@ -14,11 +14,7 @@
 #include <set>
 #include "linalg.hpp"
 
-/*! \namespace utils
- *  \brief A few helpful functions that didn't seem to be readily
- *  classified anywhere else.   
- */ 
-namespace utils{
+namespace kwantxi{
   ///Clears whitespace from front and back of string s.
   std::string trim(const std::string& s);
 

new file mode 100644
--- /dev/null
+++ b/src/vtkplot.cpp
@@ -0,0 +1,179 @@
+#include <vtkPoints.h>
+#include <vtkPointData.h>
+#include <vtkDelaunay2D.h>
+#include <vtkLookupTable.h>
+
+#include <vtkPolyData.h>
+#include <vtkPolyDataNormals.h>
+#include <vtkPolyDataMapper.h>
+
+#include <vtkWarpScalar.h>
+
+#include <vtkCamera.h>
+#include <vtkActor.h>
+#include <vtkRenderWindow.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkRenderer.h>
+
+#include <vtkInteractorStyleTerrain.h>
+
+#include <vtkSmartPointer.h>
+#include <vtkDoubleArray.h>
+
+#include <vtkGraphicsFactory.h>
+#include <vtkImagingFactory.h>
+#include <vtkWindowToImageFilter.h>
+#include <vtkPNGWriter.h>
+
+#include <string>
+using std::string;
+
+#include "include/linalg.h"
+
+typedef double (*Func3d)(double, double);
+
+void InitOffscreen(bool offscreen);
+vtkSmartPointer<vtkDelaunay2D> TriangulateTerrain(Func3d func);
+
+void PlotPoints(vtkSmartPointer<vtkDelaunay2D> delaunay, 
+                bool offscreen=true,
+                const string& filename="vtkscreenshot.png");
+
+double F(double x, double y)
+{
+  return atan(x*y*sin(y));
+}
+
+int main (int argc, char ** argv)
+{
+	auto delaunay = TriangulateTerrain(&F);	
+  PlotPoints(delaunay,true);
+	return 0;
+}
+
+void InitOffscreen(bool offscreen)
+{
+  // Graphics Factory
+  vtkSmartPointer<vtkGraphicsFactory> graphics_factory
+    = vtkGraphicsFactory::New();
+  graphics_factory->SetOffScreenOnlyMode( offscreen );
+  graphics_factory->SetUseMesaClasses( offscreen );
+
+  // Imaging Factory
+  vtkSmartPointer<vtkImagingFactory>  imaging_factory
+    = vtkImagingFactory::New();
+  imaging_factory->SetUseMesaClasses( offscreen );
+}
+
+vtkSmartPointer<vtkDelaunay2D> TriangulateTerrain(Func3d func)
+{
+
+	vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
+  vtkSmartPointer<vtkDoubleArray> scalars = vtkSmartPointer<vtkDoubleArray>::New();
+ 
+	unsigned int GridSize = 50;
+	for(double x = -3; x < 3; x += 6.0/GridSize)
+	{
+		for(double y = -3; y < 3; y += 6.0/GridSize)
+		{
+      double z = func(x,y);
+			vtkIdType idx = points->InsertNextPoint(x, y, z);
+      scalars -> InsertTuple(idx, &z);
+		}
+	}
+ 
+	//add the grid points to a polydata object
+	vtkSmartPointer<vtkPolyData> polydata = vtkPolyData::New();
+	polydata->SetPoints(points);
+  polydata->GetPointData() -> SetScalars(scalars);
+ 
+	//triangulate the grid points
+	vtkSmartPointer<vtkDelaunay2D> delaunay = vtkDelaunay2D::New();
+	delaunay->SetInput(polydata);
+	delaunay->Update();
+
+  return delaunay;
+}
+
+void PlotPoints(vtkSmartPointer<vtkDelaunay2D> delaunay, 
+                bool offscreen, 
+                const string& filename)
+{
+  InitOffscreen(offscreen);
+
+  //Normals for Gourad shading
+  vtkSmartPointer<vtkPolyDataNormals> normals = vtkPolyDataNormals::New();
+  normals -> SetInputConnection(delaunay ->GetOutputPort() );
+  normals -> SetFeatureAngle(60);
+
+  //Set the colours for the rendering
+  vtkSmartPointer<vtkLookupTable> lut = vtkLookupTable::New();
+  lut -> SetHueRange(0.66667, 0.0);
+  lut -> SetNumberOfColors(256);
+  lut -> SetRampToLinear();
+  lut -> Build();
+
+  
+	// map the contours to graphical primitives
+	vtkSmartPointer<vtkPolyDataMapper> contMapper = vtkPolyDataMapper::New();
+  contMapper->SetInput(normals -> GetOutput() );
+  contMapper->SetScalarRange(0,1);
+  contMapper->SetLookupTable(lut);
+  contMapper->ImmediateModeRenderingOn();
+
+ 
+  // create an actor for the contours
+	vtkSmartPointer<vtkActor> contActor = vtkActor::New();
+	contActor->SetMapper(contMapper);
+ 
+  // a renderer and render window
+	vtkSmartPointer<vtkRenderer> ren1 = vtkRenderer::New();
+	vtkSmartPointer<vtkRenderWindow> renWin = vtkRenderWindow::New();
+  if(offscreen)
+  {
+    renWin->SetOffScreenRendering(1);
+  }
+	renWin->AddRenderer(ren1);
+
+  // an interactor
+  vtkSmartPointer<vtkRenderWindowInteractor> iren 
+    = vtkRenderWindowInteractor::New();
+  iren->SetRenderWindow(renWin);  
+
+	// add the actors to the scene
+	ren1->SetBackground(0.0,0.0,0.0);
+  ren1->GetActiveCamera()->SetViewUp(0,0,1);
+  ren1->GetActiveCamera()->SetPosition(1, 1, 1);
+  ren1->GetActiveCamera()->SetParallelProjection(1);
+	ren1->AddActor(contActor);
+  ren1->ResetCamera();
+ 
+  // render an image (lights and cameras are created automatically)
+	renWin->Render();
+
+  if(offscreen)
+  {
+    // Write to file
+    vtkSmartPointer<vtkWindowToImageFilter> windowToImageFilter = 
+      vtkSmartPointer<vtkWindowToImageFilter>::New();
+    windowToImageFilter->SetInput(renWin);
+    windowToImageFilter->Update();
+ 
+    vtkSmartPointer<vtkPNGWriter> writer = vtkSmartPointer<vtkPNGWriter>::New();
+    writer->SetFileName(filename.c_str());
+    writer->SetInput(windowToImageFilter->GetOutput());
+    writer->Write();
+  }
+  else
+  {
+    //The style
+    vtkSmartPointer<vtkInteractorStyleTerrain> terrain_style 
+      = vtkInteractorStyleTerrain::New();
+ 
+    // Set the style
+    iren -> SetInteractorStyle(terrain_style);
+  
+    // begin mouse interaction
+    iren->Start(); 
+  }
+}