diff --git a/cspline.h b/cspline.h
index b4af29aff4bd2f94f85b5334d4187fb8e59f83d9..97ae0337152b6fe2635fe3523b1bcc85ee47ca7a 100644
--- a/cspline.h
+++ b/cspline.h
@@ -21,16 +21,16 @@
 #define WANT_STREAM
 #define WANT_MATH
 
-using namespace NEWMAT;
-using namespace std;
+
 ///////////////////////////////////////////////////////
 
 
 namespace MISCMATHS {
+
   class Cspline{
   public:
     Cspline(){}
-    Cspline(ColumnVector& pnodes,ColumnVector& pvals):
+    Cspline(NEWMAT::ColumnVector& pnodes,NEWMAT::ColumnVector& pvals):
       nodes(pnodes),
       vals(pvals),
       n(nodes.Nrows())
@@ -39,7 +39,7 @@ namespace MISCMATHS {
       fitted=true;
     }
 
-    Cspline(ColumnVector& pnodes, Matrix& pcoefs) :
+    Cspline(NEWMAT::ColumnVector& pnodes, NEWMAT::Matrix& pcoefs) :
       nodes(pnodes),
       coefs(pcoefs),
       n(nodes.Nrows())
@@ -49,23 +49,23 @@ namespace MISCMATHS {
       fitted=false;
     };
 
-    void set(ColumnVector& pnodes,ColumnVector& pvals);
-    void set(ColumnVector& pnodes, Matrix& pcoefs);
+    void set(NEWMAT::ColumnVector& pnodes,NEWMAT::ColumnVector& pvals);
+    void set(NEWMAT::ColumnVector& pnodes, NEWMAT::Matrix& pcoefs);
 
     void fit();
     float interpolate(float xx) const;
     float interpolate(float xx,int ind) const;
-    ColumnVector interpolate(const ColumnVector& x) const;
-    ColumnVector interpolate(const ColumnVector& x, const ColumnVector& indvec) const;
+    NEWMAT::ColumnVector interpolate(const NEWMAT::ColumnVector& x) const;
+    NEWMAT::ColumnVector interpolate(const NEWMAT::ColumnVector& x, const NEWMAT::ColumnVector& indvec) const;
 
   protected:
 
     bool fitted;
-    ColumnVector nodes;
-    ColumnVector vals;
-    Matrix coefs;
+    NEWMAT::ColumnVector nodes;
+    NEWMAT::ColumnVector vals;
+    NEWMAT::Matrix coefs;
     int n;
-    void diff(const ColumnVector& x, ColumnVector& dx );
+    void diff(const NEWMAT::ColumnVector& x, NEWMAT::ColumnVector& dx );
 
   };
 }
diff --git a/f2z.h b/f2z.h
index b6fa5e49093e2cfdb8d9737558e05618431e302f..c270b5dfd2b257c36ce816c817ed1406bfe2db86 100644
--- a/f2z.h
+++ b/f2z.h
@@ -14,9 +14,7 @@
 #include "armawrap/newmatap.h"
 #include "armawrap/newmatio.h"
 #include "base2z.h"
-//#include "miscmaths.h"
 
