[37] | 1 | function Z = power(X, Y) |
---|
| 2 | % .^ Array power. |
---|
| 3 | % (Quaternion overloading of standard Matlab function.) |
---|
| 4 | |
---|
| 5 | % Copyright © 2005, 2006 Stephen J. Sangwine and Nicolas Le Bihan. |
---|
| 6 | % See the file : Copyright.m for further details. |
---|
| 7 | |
---|
| 8 | error(nargchk(2, 2, nargin)), error(nargoutchk(0, 1, nargout)) |
---|
| 9 | |
---|
| 10 | % This function is implemented using two methods. For a small set of |
---|
| 11 | % integer powers, an equivalent operation is used. E.g. for Y == -1, the |
---|
| 12 | % elementwise inverse is used, for Y == 2, elementwise squaring is used. |
---|
| 13 | |
---|
| 14 | % For a power of ± 1/2, the sqrt function is used, with or without a |
---|
| 15 | % reciprocal. |
---|
| 16 | |
---|
| 17 | % For powers which are not handled in this way a general formula using |
---|
| 18 | % logarithms is used. |
---|
| 19 | |
---|
| 20 | |
---|
| 21 | if size(Y) ~= [1, 1] |
---|
| 22 | error('Quaternion .^ is not implemented for matrix exponents.'); |
---|
| 23 | else |
---|
| 24 | |
---|
| 25 | % Y is a scalar. Check for and handle the various powers that are dealt |
---|
| 26 | % with as special cases. |
---|
| 27 | |
---|
| 28 | if Y == -2 |
---|
| 29 | Z = (X .* X) .^ -1; % Use the next case recursively. |
---|
| 30 | elseif Y == -1 |
---|
| 31 | Z = conj(X) ./ modsquared(X); % I.e. elementwise inverse. If X has |
---|
| 32 | % zero norm this will result in a |
---|
| 33 | % NaN. |
---|
| 34 | elseif Y == 0 |
---|
| 35 | Z = ones(size(X)); |
---|
| 36 | elseif Y == 1 |
---|
| 37 | Z = X; |
---|
| 38 | elseif Y == 2 |
---|
| 39 | Z = X .* X; |
---|
| 40 | elseif Y == 1/2 |
---|
| 41 | Z = sqrt(X); |
---|
| 42 | elseif Y == -1/2 |
---|
| 43 | Z = sqrt(X .^ -1); |
---|
| 44 | else |
---|
| 45 | |
---|
| 46 | % The general case. The formula used here is taken from |
---|
| 47 | % A quaternion algebra tool set, Doug Sweetser, |
---|
| 48 | % http://www.theworld.com/~sweetser/quaternions/intro/tools/tools.html |
---|
| 49 | |
---|
| 50 | Z = exp(log(X) .* Y); % NB log(X) is the natural logarithm of X. |
---|
| 51 | % (Matlab convention.) |
---|
| 52 | end |
---|
| 53 | end |
---|