octave_packages/quaternion-2.0.0/@quaternion/power.m

   1 ## Copyright (C) 2010, 2011   Lukas F. Reichlin
   2 ## Copyright (c) 2011 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
   3 ##
   4 ## This program is free software: you can redistribute it and/or modify
   5 ## it under the terms of the GNU General Public License as published by
   6 ## the Free Software Foundation, either version 3 of the License, or
   7 ## (at your option) any later version.
   8 ##
   9 ## This program is distributed in the hope that it will be useful,
  10 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12 ## GNU General Public License for more details.
  13 ##
  14 ## You should have received a copy of the GNU General Public License
  15 ## along with this program.  If not, see <http://www.gnu.org/licenses/>.
  16
  17 ## -*- texinfo -*-
  18 ## Power operator of quaternions.  Used by Octave for "q.^x".
  19 ## Exponent x can be scalar or of appropriate size.
  20
  21 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
  22 ## Created: May 2010
  23 ## Version: 0.3
  24
  25 function a = power (a, b)
  26
  27   if (isa (b, "quaternion"))          # exponent is a quaternion
  28     a = exp (log (a) .* b);           # a could be real, but log doesn't care
  29   elseif (! isreal (b))
  30     error ("quaternion:invalidArgument", "quaternion: power: invalid exponent");
  31   elseif (b == -1)                    # special case for ldivide and rdivide
  32     norm2 = norm2 (a);                # a is quaternion because b is not,
  33     a.w = a.w ./ norm2;               # otherwise octave wouldn't call
  34     a.x = -a.x ./ norm2;              # quaternion's power operator.
  35     a.y = -a.y ./ norm2;
  36     a.z = -a.z ./ norm2;
  37   else                                # exponent is real
  38     na = abs (a);
  39     nv = normv (a);
  40     th = acos (a.w ./ na);
  41     nab = na.^b;
  42     snt = sin (b.*th);
  43     a.w = nab .* cos (b.*th);
  44     a.x = (a.x ./ nv) .* nab .* snt;
  45     a.y = (a.y ./ nv) .* nab .* snt;
  46     a.z = (a.z ./ nv) .* nab .* snt;
  47   endif
  48
  49 endfunction