Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/g4tools/include/tools/lina/mat4

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /externals/g4tools/include/tools/lina/mat4 (Version 11.3.0) and /externals/g4tools/include/tools/lina/mat4 (Version 9.4.p3)


  1 // Copyright (C) 2010, Guy Barrand. All rights    
  2 // See the file tools.license for terms.          
  3                                                   
  4 #ifndef tools_mat4                                
  5 #define tools_mat4                                
  6                                                   
  7 #include "mat"                                    
  8                                                   
  9 //#include <cmath>                                
 10                                                   
 11 namespace tools {                                 
 12                                                   
 13 template <class T>                                
 14 class mat4 : public mat<T,4> {                    
 15   typedef mat<T,4> parent;                        
 16   typedef mat<T,4> pr;                            
 17 public:                                           
 18 #ifdef TOOLS_MEM                                  
 19   mat4(bool a_inc = true):parent(a_inc) {}        
 20 #else                                             
 21   mat4():parent() {}                              
 22 #endif                                            
 23   mat4(const mat<T,4>& a_from):parent(a_from){    
 24   virtual ~mat4() {}                              
 25 public:                                           
 26   mat4(const mat4& a_from):parent(a_from){}       
 27   mat4& operator=(const mat4& a_from){            
 28     parent::operator=(a_from);                    
 29     return *this;                                 
 30   }                                               
 31 public:                                           
 32   mat4(const T& a_00,const T& a_01,const T& a_    
 33        const T& a_10,const T& a_11,const T& a_    
 34        const T& a_20,const T& a_21,const T& a_    
 35        const T& a_30,const T& a_31,const T& a_    
 36   {                                               
 37     set_matrix(a_00,a_01,a_02,a_03,               
 38                a_10,a_11,a_12,a_13,               
 39                a_20,a_21,a_22,a_23,               
 40                a_30,a_31,a_32,a_33);              
 41   }                                               
 42 public:                                           
 43   void set_matrix(const mat4<T>& a_m) {parent:    
 44   void set_matrix(                                
 45     const T& a_00,const T& a_01,const T& a_02,    
 46     const T& a_10,const T& a_11,const T& a_12,    
 47     const T& a_20,const T& a_21,const T& a_22,    
 48     const T& a_30,const T& a_31,const T& a_32,    
 49   {                                               
 50     //a_<R><C>                                    
 51     //pr::m_vec[R + C * 4];                       
 52     pr::m_vec[0] = a_00;pr::m_vec[4] = a_01;pr    
 53     pr::m_vec[1] = a_10;pr::m_vec[5] = a_11;pr    
 54     pr::m_vec[2] = a_20;pr::m_vec[6] = a_21;pr    
 55     pr::m_vec[3] = a_30;pr::m_vec[7] = a_31;pr    
 56   }                                               
 57                                                   
 58   void set_scale(const T& a_s) {_set_scale(a_s    
 59   void set_scale(const T& a_1,const T& a_2,con    
 60   void set_translate(const T& a_x,const T& a_y    
 61   void set_rotate(const T& a_x,const T& a_y,co    
 62     _set_rotate(a_x,a_y,a_z,a_angle,pr::m_vec,    
 63   }                                               
 64   void set_ortho(const T& a_l,const T& a_r,       
 65                  const T& a_b,const T& a_t,       
 66                  const T& a_n,const T& a_f) {     
 67     // from man glOrtho.                          
 68     T tx = -(a_r+a_l)/(a_r-a_l);                  
 69     T ty = -(a_t+a_b)/(a_t-a_b);                  
 70     T tz = -(a_f+a_n)/(a_f-a_n);                  
 71                                                   
 72     T a_00,a_01,a_02,a_03;                        
 73     T a_10,a_11,a_12,a_13;                        
 74     T a_20,a_21,a_22,a_23;                        
 75     T a_30,a_31,a_32,a_33;                        
 76     a_00 = 2/(a_r-a_l);a_01 =           0;a_02    
 77     a_10 =           0;a_11 = 2/(a_t-a_b);a_12    
 78     a_20 =           0;a_21 =           0;a_22    
 79     a_30 =           0;a_31 =           0;a_32    
 80     set_matrix(a_00,a_01,a_02,a_03,  //1 row      
 81                a_10,a_11,a_12,a_13,  //2 row      
 82                a_20,a_21,a_22,a_23,  //3 row      
 83                a_30,a_31,a_32,a_33); //4 row      
 84                                                   
 85     //NOTE : Z(x,y, z,1) = -2z/(f-n)+tz = [-2z    
 86     //       W(x,y, z,1) =  1 -> Z/W = Z          
 87     //       Z(x,y,-n,1) = -1                     
 88     //       Z(x,y,-f,1) =  1                     
 89                                                   
 90     //       Z(x,y,-(f+n)/2,1) = 0                
 91                                                   
 92     //       X(x,0, z,1) = 2x/(r-l)+tx = [2x-r    
 93     //       X(r,0, z,1) = 1                      
 94                                                   
 95     // the view direction is then (0,0,1) in t    
 96   }                                               
 97   void set_frustum(const T& a_l,const T& a_r,     
 98                    const T& a_b,const T& a_t,     
 99                    const T& a_n,const T& a_f)     
100     // from man glFrustum.                        
101     T A = (a_r+a_l)/(a_r-a_l);                    
102     T B = (a_t+a_b)/(a_t-a_b);                    
103     T C = -(a_f+a_n)/(a_f-a_n);                   
104     T D = -(2*a_f*a_n)/(a_f-a_n);                 
105                                                   
106     T a_00,a_01,a_02,a_03;                        
107     T a_10,a_11,a_12,a_13;                        
108     T a_20,a_21,a_22,a_23;                        
109     T a_30,a_31,a_32,a_33;                        
110     a_00 = (2*a_n)/(a_r-a_l);a_01 =               
111     a_10 =                 0;a_11 = (2*a_n)/(a    
112     a_20 =                 0;a_21 =               
113     a_30 =                 0;a_31 =               
114     set_matrix(a_00,a_01,a_02,a_03,  //1 row      
115                a_10,a_11,a_12,a_13,  //2 row      
116                a_20,a_21,a_22,a_23,  //3 row      
117                a_30,a_31,a_32,a_33); //4 row      
118                                                   
119     //NOTE : Z(x,y, z,1) = C*z+D = -[(f+n)*z+2    
120                                                   
121     //       Z(x,y,-n,1) = -[-fn-nn+2fn]/(f-n)    
122     //       W(x,y,-n,1) =  n                     
123     // ->  Z/W(x,y,-n,1) = -1                     
124                                                   
125     //       Z(x,y,-2fn/(f+n),1) = 0              
126     // ->  Z/W                   = 0              
127                                                   
128     //       Z(x,y,-f,1) = -[-ff-fn+2fn]/(f-n)    
129     //       W(x,y,-f,1) = f                      
130     // ->  Z/W(x,y,-f,1) = 1                      
131                                                   
132     //       X(x,0, z,1) = 2nx/(r-l)+z(r+l)/(r    
133     //       X(r,0,-n,1) = [nr-nl]/(r-l) = n      
134     //       W(r,0,-n,1) = n                      
135     // ->  X/W(r,0,-n,1) = 1                      
136                                                   
137     //       X(l,0,-n,1) = (2nl-n(r+l))/(r-l)     
138     //       W(l,0,-n,1) = n                      
139     // ->  X/W(l,0,-n,1) = -1                     
140                                                   
141     // lrbt corners are in plane z=-1 at xy=+/    
142                                                   
143     // eye ?                                      
144     // eye before ? (0,0,z,1) -> (zA,zB,zC+D,-    
145     //  z=0 -> (0,0,D=-2fn(f-n),0)                
146                                                   
147   }                                               
148                                                   
149   void get_translate(T& a_x,T& a_y,T& a_z) con    
150     a_x = pr::m_vec[12];                          
151     a_y = pr::m_vec[13];                          
152     a_z = pr::m_vec[14];                          
153   }                                               
154                                                   
155   void no_translate() {                           
156     pr::m_vec[12] = 0;                            
157     pr::m_vec[13] = 0;                            
158     pr::m_vec[14] = 0;                            
159   }                                               
160                                                   
161 //void set_translate_only(const T& a_x,const T    
162 //  pr::m_vec[12] = a_x;                          
163 //  pr::m_vec[13] = a_y;                          
164 //  pr::m_vec[14] = a_z;                          
165 //}                                               
166                                                   
167 public:                                           
168   void mul_4(T& a_x,T& a_y,T& a_z,T& a_p) cons    
169     // a_[x,y,z,p] = this * a_[x,y,z,p]           
170     //pr::m_vec[R + C * 4];                       
171     //pr::m_vec[0]= 00;pr::m_vec[4] = 01;pr::m    
172     //pr::m_vec[1]= 10;pr::m_vec[5] = 11;pr::m    
173     //pr::m_vec[2]= 20;pr::m_vec[6] = 21;pr::m    
174     //pr::m_vec[3]= 30;pr::m_vec[7] = 31;pr::m    
175     T x = pr::m_vec[0]*a_x+pr::m_vec[4]*a_y+pr    
176     T y = pr::m_vec[1]*a_x+pr::m_vec[5]*a_y+pr    
177     T z = pr::m_vec[2]*a_x+pr::m_vec[6]*a_y+pr    
178     T p = pr::m_vec[3]*a_x+pr::m_vec[7]*a_y+pr    
179     a_x = x;                                      
180     a_y = y;                                      
181     a_z = z;                                      
182     a_p = p;                                      
183   }                                               
184   void mul_4_opt(T& a_x,T& a_y,T& a_z,T& a_p,T    
185     // a_[x,y,z,p] = this * a_[x,y,z,p]           
186     //pr::m_vec[R + C * 4];                       
187     //pr::m_vec[0]= 00;pr::m_vec[4] = 01;pr::m    
188     //pr::m_vec[1]= 10;pr::m_vec[5] = 11;pr::m    
189     //pr::m_vec[2]= 20;pr::m_vec[6] = 21;pr::m    
190     //pr::m_vec[3]= 30;pr::m_vec[7] = 31;pr::m    
191     a_tmp[0] = pr::m_vec[0]*a_x+pr::m_vec[4]*a    
192     a_tmp[1] = pr::m_vec[1]*a_x+pr::m_vec[5]*a    
193     a_tmp[2] = pr::m_vec[2]*a_x+pr::m_vec[6]*a    
194     a_tmp[3] = pr::m_vec[3]*a_x+pr::m_vec[7]*a    
195     a_x = a_tmp[0];                               
196     a_y = a_tmp[1];                               
197     a_z = a_tmp[2];                               
198     a_p = a_tmp[3];                               
199   }                                               
200   void mul_3(T& a_x,T& a_y,T& a_z) const {        
201     // a_[x,y,z] = this * a_[x,y,z]               
202     //pr::m_vec[R + C * 4];                       
203     //pr::m_vec[0]= 00;pr::m_vec[4] = 01;pr::m    
204     //pr::m_vec[1]= 10;pr::m_vec[5] = 11;pr::m    
205     //pr::m_vec[2]= 20;pr::m_vec[6] = 21;pr::m    
206     //pr::m_vec[3]= 30;pr::m_vec[7] = 31;pr::m    
207     T x = pr::m_vec[0]*a_x+pr::m_vec[4]*a_y+pr    
208     T y = pr::m_vec[1]*a_x+pr::m_vec[5]*a_y+pr    
209     T z = pr::m_vec[2]*a_x+pr::m_vec[6]*a_y+pr    
210     a_x = x;                                      
211     a_y = y;                                      
212     a_z = z;                                      
213   }                                               
214   void mul_3_opt(T& a_x,T& a_y,T& a_z,T a_tmp[    
215     // a_[x,y,z] = this * a_[x,y,z]               
216     //pr::m_vec[R + C * 4];                       
217     //pr::m_vec[0]= 00;pr::m_vec[4] = 01;pr::m    
218     //pr::m_vec[1]= 10;pr::m_vec[5] = 11;pr::m    
219     //pr::m_vec[2]= 20;pr::m_vec[6] = 21;pr::m    
220     //pr::m_vec[3]= 30;pr::m_vec[7] = 31;pr::m    
221     a_tmp[0] = pr::m_vec[0]*a_x+pr::m_vec[4]*a    
222     a_tmp[1] = pr::m_vec[1]*a_x+pr::m_vec[5]*a    
223     a_tmp[2] = pr::m_vec[2]*a_x+pr::m_vec[6]*a    
224     a_x = a_tmp[0];                               
225     a_y = a_tmp[1];                               
226     a_z = a_tmp[2];                               
227   }                                               
228   void mul_2(T& a_x,T& a_y) const {               
229     // a_[x,y] = this * a_[x,y]                   
230     //pr::m_vec[R + C * 4];                       
231     //pr::m_vec[0]= 00;pr::m_vec[4] = 01;pr::m    
232     //pr::m_vec[1]= 10;pr::m_vec[5] = 11;pr::m    
233     //pr::m_vec[2]= 20;pr::m_vec[6] = 21;pr::m    
234     //pr::m_vec[3]= 30;pr::m_vec[7] = 31;pr::m    
235     T x = pr::m_vec[0]*a_x+pr::m_vec[4]*a_y+pr    
236     T y = pr::m_vec[1]*a_x+pr::m_vec[5]*a_y+pr    
237     a_x = x;                                      
238     a_y = y;                                      
239   }                                               
240                                                   
241   void mul_dir_3(T& a_x,T& a_y,T& a_z) const {    
242     // used to multiply normals.                  
243     // a_[x,y,z] = rot_part(this) * a_[x,y,z]     
244                                                   
245     T x = pr::m_vec[0]*a_x+pr::m_vec[4]*a_y+pr    
246     T y = pr::m_vec[1]*a_x+pr::m_vec[5]*a_y+pr    
247     T z = pr::m_vec[2]*a_x+pr::m_vec[6]*a_y+pr    
248     a_x = x;                                      
249     a_y = y;                                      
250     a_z = z;                                      
251   }                                               
252   void mul_dir_3_opt(T& a_x,T& a_y,T& a_z,T a_    
253     // used to multiply normals.                  
254     // a_[x,y,z] = rot_part(this) * a_[x,y,z]     
255                                                   
256     a_tmp[0] = pr::m_vec[0]*a_x+pr::m_vec[4]*a    
257     a_tmp[1] = pr::m_vec[1]*a_x+pr::m_vec[5]*a    
258     a_tmp[2] = pr::m_vec[2]*a_x+pr::m_vec[6]*a    
259     a_x = a_tmp[0];                               
260     a_y = a_tmp[1];                               
261     a_z = a_tmp[2];                               
262   }                                               
263                                                   
264   void mul_trans_3(T& a_x,T& a_y,T& a_z) const    
265     // used in sg::healpix.                       
266     // a_[x,y,z] = trans_part(this) * a_[x,y,z    
267     a_x += pr::m_vec[12];                         
268     a_y += pr::m_vec[13];                         
269     a_z += pr::m_vec[14];                         
270   }                                               
271                                                   
272   template <class VEC>                            
273   void mul_dir_3(VEC& a_dir) const {mul_dir_3(    
274                                                   
275   void mul_scale(const T& a_sx,const T& a_sy,c    
276     // this = this * mat4_scale(a_s[x,y,z]        
277     //pr::m_vec[R + C * 4];                       
278     //pr::m_vec[0]= 00;pr::m_vec[4] = 01;pr::m    
279     //pr::m_vec[1]= 10;pr::m_vec[5] = 11;pr::m    
280     //pr::m_vec[2]= 20;pr::m_vec[6] = 21;pr::m    
281     //pr::m_vec[3]= 30;pr::m_vec[7] = 31;pr::m    
282     pr::m_vec[0] *= a_sx;                         
283     pr::m_vec[1] *= a_sx;                         
284     pr::m_vec[2] *= a_sx;                         
285     pr::m_vec[3] *= a_sx;                         
286                                                   
287     pr::m_vec[4] *= a_sy;                         
288     pr::m_vec[5] *= a_sy;                         
289     pr::m_vec[6] *= a_sy;                         
290     pr::m_vec[7] *= a_sy;                         
291                                                   
292     pr::m_vec[ 8] *= a_sz;                        
293     pr::m_vec[ 9] *= a_sz;                        
294     pr::m_vec[10] *= a_sz;                        
295     pr::m_vec[11] *= a_sz;                        
296   }                                               
297   void mul_scale(const T& a_s) {                  
298     pr::m_vec[0] *= a_s;                          
299     pr::m_vec[1] *= a_s;                          
300     pr::m_vec[2] *= a_s;                          
301     pr::m_vec[3] *= a_s;                          
302                                                   
303     pr::m_vec[4] *= a_s;                          
304     pr::m_vec[5] *= a_s;                          
305     pr::m_vec[6] *= a_s;                          
306     pr::m_vec[7] *= a_s;                          
307                                                   
308     pr::m_vec[ 8] *= a_s;                         
309     pr::m_vec[ 9] *= a_s;                         
310     pr::m_vec[10] *= a_s;                         
311     pr::m_vec[11] *= a_s;                         
312   }                                               
313   void mul_translate(const T& a_x,const T& a_y    
314     pr::m_vec[12] = pr::m_vec[0]*a_x+pr::m_vec    
315     pr::m_vec[13] = pr::m_vec[1]*a_x+pr::m_vec    
316     pr::m_vec[14] = pr::m_vec[2]*a_x+pr::m_vec    
317     pr::m_vec[15] = pr::m_vec[3]*a_x+pr::m_vec    
318   }                                               
319                                                   
320   void mul_rotate(const T& a_x,const T& a_y,co    
321     T rot[16];                                    
322     _set_rotate(a_x,a_y,a_z,a_angle,rot,a_sin,    
323     parent::_mul_mtx(rot);                        
324   }                                               
325                                                   
326   void left_mul_rotate(const T& a_x,const T& a    
327     T _m[16];                                     
328     _set_rotate(a_x,a_y,a_z,a_angle,_m,a_sin,a    
329     parent::_left_mul_mtx(_m);                    
330   }                                               
331                                                   
332   void left_mul_scale(const T& a_x,const T& a_    
333     T _m[16];                                     
334     _set_scale(a_x,a_y,a_z,_m);                   
335     parent::_left_mul_mtx(_m);                    
336   }                                               
337                                                   
338   void left_mul_translate(const T& a_x,const T    
339     T _m[16];                                     
340     _set_translate(a_x,a_y,a_z,_m);               
341     parent::_left_mul_mtx(_m);                    
342   }                                               
343                                                   
344   void v00(const T& a_value){pr::m_vec[0+0*4]     
345   void v10(const T& a_value){pr::m_vec[1+0*4]     
346   void v20(const T& a_value){pr::m_vec[2+0*4]     
347   void v30(const T& a_value){pr::m_vec[3+0*4]     
348                                                   
349   void v01(const T& a_value){pr::m_vec[0+1*4]     
350   void v11(const T& a_value){pr::m_vec[1+1*4]     
351   void v21(const T& a_value){pr::m_vec[2+1*4]     
352   void v31(const T& a_value){pr::m_vec[3+1*4]     
353                                                   
354   void v02(const T& a_value){pr::m_vec[0+2*4]     
355   void v12(const T& a_value){pr::m_vec[1+2*4]     
356   void v22(const T& a_value){pr::m_vec[2+2*4]     
357   void v32(const T& a_value){pr::m_vec[3+2*4]     
358                                                   
359   void v03(const T& a_value){pr::m_vec[0+3*4]     
360   void v13(const T& a_value){pr::m_vec[1+3*4]     
361   void v23(const T& a_value){pr::m_vec[2+3*4]     
362   void v33(const T& a_value){pr::m_vec[3+3*4]     
363                                                   
364   const T& v00() const {return pr::m_vec[0+0*4    
365   const T& v10() const {return pr::m_vec[1+0*4    
366   const T& v20() const {return pr::m_vec[2+0*4    
367   const T& v30() const {return pr::m_vec[3+0*4    
368                                                   
369   const T& v01() const {return pr::m_vec[0+1*4    
370   const T& v11() const {return pr::m_vec[1+1*4    
371   const T& v21() const {return pr::m_vec[2+1*4    
372   const T& v31() const {return pr::m_vec[3+1*4    
373                                                   
374   const T& v02() const {return pr::m_vec[0+2*4    
375   const T& v12() const {return pr::m_vec[1+2*4    
376   const T& v22() const {return pr::m_vec[2+2*4    
377   const T& v32() const {return pr::m_vec[3+2*4    
378                                                   
379   const T& v03() const {return pr::m_vec[0+3*4    
380   const T& v13() const {return pr::m_vec[1+3*4    
381   const T& v23() const {return pr::m_vec[2+3*4    
382   const T& v33() const {return pr::m_vec[3+3*4    
383                                                   
384 protected:                                        
385   static void _set_translate(const T& a_x,cons    
386     v[0] = T(1);v[4] =    0;v[ 8] =    0;v[12]    
387     v[1] =    0;v[5] = T(1);v[ 9] =    0;v[13]    
388     v[2] =    0;v[6] =    0;v[10] = T(1);v[14]    
389     v[3] =    0;v[7] =    0;v[11] =    0;v[15]    
390   }                                               
391                                                   
392   static void _set_scale(const T& a_1,const T&    
393     v[0] = a_1;v[4] =   0;v[ 8] =   0;v[12] =     
394     v[1] =   0;v[5] = a_2;v[ 9] =   0;v[13] =     
395     v[2] =   0;v[6] =   0;v[10] = a_3;v[14] =     
396     v[3] =   0;v[7] =   0;v[11] =   0;v[15] =     
397   }                                               
398                                                   
399   static void _set_rotate(const T& a_x,const T    
400     //WARNING : (a_x,a_y,a_z) must be a normal    
401     // v[] = exp(-a_angle*n(a_x,a_y,a_z)*Rs).     
402     T si = a_sin(a_angle);                        
403     T co = a_cos(a_angle);                        
404     T x = a_x;                                    
405     T y = a_y;                                    
406     T z = a_z;                                    
407     T x2 = x*x;                                   
408     T y2 = y*y;                                   
409     T z2 = z*z;                                   
410     T xy = x*y;                                   
411     T xz = x*z;                                   
412     T yz = y*z;                                   
413     v[0] = x2+(1-x2)*co;v[4] = xy*(1-co)-z*si;    
414     v[1] = xy*(1-co)+z*si;v[5] = y2+(1-y2)*co;    
415     v[2] = xz*(1-co)-y*si;v[6] = yz*(1-co)+x*s    
416     v[3] =            0;v[7] =            0;v[    
417                                                   
418     // If :                                       
419     //     n =(x,y,z)                             
420     //    n2 = x2+y2+z2 = 1                       
421     //   n.E = x*E1+y*E2+z*E3                     
422     // with :                                     
423     //   E1            E2             E3          
424     //    0  0  0       0  0 -1        0  1  0    
425     //    0  0  1       0  0  0       -1  0  0    
426     //    0 -1  0       1  0  0        0  0  0    
427     //                                            
428     // R(r,c) = cos(theta)*Id(r,c)+(1-cos(thet    
429     //                                            
430     // R = exp(-theta*(n.E)) // active rotatio    
431                                                   
432   }                                               
433                                                   
434 public:                                           
435   void mul_mtx_rot_root(const T& a_00,const T&    
436                         const T& a_10,const T&    
437                         const T& a_20,const T&    
438   {                                               
439     T* tv = pr::m_vec;                            
440     //pr::m_vec[0] = 00;pr::m_vec[4] = 01;pr::    
441     //pr::m_vec[1] = 10;pr::m_vec[5] = 11;pr::    
442     //pr::m_vec[2] = 20;pr::m_vec[6] = 21;pr::    
443     //pr::m_vec[3] = 30;pr::m_vec[7] = 31;pr::    
444                                                   
445     T tv_0  = tv[0];                              
446     T tv_1  = tv[1];                              
447     T tv_2  = tv[2];                              
448     T tv_3  = tv[3];                              
449     T tv_4  = tv[4];                              
450     T tv_5  = tv[5];                              
451     T tv_6  = tv[6];                              
452     T tv_7  = tv[7];                              
453     T tv_8  = tv[8];                              
454     T tv_9  = tv[9];                              
455     T tv_10 = tv[10];                             
456     T tv_11 = tv[11];                             
457     T tv_12 = tv[12];                             
458     T tv_13 = tv[13];                             
459     T tv_14 = tv[14];                             
460     T tv_15 = tv[15];                             
461                                                   
462     T fv_0  = a_00;                               
463     T fv_1  = a_10;                               
464     T fv_2  = a_20;                               
465     //T fv_3  =    0;                             
466                                                   
467     T fv_4  = a_01;                               
468     T fv_5  = a_11;                               
469     T fv_6  = a_21;                               
470     //T fv_7  =    0;                             
471                                                   
472     T fv_8  = a_02;                               
473     T fv_9  = a_12;                               
474     T fv_10 = a_22;                               
475     //T fv_11 =    0;                             
476                                                   
477     //T fv_12 = 0;                                
478     //T fv_13 = 0;                                
479     //T fv_14 = 0;                                
480     //T fv_15 = 1;                                
481                                                   
482     tv[0] = tv_0*fv_0+tv_4*fv_1+ tv_8*fv_2;       
483     tv[1] = tv_1*fv_0+tv_5*fv_1+ tv_9*fv_2;       
484     tv[2] = tv_2*fv_0+tv_6*fv_1+tv_10*fv_2;       
485     tv[3] = tv_3*fv_0+tv_7*fv_1+tv_11*fv_2;       
486                                                   
487     tv[4] = tv_0*fv_4+tv_4*fv_5+ tv_8*fv_6;       
488     tv[5] = tv_1*fv_4+tv_5*fv_5+ tv_9*fv_6;       
489     tv[6] = tv_2*fv_4+tv_6*fv_5+tv_10*fv_6;       
490     tv[7] = tv_3*fv_4+tv_7*fv_5+tv_11*fv_6;       
491                                                   
492     tv[8]  = tv_0*fv_8+tv_4*fv_9+ tv_8*fv_10;     
493     tv[9]  = tv_1*fv_8+tv_5*fv_9+ tv_9*fv_10;     
494     tv[10] = tv_2*fv_8+tv_6*fv_9+tv_10*fv_10;     
495     tv[11] = tv_3*fv_8+tv_7*fv_9+tv_11*fv_10;     
496                                                   
497     tv[12] = tv_12;                               
498     tv[13] = tv_13;                               
499     tv[14] = tv_14;                               
500     tv[15] = tv_15;                               
501   }                                               
502                                                   
503   template <class MAT3>                           
504   void set_mat3(MAT3& a_m3) {                     
505     a_m3.set_matrix(pr::m_vec[0],pr::m_vec[4],    
506                     pr::m_vec[1],pr::m_vec[5],    
507                     pr::m_vec[2],pr::m_vec[6],    
508   }                                               
509 private:static void check_instantiation() {mat    
510 };                                                
511                                                   
512 //////////////////////////////////////////////    
513 /// common matrices : ////////////////////////    
514 //////////////////////////////////////////////    
515                                                   
516 #ifdef TOOLS_MEM                                  
517 template <class T> inline const mat4<T>& mat4_    
518 #else                                             
519 template <class T> inline const mat4<T>& mat4_    
520 #endif                                            
521                                                   
522 /*                                                
523 template <class T>                                
524 inline const mat4<T>& mat4_identity() {           
525   static mat4<T> s_v(false);                      
526   static bool s_first = true;                     
527   if(s_first) {s_v.set_identity();s_first=fals    
528   return s_v;                                     
529 }                                                 
530 */                                                
531                                                   
532 //for sf, mf :                                    
533 template <class T>                                
534 inline const T* get_data(const mat4<T>& a_v) {    
535                                                   
536 }                                                 
537                                                   
538 #include <ostream>                                
539                                                   
540 namespace tools {                                 
541                                                   
542 template <class T>                                
543 inline std::ostream& operator<<(std::ostream&     
544   const T* v = a_mtx.data();                      
545   a_out << v[0] << "," << v[4] << "," << v[ 8]    
546         << v[1] << "," << v[5] << "," << v[ 9]    
547         << v[2] << "," << v[6] << "," << v[10]    
548         << v[3] << "," << v[7] << "," << v[11]    
549   return a_out;                                   
550 }                                                 
551                                                   
552 }                                                 
553                                                   
554 #endif