[37] | 1 | function r = unit(a) |
---|
| 2 | % UNIT quaternion. Divides a quaternion by its own modulus. |
---|
| 3 | % The result is a quaternion with unit modulus. |
---|
| 4 | |
---|
| 5 | % Copyright © 2005 Stephen J. Sangwine and Nicolas Le Bihan. |
---|
| 6 | % See the file : Copyright.m for further details. |
---|
| 7 | |
---|
| 8 | error(nargchk(1, 1, nargin)), error(nargoutchk(0, 1, nargout)) |
---|
| 9 | |
---|
| 10 | m = abs(a); |
---|
| 11 | |
---|
| 12 | % Since m could be complex, the warning check below has to |
---|
| 13 | % use abs again to get a real result for comparison with eps. |
---|
| 14 | |
---|
| 15 | if any(any(abs(m) < eps)) |
---|
| 16 | warning(['At least one element has zero modulus, '... |
---|
| 17 | 'and divide by zero will occur.']); |
---|
| 18 | end |
---|
| 19 | |
---|
| 20 | r = a ./ m; |
---|
| 21 | |
---|
| 22 | % Dividing by a small modulus can result in numerical errors such that the |
---|
| 23 | % result does not have unit modulus. This is especially likely with complex |
---|
| 24 | % quaternions, so we perform a check here. |
---|
| 25 | |
---|
| 26 | % The modulus of each element of r should be either 1 or 1 + 0i. In either |
---|
| 27 | % case, subtracting ones should result in a real or complex number with |
---|
| 28 | % very small modulus. Therefore taking the modulus of the difference should |
---|
| 29 | % yield small values, which we compare with epsilon using an arbitrary |
---|
| 30 | % scale factor which we choose so that the test is not too sensitive. |
---|
| 31 | |
---|
| 32 | m = abs(abs(r) - ones(size(r))); |
---|
| 33 | |
---|
| 34 | if any(any(m > 100.*eps)) |
---|
| 35 | warning(['At least one element of the result of the unit '... |
---|
| 36 | 'function has a modulus which is not accurately 1.']); |
---|
| 37 | end |
---|
| 38 | |
---|
| 39 | % Discussion note. Unit complex pure quaternions must have real and |
---|
| 40 | % imaginary parts with moduli b and d such that b^2 - d^2 = 1. See: |
---|
| 41 | % |
---|
| 42 | % Sangwine, S. J., Biquaternion (complexified quaternion) roots of -1, |
---|
| 43 | % arXiv:math.RA/0506190, 10 June 2005, available at http://www.arxiv.org/. |
---|
| 44 | % |
---|
| 45 | % It is possible for b and d to be large, and then the difference between |
---|
| 46 | % their squares can be lost below the least significant bits. Therefore, |
---|
| 47 | % although it is possible in theory to construct a unit pure complex |
---|
| 48 | % quaternion by applying the unit function to an arbitrary pure complex |
---|
| 49 | % quaternion, it is much better to create the unit pure complex quaternion |
---|
| 50 | % accurately in the first place by choosing correct values for b and d. |
---|
| 51 | % |
---|
| 52 | % SJS 28 September 2005 |
---|