diff --git a/source/graphics/Model.cpp b/source/graphics/Model.cpp
index 2e10820f94..6832af9535 100644
--- a/source/graphics/Model.cpp
+++ b/source/graphics/Model.cpp
@@ -325,7 +325,7 @@ void CModel::ValidatePosition()
 	
 		debug_assert(m_pModelDef->GetNumBones() == m_Anim->m_AnimDef->GetNumKeys());
 	
-		m_Anim->m_AnimDef->BuildBoneMatrices(m_AnimTime,m_BoneMatrices);
+		m_Anim->m_AnimDef->BuildBoneMatrices(m_AnimTime, m_BoneMatrices, !(m_Flags & MODELFLAG_NOLOOPANIMATION));
 	
 		const CMatrix3D& transform=GetTransform();
 		for (size_t i=0;i<m_pModelDef->GetNumBones();i++) {
diff --git a/source/graphics/SkeletonAnimDef.cpp b/source/graphics/SkeletonAnimDef.cpp
index 41507e53a9..f903c67f45 100644
--- a/source/graphics/SkeletonAnimDef.cpp
+++ b/source/graphics/SkeletonAnimDef.cpp
@@ -30,7 +30,7 @@ CSkeletonAnimDef::~CSkeletonAnimDef()
 ///////////////////////////////////////////////////////////////////////////////////////////
 // BuildBoneMatrices: build matrices for all bones at the given time (in MS) in this 
 // animation
-void CSkeletonAnimDef::BuildBoneMatrices(float time,CMatrix3D* matrices) const
+void CSkeletonAnimDef::BuildBoneMatrices(float time, CMatrix3D* matrices, bool loop) const
 {
 	float fstartframe = time/m_FrameTime;
 	u32 startframe = u32(time/m_FrameTime);
@@ -41,13 +41,12 @@ void CSkeletonAnimDef::BuildBoneMatrices(float time,CMatrix3D* matrices) const
 	u32 endframe = startframe + 1;
 	endframe %= m_NumFrames; 
 
-	if (endframe == 0)
+	if (!loop && endframe == 0)
 	{
 		// This might be something like a death animation, and interpolating
 		// between the final frame and the initial frame is wrong, because they're
 		// totally different. So if we've looped around to endframe==0, just display
 		// the animation's final frame with no interpolation.
-		// (TODO: this is only sometimes valid - how can we tell the difference?)
 		for (u32 i = 0; i < m_NumKeys; i++)
 		{
 			const Key& key = GetKey(startframe, i);
@@ -60,7 +59,7 @@ void CSkeletonAnimDef::BuildBoneMatrices(float time,CMatrix3D* matrices) const
 	{
 		for (u32 i = 0; i < m_NumKeys; i++)
 		{
-			const Key& startkey=  GetKey(startframe, i);
+			const Key& startkey = GetKey(startframe, i);
 			const Key& endkey = GetKey(endframe, i);
 
 			CVector3D trans = Interpolate(startkey.m_Translation, endkey.m_Translation, deltatime);
@@ -68,8 +67,7 @@ void CSkeletonAnimDef::BuildBoneMatrices(float time,CMatrix3D* matrices) const
 			CQuaternion rot;
 			rot.Slerp(startkey.m_Rotation, endkey.m_Rotation, deltatime);
 
-			matrices[i].SetIdentity();
-			matrices[i].Rotate(rot);
+			rot.ToMatrix(matrices[i]);
 			matrices[i].Translate(trans);
 		}
 	}
diff --git a/source/graphics/SkeletonAnimDef.h b/source/graphics/SkeletonAnimDef.h
index 8d29cabbdc..997d18d506 100644
--- a/source/graphics/SkeletonAnimDef.h
+++ b/source/graphics/SkeletonAnimDef.h
@@ -61,7 +61,7 @@ public:
 	size_t GetNumFrames() const { return (size_t)m_NumFrames; }
 
 	// build matrices for all bones at the given time (in MS) in this animation
-	void BuildBoneMatrices(float time, CMatrix3D* matrices) const;
+	void BuildBoneMatrices(float time, CMatrix3D* matrices, bool loop) const;
 
 	// anim I/O functions
 	static CSkeletonAnimDef* Load(const char* filename);
diff --git a/source/graphics/UnitAnimation.cpp b/source/graphics/UnitAnimation.cpp
index 4575e7e57a..c040b8563b 100644
--- a/source/graphics/UnitAnimation.cpp
+++ b/source/graphics/UnitAnimation.cpp
@@ -81,6 +81,9 @@ void CUnitAnimation::Update(float time)
 	if (m_Looping && model->NeedsNewAnim(time))
 	{
 		m_Unit.SetRandomAnimation(m_State, !m_Looping, DesyncSpeed(m_Speed));
+		// TODO: this really ought to transition smoothly into the new animation,
+		// instead of just cutting off the end of the previous one and jumping
+		// straight into the new.
 	}
 
 	if (m_TimeToNextSync >= 0.0 && m_TimeToNextSync-time < 0.0)
diff --git a/source/maths/Matrix3D.cpp b/source/maths/Matrix3D.cpp
index eff279d52d..0147f8ef5b 100644
--- a/source/maths/Matrix3D.cpp
+++ b/source/maths/Matrix3D.cpp
@@ -218,14 +218,16 @@ void CMatrix3D::SetTranslation(const CVector3D& vector)
 //Applies a translation to the matrix
 void CMatrix3D::Translate(float x, float y, float z)
 {
-	CMatrix3D Temp;
-	Temp.SetTranslation(x,y,z);
-	Concatenate(Temp);
+	_14 += x;
+	_24 += y;
+	_34 += z;
 }
 
 void CMatrix3D::Translate(const CVector3D &vector)
 {
-	Translate(vector.X,vector.Y,vector.Z);
+	_14 += vector.X;
+	_24 += vector.Y;
+	_34 += vector.Z;
 }
 
 void CMatrix3D::Concatenate(const CMatrix3D& m)
diff --git a/source/maths/Matrix3D.h b/source/maths/Matrix3D.h
index 5b9e197fe1..5af25a2c1f 100644
--- a/source/maths/Matrix3D.h
+++ b/source/maths/Matrix3D.h
@@ -82,7 +82,9 @@ public:
 	void SetTranslation(float x, float y, float z);
 	void SetTranslation(const CVector3D& vector);
 
-	// concatenate given translation onto this matrix
+	// concatenate given translation onto this matrix. Assumes the current
+	// matrix is an affine transformation (i.e. the bottom row is [0,0,0,1])
+	// as an optimisation.
 	void Translate(float x, float y, float z);
 	void Translate(const CVector3D& vector);
 
diff --git a/source/maths/Quaternion.cpp b/source/maths/Quaternion.cpp
index 3024e69e03..a7f06e23b2 100644
--- a/source/maths/Quaternion.cpp
+++ b/source/maths/Quaternion.cpp
@@ -25,45 +25,65 @@ CQuaternion::CQuaternion(float x, float y, float z, float w)
 {
 }
 
-//quaternion addition
 CQuaternion CQuaternion::operator + (const CQuaternion &quat) const
 {
 	CQuaternion Temp;
-
 	Temp.m_W = m_W + quat.m_W;
 	Temp.m_V = m_V + quat.m_V;
-
 	return Temp;
 }
 
-//quaternion addition/assignment
 CQuaternion &CQuaternion::operator += (const CQuaternion &quat)
 {
-	m_W += quat.m_W;
-	m_V += quat.m_V;
-
-	return (*this);
+	*this = *this + quat;
+	return *this;
+}
+
+CQuaternion CQuaternion::operator - (const CQuaternion &quat) const
+{
+	CQuaternion Temp;
+	Temp.m_W = m_W - quat.m_W;
+	Temp.m_V = m_V - quat.m_V;
+	return Temp;
+}
+
+CQuaternion &CQuaternion::operator -= (const CQuaternion &quat)
+{
+	*this = *this - quat;
+	return *this;
 }
 
-//quaternion multiplication
 CQuaternion CQuaternion::operator * (const CQuaternion &quat) const
 {
 	CQuaternion Temp;
-
 	Temp.m_W = (m_W * quat.m_W) - (m_V.Dot(quat.m_V));
 	Temp.m_V = (m_V.Cross(quat.m_V)) + (quat.m_V * m_W) + (m_V * quat.m_W);
-
 	return Temp;
 }
 
-//quaternion multiplication/assignment
 CQuaternion &CQuaternion::operator *= (const CQuaternion &quat)
 {
-	(*this) = (*this) * quat;
-
-	return (*this);
+	*this = *this * quat;
+	return *this;
 }
 
+CQuaternion CQuaternion::operator * (float factor) const
+{
+	CQuaternion Temp;
+	Temp.m_W = m_W * factor;
+	Temp.m_V = m_V * factor;
+	return Temp;
+}
+
+
+float CQuaternion::Dot(const CQuaternion& quat) const
+{
+	return
+		m_V.X * quat.m_V.X +
+		m_V.Y * quat.m_V.Y +
+		m_V.Z * quat.m_V.Z +
+		m_W   * quat.m_W;
+}
 
 void CQuaternion::FromEulerAngles (float x, float y, float z)
 {
@@ -167,16 +187,13 @@ void CQuaternion::ToMatrix(CMatrix3D& result) const
 	result._44 = 1;
 }
 
-void CQuaternion::Slerp(const CQuaternion& from,const CQuaternion& to, float ratio)
+void CQuaternion::Slerp(const CQuaternion& from, const CQuaternion& to, float ratio)
 {
 	float to1[4];
 	float omega, cosom, sinom, scale0, scale1;
 	
 	// calc cosine
-	cosom = from.m_V.X * to.m_V.X  +
-			from.m_V.Y * to.m_V.Y  +
-			from.m_V.Z * to.m_V.Z  +
-			from.m_W * to.m_W;
+	cosom = from.Dot(to);
 
 
 	// adjust signs (if necessary)
@@ -220,6 +237,16 @@ void CQuaternion::Slerp(const CQuaternion& from,const CQuaternion& to, float rat
 	m_W	  = scale0 * from.m_W + scale1 * to1[3];
 }
 
+void CQuaternion::Nlerp(const CQuaternion& from, const CQuaternion& to, float ratio)
+{
+	float c = from.Dot(to);
+	if (c < 0.f)
+		*this = from - (to + from) * ratio;
+	else
+		*this = from + (to - from) * ratio;
+	Normalize();
+}
+
 ///////////////////////////////////////////////////////////////////////////////////////////////
 // FromAxisAngle: create a quaternion from axis/angle representation of a rotation
 void CQuaternion::FromAxisAngle(const CVector3D& axis, float angle)
diff --git a/source/maths/Quaternion.h b/source/maths/Quaternion.h
index 4a09667ea1..d35d7b6ec0 100644
--- a/source/maths/Quaternion.h
+++ b/source/maths/Quaternion.h
@@ -22,15 +22,18 @@ public:
 	CQuaternion();
 	CQuaternion(float x, float y, float z, float w);
 	
-	// Quaternion addition
 	CQuaternion operator + (const CQuaternion &quat) const;
-	// Quaternion addition/assignment
 	CQuaternion &operator += (const CQuaternion &quat);
 
-	// Quaternion multiplication
+	CQuaternion operator - (const CQuaternion &quat) const;
+	CQuaternion &operator -= (const CQuaternion &quat);
+
 	CQuaternion operator * (const CQuaternion &quat) const;
-	// Quaternion multiplication/assignment
 	CQuaternion &operator *= (const CQuaternion &quat);
+
+	CQuaternion operator * (float factor) const;
+
+	float Dot(const CQuaternion& quat) const;
 	
 	void FromEulerAngles (float x, float y, float z);
 	CVector3D ToEulerAngles();
@@ -42,6 +45,9 @@ public:
 	// Sphere interpolation
 	void Slerp(const CQuaternion& from, const CQuaternion& to, float ratio);
 
+	// Normalised linear interpolation
+	void Nlerp(const CQuaternion& from, const CQuaternion& to, float ratio);
+
 	// Create a quaternion from axis/angle representation of a rotation
 	void FromAxisAngle(const CVector3D& axis, float angle);