Coverage for /usr/lib/python3/dist-packages/sympy/algebras/quaternion.py: 23%

382 statements  

« prev     ^ index     » next       coverage.py v7.9.1, created at 2025-06-14 15:55 +0200

1from sympy.core.numbers import Rational 

2from sympy.core.singleton import S 

3from sympy.core.relational import is_eq 

4from sympy.functions.elementary.complexes import (conjugate, im, re, sign) 

5from sympy.functions.elementary.exponential import (exp, log as ln) 

6from sympy.functions.elementary.miscellaneous import sqrt 

7from sympy.functions.elementary.trigonometric import (acos, asin, atan2) 

8from sympy.functions.elementary.trigonometric import (cos, sin) 

9from sympy.simplify.trigsimp import trigsimp 

10from sympy.integrals.integrals import integrate 

11from sympy.matrices.dense import MutableDenseMatrix as Matrix 

12from sympy.core.sympify import sympify, _sympify 

13from sympy.core.expr import Expr 

14from sympy.core.logic import fuzzy_not, fuzzy_or 

15 

16from mpmath.libmp.libmpf import prec_to_dps 

17 

18 

19def _check_norm(elements, norm): 

20 """validate if input norm is consistent""" 

21 if norm is not None and norm.is_number: 

22 if norm.is_positive is False: 

23 raise ValueError("Input norm must be positive.") 

24 

25 numerical = all(i.is_number and i.is_real is True for i in elements) 

26 if numerical and is_eq(norm**2, sum(i**2 for i in elements)) is False: 

27 raise ValueError("Incompatible value for norm.") 

28 

29 

30def _is_extrinsic(seq): 

31 """validate seq and return True if seq is lowercase and False if uppercase""" 

32 if type(seq) != str: 

33 raise ValueError('Expected seq to be a string.') 

34 if len(seq) != 3: 

35 raise ValueError("Expected 3 axes, got `{}`.".format(seq)) 

36 

37 intrinsic = seq.isupper() 

38 extrinsic = seq.islower() 

39 if not (intrinsic or extrinsic): 

40 raise ValueError("seq must either be fully uppercase (for extrinsic " 

41 "rotations), or fully lowercase, for intrinsic " 

42 "rotations).") 

43 

44 i, j, k = seq.lower() 

45 if (i == j) or (j == k): 

46 raise ValueError("Consecutive axes must be different") 

47 

48 bad = set(seq) - set('xyzXYZ') 

49 if bad: 

50 raise ValueError("Expected axes from `seq` to be from " 

51 "['x', 'y', 'z'] or ['X', 'Y', 'Z'], " 

52 "got {}".format(''.join(bad))) 

53 

54 return extrinsic 

55 

56 

57class Quaternion(Expr): 

58 """Provides basic quaternion operations. 

59 Quaternion objects can be instantiated as Quaternion(a, b, c, d) 

60 as in (a + b*i + c*j + d*k). 

61 

62 Parameters 

63 ========== 

64 

65 norm : None or number 

66 Pre-defined quaternion norm. If a value is given, Quaternion.norm 

67 returns this pre-defined value instead of calculating the norm 

68 

69 Examples 

70 ======== 

71 

72 >>> from sympy import Quaternion 

73 >>> q = Quaternion(1, 2, 3, 4) 

74 >>> q 

75 1 + 2*i + 3*j + 4*k 

76 

77 Quaternions over complex fields can be defined as : 

78 

79 >>> from sympy import Quaternion 

80 >>> from sympy import symbols, I 

81 >>> x = symbols('x') 

82 >>> q1 = Quaternion(x, x**3, x, x**2, real_field = False) 

83 >>> q2 = Quaternion(3 + 4*I, 2 + 5*I, 0, 7 + 8*I, real_field = False) 

84 >>> q1 

85 x + x**3*i + x*j + x**2*k 

86 >>> q2 

87 (3 + 4*I) + (2 + 5*I)*i + 0*j + (7 + 8*I)*k 

88 

89 Defining symbolic unit quaternions: 

90 >>> from sympy import Quaternion 

91 >>> from sympy.abc import w, x, y, z 

92 >>> q = Quaternion(w, x, y, z, norm=1) 

93 >>> q 

94 w + x*i + y*j + z*k 

95 >>> q.norm() 

96 1 

97 

98 References 

99 ========== 

100 

101 .. [1] https://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/ 

102 .. [2] https://en.wikipedia.org/wiki/Quaternion 

103 

104 """ 

105 _op_priority = 11.0 

106 

107 is_commutative = False 

108 

109 def __new__(cls, a=0, b=0, c=0, d=0, real_field=True, norm=None): 

110 a, b, c, d = map(sympify, (a, b, c, d)) 

111 

112 if any(i.is_commutative is False for i in [a, b, c, d]): 

113 raise ValueError("arguments have to be commutative") 

114 else: 

115 obj = Expr.__new__(cls, a, b, c, d) 

116 obj._a = a 

117 obj._b = b 

118 obj._c = c 

119 obj._d = d 

120 obj._real_field = real_field 

121 obj.set_norm(norm) 

122 return obj 

123 

124 def set_norm(self, norm): 

125 """Sets norm of an already instantiated quaternion. 

126 

127 Parameters 

128 ========== 

129 

130 norm : None or number 

131 Pre-defined quaternion norm. If a value is given, Quaternion.norm 

132 returns this pre-defined value instead of calculating the norm 

133 

134 Examples 

135 ======== 

136 

137 >>> from sympy import Quaternion 

138 >>> from sympy.abc import a, b, c, d 

139 >>> q = Quaternion(a, b, c, d) 

140 >>> q.norm() 

141 sqrt(a**2 + b**2 + c**2 + d**2) 

142 

143 Setting the norm: 

144 

145 >>> q.set_norm(1) 

146 >>> q.norm() 

147 1 

148 

149 Removing set norm: 

150 

151 >>> q.set_norm(None) 

152 >>> q.norm() 

153 sqrt(a**2 + b**2 + c**2 + d**2) 

154 

155 """ 

156 norm = sympify(norm) 

157 _check_norm(self.args, norm) 

158 self._norm = norm 

159 

160 @property 

161 def a(self): 

162 return self._a 

163 

164 @property 

165 def b(self): 

166 return self._b 

167 

168 @property 

169 def c(self): 

170 return self._c 

171 

172 @property 

173 def d(self): 

174 return self._d 

175 

176 @property 

177 def real_field(self): 

178 return self._real_field 

179 

180 @property 

181 def product_matrix_left(self): 

182 r"""Returns 4 x 4 Matrix equivalent to a Hamilton product from the 

183 left. This can be useful when treating quaternion elements as column 

184 vectors. Given a quaternion $q = a + bi + cj + dk$ where a, b, c and d 

185 are real numbers, the product matrix from the left is: 

186 

187 .. math:: 

188 

189 M = \begin{bmatrix} a &-b &-c &-d \\ 

190 b & a &-d & c \\ 

191 c & d & a &-b \\ 

192 d &-c & b & a \end{bmatrix} 

193 

194 Examples 

195 ======== 

196 

197 >>> from sympy import Quaternion 

198 >>> from sympy.abc import a, b, c, d 

199 >>> q1 = Quaternion(1, 0, 0, 1) 

200 >>> q2 = Quaternion(a, b, c, d) 

201 >>> q1.product_matrix_left 

202 Matrix([ 

203 [1, 0, 0, -1], 

204 [0, 1, -1, 0], 

205 [0, 1, 1, 0], 

206 [1, 0, 0, 1]]) 

207 

208 >>> q1.product_matrix_left * q2.to_Matrix() 

209 Matrix([ 

210 [a - d], 

211 [b - c], 

212 [b + c], 

213 [a + d]]) 

214 

215 This is equivalent to: 

216 

217 >>> (q1 * q2).to_Matrix() 

218 Matrix([ 

219 [a - d], 

220 [b - c], 

221 [b + c], 

222 [a + d]]) 

223 """ 

