# HG changeset patch
# User Jordi Gutiérrez Hermoso <jordigh@gmail.com>
# Date 1218172137 18000
# Node ID 5144dd3c546835ca4db549a64cdcf2ba626a0b20
# Parent  e6d2cd3b6e773d27a089037fb80d3a7703643c0c
Almost completed adding all Doxygen docstrings.

diff --git a/Doxyfile b/Doxyfile
--- a/Doxyfile
+++ b/Doxyfile
@@ -25,7 +25,7 @@
 # The PROJECT_NAME tag is a single word (or a sequence of words surrounded 
 # by quotes) that should identify the project.
 
-PROJECT_NAME           =  
+PROJECT_NAME           =  RBF-DDM
 
 # The PROJECT_NUMBER tag can be used to enter a project or revision number. 
 # This could be handy for archiving the generated documentation or 
@@ -222,7 +222,7 @@
 # func(std::string) {}). This also make the inheritance and collaboration 
 # diagrams that involve STL classes more complete and accurate.
 
-BUILTIN_STL_SUPPORT    = NO
+BUILTIN_STL_SUPPORT    = YES
 
 # If you use Microsoft's C++/CLI language, you should set this option to YES to
 # enable parsing support.
diff --git a/diff_op.cpp b/diff_op.cpp
--- a/diff_op.cpp
+++ b/diff_op.cpp
@@ -10,39 +10,45 @@
   
 
   //************** Differential operator functions *********************
