Add a useful package (from Source forge) for octave

[CreaPhase.git] / octave_packages / quaternion-2.0.0 / @quaternion / mpower.m

diff --git a/octave_packages/quaternion-2.0.0/@quaternion/mpower.m b/octave_packages/quaternion-2.0.0/@quaternion/mpower.m
new file mode 100644 (file)
index 0000000..e8c865b
--- /dev/null
+++ b/octave_packages/quaternion-2.0.0/@quaternion/mpower.m
@@ -0,0 +1,56 @@
+## Copyright (C) 2010, 2011   Lukas F. Reichlin
+##
+## 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 3 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, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## Matrix power operator of quaternions.  Used by Octave for "q^x".
+
+## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
+## Created: May 2010
+## Version: 0.2
+
+function q = mpower (a, b)
+
+  [r, c] = size (a);
+  
+  if (r != c)
+    error ("quaternion: mpower: quaternion matrix must be square");
+  endif
+  
+  if (r == 1 && c == 1)                 # a scalar, b?
+    q = a .^ b;                         # b could be a quaternion
+  elseif (is_real_array (b) && isscalar (b) && fix (b) == b)
+    e = fix (abs (b));
+    switch (sign (b))
+      case -1                           # q^-e
+        a = inv (a);
+        q = a;
+      case 0                            # q^0
+        q = eye (r);                    # alternative: q = quaternion (eye (r))
+        return;
+      case 1;                           # q^e
+        q = a;
+    endswitch  
+    for k = 2 : e
+      q *= a;                           # improvement?: q^8 = ((q^2)^2)^2, q^9 = (((q^2)^2)^2)*q
+    endfor
+  else
+    error ("quaternion: mpower: case not implemented yet");
+    q = expm (logm (a) * b);            # don't know whether this formula is correct
+  endif
+
+  ## TODO: - q1 ^ q2
+  ##       - arrays
+
+endfunction