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 |
---|