224 return Matrix([ 

225 [self.a, -self.b, -self.c, -self.d], 

226 [self.b, self.a, -self.d, self.c], 

227 [self.c, self.d, self.a, -self.b], 

228 [self.d, -self.c, self.b, self.a]]) 

229 

230 @property 

231 def product_matrix_right(self): 

232 r"""Returns 4 x 4 Matrix equivalent to a Hamilton product from the 

233 right. This can be useful when treating quaternion elements as column 

234 vectors. Given a quaternion $q = a + bi + cj + dk$ where a, b, c and d 

235 are real numbers, the product matrix from the left is: 

236 

237 .. math:: 

238 

239 M = \begin{bmatrix} a &-b &-c &-d \\ 

240 b & a & d &-c \\ 

241 c &-d & a & b \\ 

242 d & c &-b & a \end{bmatrix} 

243 

244 

245 Examples 

246 ======== 

247 

248 >>> from sympy import Quaternion 

249 >>> from sympy.abc import a, b, c, d 

250 >>> q1 = Quaternion(a, b, c, d) 

251 >>> q2 = Quaternion(1, 0, 0, 1) 

252 >>> q2.product_matrix_right 

253 Matrix([ 

254 [1, 0, 0, -1], 

255 [0, 1, 1, 0], 

256 [0, -1, 1, 0], 

257 [1, 0, 0, 1]]) 

258 

259 Note the switched arguments: the matrix represents the quaternion on 

260 the right, but is still considered as a matrix multiplication from the 

261 left. 

262 

263 >>> q2.product_matrix_right * q1.to_Matrix() 

264 Matrix([ 

265 [ a - d], 

266 [ b + c], 

267 [-b + c], 

268 [ a + d]]) 

269 

270 This is equivalent to: 

271 

272 >>> (q1 * q2).to_Matrix() 

273 Matrix([ 

274 [ a - d], 

275 [ b + c], 

276 [-b + c], 

277 [ a + d]]) 

278 """ 

279 return Matrix([ 

280 [self.a, -self.b, -self.c, -self.d], 

281 [self.b, self.a, self.d, -self.c], 

282 [self.c, -self.d, self.a, self.b], 

283 [self.d, self.c, -self.b, self.a]]) 

284 

285 def to_Matrix(self, vector_only=False): 

286 """Returns elements of quaternion as a column vector. 

287 By default, a Matrix of length 4 is returned, with the real part as the 

288 first element. 

289 If vector_only is True, returns only imaginary part as a Matrix of 

290 length 3. 

291 

292 Parameters 

293 ========== 

294 

295 vector_only : bool 

296 If True, only imaginary part is returned. 

297 Default value: False 

298 

299 Returns 

300 ======= 

301 

302 Matrix 

303 A column vector constructed by the elements of the quaternion. 

304 

305 Examples 

306 ======== 

307 

308 >>> from sympy import Quaternion 

309 >>> from sympy.abc import a, b, c, d 

310 >>> q = Quaternion(a, b, c, d) 

311 >>> q 

312 a + b*i + c*j + d*k 

313 

314 >>> q.to_Matrix() 

315 Matrix([ 

316 [a], 

317 [b], 

318 [c], 

319 [d]]) 

320 

321 

322 >>> q.to_Matrix(vector_only=True) 

323 Matrix([ 

324 [b], 

325 [c], 

326 [d]]) 

327 

328 """ 

329 if vector_only: 

330 return Matrix(self.args[1:]) 

331 else: 

332 return Matrix(self.args) 

333 

334 @classmethod 

335 def from_Matrix(cls, elements): 

336 """Returns quaternion from elements of a column vector`. 

337 If vector_only is True, returns only imaginary part as a Matrix of 

338 length 3. 

339 

340 Parameters 

341 ========== 

342 

343 elements : Matrix, list or tuple of length 3 or 4. If length is 3, 

344 assume real part is zero. 

345 Default value: False 

346 

347 Returns 

348 ======= 

349 

350 Quaternion 

351 A quaternion created from the input elements. 

352 

353 Examples 

354 ======== 

355 

356 >>> from sympy import Quaternion 

357 >>> from sympy.abc import a, b, c, d 

358 >>> q = Quaternion.from_Matrix([a, b, c, d]) 

359 >>> q 

360 a + b*i + c*j + d*k 

361 

362 >>> q = Quaternion.from_Matrix([b, c, d]) 

363 >>> q 

364 0 + b*i + c*j + d*k 

365 

366 """ 

367 length = len(elements) 

368 if length != 3 and length != 4: 

369 raise ValueError("Input elements must have length 3 or 4, got {} " 

370 "elements".format(length)) 

371 

372 if length == 3: 

373 return Quaternion(0, *elements) 

374 else: 

375 return Quaternion(*elements) 

376 

377 @classmethod 

378 def from_euler(cls, angles, seq): 

379 """Returns quaternion equivalent to rotation represented by the Euler 

380 angles, in the sequence defined by ``seq``. 

381 

382 Parameters 

383 ========== 

384 

385 angles : list, tuple or Matrix of 3 numbers 

386 The Euler angles (in radians). 

387 seq : string of length 3 

388 Represents the sequence of rotations. 

389 For intrinsic rotations, seq must be all lowercase and its elements 

390 must be from the set ``{'x', 'y', 'z'}`` 

391 For extrinsic rotations, seq must be all uppercase and its elements 

392 must be from the set ``{'X', 'Y', 'Z'}`` 

393 

394 Returns 

395 ======= 

396 

397 Quaternion 

398 The normalized rotation quaternion calculated from the Euler angles 

399 in the given sequence. 

400 

401 Examples 

402 ======== 

403 

404 >>> from sympy import Quaternion 

405 >>> from sympy import pi 

406 >>> q = Quaternion.from_euler([pi/2, 0, 0], 'xyz') 

407 >>> q 

408 sqrt(2)/2 + sqrt(2)/2*i + 0*j + 0*k 

409 

410 >>> q = Quaternion.from_euler([0, pi/2, pi] , 'zyz') 

411 >>> q 

412 0 + (-sqrt(2)/2)*i + 0*j + sqrt(2)/2*k 

413 

414 >>> q = Quaternion.from_euler([0, pi/2, pi] , 'ZYZ') 

415 >>> q 

416 0 + sqrt(2)/2*i + 0*j + sqrt(2)/2*k 

417 

418 """ 

419 

420 if len(angles) != 3: 

421 raise ValueError("3 angles must be given.") 

422 

423 extrinsic = _is_extrinsic(seq) 

424 i, j, k = seq.lower() 

425 

426 # get elementary basis vectors 

427 ei = [1 if n == i else 0 for n in 'xyz'] 

428 ej = [1 if n == j else 0 for n in 'xyz'] 

429 ek = [1 if n == k else 0 for n in 'xyz'] 

430 

431 # calculate distinct quaternions 

432 qi = cls.from_axis_angle(ei, angles[0]) 

433 qj = cls.from_axis_angle(ej, angles[1]) 

434 qk = cls.from_axis_angle(ek, angles[2]) 

435 

436 if extrinsic: 

437 return trigsimp(qk * qj * qi) 

438 else: 

439 return trigsimp(qi * qj * qk) 

