diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 5f79b9b2fa8a2cb0ccdafd28fe8eff10f45f257d..737a4b4870f48c15ad90faa24d9670a8d1c94325 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1,4 +1,4 @@
-// $Id: GModel.cpp,v 1.36 2007-03-18 23:02:26 geuzaine Exp $
+// $Id: GModel.cpp,v 1.37 2007-04-08 23:06:53 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -121,9 +121,38 @@ void GModel::removeInvisibleElements()
     removeInvisible((*it)->lines, all);
     if((*it)->meshRep) (*it)->meshRep->destroy();
   }
+}
+
+template<class T>
+static void associateEntityWithElementVertices(GEntity *ge, std::vector<T*> &elements)
+{
+  for(unsigned int i = 0; i < elements.size(); i++)
+    for(int j = 0; j < elements[i]->getNumVertices(); j++)
+      elements[i]->getVertex(j)->setEntity(ge);
+}
 
-  // FIXME: loop over all mesh vertices in edges, faces and regions
-  // and remove all unused (nonconnected) vertices
+void GModel::associateEntityWithVertices()
+{
+  // loop on regions, then on faces, edges and vertices and store the
+  // entity pointer in the the elements' vertices (this way we
+  // associate the entity of lowest geometrical degree with each
+  // vertex)
+  for(riter it = firstRegion(); it != lastRegion(); ++it){
+    associateEntityWithElementVertices(*it, (*it)->tetrahedra);
+    associateEntityWithElementVertices(*it, (*it)->hexahedra);
+    associateEntityWithElementVertices(*it, (*it)->prisms);
+    associateEntityWithElementVertices(*it, (*it)->pyramids);
+  }
+  for(fiter it = firstFace(); it != lastFace(); ++it){
+    associateEntityWithElementVertices(*it, (*it)->triangles);
+    associateEntityWithElementVertices(*it, (*it)->quadrangles);
+  }
+  for(eiter it = firstEdge(); it != lastEdge(); ++it){
+    associateEntityWithElementVertices(*it, (*it)->lines);
+  }
+  for(viter it = firstVertex(); it != lastVertex(); ++it){
+    (*it)->mesh_vertices[0]->setEntity(*it);
+  }
 }
 
 int GModel::renumberMeshVertices()
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 4140af7a5a62764164ddaef40016e7e82911f7a6..8aa8b1f554c0c84bfca95c14ac9450ff6e618049 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -98,6 +98,9 @@ class GModel
   void remove(GEdge *e) { edges.erase(std::find(firstEdge(), lastEdge(), e)); }
   void remove(GVertex *v) { vertices.erase(std::find(firstVertex(), lastVertex(), v)); }
 
+  // loop over all vertices connected to elements and associate geo entity
+  void associateEntityWithVertices();
+
   // Renumber all the mesh vertices in a continuous sequence
   int renumberMeshVertices();
 
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 43d2145e7c66b17d9dd743877e5030a493eceb5c..17cf823d430a1c011acdf1295eba0633be3d33d4 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Mesh.cpp,v 1.13 2007-04-04 12:18:54 geuzaine Exp $
+// $Id: GModelIO_Mesh.cpp,v 1.14 2007-04-08 23:06:53 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -95,38 +95,6 @@ static void storeElementsInEntities(GModel *m,
   }
 }
 
-template<class T>
-static void associateEntityWithVertices(GEntity *ge, std::vector<T*> &elements)
-{
-  for(unsigned int i = 0; i < elements.size(); i++)
-    for(int j = 0; j < elements[i]->getNumVertices(); j++)
-      elements[i]->getVertex(j)->setEntity(ge);
-}
-
-static void associateEntityWithVertices(GModel *m)
-{
-  // loop on regions, then on faces, edges and vertices and store the
-  // entity pointer in the the elements' vertices (this way we
-  // associate the entity of lowest geometrical degree with each
-  // vertex)
-  for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
-    associateEntityWithVertices(*it, (*it)->tetrahedra);
-    associateEntityWithVertices(*it, (*it)->hexahedra);
-    associateEntityWithVertices(*it, (*it)->prisms);
-    associateEntityWithVertices(*it, (*it)->pyramids);
-  }
-  for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
-    associateEntityWithVertices(*it, (*it)->triangles);
-    associateEntityWithVertices(*it, (*it)->quadrangles);
-  }
-  for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it){
-    associateEntityWithVertices(*it, (*it)->lines);
-  }
-  for(GModel::viter it = m->firstVertex(); it != m->lastVertex(); ++it){
-    (*it)->mesh_vertices[0]->setEntity(*it);
-  }
-}
-
 static void storeVerticesInEntities(std::map<int, MVertex*> &vertices)
 {
   std::map<int, MVertex*>::const_iterator it = vertices.begin();
@@ -517,7 +485,7 @@ int GModel::readMSH(const std::string &name)
   }
 
   // associate the correct geometrical entity with each mesh vertex
