diff --git a/src/diffpy/Attributes.hpp b/src/diffpy/Attributes.hpp index 539ec67a..9f31be3c 100644 --- a/src/diffpy/Attributes.hpp +++ b/src/diffpy/Attributes.hpp @@ -38,9 +38,9 @@ class Attributes; class DoubleAttributeError : public std::runtime_error { public: - DoubleAttributeError(const std::string msg="") : + DoubleAttributeError(const std::string msg="") : std::runtime_error(msg) - { } + { } }; /// @class BaseDoubleAttribute diff --git a/src/diffpy/mathutils.ipp b/src/diffpy/mathutils.ipp index 626f4592..c95d112a 100644 --- a/src/diffpy/mathutils.ipp +++ b/src/diffpy/mathutils.ipp @@ -51,17 +51,17 @@ double cosd(double x) double xp = fmod(fabs(x), 360.0); if (remainder(xp, 60.0) == 0.0 || remainder(xp, 90.0) == 0.0) { - switch(int(round(xp))) - { - case 0: return 1.0; - case 60: - case 300: return 0.5; - case 90: - case 270: return 0.0; - case 120: - case 240: return -0.5; - case 180: return -1.0; - }; + switch(int(round(xp))) + { + case 0: return 1.0; + case 60: + case 300: return 0.5; + case 90: + case 270: return 0.0; + case 120: + case 240: return -0.5; + case 180: return -1.0; + }; } return cos(x/180.0*M_PI); } @@ -79,14 +79,14 @@ double acosd(double x) { if (remainder(x, 0.5) == 0.0) { - switch(int(round(x/0.5))) - { - case 0: return 90.0; - case 1: return 60.0; - case -1: return 120.0; - case 2: return 0.0; - case -2: return 180.0; - }; + switch(int(round(x/0.5))) + { + case 0: return 90.0; + case 1: return 60.0; + case -1: return 120.0; + case 2: return 0.0; + case -2: return 180.0; + }; } return acos(x)/M_PI*180.0; } @@ -97,14 +97,14 @@ double asind(double x) { if (remainder(x, 0.5) == 0.0) { - switch(int(round(x/0.5))) - { - case 0: return 0.0; - case 1: return 30.0; - case -1: return -30.0; - case 2: return 90.0; - case -2: return -90.0; - }; + switch(int(round(x/0.5))) + { + case 0: return 0.0; + case 1: return 30.0; + case -1: return -30.0; + case 2: return 90.0; + case -2: return -90.0; + }; } return acos(x)/M_PI*180.0; } diff --git a/src/diffpy/srreal/Lattice.cpp b/src/diffpy/srreal/Lattice.cpp index 6416d977..7f920100 100644 --- a/src/diffpy/srreal/Lattice.cpp +++ b/src/diffpy/srreal/Lattice.cpp @@ -31,7 +31,7 @@ Lattice::Lattice() Lattice::Lattice(double a0, double b0, double c0, - double alpha0, double beta0, double gamma0) + double alpha0, double beta0, double gamma0) { mbaserot = 1.0, 0.0, 0.0, @@ -43,7 +43,7 @@ Lattice::Lattice(double a0, double b0, double c0, // Public Methods ------------------------------------------------------------ void Lattice::setLatPar(double a0, double b0, double c0, - double alpha0, double beta0, double gamma0) + double alpha0, double beta0, double gamma0) { using namespace diffpy::mathutils; ma = a0; @@ -97,8 +97,8 @@ void Lattice::setLatPar(double a0, double b0, double c0, } void Lattice::setLatBase(const R3::Vector& va0, - const R3::Vector& vb0, - const R3::Vector& vc0) + const R3::Vector& vb0, + const R3::Vector& vc0) { using namespace diffpy::mathutils; mbase = va0[0], va0[1], va0[2], @@ -221,8 +221,8 @@ const R3::Matrix& Lattice::cartesianMatrix(const R3::Matrix& Ml) const const R3::Matrix& Lattice::fractionalMatrix(const R3::Matrix& Mc) const { static R3::Matrix res0, res1; - res0 = product(Mc, mrecnormbase); - res1 = product(R3::transpose(mrecnormbase), res0); + res0 = R3::product(Mc, mrecnormbase); + res1 = R3::product(R3::transpose(mrecnormbase), res0); return res1; } diff --git a/src/diffpy/srreal/PointsInSphere.cpp b/src/diffpy/srreal/PointsInSphere.cpp index 01695a3b..51d7a057 100644 --- a/src/diffpy/srreal/PointsInSphere.cpp +++ b/src/diffpy/srreal/PointsInSphere.cpp @@ -32,9 +32,9 @@ using pointsinsphere::LatticeParameters; // Constructor --------------------------------------------------------------- LatticeParameters::LatticeParameters( double _a, double _b, double _c, - double _alpha, double _beta, double _gamma ) : - a(_a), b(_b), c(_c), - alpha(_alpha), beta(_beta), gamma(_gamma) + double _alpha, double _beta, double _gamma ) : + a(_a), b(_b), c(_c), + alpha(_alpha), beta(_beta), gamma(_gamma) { update(); } @@ -86,9 +86,9 @@ LatticeParameters LatticeParameters::reciprocal() const // Constructors -------------------------------------------------------------- PointsInSphere::PointsInSphere(double rmin, double rmax, - const LatticeParameters& _latpar ) : + const LatticeParameters& _latpar ) : _Rmin(rmin), _Rmax(rmax), - latpar(_latpar), + latpar(_latpar), _m(_mno[0]), _n(_mno[1]), _o(_mno[2]) { init(); @@ -97,10 +97,10 @@ PointsInSphere::PointsInSphere(double rmin, double rmax, PointsInSphere::PointsInSphere(double rmin, double rmax, - double _a, double _b, double _c, - double _alpha, double _beta, double _gamma) : - _Rmin(rmin), _Rmax(rmax), - latpar(_a, _b, _c, _alpha, _beta, _gamma), + double _a, double _b, double _c, + double _alpha, double _beta, double _gamma) : + _Rmin(rmin), _Rmax(rmax), + latpar(_a, _b, _c, _alpha, _beta, _gamma), _m(_mno[0]), _n(_mno[1]), _o(_mno[2]) { init(); @@ -181,7 +181,7 @@ double PointsInSphere::r() const const double &a = latpar.a, &b = latpar.b, &c = latpar.c; const double &ca = latpar.ca, &cb = latpar.cb, &cg = latpar.cg; return sqrt( m()*m()*a*a + n()*n()*b*b + o()*o()*c*c - + 2*m()*n()*a*b*cg + 2*m()*o()*a*c*cb + 2*n()*o()*b*c*ca ); + + 2*m()*n()*a*b*cg + 2*m()*o()*a*c*cb + 2*n()*o()*b*c*ca ); } // Private Methods ----------------------------------------------------------- @@ -191,7 +191,7 @@ void PointsInSphere::next_m() this->_m += 1; if (finished()) { - return; + return; } // not finished here n0plane = m()*dn0dm; @@ -207,26 +207,26 @@ void PointsInSphere::next_n() { do { - this->_n += 1; - if (n() < hi_n) - { - o0line = o0plane + (n()-n0plane)*do0dn; - double RlineSquare = RplaneSquare - pow((n()-n0plane)/b2r,2); - oHalfSpan = RlineSquare > 0.0 ? sqrt(RlineSquare)*c1r : 0.0; - // parentheses improve round-off errors around [0,0,0] - double RExclSquare = (RlineSquare - RmaxSquare) + RminSquare; - oExclHalfSpan = RExclSquare > 0.0 ? sqrt(RExclSquare)*c1r : 0.0; - this->_o = int(floor(o0line - oHalfSpan)); - outside_o = int(ceil(o0line + oHalfSpan)); - hi_o = outside_o; - if (oExclHalfSpan) - { - int hole_o = int(ceil(o0line - oExclHalfSpan)); - if (fabs(hole_o-o0line) < oExclHalfSpan) hi_o = hole_o; - } - return; - } - next_m(); + this->_n += 1; + if (n() < hi_n) + { + o0line = o0plane + (n()-n0plane)*do0dn; + double RlineSquare = RplaneSquare - pow((n()-n0plane)/b2r,2); + oHalfSpan = RlineSquare > 0.0 ? sqrt(RlineSquare)*c1r : 0.0; + // parentheses improve round-off errors around [0,0,0] + double RExclSquare = (RlineSquare - RmaxSquare) + RminSquare; + oExclHalfSpan = RExclSquare > 0.0 ? sqrt(RExclSquare)*c1r : 0.0; + this->_o = int(floor(o0line - oHalfSpan)); + outside_o = int(ceil(o0line + oHalfSpan)); + hi_o = outside_o; + if (oExclHalfSpan) + { + int hole_o = int(ceil(o0line - oExclHalfSpan)); + if (fabs(hole_o-o0line) < oExclHalfSpan) hi_o = hole_o; + } + return; + } + next_m(); } while (!finished()); } @@ -236,18 +236,18 @@ void PointsInSphere::next_o() { do { - this->_o += 1; - if (o() < hi_o) - { - return; - } - if (hi_o != outside_o) - { - hi_o = outside_o; - this->_o = int( ceil(o0line+oExclHalfSpan) ) - 1; - continue; - } - next_n(); + this->_o += 1; + if (o() < hi_o) + { + return; + } + if (hi_o != outside_o) + { + hi_o = outside_o; + this->_o = int( ceil(o0line+oExclHalfSpan) ) - 1; + continue; + } + next_n(); } while (!finished()); } @@ -275,19 +275,19 @@ void PointsInSphere::init() // Constructors -------------------------------------------------------------- ReflectionsInQminQmax::ReflectionsInQminQmax(double qmin, double qmax, - const LatticeParameters& _latpar) : - _Qmin(qmin), _Qmax(qmax), - latpar(_latpar), - sph(qmin*M_1_PI/2.0, qmax*M_1_PI/2.0, latpar.reciprocal()) + const LatticeParameters& _latpar) : + _Qmin(qmin), _Qmax(qmax), + latpar(_latpar), + sph(qmin*M_1_PI/2.0, qmax*M_1_PI/2.0, latpar.reciprocal()) { } ReflectionsInQminQmax::ReflectionsInQminQmax(double qmin, double qmax, - double _a, double _b, double _c, - double _alpha, double _beta, double _gamma ) : - _Qmin(qmin), _Qmax(qmax), - latpar(_a, _b, _c, _alpha, _beta, _gamma), - sph(qmin*M_1_PI/2.0, qmax*M_1_PI/2.0, latpar.reciprocal()) + double _a, double _b, double _c, + double _alpha, double _beta, double _gamma ) : + _Qmin(qmin), _Qmax(qmax), + latpar(_a, _b, _c, _alpha, _beta, _gamma), + sph(qmin*M_1_PI/2.0, qmax*M_1_PI/2.0, latpar.reciprocal()) { } // Public Methods ------------------------------------------------------------ @@ -367,18 +367,18 @@ double ReflectionsInQminQmax::d() const // Constructors -------------------------------------------------------------- ReflectionsInDmaxDmin::ReflectionsInDmaxDmin(double dmax, double dmin, - const LatticeParameters& _latpar) : - ReflectionsInQminQmax(2.0*M_PI/dmax, 2.0*M_PI/dmin, _latpar), + const LatticeParameters& _latpar) : + ReflectionsInQminQmax(2.0*M_PI/dmax, 2.0*M_PI/dmin, _latpar), _Dmin(dmin), _Dmax(dmax) { } ReflectionsInDmaxDmin::ReflectionsInDmaxDmin(double dmax, double dmin, - double _a, double _b, double _c, - double _alpha, double _beta, double _gamma) : - ReflectionsInQminQmax(2.0*M_PI/dmax, 2.0*M_PI/dmin, - _a, _b, _c, _alpha, _beta, _gamma), - _Dmin(dmin), _Dmax(dmax) + double _a, double _b, double _c, + double _alpha, double _beta, double _gamma) : + ReflectionsInQminQmax(2.0*M_PI/dmax, 2.0*M_PI/dmin, + _a, _b, _c, _alpha, _beta, _gamma), + _Dmin(dmin), _Dmax(dmax) { } // Public Methods ------------------------------------------------------------ diff --git a/src/diffpy/srreal/PointsInSphere.hpp b/src/diffpy/srreal/PointsInSphere.hpp index 6cfd569d..68d4333c 100644 --- a/src/diffpy/srreal/PointsInSphere.hpp +++ b/src/diffpy/srreal/PointsInSphere.hpp @@ -37,7 +37,7 @@ * * ReflectionsInQminQmax ref(Qmin, Qmax, a, b, c, alpha, beta, gamma) * for (ReflectionsInQminQmax ref(Qmin, Qmax, a, b, c, alpha, beta, gamma); -* !ref.finished(); ref.next() ) +* !ref.finished(); ref.next() ) * { * // Miller indices are in ref.h(), ref.k(), ref.l() or ref.hkl() * // ref.Q() is magnitude of Q vector @@ -233,10 +233,10 @@ PointsInSphere::PointsInSphere(double rmin, double rmax, const L& lat) : template ReflectionsInQminQmax::ReflectionsInQminQmax( double _qmin, double _qmax, const L& lat) : - _Qmin(_qmin), _Qmax(_qmax), - latpar(lat.a(), lat.b(), lat.c(), + _Qmin(_qmin), _Qmax(_qmax), + latpar(lat.a(), lat.b(), lat.c(), lat.alpha(), lat.beta(), lat.gamma()), - sph(_qmin*M_1_PI/2.0, _qmax*M_1_PI/2.0, latpar.reciprocal()) + sph(_qmin*M_1_PI/2.0, _qmax*M_1_PI/2.0, latpar.reciprocal()) { } // Template Constructor for ReflectionsInDmaxDmin ---------------------------- @@ -244,12 +244,12 @@ ReflectionsInQminQmax::ReflectionsInQminQmax( template ReflectionsInDmaxDmin::ReflectionsInDmaxDmin( double dmax, double dmin, const L& lat) : - ReflectionsInQminQmax(2.0*M_PI/dmax, 2.0*M_PI/dmin, lat), - _Dmin(dmin), _Dmax(dmax) + ReflectionsInQminQmax(2.0*M_PI/dmax, 2.0*M_PI/dmin, lat), + _Dmin(dmin), _Dmax(dmax) { } } // namespace srreal } // namespace diffpy -#endif // POINTSINSPHERE_HPP_INCLUDED +#endif // POINTSINSPHERE_HPP_INCLUDED diff --git a/src/diffpy/srreal/R3linalg.cpp b/src/diffpy/srreal/R3linalg.cpp index 9bfe7d41..9bc9b253 100644 --- a/src/diffpy/srreal/R3linalg.cpp +++ b/src/diffpy/srreal/R3linalg.cpp @@ -14,10 +14,11 @@ namespace diffpy { namespace srreal { +namespace R3 { -const R3::Matrix& R3::identity() +const Matrix& identity() { - static R3::Matrix mx; + static Matrix mx; static bool mx_ready = false; if (!mx_ready) { @@ -30,26 +31,26 @@ const R3::Matrix& R3::identity() } -const R3::Matrix& R3::zeros() +const Matrix& zeros() { - static R3::Matrix mx; + static Matrix mx; static bool mx_ready = false; if (!mx_ready) mx = 0.0; return mx; } -double R3::determinant(const R3::Matrix& A) +double determinant(const Matrix& A) { - gsl_matrix* gA = gsl_matrix_alloc(R3::Ndim, R3::Ndim); - for (int i = 0; i != R3::Ndim; ++i) + gsl_matrix* gA = gsl_matrix_alloc(Ndim, Ndim); + for (int i = 0; i != Ndim; ++i) { - for (int j = 0; j != R3::Ndim; ++j) - { - gsl_matrix_set(gA, i, j, A(i,j)); - } + for (int j = 0; j != Ndim; ++j) + { + gsl_matrix_set(gA, i, j, A(i,j)); + } } - gsl_permutation* gP = gsl_permutation_alloc(R3::Ndim); + gsl_permutation* gP = gsl_permutation_alloc(Ndim); int signum; gsl_linalg_LU_decomp(gA, gP, &signum); double det = gsl_linalg_LU_det(gA, signum); @@ -59,21 +60,21 @@ double R3::determinant(const R3::Matrix& A) } -R3::Matrix R3::inverse(const R3::Matrix& A) +Matrix inverse(const Matrix& A) { - gsl_matrix* gA = gsl_matrix_alloc(R3::Ndim, R3::Ndim); - for (int i = 0; i != R3::Ndim; ++i) + gsl_matrix* gA = gsl_matrix_alloc(Ndim, Ndim); + for (int i = 0; i != Ndim; ++i) { - for (int j = 0; j != R3::Ndim; ++j) - { - gsl_matrix_set(gA, i, j, A(i,j)); - } + for (int j = 0; j != Ndim; ++j) + { + gsl_matrix_set(gA, i, j, A(i,j)); + } } - gsl_permutation* gP = gsl_permutation_alloc(R3::Ndim); + gsl_permutation* gP = gsl_permutation_alloc(Ndim); int signum; gsl_linalg_LU_decomp(gA, gP, &signum); - R3::Matrix B; - gsl_matrix_view gB = gsl_matrix_view_array(B.data(), R3::Ndim, R3::Ndim); + Matrix B; + gsl_matrix_view gB = gsl_matrix_view_array(B.data(), Ndim, Ndim); gsl_linalg_LU_invert(gA, gP, &gB.matrix); gsl_permutation_free(gP); gsl_matrix_free(gA); @@ -81,9 +82,9 @@ R3::Matrix R3::inverse(const R3::Matrix& A) } -R3::Matrix R3::transpose(const R3::Matrix& A) +Matrix transpose(const Matrix& A) { - R3::Matrix res; + Matrix res; res = A(0,0), A(1,0), A(2,0), A(0,1), A(1,1), A(2,1), A(0,2), A(1,2), A(2,2); @@ -91,6 +92,22 @@ R3::Matrix R3::transpose(const R3::Matrix& A) } +const Matrix& product(const Matrix& A, const Matrix& B) +{ + static Matrix C; + C = 0.0; + for (int i = 0; i < Ndim; ++i) { + for (int j = 0; j < Ndim; ++j) { + double& cij = C(i, j); + for (int k = 0; k < Ndim; ++k) { + C(i, j) += A(i, k) * B(k, j); + } + } + } + return C; +} + +} // namespace R3 } // namespace srreal namespace mathutils { diff --git a/src/diffpy/srreal/R3linalg.hpp b/src/diffpy/srreal/R3linalg.hpp index 26604c51..94424cd8 100644 --- a/src/diffpy/srreal/R3linalg.hpp +++ b/src/diffpy/srreal/R3linalg.hpp @@ -22,8 +22,7 @@ #include #include -#include -#include +#include #include namespace diffpy { @@ -35,7 +34,6 @@ namespace R3 { // Constants const int Ndim = 3; -using blitz::product; using ::diffpy::mathutils::SQRT_DOUBLE_EPS; // Types @@ -50,6 +48,7 @@ const Matrix& zeros(); double determinant(const Matrix& A); Matrix inverse(const Matrix& A); Matrix transpose(const Matrix& A); +const Matrix& product(const Matrix&, const Matrix&); template double norm(const V&); template double distance(const V& u, const V& v); @@ -58,7 +57,6 @@ template Vector cross(const V& u, const V& v); template const Vector& mxvecproduct(const Matrix&, const V&); template const Vector& mxvecproduct(const V&, const Matrix&); - // Template functions -------------------------------------------------------- template diff --git a/src/diffpy/srreal/scatteringfactordata.cpp b/src/diffpy/srreal/scatteringfactordata.cpp index 478e1900..70862bfd 100644 --- a/src/diffpy/srreal/scatteringfactordata.cpp +++ b/src/diffpy/srreal/scatteringfactordata.cpp @@ -27,7 +27,7 @@ #include #include #include -#include +#include #include #include diff --git a/src/tests/TestLattice.hpp b/src/tests/TestLattice.hpp index 58a168f2..49e98249 100644 --- a/src/tests/TestLattice.hpp +++ b/src/tests/TestLattice.hpp @@ -31,70 +31,70 @@ class TestLattice : public CxxTest::TestSuite void setUp() { - lattice = new Lattice(); + lattice = new Lattice(); precision = 1.0e-8; allclose = EpsilonEqual(precision); } void tearDown() { - delete lattice; + delete lattice; } void test_Lattice() { - // default lattice should be cartesian - TS_ASSERT(allclose(lattice->base(), R3::identity())); - // lattice parameters constructor - Lattice lattice1(1.0, 2.0, 3.0, 90, 90, 120); - R3::Vector va = lattice1.va(); - R3::Vector vb = lattice1.vb(); - R3::Vector vc = lattice1.vc(); + // default lattice should be cartesian + TS_ASSERT(allclose(lattice->base(), R3::identity())); + // lattice parameters constructor + Lattice lattice1(1.0, 2.0, 3.0, 90, 90, 120); + R3::Vector va = lattice1.va(); + R3::Vector vb = lattice1.vb(); + R3::Vector vc = lattice1.vc(); double adotb = R3::dot(va, vb); double adotc = R3::dot(va, vc); double bdotc = R3::dot(vb, vc); - TS_ASSERT_DELTA(1.0, R3::norm(va), precision); - TS_ASSERT_DELTA(2.0, R3::norm(vb), precision); - TS_ASSERT_DELTA(3.0, R3::norm(vc), precision); - TS_ASSERT_DELTA(-0.5*1.0*2.0, adotb, precision); - TS_ASSERT_DELTA(0.0, adotc, precision); - TS_ASSERT_DELTA(0.0, bdotc, precision); - TS_ASSERT_DELTA(sqrt(4.0/3), lattice1.ar(), precision); - TS_ASSERT_DELTA(sqrt(1.0/3), lattice1.br(), precision); - TS_ASSERT_DELTA(1.0/3, lattice1.cr(), precision); - TS_ASSERT_DELTA(90.0, lattice1.alphar(), precision); - TS_ASSERT_DELTA(90.0, lattice1.betar(), precision); - TS_ASSERT_DELTA(60.0, lattice1.gammar(), precision); - // lattice vectors constructor - va = 1.0, 1.0, 0.0; - vb = 0.0, 1.0, 1.0; - vc = 1.0, 0.0, 1.0; - Lattice lattice2(va, vb, vc); - TS_ASSERT_DELTA(sqrt(2.0), lattice2.a(), precision); - TS_ASSERT_DELTA(60.0, lattice2.alpha(), precision); + TS_ASSERT_DELTA(1.0, R3::norm(va), precision); + TS_ASSERT_DELTA(2.0, R3::norm(vb), precision); + TS_ASSERT_DELTA(3.0, R3::norm(vc), precision); + TS_ASSERT_DELTA(-0.5*1.0*2.0, adotb, precision); + TS_ASSERT_DELTA(0.0, adotc, precision); + TS_ASSERT_DELTA(0.0, bdotc, precision); + TS_ASSERT_DELTA(sqrt(4.0/3), lattice1.ar(), precision); + TS_ASSERT_DELTA(sqrt(1.0/3), lattice1.br(), precision); + TS_ASSERT_DELTA(1.0/3, lattice1.cr(), precision); + TS_ASSERT_DELTA(90.0, lattice1.alphar(), precision); + TS_ASSERT_DELTA(90.0, lattice1.betar(), precision); + TS_ASSERT_DELTA(60.0, lattice1.gammar(), precision); + // lattice vectors constructor + va = 1.0, 1.0, 0.0; + vb = 0.0, 1.0, 1.0; + vc = 1.0, 0.0, 1.0; + Lattice lattice2(va, vb, vc); + TS_ASSERT_DELTA(sqrt(2.0), lattice2.a(), precision); + TS_ASSERT_DELTA(60.0, lattice2.alpha(), precision); } void test_setLatPar() { lattice->setLatPar(1.0, 2.0, 3.0, 90, 90, 120); - R3::Matrix base_check; - base_check = sqrt(0.75), -0.5, 0.0, - 0.0, +2.0, 0.0, - 0.0, +0.0, 3.0; - TS_ASSERT(allclose(base_check, lattice->base())); - R3::Matrix recbase_check; + R3::Matrix base_check; + base_check = sqrt(0.75), -0.5, 0.0, + 0.0, +2.0, 0.0, + 0.0, +0.0, 3.0; + TS_ASSERT(allclose(base_check, lattice->base())); + R3::Matrix recbase_check; recbase_check = sqrt(4.0/3), sqrt(1.0/12), 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 1.0/3.0; - TS_ASSERT(allclose(recbase_check, lattice->recbase())); + TS_ASSERT(allclose(recbase_check, lattice->recbase())); } void test_setLatBase() { R3::Vector va, vb, vc; va = 1.0, 1.0, 0.0; - vb = 0.0, 1.0, 1.0; - vc = 1.0, 0.0, 1.0; + vb = 0.0, 1.0, 1.0; + vc = 1.0, 0.0, 1.0; lattice->setLatBase(va, vb, vc); TS_ASSERT_DELTA(sqrt(2.0), lattice->a(), precision); TS_ASSERT_DELTA(sqrt(2.0), lattice->b(), precision); @@ -102,52 +102,52 @@ class TestLattice : public CxxTest::TestSuite TS_ASSERT_DELTA(60.0, lattice->alpha(), precision); TS_ASSERT_DELTA(60.0, lattice->beta(), precision); TS_ASSERT_DELTA(60.0, lattice->gamma(), precision); - // check determinant of rotation matrix - double detRot0 = R3::determinant(lattice->baserot()); + // check determinant of rotation matrix + double detRot0 = R3::determinant(lattice->baserot()); TS_ASSERT_DELTA(1.0, detRot0, precision); - // check if rotation matrix works - R3::Matrix base_check; - base_check = va[0], va[1], va[2], - vb[0], vb[1], vb[2], - vc[0], vc[1], vc[2]; - TS_ASSERT(allclose(base_check, lattice->base())); - R3::Matrix recbase_check; - recbase_check = 0.5, -0.5, 0.5, - 0.5, 0.5, -0.5, - -0.5, 0.5, 0.5; - TS_ASSERT(allclose(recbase_check, lattice->recbase())); + // check if rotation matrix works + R3::Matrix base_check; + base_check = va[0], va[1], va[2], + vb[0], vb[1], vb[2], + vc[0], vc[1], vc[2]; + TS_ASSERT(allclose(base_check, lattice->base())); + R3::Matrix recbase_check; + recbase_check = 0.5, -0.5, 0.5, + 0.5, 0.5, -0.5, + -0.5, 0.5, 0.5; + TS_ASSERT(allclose(recbase_check, lattice->recbase())); lattice->setLatPar( lattice->a(), lattice->b(), lattice->c(), - 44.0, 66.0, 88.0 ); - TS_ASSERT(!allclose(base_check, lattice->base())); - TS_ASSERT(!allclose(recbase_check, lattice->recbase())); + 44.0, 66.0, 88.0 ); + TS_ASSERT(!allclose(base_check, lattice->base())); + TS_ASSERT(!allclose(recbase_check, lattice->recbase())); lattice->setLatPar( lattice->a(), lattice->b(), lattice->c(), - 60.0, 60.0, 60.0 ); - TS_ASSERT(allclose(base_check, lattice->base())); - TS_ASSERT(allclose(recbase_check, lattice->recbase())); + 60.0, 60.0, 60.0 ); + TS_ASSERT(allclose(base_check, lattice->base())); + TS_ASSERT(allclose(recbase_check, lattice->recbase())); } void test_dist() { - R3::Vector va, vb; - va = 1.0, 2.0, 2.0; - vb = 0.0, 0.0, 0.0; - TS_ASSERT_DELTA(3.0, + R3::Vector va, vb; + va = 1.0, 2.0, 2.0; + vb = 0.0, 0.0, 0.0; + TS_ASSERT_DELTA(3.0, lattice->distance(va, vb), precision); - lattice->setLatPar(2.0, 2.0, 2.0, 90.0, 90.0, 90.0); - TS_ASSERT_DELTA(6.0, + lattice->setLatPar(2.0, 2.0, 2.0, 90.0, 90.0, 90.0); + TS_ASSERT_DELTA(6.0, lattice->distance(va, vb), precision); } void test_angle() { - R3::Vector va, vb; - va = 1.0, 0.0, 0.0; - vb = 0.0, 1.0, 0.0; - TS_ASSERT_DELTA(90.0, - lattice->angledeg(va, vb), precision); - lattice->setLatPar(2.0, 2.0, 2.0, 90.0, 90.0, 120.0); - TS_ASSERT_DELTA(120.0, - lattice->angledeg(va, vb), precision); + R3::Vector va, vb; + va = 1.0, 0.0, 0.0; + vb = 0.0, 1.0, 0.0; + TS_ASSERT_DELTA(90.0, + lattice->angledeg(va, vb), precision); + lattice->setLatPar(2.0, 2.0, 2.0, 90.0, 90.0, 120.0); + TS_ASSERT_DELTA(120.0, + lattice->angledeg(va, vb), precision); } void test_ucvCartesian() diff --git a/src/tests/TestPointsInSphere.hpp b/src/tests/TestPointsInSphere.hpp index 98260f71..397b6959 100644 --- a/src/tests/TestPointsInSphere.hpp +++ b/src/tests/TestPointsInSphere.hpp @@ -36,12 +36,12 @@ struct vidxgroup double vijk[4]; vidxgroup(double v, const int* ijk) { - vijk[0] = v; - for (size_t i = 0; i != 3; ++i) { vijk[i+1] = ijk[i]; } + vijk[0] = v; + for (size_t i = 0; i != 3; ++i) { vijk[i+1] = ijk[i]; } } vidxgroup(double v, int i, int j, int k) { - vijk[0] = v; vijk[1] = i; vijk[2] = j; vijk[3] = k; + vijk[0] = v; vijk[1] = i; vijk[2] = j; vijk[3] = k; } }; @@ -55,14 +55,14 @@ bool operator<(const vidxgroup &x, const vidxgroup &y) bool operator==(const vidxgroup &x, const vidxgroup &y) { bool eq = (fabs(x.vijk[0] - y.vijk[0]) < eps) && - equal(x.vijk+1, x.vijk+4, y.vijk+1); + equal(x.vijk+1, x.vijk+4, y.vijk+1); return eq; } ostream& operator<<(ostream &s, const vidxgroup &x) { return s << "<" << x.vijk[0] << ";" << int(x.vijk[1]) - << ',' << int(x.vijk[2]) << ',' << int(x.vijk[3]) << '>'; + << ',' << int(x.vijk[2]) << ',' << int(x.vijk[3]) << '>'; } } // namespace @@ -84,134 +84,134 @@ class TestPointsInSphere : public CxxTest::TestSuite void setUp() { - latpar = new LatticeParameters(1, 1, 1, 90, 90, 90); + latpar = new LatticeParameters(1, 1, 1, 90, 90, 90); } void tearDown() { - delete latpar; + delete latpar; } private: int count(double Rmin, double Rmax) { - int c = 0; - PointsInSphere sph(Rmin, Rmax, *latpar); - for (; !sph.finished(); sph.next()) ++c; - return c; + int c = 0; + PointsInSphere sph(Rmin, Rmax, *latpar); + for (; !sph.finished(); sph.next()) ++c; + return c; } vector sortedPoints(double Rmin, double Rmax) { - vector ridx; - for ( PointsInSphere sph(Rmin, Rmax, *latpar); - not sph.finished(); sph.next() ) - { - ridx.push_back(vidxgroup(sph.r(), sph.mno())); - } - sort(ridx.begin(), ridx.end()); - return ridx; + vector ridx; + for ( PointsInSphere sph(Rmin, Rmax, *latpar); + not sph.finished(); sph.next() ) + { + ridx.push_back(vidxgroup(sph.r(), sph.mno())); + } + sort(ridx.begin(), ridx.end()); + return ridx; } public: void test_Cubic() { - latpar->a = latpar->b = latpar->c = 1.0; - latpar->alpha = latpar->beta = latpar->gamma = 90.0; - latpar->update(); - TS_ASSERT_EQUALS(0, count(0.0, 0.0)); - TS_ASSERT_EQUALS(0, count(eps, 0.5)); - TS_ASSERT_EQUALS(0, count(1.0 + eps, 1.1)); - TS_ASSERT_EQUALS(1, count(0.0, eps)); - TS_ASSERT_EQUALS(7, count(0.0, 1 + eps)); - TS_ASSERT_EQUALS(19, count(0.0, sqrt(2.0) + eps)); - TS_ASSERT_EQUALS(12, count(1.0 + eps, sqrt(2.0) + eps)); + latpar->a = latpar->b = latpar->c = 1.0; + latpar->alpha = latpar->beta = latpar->gamma = 90.0; + latpar->update(); + TS_ASSERT_EQUALS(0, count(0.0, 0.0)); + TS_ASSERT_EQUALS(0, count(eps, 0.5)); + TS_ASSERT_EQUALS(0, count(1.0 + eps, 1.1)); + TS_ASSERT_EQUALS(1, count(0.0, eps)); + TS_ASSERT_EQUALS(7, count(0.0, 1 + eps)); + TS_ASSERT_EQUALS(19, count(0.0, sqrt(2.0) + eps)); + TS_ASSERT_EQUALS(12, count(1.0 + eps, sqrt(2.0) + eps)); } void test_Orthorombic() { - latpar->a = 1.0; latpar->b = 2.0; latpar->c = 3.0; - latpar->alpha = latpar->beta = latpar->gamma = 90.0; - latpar->update(); - TS_ASSERT_EQUALS(3, count(0.0, 1.1)); - TS_ASSERT_EQUALS(4, count(1.9, 2.1)); - vidxgroup ep[] = { - vidxgroup(0, 0, 0, 0), - vidxgroup(1, -1, 0, 0), - vidxgroup(1, 1, 0, 0), - vidxgroup(2, -2, 0, 0), - vidxgroup(2, 0, -1, 0), - vidxgroup(2, 0, 1, 0), - vidxgroup(2, 2, 0, 0), - vidxgroup(sqrt(5.0), -1, -1, 0), - vidxgroup(sqrt(5.0), -1, 1, 0), - vidxgroup(sqrt(5.0), 1, -1, 0), - vidxgroup(sqrt(5.0), 1, 1, 0), - vidxgroup(sqrt(8.0), -2, -1, 0), - vidxgroup(sqrt(8.0), -2, 1, 0), - vidxgroup(sqrt(8.0), 2, -1, 0), - vidxgroup(sqrt(8.0), 2, 1, 0), - vidxgroup(3, -3, 0, 0), - vidxgroup(3, 0, 0, -1), - vidxgroup(3, 0, 0, 1), - vidxgroup(3, 3, 0, 0), - }; - vector exp_pts(ep, ep + sizeof(ep)/sizeof(vidxgroup)); - vector act_pts = sortedPoints(0.0, 3.0+eps); - TS_ASSERT_EQUALS(exp_pts.size(), act_pts.size()); - for (size_t i = 0; i != exp_pts.size(); ++i) - { - TS_ASSERT_EQUALS(exp_pts[i], act_pts[i]); - } + latpar->a = 1.0; latpar->b = 2.0; latpar->c = 3.0; + latpar->alpha = latpar->beta = latpar->gamma = 90.0; + latpar->update(); + TS_ASSERT_EQUALS(3, count(0.0, 1.1)); + TS_ASSERT_EQUALS(4, count(1.9, 2.1)); + vidxgroup ep[] = { + vidxgroup(0, 0, 0, 0), + vidxgroup(1, -1, 0, 0), + vidxgroup(1, 1, 0, 0), + vidxgroup(2, -2, 0, 0), + vidxgroup(2, 0, -1, 0), + vidxgroup(2, 0, 1, 0), + vidxgroup(2, 2, 0, 0), + vidxgroup(sqrt(5.0), -1, -1, 0), + vidxgroup(sqrt(5.0), -1, 1, 0), + vidxgroup(sqrt(5.0), 1, -1, 0), + vidxgroup(sqrt(5.0), 1, 1, 0), + vidxgroup(sqrt(8.0), -2, -1, 0), + vidxgroup(sqrt(8.0), -2, 1, 0), + vidxgroup(sqrt(8.0), 2, -1, 0), + vidxgroup(sqrt(8.0), 2, 1, 0), + vidxgroup(3, -3, 0, 0), + vidxgroup(3, 0, 0, -1), + vidxgroup(3, 0, 0, 1), + vidxgroup(3, 3, 0, 0), + }; + vector exp_pts(ep, ep + sizeof(ep)/sizeof(vidxgroup)); + vector act_pts = sortedPoints(0.0, 3.0+eps); + TS_ASSERT_EQUALS(exp_pts.size(), act_pts.size()); + for (size_t i = 0; i != exp_pts.size(); ++i) + { + TS_ASSERT_EQUALS(exp_pts[i], act_pts[i]); + } } void test_Hexagonal() { - latpar->a = 1.0; latpar->b = 1.0; latpar->c = 2.0; - latpar->alpha = latpar->beta = 90.0; latpar->gamma = 120.0; - latpar->update(); - TS_ASSERT_EQUALS(7, count(0.0, 1+eps)); - vidxgroup ep[] = { - vidxgroup(0, 0, 0, 0), - vidxgroup(1, -1, -1, 0), - vidxgroup(1, -1, 0, 0), - vidxgroup(1, 0, -1, 0), - vidxgroup(1, 0, 1, 0), - vidxgroup(1, 1, 0, 0), - vidxgroup(1, 1, 1, 0), - vidxgroup(sqrt(3.0), -2, -1, 0), - vidxgroup(sqrt(3.0), -1, -2, 0), - vidxgroup(sqrt(3.0), -1, 1, 0), - vidxgroup(sqrt(3.0), 1, -1, 0), - vidxgroup(sqrt(3.0), 1, 2, 0), - vidxgroup(sqrt(3.0), 2, 1, 0), - vidxgroup(2, -2, -2, 0), - vidxgroup(2, -2, 0, 0), - vidxgroup(2, 0, -2, 0), - vidxgroup(2, 0, 0, -1), - vidxgroup(2, 0, 0, 1), - vidxgroup(2, 0, 2, 0), - vidxgroup(2, 2, 0, 0), - vidxgroup(2, 2, 2, 0), - }; - vector exp_pts(ep, ep + sizeof(ep)/sizeof(vidxgroup)); - vector act_pts = sortedPoints(0.0, 2.0+eps); - TS_ASSERT_EQUALS(exp_pts.size(), act_pts.size()); - for (size_t i = 0; i != exp_pts.size(); ++i) - { - TS_ASSERT_EQUALS(exp_pts[i], act_pts[i]); - } + latpar->a = 1.0; latpar->b = 1.0; latpar->c = 2.0; + latpar->alpha = latpar->beta = 90.0; latpar->gamma = 120.0; + latpar->update(); + TS_ASSERT_EQUALS(7, count(0.0, 1+eps)); + vidxgroup ep[] = { + vidxgroup(0, 0, 0, 0), + vidxgroup(1, -1, -1, 0), + vidxgroup(1, -1, 0, 0), + vidxgroup(1, 0, -1, 0), + vidxgroup(1, 0, 1, 0), + vidxgroup(1, 1, 0, 0), + vidxgroup(1, 1, 1, 0), + vidxgroup(sqrt(3.0), -2, -1, 0), + vidxgroup(sqrt(3.0), -1, -2, 0), + vidxgroup(sqrt(3.0), -1, 1, 0), + vidxgroup(sqrt(3.0), 1, -1, 0), + vidxgroup(sqrt(3.0), 1, 2, 0), + vidxgroup(sqrt(3.0), 2, 1, 0), + vidxgroup(2, -2, -2, 0), + vidxgroup(2, -2, 0, 0), + vidxgroup(2, 0, -2, 0), + vidxgroup(2, 0, 0, -1), + vidxgroup(2, 0, 0, 1), + vidxgroup(2, 0, 2, 0), + vidxgroup(2, 2, 0, 0), + vidxgroup(2, 2, 2, 0), + }; + vector exp_pts(ep, ep + sizeof(ep)/sizeof(vidxgroup)); + vector act_pts = sortedPoints(0.0, 2.0+eps); + TS_ASSERT_EQUALS(exp_pts.size(), act_pts.size()); + for (size_t i = 0; i != exp_pts.size(); ++i) + { + TS_ASSERT_EQUALS(exp_pts[i], act_pts[i]); + } } void test_FCC() { - latpar->a = latpar->b = latpar->c = sqrt(0.5); - latpar->alpha = latpar->beta = latpar->gamma = 60.0; - latpar->update(); - TS_ASSERT_EQUALS(13, count(0.0, sqrt(0.5)+eps)); - TS_ASSERT_EQUALS(19, count(0.0, 1.0+eps)); + latpar->a = latpar->b = latpar->c = sqrt(0.5); + latpar->alpha = latpar->beta = latpar->gamma = 60.0; + latpar->update(); + TS_ASSERT_EQUALS(13, count(0.0, sqrt(0.5)+eps)); + TS_ASSERT_EQUALS(19, count(0.0, 1.0+eps)); } }; // class TestPointsInSphere diff --git a/src/tests/TestR3linalg.hpp b/src/tests/TestR3linalg.hpp index 2c34aabe..53950158 100644 --- a/src/tests/TestR3linalg.hpp +++ b/src/tests/TestR3linalg.hpp @@ -44,47 +44,47 @@ class TestR3linalg : public CxxTest::TestSuite void test_determinant() { - // default lattice should be cartesian - R3::Matrix A1, A2; - A1 = 9, 1, 9, - 6, 7, 4, - 0, 2, 9; - A2 = 9, 1, 9, - 0, 2, 9, - 6, 7, 4; - double detA = 549; - TS_ASSERT_EQUALS(detA, R3::determinant(A1)); - TS_ASSERT_EQUALS(-detA, R3::determinant(A2)); + // default lattice should be cartesian + R3::Matrix A1, A2; + A1 = 9, 1, 9, + 6, 7, 4, + 0, 2, 9; + A2 = 9, 1, 9, + 0, 2, 9, + 6, 7, 4; + double detA = 549; + TS_ASSERT_EQUALS(detA, R3::determinant(A1)); + TS_ASSERT_EQUALS(-detA, R3::determinant(A2)); } void test_inverse() { - R3::Matrix A, invA; - A = - 0.5359, -0.5904, 0.8670, - -0.0053, 0.7559, 0.2692, - -0.8926, 0.9424, 0.9692; - invA = - 0.49063197005867, 1.42323870111089, -0.83420736316541, - -0.24089965988852, 1.32489393619466, -0.15249839300481, - 0.68609361943181, 0.02249568627913, 0.41178393851247; - TS_ASSERT(allclose(invA, R3::inverse(A))); + R3::Matrix A, invA; + A = + 0.5359, -0.5904, 0.8670, + -0.0053, 0.7559, 0.2692, + -0.8926, 0.9424, 0.9692; + invA = + 0.49063197005867, 1.42323870111089, -0.83420736316541, + -0.24089965988852, 1.32489393619466, -0.15249839300481, + 0.68609361943181, 0.02249568627913, 0.41178393851247; + TS_ASSERT(allclose(invA, R3::inverse(A))); } void test_transpose() { - R3::Matrix A, Atrans; - A = - 0.5359, -0.5904, 0.8670, - -0.0053, 0.7559, 0.2692, - -0.8926, 0.9424, 0.9692; - Atrans = - 0.5359, -0.0053, -0.8926, + R3::Matrix A, Atrans; + A = + 0.5359, -0.5904, 0.8670, + -0.0053, 0.7559, 0.2692, + -0.8926, 0.9424, 0.9692; + Atrans = + 0.5359, -0.0053, -0.8926, -0.5904, 0.7559, 0.9424, 0.8670, 0.2692, 0.9692; - TS_ASSERT_EQUALS(Atrans, R3::transpose(A)); + TS_ASSERT_EQUALS(Atrans, R3::transpose(A)); } diff --git a/src/tests/testdata/PbScW25TiO3.stru b/src/tests/testdata/PbScW25TiO3.stru index 433147b3..71ba3651 100644 --- a/src/tests/testdata/PbScW25TiO3.stru +++ b/src/tests/testdata/PbScW25TiO3.stru @@ -2,11 +2,11 @@ title Pb8 (Sc16/3 W8/3)0.75 (Ti8)0.25 O24, B factors from BNL2 format pdffit scale 1.000000 sharp 0.000000, 1.000000, 0.000000 -spcgr Fm-3m +spcgr Fm-3m cell 8.072700, 8.072700, 8.072700, 90.000000, 90.000000, 90.000000 dcell 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000 ncell 1, 1, 1, 56 -atoms +atoms PB 0.25000000 0.25000000 0.25000000 1.0000 0.00000000 0.00000000 0.00000000 0.0000 0.04795025 0.04795025 0.04795025