440 

441 def to_euler(self, seq, angle_addition=True, avoid_square_root=False): 

442 r"""Returns Euler angles representing same rotation as the quaternion, 

443 in the sequence given by ``seq``. This implements the method described 

444 in [1]_. 

445 

446 For degenerate cases (gymbal lock cases), the third angle is 

447 set to zero. 

448 

449 Parameters 

450 ========== 

451 

452 seq : string of length 3 

453 Represents the sequence of rotations. 

454 For intrinsic rotations, seq must be all lowercase and its elements 

455 must be from the set ``{'x', 'y', 'z'}`` 

456 For extrinsic rotations, seq must be all uppercase and its elements 

457 must be from the set ``{'X', 'Y', 'Z'}`` 

458 

459 angle_addition : bool 

460 When True, first and third angles are given as an addition and 

461 subtraction of two simpler ``atan2`` expressions. When False, the 

462 first and third angles are each given by a single more complicated 

463 ``atan2`` expression. This equivalent expression is given by: 

464 

465 .. math:: 

466 

467 \operatorname{atan_2} (b,a) \pm \operatorname{atan_2} (d,c) = 

468 \operatorname{atan_2} (bc\pm ad, ac\mp bd) 

469 

470 Default value: True 

471 

472 avoid_square_root : bool 

473 When True, the second angle is calculated with an expression based 

474 on ``acos``, which is slightly more complicated but avoids a square 

475 root. When False, second angle is calculated with ``atan2``, which 

476 is simpler and can be better for numerical reasons (some 

477 numerical implementations of ``acos`` have problems near zero). 

478 Default value: False 

479 

480 

481 Returns 

482 ======= 

483 

484 Tuple 

485 The Euler angles calculated from the quaternion 

486 

487 Examples 

488 ======== 

489 

490 >>> from sympy import Quaternion 

491 >>> from sympy.abc import a, b, c, d 

492 >>> euler = Quaternion(a, b, c, d).to_euler('zyz') 

493 >>> euler 

494 (-atan2(-b, c) + atan2(d, a), 

495 2*atan2(sqrt(b**2 + c**2), sqrt(a**2 + d**2)), 

496 atan2(-b, c) + atan2(d, a)) 

497 

498 

499 References 

500 ========== 

501 

502 .. [1] https://doi.org/10.1371/journal.pone.0276302 

503 

504 """ 

505 if self.is_zero_quaternion(): 

506 raise ValueError('Cannot convert a quaternion with norm 0.') 

507 

508 angles = [0, 0, 0] 

509 

510 extrinsic = _is_extrinsic(seq) 

511 i, j, k = seq.lower() 

512 

513 # get index corresponding to elementary basis vectors 

514 i = 'xyz'.index(i) + 1 

515 j = 'xyz'.index(j) + 1 

516 k = 'xyz'.index(k) + 1 

517 

518 if not extrinsic: 

519 i, k = k, i 

520 

521 # check if sequence is symmetric 

522 symmetric = i == k 

523 if symmetric: 

524 k = 6 - i - j 

525 

526 # parity of the permutation 

527 sign = (i - j) * (j - k) * (k - i) // 2 

528 

529 # permutate elements 

530 elements = [self.a, self.b, self.c, self.d] 

531 a = elements[0] 

532 b = elements[i] 

533 c = elements[j] 

534 d = elements[k] * sign 

535 

536 if not symmetric: 

537 a, b, c, d = a - c, b + d, c + a, d - b 

538 

539 if avoid_square_root: 

540 if symmetric: 

541 n2 = self.norm()**2 

542 angles[1] = acos((a * a + b * b - c * c - d * d) / n2) 

543 else: 

544 n2 = 2 * self.norm()**2 

545 angles[1] = asin((c * c + d * d - a * a - b * b) / n2) 

546 else: 

547 angles[1] = 2 * atan2(sqrt(c * c + d * d), sqrt(a * a + b * b)) 

548 if not symmetric: 

549 angles[1] -= S.Pi / 2 

550 

551 # Check for singularities in numerical cases 

552 case = 0 

553 if is_eq(c, S.Zero) and is_eq(d, S.Zero): 

554 case = 1 

555 if is_eq(a, S.Zero) and is_eq(b, S.Zero): 

556 case = 2 

557 

558 if case == 0: 

559 if angle_addition: 

560 angles[0] = atan2(b, a) + atan2(d, c) 

561 angles[2] = atan2(b, a) - atan2(d, c) 

562 else: 

563 angles[0] = atan2(b*c + a*d, a*c - b*d) 

564 angles[2] = atan2(b*c - a*d, a*c + b*d) 

565 

566 else: # any degenerate case 

567 angles[2 * (not extrinsic)] = S.Zero 

568 if case == 1: 

569 angles[2 * extrinsic] = 2 * atan2(b, a) 

570 else: 

571 angles[2 * extrinsic] = 2 * atan2(d, c) 

572 angles[2 * extrinsic] *= (-1 if extrinsic else 1) 

573 

574 # for Tait-Bryan angles 

575 if not symmetric: 

576 angles[0] *= sign 

577 

578 if extrinsic: 

579 return tuple(angles[::-1]) 

580 else: 

581 return tuple(angles) 

582 

583 @classmethod 

584 def from_axis_angle(cls, vector, angle): 

585 """Returns a rotation quaternion given the axis and the angle of rotation. 

586 

587 Parameters 

588 ========== 

589 

590 vector : tuple of three numbers 

591 The vector representation of the given axis. 

592 angle : number 

593 The angle by which axis is rotated (in radians). 

594 

595 Returns 

596 ======= 

597 

598 Quaternion 

599 The normalized rotation quaternion calculated from the given axis and the angle of rotation. 

600 

601 Examples 

602 ======== 

603 

604 >>> from sympy import Quaternion 

605 >>> from sympy import pi, sqrt 

606 >>> q = Quaternion.from_axis_angle((sqrt(3)/3, sqrt(3)/3, sqrt(3)/3), 2*pi/3) 

607 >>> q 

608 1/2 + 1/2*i + 1/2*j + 1/2*k 

609 

610 """ 

611 (x, y, z) = vector 

612 norm = sqrt(x**2 + y**2 + z**2) 

613 (x, y, z) = (x / norm, y / norm, z / norm) 

614 s = sin(angle * S.Half) 

615 a = cos(angle * S.Half) 

616 b = x * s 

617 c = y * s 

618 d = z * s 

619 

620 # note that this quaternion is already normalized by construction: 

621 # c^2 + (s*x)^2 + (s*y)^2 + (s*z)^2 = c^2 + s^2*(x^2 + y^2 + z^2) = c^2 + s^2 * 1 = c^2 + s^2 = 1 

622 # so, what we return is a normalized quaternion 

623 

624 return cls(a, b, c, d) 

625 

626 @classmethod 

627 def from_rotation_matrix(cls, M): 

628 """Returns the equivalent quaternion of a matrix. The quaternion will be normalized 

629 only if the matrix is special orthogonal (orthogonal and det(M) = 1). 

630 

631 Parameters 

632 ========== 

633 

634 M : Matrix 

635 Input matrix to be converted to equivalent quaternion. M must be special 

636 orthogonal (orthogonal and det(M) = 1) for the quaternion to be normalized. 

637 

638 Returns 

639 ======= 

640 

641 Quaternion 

642 The quaternion equivalent to given matrix. 

643 

644 Examples 

645 ======== 

646 

647 >>> from sympy import Quaternion 

648 >>> from sympy import Matrix, symbols, cos, sin, trigsimp 

649 >>> x = symbols('x') 

650 >>> M = Matrix([[cos(x), -sin(x), 0], [sin(x), cos(x), 0], [0, 0, 1]]) 

651 >>> q = trigsimp(Quaternion.from_rotation_matrix(M)) 

652 >>> q 

653 sqrt(2)*sqrt(cos(x) + 1)/2 + 0*i + 0*j + sqrt(2 - 2*cos(x))*sign(sin(x))/2*k 

654 

655 """ 

