ENH: quaternion: Add power and exp functions

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/OpenFOAM/primitives/quaternion/quaternion.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,21 +58,60 @@ Foam::quaternion Foam::slerp
     const scalar t
 )
 {
-    // Calculate angle between the quaternions
+    return qa*pow((inv(qa)*qb), t);
-    scalar cosHalfTheta = qa & qb;
+}
 
-    if (mag(cosHalfTheta) >= 1)
+
+Foam::quaternion Foam::exp(const quaternion& q)
+{
+    const scalar magV = mag(q.v());
+
+    if (magV == 0)
     {
-        return qa;
+        return quaternion(1, vector::zero);
     }
 
-    scalar halfTheta = acos(cosHalfTheta);
+    const scalar expW = exp(q.w());
-    scalar sinHalfTheta = sqrt(1.0 - sqr(cosHalfTheta));
 
-    scalar wa = sin((1 - t)*halfTheta)/sinHalfTheta;
+    return quaternion
-    scalar wb = sin(t*halfTheta)/sinHalfTheta;
+    (
+        expW*cos(magV),
+        expW*sin(magV)*q.v()/magV
+    );
+}
 
-    return wa*qa + wb*qb;
+
+Foam::quaternion Foam::pow(const quaternion& q, const label power)
+{
+    const scalar magQ = mag(q);
+    const scalar magV = mag(q.v());
+
+    quaternion powq(q.v());
+
+    if (magV != 0 && magQ != 0)
+    {
+        powq /= magV;
+        powq *= power*acos(q.w()/magQ);
+    }
+
+    return pow(magQ, power)*exp(powq);
+}
+
+
+Foam::quaternion Foam::pow(const quaternion& q, const scalar power)
+{
+    const scalar magQ = mag(q);
+    const scalar magV = mag(q.v());
+
+    quaternion powq(q.v());
+
+    if (magV != 0 && magQ != 0)
+    {
+        powq /= magV;
+        powq *= power*acos(q.w()/magQ);
+    }
+
+    return pow(magQ, power)*exp(powq);
 }
 
 

src/OpenFOAM/primitives/quaternion/quaternion.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -226,6 +226,15 @@ quaternion slerp
     const scalar t
 );
 
+//- Exponent of a quaternion
+quaternion exp(const quaternion& q);
+
+//- Power of a quaternion
+quaternion pow(const quaternion& q, const label power);
+
+//- Power of a quaternion
+quaternion pow(const quaternion& q, const scalar power);
+
 //- Data associated with quaternion type are contiguous
 template<>
 inline bool contiguous<quaternion>() {return true;}