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 11.0.p4)


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