656 

657 absQ = M.det()**Rational(1, 3) 

658 

659 a = sqrt(absQ + M[0, 0] + M[1, 1] + M[2, 2]) / 2 

660 b = sqrt(absQ + M[0, 0] - M[1, 1] - M[2, 2]) / 2 

661 c = sqrt(absQ - M[0, 0] + M[1, 1] - M[2, 2]) / 2 

662 d = sqrt(absQ - M[0, 0] - M[1, 1] + M[2, 2]) / 2 

663 

664 b = b * sign(M[2, 1] - M[1, 2]) 

665 c = c * sign(M[0, 2] - M[2, 0]) 

666 d = d * sign(M[1, 0] - M[0, 1]) 

667 

668 return Quaternion(a, b, c, d) 

669 

670 def __add__(self, other): 

671 return self.add(other) 

672 

673 def __radd__(self, other): 

674 return self.add(other) 

675 

676 def __sub__(self, other): 

677 return self.add(other*-1) 

678 

679 def __mul__(self, other): 

680 return self._generic_mul(self, _sympify(other)) 

681 

682 def __rmul__(self, other): 

683 return self._generic_mul(_sympify(other), self) 

684 

685 def __pow__(self, p): 

686 return self.pow(p) 

687 

688 def __neg__(self): 

689 return Quaternion(-self._a, -self._b, -self._c, -self.d) 

690 

691 def __truediv__(self, other): 

692 return self * sympify(other)**-1 

693 

694 def __rtruediv__(self, other): 

695 return sympify(other) * self**-1 

696 

697 def _eval_Integral(self, *args): 

698 return self.integrate(*args) 

699 

700 def diff(self, *symbols, **kwargs): 

701 kwargs.setdefault('evaluate', True) 

702 return self.func(*[a.diff(*symbols, **kwargs) for a in self.args]) 

703 

704 def add(self, other): 

705 """Adds quaternions. 

706 

707 Parameters 

708 ========== 

709 

710 other : Quaternion 

711 The quaternion to add to current (self) quaternion. 

712 

713 Returns 

714 ======= 

715 

716 Quaternion 

717 The resultant quaternion after adding self to other 

718 

719 Examples 

720 ======== 

721 

722 >>> from sympy import Quaternion 

723 >>> from sympy import symbols 

724 >>> q1 = Quaternion(1, 2, 3, 4) 

725 >>> q2 = Quaternion(5, 6, 7, 8) 

726 >>> q1.add(q2) 

727 6 + 8*i + 10*j + 12*k 

728 >>> q1 + 5 

729 6 + 2*i + 3*j + 4*k 

730 >>> x = symbols('x', real = True) 

731 >>> q1.add(x) 

732 (x + 1) + 2*i + 3*j + 4*k 

733 

734 Quaternions over complex fields : 

735 

736 >>> from sympy import Quaternion 

737 >>> from sympy import I 

738 >>> q3 = Quaternion(3 + 4*I, 2 + 5*I, 0, 7 + 8*I, real_field = False) 

739 >>> q3.add(2 + 3*I) 

740 (5 + 7*I) + (2 + 5*I)*i + 0*j + (7 + 8*I)*k 

741 

742 """ 

743 q1 = self 

744 q2 = sympify(other) 

745 

746 # If q2 is a number or a SymPy expression instead of a quaternion 

747 if not isinstance(q2, Quaternion): 

748 if q1.real_field and q2.is_complex: 

749 return Quaternion(re(q2) + q1.a, im(q2) + q1.b, q1.c, q1.d) 

750 elif q2.is_commutative: 

751 return Quaternion(q1.a + q2, q1.b, q1.c, q1.d) 

752 else: 

753 raise ValueError("Only commutative expressions can be added with a Quaternion.") 

754 

755 return Quaternion(q1.a + q2.a, q1.b + q2.b, q1.c + q2.c, q1.d 

756 + q2.d) 

757 

758 def mul(self, other): 

759 """Multiplies quaternions. 

760 

761 Parameters 

762 ========== 

763 

764 other : Quaternion or symbol 

765 The quaternion to multiply to current (self) quaternion. 

766 

767 Returns 

768 ======= 

769 

770 Quaternion 

771 The resultant quaternion after multiplying self with other 

772 

773 Examples 

774 ======== 

775 

776 >>> from sympy import Quaternion 

777 >>> from sympy import symbols 

778 >>> q1 = Quaternion(1, 2, 3, 4) 

779 >>> q2 = Quaternion(5, 6, 7, 8) 

780 >>> q1.mul(q2) 

781 (-60) + 12*i + 30*j + 24*k 

782 >>> q1.mul(2) 

783 2 + 4*i + 6*j + 8*k 

784 >>> x = symbols('x', real = True) 

785 >>> q1.mul(x) 

786 x + 2*x*i + 3*x*j + 4*x*k 

787 

788 Quaternions over complex fields : 

789 

790 >>> from sympy import Quaternion 

791 >>> from sympy import I 

792 >>> q3 = Quaternion(3 + 4*I, 2 + 5*I, 0, 7 + 8*I, real_field = False) 

793 >>> q3.mul(2 + 3*I) 

794 (2 + 3*I)*(3 + 4*I) + (2 + 3*I)*(2 + 5*I)*i + 0*j + (2 + 3*I)*(7 + 8*I)*k 

795 

796 """ 

797 return self._generic_mul(self, _sympify(other)) 

798 

799 @staticmethod 

800 def _generic_mul(q1, q2): 

801 """Generic multiplication. 

802 

803 Parameters 

804 ========== 

805 

806 q1 : Quaternion or symbol 

807 q2 : Quaternion or symbol 

808 

809 It is important to note that if neither q1 nor q2 is a Quaternion, 

810 this function simply returns q1 * q2. 

811 

812 Returns 

813 ======= 

814 

815 Quaternion 

816 The resultant quaternion after multiplying q1 and q2 

817 

818 Examples 

819 ======== 

820 

821 >>> from sympy import Quaternion 

822 >>> from sympy import Symbol, S 

823 >>> q1 = Quaternion(1, 2, 3, 4) 

824 >>> q2 = Quaternion(5, 6, 7, 8) 

825 >>> Quaternion._generic_mul(q1, q2) 

826 (-60) + 12*i + 30*j + 24*k 

827 >>> Quaternion._generic_mul(q1, S(2)) 

828 2 + 4*i + 6*j + 8*k 

829 >>> x = Symbol('x', real = True) 

830 >>> Quaternion._generic_mul(q1, x) 

831 x + 2*x*i + 3*x*j + 4*x*k 

832 

833 Quaternions over complex fields : 

834 

835 >>> from sympy import I 

836 >>> q3 = Quaternion(3 + 4*I, 2 + 5*I, 0, 7 + 8*I, real_field = False) 

837 >>> Quaternion._generic_mul(q3, 2 + 3*I) 

838 (2 + 3*I)*(3 + 4*I) + (2 + 3*I)*(2 + 5*I)*i + 0*j + (2 + 3*I)*(7 + 8*I)*k 

839 

840 """ 

841 # None is a Quaternion: 

842 if not isinstance(q1, Quaternion) and not isinstance(q2, Quaternion): 

