ESyS-Particle  4.0.1
vec3.hpp
00001 
00002 //                                                         //
00003 // Copyright (c) 2003-2011 by The University of Queensland //
00004 // Earth Systems Science Computational Centre (ESSCC)      //
00005 // http://www.uq.edu.au/esscc                              //
00006 //                                                         //
00007 // Primary Business: Brisbane, Queensland, Australia       //
00008 // Licensed under the Open Software License version 3.0    //
00009 // http://www.opensource.org/licenses/osl-3.0.php          //
00010 //                                                         //
00012 
00013 
00014 #ifndef __VEC3_HPP
00015 #define __VEC3_HPP
00016 
00017 #include "Foundation/Matrix3.h"
00018 
00019 //the error...
00020 VEC3_INLINE VecErr::VecErr(const string& m):MError(m)
00021 {
00022   message.insert(0,"Vec3 "); 
00023 }
00024 
00025 // constructors
00026 VEC3_INLINE Vec3::Vec3()
00027 {
00028   data[0]=0;
00029   data[1]=0;
00030   data[2]=0;
00031 }
00032 
00033 VEC3_INLINE Vec3::Vec3(double s)
00034 {
00035   data[0]=s;
00036   data[1]=s;
00037   data[2]=s;
00038 }
00039 
00040 VEC3_INLINE Vec3::Vec3(double a,double b,double c)
00041 {
00042   data[0]=a;
00043   data[1]=b;
00044   data[2]=c;
00045 }
00046 
00047 VEC3_INLINE Vec3::Vec3(const Vec3& rhs)
00048 {
00049   data[0]=rhs.data[0];
00050   data[1]=rhs.data[1];
00051   data[2]=rhs.data[2];
00052 }
00053 
00054 // operators
00055 
00056 VEC3_INLINE Vec3& Vec3::operator=(const Vec3& rhs)
00057 {
00058   data[0]=rhs.data[0];
00059   data[1]=rhs.data[1];
00060   data[2]=rhs.data[2];
00061   return *this;
00062 }
00063 
00064 VEC3_INLINE Vec3& Vec3::operator=(double s)
00065 {
00066   data[0]=s;
00067   data[1]=s;
00068   data[2]=s;
00069   return *this;
00070 }
00071 
00072 VEC3_INLINE Vec3& Vec3::operator-=(const Vec3& rhs)
00073 {
00074   data[0]-=rhs.data[0];
00075   data[1]-=rhs.data[1];
00076   data[2]-=rhs.data[2];
00077   return *this;
00078 }
00079 
00080 VEC3_INLINE Vec3& Vec3::operator+=(const Vec3& rhs)
00081 {
00082   data[0]+=rhs.data[0];
00083   data[1]+=rhs.data[1];
00084   data[2]+=rhs.data[2];
00085   return *this;
00086 }
00087 
00088 VEC3_INLINE Vec3 Vec3::operator+(const Vec3& rhs) const
00089 {
00090   return Vec3(data[0]+rhs.data[0], data[1]+rhs.data[1], data[2]+rhs.data[2]);
00091 }
00092 
00093 VEC3_INLINE Vec3 Vec3::operator-(const Vec3& rhs) const
00094 {
00095   return Vec3(data[0]-rhs.data[0], data[1]-rhs.data[1], data[2]-rhs.data[2]);
00096 }
00097 
00098 VEC3_INLINE Vec3 Vec3::operator-() const
00099 {
00100   return Vec3( -data[0],-data[1],-data[2] );
00101 }
00102 
00103 VEC3_INLINE Vec3 Vec3::operator*(const Matrix3 &m) const
00104 {
00105   const double x = m(0,0)*data[0] + m(1,0)*data[1] + m(2,0)*data[2];
00106   const double y = m(0,1)*data[0] + m(1,1)*data[1] + m(2,1)*data[2];
00107   const double z = m(0,2)*data[0] + m(1,2)*data[1] + m(2,2)*data[2];
00108 
00109   return Vec3(x,y,z);
00110 }
00111 
00112 VEC3_INLINE double Vec3::operator*(const Vec3& rhs) const
00113 {
00114   return data[0]*rhs.data[0]+data[1]*rhs.data[1]+data[2]*rhs.data[2];
00115 }
00116 
00117 VEC3_INLINE Vec3 Vec3::operator*(double s) const 
00118 { 
00119    return Vec3(data[0]*s,data[1]*s,data[2]*s) ; 
00120 } 
00121 
00122 VEC3_INLINE Vec3& Vec3::operator*=(double rhs)
00123 {
00124   data[0]*=rhs;
00125   data[1]*=rhs;
00126   data[2]*=rhs;
00127   return *this;
00128 }
00129 
00130 VEC3_INLINE Vec3& Vec3::operator/=(double c)
00131 {
00132   data[0] /= c;
00133   data[1] /= c;
00134   data[2] /= c;
00135 
00136   return *this;
00137 }
00138 
00139 VEC3_INLINE Vec3 Vec3::operator/(double s) const 
00140 { 
00141    return Vec3(data[0]/s,data[1]/s,data[2]/s) ; 
00142 }
00143 
00144 VEC3_INLINE Vec3 Vec3::operator+(double s) const 
00145 {
00146    return Vec3(data[0]+s, data[1]+s, data[2]+s) ; 
00147 }
00148 
00149 VEC3_INLINE Vec3 Vec3::operator-(double s) const 
00150 {
00151    return Vec3(data[0]-s, data[1]-s, data[2]-s); 
00152 }
00153 
00154 VEC3_INLINE Vec3 Vec3::rotate(const Vec3 &axis, const Vec3 &axisPt) const
00155 {
00156   const double phi = axis.norm();
00157   if (phi > 0.0)
00158   {
00159     const Vec3 r = *this - axisPt;
00160     const Vec3 n = axis/phi;
00161     const double cosPhi = cos(phi);
00162     const Vec3 rotatedR =
00163       r*cosPhi + n*((dot(n, r))*(1-cosPhi)) + cross(r, n)*sin(phi);
00164     return rotatedR + axisPt;
00165   }
00166   return *this;
00167 }
00168 
00169 VEC3_INLINE Vec3 &Vec3::operator+=(double s)
00170 {
00171    data[0] += s;
00172    data[1] += s;
00173    data[2] += s;
00174    return *this; 
00175 }
00176 
00177 VEC3_INLINE Vec3 &Vec3::operator-=(double s)
00178 {
00179    data[0] -= s;
00180    data[1] -= s;
00181    data[2] -= s;
00182    return *this; 
00183 }
00184 
00185 // vector product
00186 // 9 Flops ( 6 mult, 3 sub ) 
00187 VEC3_INLINE Vec3 cross(const Vec3& lhs,const Vec3& rhs)
00188 {
00189   return Vec3(lhs.data[1]*rhs.data[2]-lhs.data[2]*rhs.data[1],
00190               lhs.data[2]*rhs.data[0]-lhs.data[0]*rhs.data[2],
00191               lhs.data[0]*rhs.data[1]-lhs.data[1]*rhs.data[0]);
00192 }
00193 
00194 //  dot product  
00195 
00196 VEC3_INLINE double dot(const Vec3& v1, const Vec3& v2)
00197 {
00198   return v1.data[0] * v2.data[0] + 
00199          v1.data[1] * v2.data[1] + 
00200          v1.data[2] * v2.data[2];
00201 }
00202 
00203 VEC3_INLINE Vec3 operator*(double f,const Vec3& rhs)
00204 {
00205   return Vec3(f*rhs.data[0], f*rhs.data[1], f*rhs.data[2]);
00206 }
00207 
00208 
00209 // euclidian norm
00210 // 6 Flops ( 3 mult, 2 add, 1 sqrt )
00211 VEC3_INLINE double Vec3::norm() const
00212 {
00213   return sqrt(data[0]*data[0]+data[1]*data[1]+data[2]*data[2]);
00214 }
00215 
00216 // square of the euclidian norm
00217 // 5 Flops ( 3 mult, 2 add)
00218 VEC3_INLINE double Vec3::norm2() const
00219 {
00220   return data[0]*data[0]+data[1]*data[1]+data[2]*data[2];
00221 }
00222 
00223 // returns unit vector in direction of the original vector
00224 // 9 Flops ( 3 mult, 2 add, 3 div, 1 sqrt ) 
00225 VEC3_INLINE Vec3 Vec3::unit() const
00226 {
00227   return (*this)/norm();
00228 }
00229 
00230 // per element min/max
00231 VEC3_INLINE Vec3 cmax(const Vec3& v1,const Vec3& v2)
00232 {
00233   Vec3 res;
00234   res.data[0]=v1.data[0]>v2.data[0] ? v1.data[0] : v2.data[0];
00235   res.data[1]=v1.data[1]>v2.data[1] ? v1.data[1] : v2.data[1];
00236   res.data[2]=v1.data[2]>v2.data[2] ? v1.data[2] : v2.data[2];
00237   return res;
00238 }
00239 
00240 VEC3_INLINE Vec3 cmin(const Vec3& v1,const Vec3& v2)
00241 {
00242   Vec3 res;
00243   res.data[0]=v1.data[0]<v2.data[0] ? v1.data[0] : v2.data[0];
00244   res.data[1]=v1.data[1]<v2.data[1] ? v1.data[1] : v2.data[1];
00245   res.data[2]=v1.data[2]<v2.data[2] ? v1.data[2] : v2.data[2];
00246   return res;
00247 }
00248 
00249 // save version, throws exception if norm()==0
00250 VEC3_INLINE Vec3 Vec3::unit_s() const
00251 {
00252   double n=norm();
00253   if(n==0) throw VecErr("norm() of data[2]ero-vector"); 
00254   Vec3 res(data[0],data[1],data[2]);
00255   return res/n;
00256 }
00257 
00258 VEC3_INLINE double Vec3::max() const
00259 {
00260   double m = ( data[0]>data[1] ? data[0] : data[1] );
00261 
00262   return ( m>data[2] ? m : data[2] );
00263 
00264 }
00265 
00266 VEC3_INLINE double Vec3::min() const
00267 {
00268   double m = ( data[0]<data[1] ? data[0] : data[1] );
00269 
00270   return ( m<data[2] ? m : data[2] );
00271 }
00272 
00273 VEC3_INLINE bool Vec3::operator==(const Vec3& V) const
00274 {
00275   return((data[0]==V.data[0])&&(data[1]==V.data[1])&&(data[2]==V.data[2]));
00276 }
00277 
00278 VEC3_INLINE bool Vec3::operator!=(const Vec3& V) const
00279 {
00280   return((data[0]!=V.data[0])||(data[1]!=V.data[1])||(data[2]!=V.data[2]));
00281 }
00282 
00283 // per component min/max
00284 
00285 VEC3_INLINE Vec3 comp_max(const Vec3& V1,const Vec3& V2)
00286 {
00287   double x=(V1.X() > V2.X()) ? V1.X() : V2.X();
00288   double y=(V1.Y() > V2.Y()) ? V1.Y() : V2.Y();
00289   double z=(V1.Z() > V2.Z()) ? V1.Z() : V2.Z();
00290  
00291   return Vec3(x,y,z);
00292 }
00293 
00294 VEC3_INLINE Vec3 comp_min(const Vec3& V1,const Vec3& V2)
00295 {
00296   double x=(V1.X() < V2.X()) ? V1.X() : V2.X();
00297   double y=(V1.Y() < V2.Y()) ? V1.Y() : V2.Y();
00298   double z=(V1.Z() < V2.Z()) ? V1.Z() : V2.Z();
00299  
00300   return Vec3(x,y,z);
00301 }
00302 
00303 // in/output
00304 
00305 VEC3_INLINE std::ostream& operator << (std::ostream& ostr,const Vec3& V)
00306 {
00307   const char delimiter = ' ';
00308   ostr
00309     << V.data[0] << delimiter
00310     << V.data[1] << delimiter
00311     << V.data[2];
00312 
00313   return ostr;
00314 }
00315 
00316 VEC3_INLINE std::istream& operator >> (std::istream& istr,Vec3& V)
00317 {
00318   istr
00319     >> V.data[0]
00320     >> V.data[1]
00321     >> V.data[2];
00322 
00323   return istr;
00324 }
00325 
00326 #endif // __VEC3_HPP