Mantid
Loading...
Searching...
No Matches
OrientedLattice.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
9
10namespace Mantid::Geometry {
13
14namespace {
15const double TWO_PI = 2. * M_PI;
16}
17
22 if (Umatrix.isRotation()) {
23 U = Umatrix;
24 UB = U * getB();
25 ModUB = UB * getModHKL();
26 } else
27 throw std::invalid_argument("U is not a proper rotation");
28}
29
39OrientedLattice::OrientedLattice(const double _a, const double _b, const double _c, const DblMatrix &Umatrix)
40 : UnitCell(_a, _b, _c) {
41 if (Umatrix.isRotation()) {
42 U = Umatrix;
43 UB = U * getB();
44 ModUB = UB * getModHKL();
45 } else
46 throw std::invalid_argument("U is not a proper rotation");
47}
48
59OrientedLattice::OrientedLattice(const double _a, const double _b, const double _c, const double _alpha,
60 const double _beta, const double _gamma, const DblMatrix &Umatrix, const int angleunit)
61 : UnitCell(_a, _b, _c, _alpha, _beta, _gamma, angleunit) {
62 if (Umatrix.isRotation()) {
63 U = Umatrix;
64 UB = U * getB();
65 ModUB = UB * getModHKL();
66 } else
67 throw std::invalid_argument("U is not a proper rotation");
68}
69
74OrientedLattice::OrientedLattice(const UnitCell &uc, const DblMatrix &Umatrix) : UnitCell(uc), U(Umatrix) {
75 if (Umatrix.isRotation()) {
76 U = Umatrix;
77 UB = U * getB();
78 ModUB = UB * getModHKL();
79 } else
80 throw std::invalid_argument("U is not a proper rotation");
81}
82
85const DblMatrix &OrientedLattice::getU() const { return U; }
86
95const DblMatrix &OrientedLattice::getUB() const { return UB; }
96
97const DblMatrix &OrientedLattice::getModUB() const { return ModUB; }
98
103void OrientedLattice::setU(const DblMatrix &newU, const bool force) {
104 // determinant ==1 or (determinant == +/-1 and force)
105 if (newU.isRotation() || (force && newU.isOrthogonal())) {
106 U = newU;
107 UB = U * getB();
108 ModUB = UB * getModHKL();
109 } else
110 throw std::invalid_argument("U is not a proper rotation");
111}
112
116 // check if determinant is close to 0. The 1e-10 value is arbitrary
117 if (std::fabs(newUB.determinant()) > 1e-10) {
118 UB = newUB;
119 DblMatrix newGstar, B;
120 newGstar = newUB.Tprime() * newUB;
121 this->recalculateFromGstar(newGstar);
122 B = this->getB();
123 B.Invert();
124 U = newUB * B;
125 } else
126 throw std::invalid_argument("determinant of UB is too close to 0");
127}
128
133 ModUB = newModUB;
134 DblMatrix UBinv, newModHKL;
135 UBinv = getUB();
136 UBinv.Invert();
137 newModHKL = UBinv * ModUB;
138 setModHKL(newModHKL);
139}
140
146 DblMatrix UBinv = this->getUB();
147 UBinv.Invert();
148 V3D out = UBinv * Q / TWO_PI; // transform back to HKL
149 return out;
150}
151
156V3D OrientedLattice::qFromHKL(const V3D &hkl) const { return UB * hkl * TWO_PI; }
157
163 V3D z(0, 0, 1);
164 DblMatrix UBinv = UB;
165 UBinv.Invert();
166 return UBinv * z;
167}
168
176 V3D x(1, 0, 0);
177 DblMatrix UBinv = UB;
178 UBinv.Invert();
179 return UBinv * x;
180}
181
196 const DblMatrix &BMatrix = this->getB();
197 V3D buVec = BMatrix * u;
198 V3D bvVec = BMatrix * v;
199 // try to make an orthonormal system
200 if (buVec.norm2() < 1e-10)
201 throw std::invalid_argument("|B.u|~0");
202 if (bvVec.norm2() < 1e-10)
203 throw std::invalid_argument("|B.v|~0");
204 buVec.normalize(); // 1st unit vector, along Bu
205 V3D bwVec = buVec.cross_prod(bvVec);
206 const auto norm = bwVec.norm();
207 if (norm < 1e-5) {
208 // 3rd unit vector, perpendicular to Bu,Bv
209 throw std::invalid_argument("u and v are parallel");
210 }
211 bwVec /= norm;
212 // 2nd unit vector, perpendicular to Bu, in the Bu,Bv plane
213 bvVec = bwVec.cross_prod(buVec);
214 DblMatrix tau(3, 3), lab(3, 3), U(3, 3);
215 /*lab = U tau
216 / 0 1 0 \ /bu[0] bv[0] bw[0]\
217 | 0 0 1 | = U |bu[1] bv[1] bw[1]|
218 \ 1 0 0 / \bu[2] bv[2] bw[2]/
219 */
220 lab[0][1] = 1.;
221 lab[1][2] = 1.;
222 lab[2][0] = 1.;
223 tau[0][0] = buVec[0];
224 tau[0][1] = bvVec[0];
225 tau[0][2] = bwVec[0];
226 tau[1][0] = buVec[1];
227 tau[1][1] = bvVec[1];
228 tau[1][2] = bwVec[1];
229 tau[2][0] = buVec[2];
230 tau[2][1] = bvVec[2];
231 tau[2][2] = bwVec[2];
232 tau.Invert();
233 U = lab * tau;
234 this->setU(U);
235 return getU();
236}
237
242void OrientedLattice::saveNexus(::NeXus::File *file, const std::string &group) const {
243 file->makeGroup(group, "NXcrystal", true);
244 file->writeData("unit_cell_a", this->a());
245 file->writeData("unit_cell_b", this->b());
246 file->writeData("unit_cell_c", this->c());
247 file->writeData("unit_cell_alpha", this->alpha());
248 file->writeData("unit_cell_beta", this->beta());
249 file->writeData("unit_cell_gamma", this->gamma());
250 // Save the UB matrix
251 std::vector<double> ub = this->UB.getVector();
252 std::vector<int> dims(2, 3); // 3x3 matrix
253 file->writeData("orientation_matrix", ub, dims);
254
255 file->closeGroup();
256}
257
262void OrientedLattice::loadNexus(::NeXus::File *file, const std::string &group) {
263 file->openGroup(group, "NXcrystal");
264 std::vector<double> ub;
265 file->readData("orientation_matrix", ub);
266 // Turn into a matrix
267 DblMatrix ubMat(ub);
268 // This will set the lattice parameters and the U matrix:
269 this->setUB(ubMat);
270 file->closeGroup();
271}
272
285bool OrientedLattice::GetUB(DblMatrix &UB, const V3D &a_dir, const V3D &b_dir, const V3D &c_dir) {
286 if (UB.numRows() != 3 || UB.numCols() != 3) {
287 throw std::invalid_argument("Find_UB(): UB matrix NULL or not 3X3");
288 }
289
290 UB.setRow(0, a_dir);
291 UB.setRow(1, b_dir);
292 UB.setRow(2, c_dir);
293 try {
294 UB.Invert();
295 } catch (...) {
296 return false;
297 }
298 return true;
299}
300
318bool OrientedLattice::GetABC(const DblMatrix &UB, V3D &a_dir, V3D &b_dir, V3D &c_dir) {
319 if (UB.numRows() != 3 || UB.numCols() != 3) {
320 throw std::invalid_argument("GetABC(): UB matrix NULL or not 3X3");
321 }
322
323 DblMatrix UB_inverse(UB);
324 try {
325 UB_inverse.Invert();
326 } catch (...) {
327 return false;
328 }
329 a_dir(UB_inverse[0][0], UB_inverse[0][1], UB_inverse[0][2]);
330 b_dir(UB_inverse[1][0], UB_inverse[1][1], UB_inverse[1][2]);
331 c_dir(UB_inverse[2][0], UB_inverse[2][1], UB_inverse[2][2]);
332
333 return true;
334}
339 UB = U * getB();
340}
341bool OrientedLattice::operator==(const OrientedLattice &other) const { return UB == other.UB; }
342bool OrientedLattice::operator!=(const OrientedLattice &other) const { return UB != other.UB; }
343} // namespace Mantid::Geometry
Class to implement UB matrix.
Kernel::V3D getuVector() const
gets a vector along beam direction when goniometers are at 0.
void setU(const Kernel::DblMatrix &newU, const bool force=true)
Sets the U matrix.
Kernel::V3D getvVector() const
gets a vector in the horizontal plane, perpendicular to the beam direction when goniometers are at 0.
OrientedLattice(const Kernel::DblMatrix &Umatrix=Kernel::DblMatrix(3, 3, true))
Default constructor.
void setModUB(const Kernel::DblMatrix &newModUB)
Sets the Modulation UB matrix.
Kernel::V3D hklFromQ(const Kernel::V3D &Q) const
Return hkl from the Q-sample coordinates.
void setUB(const Kernel::DblMatrix &newUB)
Sets the UB matrix and recalculates lattice parameters.
static bool GetABC(const Kernel::DblMatrix &UB, Kernel::V3D &a_dir, Kernel::V3D &b_dir, Kernel::V3D &c_dir)
Get the real space edge vectors a, b, c corresponding to the UB matrix.
bool operator==(const OrientedLattice &other) const
const Kernel::DblMatrix & setUFromVectors(const Kernel::V3D &u, const Kernel::V3D &v)
Create the U matrix from two vectors.
void recalculate() override
Private function, called at initialization or whenever lattice parameters are changed.
const Kernel::DblMatrix & getUB() const
Get the UB matrix.
void loadNexus(::NeXus::File *file, const std::string &group)
Load the lattice to from an open NeXus file.
void saveNexus(::NeXus::File *file, const std::string &group) const
Save the lattice to an open NeXus file.
bool operator!=(const OrientedLattice &other) const
static bool GetUB(Kernel::DblMatrix &UB, const Kernel::V3D &a_dir, const Kernel::V3D &b_dir, const Kernel::V3D &c_dir)
Get the UB matix corresponding to the real space edge vectors a, b, c.
Kernel::V3D qFromHKL(const Kernel::V3D &hkl) const
Return Q-sample coordinates from hkl.
const Kernel::DblMatrix & getModUB() const
const Kernel::DblMatrix & getU() const
Get the U matrix.
void recalculateFromGstar(const Kernel::DblMatrix &NewGstar) override
Make recalculateFrom private.
Class to implement unit cell of crystals.
Definition: UnitCell.h:44
const Kernel::DblMatrix & getB() const
Get the B-matrix.
Definition: UnitCell.cpp:758
double alpha() const
Get lattice parameter.
Definition: UnitCell.cpp:133
double c() const
Get lattice parameter.
Definition: UnitCell.cpp:128
const Kernel::DblMatrix & getModHKL() const
Get modulation vectors for satellites.
Definition: UnitCell.cpp:548
void setModHKL(double _dh1, double _dk1, double _dl1, double _dh2, double _dk2, double _dl2, double _dh3, double _dk3, double _dl3)
Set modulation vectors for satellites.
Definition: UnitCell.cpp:353
double beta() const
Get lattice parameter.
Definition: UnitCell.cpp:138
double a() const
Get lattice parameter.
Definition: UnitCell.cpp:118
virtual void recalculate()
Private function, called at initialization or whenever lattice parameters are changed.
Definition: UnitCell.cpp:766
double b() const
Get lattice parameter.
Definition: UnitCell.cpp:123
Kernel::DblMatrix B
B matrix for a right-handed coordinate system, in Busing-Levy convention.
Definition: UnitCell.h:191
double gamma() const
Get lattice parameter.
Definition: UnitCell.cpp:143
T determinant() const
Calculate the determinant.
Definition: Matrix.cpp:1048
T Invert()
LU inversion routine.
Definition: Matrix.cpp:924
Matrix< T > Tprime() const
Transpose the matrix.
Definition: Matrix.cpp:766
bool isRotation() const
Check if a matrix represents a proper rotation @ return :: true/false.
Definition: Matrix.cpp:1424
void setRow(const size_t nRow, const std::vector< T > &newRow)
Definition: Matrix.cpp:686
std::vector< T > getVector() const
Definition: Matrix.cpp:77
size_t numRows() const
Return the number of rows in the matrix.
Definition: Matrix.h:144
size_t numCols() const
Return the number of columns in the matrix.
Definition: Matrix.h:147
bool isOrthogonal() const
Check if a matrix is orthogonal.
Definition: Matrix.cpp:1444
Class for 3D vectors.
Definition: V3D.h:34
double normalize()
Make a normalized vector (return norm value)
Definition: V3D.cpp:130
constexpr V3D cross_prod(const V3D &v) const noexcept
Cross product (this * argument)
Definition: V3D.h:278
double norm() const noexcept
Definition: V3D.h:263
constexpr double norm2() const noexcept
Vector length squared.
Definition: V3D.h:265