843 return q1 * q2 

844 

845 # If q1 is a number or a SymPy expression instead of a quaternion 

846 if not isinstance(q1, Quaternion): 

847 if q2.real_field and q1.is_complex: 

848 return Quaternion(re(q1), im(q1), 0, 0) * q2 

849 elif q1.is_commutative: 

850 return Quaternion(q1 * q2.a, q1 * q2.b, q1 * q2.c, q1 * q2.d) 

851 else: 

852 raise ValueError("Only commutative expressions can be multiplied with a Quaternion.") 

853 

854 # If q2 is a number or a SymPy expression instead of a quaternion 

855 if not isinstance(q2, Quaternion): 

856 if q1.real_field and q2.is_complex: 

857 return q1 * Quaternion(re(q2), im(q2), 0, 0) 

858 elif q2.is_commutative: 

859 return Quaternion(q2 * q1.a, q2 * q1.b, q2 * q1.c, q2 * q1.d) 

860 else: 

861 raise ValueError("Only commutative expressions can be multiplied with a Quaternion.") 

862 

863 # If any of the quaternions has a fixed norm, pre-compute norm 

864 if q1._norm is None and q2._norm is None: 

865 norm = None 

866 else: 

867 norm = q1.norm() * q2.norm() 

868 

869 return Quaternion(-q1.b*q2.b - q1.c*q2.c - q1.d*q2.d + q1.a*q2.a, 

870 q1.b*q2.a + q1.c*q2.d - q1.d*q2.c + q1.a*q2.b, 

871 -q1.b*q2.d + q1.c*q2.a + q1.d*q2.b + q1.a*q2.c, 

872 q1.b*q2.c - q1.c*q2.b + q1.d*q2.a + q1.a * q2.d, 

873 norm=norm) 

874 

875 def _eval_conjugate(self): 

876 """Returns the conjugate of the quaternion.""" 

877 q = self 

878 return Quaternion(q.a, -q.b, -q.c, -q.d, norm=q._norm) 

879 

880 def norm(self): 

881 """Returns the norm of the quaternion.""" 

882 if self._norm is None: # check if norm is pre-defined 

883 q = self 

884 # trigsimp is used to simplify sin(x)^2 + cos(x)^2 (these terms 

885 # arise when from_axis_angle is used). 

886 self._norm = sqrt(trigsimp(q.a**2 + q.b**2 + q.c**2 + q.d**2)) 

887 

888 return self._norm 

889 

890 def normalize(self): 

891 """Returns the normalized form of the quaternion.""" 

892 q = self 

893 return q * (1/q.norm()) 

894 

895 def inverse(self): 

896 """Returns the inverse of the quaternion.""" 

897 q = self 

898 if not q.norm(): 

899 raise ValueError("Cannot compute inverse for a quaternion with zero norm") 

900 return conjugate(q) * (1/q.norm()**2) 

901 

902 def pow(self, p): 

903 """Finds the pth power of the quaternion. 

904 

905 Parameters 

906 ========== 

907 

908 p : int 

909 Power to be applied on quaternion. 

910 

911 Returns 

912 ======= 

913 

914 Quaternion 

915 Returns the p-th power of the current quaternion. 

916 Returns the inverse if p = -1. 

917 

918 Examples 

919 ======== 

920 

921 >>> from sympy import Quaternion 

922 >>> q = Quaternion(1, 2, 3, 4) 

923 >>> q.pow(4) 

924 668 + (-224)*i + (-336)*j + (-448)*k 

925 

926 """ 

927 p = sympify(p) 

928 q = self 

929 if p == -1: 

930 return q.inverse() 

931 res = 1 

932 

933 if not p.is_Integer: 

934 return NotImplemented 

935 

936 if p < 0: 

937 q, p = q.inverse(), -p 

938 

939 while p > 0: 

940 if p % 2 == 1: 

941 res = q * res 

942 

943 p = p//2 

944 q = q * q 

945 

946 return res 

947 

948 def exp(self): 

949 """Returns the exponential of q (e^q). 

950 

951 Returns 

952 ======= 

953 

954 Quaternion 

955 Exponential of q (e^q). 

956 

957 Examples 

958 ======== 

959 

960 >>> from sympy import Quaternion 

961 >>> q = Quaternion(1, 2, 3, 4) 

962 >>> q.exp() 

963 E*cos(sqrt(29)) 

964 + 2*sqrt(29)*E*sin(sqrt(29))/29*i 

965 + 3*sqrt(29)*E*sin(sqrt(29))/29*j 

966 + 4*sqrt(29)*E*sin(sqrt(29))/29*k 

967 

968 """ 

969 # exp(q) = e^a(cos||v|| + v/||v||*sin||v||) 

970 q = self 

971 vector_norm = sqrt(q.b**2 + q.c**2 + q.d**2) 

972 a = exp(q.a) * cos(vector_norm) 

973 b = exp(q.a) * sin(vector_norm) * q.b / vector_norm 

974 c = exp(q.a) * sin(vector_norm) * q.c / vector_norm 

975 d = exp(q.a) * sin(vector_norm) * q.d / vector_norm 

976 

977 return Quaternion(a, b, c, d) 

978 

979 def _ln(self): 

980 """Returns the natural logarithm of the quaternion (_ln(q)). 

981 

982 Examples 

983 ======== 

984 

985 >>> from sympy import Quaternion 

986 >>> q = Quaternion(1, 2, 3, 4) 

987 >>> q._ln() 

988 log(sqrt(30)) 

989 + 2*sqrt(29)*acos(sqrt(30)/30)/29*i 

990 + 3*sqrt(29)*acos(sqrt(30)/30)/29*j 

991 + 4*sqrt(29)*acos(sqrt(30)/30)/29*k 

992 

993 """ 

994 # _ln(q) = _ln||q|| + v/||v||*arccos(a/||q||) 

995 q = self 

996 vector_norm = sqrt(q.b**2 + q.c**2 + q.d**2) 

997 q_norm = q.norm() 

998 a = ln(q_norm) 

999 b = q.b * acos(q.a / q_norm) / vector_norm 

1000 c = q.c * acos(q.a / q_norm) / vector_norm 

1001 d = q.d * acos(q.a / q_norm) / vector_norm 

1002 

1003 return Quaternion(a, b, c, d) 

1004 

1005 def _eval_subs(self, *args): 

1006 elements = [i.subs(*args) for i in self.args] 

1007 norm = self._norm 

1008 try: 

1009 norm = norm.subs(*args) 

1010 except AttributeError: 

1011 pass 

1012 _check_norm(elements, norm) 

1013 return Quaternion(*elements, norm=norm) 

1014 

1015 def _eval_evalf(self, prec): 

1016 """Returns the floating point approximations (decimal numbers) of the quaternion. 

1017 

1018 Returns 

1019 ======= 

1020 

1021 Quaternion 

1022 Floating point approximations of quaternion(self) 

1023 

1024 Examples 

1025 ======== 

1026 

1027 >>> from sympy import Quaternion 

1028 >>> from sympy import sqrt 

1029 >>> q = Quaternion(1/sqrt(1), 1/sqrt(2), 1/sqrt(3), 1/sqrt(4)) 

1030 >>> q.evalf() 

1031 1.00000000000000 

1032 + 0.707106781186547*i 

1033 + 0.577350269189626*j 

1034 + 0.500000000000000*k 

1035 

1036 """ 

1037 nprec = prec_to_dps(prec) 

1038 return Quaternion(*[arg.evalf(n=nprec) for arg in self.args]) 

1039 

1040 def pow_cos_sin(self, p): 