-  diff_op::diff_op(){
-  }
-
-  double diff_op::operator()(const realfunc &f, const point &p) const{
+  double diff_op::operator()(const realfunc &f, const point &p) const
+  {
     return at(f,p);
   }
   
-  double linear_diff_op2::operator()(const realfunc &f, const point &p) const{
+  double linear_diff_op2::operator()(const realfunc &f, const point
+  &p) const
+  {
     return diff_op::operator()(f,p);
   }
 
 
   double bdry_diff_op::operator()(const realfunc &f, 
-				  const point &p) const{
+				  const point &p) const
+  {
     return at(f,p);
   }
   double bdry_diff_op::operator()(const realfunc &f, 
-				  const point &p, const vector &n) const{
+				  const point &p, const vector &n)
+				  const 
+  {
     return at(f,p,n);
   }
 
   double dirichlet_op::at(const realfunc &f, 
-			  const point &p) const{
+			  const point &p) const
+  {
     return f(p);
   }
   double dirichlet_op::at(const realfunc &f, const point &p, 
-			  const vector &n) const{
+			  const vector &n) const
+  {
     n.size();
     return f(p);
   }
 
   double neumann_op::at(const realfunc &f, const point &p, 
-			const vector &n) const{
+			const vector &n) const
+  {
     size_t dim = n.size();
     vector grad(dim);
     for(size_t i = 1; i <= dim; i++)
@@ -50,28 +56,34 @@
     return grad%n/norm(n);
   }
 
-  double Id_op::at(const realfunc &f, const point &p) const{
+  double Id_op::at(const realfunc &f, const point &p) const
+  {
     return f(p);
   }
 
-  double del1::at(const realfunc &f, const point &p) const{
+  double del1::at(const realfunc &f, const point &p) const
+  {
     return f.d(p,1);
   }
 
-  double del1::at(const realfunc &f, const point &p, size_t i) const{
+  double del1::at(const realfunc &f, const point &p, size_t i) const
+  {
     return f.d(p,i);
   }
 
-  double del2::at(const realfunc &f, const point &p) const{
+  double del2::at(const realfunc &f, const point &p) const
+  {
     return f.d2(p,1,1);
   }
 
   double del2::at(const realfunc &f, const point &p, 
-		  size_t i, size_t j) const{
+		  size_t i, size_t j) const
+  {
     return f.d2(p,i,j);
   }
 
-  double Laplacian::at(const realfunc &f, const point &p) const{
+  double Laplacian::at(const realfunc &f, const point &p) const
+  {
     size_t n = p.size();
     double result = 0;
     del2 D2;
diff --git a/func.cpp b/func.cpp
--- a/func.cpp
+++ b/func.cpp
@@ -76,7 +76,7 @@
   }
   double realfunc::at(const point& p) const{
     if(myfunc == 0){
-      no_init(__LINE__,__FILE__);
+      throw no_init(__LINE__,__FILE__);
     }
 
     return myfunc(p);
@@ -101,11 +101,12 @@
     return 0;
   }
 
-  void realfunc::no_init(int line, string file) const{
+  error_handling::badArgument
+  realfunc::no_init(int line, string file) const{
     error_handling::badArgument exc;
     exc.line = line;
     exc.file = file;
     exc.reason = "Did not assign a function pointer to a realfunc object.";
-    throw exc;
+    return exc;
   }
 }
diff --git a/include/bvp.hpp b/include/bvp.hpp
--- a/include/bvp.hpp
+++ b/include/bvp.hpp
@@ -34,7 +34,7 @@
     domain(size_t dimension, set<point> intr, 
 	   set<point> bdry, map<point, vector> ns);
 
-    /*!/brief Construct the domain from filenames.
+    /*!\brief Construct the domain from filenames.
      * 
      * Each file containing the desired points as a matrix structure
      * where the number of columns must be the dimensionality of the
@@ -212,6 +212,7 @@
 	const realfunc& f_in, 
 	const realfunc& g_in)      
       : BVP(O, L_in, B_in, f_in, g_in){};
+    //FIXME: figure out how to templatise this.
     linear_BVP2(shared_ptr<const domain> O, 
 	shared_ptr<const linear_diff_op2> L_in, 
 	shared_ptr<const bdry_diff_op> B_in, 
diff --git a/include/diff_op.hpp b/include/diff_op.hpp
--- a/include/diff_op.hpp
+++ b/include/diff_op.hpp
@@ -1,3 +1,9 @@
+/*! \file diff_op.hpp 
+ *
+ *  \brief Differential operators for boundary-value problems (BVP)
+ *  are declared here.
+ */
+
 #ifndef __DIFF_OP_HPP__
 #define __DIFF_OP_HPP__
 
@@ -28,51 +34,68 @@
    *
    */
 
-
-  //A general differential operator. All it does is evaluate the
-  //operator applied at a specified point for a real-valued function
-  //that take a vector argument. This is a pure abstract class.
+  /*! \brief A general differential operator. 
+   *
+   * All this class is evaluate the operator applied at a specified
+   * point for a real-valued function that take a vector
+   * argument. This is a pure abstract class.
+   */
   class diff_op{ 
   public:
-    diff_op();
+    ///Evaluate the function \f$f\f$ at the point \f$p\f$.
+    /** All this does is call the pure virtual at(realfunc, point)
+     * function below.
+     */
     double operator()(const realfunc &f, const point &p) const;
+    /** \brief Pure virtual function. 
+     *
+     * Actual implementation of the differential operator goes here.
+     */
     virtual double at(const realfunc &f, const point &p) const = 0;
+    ///Nothing to destroy.
     virtual ~diff_op(){};
   };
 
-  //A linear diff_op. Also pure abstract class for now.  
-  class linear_diff_op : virtual public diff_op{ 
-  public:
-    linear_diff_op(){};
-  };
+  /// A linear diff_op. Also pure abstract class for now.  
+  class linear_diff_op : virtual public diff_op{  };
 
-  //An at most 2nd-order differential operator. 
-  class diff_op2 : virtual public diff_op{
-  public:
-    diff_op2(){};
-  };
+  /// An at most 2nd-order differential operator. 
+  class diff_op2 : virtual public diff_op{ };
 
-  //The heat, wave, and Laplace's equation use this kind of
-  //operators: they're both linear and at most 2nd order. 
+  /*! \brief The heat, wave, and Laplace's equation use this kind of
+   * operators: they're both linear and at most 2nd order.
+   */
   class linear_diff_op2 : public linear_diff_op, public diff_op2{
   public:
     double operator()(const realfunc &f, const point &p) const;
   };
   
 
-  //An operator for specifying boundary conditions
+  /// An operator for specifying boundary conditions
   class bdry_diff_op : public diff_op{
   public:
+    
     double operator()(const realfunc &f, const point &p) const;
+    /// Only calls at(realfunc, point, vector) below.
     double operator()(const realfunc &f, const point &p, 
 		      const vector &n) const;
     virtual double at(const realfunc &f, 
 		      const point &p) const =0;
+    /*! \brief A boundary differential operator may need to know about
+     *  the normal vector at the boundary.
+     *
+     *  A boundary differential operator must provide evaluation at
+     *  the boundary, or explicitly ignore it. Robin and Neumann
+     *  boundary operators both make use of this normal vector.
+     *  \param f - The function to which the operator is applied.
+     *  \param p - The point at which the operator is evaluated
+     *  \param n - The normal vector for this point on the boundary.
+     */
     virtual double at(const realfunc &f, const point &p, 
 		      const vector &n) const =0;
   };
 
-  //Dirichlet boundary conditions
+  /// Dirichlet boundary conditions
   class dirichlet_op : public bdry_diff_op{
   public:
     double at(const realfunc &f, const point &p) const;
@@ -80,7 +103,7 @@
   };
 
 
-  //Neumann boundary conditions
+  /// Neumann boundary conditions
   class neumann_op : public bdry_diff_op{
   public:
     double at(const realfunc &f, const point &p, const vector &n) const;    
@@ -88,30 +111,36 @@
     double at(const realfunc &f, const point &p) const {return f(p);};
   };
 
-  //Identity operator
+  /// Identity operator
   class Id_op : public linear_diff_op2{
     double at(const realfunc &f, const point &p) const;
   }; 
 
-  //Partial wrt to some direction. By default (i.e with the overloaded
-  //at() function), wrt first coordinate. Else, last argument gives
-  //the direction to evaluate the partial.
+  /*! \brief Partial wrt to some direction. 
+   *
+   * By default (i.e with the overloaded at() function), wrt first
+   * coordinate. Else, last argument gives the direction to evaluate
+   * the partial.
+   */
   class del1 : public linear_diff_op2{
   public:
     double at(const realfunc &f, const point &p) const;
     double at(const realfunc &f, const point &p, size_t i) const;
   };
 
-  //Second partials wrt some direction. By default (i.e with the
-  //overloaded at() function), wrt first coordinate, twice. Else, last
-  //two arguments give the directions to evaluate the partials.
+  /*! \brief Second partials wrt some direction. 
+   *
+   * By default (i.e with the overloaded at() function), wrt first
+   * coordinate, twice. Else, last two arguments give the directions
+   * to evaluate the partials.
+   */
   class del2 : public linear_diff_op2{
   public:    
     double at(const realfunc &f, const point &p) const;
     double at(const realfunc &f, const point &p, size_t i, size_t j) const;
   };
 
-  //The Laplacian
+  ///The Laplacian
   class Laplacian : public linear_diff_op2{ 
   public:
     double at(const realfunc &f, const point &p) const;
diff --git a/include/func.hpp b/include/func.hpp
--- a/include/func.hpp
+++ b/include/func.hpp
@@ -1,35 +1,71 @@
+/*! \file func.hpp
+ * \brief Real-valued functions defined here
+ * 
+ * A general class for a real-valued functions from \f$\mathbb{R}^n\f$
+ * to \f$\mathbb{R}\f$ instead of using naked function pointers, meant
+ * to be a class that can be derived from. Wrappers for GSL-style
+ * function pointers are also provided.
+ */
 #ifndef __FUNC_HPP__
 #define __FUNC_HPP__
 
 #include "linalg.hpp"
+#include "error.hpp"
 #include <gsl/gsl_deriv.h>
 
 namespace bvp{
   using namespace linalg;
 
-  //A real-valued function from R^n to R.
+  /*! \brief  A real-valued function from \f$\mathbb{R}^n\f$ to
+   *  \f$\mathbb{R}\f$ .
+   *
+   * That is, a realfunc is a real-valued function that takes in
+   * linalg::point and returns doubles. This class is meant to be
+   * generic and any other real-valued function (e.g. a
+   * radial_basis_function) should derive from this one.
+   *
+   * Note that the base class provides derivatives using
+   * finite-difference approximations from the GSL. Derived classes
+   * are intended to provide their own exact derivatives instead of
+   * using the default finite-difference approximations.
+   * 
+   */
   class realfunc{
   public:
+    /// Create a useless realfunc, not assigned to anything.
     realfunc(); 
+    /// Create a realfunc from a function pointer.
     realfunc( double(*f)(const point&));
+    /// Destroys nothing, but declared virtual for derived classes.
     virtual ~realfunc(){};
 
+    /// Give a function pointer to this realfunc.
     void set_function_ptr(double (f_in)(const point &p));
     
+    /// Evaluate the function at point
     double operator()(const point& p) const;
+    /// Evaluate the function at point
     virtual double at(const point& p) const;
 
-    //Derivatives in directions k, k1, k2.
+    /// First derivative at point \f$x\f$ and in the \f$k\f$th direction.
     virtual double d(const point& x, size_t k) const; 
+    /// Second derivative (FIXME: actually implement this)
     virtual double d2(const point& x, size_t k1, size_t k2) const ; 
   protected:
+    /// Machine epsilon
     static double eps; 
+    /// Square root of epsilon
     static double sqrteps;
+    /// Cube root of epsilon
     static double root3eps;
+    /// Have already initialised epsilon above?
     static bool initialised ;
   private:
+    /// Pointer to a function that takes linalg::point and returns double
     double (*myfunc)(const point &p);
-    void no_init(int line, string file) const; //For throwing exceptions
+    /// Exception builder
+    error_handling::badArgument
+    no_init(int line, string file) const; 
   };
 
 
diff --git a/include/interpolator.hpp b/include/interpolator.hpp
--- a/include/interpolator.hpp
+++ b/include/interpolator.hpp
@@ -1,3 +1,10 @@
+/*!\file interpolator.hpp
+ *
+ * \brief An interpolator \f$ \displaystyle u(x) = \sum_{\xi \in \Xi}
+ * \lambda_\xi \varphi_\xi(x) \f$, templated by RBFs, and its helper
+ * classes are declared here.
+ * 
+ */
 #ifndef __INTERPOLATOR_HPP__
 #define __INTERPOLATOR_HPP__
 
@@ -15,6 +22,26 @@
 
   class interp_values;
 
+  /*! \brief An interpolator \f$ \displaystyle u(x) = \sum_{\xi \in
+   *  \Xi} \lambda_\xi \varphi_\xi(x) \f$, where the \f$ \varphi_\xi
+   *  \f$ are RBFs.
+   *
+   * An interpolator \f$ \displaystyle u(x) = \sum_{\xi \in \Xi}
+   * \lambda_\xi \varphi_\xi(x) \f$, where the \f$\varphi_\xi\f$ are
+   * RBFs centred at \f$\xi\f$, and \f$\Xi\f$ is a set of points on
+   * the interior and boundary of some domain \f$\Omega\f$.
+   *
+   * The interpolator can be used in various ways, initialised by
+   * either a BVP::linear_bvp2 or by interpolation data
+   * (i.e. std::map< linalg::point, double>). Once an interpolator is
+   * initialised, it becomes a bona fide bvp::realfunc, so that its
+   * evaluations and derivatives become available. For certain
+   * problems, it is convenient to be able to evaluate the
+   * interpolator or its derivatives at all points of the domain. It
+   * is also useful to be able to interpolate again on the same domain
+   * or with a slightly modified BVP. Interfaces to both of these
+   * functions are provided.
+   */
   template<typename RBF>
   class interpolator : public realfunc{
   public:
@@ -39,7 +66,11 @@
      *         at those points
      *
      *  Must provide domain information. The values of Xi must match
-     *  points on the given domain Omega.
+     *  points on the given domain Omega. The reason why the domain
+     *  might be relevant is in case it is important to specify which
+     *  points are on the interior and which are on the boundary, so
+     *  that it's possible to interpolate again by changing the values
+     *  on either the interior or the boundary.
      */
     interpolator(shared_ptr<domain> Omega, const map<point, double>& Xi);
     //@}
@@ -81,7 +112,7 @@
     interp_values operator()() const; 
     ///Full domain evaluation
     interp_values at() const; 
-    /// Full domain first erivative
+    /// Full domain first derivative
     interp_values d(size_t k) const; 
     /// Full domain second derivatives
     interp_values d2(size_t k1, size_t k2) const; 
@@ -154,11 +185,11 @@
     ///Is the interpolator ready for use?
     bool initted;
 
-    /** @name Exception throwers
+    /** @name Exception builders
      */
     //@{
-    void not_initted(int line, string file) const; 
-    void not_precomputed(int line, string file) const;
+    error_handling::badArgument not_initted(int line, string file) const; 
+    error_handling::badArgument not_precomputed(int line, string file) const;
     //@}
     
     ///Coefficients of the RBFs
@@ -235,17 +266,23 @@
     linalg::vector v;
     ///Two interp_values can be added iff these two values matches.
     size_t rbfs_hash;
-    ///Exception thrower
-    void different_rbfs(int line, string file) const;
+    ///Exception builder
+    badArgument different_rbfs(int line, string file) const;
 
     template<typename RBF> friend class interpolator;
-    friend interp_values operator*(const interp_values& v, const
-				   interp_values& w); 
+    
+    /// For comfortable syntax
+    friend interp_values operator*(double a, const interp_values& w); 
   };
 
   ///For comfortable syntax
+  //@{
   template <typename RBF>
-  interpolator<RBF> operator*(double a, const interpolator<RBF>& u);
+  interpolator<RBF> operator*(double a, const interpolator<RBF>& u)
+  {
+    return u*a;
+  }
+  //@}
 }
 
 #endif
diff --git a/include/linalg.hpp b/include/linalg.hpp
--- a/include/linalg.hpp
+++ b/include/linalg.hpp
@@ -61,7 +61,9 @@
     /// Number of columns.
     size_t cols() const;
 
-    /*! \brief Fortran-style parenthetical indexing (hence Octave-style too)
+    /*! \brief Fortran-style parenthetical indexing (hence
+     *         Octave-style too)
+     *
      *  Indices start from 1.
      *  @param i - Row number.
      *  @param j - Column number.
@@ -73,6 +75,7 @@
     const double& operator()(const size_t i, const size_t j) const;
 
     /*! \brief Indexing for vectorisation.
+     *
      *  The GSL provides limited vectorisation routines which can be
      *  useful for operating on blocks of matrices all at once.
      *  @param i - Row to be sliced.
@@ -84,6 +87,7 @@
     const vector_view operator()(const size_t  i, const slice &b) const;
 
     /*! \brief Indexing for vectorisation.
+     *
      *  The GSL provides limited vectorisation routines which can be
      *  useful for operating on blocks of matrices all at once.
      *  @param a - Slice range
@@ -118,23 +122,27 @@
      */
     matrix inv()   const;    
     /*! \brief Tranpose.
+     *
      * Returns a copy of the original matrix, transposed. The original
      * matrix is not modified.
      * \returns Transposed copy of matrix.
      */
     matrix T()     const;    
     /*! \brief Trace.
+     *
      * The trace of a square matrix.
      * \returns The trace.
      */
     double tr()    const;    
     /*! \brief Determinant.
+     *
      * The determinant computed via an LU factorisation
      * \returns The determinant.
      */
     double det()   const;    
 
     /*! \brief Solves Mv = w for v with LU factorisation. 
+     *
      * Given another vector, will return the solution vector v to the
      * equation Mv = w.
      * \param w - Another vector, must have the same size as number of
diff --git a/include/rbf.hpp b/include/rbf.hpp
--- a/include/rbf.hpp
+++ b/include/rbf.hpp
@@ -89,7 +89,7 @@
     virtual ~piecewise_smooth_rbf();
   };
   
-  /// C-infty RBFs
+  /// \f$C^\infty\f$ RBFs
   class c_infty_rbf : public radial_basis_function{
   public:
     c_infty_rbf();
@@ -108,7 +108,7 @@
 //Specific rbf's follow...
 namespace rbf{
 
-  /// r^n with n odd
+  /// \f$r^n\f$ with \f$n\f$ odd
   class piecewise_polynomial : public piecewise_smooth_rbf{
   public:
     static void set_n(size_t new_n);
@@ -136,7 +136,7 @@
 }
 
 namespace rbf{
-  /// r^n log(r) with n even
+  /// \f$r^n \log(r)\f$ with \f$n\f$ even
   class thin_plate_spline : public piecewise_smooth_rbf{
   public:
     static void set_n(size_t new_n);
@@ -160,7 +160,7 @@
 
 
 namespace rbf{
-  /// sqrt(1+(eps*r)^2) with eps > 0
+  /// \f$\sqrt{1+(\epsilon r)^2}\f$ with \f$\epsilon > 0\f$
   class multiquadric : public c_infty_rbf{
   public:
     multiquadric(){};
@@ -178,7 +178,7 @@
 }
 
 namespace rbf{
-  /// 1/sqrt(1 + (eps*r)^2) with eps > 0
+  /// \f$1/\sqrt{1 + (\epsilon r)^2}\f$ with \f$\epsilon > 0\f$
   class inverse_multiquadric : public c_infty_rbf{
   public:
     inverse_multiquadric(){};
@@ -196,7 +196,7 @@
 }
 
 namespace rbf{
-  /// 1/(1+(eps*r)^2) with eps > 0
+  /// \f$1/(1+(\epsilon r)^2)\f$ with \f$\epsilon > 0 \f$
   class inverse_quadratic : public c_infty_rbf{
   public:
     inverse_quadratic(){};
@@ -214,7 +214,7 @@
 }
 
 namespace rbf{
-  /// exp(- (eps*r)^2) with eps > 0.
+  /// \f$ e^{- (\epsilon r)^2}\f$ with \f$\epsilon > 0\f$.
   class gaussian : public c_infty_rbf{
   public:
     gaussian(){};
diff --git a/interpolator.cpp b/interpolator.cpp
--- a/interpolator.cpp
+++ b/interpolator.cpp
@@ -250,7 +250,7 @@
   void interpolator<RBF>::precompute_ev()
   {
     if(!initted){
-      not_initted(__LINE__, __FILE__);
+      throw not_initted(__LINE__, __FILE__);
     }
 
     if(!ev_precomp){
@@ -281,7 +281,7 @@
   void interpolator<RBF>::precompute_d1()
   {
     if(!initted){
-      not_initted(__LINE__, __FILE__);
+      throw not_initted(__LINE__, __FILE__);
     }
 
     //Have the RBFs been precomputed yet?
@@ -317,7 +317,7 @@
   void interpolator<RBF>::precompute_d2()
   {
     if(!initted)
-      not_initted(__LINE__, __FILE__);
+      throw not_initted(__LINE__, __FILE__);
     if(!d2_precomp){
       size_t dim = thebvp -> get_domain() -> get_dimension();
       //precompute all second derivatives
@@ -357,7 +357,7 @@
   double interpolator<RBF>::at(const point& p) const
   {
     if(!initted){
-      not_initted(__LINE__, __FILE__);
+      throw not_initted(__LINE__, __FILE__);
     }
     std::vector<size_t> alpha; //empty vector
 
@@ -377,7 +377,7 @@
   template<typename RBF>
   double interpolator<RBF>::d(const point& p, size_t k) const{
     if(!initted){
-      not_initted(__LINE__, __FILE__);
+      throw not_initted(__LINE__, __FILE__);
     }
     std::vector<size_t> alpha(k); alpha[k-1]++;
     
@@ -397,7 +397,7 @@
   template<typename RBF>
   double interpolator<RBF>::d2(const point& p, size_t k1, size_t k2) const{
     if(!initted){
-      not_initted(__LINE__, __FILE__);
+      throw not_initted(__LINE__, __FILE__);
     }
     std::vector<size_t> alpha(k1>k2?k1:k2); alpha[k1-1]++; alpha[k2-1]++;
     
@@ -417,14 +417,15 @@
   // ------------------- Whole domain evaluations ----------------
 
   template<typename RBF>
-  void interpolator<RBF>::not_precomputed(int line, string file) const
+  error_handling::badArgument
+  interpolator<RBF>::not_precomputed(int line, string file) const
   {
     badArgument exc;
     exc.reason = "Cannot do whole domain evaluations without "
       "precomputed RBFs.";
     exc.line = line;
     exc.file = file;
-    throw exc;
+    return exc;
   }
 
   template<typename RBF>
@@ -437,7 +438,7 @@
   interp_values interpolator<RBF>::at() const
   {
     if(!ev_precomp)
-      not_precomputed(__LINE__, __FILE__);
+      throw not_precomputed(__LINE__, __FILE__);
      
     std::vector<size_t> alpha;
     return interp_values(precomp_values_vec[alpha],rbfs_hash);
@@ -447,7 +448,7 @@
   interp_values interpolator<RBF>::d(size_t k) const
   {
     if(!d1_precomp)
-      not_precomputed(__LINE__, __FILE__);
+      throw not_precomputed(__LINE__, __FILE__);
     
     std::vector<size_t> alpha(k); alpha[k-1]++;
     return interp_values(precomp_values_vec[alpha],rbfs_hash);
@@ -457,7 +458,7 @@
   interp_values interpolator<RBF>::d2(size_t k1, size_t k2) const
   {
     if(!d2_precomp)
-      not_precomputed(__LINE__, __FILE__);
+      throw not_precomputed(__LINE__, __FILE__);
     
     std::vector<size_t> alpha(k1>k2?k1:k2); alpha[k1-1]++; alpha[k2-1]++;
     return interp_values(precomp_values_vec[alpha],rbfs_hash);
@@ -524,7 +525,7 @@
   template<typename RBF>
   void interpolator<RBF>::set_f(const realfunc& f){
     if(!initted){
-      not_initted(__LINE__, __FILE__);
+      throw not_initted(__LINE__, __FILE__);
     }
     thebvp -> set_f(f);
     computecoeffs();
@@ -534,7 +535,7 @@
   void interpolator<RBF>::set_g(const realfunc& g)
   {
     if(!initted){
-      not_initted(__LINE__, __FILE__);
+      throw not_initted(__LINE__, __FILE__);
     }
     thebvp -> set_g(g);
     computecoeffs();
@@ -544,7 +545,7 @@
   void interpolator<RBF>::set_f(const map<point, double>& f)
   {
     if(!initted){
-      not_initted(__LINE__, __FILE__);
+      throw not_initted(__LINE__, __FILE__);
     }
     thebvp -> set_f(f);
     computecoeffs();
@@ -554,21 +555,22 @@
   void interpolator<RBF>::set_g(const map<point, double>& g)
   {
     if(!initted){
-      not_initted(__LINE__, __FILE__);
+      throw not_initted(__LINE__, __FILE__);
     }
     thebvp -> set_g(g);
     computecoeffs();
   }
 
   template<typename RBF>
-  void interpolator<RBF>::not_initted(int line, string file) const
+  error_handling::badArgument
+  interpolator<RBF>::not_initted(int line, string file) const
   {
     badArgument exc;
     exc.reason = 
       "Interpolator can't interpolate without initialisation data.";
     exc.line = line;
     exc.file = file;
-    throw exc;
+    return exc;
   }
 
   template<typename RBF>
@@ -627,21 +629,21 @@
   }
   
   // ************************ interp_values stuff *****************
-  void interp_values::different_rbfs(int line, string file) const
+  badArgument interp_values::different_rbfs(int line, string file) const
   {
     badArgument exc;
     exc.reason = "Can't operate on interpolated values from "
       "incompatible interpolators";
     exc.line = line;
     exc.file = file;
-    throw exc;
+    return exc;
   }
 
   interp_values 
   interp_values::operator+(const interp_values& w) const
   {
     if(rbfs_hash != w.rbfs_hash)
-      different_rbfs(__LINE__, __FILE__);
+      throw different_rbfs(__LINE__, __FILE__);
     
     return interp_values(v + w.v,rbfs_hash);
   }
@@ -650,7 +652,7 @@
   interp_values::operator-(const interp_values& w) const
   {
     if(rbfs_hash != w.rbfs_hash)
-      different_rbfs(__LINE__, __FILE__);
+      throw different_rbfs(__LINE__, __FILE__);
     
     return interp_values(v - w.v,rbfs_hash);
   }
@@ -659,7 +661,7 @@
   interp_values::operator*(const interp_values& w) const
   {
     if(rbfs_hash != w.rbfs_hash)
-      different_rbfs(__LINE__, __FILE__);
+      throw different_rbfs(__LINE__, __FILE__);
     return interp_values(v * w.v,rbfs_hash);
   }
 
@@ -667,7 +669,7 @@
   interp_values::operator/(const interp_values& w) const
   {
     if(rbfs_hash != w.rbfs_hash)
-      different_rbfs(__LINE__, __FILE__);
+      throw different_rbfs(__LINE__, __FILE__);
     
     return interp_values(v / w.v,rbfs_hash);
   }
@@ -685,20 +687,9 @@
   }
 
   //*************** Non-member global friend functions ******************
-
-  template <typename RBF>
-  interpolator<RBF> operator*(double a, const interpolator<RBF>& u)
-  {
-    return u*a;
-  }  
-
-  interp_values operator*(const interp_values& v, const interp_values&
-  w)
-  {
-    if(v.rbfs_hash != w.rbfs_hash)
-      v.different_rbfs(__LINE__, __FILE__);
-    
-    return interp_values(v.v / w.v,v.rbfs_hash);
+  interp_values operator*(double a, const interp_values&  v)
+  { 
+    return interp_values(a*v.v, v.rbfs_hash);
   } 
 
   //Instantiations
@@ -709,4 +700,6 @@
   template class interpolator<inverse_multiquadric>;
   template class interpolator<inverse_quadratic>;
   template class interpolator<gaussian>;
+
+
 }
diff --git a/main.cpp b/main.cpp
--- a/main.cpp
+++ b/main.cpp
@@ -121,7 +121,7 @@
 		 size_t n, string which_interp);
 
 template<typename RBF>
-void bdry_iter(interpolator<RBF> &u, interpolator<RBF> &v, 
+void set_bdry_conds(interpolator<RBF> &u, interpolator<RBF> &v, 
 	       const boost::shared_ptr<domain> Omega);
 
 
@@ -202,14 +202,14 @@
       
       cout << "k1 done!" << endl;
 
-      bdry_iter(k1[0],k1[1],Omega);
+      set_bdry_conds(k1[0],k1[1],Omega);
            
       //k2
       f_u.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2));
       f_v.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2));
       f_h.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2));
  
-      bdry_iter(k2[0],k2[1],Omega);
+      set_bdry_conds(k2[0],k2[1],Omega);
 
       cout << "k2 done!" << endl;
 
@@ -218,7 +218,7 @@
       f_v.set_interps(u0+(k2[0]/2), v0+(k2[1]/2), h0+(k2[2]/2));
       f_h.set_interps(u0+(k2[0]/2), v0+(k2[1]/2), h0+(k2[2]/2));
       
-      bdry_iter(k3[0],k3[1],Omega);
+      set_bdry_conds(k3[0],k3[1],Omega);
       
       cout << "k3 done!" << endl;
 
@@ -227,7 +227,7 @@
       f_v.set_interps(u0+k3[0], v0+k3[1], h0+k3[2]);
       f_h.set_interps(u0+k3[0], v0+k3[1], h0+k3[2]);
       
-      bdry_iter(k4[0],k4[1],Omega);
+      set_bdry_conds(k4[0],k4[1],Omega);
 
       cout << "k4 done!" << endl;
      
@@ -237,7 +237,7 @@
       h1.set_f(h0 + (k1[2] + 2*k2[2] + 2*k3[2] + k4[2])/6);
      
       //Enforce boundary conditions iteratively
-      bdry_iter(u1,v1,Omega);
+      set_bdry_conds(u1,v1,Omega);
       
       u0 = u1;
       v0 = v1;
@@ -248,7 +248,7 @@
   }
   catch(error exc){
     utils::show_exception(exc);
-    return ;
+    return 1;
   }
 
 }
@@ -392,25 +392,20 @@
 }
 
 template<typename RBF>
-void bdry_iter(interpolator<RBF> &u, interpolator<RBF> &v, 
-	       const boost::shared_ptr<domain> Omega)
+void set_bdry_conds(interpolator<RBF> & u, interpolator<RBF> & v, 
+		    const boost::shared_ptr<domain> Omega)
 {
-  Gu_or_v<RBF_TYPE> gu(2, Omega), gv(1, Omega);
-  double err = 1;
-  do{    
-    linalg::vector u_old, u_new, v_old, v_new;
-    u_old = at_bdry(u,Omega);
-    gu.set_other(v);
-    u.set_g(gu);
-    u_new = at_bdry(u,Omega);
-
-    v_old = at_bdry(v,Omega);
-    gv.set_other(u);
-    v.set_g(gv);
-    v_new = at_bdry(v,Omega);
-    
-    double err_u = norm(u_new - u_old);
-    double err_v = norm(v_new - v_old);
-    err = (err_u > err_v? err_u : err_v);
-  }while(err > 1e-2);
+  size_t m = Omega -> get_boundary().size();
+  linalg::vector u_new(m), v_new(m);
+  {
+    set<point>::const_iterator I;
+    for(I = Omega -> get_boundary().begin(); 
+	I !=  Omega->get_boundary().end();
+	I++)
+    {
+      
+    }
+  }
+      
+	
 }