00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #ifndef GEOM_VECMAT3_OP_HH
00020 #define GEOM_VECMAT3_OP_HH
00021
00022 #include <ostream>
00023 #include "constants.hh"
00024
00025 #include <ost/geom/module_config.hh>
00026 #include <ost/geom/vec3.hh>
00027 #include <ost/geom/mat3.hh>
00028
00029 namespace geom {
00030
00031
00033 inline Real Length2(const Vec3& v)
00034 {
00035 return v[0]*v[0]+v[1]*v[1]+v[2]*v[2];
00036 }
00037
00039 inline Real Length(const Vec3& v)
00040 {
00041 return std::sqrt(Length2(v));
00042 }
00043
00045 inline bool Equal(const Vec3& v1, const Vec3& v2,
00046 Real ephilon=EPSILON)
00047 {
00048 return std::fabs(v1[0]-v2[0])<ephilon &&
00049 std::fabs(v1[1]-v2[1])<ephilon &&
00050 std::fabs(v1[2]-v2[2])<ephilon;
00051 }
00052
00054 inline bool Equal(const Mat3& m1, const Mat3& m2,
00055 Real ephilon=EPSILON)
00056 {
00057 return std::fabs(m1(0,0)-m2(0,0))<ephilon &&
00058 std::fabs(m1(0,1)-m2(0,1))<ephilon &&
00059 std::fabs(m1(0,2)-m2(0,2))<ephilon &&
00060 std::fabs(m1(1,0)-m2(1,0))<ephilon &&
00061 std::fabs(m1(1,1)-m2(1,1))<ephilon &&
00062 std::fabs(m1(1,2)-m2(1,2))<ephilon &&
00063 std::fabs(m1(2,0)-m2(2,0))<ephilon &&
00064 std::fabs(m1(2,1)-m2(2,1))<ephilon &&
00065 std::fabs(m1(2,2)-m2(2,2))<ephilon;
00066 }
00067
00069 inline Real Dot(const Vec3& v1, const Vec3& v2)
00070 {
00071 return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];
00072 }
00073
00075 inline Vec3 Normalize(const Vec3& v)
00076 {
00077 Real l=Length(v);
00078 if(l==0.0) {
00079 return v;
00080 }
00081 return v/l;
00082 }
00083
00085 inline Vec3 Cross(const Vec3& v1, const Vec3& v2)
00086 {
00087 Vec3 nrvo(v1[1]*v2[2]-v2[1]*v1[2],
00088 v1[2]*v2[0]-v2[2]*v1[0],
00089 v1[0]*v2[1]-v2[0]*v1[1]);
00090 return nrvo;
00091 }
00092
00094 inline Vec3 CompMultiply(const Vec3& v1, const Vec3& v2)
00095 {
00096 Vec3 nrvo(v1[0]*v2[0],v1[1]*v2[1],v1[2]*v2[2]);
00097 return nrvo;
00098 }
00099
00101 inline Vec3 CompDivide(const Vec3& v1, const Vec3& v2)
00102 {
00103 Vec3 nrvo(v1[0]/v2[0],v1[1]/v2[1],v1[2]/v2[2]);
00104 return nrvo;
00105 }
00106
00108 inline Vec3 operator*(const Vec3& v,const Mat3& m)
00109 {
00110 Vec3 nrvo(v[0]*m(0,0)+v[1]*m(1,0)+v[2]*m(2,0),
00111 v[0]*m(0,1)+v[1]*m(1,1)+v[2]*m(2,1),
00112 v[0]*m(0,2)+v[1]*m(1,2)+v[2]*m(2,2));
00113 return nrvo;
00114 }
00115
00117 inline Vec3 operator*(const Mat3& m, const Vec3& v)
00118 {
00119 Vec3 nrvo(v[0]*m(0,0)+v[1]*m(0,1)+v[2]*m(0,2),
00120 v[0]*m(1,0)+v[1]*m(1,1)+v[2]*m(1,2),
00121 v[0]*m(2,0)+v[1]*m(2,1)+v[2]*m(2,2));
00122 return nrvo;
00123 }
00124
00126 inline Mat3 operator*(const Mat3& m1, const Mat3& m2)
00127 {
00128 Mat3 nrvo(m1(0,0)*m2(0,0)+m1(0,1)*m2(1,0)+m1(0,2)*m2(2,0),
00129 m1(0,0)*m2(0,1)+m1(0,1)*m2(1,1)+m1(0,2)*m2(2,1),
00130 m1(0,0)*m2(0,2)+m1(0,1)*m2(1,2)+m1(0,2)*m2(2,2),
00131
00132 m1(1,0)*m2(0,0)+m1(1,1)*m2(1,0)+m1(1,2)*m2(2,0),
00133 m1(1,0)*m2(0,1)+m1(1,1)*m2(1,1)+m1(1,2)*m2(2,1),
00134 m1(1,0)*m2(0,2)+m1(1,1)*m2(1,2)+m1(1,2)*m2(2,2),
00135
00136 m1(2,0)*m2(0,0)+m1(2,1)*m2(1,0)+m1(2,2)*m2(2,0),
00137 m1(2,0)*m2(0,1)+m1(2,1)*m2(1,1)+m1(2,2)*m2(2,1),
00138 m1(2,0)*m2(0,2)+m1(2,1)*m2(1,2)+m1(2,2)*m2(2,2));
00139 return nrvo;
00140 }
00141
00142 Mat3 DLLEXPORT_OST_GEOM Invert(const Mat3& m);
00143 Mat3 DLLEXPORT_OST_GEOM Transpose(const Mat3& m);
00144
00145
00146 Real DLLEXPORT_OST_GEOM Comp(const Mat3& m, unsigned int i, unsigned int j);
00147
00148
00149 Real DLLEXPORT_OST_GEOM Minor(const Mat3& m, unsigned int i, unsigned int j);
00150
00151
00152 Real DLLEXPORT_OST_GEOM Det(const Mat3& m);
00153
00154
00155 Real DLLEXPORT_OST_GEOM Angle(const Vec3& v1, const Vec3& v2);
00156
00157
00158 Real DLLEXPORT_OST_GEOM SignedAngle(const Vec3& v1, const Vec3& v2, const Vec3& ref);
00159
00160 Mat3 DLLEXPORT_OST_GEOM EulerTransformation(Real theta, Real phi, Real xi);
00161
00162 Mat3 DLLEXPORT_OST_GEOM AxisRotation(const Vec3& axis, Real angle);
00163
00167 Vec3 DLLEXPORT_OST_GEOM OrthogonalVector(const Vec3& axis);
00168
00170 inline Real DihedralAngle(const Vec3& p1, const Vec3& p2,
00171 const Vec3& p3, const Vec3& p4) {
00172 const Vec3 r1 = p2-p1;
00173 const Vec3 r2 = p3-p2;
00174 const Vec3 r3 = p4-p3;
00175 const Vec3 r12cross = Cross(r1, r2);
00176 const Vec3 r23cross = Cross(r2, r3);
00177 return std::atan2(Dot(r1*Length(r2), r23cross), Dot(r12cross, r23cross));
00178 }
00179
00181 inline Vec3 Min(const Vec3& v1, const Vec3& v2)
00182 {
00183 Vec3 nrvo(std::min(v1[0],v2[0]),
00184 std::min(v1[1],v2[1]),
00185 std::min(v1[2],v2[2]));
00186 return nrvo;
00187 }
00188
00190 inline Vec3 Max(const Vec3& v1, const Vec3& v2)
00191 {
00192 Vec3 nrvo(std::max(v1[0],v2[0]),
00193 std::max(v1[1],v2[1]),
00194 std::max(v1[2],v2[2]));
00195 return nrvo;
00196 }
00197
00199 inline Real Distance(const Vec3& p1, const Vec3& p2)
00200 {
00201 return Length(p1-p2);
00202 }
00203
00204
00206 inline Real Distance2WithPBC(const Vec3& v1, const Vec3& v2, const Vec3& ucell_size){
00207 Vec3 v;
00208 v=v1-v2;
00209 for (int i=0; i<3; i++) {
00210 if (std::fabs(v[i])>ucell_size[i]/2.){
00211 v[i]=std::fabs(v[i])-ucell_size[i]*int(std::fabs(v[i])/ucell_size[i]+0.5);
00212 }
00213 }
00214 return Length2(v);
00215 }
00217 inline Real DistanceWithPBC(const Vec3& v1, const Vec3& v2, const Vec3& ucell_size){
00218 return sqrt(Distance2WithPBC(v1, v2, ucell_size));
00219 }
00221 Real DLLEXPORT_OST_GEOM MinDistance(const Vec3List& l1, const Vec3List& l2);
00223
00224 Real DLLEXPORT_OST_GEOM MinDistanceWithPBC(const Vec3List& l1, const Vec3List& l2, Vec3& ucell_size);
00227 std::vector<unsigned int> DLLEXPORT_OST_GEOM MinDistanceIndices(const Vec3List& l1, const Vec3List& l2);
00229 Vec3List DLLEXPORT_OST_GEOM CalculateUnitCellVectors(const Vec3& ucell_size, const Vec3& ucell_angles);
00231 Vec3 DLLEXPORT_OST_GEOM WrapVec3(const Vec3& v1,const Vec3& box_center,const Vec3& ucell_size);
00233 Vec3List DLLEXPORT_OST_GEOM WrapVec3List(const Vec3List& vl,const Vec3& box_center,const Vec3& ucell_size);
00235 Vec3 DLLEXPORT_OST_GEOM WrapVec3(const Vec3& v1,const Vec3& box_center,const Vec3& ucell_size,const Vec3& ucell_angles);
00237 Vec3List DLLEXPORT_OST_GEOM WrapVec3List(const Vec3List& vl,const Vec3& box_center,const Vec3& ucell_size,const Vec3& ucell_angles);
00238
00239 }
00240
00241 #endif