Verified Commit 5ae38229 authored by Tobias WEBER's avatar Tobias WEBER
Browse files

longitudinal fluctuations

parent f2254436
......@@ -89,7 +89,7 @@ prepare:
# -----------------------------------------------------------------------------
# Takin plugin modules
# -----------------------------------------------------------------------------
lib/skxmod.so: src/takin/takin.o src/core/skx.o src/core/fp.o src/core/heli.o src/core/magsys.o \
lib/skxmod.so: src/takin/takin.o src/core/skx.o src/core/fp.o src/core/heli.o src/core/longfluct.o src/core/magsys.o \
ext/takin/tools/monteconvo/sqwbase.o ext/tlibs2/libs/log.o
@echo "Linking Takin module $@..."
$(CXX) $(STD) $(OPT) $(DEFS) $(LIBDIRS) $(LIBDEFS) -shared -o $@ $+ -llapacke
......
/**
* Longitudinal fluctuations
* @author Tobias Weber <tweber@ill.fr>
* @date 13-jul-20
* @license GPLv2 (see 'LICENSE' file)
*/
#include "longfluct.h"
#include "constants.h"
Longfluct::Longfluct()
{
SetPinning(1,-1,0, 0,0,1);
}
/**
* generate skx satellites depending on the pinning and the up (i.e. B) vector
*/
void Longfluct::SetPinning(t_real Px, t_real Py, t_real Pz,
t_real upx, t_real upy, t_real upz)
{
m_up = tl2::make_vec<t_vec>({upx, upy, upz});
m_up /= tl2::veclen(m_up);
t_vec sat = tl2::make_vec<t_vec>({Px, Py, Pz});
sat = sat / tl2::veclen(sat) * g_kh_rlu_29K<t_real>;
t_mat rot = tl2::rotation_matrix(m_up, tl2::d2r(60.));
m_sats.clear();
m_sats.push_back(sat);
for(int i=1; i<6; ++i)
m_sats.emplace_back(tl2::prod_mv(rot, m_sats[i-1]));
}
/**
* calculates the dynamical structure factor for the longitudinal fluctuations
* using the formula from M. Garst, 10/jul/2020.
*/
Longfluct::t_real Longfluct::S_para(const Longfluct::t_vec& q, Longfluct::t_real E) const
{
if(tl2::float_equal<t_real>(m_A, 0))
return 0.;
constexpr t_real Ec2 = 0.04;
constexpr t_real kh = g_kh_rlu_29K<t_real>;
t_real S = 0.;
for(const t_vec& sat : m_sats)
{
S += g_kB<t_real> * m_T * m_A /
(
(m_inv_correl*m_inv_correl + tl2::inner(q-sat, q-sat)) / (kh*kh) +
std::pow(m_Gamma*E / Ec2, 2.)
);
}
return S;
}
/**
* Longitudinal fluctuations
* @author Tobias Weber <tweber@ill.fr>
* @date 13-jul-20
* @license GPLv2 (see 'LICENSE' file)
*/
#define USE_LAPACK
#include "tlibs2/libs/math17.h"
#include <vector>
class Longfluct
{
public:
using t_real = double;
using t_vec = ublas::vector<t_real>;
using t_mat = ublas::matrix<t_real>;
public:
Longfluct();
void SetPinning(t_real Px=1, t_real Py=-1, t_real Pz=0,
t_real upx=0, t_real upy=0, t_real upz=1);
t_real S_para(const t_vec& q, t_real E) const;
void SetT(t_real T) { m_T = T; }
void SetGamma(t_real G) { m_Gamma = G; }
void SetInvCorrel(t_real k) { m_inv_correl = k; }
void SetA(t_real A) { m_A = A; }
t_real GetT() const { return m_T; }
t_real GetGamma() const { return m_Gamma; }
t_real GetInvCorrel() const { return m_inv_correl; }
t_real GetA() const { return m_A; }
private:
// skx satellites
std::vector<t_vec> m_sats;
t_vec m_up;
t_real m_T = 29.;
t_real m_Gamma = 0.5;
t_real m_inv_correl = 0.1;
t_real m_A = 0.;
};
......@@ -18,6 +18,7 @@
#include "core/skx.h"
#include "core/fp.h"
#include "core/heli.h"
#include "core/longfluct.h"
#include "tlibs2/libs/phys.h"
#include "tlibs2/libs/str.h"
......@@ -48,6 +49,7 @@ protected:
Skx<t_real, t_cplx, SKX_ORDER> m_skx;
FP<t_real, t_cplx> m_fp;
Heli<t_real, t_cplx, HELI_ORDER> m_heli;
Longfluct m_lf;
t_real m_dErange = 2.; // E range around queried E to take into account
t_real m_dSigma = t_real(0.05);
......@@ -57,10 +59,12 @@ protected:
t_real m_dT = 29.;
t_real m_dB = 0.35;
t_real m_dcut = 0.02;
int m_iOnlyMode = -1;
int m_iPolChan = 0;
int m_iwhich_disp = 0; // 0: skx, 1: fp, 2: heli
int m_iProjNeutron = 1;
int m_ionlylf = 0;
t_vec m_vecG = tl2::make_vec<t_vec>({1,1,0});
t_vec m_vecB = tl2::make_vec<t_vec>({1,1,0});
......@@ -100,6 +104,7 @@ SqwMod::SqwMod()
m_fp.SetB(0.35);
m_heli.SetT(29.);
m_heli.SetB(0.35);
m_lf.SetT(29.);
std::vector<ublas::vector<t_cplx>> fourier_skx;
fourier_skx.reserve(_skxgs_allcomps.size()/3);
......@@ -120,6 +125,8 @@ SqwMod::SqwMod()
m_heli.SetCoords(m_vecB[0],m_vecB[1],m_vecB[2]);
m_heli.SetG(m_vecG[0], m_vecG[1], m_vecG[2]);
m_heli.SetFilterZeroWeight(true);
m_lf.SetPinning(m_vecPin[0],m_vecPin[1],m_vecPin[2], m_vecB[0],m_vecB[1],m_vecB[2]);
SqwBase::m_bOk = 1;
}
......@@ -150,6 +157,13 @@ std::tuple<std::vector<t_real>, std::vector<t_real>>
t_real SqwMod::operator()(t_real dh, t_real dk, t_real dl, t_real dE) const
{
t_vec vecq = tl2::make_vec<t_vec>({dh, dk, dl}) - m_vecG;
// longitudinal fluctuations
t_real dS_lf = m_lf.S_para(vecq, dE);
if(m_ionlylf)
return dS_lf;
std::vector<t_real> vecE, vecW[4];
// only calculate dispersion if global weight factor is not 0
if(!tl2::float_equal(m_dS0, t_real(0)))
......@@ -203,6 +217,11 @@ std::vector<SqwMod::t_var> SqwMod::GetVars() const
vecVars.push_back(SqwBase::t_var{"B_dir", "vector", vec_to_str(m_vecB)});
vecVars.push_back(SqwBase::t_var{"Pin_dir", "vector", vec_to_str(m_vecPin)});
vecVars.push_back(SqwBase::t_var{"lf_only", "int", tl2::var_to_str(m_ionlylf)});
vecVars.push_back(SqwBase::t_var{"lf_invcorrel", "real", tl2::var_to_str(m_lf.GetInvCorrel())});
vecVars.push_back(SqwBase::t_var{"lf_A", "real", tl2::var_to_str(m_lf.GetA())});
vecVars.push_back(SqwBase::t_var{"lf_gamma", "real", tl2::var_to_str(m_lf.GetGamma())});
return vecVars;
}
......@@ -236,6 +255,8 @@ void SqwMod::SetVars(const std::vector<SqwMod::t_var>& vecVars)
m_fp.SetT(m_dT);
m_heli.SetT(m_dT);
// fixed (and theo units) for skx!
m_lf.SetT(m_dT);
}
else if(strVar == "B")
{
......@@ -271,6 +292,9 @@ void SqwMod::SetVars(const std::vector<SqwMod::t_var>& vecVars)
m_heli.SetCoords(m_vecB[0],m_vecB[1],m_vecB[2]);
m_heli.SetG(m_vecG[0], m_vecG[1], m_vecG[2]);
m_lf.SetPinning(m_vecPin[0],m_vecPin[1],m_vecPin[2],
m_vecB[0],m_vecB[1],m_vecB[2]);
}
else if(strVar == "Pin_dir")
{
......@@ -278,6 +302,29 @@ void SqwMod::SetVars(const std::vector<SqwMod::t_var>& vecVars)
m_skx.SetCoords(m_vecB[0],m_vecB[1],m_vecB[2], m_vecPin[0],m_vecPin[1],m_vecPin[2]);
m_skx.SetG(m_vecG[0], m_vecG[1], m_vecG[2]);
m_lf.SetPinning(m_vecPin[0],m_vecPin[1],m_vecPin[2],
m_vecB[0],m_vecB[1],m_vecB[2]);
}
else if(strVar == "lf_only")
{
m_ionlylf = tl2::str_to_var<decltype(m_ionlylf)>(strVal);
}
else if(strVar == "lf_invcorrel")
{
t_real dInvCorrel = tl2::str_to_var<t_real>(strVal);
m_lf.SetInvCorrel(dInvCorrel);
}
else if(strVar == "lf_A")
{
t_real dA = tl2::str_to_var<t_real>(strVal);
m_lf.SetA(dA);
}
else if(strVar == "lf_gamma")
{
t_real dG = tl2::str_to_var<t_real>(strVal);
m_lf.SetGamma(dG);
}
}
}
......@@ -308,6 +355,8 @@ SqwBase* SqwMod::shallow_copy() const
pMod->m_skx = this->m_skx;
pMod->m_fp = this->m_fp;
pMod->m_heli = this->m_heli;
pMod->m_lf = this->m_lf;
pMod->m_ionlylf = this->m_ionlylf;
pMod->m_iPolChan = this->m_iPolChan;
pMod->m_vecG = this->m_vecG;
pMod->m_vecB = this->m_vecB;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment