source: Revenant/marslib/include/mars_linalg.h@ 7ef8ec4

port/mars-tycoon
Last change on this file since 7ef8ec4 was 80a6a52, checked in by Jonathan Neufeld <support@…>, 3 years ago

Get to a compile state for terrain procedural generation
  • Property mode set to 100755
File size: 40.5 KB
Line 
1#ifndef __MARSLINEARALGEBRA_H__
2#define __MARSLINEARALGEBRA_H__
3
4#include <exception>
5#include <math.h>
6
7#ifdef _DEBUG
8#include <assert.h>
9#endif
10
11#include "mars_meta.h"
12#include "mars_calc.h"
13
14namespace mars
15{
16 class LinAlgEx : public std::exception
17 {
18 };
19
20 class MathOverrides
21 {
22 public:
23 static inline int sqrt (const int x)
24 { return static_cast <int> (::sqrt (static_cast <double> (x))); }
25 static inline short sqrt (const short x)
26 { return static_cast <short> (::sqrtf (static_cast <float> (x))); }
27 static inline long sqrt (const long x)
28 { return static_cast <long> (::sqrt (static_cast <long> (x))); }
29 static inline long double sqrt (const long double x)
30 { return static_cast <long double> (::sqrt (static_cast <long double> (x))); }
31 static inline double sqrt (const double x)
32 { return static_cast <double> (::sqrt (static_cast <double> (x))); }
33 static inline float sqrt (const float x)
34 { return static_cast <float> (::sqrtf (static_cast <float> (x))); }
35
36/*
37 static inline int atan2 (const int x, const int y)
38 { return static_cast <int> (::atan2 (static_cast <double> (x), static_cast <double> (y))); }
39 static inline short atan2 (const short x, const short y)
40 { return static_cast <short> (::atan2f (static_cast <float> (x), static_cast <float> (y))); }
41 static inline long double atan2 (const long double x, const long double y)
42 { return static_cast <long double> (::atan2 (static_cast <long double> (x), static_cast <long double> (y))); }
43 static inline double atan2 (const double x, const double y)
44 { return static_cast <double> (::atan2 (static_cast <double> (x), static_cast <double> (y))); }
45 static inline float atan2 (const float x, const float y)
46 { return static_cast <float> (::atan2f (static_cast <float> (x), static_cast <float> (y))); }
47
48 static inline int acos (const int x)
49 { return static_cast <int> (::acos (static_cast <double> (x))); }
50 static inline short acos (const short x)
51 { return static_cast <short> (::acosf (static_cast <float> (x))); }
52 static inline long double acos (const long double x)
53 { return static_cast <long double> (::acos (static_cast <long double> (x))); }
54 static inline double acos (const double x)
55 { return static_cast <double> (::acos (static_cast <double> (x))); }
56 static inline float acos (const float x)
57 { return static_cast <float> (::acosf (static_cast <float> (x))); }*/
58 };
59
60 // Linear optimization technique, taken from: http://www.flipcode.com/archives/Faster_Vector_Math_Using_Templates.shtml
61 // Operator controllers:
62
63 template <typename T>
64 struct VOpAdd
65 {
66 typedef T ResultType;
67
68 template <int I, typename LVal, typename RVal>
69 static inline T evalOp (const LVal & lval, const RVal & rval)
70 { return lval.template eval<I>() + rval.template eval<I>(); }
71 };
72 template <typename T>
73 struct VOpSub
74 {
75 typedef T ResultType;
76
77 template <int I, typename LVal, typename RVal>
78 static inline T evalOp (const LVal & lval, const RVal & rval)
79 { return lval.template eval<I>() - rval.template eval<I>(); }
80 };
81 template <typename T, int D>
82 struct VOpScale
83 {
84 typedef T ResultType;
85
86 template <int I, typename LVal, typename RVal>
87 static inline T evalOp (const LVal & lval, const RVal & rval)
88 { return lval.template eval<I>() * rval.template eval<I>(); }
89 };
90 template <typename T, int D>
91 struct VOpContract
92 {
93 typedef T ResultType;
94
95 template <int I, typename LVal, typename RVal>
96 static inline T evalOp (const LVal & lval, const RVal & rval)
97 { return lval.template eval<I>() / rval.template eval<I>(); }
98 };
99 template <typename T>
100 struct VOpNegate
101 {
102 typedef T ResultType;
103
104 template <int I, typename RVal>
105 static inline T evalOp (const RVal & rval)
106 { return -rval.template eval<I>(); }
107 };
108 template <typename T, typename TNew>
109 struct VOpCast
110 {
111 typedef TNew ResultType;
112
113 template <int I, typename RVal>
114 static inline TNew evalOp (const RVal & rval)
115 { return static_cast< TNew > (rval.template eval<I>()); }
116 };
117
118 template <typename T, int D>
119 struct VExprTagDim
120 {
121 typedef T Precision;
122 typedef VOpAdd< T > OpAdd;
123 typedef VOpSub< T > OpSub;
124 typedef VOpScale< T, D > OpScale;
125 typedef VOpContract< T, D > OpContract;
126 typedef VOpNegate< T > OpNegate;
127
128 enum { DIMENSION = D };
129 };
130
131 template <typename T, typename Derived, int D>
132 class VBaseExpr : public VExprTagDim< T, D >
133 {
134 public:
135 typedef Derived ExprType;
136
137 private:
138 template <typename Functor>
139 struct Reduce
140 {
141 template <unsigned I>
142 static inline const T go (const Functor & f, const Derived & expr) { return go(f, expr, quantity<I>()); }
143
144 template <unsigned I>
145 static inline const T go (const Derived & expr) { return go(expr, quantity<I>()); }
146
147 private:
148 template <unsigned I>
149 static inline const T go (const Functor & f, const Derived & expr, quantity<I>)
150 {
151 return f.run<I-1>(expr.template eval<I-1>()) + go<I-1>(f, expr);
152 }
153 static inline const T go (const Functor &, const Derived &, quantity<0>)
154 { return static_cast <T> (0); }
155
156 template <unsigned I>
157 static inline const T go (const Derived & expr, quantity<I>)
158 {
159 return Functor::template run<I-1>(expr.template eval<I-1>()) + go<I-1>(expr);
160 }
161 static inline const T go (const Derived &, quantity<0>)
162 { return static_cast <T> (0); }
163 };
164
165 template <typename Expr>
166 struct Equivalence
167 {
168 template <unsigned I>
169 static inline bool test (const Derived & lval, const Expr & rval) { return test(lval, rval, quantity<I>()); }
170
171 private:
172 template <unsigned I>
173 static inline bool test (const Derived & lval, const Expr & rval, quantity<I>)
174 {
175 return (lval.template eval<I-1>() == rval.template eval<I-1>()) && test<I-1>(lval, rval);
176 }
177 static inline bool test (const Derived &, const Expr &, quantity<0>) { return true; }
178 };
179
180 public:
181 template <typename Functor> inline const T reduce () const
182 { return Reduce<Functor>::template go<D>(*static_cast <const Derived *> (this)); }
183
184 template <typename Functor> inline const T reduce (const Functor & f) const
185 { return Reduce<Functor>::template go<D>(f, *static_cast <const Derived *> (this)); }
186
187 template <typename Expr>
188 inline bool equals (const Expr & expr, const typename Expr::ExprType * = 0) const
189 { return Equivalence< Expr >::template test<D>(*static_cast <const Derived *> (this), expr); }
190
191 template <typename V>
192 inline bool equals (const V & v, const typename V::VectorType * = 0) const
193 { return v == *static_cast <const Derived *> (this); }
194 };
195
196 struct VRedSqFn
197 {
198 template <int I, typename T>
199 static inline T run (const T & val) { return SQ(val); }
200 };
201
202 template <typename LVal>
203 class VRedDotEFn
204 {
205 private:
206 const typename LVal::ExprType _lval;
207
208 public:
209 inline VRedDotEFn (const LVal & lval)
210 : _lval(lval) {}
211
212 template <int I, typename T>
213 inline T run (const T & val) const { return val * _lval.template eval<I>(); }
214 };
215
216 template <typename V>
217 class VRedDotVFn
218 {
219 private:
220 const typename V::VectorType _v;
221
222 public:
223 inline VRedDotVFn (const V & v)
224 : _v(v) {}
225
226 template <int I, typename T>
227 inline T run (const T & val) const { return val * _v.template element<I>(); }
228 };
229
230 // Expression builders:
231 // Argument
232 template <typename T, int D>
233 class VExpArg : public VBaseExpr<T, VExpArg< T, D >, D >
234 {
235 private:
236 const T _arg;
237
238 public:
239 inline VExpArg (const T & argg)
240 : _arg (argg) {}
241
242 template <int I>
243 inline const T eval () const
244 { return _arg; }
245 };
246
247 // Vector
248 template <typename V, int D>
249 class VExpArgVec : public VBaseExpr< typename V::Precision, VExpArgVec< typename V::VectorType, D >, D >
250 {
251 private:
252 const V _v;
253
254 public:
255 inline VExpArgVec (const V & v)
256 : _v (v) {}
257
258 template <int I>
259 inline const typename V::Precision eval () const
260 { return _v.template element<I>(); }
261 };
262
263 // Pairing
264 template <class LVal, typename RVal, class Op, int D>
265 class VExpr2 : public VBaseExpr<typename Op::ResultType, VExpr2< typename LVal::ExprType, typename RVal::ExprType, Op, D >, D >
266 {
267 private:
268 const LVal _lval;
269 const RVal _rval;
270
271 public:
272 typedef VExpr2< typename LVal::ExprType, typename RVal::ExprType, Op, D > ExprType;
273 typedef VBaseExpr<typename Op::ResultType, VExpr2< typename LVal::ExprType, typename RVal::ExprType, Op, D >, D > Base;
274 typedef typename Op::ResultType Precision;
275
276 inline VExpr2 (const LVal & lval, const RVal & rval)
277 : _lval(lval), _rval(rval) {}
278
279 template <int I>
280 inline const typename Op::ResultType eval () const
281 { return Op::template evalOp<I>(_lval, _rval); }
282
283 public:
284 // Scalar
285 inline const VExpr2< ExprType, VExpArg< Precision, D >, VOpScale< Precision, D >, D>
286 operator * (const Precision & n) const
287 { return VExpr2< ExprType, VExpArg< Precision, D >, VOpScale< Precision, D >, D> (*this, VExpArg< Precision, D > (n)); }
288
289 inline const VExpr2< ExprType, VExpArg< Precision, D >, VOpContract< Precision, D >, D>
290 operator / (const Precision & n) const
291 { return VExpr2< ExprType, VExpArg< Precision, D >, VOpContract< Precision, D >, D> (*this, VExpArg< Precision, D > (n)); }
292
293 inline const VExpr2< ExprType, VExpArg< Precision, D >, VOpAdd< Precision >, D>
294 operator + (const Precision & n) const
295 { return VExpr2< ExprType, VExpArg< Precision, D >, VOpAdd< Precision >, D> (*this, VExpArg< Precision, D > (n)); }
296
297 inline const VExpr2< ExprType, VExpArg< Precision, D >, VOpSub< Precision >, D>
298 operator - (const Precision & n) const
299 { return VExpr2< ExprType, VExpArg< Precision, D >, VOpSub< Precision >, D> (*this, VExpArg< Precision, D > (n)); }
300
301 // Vector
302 template <typename V>
303 inline const VExpr2< ExprType, VExpArgVec< typename V::VectorType, D >, VOpAdd< Precision >, D>
304 operator + (const V & v) const
305 { return VExpr2< ExprType, VExpArgVec< typename V::VectorType, D >, VOpAdd< Precision >, D> (*this, VExpArgVec< V, D > (v)); }
306
307 template <typename V>
308 inline const VExpr2< ExprType, VExpArgVec< typename V::VectorType, D >, VOpSub< Precision >, D>
309 operator - (const V & v) const
310 { return VExpr2< ExprType, VExpArgVec< typename V::VectorType, D >, VOpSub< Precision >, D> (*this, VExpArgVec< V, D > (v)); }
311
312 // Other expressions
313 template <typename RVal2>
314 inline const VExpr2< ExprType, typename RVal2::ExprType, VOpAdd< Precision >, D>
315 operator + (const RVal2 & rval) const
316 { return VExpr2< ExprType, typename RVal2::ExprType, VOpAdd< Precision >, D> (*this, rval); }
317
318 template <typename RVal2>
319 inline const VExpr2< ExprType, typename RVal2::ExprType, VOpSub< Precision >, D>
320 operator - (const RVal2 & rval) const
321 { return VExpr2< ExprType, typename RVal2::ExprType, VOpSub< Precision >, D> (*this, rval); }
322
323 template <typename RVal2>
324 inline bool operator == (const RVal2 & rval) const { return equals(rval); }
325 };
326
327 // Mono
328 template <typename RVal, typename Op, int D>
329 class VExpr1 : public VBaseExpr<typename Op::ResultType, VExpr1< typename RVal::ExprType, Op, D >, D >
330 {
331 private:
332 const RVal _rval;
333
334 public:
335 typedef VExpr1< typename RVal::ExprType, Op, D > ExprType;
336 typedef typename Op::ResultType Precision;
337
338 inline VExpr1 (const typename RVal::ExprType & rval)
339 : _rval(rval) {}
340
341 template <int I>
342 inline const typename Op::ResultType eval () const
343 { return Op::template evalOp<I>(_rval); }
344
345 public:
346
347 // Scalar
348 inline const VExpr2< ExprType, VExpArg< Precision, D >, VOpScale< Precision, D >, D>
349 operator * (const Precision & n) const
350 { return VExpr2< ExprType, VExpArg< Precision, D >, VOpScale< Precision, D >, D> (*this, VExpArg< Precision, D > (n)); }
351
352 inline const VExpr2< ExprType, VExpArg< Precision, D >, VOpContract< Precision, D >, D>
353 operator / (const Precision & n) const
354 { return VExpr2< ExprType, VExpArg< Precision, D >, VOpContract< Precision, D >, D> (*this, VExpArg< Precision, D > (n)); }
355
356 inline const VExpr2< ExprType, VExpArg< Precision, D >, VOpAdd< Precision >, D>
357 operator + (const Precision & n) const
358 { return VExpr2< ExprType, VExpArg< Precision, D >, VOpAdd< Precision >, D> (*this, VExpArg< Precision, D > (n)); }
359
360 inline const VExpr2< ExprType, VExpArg< Precision, D >, VOpSub< Precision >, D>
361 operator - (const Precision & n) const
362 { return VExpr2< ExprType, VExpArg< Precision, D >, VOpSub< Precision >, D> (*this, VExpArg< Precision, D > (n)); }
363
364 // Vector
365 template <typename V>
366 inline const VExpr2< ExprType, VExpArgVec< typename V::VectorType, D >, VOpAdd< Precision >, D>
367 operator + (const V & v) const
368 { return VExpr2< ExprType, VExpArgVec< typename V::VectorType, D >, VOpAdd< Precision >, D> (*this, VExpArgVec< V, D > (v)); }
369
370 template <typename V>
371 inline const VExpr2< ExprType, VExpArgVec< typename V::VectorType, D >, VOpSub< Precision >, D>
372 operator - (const V & v) const
373 { return VExpr2< ExprType, VExpArgVec< typename V::VectorType, D >, VOpSub< Precision >, D> (*this, VExpArgVec< V, D > (v)); }
374
375 // Other expressions
376 template <typename RVal2>
377 inline const VExpr2< ExprType, typename RVal2::ExprType, VOpAdd< Precision >, D>
378 operator + (const RVal2 & rval) const
379 { return VExpr2< ExprType, typename RVal2::ExprType, VOpAdd< Precision >, D> (*this, rval); }
380
381 template <typename RVal2>
382 inline const VExpr2< ExprType, typename RVal2::ExprType, VOpSub< Precision >, D>
383 operator - (const RVal2 & rval) const
384 { return VExpr2< ExprType, typename RVal2::ExprType, VOpSub< Precision >, D> (*this, rval); }
385
386 template <typename RVal2>
387 inline bool operator == (const RVal2 & rval) const { return equals(rval); }
388 };
389
390 template <typename T, typename V, typename Desc, int D >
391 struct VectorDefsTag : public VExprTagDim< T, D >
392 {
393 typedef VectorDefsTag<T , V, Desc, D> VectorTagType;
394 typedef V VectorType;
395 typedef VExpArg< T, D > Arg;
396 typedef VExpArgVec< V, D > ArgVec;
397
398 template <int I>
399 inline const T element () const
400 {
401#ifdef _DEBUG
402 int i = I;
403 assert(I < D);
404#endif
405 return reinterpret_cast <const T *> (static_cast <const Desc *> (this)) [I];
406 }
407 };
408
409 template <typename T, typename V, typename Desc, int D >
410 class ArithmeticVector : public VectorDefsTag< T, V, Desc, D >
411 {
412 public:
413 //typedef ArithmeticVector< T, V, Desc, D > ArithmeticVectorType;
414 typedef VectorDefsTag< T, V, Desc, D > Base;
415 typedef VExprTagDim< T, D > ExprTag;
416
417 using typename Base::Arg;
418 using typename Base::ArgVec;
419 using typename ExprTag::OpNegate;
420 using typename ExprTag::OpScale;
421 using typename ExprTag::OpContract;
422 using typename ExprTag::OpAdd;
423 using typename ExprTag::OpSub;
424
425 typedef VExpr1< ArgVec, OpNegate, D > ExpNeg;
426 typedef VExpr2< ArgVec, Arg, OpScale, D > ExpScale;
427 typedef VExpr2< ArgVec, Arg, OpContract, D > ExpContract;
428 typedef VExpr2< ArgVec, Arg, OpAdd, D > ExpAddScalar;
429 typedef VExpr2< ArgVec, ArgVec, OpAdd, D > ExpAddVector;
430 typedef VExpr2< ArgVec, Arg, OpSub, D > ExpSubScalar;
431 typedef VExpr2< ArgVec, ArgVec, OpSub, D > ExpSubVector;
432
433 template <typename Expr>
434 struct ExpOpExp
435 {
436 typedef VExpr2< ArgVec, typename Expr::ExprType, OpAdd, D > Add;
437 typedef VExpr2< ArgVec, typename Expr::ExprType, OpSub, D > Sub;
438 };
439
440 // Negation
441 inline ExpNeg neg () const
442 { return ExpNeg (ArgVec(* static_cast <const V *> (this))); }
443
444 // Scaling
445 inline ExpScale scale (const T & n) const
446 { return ExpScale (ArgVec(* static_cast <const V *> (this)), Arg(n)); }
447
448 // Contraction
449 inline ExpContract contract (const T & n) const
450 { return ExpContract (ArgVec(* static_cast <const V *> (this)), Arg(n)); }
451
452 // Add scalar
453 inline ExpAddScalar add (const T & n) const
454 { return ExpAddScalar (ArgVec(* static_cast <const V *> (this)), Arg(n)); }
455
456 // Add vector
457 inline ExpAddVector add (const V & v) const
458 { return ExpAddVector (ArgVec(* static_cast <const V *> (this)), ArgVec(v)); }
459
460 // Subtract scalar
461 inline ExpSubScalar sub (const T & n) const
462 { return ExpSubScalar (ArgVec(* static_cast <const V *> (this)), Arg(n)); }
463
464 // Subtract vector
465 inline ExpSubVector sub (const V & v) const
466 { return ExpSubVector (ArgVec(* static_cast <const V *> (this)), ArgVec(v)); }
467 };
468
469 template <typename T, typename V>
470 struct DescV2 : public ArithmeticVector< T, V, DescV2< T, V >, 2 >
471 {
472 typedef ArithmeticVector< T, V, DescV2< T, V >, 2 > Base;
473
474 T x, y;
475
476 inline DescV2 ()
477 : x(0), y(0) {}
478 inline DescV2 (const T x, const T y)
479 : x(x), y(y) {}
480
481 template <typename E, typename Vv>
482 inline DescV2 (const DescV2 <E, Vv> & copy)
483 : x(static_cast <T> (copy.x)),
484 y(static_cast <T> (copy.y)) {}
485
486 template <typename Exp>
487 inline V & operator = (const Exp & expr) const
488 {
489 x = expr.template eval<0>();
490 y = expr.template eval<1>();
491 return * static_cast <const V *> (this);
492 }
493 };
494
495 template <typename T, typename V>
496 struct DescV3 : public ArithmeticVector< T, V, DescV3< T, V >, 3 >
497 {
498 typedef ArithmeticVector< T, V, DescV3< T, V >, 3 > Base;
499
500 T x, y, z;
501
502 inline DescV3 ()
503 : x(0), y(0), z(0) {}
504 inline DescV3 (const T x, const T y, const T z)
505 : x(x), y(y), z(z) {}
506
507 template <typename E, typename Vv>
508 inline DescV3 (const DescV3 <E, Vv> & copy)
509 : x(static_cast <T> (copy.x)),
510 y(static_cast <T> (copy.y)),
511 z(static_cast <T> (copy.z)) {}
512
513 template <typename Exp>
514 inline V & operator = (const Exp & expr) const
515 {
516 x = expr.template eval<0>();
517 y = expr.template eval<1>();
518 z = expr.template eval<2>();
519 return * static_cast <const V *> (this);
520 }
521 };
522
523 template <typename T, typename V>
524 struct DescQ : public VectorDefsTag< T, V, DescQ< T, V >, 3 >
525 {
526 typedef VectorDefsTag< T, V, DescQ< T, V >, 3 > Base;
527
528 T x, y, z, w;
529
530 inline DescQ ()
531 : x(0), y(0), z(0), w(0) {}
532 inline DescQ (const T x, const T y, const T z, const T w)
533 : x(x), y(y), z(z), w(w) {}
534
535 template <typename Vv>
536 inline DescQ (const DescV3< T, Vv> & v)
537 : x(v.x), y(v.y), z(v.z), w(0) {}
538
539 template <typename Vv>
540 inline DescQ (const T s, const DescV3< T, Vv> & v)
541 : x(v.x), y(v.y), z(v.z), w(s) {}
542
543 template <typename E>
544 inline DescQ (const DescQ <E, V> & copy)
545 : x(static_cast <T> (copy.x)),
546 y(static_cast <T> (copy.y)),
547 z(static_cast <T> (copy.z)),
548 w(static_cast <T> (copy.w)) {}
549 };
550
551 template <typename T> class vector3D;
552 template <typename T> class SphericalCoords;
553 template <typename T> class PolarCoords;
554 template <typename T> class CylindricalCoords;
555
556 // *** Vector 2D class template ***
557 template <typename T>
558 class vector2D : public DescV2< T, vector2D< T > >, private MathOverrides
559 {
560 private:
561 typedef DescV2< T, vector2D< T > > Base;
562
563 public:
564 template <typename Exp>
565 inline vector2D (const Exp & expr, const typename Exp::ExprType * dummy = 0)
566 : DescV2< T, vector2D< T > >(expr.template eval<0>(), expr.template eval<1>()) {}
567
568 inline vector2D ()
569 : DescV2< T, vector2D< T > >() {}
570 inline vector2D (const T x, const T y)
571 : DescV2< T, vector2D< T > >(x, y) {}
572 template <typename E>
573 inline vector2D (const vector2D <E> & other)
574 : DescV2< T, vector2D< T > >(other) {}
575
576 // BEGIN ArithmeticVector INHERITANCE {
577 // Binary expression mappings
578 template <typename Expr>
579 inline typename Base::template ExpOpExp<typename Expr::ExprType>::Add operator + (const Expr & rval) const
580 { return typename Base::template ExpOpExp<typename Expr::ExprType>::Add (typename Base::ArgVec(*this), rval); }
581
582 template <typename Expr>
583 inline typename Base::template ExpOpExp<typename Expr::ExprType>::Sub operator - (const Expr & rval) const
584 { return typename Base::template ExpOpExp<typename Expr::ExprType>::Sub (typename Base::ArgVec(*this), rval); }
585
586 template <typename Expr>
587 inline typename enable_ifb< has_def< typename Expr::ExprType >::value, bool >::Type operator == (const Expr & rval) const
588 { return rval.equals(*this); }
589
590 // Import dependent types
591 using typename Base::ExpNeg;
592 using typename Base::ExpScale;
593 using typename Base::ExpContract;
594 using typename Base::ExpAddScalar;
595 using typename Base::ExpAddVector;
596 using typename Base::ExpSubScalar;
597 using typename Base::ExpSubVector;
598
599 using typename Base::VectorType;
600
601 // Basic binary mappings
602 inline ExpNeg operator - () const { return this->neg(); }
603 inline ExpScale operator * (const T & n) const { return this->scale(n); }
604 inline ExpContract operator / (const T & n) const { return this->contract(n); }
605 inline ExpAddScalar operator + (const T & n) const { return this->add(n); }
606 inline ExpAddVector operator + (const VectorType & v) const { return this->add(v); }
607 inline ExpSubScalar operator - (const T & n) const { return this->sub(n); }
608 inline ExpSubVector operator - (const VectorType & v) const { return this->sub(v); }
609 // } END ArithmeticVector INHERITANCE
610
611 // Dot-product
612 inline T operator * (const vector2D<T> & v) const
613 { return this->x * v.x + this->y * v.y; }
614
615 // Maps to 2D-cross product below
616 inline const vector2D <T> operator & (const vector3D <T> & v) const
617 { return operator & (v.z); }
618
619 // 2D-cross product (perpendicular)
620 inline const vector2D <T> operator & (const T s) const
621 { return vector2D<T> (this->y * s, -this->x * s); }
622
623 // Accumulate scalar
624 inline const vector2D <T> & operator += (const T & n)
625 {
626 this->x += n;
627 this->y += n;
628 return *this;
629 }
630
631 // Accumulate vector
632 inline const vector2D <T> & operator += (const vector2D<T> & v)
633 {
634 this->x += v.x;
635 this->y += v.y;
636
637 return *this;
638 }
639
640 // Decrement vector
641 inline const vector2D <T> & operator -= (const vector2D<T> & v)
642 {
643 this->x -= v.x;
644 this->y -= v.y;
645
646 return *this;
647 }
648 inline T magSQ () const { return static_cast <T> (SQ(this->x) + SQ(this->y)); }
649 inline T magnitude () const { return static_cast <T> (sqrt(magSQ())); }
650 inline vector2D <T> normalize () const
651 {
652 const T m = magnitude();
653
654 if (m == 0)
655 return vector2D <T> (0,0);
656 else
657 return vector2D<T> (
658 this->x / m,
659 this->y / m
660 );
661 }
662
663 template <typename J>
664 inline bool operator != (const vector2D <J> & other)
665 {
666 return this->x != other.x || this->y != other.y;
667 }
668 template <typename J>
669 inline bool operator == (const vector2D <J> & other)
670 {
671 return this->x == other.x && this->y == other.y;
672 }
673
674 // Assignment
675 template <typename J>
676 inline vector2D <T> & operator = (const vector2D <J> & copy)
677 {
678 this->x = static_cast <T> (copy.x);
679 this->y = static_cast <T> (copy.y);
680
681 return *this;
682 }
683
684 // Scale vector
685 inline vector2D <T> & operator *= (const T factor)
686 {
687 this->x *= factor;
688 this->y *= factor;
689 return *this;
690 }
691
692 // Contract vector
693 inline vector3D <T> & operator /= (const T factor)
694 {
695 this->x /= factor;
696 this->y /= factor;
697 return *this;
698 }
699
700 // Length comparisons
701 inline bool operator < (const T len) const
702 { return magSQ() < SQ(len); }
703 inline bool operator > (const T len) const
704 { return magSQ() > SQ(len); }
705 inline bool operator <= (const T len) const
706 { return magSQ() <= SQ(len); }
707 inline bool operator >= (const T len) const
708 { return magSQ() >= SQ(len); }
709
710 inline bool operator == (const vector2D< T > & b) const
711 { return this->x == b.x && this->y == b.y; }
712 };
713
714 // *** Vector 3D class template ***
715 template <typename T>
716 class vector3D : public DescV3< T, vector3D< T > >, private MathOverrides
717 {
718 private:
719 typedef DescV3< T, vector3D< T > > Base;
720
721 public:
722 template <typename Exp>
723 inline vector3D (const Exp & expr, const typename Exp::ExprType * dummy = 0)
724 : DescV3< T, vector3D< T > >(expr.template eval<0>(), expr.template eval<1>(), expr.template eval<2>()) {}
725
726 template <typename E>
727 inline vector3D (const SphericalCoords< E > & sphc)
728 : DescV3< T, vector3D< T > >(
729 static_cast <T> (sphc.r * sin(sphc.zenith) * cos(sphc.azimuth)),
730 static_cast <T> (sphc.r * sin(sphc.zenith) * sin(sphc.azimuth)),
731 static_cast <T> (sphc.r * cos(sphc.zenith))
732 ) {}
733
734 inline vector3D ()
735 : DescV3< T, vector3D< T > >() {}
736 inline vector3D (const T x, const T y, const T z)
737 : DescV3< T, vector3D< T > >(x, y, z) {}
738
739 template <typename E>
740 inline vector3D (const vector3D <E> & other)
741 : DescV3< T, vector3D< T > > (other) {}
742
743 // Spherical adapter
744 inline vector3D operator = (const SphericalCoords <T> & sphc)
745 {
746 this->x = sphc.r * sin(sphc.zenith) * cos(sphc.azimuth);
747 this->y = sphc.r * sin(sphc.zenith) * sin(sphc.azimuth);
748 this->z = sphc.r * cos(sphc.zenith);
749 }
750 template <typename E>
751 inline vector3D<E> & cast ()
752 {
753 return vector3D <E> (
754 static_cast<T> (this->other.x),
755 static_cast<T> (this->other.y),
756 static_cast<T> (this->other.z)
757 );
758 }
759
760 // BEGIN ArithmeticVector INHERITANCE {
761 // Binary expression mappings
762 template <typename Expr>
763 inline typename Base::template ExpOpExp<typename Expr::ExprType>::Add operator + (const Expr & rval) const
764 { return typename Base::template ExpOpExp<typename Expr::ExprType>::Add (typename Base::ArgVec(*this), rval); }
765
766 template <typename Expr>
767 inline typename Base::template ExpOpExp<typename Expr::ExprType>::Sub operator - (const Expr & rval) const
768 { return typename Base::template ExpOpExp<typename Expr::ExprType>::Sub (typename Base::ArgVec(*this), rval); }
769
770 template <typename Expr>
771 inline typename enable_ifb< has_def< typename Expr::ExprType >::value, bool >::Type operator == (const Expr & rval) const
772 { return rval.equals(*this); }
773
774 // Import dependent types
775 using typename Base::ExpNeg;
776 using typename Base::ExpScale;
777 using typename Base::ExpContract;
778 using typename Base::ExpAddScalar;
779 using typename Base::ExpAddVector;
780 using typename Base::ExpSubScalar;
781 using typename Base::ExpSubVector;
782
783 using typename Base::VectorType;
784
785 // Basic binary mappings
786 inline ExpNeg operator - () const { return this->neg(); }
787 inline ExpScale operator * (const T & n) const { return this->scale(n); }
788 inline ExpContract operator / (const T & n) const { return this->contract(n); }
789 inline ExpAddScalar operator + (const T & n) const { return this->add(n); }
790 inline ExpAddVector operator + (const VectorType & v) const { return this->add(v); }
791 inline ExpSubScalar operator - (const T & n) const { return this->sub(n); }
792 inline ExpSubVector operator - (const VectorType & v) const { return this->sub(v); }
793 // } END ArithmeticVector INHERITANCE
794
795 // Shortcut mappings
796 inline T operator * (const vector3D<T> & v) const { return dot(v); }
797 inline bool operator || (const vector3D<T> & v) const { return orthogonal(v); }
798 inline const vector3D <T> operator & (const vector3D <T> & v) const { return cross(v); }
799
800 inline bool operator == (const vector3D< T > & b) const { return equals(b); }
801 template <typename J> inline bool operator != (const vector3D <J> & other) const { return nequal(other); }
802 template <typename J> inline bool operator == (const vector3D <J> & other) const { return equals(other); }
803
804 inline bool operator < (const T len) const { return lesser(len); }
805 inline bool operator > (const T len) const { return greater(len); }
806 inline bool operator <= (const T len) const { return lessequal(len); }
807 inline bool operator >= (const T len) const { return greatequal(len); }
808
809 // Vector2 adapter
810 inline operator vector2D <T> () const { return vector2D <T> (this->x, this->y); }
811
812 // 3D-specific operations
813 inline T dot (const vector3D<T> & v) const
814 { return this->x * v.x + this->y * v.y + this->z * v.z; }
815 inline bool orthogonal (const vector3D<T> & v) const
816 { return Matrix3D(vector3D<T>(1,1,1), *this, v) == 0; }
817 inline const vector3D <T> cross (const vector3D <T> & v) const
818 {
819 return vector3D <T> (
820 this->y*v.z - this->z*v.y,
821 this->z*v.x - this->x*v.z,
822 this->x*v.y - this->y*v.x
823 );
824 }
825
826 // Length
827 inline T magSQ () const { return static_cast <T> (SQ(this->x) + SQ(this->y) + SQ(this->z)); }
828 inline T magnitude () const { return static_cast <T> (sqrt(magSQ())); }
829 inline vector3D <T> normalize () const
830 {
831 const T m = magnitude();
832
833 if (m == 0)
834 return vector3D <T> (0,0,0);
835 else
836 return vector3D<T> (
837 this->x / m,
838 this->y / m,
839 this->z / m
840 );
841 }
842
843 // Equivalence
844 inline bool equals (const vector3D< T > & b) const
845 { return this->x == b.x && this->y == b.y && this->z == b.z; }
846
847 template <typename J>
848 inline bool nequal (const vector3D <J> & other) const
849 { return this->x != other.x || this->y != other.y || this->z != other.z; }
850
851 template <typename J>
852 inline bool equals (const vector3D <J> & other) const
853 { return this->x == other.x && this->y == other.y && this->z == other.z; }
854
855 // Length comparisons
856 inline bool lesser (const T len) const
857 { return magSQ() < SQ(len); }
858 inline bool greater (const T len) const
859 { return magSQ() > SQ(len); }
860 inline bool lessequal (const T len) const
861 { return magSQ() <= SQ(len); }
862 inline bool greatequal (const T len) const
863 { return magSQ() >= SQ(len); }
864
865 inline vector2D <T> v2D() const { return vector2D <T> (this->x, this->y); }
866
867 // Accumulate scalar
868 inline const vector3D <T> & operator += (const T & n)
869 {
870 this->x += n;
871 this->y += n;
872 this->z += n;
873 }
874
875 // Accumulate vector
876 inline const vector3D <T> & operator += (const vector3D<T> & v)
877 {
878 this->x += v.x;
879 this->y += v.y;
880 this->z += v.z;
881
882 return *this;
883 }
884
885 // Decrement by vector
886 inline const vector3D <T> & operator -= (const vector3D<T> & v)
887 {
888 this->x -= v.x;
889 this->y -= v.y;
890 this->z -= v.z;
891
892 return *this;
893 }
894
895 // Assignment
896 template <typename J>
897 inline vector3D <T> & operator = (const vector3D <J> & copy)
898 {
899 this->x = static_cast <T> (copy.x);
900 this->y = static_cast <T> (copy.y);
901 this->z = static_cast <T> (copy.z);
902
903 return *this;
904 }
905
906 // Scale vector
907 inline vector3D <T> & operator *= (const T factor)
908 {
909 this->x *= factor;
910 this->y *= factor;
911 this->z *= factor;
912 return *this;
913 }
914
915 // Contract vector
916 inline vector3D <T> & operator /= (const T factor)
917 {
918 this->x /= factor;
919 this->y /= factor;
920 this->z /= factor;
921 return *this;
922 }
923 };
924 // Swap-operator mappings
925 template <typename V>
926 inline typename V::ExpScale operator * (const typename V::Precision f, const V & v)
927 { return v.scale(f); }
928
929 template <typename T, typename RVal>
930 inline VExpr2< typename RVal::ExprType, VExpArg< T, RVal::DIMENSION >, VOpScale< T, RVal::DIMENSION >, RVal::DIMENSION >
931 operator * (const T f, const RVal & rval)
932 { return VExpr2< typename RVal::ExprType, VExpArg< T, RVal::DIMENSION >, VOpScale< T, RVal::DIMENSION >, RVal::DIMENSION > (rval, VExpArg< T, RVal::DIMENSION > (f)); }
933
934 // *** Generic mappings ***
935 // ** Expressions **
936 template <typename E>
937 inline typename E::ExprType::Precision
938 MAG(const E & expr) { return MathOverrides::sqrt(expr.template reduce<VRedSqFn>()); }
939
940 template <typename E>
941 inline typename E::ExprType::Precision
942 MAGSQ(const E & expr) { return expr.template reduce<VRedSqFn>(); }
943
944 template <typename V>
945 inline typename V::VectorType::Precision
946 DOT(const V & a, const V & b) { return a * b; }
947
948 template <typename V, typename Expr>
949 inline typename enable_ifb< andb< has_def< typename V::VectorType>::value, has_def< typename Expr::ExprType >::value >::value, typename V::Precision >::Type
950 DOT(const V & a, const Expr & b) { return b.template reduce< VRedDotVFn< V > > (VRedDotVFn< V >(a)); }
951
952 template <typename Expr, typename V>
953 inline typename enable_ifb< andb< has_def< typename Expr::ExprType >::value, has_def< typename V::VectorType >::value >::value, typename Expr::Precision >::Type
954 DOT(const Expr & a, const V & b) { return a.template reduce< VRedDotVFn< V > > (VRedDotVFn< V >(b)); }
955
956 template <typename LVal, typename RVal>
957 inline typename enable_ifb< andb< has_def< typename LVal::ExprType >::value, has_def< typename RVal::ExprType >::value >::value, typename LVal::Precision >::Type
958 DOT(const LVal & a, const RVal & b) { return a.template reduce<VRedDotEFn< RVal > > (VRedDotEFn< RVal > (b)); }
959
960 template <typename A, typename B>
961 inline typename enable_ifb< andb< has_def<typename A::VExprTagType>::value, has_def<typename B::VExprTagType>::value >::value, typename A::Precision>::Type
962 ANGLE (const A & a, const B & b) { return acosf (DOT(a,b) / (MAG(a) * MAG(b))); }
963
964 template <typename V>
965 inline typename V::VectorType::Precision
966 MAG(const V & v) { return v.magnitude(); }
967
968 template <typename V>
969 inline typename V::VectorType::Precision
970 MAGSQ(const V & v) { return v.magSQ(); }
971
972 template <typename V>
973 inline typename V::VectorType
974 U(const V & a) { return a.normalize(); }
975
976 template <typename Expr, typename V>
977 inline typename enable_ifb< andb< has_def< typename Expr::ExprType >::value, has_def< typename V::VectorType >::value >::value >::Type
978 U(const Expr & a) { return V(a.normalize()); }
979
980 template < typename TNew, typename E >
981 inline VExpr1< E, VOpCast< typename E::ExprType::Precision, TNew >, E::ExprType::DIMENSION >
982 CAST(const E & expr) { return VExpr1< E, VOpCast< typename E::ExprType::Precision, TNew >, E::ExprType::DIMENSION > (expr); }
983
984 typedef vector2D<float> vector2Df;
985 typedef vector2D<double> vector2Dd;
986 typedef vector3D<float> vector3Df;
987 typedef vector3D<double> vector3Dd;
988
989 // Spherical coordinates, definition taken from http://mathworld.wolfram.com/SphericalCoordinates.html
990 template <typename T>
991 class SphericalCoords : private MathOverrides
992 {
993 public:
994 typedef T Precision;
995
996 T r;
997 float
998 zenith, // The circle with a vertical bar through it, thus represents the polar angle
999 azimuth; // Looks like "theta", circle with horizontal bar through it, thus represents the angle in the x/y plane
1000
1001 inline SphericalCoords () : r(0), zenith(0), azimuth(0) {}
1002
1003 inline SphericalCoords (const vector3D <T> & v)
1004 : r(static_cast< T > (sqrt(SQ(v.x) + SQ(v.y) + SQ(v.z)))),
1005 zenith(atan2f(static_cast< float > (v.y), static_cast< float > (v.x))),
1006 azimuth(0)
1007 {
1008 if (r != 0)
1009 azimuth = acos(static_cast< float > (v.z) / static_cast< float > (r));
1010 }
1011 inline SphericalCoords (const T r, const float zenith, const float azimuth)
1012 : r(r), zenith(zenith), azimuth(azimuth) {}
1013
1014 inline SphericalCoords <T> & operator = (const vector3D <T> & v)
1015 {
1016 r = sqrt(SQ(v.x) + SQ(v.y) + SQ(v.z));
1017
1018 azimuth = atan2(
1019 static_cast< float > (v.y),
1020 static_cast< float > (v.x)
1021 );
1022
1023 if (r != 0)
1024 zenith = acos(
1025 static_cast< float > (v.z)
1026 /
1027 static_cast< float > (r)
1028 );
1029
1030 return *this;
1031 }
1032 inline const vector3D <T> operator + (const vector3D <T> rval) const { return static_cast <const vector3D<T> > (*this) + rval; }
1033 inline operator const vector3D <T> () const
1034 {
1035 return vector3D <T>
1036 (
1037 static_cast< T > (static_cast< float > (r) * sin(zenith) * cos(azimuth)),
1038 static_cast< T > (static_cast< float > (r) * sin(zenith) * sin(azimuth)),
1039 static_cast< T > (static_cast< float > (r) * cos(zenith))
1040 );
1041 }
1042 inline SphericalCoords <T> & operator = (const CylindricalCoords <T> & c)
1043 {
1044 r = static_cast< T > (sqrt(SQ(c.p) + SQ(c.z)));
1045 zenith = atan2(c.p, c.z);
1046 azimuth = c.azimuth;
1047 return *this;
1048 }
1049 inline operator const CylindricalCoords <T> () const
1050 {
1051 return CylindricalCoords <T>
1052 (
1053 static_cast< T > (static_cast< float > (r) * sin(zenith)),
1054 azimuth,
1055 static_cast< T > (static_cast< float > (r) * cos(zenith))
1056 );
1057 }
1058 };
1059
1060 template <typename T>
1061 class PolarCoords
1062 {
1063 private:
1064 inline PolarCoords <T> & assign (const T x, const T y)
1065 {
1066 p = sqrt(SQ(x) + SQ(y));
1067 azimuth = atan2(y, x);
1068 return *this;
1069 }
1070
1071 public:
1072 T p;
1073 float azimuth;
1074
1075 inline PolarCoords () : p(0), azimuth(0) {}
1076 inline PolarCoords (const vector2D <T> & v) : p(sqrt(SQ(v.x)) + sqrt(SQ(v.y))), azimuth(atan2(v.y, v.x)) {}
1077
1078 inline PolarCoords (const T p, const float azimuth)
1079 : p(p), azimuth(azimuth) {}
1080
1081 inline PolarCoords <T> & operator = (const vector2D <T> & v)
1082 { return assign(v.x, v.y); }
1083 inline PolarCoords <T> & operator = (const vector3D <T> & v)
1084 { return assign(v.x, v.y); }
1085
1086 template <typename J>
1087 inline operator const vector2D <J> () const
1088 {
1089 return vector2D <J> (
1090 static_cast <J> (static_cast <double> (p) * cos(azimuth)),
1091 static_cast <J> (static_cast <double> (p) * sin(azimuth))
1092 );
1093 }
1094 inline PolarCoords <T> & operator = (const SphericalCoords <T> & s)
1095 {
1096 p = s.r * sin(s.zenith);
1097 azimuth = s.azimuth;
1098 return *this;
1099 }
1100 };
1101
1102 template <typename T>
1103 class CylindricalCoords : public PolarCoords <T>
1104 {
1105 public:
1106 T z;
1107
1108 inline CylindricalCoords () : PolarCoords<T> (0, 0), z(0) {}
1109 inline CylindricalCoords (const vector3D <T> & v) : PolarCoords<T> (v), z(v.z) {}
1110
1111 inline CylindricalCoords (const T p, const float azimuth, const T z)
1112 : PolarCoords<T>(p, azimuth), z(z) {}
1113
1114 inline CylindricalCoords <T> & operator = (const vector3D <T> & v)
1115 {
1116 PolarCoords <T>::operator = (v);
1117 z = v.z;
1118 return *this;
1119 }
1120 inline CylindricalCoords <T> & operator = (const SphericalCoords <T> & s)
1121 {
1122 PolarCoords <T>::operator = (s);
1123 z = s.r * cos(s.zenith);
1124 return *this;
1125 }
1126 template <typename J>
1127 inline operator const vector3D <J> () const
1128 {
1129 return vector3D <J> (
1130 static_cast <J> (static_cast <double> (this->p) * cos(this->azimuth)),
1131 static_cast <J> (static_cast <double> (this->p) * sin(this->azimuth)),
1132 static_cast <J> (z)
1133 );
1134 }
1135 inline operator const SphericalCoords <T> () const
1136 {
1137 return SphericalCoords <T>
1138 (
1139 sqrt(SQ(this->p) + SQ(z)),
1140 atan2(this->p, z),
1141 this->azimuth
1142 );
1143 }
1144 };
1145
1146 template <typename T>
1147 class Matrix3D
1148 {
1149 public:
1150 typedef T Precision;
1151
1152 T values [3][3];
1153
1154 inline Matrix3D (const T a1, const T b2, const T c3)
1155 {
1156 values[0][0] = a1; values[0][1] = 0; values[0][2] = 0;
1157 values[1][0] = 0; values[1][1] = b2; values[1][2] = 0;
1158 values[2][0] = 0; values[2][1] = 0; values[2][2] = c3;
1159 }
1160 inline Matrix3D (
1161 const T a1, const T b1, const T c1,
1162 const T a2, const T b2, const T c2,
1163 const T a3, const T b3, const T c3
1164 )
1165 {
1166 values[0][0] = a1; values[0][1] = b1; values[0][2] = c1;
1167 values[1][0] = a2; values[1][1] = b2; values[1][2] = c2;
1168 values[2][0] = a3; values[2][1] = b3; values[2][2] = c3;
1169 }
1170 inline Matrix3D (
1171 const vector3D<T> & a, vector3D<T> b, vector3D<T> c
1172 )
1173 {
1174 values[0][0] = a.x; values[0][1] = a.y; values[0][2] = a.z;
1175 values[1][0] = b.x; values[1][1] = b.y; values[1][2] = b.z;
1176 values[2][0] = c.x; values[2][1] = c.y; values[2][2] = c.z;
1177 }
1178
1179 inline operator T ()
1180 {
1181 return
1182 values[0][0] * values[1][1] * values[2][2] -
1183 values[0][0] * values[1][2] * values[2][1] +
1184 values[0][1] * values[1][2] * values[2][0] -
1185 values[0][1] * values[1][0] * values[2][2] +
1186 values[0][2] * values[1][0] * values[2][1] -
1187 values[0][2] * values[1][1] * values[2][0];
1188 }
1189 inline mars::vector3D <T> operator * (const mars::vector3D <T> v) const
1190 {
1191 return vector3D<T> (
1192 v.x * values[0][0] + v.y * values[0][1] + v.z * values[0][2],
1193 v.x * values[1][0] + v.y * values[1][1] + v.z * values[1][2],
1194 v.x * values[2][0] + v.y * values[2][1] + v.z * values[2][2]
1195 );
1196 }
1197 inline mars::Matrix3D <T> transpose () const
1198 {
1199 return Matrix3D (
1200 values[0][0], values[1][0], values[2][0],
1201 values[0][1], values[1][1], values[2][1],
1202 values[0][2], values[1][2], values[2][2]
1203 );
1204 }
1205 static inline mars::Matrix3D <T> rotateZ (float theta)
1206 {
1207 return Matrix3D <T> (
1208 cos(theta), -sin(theta), 0,
1209 sin(theta), cos(theta), 0,
1210 0, 0, 1
1211 );
1212 }
1213 };
1214
1215 template <typename T>
1216 class Quaternion : public DescQ< T, Quaternion< T > >
1217 {
1218 public:
1219 inline Quaternion (const T w, const T x, const T y, const T z)
1220 : DescQ< T, Quaternion< T > >(x, y, z, w) {}
1221 inline Quaternion (const vector3D <T> & vec)
1222 : DescQ< T, Quaternion< T > >(vec) {}
1223 inline Quaternion (const T s, const vector3D <T> & v)
1224 : DescQ< T, Quaternion< T > >(s, v) {}
1225 inline Quaternion ()
1226 : DescQ< T, Quaternion< T > >() {}
1227
1228 static
1229 inline Quaternion <T>
1230 rotation (const T theta, const vector3D <T> & u)
1231 {
1232 return Quaternion <T> (
1233 cos(theta / 2),
1234 u * sin(theta / 2)
1235 );
1236 }
1237
1238 inline const Quaternion <T> operator * (const Quaternion <T> & q) const
1239 {
1240 const T s1 = this->w, s2 = q.w;
1241 const vector3D <T> & v1 = this->v(), v2 = q.v();
1242
1243 return Quaternion <T> (
1244 s1 * s2 - (v1 * v2),
1245 (s1 * v2) + (s2 * v1) + (v1 & v2)
1246 );
1247 }
1248 inline const Quaternion <T> operator - () const
1249 {
1250 return Quaternion <T> (this->w, -this->x, -this->y, -this->z);
1251 }
1252 inline const Quaternion <T> operator ! () const
1253 {
1254 const Quaternion <T> qq = operator - ();
1255 return qq / (operator * (qq));
1256 }
1257 inline const Quaternion <T> operator / (const Quaternion <T> & q) const
1258 {
1259 const T divisor = (q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w);
1260 return Quaternion <T> (
1261 (q.x*this->x + q.y*this->y + q.z*this->z + q.w*this->w) / divisor,
1262 (this->x*q.w - this->y*q.z + this->z*q.y - this->w*q.x) / divisor,
1263 (this->x*q.z + this->y*q.w - this->z*q.x - this->w*q.y) / divisor,
1264 (this->y*q.x + this->z*q.w - this->x*q.y - this->w*q.z) / divisor
1265 );
1266 }
1267 inline const Quaternion <T> operator & (const Quaternion <T> & q) const
1268 {
1269 return operator * (q * operator ! ());
1270 }
1271 inline operator vector3D <T> () const
1272 {
1273 return this->v();
1274 }
1275 };
1276}
1277
1278#endif
Note: See TracBrowser for help on using the repository browser.