diff --git a/vecmat3.cc b/vecmat3.cc
index 6357666..e0f0e0d 100644
--- a/vecmat3.cc
+++ b/vecmat3.cc
@@ -69,6 +69,23 @@ m[2][2]=q[2]*rhs.q[2];
 return m;
 }
 
+template<typename T>
+void Vec3<T>::addouter(Mat3<T> &m, const Vec3<T> &rhs,  const T weight) const
+{
+m[0][0]+=weight* q[0]*rhs.q[0];
+m[0][1]+=weight* q[0]*rhs.q[1];
+m[0][2]+=weight* q[0]*rhs.q[2];
+
+m[1][0]+=weight* q[1]*rhs.q[0];
+m[1][1]+=weight* q[1]*rhs.q[1];
+m[1][2]+=weight* q[1]*rhs.q[2];
+
+m[2][0]+=weight* q[2]*rhs.q[0];
+m[2][1]+=weight* q[2]*rhs.q[1];
+m[2][2]+=weight* q[2]*rhs.q[2];
+
+}
+
 
 template<typename T>
 void Vec3<T>::inertia(Mat3<T> &m, const T weight) const
diff --git a/vecmat3.h b/vecmat3.h
index b1b2eda..44670b4 100644
--- a/vecmat3.h
+++ b/vecmat3.h
@@ -96,6 +96,7 @@ public:
 	const Vec3 operator*(const Mat3<T> &rhs) const;
 	const Vec3 timesT(const Mat3<T> &rhs) const; //with transpose
 	Mat3<T> outer(const Vec3 &rhs) const; //tensor product
+	void addouter(Mat3<T> &m, const Vec3 &rhs, const T weight) const; //tensor product
 	void inertia(Mat3<T> &itensor, const T weight) const; //contribution to inertia tensor
 	void randomize(const T x);