1041 """Computes the pth power in the cos-sin form. 

1042 

1043 Parameters 

1044 ========== 

1045 

1046 p : int 

1047 Power to be applied on quaternion. 

1048 

1049 Returns 

1050 ======= 

1051 

1052 Quaternion 

1053 The p-th power in the cos-sin form. 

1054 

1055 Examples 

1056 ======== 

1057 

1058 >>> from sympy import Quaternion 

1059 >>> q = Quaternion(1, 2, 3, 4) 

1060 >>> q.pow_cos_sin(4) 

1061 900*cos(4*acos(sqrt(30)/30)) 

1062 + 1800*sqrt(29)*sin(4*acos(sqrt(30)/30))/29*i 

1063 + 2700*sqrt(29)*sin(4*acos(sqrt(30)/30))/29*j 

1064 + 3600*sqrt(29)*sin(4*acos(sqrt(30)/30))/29*k 

1065 

1066 """ 

1067 # q = ||q||*(cos(a) + u*sin(a)) 

1068 # q^p = ||q||^p * (cos(p*a) + u*sin(p*a)) 

1069 

1070 q = self 

1071 (v, angle) = q.to_axis_angle() 

1072 q2 = Quaternion.from_axis_angle(v, p * angle) 

1073 return q2 * (q.norm()**p) 

1074 

1075 def integrate(self, *args): 

1076 """Computes integration of quaternion. 

1077 

1078 Returns 

1079 ======= 

1080 

1081 Quaternion 

1082 Integration of the quaternion(self) with the given variable. 

1083 

1084 Examples 

1085 ======== 

1086 

1087 Indefinite Integral of quaternion : 

1088 

1089 >>> from sympy import Quaternion 

1090 >>> from sympy.abc import x 

1091 >>> q = Quaternion(1, 2, 3, 4) 

1092 >>> q.integrate(x) 

1093 x + 2*x*i + 3*x*j + 4*x*k 

1094 

1095 Definite integral of quaternion : 

1096 

1097 >>> from sympy import Quaternion 

1098 >>> from sympy.abc import x 

1099 >>> q = Quaternion(1, 2, 3, 4) 

1100 >>> q.integrate((x, 1, 5)) 

1101 4 + 8*i + 12*j + 16*k 

1102 

1103 """ 

1104 # TODO: is this expression correct? 

1105 return Quaternion(integrate(self.a, *args), integrate(self.b, *args), 

1106 integrate(self.c, *args), integrate(self.d, *args)) 

1107 

1108 @staticmethod 

1109 def rotate_point(pin, r): 

1110 """Returns the coordinates of the point pin(a 3 tuple) after rotation. 

1111 

1112 Parameters 

1113 ========== 

1114 

1115 pin : tuple 

1116 A 3-element tuple of coordinates of a point which needs to be 

1117 rotated. 

1118 r : Quaternion or tuple 

1119 Axis and angle of rotation. 

1120 

1121 It's important to note that when r is a tuple, it must be of the form 

1122 (axis, angle) 

1123 

1124 Returns 

1125 ======= 

1126 

1127 tuple 

1128 The coordinates of the point after rotation. 

1129 

1130 Examples 

1131 ======== 

1132 

1133 >>> from sympy import Quaternion 

1134 >>> from sympy import symbols, trigsimp, cos, sin 

1135 >>> x = symbols('x') 

1136 >>> q = Quaternion(cos(x/2), 0, 0, sin(x/2)) 

1137 >>> trigsimp(Quaternion.rotate_point((1, 1, 1), q)) 

1138 (sqrt(2)*cos(x + pi/4), sqrt(2)*sin(x + pi/4), 1) 

1139 >>> (axis, angle) = q.to_axis_angle() 

1140 >>> trigsimp(Quaternion.rotate_point((1, 1, 1), (axis, angle))) 

1141 (sqrt(2)*cos(x + pi/4), sqrt(2)*sin(x + pi/4), 1) 

1142 

1143 """ 

1144 if isinstance(r, tuple): 

1145 # if r is of the form (vector, angle) 

1146 q = Quaternion.from_axis_angle(r[0], r[1]) 

1147 else: 

1148 # if r is a quaternion 

1149 q = r.normalize() 

1150 pout = q * Quaternion(0, pin[0], pin[1], pin[2]) * conjugate(q) 

1151 return (pout.b, pout.c, pout.d) 

1152 

1153 def to_axis_angle(self): 

1154 """Returns the axis and angle of rotation of a quaternion. 

1155 

1156 Returns 

1157 ======= 

1158 

1159 tuple 

1160 Tuple of (axis, angle) 

1161 

1162 Examples 

1163 ======== 

1164 

1165 >>> from sympy import Quaternion 

1166 >>> q = Quaternion(1, 1, 1, 1) 

1167 >>> (axis, angle) = q.to_axis_angle() 

1168 >>> axis 

1169 (sqrt(3)/3, sqrt(3)/3, sqrt(3)/3) 

1170 >>> angle 

1171 2*pi/3 

1172 

1173 """ 

1174 q = self 

1175 if q.a.is_negative: 

1176 q = q * -1 

1177 

1178 q = q.normalize() 

1179 angle = trigsimp(2 * acos(q.a)) 

1180 

1181 # Since quaternion is normalised, q.a is less than 1. 

1182 s = sqrt(1 - q.a*q.a) 

1183 

1184 x = trigsimp(q.b / s) 

1185 y = trigsimp(q.c / s) 

1186 z = trigsimp(q.d / s) 

1187 

1188 v = (x, y, z) 

1189 t = (v, angle) 

1190 

1191 return t 

1192 

1193 def to_rotation_matrix(self, v=None, homogeneous=True): 

1194 """Returns the equivalent rotation transformation matrix of the quaternion 

1195 which represents rotation about the origin if v is not passed. 

1196 

1197 Parameters 

1198 ========== 

1199 

1200 v : tuple or None 

1201 Default value: None 

1202 homogeneous : bool 

1203 When True, gives an expression that may be more efficient for 

1204 symbolic calculations but less so for direct evaluation. Both 

1205 formulas are mathematically equivalent. 

1206 Default value: True 

1207 

1208 Returns 

1209 ======= 

1210 

1211 tuple 

1212 Returns the equivalent rotation transformation matrix of the quaternion 

1213 which represents rotation about the origin if v is not passed. 

1214 

1215 Examples 

1216 ======== 

1217 

1218 >>> from sympy import Quaternion 

1219 >>> from sympy import symbols, trigsimp, cos, sin 

1220 >>> x = symbols('x') 

1221 >>> q = Quaternion(cos(x/2), 0, 0, sin(x/2)) 

1222 >>> trigsimp(q.to_rotation_matrix()) 

1223 Matrix([ 

1224 [cos(x), -sin(x), 0], 

1225 [sin(x), cos(x), 0], 

1226 [ 0, 0, 1]]) 

1227 

1228 Generates a 4x4 transformation matrix (used for rotation about a point 

1229 other than the origin) if the point(v) is passed as an argument. 

1230 """ 

1231 

1232 q = self 

1233 s = q.norm()**-2 

1234 

1235 # diagonal elements are different according to parameter normal 

1236 if homogeneous: 

1237 m00 = s*(q.a**2 + q.b**2 - q.c**2 - q.d**2) 

1238 m11 = s*(q.a**2 - q.b**2 + q.c**2 - q.d**2) 

1239 m22 = s*(q.a**2 - q.b**2 - q.c**2 + q.d**2) 

1240 else: 

1241 m00 = 1 - 2*s*(q.c**2 + q.d**2) 

