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