070c1185af80a553a5ea5e71e2cb0de18b1c9e80
[blender.git] / extern / ode / dist / include / ode / odemath.h
1 /*************************************************************************
2  *                                                                       *
3  * Open Dynamics Engine, Copyright (C) 2001,2002 Russell L. Smith.       *
4  * All rights reserved.  Email: russ@q12.org   Web: www.q12.org          *
5  *                                                                       *
6  * This library is free software; you can redistribute it and/or         *
7  * modify it under the terms of EITHER:                                  *
8  *   (1) The GNU Lesser General Public License as published by the Free  *
9  *       Software Foundation; either version 2.1 of the License, or (at  *
10  *       your option) any later version. The text of the GNU Lesser      *
11  *       General Public License is included with this library in the     *
12  *       file LICENSE.TXT.                                               *
13  *   (2) The BSD-style license that is included with this library in     *
14  *       the file LICENSE-BSD.TXT.                                       *
15  *                                                                       *
16  * This library is distributed in the hope that it will be useful,       *
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the files    *
19  * LICENSE.TXT and LICENSE-BSD.TXT for more details.                     *
20  *                                                                       *
21  *************************************************************************/
22
23 #ifndef _ODE_ODEMATH_H_
24 #define _ODE_ODEMATH_H_
25
26 #include <ode/common.h>
27
28 #ifdef __cplusplus
29 extern "C" {
30 #endif
31
32
33 /* 3-way dot product. dDOTpq means that elements of `a' and `b' are spaced
34  * p and q indexes apart respectively. dDOT() means dDOT11.
35  */
36
37 #ifdef __cplusplus
38 inline dReal dDOT (const dReal *a, const dReal *b)
39   { return ((a)[0]*(b)[0] + (a)[1]*(b)[1] + (a)[2]*(b)[2]); }
40 inline dReal dDOT14(const dReal *a, const dReal *b)
41   { return ((a)[0]*(b)[0] + (a)[1]*(b)[4] + (a)[2]*(b)[8]); }
42 inline dReal dDOT41(const dReal *a, const dReal *b)
43   { return ((a)[0]*(b)[0] + (a)[4]*(b)[1] + (a)[8]*(b)[2]); }
44 inline dReal dDOT44(const dReal *a, const dReal *b)
45   { return ((a)[0]*(b)[0] + (a)[4]*(b)[4] + (a)[8]*(b)[8]); }
46 #else
47 #define dDOT(a,b)   ((a)[0]*(b)[0] + (a)[1]*(b)[1] + (a)[2]*(b)[2])
48 #define dDOT14(a,b) ((a)[0]*(b)[0] + (a)[1]*(b)[4] + (a)[2]*(b)[8])
49 #define dDOT41(a,b) ((a)[0]*(b)[0] + (a)[4]*(b)[1] + (a)[8]*(b)[2])
50 #define dDOT44(a,b) ((a)[0]*(b)[0] + (a)[4]*(b)[4] + (a)[8]*(b)[8])
51 #endif
52
53
54 /* cross product, set a = b x c. dCROSSpqr means that elements of `a', `b'
55  * and `c' are spaced p, q and r indexes apart respectively.
56  * dCROSS() means dCROSS111. `op' is normally `=', but you can set it to
57  * +=, -= etc to get other effects.
58  */
59
60 #define dCROSS(a,op,b,c) \
61   (a)[0] op ((b)[1]*(c)[2] - (b)[2]*(c)[1]); \
62   (a)[1] op ((b)[2]*(c)[0] - (b)[0]*(c)[2]); \
63   (a)[2] op ((b)[0]*(c)[1] - (b)[1]*(c)[0]);
64 #define dCROSSpqr(a,op,b,c,p,q,r) \
65   (a)[  0] op ((b)[  q]*(c)[2*r] - (b)[2*q]*(c)[  r]); \
66   (a)[  p] op ((b)[2*q]*(c)[  0] - (b)[  0]*(c)[2*r]); \
67   (a)[2*p] op ((b)[  0]*(c)[  r] - (b)[  q]*(c)[  0]);
68 #define dCROSS114(a,op,b,c) dCROSSpqr(a,op,b,c,1,1,4)
69 #define dCROSS141(a,op,b,c) dCROSSpqr(a,op,b,c,1,4,1)
70 #define dCROSS144(a,op,b,c) dCROSSpqr(a,op,b,c,1,4,4)
71 #define dCROSS411(a,op,b,c) dCROSSpqr(a,op,b,c,4,1,1)
72 #define dCROSS414(a,op,b,c) dCROSSpqr(a,op,b,c,4,1,4)
73 #define dCROSS441(a,op,b,c) dCROSSpqr(a,op,b,c,4,4,1)
74 #define dCROSS444(a,op,b,c) dCROSSpqr(a,op,b,c,4,4,4)
75
76
77 /* set a 3x3 submatrix of A to a matrix such that submatrix(A)*b = a x b.
78  * A is stored by rows, and has `skip' elements per row. the matrix is
79  * assumed to be already zero, so this does not write zero elements!
80  * if (plus,minus) is (+,-) then a positive version will be written.
81  * if (plus,minus) is (-,+) then a negative version will be written.
82  */
83
84 #define dCROSSMAT(A,a,skip,plus,minus) \
85   (A)[1] = minus (a)[2]; \
86   (A)[2] = plus (a)[1]; \
87   (A)[(skip)+0] = plus (a)[2]; \
88   (A)[(skip)+2] = minus (a)[0]; \
89   (A)[2*(skip)+0] = minus (a)[1]; \
90   (A)[2*(skip)+1] = plus (a)[0];
91
92
93 /* compute the distance between two 3-vectors (oops, C++!) */
94 #ifdef __cplusplus
95 inline dReal dDISTANCE (const dVector3 a, const dVector3 b)
96   { return dSqrt( (a[0]-b[0])*(a[0]-b[0]) + (a[1]-b[1])*(a[1]-b[1]) +
97                    (a[2]-b[2])*(a[2]-b[2]) ); }
98 #else
99 #define dDISTANCE(a,b) \
100  (dSqrt( ((a)[0]-(b)[0])*((a)[0]-(b)[0]) + ((a)[1]-(b)[1])*((a)[1]-(b)[1]) + \
101          ((a)[2]-(b)[2])*((a)[2]-(b)[2]) ))
102 #endif
103
104
105 /* normalize 3x1 and 4x1 vectors (i.e. scale them to unit length) */
106 void dNormalize3 (dVector3 a);
107 void dNormalize4 (dVector4 a);
108
109
110 /* given a unit length "normal" vector n, generate vectors p and q vectors
111  * that are an orthonormal basis for the plane space perpendicular to n.
112  * i.e. this makes p,q such that n,p,q are all perpendicular to each other.
113  * q will equal n x p. if n is not unit length then p will be unit length but
114  * q wont be.
115  */
116
117 void dPlaneSpace (const dVector3 n, dVector3 p, dVector3 q);
118
119
120 /* special case matrix multipication, with operator selection */
121
122 #define dMULTIPLYOP0_331(A,op,B,C) \
123   (A)[0] op dDOT((B),(C)); \
124   (A)[1] op dDOT((B+4),(C)); \
125   (A)[2] op dDOT((B+8),(C));
126 #define dMULTIPLYOP1_331(A,op,B,C) \
127   (A)[0] op dDOT41((B),(C)); \
128   (A)[1] op dDOT41((B+1),(C)); \
129   (A)[2] op dDOT41((B+2),(C));
130 #define dMULTIPLYOP0_133(A,op,B,C) \
131   (A)[0] op dDOT14((B),(C)); \
132   (A)[1] op dDOT14((B),(C+1)); \
133   (A)[2] op dDOT14((B),(C+2));
134 #define dMULTIPLYOP0_333(A,op,B,C) \
135   (A)[0] op dDOT14((B),(C)); \
136   (A)[1] op dDOT14((B),(C+1)); \
137   (A)[2] op dDOT14((B),(C+2)); \
138   (A)[4] op dDOT14((B+4),(C)); \
139   (A)[5] op dDOT14((B+4),(C+1)); \
140   (A)[6] op dDOT14((B+4),(C+2)); \
141   (A)[8] op dDOT14((B+8),(C)); \
142   (A)[9] op dDOT14((B+8),(C+1)); \
143   (A)[10] op dDOT14((B+8),(C+2));
144 #define dMULTIPLYOP1_333(A,op,B,C) \
145   (A)[0] op dDOT44((B),(C)); \
146   (A)[1] op dDOT44((B),(C+1)); \
147   (A)[2] op dDOT44((B),(C+2)); \
148   (A)[4] op dDOT44((B+1),(C)); \
149   (A)[5] op dDOT44((B+1),(C+1)); \
150   (A)[6] op dDOT44((B+1),(C+2)); \
151   (A)[8] op dDOT44((B+2),(C)); \
152   (A)[9] op dDOT44((B+2),(C+1)); \
153   (A)[10] op dDOT44((B+2),(C+2));
154 #define dMULTIPLYOP2_333(A,op,B,C) \
155   (A)[0] op dDOT((B),(C)); \
156   (A)[1] op dDOT((B),(C+4)); \
157   (A)[2] op dDOT((B),(C+8)); \
158   (A)[4] op dDOT((B+4),(C)); \
159   (A)[5] op dDOT((B+4),(C+4)); \
160   (A)[6] op dDOT((B+4),(C+8)); \
161   (A)[8] op dDOT((B+8),(C)); \
162   (A)[9] op dDOT((B+8),(C+4)); \
163   (A)[10] op dDOT((B+8),(C+8));
164
165 #ifdef __cplusplus
166
167 inline void dMULTIPLY0_331(dReal *A, const dReal *B, const dReal *C)
168   { dMULTIPLYOP0_331(A,=,B,C) }
169 inline void dMULTIPLY1_331(dReal *A, const dReal *B, const dReal *C)
170   { dMULTIPLYOP1_331(A,=,B,C) }
171 inline void dMULTIPLY0_133(dReal *A, const dReal *B, const dReal *C)
172   { dMULTIPLYOP0_133(A,=,B,C) }
173 inline void dMULTIPLY0_333(dReal *A, const dReal *B, const dReal *C)
174   { dMULTIPLYOP0_333(A,=,B,C) }
175 inline void dMULTIPLY1_333(dReal *A, const dReal *B, const dReal *C)
176   { dMULTIPLYOP1_333(A,=,B,C) }
177 inline void dMULTIPLY2_333(dReal *A, const dReal *B, const dReal *C)
178   { dMULTIPLYOP2_333(A,=,B,C) }
179
180 inline void dMULTIPLYADD0_331(dReal *A, const dReal *B, const dReal *C)
181   { dMULTIPLYOP0_331(A,+=,B,C) }
182 inline void dMULTIPLYADD1_331(dReal *A, const dReal *B, const dReal *C)
183   { dMULTIPLYOP1_331(A,+=,B,C) }
184 inline void dMULTIPLYADD0_133(dReal *A, const dReal *B, const dReal *C)
185   { dMULTIPLYOP0_133(A,+=,B,C) }
186 inline void dMULTIPLYADD0_333(dReal *A, const dReal *B, const dReal *C)
187   { dMULTIPLYOP0_333(A,+=,B,C) }
188 inline void dMULTIPLYADD1_333(dReal *A, const dReal *B, const dReal *C)
189   { dMULTIPLYOP1_333(A,+=,B,C) }
190 inline void dMULTIPLYADD2_333(dReal *A, const dReal *B, const dReal *C)
191   { dMULTIPLYOP2_333(A,+=,B,C) }
192
193 #else
194
195 #define dMULTIPLY0_331(A,B,C) dMULTIPLYOP0_331(A,=,B,C)
196 #define dMULTIPLY1_331(A,B,C) dMULTIPLYOP1_331(A,=,B,C)
197 #define dMULTIPLY0_133(A,B,C) dMULTIPLYOP0_133(A,=,B,C)
198 #define dMULTIPLY0_333(A,B,C) dMULTIPLYOP0_333(A,=,B,C)
199 #define dMULTIPLY1_333(A,B,C) dMULTIPLYOP1_333(A,=,B,C)
200 #define dMULTIPLY2_333(A,B,C) dMULTIPLYOP2_333(A,=,B,C)
201
202 #define dMULTIPLYADD0_331(A,B,C) dMULTIPLYOP0_331(A,+=,B,C)
203 #define dMULTIPLYADD1_331(A,B,C) dMULTIPLYOP1_331(A,+=,B,C)
204 #define dMULTIPLYADD0_133(A,B,C) dMULTIPLYOP0_133(A,+=,B,C)
205 #define dMULTIPLYADD0_333(A,B,C) dMULTIPLYOP0_333(A,+=,B,C)
206 #define dMULTIPLYADD1_333(A,B,C) dMULTIPLYOP1_333(A,+=,B,C)
207 #define dMULTIPLYADD2_333(A,B,C) dMULTIPLYOP2_333(A,+=,B,C)
208
209 #endif
210
211
212 #ifdef __cplusplus
213 }
214 #endif
215
216 #endif