1242 m11 = 1 - 2*s*(q.b**2 + q.d**2) 

1243 m22 = 1 - 2*s*(q.b**2 + q.c**2) 

1244 

1245 m01 = 2*s*(q.b*q.c - q.d*q.a) 

1246 m02 = 2*s*(q.b*q.d + q.c*q.a) 

1247 

1248 m10 = 2*s*(q.b*q.c + q.d*q.a) 

1249 m12 = 2*s*(q.c*q.d - q.b*q.a) 

1250 

1251 m20 = 2*s*(q.b*q.d - q.c*q.a) 

1252 m21 = 2*s*(q.c*q.d + q.b*q.a) 

1253 

1254 if not v: 

1255 return Matrix([[m00, m01, m02], [m10, m11, m12], [m20, m21, m22]]) 

1256 

1257 else: 

1258 (x, y, z) = v 

1259 

1260 m03 = x - x*m00 - y*m01 - z*m02 

1261 m13 = y - x*m10 - y*m11 - z*m12 

1262 m23 = z - x*m20 - y*m21 - z*m22 

1263 m30 = m31 = m32 = 0 

1264 m33 = 1 

1265 

1266 return Matrix([[m00, m01, m02, m03], [m10, m11, m12, m13], 

1267 [m20, m21, m22, m23], [m30, m31, m32, m33]]) 

1268 

1269 def scalar_part(self): 

1270 r"""Returns scalar part($\mathbf{S}(q)$) of the quaternion q. 

1271 

1272 Explanation 

1273 =========== 

1274 

1275 Given a quaternion $q = a + bi + cj + dk$, returns $\mathbf{S}(q) = a$. 

1276 

1277 Examples 

1278 ======== 

1279 

1280 >>> from sympy.algebras.quaternion import Quaternion 

1281 >>> q = Quaternion(4, 8, 13, 12) 

1282 >>> q.scalar_part() 

1283 4 

1284 

1285 """ 

1286 

1287 return self.a 

1288 

1289 def vector_part(self): 

1290 r""" 

1291 Returns vector part($\mathbf{V}(q)$) of the quaternion q. 

1292 

1293 Explanation 

1294 =========== 

1295 

1296 Given a quaternion $q = a + bi + cj + dk$, returns $\mathbf{V}(q) = bi + cj + dk$. 

1297 

1298 Examples 

1299 ======== 

1300 

1301 >>> from sympy.algebras.quaternion import Quaternion 

1302 >>> q = Quaternion(1, 1, 1, 1) 

1303 >>> q.vector_part() 

1304 0 + 1*i + 1*j + 1*k 

1305 

1306 >>> q = Quaternion(4, 8, 13, 12) 

1307 >>> q.vector_part() 

1308 0 + 8*i + 13*j + 12*k 

1309 

1310 """ 

1311 

1312 return Quaternion(0, self.b, self.c, self.d) 

1313 

1314 def axis(self): 

1315 r""" 

1316 Returns the axis($\mathbf{Ax}(q)$) of the quaternion. 

1317 

1318 Explanation 

1319 =========== 

1320 

1321 Given a quaternion $q = a + bi + cj + dk$, returns $\mathbf{Ax}(q)$ i.e., the versor of the vector part of that quaternion 

1322 equal to $\mathbf{U}[\mathbf{V}(q)]$. 

1323 The axis is always an imaginary unit with square equal to $-1 + 0i + 0j + 0k$. 

1324 

1325 Examples 

1326 ======== 

1327 

1328 >>> from sympy.algebras.quaternion import Quaternion 

1329 >>> q = Quaternion(1, 1, 1, 1) 

1330 >>> q.axis() 

1331 0 + sqrt(3)/3*i + sqrt(3)/3*j + sqrt(3)/3*k 

1332 

1333 See Also 

1334 ======== 

1335 

1336 vector_part 

1337 

1338 """ 

1339 axis = self.vector_part().normalize() 

1340 

1341 return Quaternion(0, axis.b, axis.c, axis.d) 

1342 

1343 def is_pure(self): 

1344 """ 

1345 Returns true if the quaternion is pure, false if the quaternion is not pure 

1346 or returns none if it is unknown. 

1347 

1348 Explanation 

1349 =========== 

1350 

1351 A pure quaternion (also a vector quaternion) is a quaternion with scalar 

1352 part equal to 0. 

1353 

1354 Examples 

1355 ======== 

1356 

1357 >>> from sympy.algebras.quaternion import Quaternion 

1358 >>> q = Quaternion(0, 8, 13, 12) 

1359 >>> q.is_pure() 

1360 True 

1361 

1362 See Also 

1363 ======== 

1364 scalar_part 

1365 

1366 """ 

1367 

1368 return self.a.is_zero 

1369 

1370 def is_zero_quaternion(self): 

1371 """ 

1372 Returns true if the quaternion is a zero quaternion or false if it is not a zero quaternion 

1373 and None if the value is unknown. 

1374 

1375 Explanation 

1376 =========== 

1377 

1378 A zero quaternion is a quaternion with both scalar part and 

1379 vector part equal to 0. 

1380 

1381 Examples 

1382 ======== 

1383 

1384 >>> from sympy.algebras.quaternion import Quaternion 

1385 >>> q = Quaternion(1, 0, 0, 0) 

1386 >>> q.is_zero_quaternion() 

1387 False 

1388 

1389 >>> q = Quaternion(0, 0, 0, 0) 

1390 >>> q.is_zero_quaternion() 

1391 True 

1392 

1393 See Also 

1394 ======== 

1395 scalar_part 

1396 vector_part 

1397 

1398 """ 

1399 

1400 return self.norm().is_zero 

1401 

1402 def angle(self): 

1403 r""" 

1404 Returns the angle of the quaternion measured in the real-axis plane. 

1405 

1406 Explanation 

1407 =========== 

1408 

1409 Given a quaternion $q = a + bi + cj + dk$ where a, b, c and d 

1410 are real numbers, returns the angle of the quaternion given by 

1411 

1412 .. math:: 

1413 angle := atan2(\sqrt{b^2 + c^2 + d^2}, {a}) 

1414 

1415 Examples 

1416 ======== 

1417 

1418 >>> from sympy.algebras.quaternion import Quaternion 

1419 >>> q = Quaternion(1, 4, 4, 4) 

1420 >>> q.angle() 

1421 atan(4*sqrt(3)) 

1422 

1423 """ 

1424 

1425 return atan2(self.vector_part().norm(), self.scalar_part()) 

1426 

1427 

1428 def arc_coplanar(self, other): 

1429 """ 

1430 Returns True if the transformation arcs represented by the input quaternions happen in the same plane. 

1431 

1432 Explanation 

1433 =========== 

1434 

1435 Two quaternions are said to be coplanar (in this arc sense) when their axes are parallel. 

1436 The plane of a quaternion is the one normal to its axis. 

1437 

1438 Parameters 

1439 ========== 

1440 

1441 other : a Quaternion 

1442 

1443 Returns 

1444 ======= 

1445 

1446 True : if the planes of the two quaternions are the same, apart from its orientation/sign. 

1447 False : if the planes of the two quaternions are not the same, apart from its orientation/sign. 

1448 None : if plane of either of the quaternion is unknown. 

1449 

1450 Examples 

1451 ======== 

1452 

1453 >>> from sympy.algebras.quaternion import Quaternion 

1454 >>> q1 = Quaternion(1, 4, 4, 4) 

1455 >>> q2 = Quaternion(3, 8, 8, 8) 

1456 >>> Quaternion.arc_coplanar(q1, q2) 

1457 True 

1458 

1459 >>> q1 = Quaternion(2, 8, 13, 12) 

1460 >>> Quaternion.arc_coplanar(q1, q2) 

1461 False 

1462 

1463 See Also 

1464 ======== 

1465 

1466 vector_coplanar 

1467 is_pure 

1468 

1469 """ 