-  associateEntityWithVertices(this);
+  associateEntityWithVertices();
 
   // special case for geometry vertices: now that the correct
   // geometrical entity has been associated with the vertices, we
@@ -1125,7 +1093,7 @@ int GModel::readVRML(const std::string &name)
 
   for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++) 
     storeElementsInEntities(this, elements[i]);
-  associateEntityWithVertices(this);
+  associateEntityWithVertices();
   storeVerticesInEntities(allVertexVector);
 
   fclose(fp);
@@ -1328,7 +1296,7 @@ int GModel::readUNV(const std::string &name)
   
   for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++) 
     storeElementsInEntities(this, elements[i]);
-  associateEntityWithVertices(this);
+  associateEntityWithVertices();
   storeVerticesInEntities(vertexMap);
   for(int i = 0; i < 4; i++)  
     storePhysicalTagsInEntities(this, i, physicals[i]);
@@ -1494,7 +1462,7 @@ int GModel::readMESH(const std::string &name)
 
   for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++) 
     storeElementsInEntities(this, elements[i]);
-  associateEntityWithVertices(this);
+  associateEntityWithVertices();
   storeVerticesInEntities(vertexVector);
 
   fclose(fp);
@@ -1831,7 +1799,7 @@ int GModel::readBDF(const std::string &name)
   
   for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++) 
     storeElementsInEntities(this, elements[i]);
-  associateEntityWithVertices(this);
+  associateEntityWithVertices();
   storeVerticesInEntities(vertexMap);
 
   fclose(fp);
diff --git a/Geo/Makefile b/Geo/Makefile
index 79025c85970e2d03870d2d9702585fcfa9171f8e..52d62ed92fdc58063e90a2a9e66622a7e670f7e4 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.139 2007-03-02 09:20:21 remacle Exp $
+# $Id: Makefile,v 1.140 2007-04-08 23:06:53 geuzaine Exp $
 #
 # Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 #
@@ -43,9 +43,7 @@ SRC = GEntity.cpp\
 	GeoStringInterface.cpp GeoInterpolation.cpp\
       findLinks.cpp\
       MVertex.cpp \
-      MElement.cpp \
-      SVector3.cpp\
-      SBoundingBox3d.cpp
+      MElement.cpp
 
 OBJ = ${SRC:.cpp=.o}
 
@@ -267,5 +265,3 @@ MElement.o: MElement.cpp MElement.h ../Common/GmshDefines.h MVertex.h \
   SPoint3.h MEdge.h SVector3.h ../Common/Hash.h MFace.h \
   ../Numeric/Numeric.h ../Common/Context.h ../DataStr/List.h GEntity.h \
   Range.h SBoundingBox3d.h ../Common/Message.h
