Back to home page

EIC code displayed by LXR

 
 

    


Warning, /include/Geant4/tools/lina/vec3 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_vec3
0005 #define tools_vec3
0006 
0007 #include <cstddef> //size_t
0008 
0009 #ifdef TOOLS_MEM
0010 #include "../mem"
0011 #endif
0012 
0013 namespace tools {
0014 
0015 template <class T>
0016 class vec3 {
0017 #ifdef TOOLS_MEM
0018   static const std::string& s_class() {
0019     static const std::string s_v("tools::vec3");
0020     return s_v;
0021   }
0022 #endif
0023 public:
0024   typedef T elem_t;
0025   unsigned int dimension() const {return 3;}
0026 public:
0027   vec3(){
0028 #ifdef TOOLS_MEM
0029     mem::increment(s_class().c_str());
0030 #endif
0031     m_data[0] = T();
0032     m_data[1] = T();
0033     m_data[2] = T();
0034   }
0035   vec3(const T a_vec[3]) {
0036 #ifdef TOOLS_MEM
0037     mem::increment(s_class().c_str());
0038 #endif
0039     m_data[0] = a_vec[0];
0040     m_data[1] = a_vec[1];
0041     m_data[2] = a_vec[2];
0042   }
0043   vec3(const T& a0,const T& a1,const T& a2
0044 #ifdef TOOLS_MEM
0045        ,bool a_inc = true
0046 #endif
0047        ) {
0048 #ifdef TOOLS_MEM
0049     if(a_inc) mem::increment(s_class().c_str());
0050 #endif
0051     m_data[0] = a0;
0052     m_data[1] = a1;
0053     m_data[2] = a2;
0054   }
0055   virtual ~vec3() {
0056 #ifdef TOOLS_MEM
0057     mem::decrement(s_class().c_str());
0058 #endif
0059   }
0060 public:
0061   vec3(const vec3& a_from){
0062 #ifdef TOOLS_MEM
0063     mem::increment(s_class().c_str());
0064 #endif
0065     m_data[0] = a_from.m_data[0];
0066     m_data[1] = a_from.m_data[1];
0067     m_data[2] = a_from.m_data[2];
0068   }
0069   vec3& operator=(const vec3& a_from) {
0070     m_data[0] = a_from.m_data[0];
0071     m_data[1] = a_from.m_data[1];
0072     m_data[2] = a_from.m_data[2];
0073     return *this;
0074   }
0075 public:
0076   const T& v0() const { return m_data[0];}
0077   const T& v1() const { return m_data[1];}
0078   const T& v2() const { return m_data[2];}
0079 
0080   void v0(const T& a_value) { m_data[0] = a_value;}
0081   void v1(const T& a_value) { m_data[1] = a_value;}
0082   void v2(const T& a_value) { m_data[2] = a_value;}
0083 
0084   const T& x() const {return m_data[0];}
0085   const T& y() const {return m_data[1];}
0086   const T& z() const {return m_data[2];}
0087 
0088   T& x() {return m_data[0];}
0089   T& y() {return m_data[1];}
0090   T& z() {return m_data[2];}
0091 
0092   void set_value(const T& a0,const T& a1,const T& a2) {
0093     m_data[0] = a0;
0094     m_data[1] = a1;
0095     m_data[2] = a2;
0096   }
0097   void set_value(const T aV[3]) {
0098     m_data[0] = aV[0];
0099     m_data[1] = aV[1];
0100     m_data[2] = aV[2];
0101   }
0102   void value(T& a0,T& a1,T& a2) const {
0103     a0 = m_data[0];
0104     a1 = m_data[1];
0105     a2 = m_data[2];
0106   }
0107 
0108   //bool set_value(unsigned int a_index,const T& a_value) {
0109   //  if(a_index>=3) return false;
0110   //  m_[a_index] = a_value;
0111   //  return true;
0112   //}
0113 
0114   T length(T(*a_sqrt)(T)) const {
0115     return a_sqrt(m_data[0]*m_data[0]+m_data[1]*m_data[1]+m_data[2]*m_data[2]);
0116   }
0117 
0118   T normalize(T(*a_sqrt)(T)) {
0119     T norme = length(a_sqrt);
0120     if(norme==T()) return T();
0121     divide(norme);
0122     return norme;
0123   }
0124 
0125   T dot(const vec3& aV) const {
0126     return (m_data[0] * aV.m_data[0] +
0127             m_data[1] * aV.m_data[1] +
0128             m_data[2] * aV.m_data[2]);
0129   }
0130 
0131   void cross(const vec3<T>& aV,vec3<T>& a_value) const {
0132     a_value.set_value(m_data[1] * aV.m_data[2] - m_data[2] * aV.m_data[1],
0133                       m_data[2] * aV.m_data[0] - m_data[0] * aV.m_data[2],
0134                       m_data[0] * aV.m_data[1] - m_data[1] * aV.m_data[0]);
0135   }
0136 
0137   bool equal(const vec3& aV) const {
0138     if(m_data[0]!=aV.m_data[0]) return false;
0139     if(m_data[1]!=aV.m_data[1]) return false;
0140     if(m_data[2]!=aV.m_data[2]) return false;
0141     return true;
0142   }
0143 
0144   template <class PREC>
0145   bool equal_prec(const vec3& a_v,PREC a_prec,PREC(*a_fabs)(const T&)) const {
0146     if(&a_v==this) return true;
0147     for(unsigned int index=0;index<3;index++) {
0148       T diff = m_data[index]-a_v.m_data[index];
0149       if(a_fabs(diff)>=a_prec) return false;
0150     }
0151     return true;
0152   }
0153 
0154   vec3<T> _cross(const vec3<T>& aV) const {
0155     //not effective.
0156     return vec3<T>(m_data[1] * aV.m_data[2] - m_data[2] * aV.m_data[1],
0157                    m_data[2] * aV.m_data[0] - m_data[0] * aV.m_data[2],
0158                    m_data[0] * aV.m_data[1] - m_data[1] * aV.m_data[0]);
0159   }
0160 
0161   bool divide(const T& a_T) {
0162     if(a_T==T()) return false;
0163     m_data[0] /= a_T;
0164     m_data[1] /= a_T;
0165     m_data[2] /= a_T;
0166     return true;
0167   }
0168 
0169   void multiply(const T& a_T) {
0170     m_data[0] *= a_T;
0171     m_data[1] *= a_T;
0172     m_data[2] *= a_T;
0173   }
0174 
0175   void add(const vec3& a_v) {
0176     m_data[0] += a_v.m_data[0];
0177     m_data[1] += a_v.m_data[1];
0178     m_data[2] += a_v.m_data[2];
0179   }
0180 
0181   void add(const T& a0,const T& a1,const T& a2) {
0182     m_data[0] += a0;
0183     m_data[1] += a1;
0184     m_data[2] += a2;
0185   }
0186 
0187   void subtract(const vec3& a_v) {
0188     m_data[0] -= a_v.m_data[0];
0189     m_data[1] -= a_v.m_data[1];
0190     m_data[2] -= a_v.m_data[2];
0191   }
0192 
0193   void subtract(const T& a0,const T& a1,const T& a2) {
0194     m_data[0] -= a0;
0195     m_data[1] -= a1;
0196     m_data[2] -= a2;
0197   }
0198 
0199   bool cos_angle(const vec3& a_v,T& a_cos,T(*a_sqrt)(T)) const {
0200     //WARNING : if ret false, a_cos is not set.
0201     T this_length = length(a_sqrt);
0202     if(this_length==T()) return false;
0203     T a_v_length = a_v.length(a_sqrt);
0204     if(a_v_length==T()) return false;
0205     a_cos = dot(a_v)/(this_length*a_v_length);
0206     return true;
0207   }
0208 
0209   bool theta_phi(T& a_theta,T& a_phi,T(*a_sqrt)(T),T(*a_atan2)(T,T)) const {
0210     //WARNING : if ret false, a_theta, a_phi are not set.
0211     if(length(a_sqrt)==T()) return false;
0212     a_phi = a_atan2(m_data[1],m_data[0]);
0213     T xy = a_sqrt(m_data[0]*m_data[0]+m_data[1]*m_data[1]);
0214     a_theta = a_atan2(xy,m_data[2]);
0215     return true;
0216   }
0217 
0218 public: //operators
0219   T& operator[](size_t a_index) {
0220     //WARNING : no check on a_index.
0221     return m_data[a_index];
0222   }
0223   const T& operator[](size_t a_index) const {
0224     //WARNING : no check on a_index.
0225     return m_data[a_index];
0226   }
0227 
0228   vec3& operator*=(const T& a_v) {
0229     m_data[0] *= a_v;
0230     m_data[1] *= a_v;
0231     m_data[2] *= a_v;
0232     return *this;
0233   }
0234 
0235   vec3 operator+(const vec3& a_v) const {
0236     return vec3(m_data[0]+a_v.m_data[0],
0237                 m_data[1]+a_v.m_data[1],
0238                 m_data[2]+a_v.m_data[2]);
0239   }
0240 
0241   vec3 operator-(const vec3& a_v) const {
0242     return vec3(m_data[0]-a_v.m_data[0],
0243                 m_data[1]-a_v.m_data[1],
0244                 m_data[2]-a_v.m_data[2]);
0245   }
0246 
0247   vec3 operator*(const T& a_v) const {
0248     return vec3(m_data[0]*a_v,
0249                 m_data[1]*a_v,
0250                 m_data[2]*a_v);
0251   }
0252 
0253   vec3 operator/(const T& a_v) const {
0254     if(a_v==T()) return vec3();
0255     return vec3(m_data[0]/a_v,
0256                 m_data[1]/a_v,
0257                 m_data[2]/a_v);
0258   }
0259 
0260   bool operator==(const vec3& a_v) const {return equal(a_v);}
0261   bool operator!=(const vec3& a_v) const {return !operator==(a_v);}
0262 
0263 public: //for tools/sg/sf_vec
0264   typedef unsigned int size_type;
0265   size_type size() const {return 3;}
0266   const T* data() const {return m_data;}
0267   size_type data_size() const {return 3;} //for eqT.
0268 public: //for iv2sg
0269   const T* getValue() const {return m_data;}
0270   void setValue(const T& a0,const T& a1,const T& a2) {
0271     m_data[0] = a0;
0272     m_data[1] = a1;
0273     m_data[2] = a2;
0274   }
0275   void getValue(T& a0,T& a1,T& a2) const {
0276     a0 = m_data[0];
0277     a1 = m_data[1];
0278     a2 = m_data[2];
0279   }
0280   void setValue(const vec3& a_v) {
0281     m_data[0] = a_v.m_data[0];
0282     m_data[1] = a_v.m_data[1];
0283     m_data[2] = a_v.m_data[2];
0284   }
0285   void setValue(const T aV[3]) {
0286     m_data[0] = aV[0];
0287     m_data[1] = aV[1];
0288     m_data[2] = aV[2];
0289   }
0290 
0291   vec3& setValue(const vec3& a_bary,
0292                  const vec3& a_v0,const vec3& a_v1,const vec3& a_v2) {
0293     m_data[0] = a_bary[0]*a_v0[0]+a_bary[1]*a_v1[0]+a_bary[2]*a_v2[0];
0294     m_data[1] = a_bary[0]*a_v0[1]+a_bary[1]*a_v1[1]+a_bary[2]*a_v2[1];
0295     m_data[2] = a_bary[0]*a_v0[2]+a_bary[1]*a_v1[2]+a_bary[2]*a_v2[2];
0296     return *this;
0297   }
0298 
0299 public:
0300 #if defined(TOOLS_MEM) && !defined(TOOLS_MEM_ATEXIT)
0301   static const vec3<T>& s_x() {static const vec3<T> s_v(1,0,0,false);return s_v;}
0302   static const vec3<T>& s_y() {static const vec3<T> s_v(0,1,0,false);return s_v;}
0303   static const vec3<T>& s_z() {static const vec3<T> s_v(0,0,1,false);return s_v;}
0304 #else
0305   static const vec3<T>& s_x() {static const vec3<T> s_v(1,0,0);return s_v;}
0306   static const vec3<T>& s_y() {static const vec3<T> s_v(0,1,0);return s_v;}
0307   static const vec3<T>& s_z() {static const vec3<T> s_v(0,0,1);return s_v;}
0308 #endif
0309 protected:
0310   T m_data[3];
0311 
0312 private:static void check_instantiation() {vec3<float> v;}
0313 };
0314 
0315 //for sf, mf :
0316 template <class T>
0317 inline const T* get_data(const vec3<T>& a_v) {return a_v.data();}
0318 
0319 template <class T>
0320 inline void get_normal(const vec3<T>& a_p0,const vec3<T>& a_p1,const vec3<T>& a_p2,vec3<T>& a_nm,
0321                        vec3<T>& a_tmp_1,vec3<T>& a_tmp_2,T(*a_sqrt)(T)) {
0322   // Used to optimize sg::bin().
0323   //(a_p1-a_p0).cross(a_p2-a_p1,a_nm);
0324 
0325   a_tmp_1 = a_p1;
0326   a_tmp_1.subtract(a_p0);
0327 
0328   a_tmp_2 = a_p2;
0329   a_tmp_2.subtract(a_p1);
0330 
0331   a_tmp_1.cross(a_tmp_2,a_nm);
0332 
0333   a_nm.normalize(a_sqrt);
0334 }
0335 
0336 /*
0337 template <class VEC3>
0338 inline void get_normal(const VEC3& a_p0,const VEC3& a_p1,const VEC3& a_p2,VEC3& a_nm) {
0339   VEC3 tmp1,tmp2;
0340   get_normal(a_p0,a_p1,a_p2,a_nm,tmp1,tmp2);
0341 }
0342 */
0343 template <class VEC3>
0344 inline void direction(const VEC3& a_p0,const VEC3& a_p1,const VEC3& a_p2,VEC3& a_value) {
0345   // Orientation is computed by taking (p1 - p0) x (p2 - p0)
0346   VEC3 P = a_p1;
0347   P.subtract(a_p0);
0348   VEC3 P2 = a_p2;
0349   P2.subtract(a_p0);
0350   P.cross(P2,a_value);
0351 }
0352 
0353 template <class VEC3>
0354 inline void area(const VEC3& a_p0,const VEC3& a_p1,const VEC3& a_p2,typename VEC3::elem_t& a_value,
0355                  VEC3& a_tmp_1,VEC3& a_tmp_2,VEC3& a_tmp_3) {
0356   // area of the triangle (a_p0,a_p1,a_p2)
0357   typedef typename VEC3::elem_t T;
0358 
0359   a_tmp_1 = a_p1;
0360   a_tmp_1.subtract(a_p0);
0361 
0362   a_tmp_2 = a_p2;
0363   a_tmp_2.subtract(a_p1);
0364 
0365   a_tmp_1.cross(a_tmp_2,a_tmp_3);
0366 
0367   a_value = a_tmp_3.length()/T(2);
0368 }
0369 
0370 template <class VEC3>
0371 inline void area(const VEC3& a_p0,const VEC3& a_p1,const VEC3& a_p2,typename VEC3::elem_t& a_value) {
0372   VEC3 tmp1,tmp2,tmp3;
0373   area(a_p0,a_p1,a_p2,a_value,tmp1,tmp2,tmp3);
0374 }
0375 
0376 template <class T>
0377 inline void direction(const T& a_0_x,const T& a_0_y,const T& a_0_z,
0378                       const T& a_1_x,const T& a_1_y,const T& a_1_z,
0379                       const T& a_2_x,const T& a_2_y,const T& a_2_z,vec3<T>& a_value) {
0380   direction(vec3<T>(a_0_x,a_0_y,a_0_z),
0381             vec3<T>(a_1_x,a_1_y,a_1_z),
0382             vec3<T>(a_2_x,a_2_y,a_2_z),a_value);
0383 }
0384 
0385 }
0386 
0387 #include <ostream>
0388 
0389 namespace tools {
0390 
0391 // for sf_vec::dump().
0392 template <class T>
0393 inline std::ostream& operator<<(std::ostream& a_out,const vec3<T>& a_this){
0394   a_out << "x = " << a_this.v0()
0395         << ",y = " << a_this.v1()
0396         << ",z = " << a_this.v2();
0397   return a_out;
0398 }
0399 
0400 }
0401 
0402 #endif