C function for energy distribution spectral density

This commit is contained in:
Dominik Demuth 2023-04-17 20:17:54 +02:00
parent 4a9d50e8a4
commit 2042148d0f
2 changed files with 60 additions and 12 deletions

View File

@ -1,6 +1,8 @@
/* integrands used in quadrature integration with scipy's LowLevelCallables */
#include <math.h>
#define KB 8.617333262145179e-05
/* FFHS functions */
double ffhsSD(double x, void *user_data) {
double *c = (double *)user_data;
@ -18,7 +20,7 @@ double ffhsSD(double x, void *user_data) {
/* log-gaussian functions */
double logNormalDist(double tau, double tau0, double sigma) {
return exp(- pow((log(tau/tau0) / sigma), 2) / 2) / sqrt(2*M_PI)/sigma;
return exp(- pow((log(tau/tau0) / sigma), 2) / 2.) / sqrt(2*M_PI)/sigma;
}
double logGaussianSD_high(double u, void *user_data) {
@ -61,3 +63,26 @@ double logGaussianCorrelation(double x, void *user_data) {
return dist * exp(-t/uu);
}
double normalDist(double x, double x0, double sigma) {
return exp(- pow((x-x0) / sigma, 2) / 2.) / sqrt(2 * M_PI) / sigma;
}
double rate(double tau0, double ea, double t) {
return exp(-ea / t / KB) / tau0;
}
double energyDist_SD(double x, void *user_data) {
double *c = (double *)user_data;
double omega = c[0];
double tau0 = c[1];
double e_m = c[2];
double e_b = c[3];
double temp = c[4];
double r = rate(tau0, x, temp);
return r/(pow(r, 2) + pow(omega, 2)) * normalDist(x, e_m, e_b);
}

View File

@ -1,13 +1,18 @@
import ctypes
from itertools import product
import numpy as np
from scipy import LowLevelCallable
from scipy.integrate import quad, simps as simpson
from .base import Distribution
from ..lib.utils import ArrayLike
from ..utils.constants import kB
from .helper import HAS_C_FUNCS, lib
# noinspection PyMethodOverriding
class EnergyBarriers(Distribution):
name = 'Energy barriers'
parameter = [r'\tau_{0}', r'E_{m}', r'\Delta E']
@ -51,7 +56,7 @@ class EnergyBarriers(Distribution):
omega = np.atleast_1d(omega)
temperature = np.atleast_1d(temperature)
e_axis = np.linspace(max(0, e_m-50*e_b), e_m+50*e_b, num=5001)
e_axis = np.linspace(max(0., e_m-50*e_b), e_m+50*e_b, num=5001)
ret_val = []
for o, tt in product(omega, temperature):
ret_val.append(simpson(_integrand_freq_real(e_axis, o, tau0, e_m, e_b, tt), e_axis) -
@ -60,23 +65,41 @@ class EnergyBarriers(Distribution):
return np.array(ret_val)
@staticmethod
def specdens(omega, temperature, *args):
def specdens(omega, temperature, tau0: float, e_m: float, e_b: float):
# in contrast to other spectral densities, it's omega and temperature
tau0, e_m, e_b = args
def integrand(e_a, w, t0, mu, sigma, t):
r = EnergyBarriers.rate(t0, e_a, t)
return r/(r**2 + w**2) * EnergyBarriers.energydistribution(e_a, mu, sigma)
omega = np.atleast_1d(omega)
temperature = np.atleast_1d(temperature)
e_axis = np.linspace(max(0, e_m-50*e_b), e_m+50*e_b, num=5001)
if HAS_C_FUNCS:
ret_val = EnergyBarriers.spec_dens_c(omega, temperature, tau0, e_m, e_b)
else:
ret_val = EnergyBarriers.spec_dens_py(omega, temperature, tau0, e_m, e_b)
ret_val = np.array([simpson(integrand(e_axis, o, tau0, e_m, e_b, tt), e_axis)
for o in omega for tt in temperature])
return ret_val.squeeze()
return ret_val
@staticmethod
def spec_dens_c(omega: np.ndarray, temperature: np.ndarray, tau0: float, e_m: float, e_b: float) -> np.ndarray:
res = []
for o, t in product(omega, temperature):
c = (ctypes.c_double * 5)(o, tau0, e_m, e_b, t)
user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p)
area = quad(LowLevelCallable(lib.energyDist_SD, user_data), 0, np.infty, epsabs=1e-10)[0]
res.append(area)
return np.array(res)
@staticmethod
def spec_dens_py(omega: np.ndarray, temperature: np.ndarray, tau0: float, e_m: float, e_b: float) -> np.ndarray:
def integrand(e_a, w, t0, mu, sigma, t):
r = EnergyBarriers.rate(t0, e_a, t)
return r/(r**2 + w**2) * EnergyBarriers.energydistribution(e_a, mu, sigma)
e_axis = np.linspace(max(0., e_m-50*e_b), e_m+50*e_b, num=5001)
ret_val = [simpson(integrand(e_axis, o, tau0, e_m, e_b, tt), e_axis) for o in omega for tt in temperature]
return np.array(ret_val)
@staticmethod
def mean(*args):