-SVector3.o: SVector3.cpp SVector3.h SPoint3.h
-SBoundingBox3d.o: SBoundingBox3d.cpp SBoundingBox3d.h SPoint3.h
diff --git a/Geo/SBoundingBox3d.cpp b/Geo/SBoundingBox3d.cpp
deleted file mode 100644
index 0a7261c4acada0be83b9e7418e097a1554d5783b..0000000000000000000000000000000000000000
--- a/Geo/SBoundingBox3d.cpp
+++ /dev/null
@@ -1,117 +0,0 @@
-// $Id: SBoundingBox3d.cpp,v 1.5 2007-01-10 13:48:46 geuzaine Exp $
-//
-// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
-//
-// This program is free software; you can redistribute it and/or modify
-// it under the terms of the GNU General Public License as published by
-// the Free Software Foundation; either version 2 of the License, or
-// (at your option) any later version.
-//
-// This program is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with this program; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
-// USA.
-// 
-// Please report all bugs and problems to <gmsh@geuz.org>.
-
-#include "SBoundingBox3d.h"
-#include <float.h>
-
-SBoundingBox3d::SBoundingBox3d()
-: MinPt(DBL_MAX,DBL_MAX,DBL_MAX), MaxPt(-DBL_MAX,-DBL_MAX,-DBL_MAX)
-{}
-
-SBoundingBox3d::SBoundingBox3d(const SPoint3 & pt)
-: MinPt(pt), MaxPt(pt)
-{}
-
-void SBoundingBox3d::reset()
-{
-  MinPt = SPoint3(DBL_MAX,DBL_MAX,DBL_MAX);
-  MaxPt = SPoint3(-DBL_MAX,-DBL_MAX,-DBL_MAX);
-}
-
-bool SBoundingBox3d::empty()
-{
-  if(MinPt.x() == DBL_MAX || MinPt.y() == DBL_MAX || MinPt.z() == DBL_MAX ||
-     MaxPt.x() == -DBL_MAX || MaxPt.y() == -DBL_MAX || MaxPt.z() == -DBL_MAX) 
-    return true;
-  return false;
-}
-
-void SBoundingBox3d::operator+=(const SPoint3 &pt)
-// note: it is possible for pt[i] to be both > MaxPt[i] and < MinPt[i]
-// the first point always will be both
-{
-  if(pt[0] < MinPt[0])
-    MinPt[0] = pt[0];
-  if (pt[0] > MaxPt[0])
-    MaxPt[0] = pt[0];
-
-  if(pt[1] < MinPt[1])
-    MinPt[1] = pt[1];
-  if (pt[1] > MaxPt[1])
-    MaxPt[1] = pt[1];
-
-  if(pt[2] < MinPt[2])
-    MinPt[2] = pt[2];
-  if (pt[2] > MaxPt[2])
-    MaxPt[2] = pt[2];
-}
-
-void SBoundingBox3d::operator+=(const SBoundingBox3d &box)
-{
-  (*this)+=box.MinPt;
-  (*this)+=box.MaxPt;
-}
-
-SPoint3 SBoundingBox3d::min() const
-{ return MinPt; }
-
-SPoint3 SBoundingBox3d::max() const
-{ return MaxPt; }
-
-SPoint3 SBoundingBox3d::center() const
-{ return (MinPt+MaxPt)*.5; }
-
-void SBoundingBox3d::operator*=(double scale)
-{
-  SPoint3 center = (MinPt + MaxPt)*.5;
-  MaxPt -= center;
-  MinPt -= center;
-  MaxPt *= scale;
-  MinPt *= scale;
-  MaxPt += center;
-  MinPt += center;
-}
-
-void SBoundingBox3d::scale(double sx, double sy, double sz)
-{
-  SPoint3 center = (MinPt + MaxPt)*.5;
-  MaxPt -= center;
-  MinPt -= center;
-  MaxPt[0] *= sx;   MaxPt[1] *= sy; MaxPt[2] *= sz;
-  MinPt[0] *= sx;   MinPt[1] *= sy; MinPt[2] *= sz;
-  MaxPt += center;
-  MinPt += center;
-}
-
-void SBoundingBox3d::makeCube()
-{
-  SPoint3 len = MaxPt-MinPt;
-  double scales[3];
-  double max=-1.0;
-
-  for(int i = 0; i < 3; i++)
-    max = len[i] > max ? len[i] : max;
-
-  for(int j = 0; j < 3; j++)
-    scales[j] = max/len[j];
-
-  scale(scales[0],scales[1],scales[2]);
-}
diff --git a/Geo/SBoundingBox3d.h b/Geo/SBoundingBox3d.h
index 462c7bcc9ce8c41f7e19efd9ed58d49d0417c1e8..003d24c9fef4edc90c556f40f9b273d56657366c 100644
--- a/Geo/SBoundingBox3d.h
+++ b/Geo/SBoundingBox3d.h
@@ -20,26 +20,89 @@
 // 
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
+#include <float.h>
 #include "SPoint3.h"
 
 // A bounding box class - add points and it grows to be the bounding
 // box of the point set
 class SBoundingBox3d {
  public:
-  SBoundingBox3d();
-  SBoundingBox3d(const SPoint3 &);
-  bool empty();
-  void reset();
-  void operator+=(const SPoint3 &pt);
-  void operator+=(const SBoundingBox3d &pt);
-  void operator*=(double scale);
-  void scale(double, double, double);
-  SPoint3 min() const;
-  SPoint3 max() const;
-  SPoint3 center() const;
-  void makeCube();
+  SBoundingBox3d() 
+    : MinPt(DBL_MAX,DBL_MAX,DBL_MAX), MaxPt(-DBL_MAX,-DBL_MAX,-DBL_MAX) {}
+  SBoundingBox3d(const SPoint3 &pt)
+    : MinPt(pt), MaxPt(pt) {}
+  bool empty()
+  {
+    if(MinPt.x() == DBL_MAX || MinPt.y() == DBL_MAX || MinPt.z() == DBL_MAX ||
+       MaxPt.x() == -DBL_MAX || MaxPt.y() == -DBL_MAX || MaxPt.z() == -DBL_MAX) 
+      return true;
+    return false;
+  }
+  void reset()
+  {
+    MinPt = SPoint3(DBL_MAX,DBL_MAX,DBL_MAX);
+    MaxPt = SPoint3(-DBL_MAX,-DBL_MAX,-DBL_MAX);
+  }
+  void operator+=(const SPoint3 &pt)
+  {
+    // note: it is possible for pt[i] to be both > MaxPt[i] and < MinPt[i]
+    // the first point always will be both
+    if(pt[0] < MinPt[0])
+      MinPt[0] = pt[0];
+    if (pt[0] > MaxPt[0])
+      MaxPt[0] = pt[0];
+    
+    if(pt[1] < MinPt[1])
+      MinPt[1] = pt[1];
+    if (pt[1] > MaxPt[1])
+      MaxPt[1] = pt[1];
+    
+    if(pt[2] < MinPt[2])
+      MinPt[2] = pt[2];
+    if (pt[2] > MaxPt[2])
+      MaxPt[2] = pt[2];
+  }
+  void operator+=(const SBoundingBox3d &box)
+  {
+    (*this) += box.MinPt;
+    (*this) += box.MaxPt;
+  }
+  void operator*=(double scale)
+  {
+    SPoint3 center = (MinPt + MaxPt) * .5;
+    MaxPt -= center;
+    MinPt -= center;
+    MaxPt *= scale;
+    MinPt *= scale;
+    MaxPt += center;
+    MinPt += center;
+  }
+  void scale(double sx, double sy, double sz)
+  {
+    SPoint3 center = (MinPt + MaxPt) * .5;
+    MaxPt -= center;
+    MinPt -= center;
+    MaxPt[0] *= sx; MaxPt[1] *= sy; MaxPt[2] *= sz;
+    MinPt[0] *= sx; MinPt[1] *= sy; MinPt[2] *= sz;
+    MaxPt += center;
+    MinPt += center;
+  }
+  SPoint3 min() const { return MinPt; }
+  SPoint3 max() const { return MaxPt; }
+  SPoint3 center() const { return (MinPt + MaxPt) * .5; }
+  void makeCube()
+  {
+    SPoint3 len = MaxPt - MinPt;
+    double scales[3];
+    double max = -1.0;
+    for(int i = 0; i < 3; i++)
+      max = len[i] > max ? len[i] : max;
+    for(int j = 0; j < 3; j++)
+      scales[j] = max/len[j];
+    scale(scales[0],scales[1],scales[2]);
+  }
  private:
-  SPoint3 MinPt,MaxPt;
+  SPoint3 MinPt, MaxPt;
 };
 
 #endif