-using namespace NEWMAT;
 
 namespace MISCMATHS {
 
@@ -28,9 +26,9 @@ namespace MISCMATHS {
 
       float convert(float f, int d1, int d2);
 
-      static void ComputeFStats(const ColumnVector& p_fs, int p_dof1, int p_dof2, ColumnVector& p_zs);
-      static void ComputeFStats(const ColumnVector& p_fs, int p_dof1, const ColumnVector& p_dof2, ColumnVector& p_zs);
-      static void ComputeFStats(const ColumnVector& p_fs, const ColumnVector& p_dof1, const ColumnVector& p_dof2, ColumnVector& p_zs);
+      static void ComputeFStats(const NEWMAT::ColumnVector& p_fs, int p_dof1, int p_dof2, NEWMAT::ColumnVector& p_zs);
+      static void ComputeFStats(const NEWMAT::ColumnVector& p_fs, int p_dof1, const NEWMAT::ColumnVector& p_dof2, NEWMAT::ColumnVector& p_zs);
+      static void ComputeFStats(const NEWMAT::ColumnVector& p_fs, const NEWMAT::ColumnVector& p_dof1, const NEWMAT::ColumnVector& p_dof2, NEWMAT::ColumnVector& p_zs);
 
     private:
       F2z() : Base2z()
diff --git a/histogram.h b/histogram.h
index 8ddbd0ec5b829592659a5109377a0a2f174723cf..a109e8e6916878a4cba907b6b33dc4c53ec40f80 100644
--- a/histogram.h
+++ b/histogram.h
@@ -18,8 +18,6 @@
 #include "armawrap/newmatio.h"
 #include "miscmaths.h"
 
-using namespace NEWMAT;
-
 namespace MISCMATHS {
 
   class Histogram
@@ -33,35 +31,35 @@ namespace MISCMATHS {
 
       Histogram(const Histogram& in){ *this=in;}
 
-      Histogram(const ColumnVector& psourceData, int numBins)
+      Histogram(const NEWMAT::ColumnVector& psourceData, int numBins)
 	: sourceData(psourceData), calcRange(true), bins(numBins){}
 
-      Histogram(const ColumnVector& psourceData, float phistMin, float phistMax, int numBins)
+      Histogram(const NEWMAT::ColumnVector& psourceData, float phistMin, float phistMax, int numBins)
 	: sourceData(psourceData), calcRange(false), histMin(phistMin), histMax(phistMax), bins(numBins){}
 
-      void set(const ColumnVector& psourceData, int numBins) {
+      void set(const NEWMAT::ColumnVector& psourceData, int numBins) {
 	sourceData=psourceData; calcRange=true; bins=numBins;
       }
 
-      void set(const ColumnVector& psourceData, float phistMin, float phistMax, int numBins) {
+      void set(const NEWMAT::ColumnVector& psourceData, float phistMin, float phistMax, int numBins) {
 	sourceData=psourceData; calcRange=false; histMin=phistMin; histMax=phistMax; bins=numBins;
       }
 
       void generate();
-      void generate(ColumnVector exclusion_values);
-      ColumnVector generateCDF();
+      void generate(NEWMAT::ColumnVector exclusion_values);
+      NEWMAT::ColumnVector generateCDF();
 
       float getHistMin() const {return histMin;}
       float getHistMax() const {return histMax;}
       void setHistMax(float phistMax) {histMax = phistMax;}
       void setHistMin(float phistMin) {histMin = phistMin;}
-      void setexclusion(ColumnVector exclusion_values) {exclusion =exclusion_values;}
+      void setexclusion(NEWMAT::ColumnVector exclusion_values) {exclusion =exclusion_values;}
       void smooth();
 
       int integrateAll() {return sourceData.Nrows();}
 
-      const ColumnVector& getData() {return histogram;}
-      void setData(const ColumnVector& phist) { histogram=phist;}
+      const NEWMAT::ColumnVector& getData() {return histogram;}
+      void setData(const NEWMAT::ColumnVector& phist) { histogram=phist;}
 
       int integrateToInf(float value) const { return integrate(value, histMax); }
       int integrateFromInf(float value) const { return integrate(histMin, value); }
@@ -76,16 +74,16 @@ namespace MISCMATHS {
       float getPercentile(float perc);
 
       inline int getNumBins() const {return bins;}
-      inline ColumnVector getCDF() const {return CDF;}
-      inline ColumnVector getsourceData()const {return sourceData;}
+      inline NEWMAT::ColumnVector getCDF() const {return CDF;}
+      inline NEWMAT::ColumnVector getsourceData()const {return sourceData;}
     protected:
 
     private:
 
-      ColumnVector sourceData;
-      ColumnVector histogram;
-      ColumnVector exclusion;
-      ColumnVector CDF;
+      NEWMAT::ColumnVector sourceData;
+      NEWMAT::ColumnVector histogram;
+      NEWMAT::ColumnVector exclusion;
+      NEWMAT::ColumnVector CDF;
 
       bool calcRange;
 
@@ -107,7 +105,7 @@ namespace MISCMATHS {
       return (bin*(histMax-histMin))/(float)bins + histMin;
     }
 
-  inline ColumnVector Histogram::generateCDF()
+  inline NEWMAT::ColumnVector Histogram::generateCDF()
     {
 
 
diff --git a/kernel.h b/kernel.h
index 43e305f032c7f58e06da113dc3efed199a22ff44..131084e48a76b4b7cf803fae684a1c729e597431 100644
--- a/kernel.h
+++ b/kernel.h
@@ -17,8 +17,6 @@
 #include <cmath>
 #include "armawrap/newmat.h"
 
-using namespace NEWMAT;
-using namespace std;
 
 namespace MISCMATHS {
 
@@ -33,9 +31,9 @@ namespace MISCMATHS {
       int p_widthx;
       int p_widthy;
       int p_widthz;
-      ColumnVector p_kernelx;
-      ColumnVector p_kernely;
-      ColumnVector p_kernelz;
+      NEWMAT::ColumnVector p_kernelx;
+      NEWMAT::ColumnVector p_kernely;
+      NEWMAT::ColumnVector p_kernelz;
 
       // stop all forms of creation except the constructors below
       kernelstorage();
@@ -48,8 +46,8 @@ namespace MISCMATHS {
       float *storey;
       float *storez;
 
-      kernelstorage(const ColumnVector& kx, const ColumnVector& ky,
-		    const ColumnVector& kz, int wx, int wy, int wz)
+      kernelstorage(const NEWMAT::ColumnVector& kx, const NEWMAT::ColumnVector& ky,
+		    const NEWMAT::ColumnVector& kz, int wx, int wy, int wz)
         {
 	  p_kernelx = kx; p_kernely = ky; p_kernelz = kz;
 	  p_widthx = wx;  p_widthy = wy;  p_widthz = wz;
@@ -92,9 +90,9 @@ namespace MISCMATHS {
       int widthx() const { return p_widthx; }
       int widthy() const { return p_widthy; }
       int widthz() const { return p_widthz; }
-      const ColumnVector& kernelx() const { return p_kernelx; }
-      const ColumnVector& kernely() const { return p_kernely; }
-      const ColumnVector& kernelz() const { return p_kernelz; }
+      const NEWMAT::ColumnVector& kernelx() const { return p_kernelx; }
+      const NEWMAT::ColumnVector& kernely() const { return p_kernely; }
+      const NEWMAT::ColumnVector& kernelz() const { return p_kernelz; }
 
     };
 
@@ -132,8 +130,8 @@ namespace MISCMATHS {
       }
 
 
-      void setkernel(const ColumnVector& kx, const ColumnVector& ky,
-		     const ColumnVector& kz, int wx, int wy, int wz)
+      void setkernel(const NEWMAT::ColumnVector& kx, const NEWMAT::ColumnVector& ky,
+		     const NEWMAT::ColumnVector& kz, int wx, int wy, int wz)
       {
 	// see if already in list:
 	storedkernel = new kernelstorage(kx,ky,kz,wx,wy,wz);
@@ -158,21 +156,21 @@ namespace MISCMATHS {
 
   //////// Support functions /////////
 
-  float kernelval(float x, int w, const ColumnVector& kernel);
+  float kernelval(float x, int w, const NEWMAT::ColumnVector& kernel);
   float sincfn(float x);
   float hanning(float x, int w);
   float blackman(float x, int w);
   float rectangular(float x, int w);
-  ColumnVector sinckernel1D(const string& sincwindowtype, int w, int n);
-  kernel sinckernel(const string& sincwindowtype, int w, int nstore);
-  kernel sinckernel(const string& sincwindowtype,
+  NEWMAT::ColumnVector sinckernel1D(const std::string& sincwindowtype, int w, int n);
+  kernel sinckernel(const std::string& sincwindowtype, int w, int nstore);
+  kernel sinckernel(const std::string& sincwindowtype,
 		    int wx, int wy, int wz, int nstore);
-  float extrapolate_1d(const ColumnVector& data, const int index);
-  float interpolate_1d(const ColumnVector& data, const float index);
-  float kernelinterpolation_1d(const ColumnVector& data, float index, const ColumnVector& userkernel, int width);
-  float kernelinterpolation_1d(const ColumnVector& data, float index);
-  float kernelinterpolation_1d(RowVector data, float index);
-  float hermiteinterpolation_1d(const ColumnVector& data, int p1, int p4, float t);
+  float extrapolate_1d(const NEWMAT::ColumnVector& data, const int index);
+  float interpolate_1d(const NEWMAT::ColumnVector& data, const float index);
+  float kernelinterpolation_1d(const NEWMAT::ColumnVector& data, float index, const NEWMAT::ColumnVector& userkernel, int width);
+  float kernelinterpolation_1d(const NEWMAT::ColumnVector& data, float index);
+  float kernelinterpolation_1d(NEWMAT::RowVector data, float index);
+  float hermiteinterpolation_1d(const NEWMAT::ColumnVector& data, int p1, int p4, float t);
 }
 
 #endif
diff --git a/minimize.h b/minimize.h
index 8201c8e5e2640af4cf27e86a9d3a686e8d11e0b3..5039346cc84c12c7180ac80bc23ef9592485ee40 100644
--- a/minimize.h
+++ b/minimize.h
@@ -20,12 +20,6 @@
 #include "miscmaths.h"
 
 
-#define WANT_STREAM
-#define WANT_MATH
-
-using namespace MISCMATHS;
-using namespace NEWMAT;
-using namespace std;
 ///////////////////////////////////////////////////////
 
 //fminsearch.m
@@ -35,7 +29,7 @@ namespace MISCMATHS {
 class pair_comparer
 {
 public:
-  bool operator()(const pair<float,ColumnVector>& p1,const pair<float,ColumnVector>& p2) const
+  bool operator()(const pair<float,NEWMAT::ColumnVector>& p1,const pair<float,NEWMAT::ColumnVector>& p2) const
   {
     return p1.first < p2.first;
   }
@@ -44,35 +38,35 @@ public:
   class EvalFunction;
   class gEvalFunction;
 
-float diff1(const ColumnVector& x, const EvalFunction& func, int i,float h,int errorord=4);// finite diff derivative
+float diff1(const NEWMAT::ColumnVector& x, const EvalFunction& func, int i,float h,int errorord=4);// finite diff derivative
 
-float diff2(const ColumnVector& x, const EvalFunction& func, int i,float h,int errorord=4);// finite diff 2nd derivative
+float diff2(const NEWMAT::ColumnVector& x, const EvalFunction& func, int i,float h,int errorord=4);// finite diff 2nd derivative
 
-float diff2(const ColumnVector& x, const EvalFunction& func, int i,int j,float h,int errorord=4);// finite diff cross derivative
+float diff2(const NEWMAT::ColumnVector& x, const EvalFunction& func, int i,int j,float h,int errorord=4);// finite diff cross derivative
 
-ReturnMatrix gradient(const ColumnVector& x, const EvalFunction& func,float h,int errorord=4);// finite diff derivative vector
+NEWMAT::ReturnMatrix gradient(const NEWMAT::ColumnVector& x, const EvalFunction& func,float h,int errorord=4);// finite diff derivative vector
 
-ReturnMatrix hessian(const ColumnVector& x, const EvalFunction& func,float h,int errorord=4);// finite diff hessian
+NEWMAT::ReturnMatrix hessian(const NEWMAT::ColumnVector& x, const EvalFunction& func,float h,int errorord=4);// finite diff hessian
 
-void minsearch(ColumnVector& x, const EvalFunction& func, ColumnVector& paramstovary);
+void minsearch(NEWMAT::ColumnVector& x, const EvalFunction& func, NEWMAT::ColumnVector& paramstovary);
 
-void scg(ColumnVector& x, const gEvalFunction& func, ColumnVector& paramstovary, float tol = 0.0000001, float eps=1e-16, int niters=500);
+void scg(NEWMAT::ColumnVector& x, const gEvalFunction& func, NEWMAT::ColumnVector& paramstovary, float tol = 0.0000001, float eps=1e-16, int niters=500);
 
 class EvalFunction
 {//Function where gradient is not analytic (or you are too lazy to work it out) (required for fminsearch)
 public:
   EvalFunction(){}
-  virtual float evaluate(const ColumnVector& x) const = 0; //evaluate the function
+  virtual float evaluate(const NEWMAT::ColumnVector& x) const = 0; //evaluate the function
   virtual ~EvalFunction(){};
 
-  virtual void minimize(ColumnVector& x)
+  virtual void minimize(NEWMAT::ColumnVector& x)
   {
-    ColumnVector paramstovary(x.Nrows());
+    NEWMAT::ColumnVector paramstovary(x.Nrows());
     paramstovary = 1;
     minsearch(x,*this,paramstovary);
   }
 
-  virtual void minimize(ColumnVector& x, ColumnVector& paramstovary)
+  virtual void minimize(NEWMAT::ColumnVector& x, NEWMAT::ColumnVector& paramstovary)
   {
     minsearch(x,*this,paramstovary);
   }
@@ -88,17 +82,17 @@ public:
   gEvalFunction() : EvalFunction(){}
   // evaluate is inherited from EvalFunction
 
-  virtual ReturnMatrix g_evaluate(const ColumnVector& x) const = 0; //evaluate the gradient
+  virtual NEWMAT::ReturnMatrix g_evaluate(const NEWMAT::ColumnVector& x) const = 0; //evaluate the gradient
   virtual ~gEvalFunction(){};
 
-  virtual void minimize(ColumnVector& x)
+  virtual void minimize(NEWMAT::ColumnVector& x)
   {
-    ColumnVector paramstovary(x.Nrows());
+    NEWMAT::ColumnVector paramstovary(x.Nrows());
     paramstovary = 1;
     scg(x,*this,paramstovary);
   }
 
-  virtual void minimize(ColumnVector& x, ColumnVector& paramstovary)
+  virtual void minimize(NEWMAT::ColumnVector& x, NEWMAT::ColumnVector& paramstovary)
   {
     scg(x,*this,paramstovary);
   }
diff --git a/miscmaths.h b/miscmaths.h
index 11717b708442d7387637426125285950ef054338..97d1839585dc18d00acd64e3d3ded488054ab6d0 100644
--- a/miscmaths.h
+++ b/miscmaths.h
@@ -30,12 +30,9 @@
 
 //#pragma interface
 
-using namespace NEWMAT;
-using namespace NiftiIO;
-using namespace std;
-
 namespace MISCMATHS {
 
+
 #ifndef M_PI
 #define M_PI            3.14159265358979323846
 #endif
@@ -44,40 +41,40 @@ namespace MISCMATHS {
 #define LOGOUT(t) LogSingleton::getInstance().str()<<#t "="<<t<<endl;
 
   // IO/string stuff
-  template <class T> string num2str(T n, int width=-1);
+  template <class T> std::string num2str(T n, int width=-1);
 
 #if defined(_MSC_VER) && (_MSC_VER < 1300)
-  template <class T> string num2str(T n) { return num2str(n -1); }
+  template <class T> std::string num2str(T n) { return num2str(n -1); }
 #endif
 
-  string size(const Matrix& mat);
-  bool isNumber(const string& str);
-  ReturnMatrix read_ascii_matrix(const string& filename, int nrows, int ncols);
-  ReturnMatrix read_ascii_matrix(int nrows, int ncols, const string& filename);
-  ReturnMatrix read_ascii_matrix(const string& filename);
-  ReturnMatrix read_vest(string p_fname);
-  int read_binary_matrix(Matrix& mres, const string& filename);
-  ReturnMatrix read_binary_matrix(const string& filename);
-  ReturnMatrix read_matrix(const string& filename);
-
-  int write_ascii_matrix(const Matrix& mat, const string& filename,
+  std::string size(const NEWMAT::Matrix& mat);
+  bool isNumber(const std::string& str);
+  NEWMAT::ReturnMatrix read_ascii_matrix(const std::string& filename, int nrows, int ncols);
+  NEWMAT::ReturnMatrix read_ascii_matrix(int nrows, int ncols, const std::string& filename);
+  NEWMAT::ReturnMatrix read_ascii_matrix(const std::string& filename);
+  NEWMAT::ReturnMatrix read_vest(std::string p_fname);
+  int read_binary_matrix(NEWMAT::Matrix& mres, const std::string& filename);
+  NEWMAT::ReturnMatrix read_binary_matrix(const std::string& filename);
+  NEWMAT::ReturnMatrix read_matrix(const std::string& filename);
+
+  int write_ascii_matrix(const NEWMAT::Matrix& mat, const std::string& filename,
 			 int precision=-1);
-  int write_ascii_matrix(const string& filename, const Matrix& mat,
+  int write_ascii_matrix(const std::string& filename, const NEWMAT::Matrix& mat,
 			 int precision=-1);
-  int write_vest(const Matrix& x, string p_fname, int precision=-1);
-  int write_vest(string p_fname, const Matrix& x, int precision=-1);
-  int write_binary_matrix(const Matrix& mat, const string& filename);
+  int write_vest(const NEWMAT::Matrix& x, std::string p_fname, int precision=-1);
+  int write_vest(std::string p_fname, const NEWMAT::Matrix& x, int precision=-1);
+  int write_binary_matrix(const NEWMAT::Matrix& mat, const std::string& filename);
 
   // more basic IO calls
-  string skip_alpha(ifstream& fs);
-  ReturnMatrix read_ascii_matrix(ifstream& fs, int nrows, int ncols);
-  ReturnMatrix read_ascii_matrix(int nrows, int ncols, ifstream& fs);
-  ReturnMatrix read_ascii_matrix(ifstream& fs);
-  int read_binary_matrix(Matrix& mres, ifstream& fs);
-  ReturnMatrix read_binary_matrix(ifstream& fs);
-  int write_ascii_matrix(const Matrix& mat, ofstream& fs, int precision=-1);
-  int write_ascii_matrix(ofstream& fs, const Matrix& mat, int precision=-1);
-  int write_binary_matrix(const Matrix& mat, ofstream& fs);
+  std::string skip_alpha(std::ifstream& fs);
+  NEWMAT::ReturnMatrix read_ascii_matrix(std::ifstream& fs, int nrows, int ncols);
+  NEWMAT::ReturnMatrix read_ascii_matrix(int nrows, int ncols, std::ifstream& fs);
+  NEWMAT::ReturnMatrix read_ascii_matrix(std::ifstream& fs);
+  int read_binary_matrix(NEWMAT::Matrix& mres, std::ifstream& fs);
+  NEWMAT::ReturnMatrix read_binary_matrix(std::ifstream& fs);
+  int write_ascii_matrix(const NEWMAT::Matrix& mat, std::ofstream& fs, int precision=-1);
+  int write_ascii_matrix(std::ofstream& fs, const NEWMAT::Matrix& mat, int precision=-1);
+  int write_binary_matrix(const NEWMAT::Matrix& mat, std::ofstream& fs);
 
   // General maths
 
@@ -125,83 +122,83 @@ namespace MISCMATHS {
   template<class T>
    inline T Sqr(const T& x) { return x*x; }
 
-  ColumnVector cross(const ColumnVector& a, const ColumnVector& b);
-  ColumnVector cross(const double *a, const double *b);
+  NEWMAT::ColumnVector cross(const NEWMAT::ColumnVector& a, const NEWMAT::ColumnVector& b);
+  NEWMAT::ColumnVector cross(const double *a, const double *b);
 
-  inline float dot(const ColumnVector& a, const ColumnVector& b)
+  inline float dot(const NEWMAT::ColumnVector& a, const NEWMAT::ColumnVector& b)
     { return Sum(SP(a,b)); }
 
-  double norm2(const ColumnVector& x);
+  double norm2(const NEWMAT::ColumnVector& x);
   double norm2sq(double a, double b, double c);
   float norm2sq(float a, float b, float c);
 
-  ColumnVector seq(const int num);
+  NEWMAT::ColumnVector seq(const int num);
 
-  int diag(Matrix& m, const float diagvals[]);
-  int diag(Matrix& m, const ColumnVector& diagvals);
-  int diag(DiagonalMatrix& m, const ColumnVector& diagvals);
-  ReturnMatrix diag(const Matrix& mat);
-  ReturnMatrix pinv(const Matrix& mat);
-  int rank(const Matrix& X);
-  ReturnMatrix sqrtaff(const Matrix& mat);
+  int diag(NEWMAT::Matrix& m, const float diagvals[]);
+  int diag(NEWMAT::Matrix& m, const NEWMAT::ColumnVector& diagvals);
+  int diag(NEWMAT::DiagonalMatrix& m, const NEWMAT::ColumnVector& diagvals);
+  NEWMAT::ReturnMatrix diag(const NEWMAT::Matrix& mat);
+  NEWMAT::ReturnMatrix pinv(const NEWMAT::Matrix& mat);
+  int rank(const NEWMAT::Matrix& X);
+  NEWMAT::ReturnMatrix sqrtaff(const NEWMAT::Matrix& mat);
 
   // in the following mode = "new2old" or "old2new" (see .cc for more info)
-  vector<int> get_sortindex(const Matrix& vals, const string& mode, int col=1);
-  Matrix apply_sortindex(const Matrix& vals, vector<int> sidx, const string& mode);
-
-  void reshape(Matrix& r, const Matrix& m, int nrows, int ncols);
-  ReturnMatrix reshape(const Matrix& m, int nrows, int ncols);
-  int addrow(Matrix& m, int ncols);
-
-  int construct_rotmat_euler(const ColumnVector& params, int n, Matrix& aff);
-  int construct_rotmat_euler(const ColumnVector& params, int n, Matrix& aff,
-		   const ColumnVector& centre);
-  int construct_rotmat_quat(const ColumnVector& params, int n, Matrix& aff);
-  int construct_rotmat_quat(const ColumnVector& params, int n, Matrix& aff,
-		   const ColumnVector& centre);
-  int make_rot(const ColumnVector& angl, const ColumnVector& centre,
-	       Matrix& rot);
-
-  int getrotaxis(ColumnVector& axis, const Matrix& rotmat);
-  int rotmat2euler(ColumnVector& angles, const Matrix& rotmat);
-  int rotmat2quat(ColumnVector& quaternion, const Matrix& rotmat);
-  int decompose_aff(ColumnVector& params, const Matrix& affmat,
-		    int (*rotmat2params)(ColumnVector& , const Matrix& ));
-  int decompose_aff(ColumnVector& params, const Matrix& affmat,
-		    const ColumnVector& centre,
-		    int (*rotmat2params)(ColumnVector& , const Matrix& ));
-  int compose_aff(const ColumnVector& params, int n, const ColumnVector& centre,
-		  Matrix& aff,
-		  int (*params2rotmat)(const ColumnVector& , int , Matrix& ,
-				       const ColumnVector& ) );
-  float rms_deviation(const Matrix& affmat1, const Matrix& affmat2,
-		      const ColumnVector& centre, const float rmax);
-  float rms_deviation(const Matrix& affmat1, const Matrix& affmat2,
+  std::vector<int> get_sortindex(const NEWMAT::Matrix& vals, const std::string& mode, int col=1);
+  NEWMAT::Matrix apply_sortindex(const NEWMAT::Matrix& vals, std::vector<int> sidx, const std::string& mode);
+
+  void reshape(NEWMAT::Matrix& r, const NEWMAT::Matrix& m, int nrows, int ncols);
+  NEWMAT::ReturnMatrix reshape(const NEWMAT::Matrix& m, int nrows, int ncols);
+  int addrow(NEWMAT::Matrix& m, int ncols);
+
+  int construct_rotmat_euler(const NEWMAT::ColumnVector& params, int n, NEWMAT::Matrix& aff);
+  int construct_rotmat_euler(const NEWMAT::ColumnVector& params, int n, NEWMAT::Matrix& aff,
+		   const NEWMAT::ColumnVector& centre);
+  int construct_rotmat_quat(const NEWMAT::ColumnVector& params, int n, NEWMAT::Matrix& aff);
+  int construct_rotmat_quat(const NEWMAT::ColumnVector& params, int n, NEWMAT::Matrix& aff,
+		   const NEWMAT::ColumnVector& centre);
+  int make_rot(const NEWMAT::ColumnVector& angl, const NEWMAT::ColumnVector& centre,
+	       NEWMAT::Matrix& rot);
+
+  int getrotaxis(NEWMAT::ColumnVector& axis, const NEWMAT::Matrix& rotmat);
+  int rotmat2euler(NEWMAT::ColumnVector& angles, const NEWMAT::Matrix& rotmat);
+  int rotmat2quat(NEWMAT::ColumnVector& quaternion, const NEWMAT::Matrix& rotmat);
+  int decompose_aff(NEWMAT::ColumnVector& params, const NEWMAT::Matrix& affmat,
+		    int (*rotmat2params)(NEWMAT::ColumnVector& , const NEWMAT::Matrix& ));
+  int decompose_aff(NEWMAT::ColumnVector& params, const NEWMAT::Matrix& affmat,
+		    const NEWMAT::ColumnVector& centre,
+                    int (*rotmat2params)(NEWMAT::ColumnVector& , const NEWMAT::Matrix& ));
+  int compose_aff(const NEWMAT::ColumnVector& params, int n, const NEWMAT::ColumnVector& centre,
+		  NEWMAT::Matrix& aff,
+		  int (*params2rotmat)(const NEWMAT::ColumnVector& , int , NEWMAT::Matrix& ,
+				       const NEWMAT::ColumnVector& ) );
+  float rms_deviation(const NEWMAT::Matrix& affmat1, const NEWMAT::Matrix& affmat2,
+		      const NEWMAT::ColumnVector& centre, const float rmax);
+  float rms_deviation(const NEWMAT::Matrix& affmat1, const NEWMAT::Matrix& affmat2,
 		      const float rmax=80.0);
 
-  void get_axis_orientations(const Matrix& sform_mat, int sform_code,
-			     const Matrix& qform_mat, int qform_code,
+  void get_axis_orientations(const NEWMAT::Matrix& sform_mat, int sform_code,
+			     const NEWMAT::Matrix& qform_mat, int qform_code,
 			     int& icode, int& jcode, int& kcode);
 
   // 1D lookup table with linear interpolation
-  float interp1(const ColumnVector& x, const ColumnVector& y, float xi);
+  float interp1(const NEWMAT::ColumnVector& x, const NEWMAT::ColumnVector& y, float xi);
 
-  float quantile(const ColumnVector& in, int which);
-  float percentile(const ColumnVector& in, float p);
-  inline float median(const ColumnVector& in){ return quantile(in,2);}
-  inline float iqr(const ColumnVector &in) { return quantile(in,3) - quantile(in,1); }
+  float quantile(const NEWMAT::ColumnVector& in, int which);
+  float percentile(const NEWMAT::ColumnVector& in, float p);
+  inline float median(const NEWMAT::ColumnVector& in){ return quantile(in,2);}
+  inline float iqr(const NEWMAT::ColumnVector &in) { return quantile(in,3) - quantile(in,1); }
 
-  ReturnMatrix quantile(const Matrix& in, int which);
-  ReturnMatrix percentile(const Matrix& in, float p);
-  inline ReturnMatrix median(const Matrix& in){ return quantile(in,2);}
-  inline ReturnMatrix iqr(const Matrix& in){ Matrix res = quantile(in,3) - quantile(in,1); res.Release(); return res;}
+  NEWMAT::ReturnMatrix quantile(const NEWMAT::Matrix& in, int which);
+  NEWMAT::ReturnMatrix percentile(const NEWMAT::Matrix& in, float p);
+  inline NEWMAT::ReturnMatrix median(const NEWMAT::Matrix& in){ return quantile(in,2);}
+  inline NEWMAT::ReturnMatrix iqr(const NEWMAT::Matrix& in){ NEWMAT::Matrix res = quantile(in,3) - quantile(in,1); res.Release(); return res;}
 
-  void cart2sph(const ColumnVector& dir, float& th, float& ph);// cartesian to sperical polar coordinates
-  void cart2sph(const Matrix& dir,ColumnVector& th,ColumnVector& ph);//ditto
-  void cart2sph(const vector<ColumnVector>& dir,ColumnVector& th,ColumnVector& ph);// same but in a vector
+  void cart2sph(const NEWMAT::ColumnVector& dir, float& th, float& ph);// cartesian to sperical polar coordinates
+  void cart2sph(const NEWMAT::Matrix& dir,NEWMAT::ColumnVector& th,NEWMAT::ColumnVector& ph);//ditto
+  void cart2sph(const std::vector<NEWMAT::ColumnVector>& dir,NEWMAT::ColumnVector& th,NEWMAT::ColumnVector& ph);// same but in a vector
 
   // geometry function
-  inline float point_plane_distance(const ColumnVector& X,const ColumnVector& P){//plane defined by a,b,c,d with a^2+b^2+c^2=1
+  inline float point_plane_distance(const NEWMAT::ColumnVector& X,const NEWMAT::ColumnVector& P){//plane defined by a,b,c,d with a^2+b^2+c^2=1
     return( dot(X,P.SubMatrix(1,3,1,1))+P(4) );
   }
 
@@ -210,75 +207,75 @@ namespace MISCMATHS {
 
   // Auto-correlation function estimate of columns of p_ts
   // gives unbiased estimate - scales the raw correlation by 1/(N-abs(lags))
-  void xcorr(const Matrix& p_ts, Matrix& ret, int lag = 0, int p_zeropad = 0);
-  ReturnMatrix xcorr(const Matrix& p_ts, int lag = 0, int p_zeropad = 0);
+  void xcorr(const NEWMAT::Matrix& p_ts, NEWMAT::Matrix& ret, int lag = 0, int p_zeropad = 0);
+  NEWMAT::ReturnMatrix xcorr(const NEWMAT::Matrix& p_ts, int lag = 0, int p_zeropad = 0);
 
   // removes trend from columns of p_ts
   // if p_level==0 it just removes the mean
   // if p_level==1 it removes linear trend
   // if p_level==2 it removes quadratic trend
-  void detrend(Matrix& p_ts, int p_level=1);
-
-  ReturnMatrix zeros(const int dim1, const int dim2 = -1);
-  ReturnMatrix ones(const int dim1, const int dim2 = -1);
-  ReturnMatrix repmat(const Matrix& mat, const int rows = 1, const int cols = 1);
-  ReturnMatrix dist2(const Matrix& mat1, const Matrix& mat2);
-  ReturnMatrix abs(const Matrix& mat);
-  void abs_econ(Matrix& mat);
-  ReturnMatrix sqrt(const Matrix& mat);
-  void sqrt_econ(Matrix& mat);
-  ReturnMatrix sqrtm(const Matrix& mat);
-  ReturnMatrix log(const Matrix& mat);
-  void log_econ(Matrix& mat);
-  ReturnMatrix exp(const Matrix& mat);
-  void exp_econ(Matrix& mat);
-  ReturnMatrix expm(const Matrix& mat);
-  ReturnMatrix tanh(const Matrix& mat);
-  void tanh_econ(Matrix& mat);
-  ReturnMatrix pow(const Matrix& mat, const double exp);
-  void pow_econ(Matrix& mat, const double exp);
-  ReturnMatrix sum(const Matrix& mat, const int dim = 1);
-  ReturnMatrix mean(const Matrix& mat, const int dim = 1);
-  ReturnMatrix mean(const Matrix& mat, const RowVector& weights, const int dim=1);
-  ReturnMatrix var(const Matrix& mat, const int dim = 1);
-  ReturnMatrix max(const Matrix& mat);
-  ReturnMatrix max(const Matrix& mat,ColumnVector& index);
-  ReturnMatrix min(const Matrix& mat);
-  ReturnMatrix gt(const Matrix& mat1,const Matrix& mat2);
-  ReturnMatrix lt(const Matrix& mat1,const Matrix& mat2);
-  ReturnMatrix geqt(const Matrix& mat1,const Matrix& mat2);
-  ReturnMatrix geqt(const Matrix& mat1,const float a);
-  ReturnMatrix leqt(const Matrix& mat1,const Matrix& mat2);
-  ReturnMatrix eq(const Matrix& mat1,const Matrix& mat2);
-  ReturnMatrix neq(const Matrix& mat1,const Matrix& mat2);
-  ReturnMatrix SD(const Matrix& mat1,const Matrix& mat2); // Schur (element-wise) divide
-  void SD_econ(Matrix& mat1,const Matrix& mat2); // Schur (element-wise) divide
-  void SP_econ(Matrix& mat1,const Matrix& mat2); // Schur (element-wise) divide
-
-  ReturnMatrix vox_to_vox(const ColumnVector& xyz1,const ColumnVector& dims1,const ColumnVector& dims2,const Matrix& xfm);
-  ReturnMatrix mni_to_imgvox(const ColumnVector& mni,const ColumnVector& mni_origin,const Matrix& mni2img, const ColumnVector& img_dims);
-
-  void remmean_econ(Matrix& mat, const int dim = 1);
-  void remmean(Matrix& mat, Matrix& Mean, const int dim = 1);
-  void remmean(const Matrix& mat, Matrix& demeanedmat, Matrix& Mean,  const int dim = 1);
-  ReturnMatrix remmean(const Matrix& mat, const int dim = 1);
-
-  ReturnMatrix stdev(const Matrix& mat, const int dim = 1);
-  ReturnMatrix cov(const Matrix& mat, const bool sampleCovariance = false, const int econ=20000);
-  ReturnMatrix cov_r(const Matrix& mat, const bool sampleCovariance = false, const int econ=20000);
-  ReturnMatrix cov_r(const Matrix& data, const Matrix& weights, int econ=20000);
-
-
-
-  ReturnMatrix oldcov(const Matrix& mat, const bool norm = false);
-  ReturnMatrix corrcoef(const Matrix& mat, const bool norm = false);
-  void symm_orth(Matrix &Mat);
-  void powerspectrum(const Matrix &Mat1, Matrix &Result, bool useLog);
-  void element_mod_n(Matrix& Mat,double n); //represent each element in modulo n (useful for wrapping phases (n=2*M_PI))
+  void detrend(NEWMAT::Matrix& p_ts, int p_level=1);
+
+  NEWMAT::ReturnMatrix zeros(const int dim1, const int dim2 = -1);
+  NEWMAT::ReturnMatrix ones(const int dim1, const int dim2 = -1);
+  NEWMAT::ReturnMatrix repmat(const NEWMAT::Matrix& mat, const int rows = 1, const int cols = 1);
+  NEWMAT::ReturnMatrix dist2(const NEWMAT::Matrix& mat1, const NEWMAT::Matrix& mat2);
+  NEWMAT::ReturnMatrix abs(const NEWMAT::Matrix& mat);
+  void abs_econ(NEWMAT::Matrix& mat);
+  NEWMAT::ReturnMatrix sqrt(const NEWMAT::Matrix& mat);
+  void sqrt_econ(NEWMAT::Matrix& mat);
+  NEWMAT::ReturnMatrix sqrtm(const NEWMAT::Matrix& mat);
+  NEWMAT::ReturnMatrix log(const NEWMAT::Matrix& mat);
+  void log_econ(NEWMAT::Matrix& mat);
+  NEWMAT::ReturnMatrix exp(const NEWMAT::Matrix& mat);
+  void exp_econ(NEWMAT::Matrix& mat);
+  NEWMAT::ReturnMatrix expm(const NEWMAT::Matrix& mat);
+  NEWMAT::ReturnMatrix tanh(const NEWMAT::Matrix& mat);
+  void tanh_econ(NEWMAT::Matrix& mat);
+  NEWMAT::ReturnMatrix pow(const NEWMAT::Matrix& mat, const double exp);
+  void pow_econ(NEWMAT::Matrix& mat, const double exp);
+  NEWMAT::ReturnMatrix sum(const NEWMAT::Matrix& mat, const int dim = 1);
+  NEWMAT::ReturnMatrix mean(const NEWMAT::Matrix& mat, const int dim = 1);
+  NEWMAT::ReturnMatrix mean(const NEWMAT::Matrix& mat, const NEWMAT::RowVector& weights, const int dim=1);
+  NEWMAT::ReturnMatrix var(const NEWMAT::Matrix& mat, const int dim = 1);
+  NEWMAT::ReturnMatrix max(const NEWMAT::Matrix& mat);
+  NEWMAT::ReturnMatrix max(const NEWMAT::Matrix& mat,NEWMAT::ColumnVector& index);
+  NEWMAT::ReturnMatrix min(const NEWMAT::Matrix& mat);
+  NEWMAT::ReturnMatrix gt(const NEWMAT::Matrix& mat1,const NEWMAT::Matrix& mat2);
+  NEWMAT::ReturnMatrix lt(const NEWMAT::Matrix& mat1,const NEWMAT::Matrix& mat2);
+  NEWMAT::ReturnMatrix geqt(const NEWMAT::Matrix& mat1,const NEWMAT::Matrix& mat2);
+  NEWMAT::ReturnMatrix geqt(const NEWMAT::Matrix& mat1,const float a);
+  NEWMAT::ReturnMatrix leqt(const NEWMAT::Matrix& mat1,const NEWMAT::Matrix& mat2);
+  NEWMAT::ReturnMatrix eq(const NEWMAT::Matrix& mat1,const NEWMAT::Matrix& mat2);
+  NEWMAT::ReturnMatrix neq(const NEWMAT::Matrix& mat1,const NEWMAT::Matrix& mat2);
+  NEWMAT::ReturnMatrix SD(const NEWMAT::Matrix& mat1,const NEWMAT::Matrix& mat2); // Schur (element-wise) divide
+  void SD_econ(NEWMAT::Matrix& mat1,const NEWMAT::Matrix& mat2); // Schur (element-wise) divide
+  void SP_econ(NEWMAT::Matrix& mat1,const NEWMAT::Matrix& mat2); // Schur (element-wise) divide
+
+  NEWMAT::ReturnMatrix vox_to_vox(const NEWMAT::ColumnVector& xyz1,const NEWMAT::ColumnVector& dims1,const NEWMAT::ColumnVector& dims2,const NEWMAT::Matrix& xfm);
+  NEWMAT::ReturnMatrix mni_to_imgvox(const NEWMAT::ColumnVector& mni,const NEWMAT::ColumnVector& mni_origin,const NEWMAT::Matrix& mni2img, const NEWMAT::ColumnVector& img_dims);
+
+  void remmean_econ(NEWMAT::Matrix& mat, const int dim = 1);
+  void remmean(NEWMAT::Matrix& mat, NEWMAT::Matrix& Mean, const int dim = 1);
+  void remmean(const NEWMAT::Matrix& mat, NEWMAT::Matrix& demeanedmat, NEWMAT::Matrix& Mean,  const int dim = 1);
+  NEWMAT::ReturnMatrix remmean(const NEWMAT::Matrix& mat, const int dim = 1);
+
+  NEWMAT::ReturnMatrix stdev(const NEWMAT::Matrix& mat, const int dim = 1);
+  NEWMAT::ReturnMatrix cov(const NEWMAT::Matrix& mat, const bool sampleCovariance = false, const int econ=20000);
+  NEWMAT::ReturnMatrix cov_r(const NEWMAT::Matrix& mat, const bool sampleCovariance = false, const int econ=20000);
+  NEWMAT::ReturnMatrix cov_r(const NEWMAT::Matrix& data, const NEWMAT::Matrix& weights, int econ=20000);
+
+
+
+  NEWMAT::ReturnMatrix oldcov(const NEWMAT::Matrix& mat, const bool norm = false);
+  NEWMAT::ReturnMatrix corrcoef(const NEWMAT::Matrix& mat, const bool norm = false);
+  void symm_orth(NEWMAT::Matrix &Mat);
+  void powerspectrum(const NEWMAT::Matrix &Mat1, NEWMAT::Matrix &Result, bool useLog);
+  void element_mod_n(NEWMAT::Matrix& Mat,double n); //represent each element in modulo n (useful for wrapping phases (n=2*M_PI))
 
   // matlab-like flip function
-  ReturnMatrix flipud(const Matrix& mat);
-  ReturnMatrix fliplr(const Matrix& mat);
+  NEWMAT::ReturnMatrix flipud(const NEWMAT::Matrix& mat);
+  NEWMAT::ReturnMatrix fliplr(const NEWMAT::Matrix& mat);
 
   // ols
   // data is t x v
@@ -286,24 +283,24 @@ namespace MISCMATHS {
   // tc is cons x ev (contrast matrix)
   // cope and varcope will be cons x v
   // but will be resized if they are wrong
-  void ols(const Matrix& data,const Matrix& des,const Matrix& tc, Matrix& cope,Matrix& varcope);
-  float ols_dof(const Matrix& des);
+  void ols(const NEWMAT::Matrix& data,const NEWMAT::Matrix& des,const NEWMAT::Matrix& tc, NEWMAT::Matrix& cope,NEWMAT::Matrix& varcope);
+  float ols_dof(const NEWMAT::Matrix& des);
 
 
   // Conjugate Gradient methods to solve for x in:   A * x = b
   // A must be symmetric and positive definite
-  int conjgrad(ColumnVector& x, const Matrix& A, const ColumnVector& b,
+  int conjgrad(NEWMAT::ColumnVector& x, const NEWMAT::Matrix& A, const NEWMAT::ColumnVector& b,
 	       int maxit=3);
   // allow specification of reltol = relative tolerance of residual error
   //  (stops when error < reltol * initial error)
-  int conjgrad(ColumnVector& x, const Matrix& A, const ColumnVector& b,
+  int conjgrad(NEWMAT::ColumnVector& x, const NEWMAT::Matrix& A, const NEWMAT::ColumnVector& b,
 	       int maxit, float reltol);
 
-  float csevl(const float x, const ColumnVector& cs, const int n);
+  float csevl(const float x, const NEWMAT::ColumnVector& cs, const int n);
   float digamma(const float x);
-  void glm_vb(const Matrix& X, const ColumnVector& Y, ColumnVector& B, SymmetricMatrix& ilambda_B, int niters=20);
+  void glm_vb(const NEWMAT::Matrix& X, const NEWMAT::ColumnVector& Y, NEWMAT::ColumnVector& B, NEWMAT::SymmetricMatrix& ilambda_B, int niters=20);
 
-  vector<float> ColumnVector2vector(const ColumnVector& col);
+  vector<float> ColumnVector2vector(const NEWMAT::ColumnVector& col);
 
   ///////////////////////////////////////////////////////////////////////////
   // Uninteresting byte swapping functions
@@ -318,9 +315,9 @@ namespace MISCMATHS {
 
   // TEMPLATE DEFINITIONS //
 
-  template<class t> ReturnMatrix vector2ColumnVector(const vector<t>& vec)
+  template<class t> NEWMAT::ReturnMatrix vector2ColumnVector(const std::vector<t>& vec)
   {
-    ColumnVector col(vec.size());
+    NEWMAT::ColumnVector col(vec.size());
 
     for(unsigned int c = 0; c < vec.size(); c++)
       col(c+1) = vec[c];
@@ -329,22 +326,22 @@ namespace MISCMATHS {
     return col;
   }
 
-  template<class t> void write_vector(const string& fname, const vector<t>& vec)
+  template<class t> void write_vector(const std::string& fname, const std::vector<t>& vec)
   {
-    ofstream out;
+    std::ofstream out;
     out.open(fname.c_str(), ios::out);
     copy(vec.begin(), vec.end(), ostream_iterator<t>(out, " "));
   }
 
-  template<class t> void write_vector(const vector<t>& vec, const string& fname)
+  template<class t> void write_vector(const std::vector<t>& vec, const std::string& fname)
   {
     write_vector(fname,vec);
   }
 
   template <class T>
-  string num2str(T n, int width)
+  std::string num2str(T n, int width)
   {
-    ostringstream os;
+    std::ostringstream os;
     if (width>0) {
       os.fill('0');
       os.width(width);
diff --git a/miscprob.h b/miscprob.h
index 59bdbc4ae7a8d45057194a9ba5582f3156e05a8a..c69e2860cf7d176c57f027aee49717e0b6ed3102 100644
--- a/miscprob.h
+++ b/miscprob.h
@@ -16,48 +16,46 @@
 #include "cprob/libprob.h"
 #include "stdlib.h"
 
-using namespace NEWMAT;
-
 namespace MISCMATHS {
 
 //   ReturnMatrix betarnd(const int dim1, const int dim2,
 // 		       const float a, const float b);
 
-  ReturnMatrix betapdf(const RowVector& vals,
+  NEWMAT::ReturnMatrix betapdf(const NEWMAT::RowVector& vals,
 		       const float a, const float b);
 
-  ReturnMatrix unifrnd(const int dim1 = 1, const int dim2 = -1,
+  NEWMAT::ReturnMatrix unifrnd(const int dim1 = 1, const int dim2 = -1,
 		       const float start = 0, const float end = 1);
 
-  ReturnMatrix normrnd(const int dim1 = 1, const int dim2 = -1,
+  NEWMAT::ReturnMatrix normrnd(const int dim1 = 1, const int dim2 = -1,
 		       const float mu = 0, const float sigma = 1);
 
   // returns nsamps*nparams matrix:
-  ReturnMatrix mvnrnd(const RowVector& mu, const SymmetricMatrix& covar, int nsamp = 1);
+  NEWMAT::ReturnMatrix mvnrnd(const NEWMAT::RowVector& mu, const NEWMAT::SymmetricMatrix& covar, int nsamp = 1);
 
-  float mvnpdf(const RowVector& vals, const RowVector& mu, const SymmetricMatrix& covar);
+  float mvnpdf(const NEWMAT::RowVector& vals, const NEWMAT::RowVector& mu, const NEWMAT::SymmetricMatrix& covar);
 
-  float bvnpdf(const RowVector& vals, const RowVector& mu, const SymmetricMatrix& covar);
+  float bvnpdf(const NEWMAT::RowVector& vals, const NEWMAT::RowVector& mu, const NEWMAT::SymmetricMatrix& covar);
 
   float normpdf(const float val, const float mu = 0, const float var = 1);
   float lognormpdf(const float val, const float mu = 0, const float var = 1);
 
-  ReturnMatrix normpdf(const RowVector& vals, const float mu = 0, const float var = 1);
+  NEWMAT::ReturnMatrix normpdf(const NEWMAT::RowVector& vals, const float mu = 0, const float var = 1);
 
-  ReturnMatrix normpdf(const RowVector& vals, const RowVector& mus,
-		       const RowVector& vars);
+  NEWMAT::ReturnMatrix normpdf(const NEWMAT::RowVector& vals, const NEWMAT::RowVector& mus,
+		       const NEWMAT::RowVector& vars);
 
-  ReturnMatrix normcdf(const RowVector& vals, const float mu = 0, const float var = 1);
+  NEWMAT::ReturnMatrix normcdf(const NEWMAT::RowVector& vals, const float mu = 0, const float var = 1);
 
-  ReturnMatrix gammapdf(const RowVector& vals, const float mu = 0, const float var = 1);
+  NEWMAT::ReturnMatrix gammapdf(const NEWMAT::RowVector& vals, const float mu = 0, const float var = 1);
 
-  ReturnMatrix gammacdf(const RowVector& vals, const float mu = 0, const float var = 1);
+  NEWMAT::ReturnMatrix gammacdf(const NEWMAT::RowVector& vals, const float mu = 0, const float var = 1);
 
-//   ReturnMatrix gammarnd(const int dim1, const int dim2,
+//   NEWMAT::ReturnMatrix gammarnd(const int dim1, const int dim2,
 // 			const float a, const float b);
 
   // returns n! * n matrix of all possible permutations
-  ReturnMatrix perms(const int n);
+  NEWMAT::ReturnMatrix perms(const int n);
 
 
   class Mvnormrandm
@@ -65,42 +63,42 @@ namespace MISCMATHS {
     public:
       Mvnormrandm(){}
 
-      Mvnormrandm(const RowVector& pmu, const SymmetricMatrix& pcovar) :
+      Mvnormrandm(const NEWMAT::RowVector& pmu, const NEWMAT::SymmetricMatrix& pcovar) :
 	mu(pmu),
 	covar(pcovar)
 	{
-	  Matrix eig_vec;
-	  DiagonalMatrix eig_val;
+	  NEWMAT::Matrix eig_vec;
+	  NEWMAT::DiagonalMatrix eig_val;
 	  EigenValues(covar,eig_val,eig_vec);
 
 	  covarw = sqrt(eig_val)*eig_vec.t();
 	}
 
-      ReturnMatrix next(int nsamp = 1) const
+      NEWMAT::ReturnMatrix next(int nsamp = 1) const
 	{
-	  Matrix ret = ones(nsamp, 1)*mu + normrnd(nsamp,mu.Ncols())*covarw;
+	  NEWMAT::Matrix ret = ones(nsamp, 1)*mu + normrnd(nsamp,mu.Ncols())*covarw;
 	  ret.Release();
 	  return ret;
 	}
 
-      ReturnMatrix next(const RowVector& pmu, int nsamp = 1)
+      NEWMAT::ReturnMatrix next(const NEWMAT::RowVector& pmu, int nsamp = 1)
 	{
 	  mu=pmu;
 
-	  Matrix ret = ones(nsamp, 1)*mu + normrnd(nsamp,mu.Ncols())*covarw;
+	  NEWMAT::Matrix ret = ones(nsamp, 1)*mu + normrnd(nsamp,mu.Ncols())*covarw;
 	  ret.Release();
 	  return ret;
 	}
 
-      void setcovar(const SymmetricMatrix& pcovar)
+      void setcovar(const NEWMAT::SymmetricMatrix& pcovar)
 	{
 	  covar=pcovar;
 
 	  mu.ReSize(covar.Nrows());
 	  mu=0;
 
-	  Matrix eig_vec;
-	  DiagonalMatrix eig_val;
+	  NEWMAT::Matrix eig_vec;
+	  NEWMAT::DiagonalMatrix eig_val;
 	  EigenValues(covar,eig_val,eig_vec);
 
 	  covarw = sqrt(eig_val)*eig_vec.t();
@@ -108,10 +106,10 @@ namespace MISCMATHS {
 
     private:
 
-      RowVector mu;
-      SymmetricMatrix covar;
+      NEWMAT::RowVector mu;
+      NEWMAT::SymmetricMatrix covar;
 
-      Matrix covarw;
+      NEWMAT::Matrix covarw;
 
     };
 }
diff --git a/optimise.h b/optimise.h
index e3a0de28cb12e6f2d733cdc08f9499d6b8ac8c2d..998e6b6fae9f6ac9eaf792e4f85a914589ccb76c 100644
--- a/optimise.h
+++ b/optimise.h
@@ -16,19 +16,17 @@
 #include "armawrap/newmatap.h"
 #include "string"
 
-using namespace NEWMAT;
-
 namespace MISCMATHS {
 
-float optimise1d(ColumnVector &pt, const ColumnVector dir,
-		const ColumnVector tol, int &iterations_done,
-		float (*func)(const ColumnVector &), int max_iter,
+  float optimise1d(NEWMAT::ColumnVector &pt, const NEWMAT::ColumnVector dir,
+		const NEWMAT::ColumnVector tol, int &iterations_done,
+		float (*func)(const NEWMAT::ColumnVector &), int max_iter,
 		float &init_value, float boundguess);
 
 
- float optimise(ColumnVector &pt, int numopt, const ColumnVector &tol,
-		float (*func)(const ColumnVector &), int &iterations_done,
-		int max_iter, const ColumnVector& boundguess,
+ float optimise(NEWMAT::ColumnVector &pt, int numopt, const NEWMAT::ColumnVector &tol,
+		float (*func)(const NEWMAT::ColumnVector &), int &iterations_done,
+		int max_iter, const NEWMAT::ColumnVector& boundguess,
 		const std::string type="brent");
 
 }
diff --git a/rungekutta.h b/rungekutta.h
index 3af3ba62b6ce9932d6007ed4f1fd9a47c194d01c..7e229f740c139eb4f07016697f80759c7e88f13a 100644
--- a/rungekutta.h
+++ b/rungekutta.h
@@ -17,8 +17,6 @@
 #include "armawrap/newmatap.h"
 #include "armawrap/newmatio.h"
 
-using namespace NEWMAT;
-
 namespace MISCMATHS {
 
 class Derivative
@@ -29,20 +27,20 @@ public:
   // x is time point to evaluate at
   // y is state variables
   // paramvalues are "constants" in the diff eqn
-  virtual const ColumnVector& evaluate(float x,const ColumnVector& y,const ColumnVector& paramvalues) const = 0;
+  virtual const NEWMAT::ColumnVector& evaluate(float x,const NEWMAT::ColumnVector& y,const NEWMAT::ColumnVector& paramvalues) const = 0;
 
   virtual ~Derivative(){};
 
 protected:
   int ny;
-  mutable ColumnVector dy;
+  mutable NEWMAT::ColumnVector dy;
 };
 
-void rk(ColumnVector& ret, const ColumnVector& y, const ColumnVector& dy, float x, float h, const Derivative& deriv,const ColumnVector& paramvalues);
+void rk(NEWMAT::ColumnVector& ret, const NEWMAT::ColumnVector& y, const NEWMAT::ColumnVector& dy, float x, float h, const Derivative& deriv,const NEWMAT::ColumnVector& paramvalues);
 
-void rkqc(ColumnVector& y, float& x, float& hnext, ColumnVector& dy, float htry, float eps, const Derivative& deriv,const ColumnVector& paramvalues);
+void rkqc(NEWMAT::ColumnVector& y, float& x, float& hnext, NEWMAT::ColumnVector& dy, float htry, float eps, const Derivative& deriv,const NEWMAT::ColumnVector& paramvalues);
 
-void runge_kutta(Matrix& yp, ColumnVector& xp, ColumnVector& hp, const ColumnVector& ystart, float x1, float x2, float eps, float hmin, const Derivative& deriv,const ColumnVector& paramvalues);
+void runge_kutta(NEWMAT::Matrix& yp, NEWMAT::ColumnVector& xp, NEWMAT::ColumnVector& hp, const NEWMAT::ColumnVector& ystart, float x1, float x2, float eps, float hmin, const Derivative& deriv,const NEWMAT::ColumnVector& paramvalues);
 
 }
 
diff --git a/sparse_matrix.h b/sparse_matrix.h
index a8ffac0c280802ab57f310344df5ef546026c119..e2c4f652a9aaaf73cdfbb6fbffead0afad82c583 100644
--- a/sparse_matrix.h
+++ b/sparse_matrix.h
@@ -16,16 +16,13 @@
 #include "armawrap/newmat.h"
 #include "armawrap/newmatio.h"
 
-using namespace NEWMAT;
-using namespace std;
-
 namespace MISCMATHS {
 
   class SparseMatrix
     {
     public:
 
-      typedef map<int,double> Row;
+      typedef std::map<int,double> Row;
 
       SparseMatrix() : nrows(0), ncols(0) {}
 
@@ -45,12 +42,12 @@ namespace MISCMATHS {
 	  return *this;
 	}
 
-      SparseMatrix(const Matrix& pmatin)
+      SparseMatrix(const NEWMAT::Matrix& pmatin)
 	{
 	  operator=(pmatin);
 	}
 
-      const SparseMatrix& operator=(const Matrix& pmatin);
+      const SparseMatrix& operator=(const NEWMAT::Matrix& pmatin);
 
       //      void ReSize(int pnrows, int pncols)
       void ReSize(int pnrows, int pncols);
@@ -62,16 +59,16 @@ namespace MISCMATHS {
 
       void transpose(SparseMatrix& ret);
 
-      ReturnMatrix RowAsColumn(int r) const;
+      NEWMAT::ReturnMatrix RowAsColumn(int r) const;
 
       int maxnonzerosinrow() const;
 
-      void permute(const ColumnVector& p, SparseMatrix& pA);
+      void permute(const NEWMAT::ColumnVector& p, SparseMatrix& pA);
 
       const double operator()(int x, int y) const
 	{
 	  double ret = 0.0;
-	  map<int,double>::const_iterator it=data[x-1].find(y-1);
+      std::map<int,double>::const_iterator it=data[x-1].find(y-1);
 	  if(it != data[x-1].end())
 	    ret = (*it).second;
 
@@ -111,7 +108,7 @@ namespace MISCMATHS {
 
       const Row& row(int r) const { return data[r-1]; }
 
-      ReturnMatrix AsMatrix() const;
+      NEWMAT::ReturnMatrix AsMatrix() const;
 
       int Nrows() const { return nrows; }
       int Ncols() const { return ncols; }
@@ -128,20 +125,20 @@ namespace MISCMATHS {
       int nrows;
       int ncols;
 
-      vector<map<int,double> > data;
+      std::vector<map<int,double> > data;
 
     };
 
   void multiply(const SparseMatrix& lm, const SparseMatrix& rm, SparseMatrix& ret);
-  void multiply(const DiagonalMatrix& lm, const SparseMatrix& rm, SparseMatrix& ret);
+  void multiply(const NEWMAT::DiagonalMatrix& lm, const SparseMatrix& rm, SparseMatrix& ret);
 
-  void multiply(const SparseMatrix& lm, const ColumnVector& rm, ColumnVector& ret);
+  void multiply(const SparseMatrix& lm, const NEWMAT::ColumnVector& rm, NEWMAT::ColumnVector& ret);
 
-  void multiply(const SparseMatrix& lm, const SparseMatrix::Row& rm, ColumnVector& ret);
+  void multiply(const SparseMatrix& lm, const SparseMatrix::Row& rm, NEWMAT::ColumnVector& ret);
 
   void add(const SparseMatrix& lm, const SparseMatrix& rm, SparseMatrix& ret);
 
-  void colvectosparserow(const ColumnVector& col, SparseMatrix::Row& row);
+  void colvectosparserow(const NEWMAT::ColumnVector& col, SparseMatrix::Row& row);
 
   void vertconcat(const SparseMatrix& A, const SparseMatrix& B, SparseMatrix& ret);
 
diff --git a/sparsefn.h b/sparsefn.h
index 5f7a887e3f9e2395a82b5420d33b9b997daa6e7f..51054aae018fcd23462c631bb89da3b6e67e157a 100644
--- a/sparsefn.h
+++ b/sparsefn.h
@@ -18,26 +18,24 @@
 #include "sparse_matrix.h"
 #include "armawrap/newmat.h"
 
-using namespace NEWMAT;
-
 namespace MISCMATHS {
 
-  float quadratic(const ColumnVector& m, const SparseMatrix& C);
+  float quadratic(const NEWMAT::ColumnVector& m, const SparseMatrix& C);
   void speye(int n, SparseMatrix& ret);
   void chol(const SparseMatrix& A, SparseMatrix& U, SparseMatrix& L);
   void inv(const SparseMatrix& U, const SparseMatrix& L, SparseMatrix& ret);
   void solvefortracex(const SparseMatrix& U, const SparseMatrix& L, const SparseMatrix& b1, const SparseMatrix& b2, float& tr1, float& tr2);
-  void solveforx(const SparseMatrix& U, const SparseMatrix& L, const ColumnVector& b, ColumnVector& x);
-  void solveforx(const SparseMatrix& A, const ColumnVector& b, ColumnVector& x, float tol = 0.001, int kmax = 500);
-  void solveforx(const SparseMatrix& A, const ColumnVector& b, SparseMatrix& x);
+  void solveforx(const SparseMatrix& U, const SparseMatrix& L, const NEWMAT::ColumnVector& b, NEWMAT::ColumnVector& x);
+  void solveforx(const SparseMatrix& A, const NEWMAT::ColumnVector& b, NEWMAT::ColumnVector& x, float tol = 0.001, int kmax = 500);
+  void solveforx(const SparseMatrix& A, const NEWMAT::ColumnVector& b, SparseMatrix& x);
   void solveforx(const SparseMatrix& A, const SparseMatrix& b, SparseMatrix& x);
   float solvefortracex(const SparseMatrix& A, const SparseMatrix& b, SparseMatrix& x, int nsamps = 50, float tol = 0.001);
-  void solve(const SparseMatrix& A, const Matrix& b, SparseMatrix& x);
+  void solve(const SparseMatrix& A, const NEWMAT::Matrix& b, SparseMatrix& x);
   void addto(SparseMatrix& A, const SparseMatrix& B, float S);
   void symmetric_addto(SparseMatrix& A, const SparseMatrix& B, float S);
   void addto(const SparseMatrix::Row& A, const SparseMatrix::Row& B, float S);
-  void addto(SparseMatrix& A, const Matrix& B);
-  void cov(const ColumnVector& A, SparseMatrix& ret);
+  void addto(SparseMatrix& A, const NEWMAT::Matrix& B);
+  void cov(const NEWMAT::ColumnVector& A, SparseMatrix& ret);
 }
 
 #endif
diff --git a/t2z.h b/t2z.h
index a5bf6f1e791e667792d4909cf256705e11ceacfe..963874db4f44ff5239f113b1e8b8e848ddb8a2ef 100644
--- a/t2z.h
+++ b/t2z.h
@@ -15,8 +15,6 @@
 #include "armawrap/newmatio.h"
 #include "base2z.h"
 
-using namespace NEWMAT;
-
 namespace MISCMATHS {
 
   class T2z : public Base2z
@@ -28,9 +26,9 @@ namespace MISCMATHS {
       float convert(float t, int dof,double *newp=NULL);
       float converttologp(float t, int dof);
 
-      static void ComputePs(const ColumnVector& p_vars, const ColumnVector& p_cbs, int p_dof, ColumnVector& p_ps);
-      static void ComputeZStats(const ColumnVector& p_vars, const ColumnVector& p_cbs, int p_dof, ColumnVector& p_zs);
-      static void ComputeZStats(const ColumnVector& p_vars, const ColumnVector& p_cbs, const ColumnVector& p_dof, ColumnVector& p_zs);
+      static void ComputePs(const NEWMAT::ColumnVector& p_vars, const NEWMAT::ColumnVector& p_cbs, int p_dof, NEWMAT::ColumnVector& p_ps);
+      static void ComputeZStats(const NEWMAT::ColumnVector& p_vars, const NEWMAT::ColumnVector& p_cbs, int p_dof, NEWMAT::ColumnVector& p_zs);
+      static void ComputeZStats(const NEWMAT::ColumnVector& p_vars, const NEWMAT::ColumnVector& p_cbs, const NEWMAT::ColumnVector& p_dof, NEWMAT::ColumnVector& p_zs);
 
     private:
       T2z() : Base2z()