Skip to content
GitLab
Commit ce57b567 authored by Tobias WEBER's avatar Tobias WEBER
Browse files

added fixes from https://github.com/t-weber/tlibs

parent 5425851e
Hide whitespace changes
Inline Side-by-side
LITERATURE
View file @ ce57b567
......@@ -8,6 +8,7 @@ Neutrons
* (Shirane 2002), G. Shirane et al., "Neutron Scattering with a Triple-Axis Spectrometer", 2002, ISBN: 978-0521411264.
* (Squires 2012), G. L. Squires, "Thermal Neutron Scattering", 2012, ISBN: 9781139107808.
* (ILL Neutron Data Booklet), ISBN: 0-9704143-7-4, 2003, https://www.ill.eu/about-the-ill/documentation/
-------------------------------------------------------------------------------
......
math/geo.h
View file @ ce57b567
......@@ -10,7 +10,6 @@
#include "../helper/flags.h"
#include "../helper/exception.h"
//#include "quat.h"
#include "linalg.h"
#include "linalg2.h"
#include "stat.h"
......@@ -37,6 +36,7 @@ public:
using t_vec = ublas::vector<T>;
using t_mat = ublas::matrix<T>;
protected:
bool m_bValid = 0;
t_vec m_vecX0;
......@@ -44,6 +44,7 @@ protected:
t_vec m_vecNorm;
T m_d;
public:
/**
* plane from a point and a normal
......@@ -258,21 +259,27 @@ public:
};
//------------------------------------------------------------------------------
template<typename T> class Line
{
public:
using t_vec = ublas::vector<T>;
using t_mat = ublas::matrix<T>;
protected:
t_vec m_vecX0;
t_vec m_vecDir;
public:
Line() {}
Line(const t_vec& vec0, const t_vec& dir)
: m_vecX0(vec0), m_vecDir(dir)
{}
......@@ -291,6 +298,25 @@ public:
const t_vec& GetDir() const { return m_vecDir; }
/**
* distance to a point
*/
T GetDist(const t_vec& vecPt) const
{
const t_vec& vecX0 = GetX0();
t_vec vecDir = GetDir() / veclen(GetDir());
// shift everything so that line goes through the origin
t_vec vecPtShift = vecPt - vecX0;
// project point on direction vector
t_vec vecClosestPt = inner(vecPtShift, vecDir) * vecDir;
// distance between point and projected point
return veclen(vecClosestPt - vecPtShift);
}
/**
* distance to line l1
*/
......@@ -299,14 +325,14 @@ public:
const Line<T>& l0 = *this;
// vector normal to both directions defining the distance line
t_vec vecNorm = cross_3(l0.GetDir(), l1.GetDir());
t_vec vecNorm = cross_3<t_vec>(l0.GetDir(), l1.GetDir());
T tlenNorm = veclen(vecNorm);
t_vec vec01 = l1.GetX0() - l0.GetX0();
// if the lines are parallel, any point (e.g. the X0s) can be used
// if the lines are parallel, any point (e.g. the x0s) can be used
if(float_equal(tlenNorm, T(0)))
return veclen(vec01);
return GetDist(l1.GetX0());
// project x0_1 - x0_0 onto vecNorm
T tdot = std::abs(inner(vec01, vecNorm));
......@@ -314,26 +340,6 @@ public:
}
/**
* distance to a point
*/
T GetDist(const t_vec& vecPt) const
{
const t_vec& vecX0 = GetX0();
const t_vec& vecDir = GetDir();
T tlenDir = veclen(vecDir);
// area of parallelogram spanned by vecDir and vecPt-vecX0
// == ||vecDir||*||vecPt-vecX0||*sin(th)
T tArea = veclen(cross_3(vecDir, vecPt-vecX0));
// length of perpendicular line from vecPt, also by sine theorem
// == ||vecPt-vecX0||*sin(th)
return tArea / tlenDir;
}
bool IsParallel(const Line<T>& line, T eps = tl::get_epsilon<T>()) const
{
return vec_is_collinear<t_vec>(GetDir(), line.GetDir(), eps);
......@@ -420,7 +426,7 @@ public:
const t_vec xp2 = plane.GetX0() + plane.GetDir1();
t_mat matDenom(N+1,N+1);
matDenom(0,0) = 1; matDenom(0,1) = 1; matDenom(0,2) = 1; matDenom(0,3) = 0;
matDenom(0,0) = 1; matDenom(0,1) = 1; matDenom(0,2) = 1; matDenom(0,3) = 0;
matDenom(1,0) = xp0[0]; matDenom(1,1) = xp1[0]; matDenom(1,2) = xp2[0]; matDenom(1,3) = dirl[0];
matDenom(2,0) = xp0[1]; matDenom(2,1) = xp1[1]; matDenom(2,2) = xp2[1]; matDenom(2,3) = dirl[1];
matDenom(3,0) = xp0[2]; matDenom(3,1) = xp1[2]; matDenom(3,2) = xp2[2]; matDenom(3,3) = dirl[2];
......@@ -430,7 +436,7 @@ public:
return false;
t_mat matNum(N+1,N+1);
matNum(0,0) = 1; matNum(0,1) = 1; matNum(0,2) = 1; matNum(0,3) = 1;
matNum(0,0) = 1; matNum(0,1) = 1; matNum(0,2) = 1; matNum(0,3) = 1;
matNum(1,0) = xp0[0]; matNum(1,1) = xp1[0]; matNum(1,2) = xp2[0]; matNum(1,3) = posl[0];
matNum(2,0) = xp0[1]; matNum(2,1) = xp1[1]; matNum(2,2) = xp2[1]; matNum(2,3) = posl[1];
matNum(3,0) = xp0[2]; matNum(3,1) = xp1[2]; matNum(3,2) = xp2[2]; matNum(3,3) = posl[2];
......@@ -448,19 +454,27 @@ public:
*
* pos0 + t0*dir0 = pos1 + t1*dir1
* pos0 - pos1 = t1*dir1 - t0*dir0
* exact: b = Mx -> M^(-1)*b = x
* approx: M^t b = M^t M x -> (M^t M)^(-1) * M^t b = x
* pos0 - pos1 = (-dir0 dir1) * (t0 t1)^t
* exact solution for the t params: (-dir0 dir1)^(-1) * (pos0 - pos1) = (t0 t1)^t
*
* generally:
* exact: b = M t -> M^(-1)*b = t
* approx.: M^t b = M^t M t -> (M^t M)^(-1) * M^t b = t
*/
bool intersect(const Line<T>& line, T& t, T eps = tl::get_epsilon<T>()) const
bool intersect(const Line<T>& line1, T& t0, T eps = tl::get_epsilon<T>(), T *pt1=nullptr) const
{
if(IsParallel(line, eps))
if(IsParallel(line1, eps))
{
t0 = 0;
if(pt1) *pt1 = 0;
return false;
}
const t_vec& pos0 = this->GetX0();
const t_vec& pos1 = line.GetX0();
const t_vec& pos1 = line1.GetX0();
const t_vec& dir0 = this->GetDir();
const t_vec& dir1 = line.GetDir();
const t_vec& dir1 = line1.GetDir();
const t_vec pos = pos0-pos1;
t_mat M = column_matrix({-dir0, dir1});
......@@ -473,9 +487,17 @@ public:
t_vec Mtb = prod_mv(Mt, pos);
t_vec params = prod_mv(MtMinv, Mtb);
t = params[0];
return true;
// get parameters
t0 = params[0];
T t1 = params[1];
if(pt1) *pt1 = t1;
// get closest points between the two lines
t_vec vecInters0 = (*this)(t0);
t_vec vecInters1 = line1(t1);
return veclen(vecInters1-vecInters0) <= eps;
}
......@@ -501,6 +523,7 @@ public:
return true;
}
/**
* middle perpendicular plane (in 3d)
*/
......@@ -586,9 +609,6 @@ bool intersect_line_poly(const Line<T>& line,
// is intersection point within polygon?
const t_vec vecFaceCentre = mean_value(vertPoly);
//std::ostringstream ostr;
//ostr << "intersect: " << vecIntersect << ", face centre: " << vecFaceCentre << "\n";
for(std::size_t iVert = 0; iVert < vertPoly.size(); ++iVert)
{
std::size_t iNextVert = iVert < (vertPoly.size()-1) ? iVert+1 : 0;
......@@ -604,9 +624,6 @@ bool intersect_line_poly(const Line<T>& line,
if(planeEdge.GetDist(vecIntersect) < -eps)
return false;
//ostr << "Face " << iVert << ": " << vecEdgeCentre << ", " << vecNorm << ", "
// << planeEdge.GetDist(vecIntersect) << "\n";
}
return true;
......@@ -740,7 +757,6 @@ t_vec get_face_normal(const t_cont<t_vec>& vecVerts, t_vec vecCentre,
if(inner(vecNorm, vecCentreToFace) < T(0))
vecNorm = -vecNorm;
//tl::log_debug("vec = ", vecCentreToFace, ",\tnorm = ", vecNorm);
return vecNorm;
}
......@@ -799,7 +815,7 @@ t_cont<t_cont<t_vec>> verts_to_polyhedron(
t_cont<t_vec>& vecPoly = vecPolys[iPlane];
// plane already in set?
if(plane.IsOnPlane(vert0, T(10)*eps) &&
if(plane.IsOnPlane(vert0, T(10)*eps) &&
plane.IsOnPlane(vert1, T(10)*eps) &&
plane.IsOnPlane(vert2, T(10)*eps))
{
......@@ -864,6 +880,7 @@ public:
using t_vec = ublas::vector<T>;
using t_mat = ublas::matrix<T>;
protected:
// x^T Q x + r x + s = 0
t_mat m_Q = zero_m<t_mat>(3,3);
......@@ -873,26 +890,36 @@ protected:
t_vec m_vecOffs = zero_v<t_vec>(3);
bool m_bQSymm = 1;
protected:
void CheckSymm()
{
m_bQSymm = is_symmetric(m_Q, std::cbrt(get_epsilon<T>()));
}
public:
Quadric() {}
Quadric(std::size_t iDim)
: m_Q(zero_m<t_mat>(iDim,iDim)), m_r(zero_v<t_vec>(iDim))
{ CheckSymm(); }
Quadric(const t_mat& Q) : m_Q(Q)
{ CheckSymm(); }
Quadric(const t_mat& Q, const t_vec& r, T s)
: m_Q(Q), m_r(r), m_s(s)
{ CheckSymm(); }
~Quadric() {}
void SetDim(std::size_t iDim) { m_Q.resize(iDim, iDim, 1); }
const Quadric<T>& operator=(const Quadric<T>& quad)
{
this->m_Q = quad.m_Q;
......@@ -904,6 +931,7 @@ public:
return *this;
}
Quadric<T>& operator=(Quadric<T>&& quad)
{
this->m_Q = std::move(quad.m_Q);
......@@ -915,20 +943,25 @@ public:
return *this;
}
Quadric(const Quadric<T>& quad) { *this = quad; }
Quadric(Quadric<T>&& quad) { *this = quad; }
void SetOffset(const t_vec& vec) { m_vecOffs = vec; }
const t_vec& GetOffset() const { return m_vecOffs; }
const t_mat& GetQ() const { return m_Q; }
const t_vec& GetR() const { return m_r; }
T GetS() const { return m_s; }
void SetQ(const t_mat& Q) { m_Q = Q; CheckSymm(); }
void SetR(const t_vec& r) { m_r = r; }
void SetS(T s) { m_s = s; }
T operator()(const t_vec& _x) const
{
t_vec x = _x-m_vecOffs;
......@@ -940,6 +973,7 @@ public:
return dQ + dR + m_s;
}
/**
* remove column and row iIdx
*/
......@@ -950,6 +984,7 @@ public:
m_vecOffs = remove_elem(m_vecOffs, iIdx);
}
void transform(const t_mat& S)
{
m_Q = tl::transform<t_mat>(m_Q, S, 1);
......
phys/mag.h
View file @ ce57b567
......@@ -43,7 +43,6 @@ T ferromag(const t_cont<t_vec>& vecNeighbours, const t_cont<std::complex<T>>& ve
return T(2)*tS*(J0 - J).real();
}
template<class t_vec = ublas::vector<double>,
typename T = typename t_vec::value_type,
typename t_cont = std::initializer_list<std::tuple<t_vec, std::complex<T>>>>
......@@ -56,10 +55,9 @@ T ferromag(const t_cont& lstNeighbours, const ublas::vector<T>& vecq, T tS)
// ----------------------------------------------------------------------------
/**
* Magnetic form factors
* @desc see: ILL Neutron Data Booklet sec. 2.5-1 (p. 60)
* @desc see: (ILL Neutron Data Booklet), sec. 2.5-1 (p. 60)
* @desc also see: https://www.ill.eu/sites/ccsl/ffacts/
*/
template<class T=double, template<class...> class t_vec=std::initializer_list>
......@@ -75,7 +73,6 @@ T j0_avg(T Q, const t_vec<T>& A, const t_vec<T>& a)
return tl::formfact<T, std::vector>(Q, vecA, vecB, c);
}
template<class T=double, template<class...> class t_vec=std::initializer_list>
T j2_avg(T Q, const t_vec<T>& A, const t_vec<T>& a)
{
......@@ -83,7 +80,6 @@ T j2_avg(T Q, const t_vec<T>& A, const t_vec<T>& a)
return j0_avg<T, t_vec>(Q, A, a) * s * s;
}
template<class T=double, template<class...> class t_vec=std::initializer_list>
T mag_formfact(T Q, T L, T S,
const t_vec<T>& A0, const t_vec<T>& a0,
......@@ -92,7 +88,6 @@ T mag_formfact(T Q, T L, T S,
return (L+T(2)*S) * j0_avg<T, t_vec>(Q, A0, a0) * L * j2_avg<T, t_vec>(Q, A2, a2);
}
/**
* form factor for transition metals (d orbitals, weak LS, spin-only)
* @desc see: (Squires 2012), p. 138
......@@ -109,7 +104,6 @@ std::tuple<T,T,T> mag_formfact_d(T Q, T g,
return std::tuple<T,T,T>(j0, j2, j0 + (T(1)-T(2)/g)*j2);
}
/**
* form factor for rare earths (f orbitals, strong LS, jj)
* @desc see: (Squires 2012), p. 139
......
phys/neutrons.h
View file @ ce57b567
......@@ -45,8 +45,8 @@ template<typename T=double> T get_E2KSQ()
template<class T=double> T t_E2KSQ = T(1)/t_KSQ2E<T>;
#endif
// --------------------------------------------------------------------------------
// --------------------------------------------------------------------------------
// --------------------------------------------------------------------------------
......@@ -59,7 +59,6 @@ t_momentum<Sys,Y> lam2p(const t_length<Sys,Y>& lam)
return get_h<Y>() / lam;
}
template<class Sys, class Y>
t_length<Sys,Y> p2lam(const t_momentum<Sys,Y>& p)
{
......@@ -74,35 +73,30 @@ t_length<Sys,Y> k2lam(const t_wavenumber<Sys,Y>& k)
return Y(2.)*get_pi<Y>() / k;
}
template<class Sys, class Y>
t_wavenumber<Sys,Y> lam2k(const t_length<Sys,Y>& lam)
{
return Y(2.)*get_pi<Y>() / lam;
}
template<class Sys, class Y>
t_momentum<Sys,Y> k2p(const t_wavenumber<Sys,Y>& k)
{
return get_hbar<Y>()*k;
}
template<class Sys, class Y>
t_wavenumber<Sys,Y> p2k(const t_momentum<Sys,Y>& p)
{
return p/get_hbar<Y>();
}
template<class Sys, class Y>
t_velocity<Sys,Y> k2v(const t_wavenumber<Sys,Y>& k)
{
return k2p(k) / get_m_n<Y>();
}
template<class Sys, class Y>
t_wavenumber<Sys,Y> v2k(const t_velocity<Sys,Y>& v)
{
......@@ -122,14 +116,12 @@ t_energy<Sys,Y> omega2E(const t_freq<Sys,Y>& omega)
return get_hbar<Y>() * omega;
}
template<class Sys, class Y>
t_freq<Sys,Y> E2omega(const t_energy<Sys,Y>& en)
{
return en / get_hbar<Y>();
}
template<class Sys, class Y>
t_energy<Sys,Y> k2E_direct(const t_wavenumber<Sys,Y>& k)
{
......@@ -138,7 +130,6 @@ t_energy<Sys,Y> k2E_direct(const t_wavenumber<Sys,Y>& k)
return E;
}
template<class Sys, class Y>
t_wavenumber<Sys,Y> E2k_direct(const t_energy<Sys,Y>& _E, bool &bImag)
{
......@@ -153,7 +144,6 @@ t_wavenumber<Sys,Y> E2k_direct(const t_energy<Sys,Y>& _E, bool &bImag)
// --------------------------------------------------------------------------------
// --------------------------------------------------------------------------------
// indirect calculations using conversion factors for numerical stability
......@@ -192,7 +182,6 @@ t_length<Sys,Y> bragg_real_lam(const t_length<Sys,Y>& d,
return Y(2.)*d/n * units::sin(twotheta/Y(2.));
}
template<class Sys, class Y>
t_length<Sys,Y> bragg_real_d(const t_length<Sys,Y>& lam,
const t_angle<Sys,Y>& twotheta, Y n = Y(1))
......@@ -200,7 +189,6 @@ t_length<Sys,Y> bragg_real_d(const t_length<Sys,Y>& lam,
return n * lam / (Y(2.)*units::sin(twotheta/Y(2.)));
}
template<class Sys, class Y>
t_angle<Sys,Y> bragg_real_twotheta(const t_length<Sys,Y>& d,
const t_length<Sys,Y>& lam, Y n = Y(1))
......@@ -225,7 +213,6 @@ t_angle<Sys,Y> bragg_recip_twotheta(const t_wavenumber<Sys,Y>& G,
return units::asin(dS) * Y(2);
}
template<class Sys, class Y>
t_wavenumber<Sys,Y> bragg_recip_G(const t_length<Sys,Y>& lam,
const t_angle<Sys,Y>& twotheta, Y n = Y(1))
......@@ -233,13 +220,11 @@ t_wavenumber<Sys,Y> bragg_recip_G(const t_length<Sys,Y>& lam,
return Y(4)*get_pi<Y>() / (n*lam) * units::sin(twotheta/Y(2));
}
template<class Sys, class Y>
t_wavenumber<Sys,Y> bragg_recip_Q(const t_length<Sys,Y>& lam,
const t_angle<Sys,Y>& twotheta, Y n = Y(1))
{ return bragg_recip_G<Sys,Y>(lam,twotheta,n); }
template<class Sys, class Y>
t_length<Sys,Y> bragg_recip_lam(const t_wavenumber<Sys,Y>& G,
const t_angle<Sys,Y>& twotheta, Y n = Y(1))
......@@ -258,7 +243,6 @@ t_wavenumber<Sys,Y> bragg_recip_G(const t_wavenumber<Sys,Y>& k,
return Y(2)*k / n * units::sin(twotheta/Y(2));
}
template<class Sys, class Y>
t_wavenumber<Sys,Y> bragg_recip_k(const t_wavenumber<Sys,Y>& G,
const t_angle<Sys,Y>& twotheta, Y n = Y(1))
......@@ -266,7 +250,6 @@ t_wavenumber<Sys,Y> bragg_recip_k(const t_wavenumber<Sys,Y>& G,
return n*G / (Y(2) * units::sin(twotheta/Y(2)));
}
template<class Sys, class Y>
t_angle<Sys,Y> bragg_recip_twotheta(const t_wavenumber<Sys,Y>& G,
const t_wavenumber<Sys,Y>& k, Y n = Y(1))
......@@ -286,7 +269,6 @@ t_length<Sys,Y> G2d(const t_wavenumber<Sys,Y>& G)
return Y(2.)*get_pi<Y>() / G;
}
template<class Sys, class Y>
t_wavenumber<Sys,Y> d2G(const t_length<Sys,Y>& d)
{
......@@ -296,11 +278,10 @@ t_wavenumber<Sys,Y> d2G(const t_length<Sys,Y>& d)
// --------------------------------------------------------------------------------
// --------------------------------------------------------------------------------
/**
* differentiated Bragg equation:
* n lam = 2d sin(th) | diff
* n lam = 2d sin(th) | diff
* n dlam = 2dd sin(th) + 2d cos(th) dth | / Bragg equ
* dlam/lam = dd/d + cos(th)/sin(th) dth
*
......@@ -321,64 +302,88 @@ Y bragg_diff(Y dDoverD, const t_angle<Sys,Y>& theta, Y dTheta)
// --------------------------------------------------------------------------------
// see e.g. ILL blue book sec. 2.6-2
/**
* kinematic plane
* see e.g. (ILL Neutron Data Booklet), sec. 2.6-2
*
* Q_vec = ki_vec - kf_vec
* Q^2 = ki^2 + kf^2 - 2ki kf cos 2th | * hbar^2 / (2 mn)
*
* using:
* Ei = hbar^2 ki^2 / (2 mn)
*
* Q^2 * hbar^2 / (2 mn) = Ei + Ef - 2 ki kf cos(2th) * hbar^2 / (2 mn)
*
* using:
* ki^2 = 2 mn Ei / hbar^2
*
* Q^2 = [Ei + Ef - 2 sqrt(Ei) sqrt(Ef) cos 2th] * 2 mn / hbar^2
*
* using:
* dE = Ei - Ef
* Ef = Ei - dE
*
* Q^2 = [2 Ei - dE - 2 sqrt(Ei (Ei - dE)) cos 2th] * 2 mn / hbar^2
*/
template<class Sys, class Y>
t_wavenumber<Sys,Y> kinematic_plane(bool bFixedKi,
const t_energy<Sys,Y>& EiEf, const t_energy<Sys,Y>& DeltaE,
const t_angle<Sys,Y>& twotheta)
{
const t_energy<Sys,Y> dE = DeltaE;
t_energy<Sys,Y> dE = DeltaE;
if(bFixedKi)
dE = -dE;
auto c = Y(2.)*get_m_n<Y>() / (get_hbar<Y>()*get_hbar<Y>());
t_wavenumber<Sys,Y> Q =
units::sqrt(Y(2.)*get_m_n<Y>() / get_hbar<Y>()) *
units::sqrt(c *
(Y(2.)*EiEf + dE - Y(2.)*units::cos(twotheta) *
units::sqrt(EiEf*(EiEf + dE)));
units::sqrt(EiEf*(EiEf + dE))));
return Q;
}
/**
* kinematic plane
* see e.g. (ILL Neutron Data Booklet), sec. 2.6-2
*
* solving the above equation for dE using sage:
* Q, Ei, dE, ctt, c = var("Q, Ei, dE, ctt, c")
* equ = (Q^2 -2*Ei*c + dE*c)^2 == 4*Ei*(Ei-dE)*c^2*ctt^2
* equ.solve(dE)
*/
template<class Sys, class Y>
t_energy<Sys,Y> kinematic_plane(bool bFixedKi, bool bBranch,
const t_energy<Sys,Y>& EiEf, const t_wavenumber<Sys,Y>& Q,
const t_angle<Sys,Y>& twotheta)
{
auto c = Y(2.)*get_m_n<Y>() / (get_hbar<Y>()*get_hbar<Y>());
Y ctt = units::cos(twotheta);
Y c2tt = units::cos(Y(2.)*twotheta);
auto c2 = c*c;
Y dSign = Y(-1.);
if(bBranch)
dSign = Y(1.);
auto EiEf2 = EiEf*EiEf;
Y dSignFixedKf = Y(1.);
if(bFixedKi)
dSignFixedKf = Y(-1.);
Y ctt = units::cos(twotheta);
Y ctt2 = ctt*ctt;
using t_sqrt_rt = decltype(c*c*EiEf*ctt);
using t_rt = decltype(t_sqrt_rt()*t_sqrt_rt());
t_rt rt = c*c*c*c * (-EiEf*EiEf)*ctt*ctt
+ c*c*c*c*EiEf*EiEf*ctt*ctt*c2tt
+ Y(2.)*c*c*c*EiEf*Q*Q*ctt*ctt;
Y dSign = bBranch ? Y(1.) : Y(-1.);
Y dSignFixedKf = bFixedKi ? Y(-1.) : Y(1.);
t_energy<Sys,Y> E =
Y(1.)/(c*c)*(dSignFixedKf*Y(2.)*c*c*EiEf*ctt*ctt
- dSignFixedKf*Y(2.)*c*c*EiEf
+ dSign*std::sqrt(Y(2.)) * my_units_sqrt<t_sqrt_rt>(rt)