Orthogonal plane decomposition
[p, m] = opd(q, a, b)
opd(q, a ,b) decomposes a quaternion array element-by-element q into two components in orthogonal planes defined by a and b using the formulae p = ½(q + aqb) and m = ½(q - aqb). The third parameter is optional, if omitted it is set equal to the second. If a is a pure quaternion, and b is omitted, p is parallel to a, and m is in the plane normal to a. (Parallel/perpendicular decomposition.) Other cases are determined by the values of the second and third parameters. Coxeter's paper is recommended as a reference for details of the geometry.
>> q = randq; >> [p,m] = opd(q, randv) p = 2.776e-17 - 0.1426 * I - 0.8066 * J - 0.07652 * K m = -0.4503 + 0.3269 * I - 0.0668 * J + 0.09476 * K >> scalar_product(p,m) ans = -3.9899e-17 % p and m are orthogonal.