diff --git a/Geo/SPoint3.h b/Geo/SPoint3.h
index 8de312ed774eb0163c18a995b419730372121555..deaeb02437a26f173b2279cc5f929078db35507e 100644
--- a/Geo/SPoint3.h
+++ b/Geo/SPoint3.h
@@ -28,9 +28,9 @@ class SPoint3 {
   double P[3];
  public:
   SPoint3() {}
-  SPoint3(double x, double y, double z) {P[0] = x; P[1] = y; P[2] = z;}
-  SPoint3(const double *p) {P[0] = p[0]; P[1] = p[1]; P[2] = p[2];}
-  SPoint3(const SPoint3 &pt) {P[0] = pt.P[0]; P[1] = pt.P[1]; P[2] = pt.P[2]; }
+  SPoint3(double x, double y, double z) { P[0] = x; P[1] = y; P[2] = z; }
+  SPoint3(const double *p) { P[0] = p[0]; P[1] = p[1]; P[2] = p[2]; }
+  SPoint3(const SPoint3 &pt) { P[0] = pt.P[0]; P[1] = pt.P[1]; P[2] = pt.P[2]; }
   virtual ~SPoint3() {}
   void setPosition(double xx, double yy, double zz);
   void getPosition(double *xx, double *yy, double *zz) const;
@@ -49,16 +49,16 @@ class SPoint3 {
 };
 
 inline SPoint3 operator + (const SPoint3 &a, const SPoint3 &b)
-{ return SPoint3(a.x()+b.x(),a.y()+b.y(),a.z()+b.z()); }
+{ return SPoint3(a.x() + b.x(), a.y() + b.y(), a.z() + b.z()); }
 
 inline SPoint3 operator - (const SPoint3 &a, const SPoint3 &b)
-{ return SPoint3(a.x()-b.x(),a.y()-b.y(),a.z()-b.z()); }
+{ return SPoint3(a.x() - b.x(), a.y() - b.y(), a.z() - b.z()); }
 
 inline void SPoint3::setPosition(double xx, double yy, double zz)
-{ P[0] = xx;  P[1] = yy;  P[2] = zz; }
+{ P[0] = xx; P[1] = yy; P[2] = zz; }
 
 inline void SPoint3::getPosition(double *xx, double *yy, double *zz) const
-{ *xx = P[0];  *yy = P[1];  *zz = P[2]; }
+{ *xx = P[0]; *yy = P[1]; *zz = P[2]; }
 
 inline void SPoint3::position(double *p) const
 { p[0] = P[0]; p[1] = P[1]; p[2] = P[2]; }
@@ -73,19 +73,19 @@ inline double SPoint3::z(void) const
 { return P[2]; }
 
 inline SPoint3 & SPoint3::operator=(const SPoint3 &p)
-{ P[0] = p.P[0]; P[1]=p.P[1]; P[2]=p.P[2]; return *this; }
+{ P[0] = p.P[0]; P[1] = p.P[1]; P[2] = p.P[2]; return *this; }
 
 inline void SPoint3::operator+=(const SPoint3 &p)
-{ P[0] += p.P[0]; P[1] += p.P[1]; P[2] += p.P[2];}
+{ P[0] += p.P[0]; P[1] += p.P[1]; P[2] += p.P[2]; }
 
 inline void SPoint3::operator-=(const SPoint3 &p)
-{ P[0] -= p.P[0]; P[1] -= p.P[1]; P[2] -= p.P[2];}
+{ P[0] -= p.P[0]; P[1] -= p.P[1]; P[2] -= p.P[2]; }
 
 inline void SPoint3::operator*=(double mult)
 { P[0] *= mult; P[1] *= mult; P[2] *= mult; }
 
 inline SPoint3 SPoint3::operator*(double mult)
-{ return SPoint3(P[0]*mult, P[1]*mult, P[2] *= mult); }
+{ return SPoint3(P[0] * mult, P[1] * mult, P[2] * mult); }
 
 inline double &SPoint3::operator[](int i)
 { return P[i]; }
diff --git a/Geo/SVector3.cpp b/Geo/SVector3.cpp
deleted file mode 100644
index b2c8846fdbffe404ec52f1300833655ebbe16c79..0000000000000000000000000000000000000000
--- a/Geo/SVector3.cpp
+++ /dev/null
@@ -1,122 +0,0 @@
-// $Id: SVector3.cpp,v 1.5 2006-11-27 22:22:14 geuzaine Exp $
-//
-// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
-//
-// This program is free software; you can redistribute it and/or modify
-// it under the terms of the GNU General Public License as published by
-// the Free Software Foundation; either version 2 of the License, or
-// (at your option) any later version.
-//
-// This program is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with this program; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
-// USA.
-// 
-// Please report all bugs and problems to <gmsh@geuz.org>.
-
-#include "SVector3.h"
-#include "math.h"
-
-
-SVector3 crossprod(const SVector3 &a, const SVector3 &b)
-{
-  return SVector3(a.y()*b.z()-b.y()*a.z(),
-		 -(a.x()*b.z()-b.x()*a.z()),
-		 a.x()*b.y()-b.x()*a.y());
-}
-
-double dot(const SVector3 &a, const SVector3 &b)
-{
-  return( a.x()*b.x()+a.y()*b.y()+a.z()*b.z() );
-}
-
-double angle(const SVector3 &a, const SVector3 &b, const SVector3 &n)
-{
-  double ddot = dot(a,b)/(norm(a)*norm(b));
-  if(ddot>1.0)
-    ddot = 1.0;
-  if(ddot<-1.0)
-    ddot = -1.0;
-  double dang = acos(ddot);
-  // check if we just found the angle or 2*PI-angle
-  SVector3 check = crossprod(a,b);
-  // ***** check this below ******
-  if( norm(check) > 0.0000001*(norm(a)+norm(b))){  // check if a,b are colinear
-    double dir = dot(check,n);
-    dang = dir < 0 ? 2*M_PI-dang : dang;
-  } else {
-    if( ddot > 0)
-      dang = 2*M_PI;
-  }
-  return dang;
-
-}
-
-double norm(const SVector3 &v)
-{ return sqrt(dot(v,v)); }
-
-double SVector3::normalize()
-{ 
-  double n = norm(*this);
-  if(n){
-    P[0] /= n; P[1]/= n; P[2] /= n;
-  }
-  return n;
-}
-
-SVector3 operator*(double m,const SVector3 &v)
-{
-  return SVector3(v[0]*m,v[1]*m,v[2]*m);
-}
-
-SVector3 operator*(const SVector3 &v, double m)
-{
-  return SVector3(v[0]*m,v[1]*m,v[2]*m);
-}
-
-SVector3 operator*(const SVector3 &v1, const SVector3 &v2)
-{
-  return SVector3(v1[0]*v2[0],v1[1]*v2[1],v1[2]*v2[2]);
-}
-
-SVector3 operator+(const SVector3 &a,const SVector3 &b)
-{
-  return SVector3(a[0]+b[0],a[1]+b[1],a[2]+b[2]);
-}
-
-SVector3 operator-(const SVector3 &a,const SVector3 &b)
-{
-  return SVector3(a[0]-b[0],a[1]-b[1],a[2]-b[2]);
-}
-
-SVector3 & SVector3::operator += (const SVector3 &a)
-{ 
-  for(int i = 0; i < 3; i++)
-    P[i] += a[i];
-  return *this;
-}
-
-SVector3 & SVector3::operator -= (const SVector3 &a)
-{ 
-  for(int i = 0; i < 3; i++)
-    P[i] -= a[i];
-  return *this;
-}
-
-SVector3 & SVector3::operator *= (const SVector3 &a)
-{ 
-  for(int i = 0; i < 3; i++)
-    P[i] *= a[i];
-  return *this;
-}
-
-SVector3 & SVector3::operator = (double v)
-{
-  P[0] = v; P[1] = v; P[2] = v;
-  return *this;
-}
diff --git a/Geo/SVector3.h b/Geo/SVector3.h
index e632026ff3def17e69b1b32496d82584e2c7760c..6bf6ecec7b2f8df97b2b8b661c56c05464e2f52b 100644
--- a/Geo/SVector3.h
+++ b/Geo/SVector3.h
@@ -27,50 +27,78 @@ class SVector3 {
  public:
   SVector3() {}
   // Construct from 2 SPoints, vector from p1 to p2
-  SVector3(const SPoint3 &p1, const SPoint3 &p2) : P(p2-p1) {}
+  SVector3(const SPoint3 &p1, const SPoint3 &p2) : P(p2 - p1) {}
   // Construct from a single SPoint, vector from origin to p1
   SVector3(const SPoint3 &p1) : P(p1) {}
-  SVector3(double x, double y, double z) : P(x,y,z) {}
-  SVector3(double v) : P(v,v,v) {}
+  SVector3(double x, double y, double z) : P(x, y, z) {}
+  SVector3(double v) : P(v, v, v) {}
   SVector3(const double *array) : P(array) {}
   inline double x(void) const { return P.x(); }
   inline double y(void) const { return P.y(); }
   inline double z(void) const { return P.z(); }
-  double normalize();
+  inline double norm() { return sqrt(P[0] * P[0] + P[1] * P[1] + P[2] * P[2]); }
+  double normalize() 
+  { 
+    double n = norm(); if(n){ P[0] /= n; P[1] /= n; P[2] /= n; }
+    return n;
+  }
   // why both [] and (), why not
-  double &operator[](int);
-  double operator[](int) const;
-  double &operator()(int);
-  double operator()(int) const;
-  SVector3 & operator += (const SVector3 &);
-  SVector3 & operator -= (const SVector3 &);
-  SVector3 & operator *= (const SVector3 &);
-  SVector3 & operator = (double);
+  double &operator[](int i){ return P[i]; }
+  double operator[](int i) const { return P[i]; }
+  double &operator()(int i){ return P[i]; }
+  double operator()(int i) const { return P[i]; }
+  SVector3 & operator += (const SVector3 &a)
+  {
+    for(int i = 0; i < 3; i++)
+      P[i] += a[i];
+    return *this;
+  }
+  SVector3 & operator -= (const SVector3 &a)
+  { 
+    for(int i = 0; i < 3; i++)
+      P[i] -= a[i];
+    return *this;
+  }
+  SVector3 & operator *= (const SVector3 &a)
+  { 
+    for(int i = 0; i < 3; i++)
+      P[i] *= a[i];
+    return *this;
+  }
+  SVector3 & operator = (double v)
+  {
+    P[0] = v; P[1] = v; P[2] = v;
+    return *this;
+  }
   operator double *() { return P; }
  protected:
   SPoint3 P;
 };
 
-double norm(const SVector3 &v);
-double dot(const SVector3 &a, const SVector3 &b);
-SVector3 crossprod(const SVector3 &a, const SVector3 &b);
-double angle(const SVector3 &a, const SVector3 &b, const SVector3 &n);
-SVector3 operator*(double m,const SVector3 &v);
-SVector3 operator*(const SVector3 &v, double m);
-SVector3 operator*(const SVector3 &v1, const SVector3 &v2);
-SVector3 operator+(const SVector3 &a,const SVector3 &b);
-SVector3 operator-(const SVector3 &a,const SVector3 &b);
+inline double dot(const SVector3 &a, const SVector3 &b)
+{ return a.x() * b.x() + a.y() * b.y() + a.z() * b.z(); }
 
-inline double &SVector3::operator[](int i)
-{ return P[i]; }
+inline double norm(const SVector3 &v)
+{ return sqrt(dot(v, v)); }
 
-inline double SVector3::operator[](int i) const
-{ return P[i]; }
+inline SVector3 crossprod(const SVector3 &a, const SVector3 &b)
+{ return SVector3(a.y() * b.z() - b.y() * a.z(), 
+		  -(a.x() * b.z() - b.x() * a.z()), 
+		  a.x() * b.y() - b.x() * a.y()); }
 
-inline double &SVector3::operator()(int i)
-{ return P[i]; }
+inline SVector3 operator*(double m,const SVector3 &v)
+{ return SVector3(v[0] * m, v[1] * m, v[2] * m); }
 
-inline double SVector3::operator()(int i) const
-{ return P[i]; }
+inline SVector3 operator*(const SVector3 &v, double m)
+{ return SVector3(v[0] * m, v[1] * m, v[2] * m); }
+
+inline SVector3 operator*(const SVector3 &v1, const SVector3 &v2)
+{ return SVector3(v1[0] * v2[0], v1[1] * v2[1], v1[2] * v2[2]); }
+
+inline SVector3 operator+(const SVector3 &a,const SVector3 &b)
+{ return SVector3(a[0] + b[0], a[1] + b[1], a[2] + b[2]); }
+
+inline SVector3 operator-(const SVector3 &a,const SVector3 &b)
+{ return SVector3(a[0] - b[0], a[1] - b[1], a[2] - b[2]); }
 
 #endif