1470 if (self.is_zero_quaternion()) or (other.is_zero_quaternion()): 

1471 raise ValueError('Neither of the given quaternions can be 0') 

1472 

1473 return fuzzy_or([(self.axis() - other.axis()).is_zero_quaternion(), (self.axis() + other.axis()).is_zero_quaternion()]) 

1474 

1475 @classmethod 

1476 def vector_coplanar(cls, q1, q2, q3): 

1477 r""" 

1478 Returns True if the axis of the pure quaternions seen as 3D vectors 

1479 q1, q2, and q3 are coplanar. 

1480 

1481 Explanation 

1482 =========== 

1483 

1484 Three pure quaternions are vector coplanar if the quaternions seen as 3D vectors are coplanar. 

1485 

1486 Parameters 

1487 ========== 

1488 

1489 q1 

1490 A pure Quaternion. 

1491 q2 

1492 A pure Quaternion. 

1493 q3 

1494 A pure Quaternion. 

1495 

1496 Returns 

1497 ======= 

1498 

1499 True : if the axis of the pure quaternions seen as 3D vectors 

1500 q1, q2, and q3 are coplanar. 

1501 False : if the axis of the pure quaternions seen as 3D vectors 

1502 q1, q2, and q3 are not coplanar. 

1503 None : if the axis of the pure quaternions seen as 3D vectors 

1504 q1, q2, and q3 are coplanar is unknown. 

1505 

1506 Examples 

1507 ======== 

1508 

1509 >>> from sympy.algebras.quaternion import Quaternion 

1510 >>> q1 = Quaternion(0, 4, 4, 4) 

1511 >>> q2 = Quaternion(0, 8, 8, 8) 

1512 >>> q3 = Quaternion(0, 24, 24, 24) 

1513 >>> Quaternion.vector_coplanar(q1, q2, q3) 

1514 True 

1515 

1516 >>> q1 = Quaternion(0, 8, 16, 8) 

1517 >>> q2 = Quaternion(0, 8, 3, 12) 

1518 >>> Quaternion.vector_coplanar(q1, q2, q3) 

1519 False 

1520 

1521 See Also 

1522 ======== 

1523 

1524 axis 

1525 is_pure 

1526 

1527 """ 

1528 

1529 if fuzzy_not(q1.is_pure()) or fuzzy_not(q2.is_pure()) or fuzzy_not(q3.is_pure()): 

1530 raise ValueError('The given quaternions must be pure') 

1531 

1532 M = Matrix([[q1.b, q1.c, q1.d], [q2.b, q2.c, q2.d], [q3.b, q3.c, q3.d]]).det() 

1533 return M.is_zero 

1534 

1535 def parallel(self, other): 

1536 """ 

1537 Returns True if the two pure quaternions seen as 3D vectors are parallel. 

1538 

1539 Explanation 

1540 =========== 

1541 

1542 Two pure quaternions are called parallel when their vector product is commutative which 

1543 implies that the quaternions seen as 3D vectors have same direction. 

1544 

1545 Parameters 

1546 ========== 

1547 

1548 other : a Quaternion 

1549 

1550 Returns 

1551 ======= 

1552 

1553 True : if the two pure quaternions seen as 3D vectors are parallel. 

1554 False : if the two pure quaternions seen as 3D vectors are not parallel. 

1555 None : if the two pure quaternions seen as 3D vectors are parallel is unknown. 

1556 

1557 Examples 

1558 ======== 

1559 

1560 >>> from sympy.algebras.quaternion import Quaternion 

1561 >>> q = Quaternion(0, 4, 4, 4) 

1562 >>> q1 = Quaternion(0, 8, 8, 8) 

1563 >>> q.parallel(q1) 

1564 True 

1565 

1566 >>> q1 = Quaternion(0, 8, 13, 12) 

1567 >>> q.parallel(q1) 

1568 False 

1569 

1570 """ 

1571 

1572 if fuzzy_not(self.is_pure()) or fuzzy_not(other.is_pure()): 

1573 raise ValueError('The provided quaternions must be pure') 

1574 

1575 return (self*other - other*self).is_zero_quaternion() 

1576 

1577 def orthogonal(self, other): 

1578 """ 

1579 Returns the orthogonality of two quaternions. 

1580 

1581 Explanation 

1582 =========== 

1583 

1584 Two pure quaternions are called orthogonal when their product is anti-commutative. 

1585 

1586 Parameters 

1587 ========== 

1588 

1589 other : a Quaternion 

1590 

1591 Returns 

1592 ======= 

1593 

1594 True : if the two pure quaternions seen as 3D vectors are orthogonal. 

1595 False : if the two pure quaternions seen as 3D vectors are not orthogonal. 

1596 None : if the two pure quaternions seen as 3D vectors are orthogonal is unknown. 

1597 

1598 Examples 

1599 ======== 

1600 

1601 >>> from sympy.algebras.quaternion import Quaternion 

1602 >>> q = Quaternion(0, 4, 4, 4) 

1603 >>> q1 = Quaternion(0, 8, 8, 8) 

1604 >>> q.orthogonal(q1) 

1605 False 

1606 

1607 >>> q1 = Quaternion(0, 2, 2, 0) 

1608 >>> q = Quaternion(0, 2, -2, 0) 

1609 >>> q.orthogonal(q1) 

1610 True 

1611 

1612 """ 

1613 

1614 if fuzzy_not(self.is_pure()) or fuzzy_not(other.is_pure()): 

1615 raise ValueError('The given quaternions must be pure') 

1616 

1617 return (self*other + other*self).is_zero_quaternion() 

1618 

1619 def index_vector(self): 

1620 r""" 

1621 Returns the index vector of the quaternion. 

1622 

1623 Explanation 

1624 =========== 

1625 

1626 Index vector is given by $\mathbf{T}(q)$ multiplied by $\mathbf{Ax}(q)$ where $\mathbf{Ax}(q)$ is the axis of the quaternion q, 

1627 and mod(q) is the $\mathbf{T}(q)$ (magnitude) of the quaternion. 

1628 

1629 Returns 

1630 ======= 

1631 

1632 Quaternion: representing index vector of the provided quaternion. 

1633 

1634 Examples 

1635 ======== 

1636 

1637 >>> from sympy.algebras.quaternion import Quaternion 

1638 >>> q = Quaternion(2, 4, 2, 4) 

1639 >>> q.index_vector() 

1640 0 + 4*sqrt(10)/3*i + 2*sqrt(10)/3*j + 4*sqrt(10)/3*k 

1641 

1642 See Also 

1643 ======== 

1644 

1645 axis 

1646 norm 

1647 

1648 """ 

1649 

1650 return self.norm() * self.axis() 

1651 

1652 def mensor(self): 

1653 """ 

1654 Returns the natural logarithm of the norm(magnitude) of the quaternion. 

1655 

1656 Examples 

1657 ======== 

1658 

1659 >>> from sympy.algebras.quaternion import Quaternion 

1660 >>> q = Quaternion(2, 4, 2, 4) 

1661 >>> q.mensor() 

1662 log(2*sqrt(10)) 

1663 >>> q.norm() 

1664 2*sqrt(10) 

1665 

1666 See Also 

1667 ======== 

1668 

1669 norm 

1670 

1671 """ 

1672 

1673 return ln(self.norm())