Back to home page

EIC code displayed by LXR

 
 

    


Warning, /include/Geant4/tools/lina/mat4 is written in an unsupported language. File is not indexed.

0001 // Copyright (C) 2010, Guy Barrand. All rights reserved.
0002 // See the file tools.license for terms.
0003 
0004 #ifndef tools_mat4
0005 #define tools_mat4
0006 
0007 #include "mat"
0008 
0009 //#include <cmath>
0010 
0011 namespace tools {
0012 
0013 template <class T>
0014 class mat4 : public mat<T,4> {
0015   typedef mat<T,4> parent;
0016   typedef mat<T,4> pr;
0017 public:
0018 #ifdef TOOLS_MEM
0019   mat4(bool a_inc = true):parent(a_inc) {}
0020 #else
0021   mat4():parent() {}
0022 #endif
0023   mat4(const mat<T,4>& a_from):parent(a_from){}
0024   virtual ~mat4() {}
0025 public:
0026   mat4(const mat4& a_from):parent(a_from){}
0027   mat4& operator=(const mat4& a_from){
0028     parent::operator=(a_from);
0029     return *this;
0030   }
0031 public:
0032   mat4(const T& a_00,const T& a_01,const T& a_02,const T& a_03, //first  row
0033        const T& a_10,const T& a_11,const T& a_12,const T& a_13, //second row
0034        const T& a_20,const T& a_21,const T& a_22,const T& a_23, //third  row
0035        const T& a_30,const T& a_31,const T& a_32,const T& a_33) //fourth row
0036   {
0037     set_matrix(a_00,a_01,a_02,a_03,
0038                a_10,a_11,a_12,a_13,
0039                a_20,a_21,a_22,a_23,
0040                a_30,a_31,a_32,a_33);
0041   }
0042 public:
0043   void set_matrix(const mat4<T>& a_m) {parent::set_matrix(a_m);}
0044   void set_matrix(
0045     const T& a_00,const T& a_01,const T& a_02,const T& a_03, //1 row
0046     const T& a_10,const T& a_11,const T& a_12,const T& a_13, //2 row
0047     const T& a_20,const T& a_21,const T& a_22,const T& a_23, //3 row
0048     const T& a_30,const T& a_31,const T& a_32,const T& a_33) //4 row
0049   {
0050     //a_<R><C>
0051     //pr::m_vec[R + C * 4];
0052     pr::m_vec[0] = a_00;pr::m_vec[4] = a_01;pr::m_vec[ 8] = a_02;pr::m_vec[12] = a_03;
0053     pr::m_vec[1] = a_10;pr::m_vec[5] = a_11;pr::m_vec[ 9] = a_12;pr::m_vec[13] = a_13;
0054     pr::m_vec[2] = a_20;pr::m_vec[6] = a_21;pr::m_vec[10] = a_22;pr::m_vec[14] = a_23;
0055     pr::m_vec[3] = a_30;pr::m_vec[7] = a_31;pr::m_vec[11] = a_32;pr::m_vec[15] = a_33;
0056   }
0057 
0058   void set_scale(const T& a_s) {_set_scale(a_s,a_s,a_s,pr::m_vec);}
0059   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);}
0060   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);}
0061   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)) {
0062     _set_rotate(a_x,a_y,a_z,a_angle,pr::m_vec,a_sin,a_cos);
0063   }
0064   void set_ortho(const T& a_l,const T& a_r,   //left,right
0065                  const T& a_b,const T& a_t,   //bottom,top
0066                  const T& a_n,const T& a_f) { //znear,zfar
0067     // from man glOrtho.
0068     T tx = -(a_r+a_l)/(a_r-a_l);
0069     T ty = -(a_t+a_b)/(a_t-a_b);
0070     T tz = -(a_f+a_n)/(a_f-a_n);
0071 
0072     T a_00,a_01,a_02,a_03;
0073     T a_10,a_11,a_12,a_13;
0074     T a_20,a_21,a_22,a_23;
0075     T a_30,a_31,a_32,a_33;
0076     a_00 = 2/(a_r-a_l);a_01 =           0;a_02 =            0;a_03 = tx;
0077     a_10 =           0;a_11 = 2/(a_t-a_b);a_12 =            0;a_13 = ty;
0078     a_20 =           0;a_21 =           0;a_22 = -2/(a_f-a_n);a_23 = tz;
0079     a_30 =           0;a_31 =           0;a_32 =            0;a_33 =  1;
0080     set_matrix(a_00,a_01,a_02,a_03,  //1 row
0081                a_10,a_11,a_12,a_13,  //2 row
0082                a_20,a_21,a_22,a_23,  //3 row
0083                a_30,a_31,a_32,a_33); //4 row
0084 
0085     //NOTE : Z(x,y, z,1) = -2z/(f-n)+tz = [-2z-f-n]/(f-n)
0086     //       W(x,y, z,1) =  1 -> Z/W = Z
0087     //       Z(x,y,-n,1) = -1
0088     //       Z(x,y,-f,1) =  1
0089 
0090     //       Z(x,y,-(f+n)/2,1) = 0
0091 
0092     //       X(x,0, z,1) = 2x/(r-l)+tx = [2x-r-l]/(r-l)
0093     //       X(r,0, z,1) = 1
0094 
0095     // the view direction is then (0,0,1) in the final projection.
0096   }
0097   void set_frustum(const T& a_l,const T& a_r,   //left,right
0098                    const T& a_b,const T& a_t,   //bottom,top
0099                    const T& a_n,const T& a_f) { //znear,zfar
0100     // from man glFrustum.
0101     T A = (a_r+a_l)/(a_r-a_l);
0102     T B = (a_t+a_b)/(a_t-a_b);
0103     T C = -(a_f+a_n)/(a_f-a_n);
0104     T D = -(2*a_f*a_n)/(a_f-a_n);
0105 
0106     T a_00,a_01,a_02,a_03;
0107     T a_10,a_11,a_12,a_13;
0108     T a_20,a_21,a_22,a_23;
0109     T a_30,a_31,a_32,a_33;
0110     a_00 = (2*a_n)/(a_r-a_l);a_01 =                 0;a_02 =  A;a_03 = 0;
0111     a_10 =                 0;a_11 = (2*a_n)/(a_t-a_b);a_12 =  B;a_13 = 0;
0112     a_20 =                 0;a_21 =                 0;a_22 =  C;a_23 = D;
0113     a_30 =                 0;a_31 =                 0;a_32 = -1;a_33 = 0;
0114     set_matrix(a_00,a_01,a_02,a_03,  //1 row
0115                a_10,a_11,a_12,a_13,  //2 row
0116                a_20,a_21,a_22,a_23,  //3 row
0117                a_30,a_31,a_32,a_33); //4 row
0118 
0119     //NOTE : Z(x,y, z,1) = C*z+D = -[(f+n)*z+2*f*n]/(f-n)
0120 
0121     //       Z(x,y,-n,1) = -[-fn-nn+2fn]/(f-n) = -[fn-nn]/(f-n) = -n
0122     //       W(x,y,-n,1) =  n
0123     // ->  Z/W(x,y,-n,1) = -1
0124 
0125     //       Z(x,y,-2fn/(f+n),1) = 0
0126     // ->  Z/W                   = 0
0127 
0128     //       Z(x,y,-f,1) = -[-ff-fn+2fn]/(f-n) = -[fn-ff]/(f-n) = f
0129     //       W(x,y,-f,1) = f
0130     // ->  Z/W(x,y,-f,1) = 1
0131 
0132     //       X(x,0, z,1) = 2nx/(r-l)+z(r+l)/(r-l) = [2nx+zr+zl]/(r-l)
0133     //       X(r,0,-n,1) = [nr-nl]/(r-l) = n
0134     //       W(r,0,-n,1) = n
0135     // ->  X/W(r,0,-n,1) = 1
0136 
0137     //       X(l,0,-n,1) = (2nl-n(r+l))/(r-l) = -n
0138     //       W(l,0,-n,1) = n
0139     // ->  X/W(l,0,-n,1) = -1
0140 
0141     // lrbt corners are in plane z=-1 at xy=+/-1.
0142 
0143     // eye ?
0144     // eye before ? (0,0,z,1) -> (zA,zB,zC+D,-z)  /W -> (-A,-B,-(C+D/z),1)
0145     //  z=0 -> (0,0,D=-2fn(f-n),0)
0146 
0147   }
0148 
0149   void get_translate(T& a_x,T& a_y,T& a_z) const {
0150     a_x = pr::m_vec[12];
0151     a_y = pr::m_vec[13];
0152     a_z = pr::m_vec[14];
0153   }
0154 
0155   void no_translate() {
0156     pr::m_vec[12] = 0;
0157     pr::m_vec[13] = 0;
0158     pr::m_vec[14] = 0;
0159   }
0160 
0161 //void set_translate_only(const T& a_x,const T& a_y,const T& a_z) {
0162 //  pr::m_vec[12] = a_x;
0163 //  pr::m_vec[13] = a_y;
0164 //  pr::m_vec[14] = a_z;
0165 //}
0166 
0167 public:
0168   void mul_4(T& a_x,T& a_y,T& a_z,T& a_p) const {
0169     // a_[x,y,z,p] = this * a_[x,y,z,p]
0170     //pr::m_vec[R + C * 4];
0171     //pr::m_vec[0]= 00;pr::m_vec[4] = 01;pr::m_vec[ 8] = 02;pr::m_vec[12] = 03;
0172     //pr::m_vec[1]= 10;pr::m_vec[5] = 11;pr::m_vec[ 9] = 12;pr::m_vec[13] = 13;
0173     //pr::m_vec[2]= 20;pr::m_vec[6] = 21;pr::m_vec[10] = 22;pr::m_vec[14] = 23;
0174     //pr::m_vec[3]= 30;pr::m_vec[7] = 31;pr::m_vec[11] = 32;pr::m_vec[15] = 33;
0175     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;
0176     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;
0177     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;
0178     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;
0179     a_x = x;
0180     a_y = y;
0181     a_z = z;
0182     a_p = p;
0183   }
0184   void mul_4_opt(T& a_x,T& a_y,T& a_z,T& a_p,T a_tmp[4]) const {
0185     // a_[x,y,z,p] = this * a_[x,y,z,p]
0186     //pr::m_vec[R + C * 4];
0187     //pr::m_vec[0]= 00;pr::m_vec[4] = 01;pr::m_vec[ 8] = 02;pr::m_vec[12] = 03;
0188     //pr::m_vec[1]= 10;pr::m_vec[5] = 11;pr::m_vec[ 9] = 12;pr::m_vec[13] = 13;
0189     //pr::m_vec[2]= 20;pr::m_vec[6] = 21;pr::m_vec[10] = 22;pr::m_vec[14] = 23;
0190     //pr::m_vec[3]= 30;pr::m_vec[7] = 31;pr::m_vec[11] = 32;pr::m_vec[15] = 33;
0191     a_tmp[0] = pr::m_vec[0]*a_x+pr::m_vec[4]*a_y+pr::m_vec[ 8]*a_z+pr::m_vec[12]*a_p;
0192     a_tmp[1] = pr::m_vec[1]*a_x+pr::m_vec[5]*a_y+pr::m_vec[ 9]*a_z+pr::m_vec[13]*a_p;
0193     a_tmp[2] = pr::m_vec[2]*a_x+pr::m_vec[6]*a_y+pr::m_vec[10]*a_z+pr::m_vec[14]*a_p;
0194     a_tmp[3] = pr::m_vec[3]*a_x+pr::m_vec[7]*a_y+pr::m_vec[11]*a_z+pr::m_vec[15]*a_p;
0195     a_x = a_tmp[0];
0196     a_y = a_tmp[1];
0197     a_z = a_tmp[2];
0198     a_p = a_tmp[3];
0199   }
0200   void mul_3(T& a_x,T& a_y,T& a_z) const {
0201     // a_[x,y,z] = this * a_[x,y,z]
0202     //pr::m_vec[R + C * 4];
0203     //pr::m_vec[0]= 00;pr::m_vec[4] = 01;pr::m_vec[ 8] = 02;pr::m_vec[12] = 03;
0204     //pr::m_vec[1]= 10;pr::m_vec[5] = 11;pr::m_vec[ 9] = 12;pr::m_vec[13] = 13;
0205     //pr::m_vec[2]= 20;pr::m_vec[6] = 21;pr::m_vec[10] = 22;pr::m_vec[14] = 23;
0206     //pr::m_vec[3]= 30;pr::m_vec[7] = 31;pr::m_vec[11] = 32;pr::m_vec[15] = 33;
0207     T x = pr::m_vec[0]*a_x+pr::m_vec[4]*a_y+pr::m_vec[ 8]*a_z+pr::m_vec[12];
0208     T y = pr::m_vec[1]*a_x+pr::m_vec[5]*a_y+pr::m_vec[ 9]*a_z+pr::m_vec[13];
0209     T z = pr::m_vec[2]*a_x+pr::m_vec[6]*a_y+pr::m_vec[10]*a_z+pr::m_vec[14];
0210     a_x = x;
0211     a_y = y;
0212     a_z = z;
0213   }
0214   void mul_3_opt(T& a_x,T& a_y,T& a_z,T a_tmp[3]) const {
0215     // a_[x,y,z] = this * a_[x,y,z]
0216     //pr::m_vec[R + C * 4];
0217     //pr::m_vec[0]= 00;pr::m_vec[4] = 01;pr::m_vec[ 8] = 02;pr::m_vec[12] = 03;
0218     //pr::m_vec[1]= 10;pr::m_vec[5] = 11;pr::m_vec[ 9] = 12;pr::m_vec[13] = 13;
0219     //pr::m_vec[2]= 20;pr::m_vec[6] = 21;pr::m_vec[10] = 22;pr::m_vec[14] = 23;
0220     //pr::m_vec[3]= 30;pr::m_vec[7] = 31;pr::m_vec[11] = 32;pr::m_vec[15] = 33;
0221     a_tmp[0] = pr::m_vec[0]*a_x+pr::m_vec[4]*a_y+pr::m_vec[ 8]*a_z+pr::m_vec[12];
0222     a_tmp[1] = pr::m_vec[1]*a_x+pr::m_vec[5]*a_y+pr::m_vec[ 9]*a_z+pr::m_vec[13];
0223     a_tmp[2] = pr::m_vec[2]*a_x+pr::m_vec[6]*a_y+pr::m_vec[10]*a_z+pr::m_vec[14];
0224     a_x = a_tmp[0];
0225     a_y = a_tmp[1];
0226     a_z = a_tmp[2];
0227   }
0228   void mul_2(T& a_x,T& a_y) const {
0229     // a_[x,y] = this * a_[x,y]
0230     //pr::m_vec[R + C * 4];
0231     //pr::m_vec[0]= 00;pr::m_vec[4] = 01;pr::m_vec[ 8] = 02;pr::m_vec[12] = 03;
0232     //pr::m_vec[1]= 10;pr::m_vec[5] = 11;pr::m_vec[ 9] = 12;pr::m_vec[13] = 13;
0233     //pr::m_vec[2]= 20;pr::m_vec[6] = 21;pr::m_vec[10] = 22;pr::m_vec[14] = 23;
0234     //pr::m_vec[3]= 30;pr::m_vec[7] = 31;pr::m_vec[11] = 32;pr::m_vec[15] = 33;
0235     T x = pr::m_vec[0]*a_x+pr::m_vec[4]*a_y+pr::m_vec[12];
0236     T y = pr::m_vec[1]*a_x+pr::m_vec[5]*a_y+pr::m_vec[13];
0237     a_x = x;
0238     a_y = y;
0239   }
0240 
0241   void mul_dir_3(T& a_x,T& a_y,T& a_z) const {
0242     // used to multiply normals.
0243     // a_[x,y,z] = rot_part(this) * a_[x,y,z]
0244 
0245     T x = pr::m_vec[0]*a_x+pr::m_vec[4]*a_y+pr::m_vec[ 8]*a_z;
0246     T y = pr::m_vec[1]*a_x+pr::m_vec[5]*a_y+pr::m_vec[ 9]*a_z;
0247     T z = pr::m_vec[2]*a_x+pr::m_vec[6]*a_y+pr::m_vec[10]*a_z;
0248     a_x = x;
0249     a_y = y;
0250     a_z = z;
0251   }
0252   void mul_dir_3_opt(T& a_x,T& a_y,T& a_z,T a_tmp[3]) const {
0253     // used to multiply normals.
0254     // a_[x,y,z] = rot_part(this) * a_[x,y,z]
0255 
0256     a_tmp[0] = pr::m_vec[0]*a_x+pr::m_vec[4]*a_y+pr::m_vec[ 8]*a_z;
0257     a_tmp[1] = pr::m_vec[1]*a_x+pr::m_vec[5]*a_y+pr::m_vec[ 9]*a_z;
0258     a_tmp[2] = pr::m_vec[2]*a_x+pr::m_vec[6]*a_y+pr::m_vec[10]*a_z;
0259     a_x = a_tmp[0];
0260     a_y = a_tmp[1];
0261     a_z = a_tmp[2];
0262   }
0263 
0264   void mul_trans_3(T& a_x,T& a_y,T& a_z) const {
0265     // used in sg::healpix.
0266     // a_[x,y,z] = trans_part(this) * a_[x,y,z]
0267     a_x += pr::m_vec[12];
0268     a_y += pr::m_vec[13];
0269     a_z += pr::m_vec[14];
0270   }
0271 
0272   template <class VEC>
0273   void mul_dir_3(VEC& a_dir) const {mul_dir_3(a_dir[0],a_dir[1],a_dir[2]);}
0274 
0275   void mul_scale(const T& a_sx,const T& a_sy,const T& a_sz) {
0276     // this = this * mat4_scale(a_s[x,y,z]
0277     //pr::m_vec[R + C * 4];
0278     //pr::m_vec[0]= 00;pr::m_vec[4] = 01;pr::m_vec[ 8] = 02;pr::m_vec[12] = 03;
0279     //pr::m_vec[1]= 10;pr::m_vec[5] = 11;pr::m_vec[ 9] = 12;pr::m_vec[13] = 13;
0280     //pr::m_vec[2]= 20;pr::m_vec[6] = 21;pr::m_vec[10] = 22;pr::m_vec[14] = 23;
0281     //pr::m_vec[3]= 30;pr::m_vec[7] = 31;pr::m_vec[11] = 32;pr::m_vec[15] = 33;
0282     pr::m_vec[0] *= a_sx;
0283     pr::m_vec[1] *= a_sx;
0284     pr::m_vec[2] *= a_sx;
0285     pr::m_vec[3] *= a_sx;
0286 
0287     pr::m_vec[4] *= a_sy;
0288     pr::m_vec[5] *= a_sy;
0289     pr::m_vec[6] *= a_sy;
0290     pr::m_vec[7] *= a_sy;
0291 
0292     pr::m_vec[ 8] *= a_sz;
0293     pr::m_vec[ 9] *= a_sz;
0294     pr::m_vec[10] *= a_sz;
0295     pr::m_vec[11] *= a_sz;
0296   }
0297   void mul_scale(const T& a_s) {
0298     pr::m_vec[0] *= a_s;
0299     pr::m_vec[1] *= a_s;
0300     pr::m_vec[2] *= a_s;
0301     pr::m_vec[3] *= a_s;
0302 
0303     pr::m_vec[4] *= a_s;
0304     pr::m_vec[5] *= a_s;
0305     pr::m_vec[6] *= a_s;
0306     pr::m_vec[7] *= a_s;
0307 
0308     pr::m_vec[ 8] *= a_s;
0309     pr::m_vec[ 9] *= a_s;
0310     pr::m_vec[10] *= a_s;
0311     pr::m_vec[11] *= a_s;
0312   }
0313   void mul_translate(const T& a_x,const T& a_y,const T& a_z) {
0314     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];
0315     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];
0316     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];
0317     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];
0318   }
0319 
0320   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)) {
0321     T rot[16];
0322     _set_rotate(a_x,a_y,a_z,a_angle,rot,a_sin,a_cos);
0323     parent::_mul_mtx(rot);
0324   }
0325 
0326   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)) {
0327     T _m[16];
0328     _set_rotate(a_x,a_y,a_z,a_angle,_m,a_sin,a_cos);
0329     parent::_left_mul_mtx(_m);
0330   }
0331 
0332   void left_mul_scale(const T& a_x,const T& a_y,const T& a_z) {
0333     T _m[16];
0334     _set_scale(a_x,a_y,a_z,_m);
0335     parent::_left_mul_mtx(_m);
0336   }
0337 
0338   void left_mul_translate(const T& a_x,const T& a_y,const T& a_z) {
0339     T _m[16];
0340     _set_translate(a_x,a_y,a_z,_m);
0341     parent::_left_mul_mtx(_m);
0342   }
0343 
0344   void v00(const T& a_value){pr::m_vec[0+0*4] = a_value;}
0345   void v10(const T& a_value){pr::m_vec[1+0*4] = a_value;}
0346   void v20(const T& a_value){pr::m_vec[2+0*4] = a_value;}
0347   void v30(const T& a_value){pr::m_vec[3+0*4] = a_value;}
0348 
0349   void v01(const T& a_value){pr::m_vec[0+1*4] = a_value;}
0350   void v11(const T& a_value){pr::m_vec[1+1*4] = a_value;}
0351   void v21(const T& a_value){pr::m_vec[2+1*4] = a_value;}
0352   void v31(const T& a_value){pr::m_vec[3+1*4] = a_value;}
0353 
0354   void v02(const T& a_value){pr::m_vec[0+2*4] = a_value;}
0355   void v12(const T& a_value){pr::m_vec[1+2*4] = a_value;}
0356   void v22(const T& a_value){pr::m_vec[2+2*4] = a_value;}
0357   void v32(const T& a_value){pr::m_vec[3+2*4] = a_value;}
0358 
0359   void v03(const T& a_value){pr::m_vec[0+3*4] = a_value;}
0360   void v13(const T& a_value){pr::m_vec[1+3*4] = a_value;}
0361   void v23(const T& a_value){pr::m_vec[2+3*4] = a_value;}
0362   void v33(const T& a_value){pr::m_vec[3+3*4] = a_value;}
0363 
0364   const T& v00() const {return pr::m_vec[0+0*4];}
0365   const T& v10() const {return pr::m_vec[1+0*4];}
0366   const T& v20() const {return pr::m_vec[2+0*4];}
0367   const T& v30() const {return pr::m_vec[3+0*4];}
0368 
0369   const T& v01() const {return pr::m_vec[0+1*4];}
0370   const T& v11() const {return pr::m_vec[1+1*4];}
0371   const T& v21() const {return pr::m_vec[2+1*4];}
0372   const T& v31() const {return pr::m_vec[3+1*4];}
0373 
0374   const T& v02() const {return pr::m_vec[0+2*4];}
0375   const T& v12() const {return pr::m_vec[1+2*4];}
0376   const T& v22() const {return pr::m_vec[2+2*4];}
0377   const T& v32() const {return pr::m_vec[3+2*4];}
0378 
0379   const T& v03() const {return pr::m_vec[0+3*4];}
0380   const T& v13() const {return pr::m_vec[1+3*4];}
0381   const T& v23() const {return pr::m_vec[2+3*4];}
0382   const T& v33() const {return pr::m_vec[3+3*4];}
0383 
0384 protected:
0385   static void _set_translate(const T& a_x,const T& a_y,const T& a_z,T v[]) {
0386     v[0] = T(1);v[4] =    0;v[ 8] =    0;v[12] = a_x;
0387     v[1] =    0;v[5] = T(1);v[ 9] =    0;v[13] = a_y;
0388     v[2] =    0;v[6] =    0;v[10] = T(1);v[14] = a_z;
0389     v[3] =    0;v[7] =    0;v[11] =    0;v[15] = T(1);
0390   }
0391 
0392   static void _set_scale(const T& a_1,const T& a_2,const T& a_3,T v[]) {
0393     v[0] = a_1;v[4] =   0;v[ 8] =   0;v[12] = 0;
0394     v[1] =   0;v[5] = a_2;v[ 9] =   0;v[13] = 0;
0395     v[2] =   0;v[6] =   0;v[10] = a_3;v[14] = 0;
0396     v[3] =   0;v[7] =   0;v[11] =   0;v[15] = T(1);
0397   }
0398 
0399   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)) {
0400     //WARNING : (a_x,a_y,a_z) must be a normalized vector.
0401     // 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).
0402     T si = a_sin(a_angle);
0403     T co = a_cos(a_angle);
0404     T x = a_x;
0405     T y = a_y;
0406     T z = a_z;
0407     T x2 = x*x;
0408     T y2 = y*y;
0409     T z2 = z*z;
0410     T xy = x*y;
0411     T xz = x*z;
0412     T yz = y*z;
0413     v[0] = x2+(1-x2)*co;v[4] = xy*(1-co)-z*si;v[ 8] = xz*(1-co)+y*si;v[12] = 0;
0414     v[1] = xy*(1-co)+z*si;v[5] = y2+(1-y2)*co;v[ 9] = yz*(1-co)-x*si;v[13] = 0;
0415     v[2] = xz*(1-co)-y*si;v[6] = yz*(1-co)+x*si;v[10] = z2+(1-z2)*co;v[14] = 0;
0416     v[3] =            0;v[7] =            0;v[11] =            0;v[15] = 1;
0417 
0418     // If :
0419     //     n =(x,y,z)
0420     //    n2 = x2+y2+z2 = 1
0421     //   n.E = x*E1+y*E2+z*E3
0422     // with :
0423     //   E1            E2             E3
0424     //    0  0  0       0  0 -1        0  1  0
0425     //    0  0  1       0  0  0       -1  0  0
0426     //    0 -1  0       1  0  0        0  0  0
0427     //
0428     // R(r,c) = cos(theta)*Id(r,c)+(1-cos(theta))*nr*nc-sin(theta)*(n.E)(r,c)
0429     //
0430     // R = exp(-theta*(n.E)) // active rotation (it models the rotation of an object by theta along n).
0431 
0432   }
0433 
0434 public:
0435   void mul_mtx_rot_root(const T& a_00,const T& a_01,const T& a_02, //1 row
0436                         const T& a_10,const T& a_11,const T& a_12, //2 row
0437                         const T& a_20,const T& a_21,const T& a_22) //3 row
0438   {
0439     T* tv = pr::m_vec;
0440     //pr::m_vec[0] = 00;pr::m_vec[4] = 01;pr::m_vec[ 8] = 02;pr::m_vec[12] = 03;
0441     //pr::m_vec[1] = 10;pr::m_vec[5] = 11;pr::m_vec[ 9] = 12;pr::m_vec[13] = 13;
0442     //pr::m_vec[2] = 20;pr::m_vec[6] = 21;pr::m_vec[10] = 22;pr::m_vec[14] = 23;
0443     //pr::m_vec[3] = 30;pr::m_vec[7] = 31;pr::m_vec[11] = 32;pr::m_vec[15] = 33;
0444 
0445     T tv_0  = tv[0];
0446     T tv_1  = tv[1];
0447     T tv_2  = tv[2];
0448     T tv_3  = tv[3];
0449     T tv_4  = tv[4];
0450     T tv_5  = tv[5];
0451     T tv_6  = tv[6];
0452     T tv_7  = tv[7];
0453     T tv_8  = tv[8];
0454     T tv_9  = tv[9];
0455     T tv_10 = tv[10];
0456     T tv_11 = tv[11];
0457     T tv_12 = tv[12];
0458     T tv_13 = tv[13];
0459     T tv_14 = tv[14];
0460     T tv_15 = tv[15];
0461 
0462     T fv_0  = a_00;
0463     T fv_1  = a_10;
0464     T fv_2  = a_20;
0465     //T fv_3  =    0;
0466 
0467     T fv_4  = a_01;
0468     T fv_5  = a_11;
0469     T fv_6  = a_21;
0470     //T fv_7  =    0;
0471 
0472     T fv_8  = a_02;
0473     T fv_9  = a_12;
0474     T fv_10 = a_22;
0475     //T fv_11 =    0;
0476 
0477     //T fv_12 = 0;
0478     //T fv_13 = 0;
0479     //T fv_14 = 0;
0480     //T fv_15 = 1;
0481 
0482     tv[0] = tv_0*fv_0+tv_4*fv_1+ tv_8*fv_2;
0483     tv[1] = tv_1*fv_0+tv_5*fv_1+ tv_9*fv_2;
0484     tv[2] = tv_2*fv_0+tv_6*fv_1+tv_10*fv_2;
0485     tv[3] = tv_3*fv_0+tv_7*fv_1+tv_11*fv_2;
0486 
0487     tv[4] = tv_0*fv_4+tv_4*fv_5+ tv_8*fv_6;
0488     tv[5] = tv_1*fv_4+tv_5*fv_5+ tv_9*fv_6;
0489     tv[6] = tv_2*fv_4+tv_6*fv_5+tv_10*fv_6;
0490     tv[7] = tv_3*fv_4+tv_7*fv_5+tv_11*fv_6;
0491 
0492     tv[8]  = tv_0*fv_8+tv_4*fv_9+ tv_8*fv_10;
0493     tv[9]  = tv_1*fv_8+tv_5*fv_9+ tv_9*fv_10;
0494     tv[10] = tv_2*fv_8+tv_6*fv_9+tv_10*fv_10;
0495     tv[11] = tv_3*fv_8+tv_7*fv_9+tv_11*fv_10;
0496 
0497     tv[12] = tv_12;
0498     tv[13] = tv_13;
0499     tv[14] = tv_14;
0500     tv[15] = tv_15;
0501   }
0502 
0503   template <class MAT3>
0504   void set_mat3(MAT3& a_m3) {
0505     a_m3.set_matrix(pr::m_vec[0],pr::m_vec[4],pr::m_vec[ 8],   //1 row
0506                     pr::m_vec[1],pr::m_vec[5],pr::m_vec[ 9],   //2 row
0507                     pr::m_vec[2],pr::m_vec[6],pr::m_vec[10]);  //3 row
0508   }
0509 private:static void check_instantiation() {mat4<float> dummy;}
0510 };
0511 
0512 ////////////////////////////////////////////////
0513 /// common matrices : //////////////////////////
0514 ////////////////////////////////////////////////
0515 
0516 #ifdef TOOLS_MEM
0517 template <class T> inline const mat4<T>& mat4_zero() {static const mat4<T> s_v(false);return s_v;}
0518 #else
0519 template <class T> inline const mat4<T>& mat4_zero() {static const mat4<T> s_v;return s_v;}
0520 #endif
0521 
0522 /*
0523 template <class T>
0524 inline const mat4<T>& mat4_identity() {
0525   static mat4<T> s_v(false);
0526   static bool s_first = true;
0527   if(s_first) {s_v.set_identity();s_first=false;}
0528   return s_v;
0529 }
0530 */
0531 
0532 //for sf, mf :
0533 template <class T>
0534 inline const T* get_data(const mat4<T>& a_v) {return a_v.data();}
0535 
0536 }
0537 
0538 #include <ostream>
0539 
0540 namespace tools {
0541 
0542 template <class T>
0543 inline std::ostream& operator<<(std::ostream& a_out,const mat4<T>& a_mtx){
0544   const T* v = a_mtx.data();
0545   a_out << v[0] << "," << v[4] << "," << v[ 8] << "," << v[12] << std::endl
0546         << v[1] << "," << v[5] << "," << v[ 9] << "," << v[13] << std::endl
0547         << v[2] << "," << v[6] << "," << v[10] << "," << v[14] << std::endl
0548         << v[3] << "," << v[7] << "," << v[11] << "," << v[15] << std::endl;
0549   return a_out;
0550 }
0551 
0552 }
0553 
0554 #endif