From bcc1597e30dd20b7713e2abe9eb658cd92b085b0 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Fri, 9 May 2025 16:02:04 +0200 Subject: [PATCH] add LF Gaussian/Lorentzian testing code. --- src/external/LF_GL/CMakeLists.txt | 65 ++++++ src/external/LF_GL/main.cpp | 329 ++++++++++++++++++++++++++++++ 2 files changed, 394 insertions(+) create mode 100644 src/external/LF_GL/CMakeLists.txt create mode 100644 src/external/LF_GL/main.cpp diff --git a/src/external/LF_GL/CMakeLists.txt b/src/external/LF_GL/CMakeLists.txt new file mode 100644 index 00000000..22573406 --- /dev/null +++ b/src/external/LF_GL/CMakeLists.txt @@ -0,0 +1,65 @@ +# - LF GL --------------------------------------------------------------------- + +cmake_minimum_required(VERSION 3.17) + +project(lf_gl VERSION 0.9 LANGUAGES C CXX) + +#--- set a default build type if none was specified --------------------------- +set(default_build_type "Release") + +if (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) + message(STATUS "Setting build type to '${default_build_type}' as none was specified.") + set(CMAKE_BUILD_TYPE "${default_build_type}" CACHE + STRING "Choose the type of build." FORCE) + # Set the possible values of build type for cmake-gui + set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS + "Debug" "Release" "MinSizeRel" "RelWithDebInfo") +endif () + +#--- check for boost ---------------------------------------------------------- +find_package(Boost REQUIRED + COMPONENTS + system + filesystem +) +message(STATUS "Boost libs: ${Boost_LIBRARIES}") + +#--- check for gsl ------------------------------------------------------------ +find_package(GSL REQUIRED) + +#--- write summary of the installation +cmake_host_system_information(RESULT PROCESSOR QUERY PROCESSOR_DESCRIPTION) + +message("") +message("|-----------------------------------------------------------------------|") +message("| |") +message("| Summary |") +message("| |") +message("|-----------------------------------------------------------------------|") +message("") +message(" System: ${CMAKE_HOST_SYSTEM_NAME} ${CMAKE_SYSTEM_PROCESSOR} - ${CMAKE_HOST_SYSTEM_VERSION}") +message(" Processor: ${PROCESSOR} (${CMAKE_SYSTEM_PROCESSOR})") +message(" ----------") +message("") +message(" lf_gl Version: ${lf_gl_VERSION}") +message(" --------------") +message("") +message(" Build Type: ${CMAKE_BUILD_TYPE}") +message(" -----------") +message("") +message("-------------------------------------------------------------------------") +message("") +message(" GSL found in ${GSL_INCLUDE_DIRS}, Version: ${GSL_VERSION}") +message(" BOOST found in ${Boost_INCLUDE_DIRS}, Version: ${Boost_VERSION}") +message("") +message("-------------------------------------------------------------------------") +message("") + +add_executable(lf_gl main.cpp) +set_property(TARGET lf_gl PROPERTY CXX_STANDARD 20) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +target_include_directories(lf_gl + BEFORE PRIVATE + $ +) +target_link_libraries(lf_gl -lm GSL::gsl) diff --git a/src/external/LF_GL/main.cpp b/src/external/LF_GL/main.cpp new file mode 100644 index 00000000..7d0d1d60 --- /dev/null +++ b/src/external/LF_GL/main.cpp @@ -0,0 +1,329 @@ +#include +#include +#include +#include +#include +#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; + +// ---------------------------------------------------------------------------- +void lf_gl_syntax() +{ + std::cout << std::endl; + std::cout << "usage lg_gl [-p field width hopp -g [1,0] -o flnOut] | [-h]" << std::endl; + std::cout << " -p: field in (G)" << std::endl; + std::cout << " width in (1/us)" << std::endl; + std::cout << " hopp in (1/us)" << std::endl; + std::cout << " -g: 1=Gaussian LF data, 0=Gaussian averaged -> Lorentz LF data" << std::endl; + std::cout << " -o: flnOut = output file name." << std::endl; + std::cout << " -h: this help" << std::endl; + std::cout << std::endl; +} + +// ---------------------------------------------------------------------------- +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], + 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 << "# t (us), pol" << std::endl; + for (unsigned int i=0; i 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) { + std::cout << "**ERROR** in Gaussian LF calculation" << std::endl; + return 2; + } + } else { + std::cout << "not yet implemented." << std::endl; + } + + lf_gl_write(flnOut, field, param, tt, pol); + + return 0; +}