Figure 2: A partial listing of the tnVector class

template<int Dim> class tnVector {
public:
     
   /* Note:
    *   all functions are considered inline according 
    *   to C++ specifications.
    */
   tnVector() {}
   tnVector(const double x, const double y) { setValue(x, y); }
   tnVector(const double x, const double y, const double z)
      { setValue(x,y,z); }
     
   double operator[](const int ix) const { return _v[ix]; } 
   double& operator[](const int ix) { return _v[ix]; }
      
   double lengthSquared() const { return dot(*this); }
     
   double length() const { return sqrt(dot(*this)) }; 
   double distanceSquared(const tnVector<Dim>& p) const {
      double t, a(0.0);
      for (int ix = 0; ix < Dim; ++ix)
         { t = _v[ix] - p[ix]; a += t * t; }
      return a;
   }
     
     
   void increment(const tnVector<Dim>& x, const double s = 1.0)
      { for (int ix = 0; ix<Dim; ++ix) _v[ix] += x._v[ix] * s; }
     
   // Math operators
   void operator+=(const tnVector<Dim>& p)
      { for (int ix = 0; ix < Dim; ++ix) _v[ix] += p._v[ix]; }
     
   void operator-=(const tnVector<Dim>& p)
      { for (int ix = 0; ix < Dim; ++ix) _v[ix] -= p._v[ix]; }
     
   void operator+=(const double s)
      { for (int ix=0; ix<Dim; ++ix) _v[ix] += s;}
   void operator-=(const double s)
      { for (int ix=0; ix<Dim; ++ix) _v[ix] -= s;}
   void operator*=(const double s)
      { for (int ix=0; ix<Dim; ++ix) _v[ix] *= s;}
   void operator/=(const double s)
      { for (int ix=0; ix<Dim; ++ix) _v[ix] /= s;}
     
   void setDifference(const tnVector<Dim>& x,
                      const tnVector<Dim>& y) {
     for (int ix=0; ix<Dim; ++ix)
        _v[ix] = x._v[ix] - y._v[ix];
   }
     
   void setSum(const tnVector<Dim>& x, const tnVector<Dim>& y) {
      for (int ix=0; ix<Dim; ++ix) _v[ix] = x._v[ix] + y._v[ix];
   }
     
   // Set to weighted sum of vectors: self = x * w_x + y * w_y. 
   void setWeightedSum(const tnVector<Dim>& x,
                       const double w_x,
                       const tnVector<Dim>& y,
                       const double w_y = 1.0) {
      for(int ix=0; ix<Dim; ++ix)
         _v[ix] = x._v[ix]*w_x + y._v[ix] * w_y; }
     
   // Dimension-specific functions
   void setValue(const double x, const double y); 
   double directedAngle(const tnVector<2>&) const;
     
   void setValue(const double, const double, const double); 
   void setCross(const tnVector<3>& x, const tnVector<3>& y); 
   double tripleProduct(const tnVector<3>&,
                        const tnVector<3>&) const;
     
protected:
   double _v[Dim]; // Storage for vector elements.
     
};  // End tnVector
     
// Inline definitions of dimension-specific functions
     
inline void tnVector<2>
::setValue(const double x, const double y) 
{ _v[0] = x; _v[1] = y; }
     
inline void tnVector<3>
::setValue(const double x, const double y, const double z) 
{ _v[0] = x; _v[1] = y; _v[2] = z; }
     
inline void tnVector<3>
::setCross(const tnVector<3>& x, const tnVector<3>& y) {
   _v[0] = x._v[1] * y._v[2] - x._v[2] * y._v[1]; 
   _v[1] = x._v[2] * y._v[0] - x._v[0] * y._v[2]; 
   _v[2] = x._v[0] * y._v[1] - x._v[1] * y._v[0];
}
     
inline double tnVector<2>
::directedAngle(const tnVector<2>& b) const {
   return atan2(_v[0]*b._v[1] - _v[1]*b._v[0],
                _v[0]*b._v[0] + _v[1]*b._v[1]); }
     
inline double tnVector<3>
::tripleProduct(const tnVector<3>& b,
                const tnVector<3>& c) const {
   return
      _v[0] * (b._v[1] * c._v[2] - b._v[2] * c._v[1]) + 
      _v[1] * (b._v[2] * c._v[0] - b._v[0] * c._v[2]) + 
      _v[2] * (b._v[0] * c._v[1] - b._v[1] * c._v[0]);
}