add LF Gaussian/Lorentzian testing code.
This commit is contained in:
parent
9fa90c24fb
commit
bcc1597e30
65
src/external/LF_GL/CMakeLists.txt
vendored
Normal file
65
src/external/LF_GL/CMakeLists.txt
vendored
Normal file
@ -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
|
||||
$<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}>
|
||||
)
|
||||
target_link_libraries(lf_gl -lm GSL::gsl)
|
||||
329
src/external/LF_GL/main.cpp
vendored
Normal file
329
src/external/LF_GL/main.cpp
vendored
Normal file
@ -0,0 +1,329 @@
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <string>
|
||||
#include <cstring>
|
||||
#include <vector>
|
||||
#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;
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
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<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],
|
||||
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 << "# t (us), pol" << std::endl;
|
||||
for (unsigned int i=0; i<tt.size(); i++) {
|
||||
fout << tt[i] << ", " << pol[i] << std::endl;
|
||||
}
|
||||
|
||||
fout.close();
|
||||
return 0;
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
double param[3] = {0.0, 0.0, 0.0};
|
||||
double field{0.0};
|
||||
bool gaussian_only = false;
|
||||
std::string flnOut="";
|
||||
|
||||
if (argc < 9) {
|
||||
lf_gl_syntax();
|
||||
return 0;
|
||||
}
|
||||
double dval;
|
||||
for (int i=1; i<argc; i++) {
|
||||
if (!strcmp(argv[i], "-p")) {
|
||||
try {
|
||||
dval = std::stod(argv[i+1]);
|
||||
} catch (...) {
|
||||
std::cout << std::endl;
|
||||
std::cout << "**ERROR** in handling of field." << std::endl;
|
||||
lf_gl_syntax();
|
||||
return 1;
|
||||
}
|
||||
param[0] = dval;
|
||||
try {
|
||||
dval = std::stod(argv[i+2]);
|
||||
} catch (...) {
|
||||
std::cout << std::endl;
|
||||
std::cout << "**ERROR** in handling of width." << std::endl;
|
||||
lf_gl_syntax();
|
||||
return 1;
|
||||
}
|
||||
param[1] = dval;
|
||||
try {
|
||||
dval = std::stod(argv[i+3]);
|
||||
} catch (...) {
|
||||
std::cout << std::endl;
|
||||
std::cout << "**ERROR** in handling of hopp." << std::endl;
|
||||
lf_gl_syntax();
|
||||
return 1;
|
||||
}
|
||||
param[2] = dval;
|
||||
i += 3;
|
||||
} else if (!strcmp(argv[i], "-g")) {
|
||||
if (!strcmp(argv[i+1], "1")) {
|
||||
gaussian_only = true;
|
||||
} else if (!strcmp(argv[i+1], "0")) {
|
||||
gaussian_only = false;
|
||||
} else {
|
||||
std::cout << std::endl;
|
||||
std::cout << "**ERROR** in handling option -g." << std::endl;
|
||||
lf_gl_syntax();
|
||||
return 1;
|
||||
}
|
||||
i++;
|
||||
} else if (!strcmp(argv[i], "-o")) {
|
||||
flnOut = std::string(argv[i+1]);
|
||||
}
|
||||
}
|
||||
|
||||
if (flnOut == "") {
|
||||
std::cout << std::endl;
|
||||
std::cout << "**ERROR** missing output file name." << std::endl;
|
||||
lf_gl_syntax();
|
||||
return 1;
|
||||
}
|
||||
|
||||
std::cout << "param: " << param[0] << ", " << param[1] << ", " << param[2] << std::endl;
|
||||
if (gaussian_only)
|
||||
std::cout << "Gaussian LF only" << std::endl;
|
||||
else
|
||||
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) {
|
||||
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;
|
||||
}
|
||||
Loading…
x
Reference in New Issue
Block a user