LMU/musrfit
LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 3 Activity

proper class for GKT LF.

Browse Source
This commit is contained in:
suter_a 2025-05-12 10:58:53 +02:00
parent bcc1597e30
commit 819d209863
4 changed files with 342 additions and 207 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
src/external/LF_GL/CMakeLists.txt vendored
View File

@ -55,7 +55,10 @@ message("")
message("-------------------------------------------------------------------------")
message("")
add_executable(lf_gl main.cpp)
add_executable(lf_gl
main.cpp
PGKT_LF.cpp
)
set_property(TARGET lf_gl PROPERTY CXX_STANDARD 20)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
target_include_directories(lf_gl

265
src/external/LF_GL/PGKT_LF.cpp vendored Normal file
View File

@ -0,0 +1,265 @@
/***************************************************************************
PGKT_LF.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2025 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#include <iostream>
#include <cmath>
#include <numbers>
#include "PGKT_LF.h"
//-----------------------------------------------------------------------------
/**
* <p>CTOR
*
* @param param param[0]=field (G), param[1]=width (1/us), param[2]=hopp (1/us)
* @param tt time vector in (ns)
* @param pol polarization vector
*/
PGKT_LF::PGKT_LF(std::vector<double> &param, std::vector<double> &tt, std::vector<double> &pol) :
fParam(param), fTime(tt), fPol(pol)
{
fParam[0] *= fGammaMu;
if (DynamicGaussKTLF() == 0)
fValid = true;
}
//-----------------------------------------------------------------------------
/**
* <p>Calculate GKT LF polarization values.
*
* @return 0 on success, >0 otherwise
*/
int PGKT_LF::DynamicGaussKTLF()
{
bool useKeren{false};
// check if there is an empty time vector, if yes create one, otherwise use the given one
if (fTime.size() == 0) {
double t = 0.0;
unsigned int i=0;
do {
t = i++ * 1.0e-4; // 100ps steps
fTime.push_back(t);
} while (t < 10.0);
}
fPol.resize(fTime.size());
if (fParam[2]/fParam[1] > 5.0) // hopp/width = nu/Delta > 5.0
useKeren=true;
unsigned int i{0};
if (useKeren) {
double wL = fParam[0];
double wL2 = wL*wL;
double nu2 = fParam[2]*fParam[2];
double Gamma_t{0.0};
for (auto t : fTime) {
Gamma_t = 2.0*fParam[1]*fParam[1]/((wL2+nu2)*(wL2+nu2))* ((wL2+nu2)*fParam[2]*t
+ (wL2-nu2)*(1.0 - exp(-fParam[2]*t)*cos(wL*t))
- 2.0*fParam[2]*wL*exp(-fParam[2]*t)*sin(wL*t));
fPol[i++] = exp(-Gamma_t);
}
} else {
CalculateDynKTLF();
for (auto t : fTime) {
fPol[i++] = GetDynKTLFValue(t);
}
}
return 0;
}
//-----------------------------------------------------------------------------
/**
* <p>Calculate GKT LF polarization values, where the Keren approximation fails.
*/
void PGKT_LF::CalculateDynKTLF()
{
const double Tmax = 20.0; // 20 us
unsigned int N = static_cast<unsigned int>(16.0*Tmax*fParam[0]);
// check if width is very high
if (fParam[1] > 0.1) {
double tmin = 20.0;
tmin = fabs(sqrt(3.0)/fParam[1]);
unsigned int Nrate = static_cast<unsigned int>(25.0 * Tmax / tmin);
if (Nrate > N) {
N = Nrate;
}
}
if (N < 300) // if too few points, i.e. hopp very small, take 300 points
N = 300;
if (N>1e6) // make sure that N is not too large to prevent memory overflow
N = 1e6;
fDynLFFuncValue.resize(N);
CalculateGaussLFIntegral();
// calculate the P^(0)(t) exp(-nu t) vector
std::vector<double> p0exp(N);
double t = 0.0;
double dt = Tmax/N;
double nu = fParam[0] / (2.0*std::numbers::pi_v<double>);
fDynLFdt = dt; // keep it since it is needed in GetDynKTLFValue()
for (unsigned int i=0; i<N; i++) {
if (nu < 0.02) { // if smaller 20kHz ~ 0.27G use zero field formula
double sigma_t_2 = t*t*fParam[1]*fParam[1];
p0exp[i] = 0.333333333333333 * (1.0 + 2.0*(1.0 - sigma_t_2)*exp(-0.5*sigma_t_2));
} else if (fParam[1]/nu > 79.5775) { // check if Delta/w0 > 500.0, in which case the ZF formula is used
double sigma_t_2 = t*t*fParam[1]*fParam[1];
p0exp[i] = 0.333333333333333 * (1.0 + 2.0*(1.0 - sigma_t_2)*exp(-0.5*sigma_t_2));
} else {
double width = fParam[1];
double w0 = fParam[0];
p0exp[i] = 1.0 - 2.0*pow(width/w0,2.0)*(1.0 - exp(-0.5*pow(width*t, 2.0))*cos(w0*t)) + GetLFIntegralValue(t);
}
p0exp[i] *= exp(-fParam[2]*t);
t += dt;
}
// solve the volterra equation (trapezoid integration)
fDynLFFuncValue[0]=p0exp[0];
double sum;
double a;
double preFactor = dt*fParam[2];
for (unsigned int i=1; i<N; i++) {
sum = p0exp[i];
sum += 0.5*preFactor*p0exp[i]*fDynLFFuncValue[0];
for (unsigned int j=1; j<i; j++) {
sum += preFactor*p0exp[i-j]*fDynLFFuncValue[j];
}
a = 1.0-0.5*preFactor*p0exp[0];
fDynLFFuncValue[i]=sum/a;
}
}
//-----------------------------------------------------------------------------
/**
* <p>Calculate static LF integral.
*/
void PGKT_LF::CalculateGaussLFIntegral()
{
// fParam[0] = omega (field), fParam[1] = width
double nu = fParam[0] / (2.0*std::numbers::pi_v<double>);
if (fParam[0] == 0.0) // field == 0
return;
double dt=0.001; // all times in usec
double t, ft;
double w0 = fParam[0];
double width = fParam[1];
double preFactor = 2.0*pow(width, 4.0) / pow(w0, 3.0);
// check if the time resolution needs to be increased
const unsigned int samplingPerPeriod = 20;
const unsigned int samplingOnExp = 3000;
if ((width <= w0) && (1.0/nu < 20.0)) { // makes sure that the frequency sampling is fine enough
if (1.0/nu/samplingPerPeriod < 0.001) {
dt = 1.0/nu/samplingPerPeriod;
}
} else if ((width > w0) && (width <= 10.0)) {
if (width/w0 > 10.0) {
dt = 0.00005;
}
} else if ((width > w0) && (width > 10.0)) { // makes sure there is a fine enough sampling for large Delta's
if (1.0/width/samplingOnExp < 0.001) {
dt = 1.0/width/samplingOnExp;
}
}
fSamplingTime = dt;
fLFIntegral.clear();
// calculate integral
t = 0.0;
fLFIntegral.push_back(0.0); // start value of the integral
ft = 0.0;
double step = 0.0, lastft = 1.0, diff = 0.0;
do {
t += dt;
step = 0.5*dt*preFactor*(exp(-0.5*pow(width * (t-dt), 2.0))*sin(w0*(t-dt))+
exp(-0.5*pow(width * t, 2.0))*sin(w0*t));
ft += step;
diff = fabs(fabs(lastft)-fabs(ft));
lastft = ft;
fLFIntegral.push_back(ft);
} while ((t <= 20.0) && (diff > 1.0e-10));
}
//-----------------------------------------------------------------------------
/**
* <p>Gets value of the non-analytic static LF integral.
*
* @param t time in (usec)
*
* @return interpolated LF integral value
*/
double PGKT_LF::GetLFIntegralValue(const double t) const
{
unsigned int idx = static_cast<unsigned int>(t/fSamplingTime);
if (idx + 2 > fLFIntegral.size())
return fLFIntegral.back();
// linearly interpolate between the two relevant function bins
double df = (fLFIntegral[idx+1]-fLFIntegral[idx])*(t/fSamplingTime-static_cast<double>(idx));
return fLFIntegral[idx]+df;
}
//-----------------------------------------------------------------------------
/**
* <p><p>Gets value of the dynamic Kubo-Toyabe LF function.
*
* @param t time (usec)
*
* @return interpolated LF value.
*/
double PGKT_LF::GetDynKTLFValue(const double t) const
{
unsigned int idx = static_cast<unsigned int>(t/fDynLFdt);
if (idx + 2 > fDynLFFuncValue.size())
return fDynLFFuncValue.back();
// linearly interpolate between the two relvant function bins
double df = (fDynLFFuncValue[idx+1]-fDynLFFuncValue[idx])*(t/fDynLFdt-static_cast<double>(idx));
return fDynLFFuncValue[idx]+df;
}

64
src/external/LF_GL/PGKT_LF.h vendored Normal file
View File

@ -0,0 +1,64 @@
/***************************************************************************
PGKT_LF.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2025 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef _PGKT_LF_
#define _PGKT_LF_
#include <vector>
class PGKT_LF
{
public:
PGKT_LF(std::vector<double> &param, std::vector<double> &tt, std::vector<double> &pol);
int DynamicGaussKTLF();
bool IsValid() { return fValid; }
std::vector<double> GetTime() { return fTime; }
std::vector<double> GetPol() { return fPol; }
private:
const double fGammaMu=8.5161545577e-2;
std::vector<double> fDynLFFuncValue;
std::vector<double> fLFIntegral;
double fSamplingTime{0.0001};
double fDynLFdt{0.0001};
bool fValid{false};
std::vector<double> fParam;
std::vector<double> fTime;
std::vector<double> fPol;
void CalculateDynKTLF();
void CalculateGaussLFIntegral();
double GetLFIntegralValue(const double t) const;
double GetDynKTLFValue(const double t) const;
};
#endif // _PGKT_LF_

215
src/external/LF_GL/main.cpp vendored
View File

@ -6,12 +6,7 @@
#include <cmath>
#include <numbers>
// global variable ------------------------------------------------------------
const double gamma_mu=8.5161545577e-2;
std::vector<double> gDynLFFuncValue;
std::vector<double> gLFIntegral;
double gSamplingTime = 0.0001;
double gDynLFdt = 0.0001;
#include "PGKT_LF.h"
// ----------------------------------------------------------------------------
void lf_gl_syntax()
@ -28,203 +23,12 @@ void lf_gl_syntax()
}
// ----------------------------------------------------------------------------
void lf_calc_int(const double param[3])
{
// param[0] = omega (field), param[1] = width
double nu = param[0] / (2.0*std::numbers::pi_v<double>);
if (param[0] == 0.0) // field == 0
return;
double dt=0.001; // all times in usec
double t, ft;
double w0 = param[0];
double width = param[1];
double preFactor = 2.0*pow(width, 4.0) / pow(w0, 3.0);
// check if the time resolution needs to be increased
const unsigned int samplingPerPeriod = 20;
const unsigned int samplingOnExp = 3000;
if ((width <= w0) && (1.0/nu < 20.0)) { // makes sure that the frequency sampling is fine enough
if (1.0/nu/samplingPerPeriod < 0.001) {
dt = 1.0/nu/samplingPerPeriod;
}
} else if ((width > w0) && (width <= 10.0)) {
if (width/w0 > 10.0) {
dt = 0.00005;
}
} else if ((width > w0) && (width > 10.0)) { // makes sure there is a fine enough sampling for large Delta's
if (1.0/width/samplingOnExp < 0.001) {
dt = 1.0/width/samplingOnExp;
}
}
gSamplingTime = dt;
gLFIntegral.clear();
// calculate integral
t = 0.0;
gLFIntegral.push_back(0.0); // start value of the integral
ft = 0.0;
double step = 0.0, lastft = 1.0, diff = 0.0;
do {
t += dt;
step = 0.5*dt*preFactor*(exp(-0.5*pow(width * (t-dt), 2.0))*sin(w0*(t-dt))+
exp(-0.5*pow(width * t, 2.0))*sin(w0*t));
ft += step;
diff = fabs(fabs(lastft)-fabs(ft));
lastft = ft;
gLFIntegral.push_back(ft);
} while ((t <= 20.0) && (diff > 1.0e-10));
}
// ----------------------------------------------------------------------------
double lf_getLFIntegralValue(const double t)
{
unsigned int idx = static_cast<unsigned int>(t/gSamplingTime);
if (idx + 2 > gLFIntegral.size())
return gLFIntegral.back();
// linearly interpolate between the two relevant function bins
double df = (gLFIntegral[idx+1]-gLFIntegral[idx])*(t/gSamplingTime-static_cast<double>(idx));
return gLFIntegral[idx]+df;
}
// ----------------------------------------------------------------------------
void lf_dyn_kt(const double param[3])
{
const double Tmax = 20.0; // 20 us
unsigned int N = static_cast<unsigned int>(16.0*Tmax*param[0]);
// check if width is very high
if (param[1] > 0.1) {
double tmin = 20.0;
tmin = fabs(sqrt(3.0)/param[1]);
unsigned int Nrate = static_cast<unsigned int>(25.0 * Tmax / tmin);
if (Nrate > N) {
N = Nrate;
}
}
if (N < 300) // if too few points, i.e. hopp very small, take 300 points
N = 300;
if (N>1e6) // make sure that N is not too large to prevent memory overflow
N = 1e6;
gDynLFFuncValue.resize(N);
lf_calc_int(param);
// calculate the P^(0)(t) exp(-nu t) vector
std::vector<double> p0exp(N);
double t = 0.0;
double dt = Tmax/N;
double nu = param[0] / (2.0*std::numbers::pi_v<double>);
gDynLFdt = dt; // keep it since it is needed in lf_getDynKTLFValue()
for (unsigned int i=0; i<N; i++) {
if (nu < 0.02) { // if smaller 20kHz ~ 0.27G use zero field formula
double sigma_t_2 = t*t*param[1]*param[1];
p0exp[i] = 0.333333333333333 * (1.0 + 2.0*(1.0 - sigma_t_2)*exp(-0.5*sigma_t_2));
} else if (param[1]/nu > 79.5775) { // check if Delta/w0 > 500.0, in which case the ZF formula is used
double sigma_t_2 = t*t*param[1]*param[1];
p0exp[i] = 0.333333333333333 * (1.0 + 2.0*(1.0 - sigma_t_2)*exp(-0.5*sigma_t_2));
} else {
double width = param[1];
double w0 = param[0];
p0exp[i] = 1.0 - 2.0*pow(width/w0,2.0)*(1.0 - exp(-0.5*pow(width*t, 2.0))*cos(w0*t)) + lf_getLFIntegralValue(t);
}
p0exp[i] *= exp(-param[2]*t);
t += dt;
}
// solve the volterra equation (trapezoid integration)
gDynLFFuncValue[0]=p0exp[0];
double sum;
double a;
double preFactor = dt*param[2];
for (unsigned int i=1; i<N; i++) {
sum = p0exp[i];
sum += 0.5*preFactor*p0exp[i]*gDynLFFuncValue[0];
for (unsigned int j=1; j<i; j++) {
sum += preFactor*p0exp[i-j]*gDynLFFuncValue[j];
}
a = 1.0-0.5*preFactor*p0exp[0];
gDynLFFuncValue[i]=sum/a;
}
}
// ----------------------------------------------------------------------------
double lf_getDynKTLFValue(const double t)
{
unsigned int idx = static_cast<unsigned int>(t/gDynLFdt);
if (idx + 2 > gDynLFFuncValue.size())
return gDynLFFuncValue.back();
// linearly interpolate between the two relvant function bins
double df = (gDynLFFuncValue[idx+1]-gDynLFFuncValue[idx])*(t/gDynLFdt-static_cast<double>(idx));
return gDynLFFuncValue[idx]+df;
}
// ----------------------------------------------------------------------------
int lf_gl_calc_gaussian(const double param[3], std::vector<double> &tt, std::vector<double> &pol)
{
bool useKeren{false};
// check if there is an empty time vector, if yes create one, otherwise use the given one
if (tt.size() == 0) {
double t = 0.0;
unsigned int i=0;
do {
t = i++ * 1.0e-4; // 100ps steps
tt.push_back(t);
} while (t < 10.0);
}
pol.resize(tt.size());
std::cout << "debug> tt.size()=" << tt.size() << ", tt[0]=" << tt[0] << ", tt[end]=" << tt[tt.size()-1] << std::endl;
if (param[2]/param[1] > 5.0) // hopp/width = nu/Delta > 5.0
useKeren=true;
unsigned int i{0};
if (useKeren) {
double wL = param[0];
double wL2 = wL*wL;
double nu2 = param[2]*param[2];
double Gamma_t{0.0};
for (auto t : tt) {
Gamma_t = 2.0*param[1]*param[1]/((wL2+nu2)*(wL2+nu2))* ((wL2+nu2)*param[2]*t
+ (wL2-nu2)*(1.0 - exp(-param[2]*t)*cos(wL*t))
- 2.0*param[2]*wL*exp(-param[2]*t)*sin(wL*t));
pol[i++] = exp(-Gamma_t);
}
} else {
lf_dyn_kt(param);
for (auto t : tt) {
pol[i++] = lf_getDynKTLFValue(t);
}
}
return 0;
}
// ----------------------------------------------------------------------------
int lf_gl_write(const std::string fln, const double field, const double param[3],
int lf_gl_write(const std::string fln, const std::vector<double> &param,
const std::vector<double> &tt, const std::vector<double> &pol)
{
std::ofstream fout(fln, std::ofstream::out);
fout << "# field=" << field << " (G), width=" << param[1] << " (1/us), hopp=" << param[2] << " (1/us)" << std::endl;
fout << "# field=" << param[0] << " (G), width=" << param[1] << " (1/us), hopp=" << param[2] << " (1/us)" << std::endl;
fout << "# t (us), pol" << std::endl;
for (unsigned int i=0; i<tt.size(); i++) {
fout << tt[i] << ", " << pol[i] << std::endl;
@ -237,7 +41,7 @@ int lf_gl_write(const std::string fln, const double field, const double param[3]
// ----------------------------------------------------------------------------
int main(int argc, char *argv[])
{
double param[3] = {0.0, 0.0, 0.0};
std::vector<double> param{0.0, 0.0, 0.0};
double field{0.0};
bool gaussian_only = false;
std::string flnOut="";
@ -308,22 +112,21 @@ int main(int argc, char *argv[])
std::cout << "Gaussian averaged -> Lorentz LF" << std::endl;
std::cout << "flnOut: " << flnOut << std::endl;
// field -> omega
field = param[0];
param[0] *= gamma_mu;
std::vector<double> tt;
std::vector<double> pol;
if (gaussian_only) {
if (lf_gl_calc_gaussian(param, tt, pol) != 0) {
PGKT_LF gkt_lf(param, tt, pol);
if (!gkt_lf.IsValid()) {
std::cout << "**ERROR** in Gaussian LF calculation" << std::endl;
return 2;
}
tt = gkt_lf.GetTime();
pol = gkt_lf.GetPol();
} else {
std::cout << "not yet implemented." << std::endl;
}
lf_gl_write(flnOut, field, param, tt, pol);
lf_gl_write(flnOut, param, tt, pol);
return 0;
}