1## Copyright (C) 2010-2015   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.4
24
25function a = power (a, b)
26
27  if (nargin != 2)
28    error ("quaternion: power: this is a binary operator");
29  endif
30
31  if (isa (b, "quaternion"))          # exponent is a quaternion
32    a = exp (log (a) .* b);           # a could be real, but log doesn't care
33  elseif (! isreal (b))
34    error ("quaternion:invalidArgument", "quaternion: power: invalid exponent");
35  elseif (b == -1)                    # special case for ldivide and rdivide
36    norm2 = norm2 (a);                # a is quaternion because b is not,
37    a.w = a.w ./ norm2;               # otherwise octave wouldn't call
38    a.x = -a.x ./ norm2;              # quaternion's power operator.
39    a.y = -a.y ./ norm2;
40    a.z = -a.z ./ norm2;
41  else                                # exponent is real
42    na = abs (a);
43    nv = normv (a);
44    th = acos (a.w ./ na);
45    nab = na.^b;
46    snt = sin (b.*th);
47    a.w = nab .* cos (b.*th);
48    a.x = (a.x ./ nv) .* nab .* snt;
49    a.y = (a.y ./ nv) .* nab .* snt;
50    a.z = (a.z ./ nv) .* nab .* snt;
51  endif
52
53endfunction
54