From 2454e018856db2c6ae8cdbbe83de6a31f9584169 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Mon, 12 May 2025 10:58:53 +0200 Subject: [PATCH] proper class for GKT LF. --- src/external/LF_GL/CMakeLists.txt | 5 +- src/external/LF_GL/PGKT_LF.cpp | 265 ++++++++++++++++++++++++++++++ src/external/LF_GL/PGKT_LF.h | 64 ++++++++ src/external/LF_GL/main.cpp | 215 +----------------------- 4 files changed, 342 insertions(+), 207 deletions(-) create mode 100644 src/external/LF_GL/PGKT_LF.cpp create mode 100644 src/external/LF_GL/PGKT_LF.h diff --git a/src/external/LF_GL/CMakeLists.txt b/src/external/LF_GL/CMakeLists.txt index 22573406..8016f7dd 100644 --- a/src/external/LF_GL/CMakeLists.txt +++ b/src/external/LF_GL/CMakeLists.txt @@ -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 diff --git a/src/external/LF_GL/PGKT_LF.cpp b/src/external/LF_GL/PGKT_LF.cpp new file mode 100644 index 00000000..01c3580d --- /dev/null +++ b/src/external/LF_GL/PGKT_LF.cpp @@ -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 +#include +#include + +#include "PGKT_LF.h" + +//----------------------------------------------------------------------------- +/** + *

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 ¶m, std::vector &tt, std::vector &pol) : + fParam(param), fTime(tt), fPol(pol) +{ + fParam[0] *= fGammaMu; + + if (DynamicGaussKTLF() == 0) + fValid = true; +} + +//----------------------------------------------------------------------------- +/** + *

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; +} + +//----------------------------------------------------------------------------- +/** + *

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(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(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 p0exp(N); + double t = 0.0; + double dt = Tmax/N; + double nu = fParam[0] / (2.0*std::numbers::pi_v); + fDynLFdt = dt; // keep it since it is needed in GetDynKTLFValue() + for (unsigned int i=0; i 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; iCalculate static LF integral. + */ +void PGKT_LF::CalculateGaussLFIntegral() +{ + // fParam[0] = omega (field), fParam[1] = width + double nu = fParam[0] / (2.0*std::numbers::pi_v); + + 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)); +} + +//----------------------------------------------------------------------------- +/** + *

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(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(idx)); + + return fLFIntegral[idx]+df; +} + +//----------------------------------------------------------------------------- +/** + *

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(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(idx)); + + return fDynLFFuncValue[idx]+df; +} diff --git a/src/external/LF_GL/PGKT_LF.h b/src/external/LF_GL/PGKT_LF.h new file mode 100644 index 00000000..47f9c9ba --- /dev/null +++ b/src/external/LF_GL/PGKT_LF.h @@ -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 + +class PGKT_LF +{ + public: + PGKT_LF(std::vector ¶m, std::vector &tt, std::vector &pol); + + int DynamicGaussKTLF(); + bool IsValid() { return fValid; } + + std::vector GetTime() { return fTime; } + std::vector GetPol() { return fPol; } + + private: + const double fGammaMu=8.5161545577e-2; + std::vector fDynLFFuncValue; + std::vector fLFIntegral; + double fSamplingTime{0.0001}; + double fDynLFdt{0.0001}; + + bool fValid{false}; + std::vector fParam; + std::vector fTime; + std::vector fPol; + + void CalculateDynKTLF(); + void CalculateGaussLFIntegral(); + double GetLFIntegralValue(const double t) const; + double GetDynKTLFValue(const double t) const; +}; + +#endif // _PGKT_LF_ diff --git a/src/external/LF_GL/main.cpp b/src/external/LF_GL/main.cpp index 7d0d1d60..5ee0113c 100644 --- a/src/external/LF_GL/main.cpp +++ b/src/external/LF_GL/main.cpp @@ -6,12 +6,7 @@ #include #include -// global variable ------------------------------------------------------------ -const double gamma_mu=8.5161545577e-2; -std::vector gDynLFFuncValue; -std::vector 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); - - 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(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(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(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(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 p0exp(N); - double t = 0.0; - double dt = Tmax/N; - double nu = param[0] / (2.0*std::numbers::pi_v); - gDynLFdt = dt; // keep it since it is needed in lf_getDynKTLFValue() - for (unsigned int i=0; i 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(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(idx)); - - return gDynLFFuncValue[idx]+df; -} - -// ---------------------------------------------------------------------------- -int lf_gl_calc_gaussian(const double param[3], std::vector &tt, std::vector &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 ¶m, const std::vector &tt, const std::vector &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 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 tt; std::vector 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; }