11#include <boost/version.hpp>
14#if BOOST_VERSION < 106700
15#include <boost/math/common_factor.hpp>
16using boost::math::gcd;
18#include <boost/integer/common_factor.hpp>
19using boost::integer::gcd;
22#include <nexus/NeXusFile.hpp>
32double nearInt(
double val,
double eps,
double mult)
noexcept {
37 if (std::abs(val - std::round(val)) > eps) {
38 mult *= std::ceil(val / eps) * eps / val;
57void V3D::spherical(
const double R,
const double theta,
const double phi)
noexcept {
58 constexpr double deg2rad = M_PI / 180.0;
71 m_pt[2] = R * cos(polar);
72 const double ct = R * sin(polar);
73 m_pt[0] = ct * cos(azimuth);
74 m_pt[1] = ct * sin(azimuth);
96 m_pt[1] = R * cos(polar);
97 const double ct = R * sin(polar);
98 m_pt[0] = ct * cos(azimuth);
99 m_pt[2] = ct * sin(azimuth);
118 constexpr double rad2deg = 180.0 / M_PI;
122 theta = acos(m_pt[2] / R) *
rad2deg;
123 phi = atan2(m_pt[1], m_pt[0]) *
rad2deg;
131 const double ND(
norm());
133 throw std::runtime_error(
"Unable to normalize a zero length vector.");
141 for (
auto &p :
m_pt) {
151 const double R = distance(v);
153 const double zOffset = m_pt[2] - v.m_pt[2];
154 return acos(zOffset / R);
166 const double ratio = this->
cosAngle(v);
169 else if (ratio <= -1.0)
181 const double n1 =
norm();
182 const double n2 = v.
norm();
183 if (n1 == 0. || n2 == 0.) {
184 throw std::runtime_error(
"Cannot calculate an angle when one of the vectors has zero length.");
201 for (
int i = 0; i < 3; i++) {
206 const double det = T.
Invert();
207 if (std::abs(det) < 1e-13)
218 const double xold(m_pt[0]), yold(m_pt[1]), zold(m_pt[2]);
219 for (
size_t i = 0; i < m_pt.size(); ++i) {
220 m_pt[i] = A[i][0] * xold + A[i][1] * yold + A[i][2] * zold;
231 const V3D &Av = *
this;
232 const V3D Tmp((Bv - Av).cross_prod(Cv - Av));
242 for (
const double p : m_pt) {
252 const auto l = norm();
266 double max = m_pt[0] * m_pt[0];
268 double u2 = m_pt[1] * m_pt[1];
269 int idx = (m_pt[0] > 0) ? 1 : -1;
272 idx = (m_pt[1] > 0) ? 2 : -2;
275 u2 = m_pt[2] * m_pt[2];
278 idx = (m_pt[2] > 0) ? 3 : -3;
299 if (vectors.size() != 2)
300 throw std::invalid_argument(
"makeVectorsOrthogonal() only works with 2 vectors");
305 std::vector<V3D> out;
307 out.emplace_back(v0);
315 out.emplace_back(v1);
319 out.emplace_back(v2);
336 OX << m_pt[0] <<
" " << m_pt[1] <<
" " << m_pt[2];
341 std::ostringstream mess;
349 std::istringstream mess(str);
366 std::getline(IX, in);
367 size_t i = in.find_first_of(
'[');
368 if (i == std::string::npos)
369 throw std::runtime_error(
"Wrong format for V3D input: " + in);
370 size_t j = in.find_last_of(
']');
371 if (j == std::string::npos || j < i + 6)
372 throw std::runtime_error(
"Wrong format for V3D input: " + in);
374 size_t c1 = in.find_first_of(
',');
375 size_t c2 = in.find_first_of(
',', c1 + 1);
376 if (c1 == std::string::npos || c2 == std::string::npos)
377 throw std::runtime_error(
"Wrong format for V3D input: [" + in +
"]");
379 m_pt[0] = std::stod(in.substr(i + 1, c1 - i - 1));
380 m_pt[1] = std::stod(in.substr(c1 + 1, c2 - c1 - 1));
381 m_pt[2] = std::stod(in.substr(c2 + 1, j - c2 - 1));
413 file->makeData(name, ::NeXus::FLOAT64, 3,
true);
414 file->putData(
m_pt.data());
424 std::vector<double> data;
425 file->readData(name, data);
426 if (data.size() != 3)
427 throw std::runtime_error(
"Unexpected data size when reading a V3D NXS field '" + name +
"'. Expected 3.");
428 std::copy(data.cbegin(), data.cend(),
m_pt.begin());
434 if (eps < FLT_EPSILON)
439 double ax = std::abs(
m_pt[0]);
440 double ay = std::abs(
m_pt[1]);
441 double az = std::abs(
m_pt[2]);
443 double amax = (ax > ay) ? ax : ay;
444 amax = (az > amax) ? az : amax;
445 if (amax < FLT_EPSILON)
446 throw(std::invalid_argument(
"vector length is less then accuracy requested"));
462 mult = nearInt(ax, eps, mult);
463 mult = nearInt(ay, eps, mult);
464 mult = nearInt(az, eps, mult);
466 const size_t iax = std::lround(ax * mult / eps);
467 const size_t iay = std::lround(ay * mult / eps);
468 const size_t iaz = std::lround(az * mult / eps);
470 const size_t div = gcd(iax, gcd(iay, iaz));
471 mult /= (
static_cast<double>(div) * eps);
484 const double mag_sq_1 = v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2];
485 const double mag_sq_2 = v2[0] * v2[0] + v2[1] * v2[1] + v2[2] * v2[2];
486 return (mag_sq_1 < mag_sq_2);
496 const double conversionFactor = inDegrees ? 180. / M_PI : 1.0;
497 const double divisor = this->
norm();
499 throw std::runtime_error(
"Cannot calculate direction angles for zero length vector");
501 return V3D(conversionFactor * acos(
m_pt[0] / divisor), conversionFactor * acos(
m_pt[1] / divisor),
502 conversionFactor * acos(
m_pt[2] / divisor));
511 if (abs(
static_cast<int>(
m_pt[0])) > MaxOrder)
512 MaxOrder = abs(
static_cast<int>(
m_pt[0]));
513 if (abs(
static_cast<int>(
m_pt[1])) > MaxOrder)
514 MaxOrder = abs(
static_cast<int>(
m_pt[1]));
515 if (abs(
static_cast<int>(
m_pt[2])) > MaxOrder)
516 MaxOrder = abs(
static_cast<int>(
m_pt[2]));
T Invert()
LU inversion routine.
void rotate(V3D &) const
Rotate a vector.
void setRotation(const double deg)
Set the rotation (both don't change rotation axis)
void loadNexus(::NeXus::File *file, const std::string &name)
Load the object from an open NeXus file.
void spherical_rad(const double R, const double polar, const double azimuth) noexcept
Sets the vector position based on spherical coordinates, in radians.
std::array< double, 3 > m_pt
static bool compareMagnitude(const Kernel::V3D &v1, const Kernel::V3D &v2)
Convenience method for sorting list of V3D objects based on magnitude.
bool coLinear(const V3D &, const V3D &) const noexcept
Determines if this,B,C are colinear.
static std::vector< V3D > makeVectorsOrthogonal(const std::vector< V3D > &vectors)
Take a list of 2 vectors and makes a 3D orthogonal system out of them The first vector i0 is taken as...
void azimuth_polar_SNS(const double R, const double azimuth, const double polar) noexcept
Sets the vector position based on azimuth and polar angle, in RADIANS, in the SNS instrument coordina...
V3D absoluteValue() const
Absolute value.
double zenith(const V3D &) const noexcept
Zenith (theta) angle between this and another vector.
constexpr double scalar_prod(const V3D &v) const noexcept
Calculates the cross product.
void printSelf(std::ostream &) const
Prints a text representation of itself in format "[x,y,z]".
V3D & operator*=(const V3D &v) noexcept
Self-Inner product.
V3D & operator/=(const V3D &v) noexcept
Self-Inner division.
double normalize()
Make a normalized vector (return norm value)
constexpr V3D cross_prod(const V3D &v) const noexcept
Cross product (this * argument)
std::string toString() const
int maxCoeff()
Maximum absolute integer value.
V3D directionAngles(bool inDegrees=true) const
Direction angles.
void spherical(const double R, const double theta, const double phi) noexcept
Sets the vector position based on spherical coordinates.
void read(std::istream &)
Read data from a stream.
double angle(const V3D &) const
Angle between this and another vector.
void rotate(const Matrix< double > &) noexcept
Rotate a point by a matrix.
double norm() const noexcept
int reBase(const V3D &, const V3D &, const V3D &) noexcept
rebase to new basis vector
void fromString(const std::string &str)
Sets the vector using a string.
bool unitVector(const double tolerance=Kernel::Tolerance) const noexcept
double hklError() const
Calculates the error in hkl.
void readPrinted(std::istream &)
Read data from a stream in the format returned by printSelf ("[x,y,z]").
void write(std::ostream &) const
Write out the point values.
void getSpherical(double &R, double &theta, double &phi) const noexcept
Return the vector's position in spherical coordinates.
double toMillerIndexes(double eps=1.e-3)
transform vector into form, used to describe directions in crystallogaphical coodinate system
double cosAngle(const V3D &) const
cos(Angle) between this and another vector
void round() noexcept
Round each component to the nearest integer.
int masterDir(const double Tol=1e-3) const noexcept
Determine if there is a master direction.
void saveNexus(::NeXus::File *file, const std::string &name) const
Save the object to an open NeXus file.
bool nullVector(const double tolerance=1e-3) const noexcept
Determine if the point is null.
constexpr double deg2rad
Defines units/enum for Crystal work.
MANTID_KERNEL_DLL std::ostream & operator<<(std::ostream &, CPUTimer &)
Convenience function to provide for easier debug printing.
MANTID_KERNEL_DLL std::istream & operator>>(std::istream &, Interpolation &)
Reads in parameter value.
constexpr double Tolerance
Standard tolerance value.
MANTID_KERNEL_DLL V3D normalize(V3D v)
Normalizes a V3D.