Commit 8d4c40c6 authored by Jan Oliver Oelerich's avatar Jan Oliver Oelerich

Some improvements, fixes, and integration of GSL.

parent 06529fb5
......@@ -70,6 +70,11 @@ find_package(LibConfig REQUIRED)
include_directories(${LIBCONFIG_INCLUDE_DIR})
set(LIBS ${LIBS} ${LIBCONFIG_LIBRARIES})
# look for GSL
find_package(GSL REQUIRED)
include_directories(${GSL_INCLUDE_DIRS})
set(LIBS ${LIBS} ${GSL_LIBRARIES})
# look for NetCDF
find_package(NetCDF REQUIRED)
include_directories(${NETCDF_INCLUDE_DIR})
......
......@@ -204,9 +204,9 @@ void GridManager::generateGrids() {
_energy_loss_grid[1] = p.plasmonsMaxEnergy();
} else if(p.plasmonsEnabled()) {
auto num_plasmon_energy_bins = (unsigned int) ceil(p.plasmonsMaxEnergy() * p.plasmonsEnergyGridDensity()) + 1;
_energy_loss_grid.resize(num_plasmon_energy_bins,0);
for(unsigned int ibin = 0; ibin < num_plasmon_energy_bins; ++ibin) {
_energy_loss_grid[ibin] = ibin / (double) num_plasmon_energy_bins * p.plasmonsMaxEnergy();
_energy_loss_grid[ibin] = (float)ibin / num_plasmon_energy_bins * p.plasmonsMaxEnergy();
}
} else {
_energy_loss_grid.resize(1,0);
......
......@@ -482,7 +482,7 @@ namespace stemsalabim {
* energy loss bins, i.e., 0-th entry is energy loss between 0 and the energy of index 1.
* @return vector with energy loss bins
*/
const std::vector<double> &energyLossGrid() const {
const std::vector<float> &energyLossGrid() const {
return _energy_loss_grid;
}
......@@ -588,7 +588,7 @@ namespace stemsalabim {
std::vector<double> _cbed_x_grid;
std::vector<double> _cbed_y_grid;
std::vector<double> _energy_loss_grid;
std::vector<float> _energy_loss_grid;
unsigned int _sampling_x;
unsigned int _sampling_y;
......
......@@ -26,6 +26,10 @@
#include <mutex>
#include <algorithm>
#include <cmath>
#include <gsl/gsl_integration.h>
#include "../utilities/algorithms.hpp"
#include "../classes/GridManager.hpp"
#include "../classes/IO.hpp"
......@@ -34,13 +38,12 @@ namespace stemsalabim { namespace atomic {
class Plasmons {
public:
static Plasmons &
getInstance() {
static Plasmons &getInstance() {
static Plasmons instance;
return instance;
}
void populateCache(const std::shared_ptr<GridManager> & gridman, double electron_density, double atom_density) {
void populateCache(const std::shared_ptr<GridManager> &gridman, double electron_density, double atom_density) {
if(_cache_populated)
return;
......@@ -54,13 +57,13 @@ namespace stemsalabim { namespace atomic {
unsigned int ly = gridman->samplingY();
double sdense2_x = pow(1. / 3. * gridman->densityX(), 2);
double sdense2_y = pow(1. / 3. * gridman->densityY(), 2);
double energy_delta = gridman->energyLossGrid()[1];
int fine_grid_num = 10;
if(p.plasmonsSimple())
fine_grid_num = 1000;
gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000000);
double error, result;
gsl_function F;
F.function = algorithms::unwrap;
_cache.resize(gridman->energyLossGrid().size()-1);
_cache.resize(gridman->energyLossGrid().size() - 1);
//unsigned int energy_grid_index = 0;
for(unsigned int eii = 0; eii < _cache.size(); ++eii) {
......@@ -72,33 +75,42 @@ namespace stemsalabim { namespace atomic {
for(unsigned int iy = 0; iy < ly; iy++) {
double theta = sqrt(pow(gridman->kx(ix), 2) + pow(gridman->ky(iy), 2)) * p.wavelength();
if(!p.bandwidthLimiting() || pow(gridman->kx(ix), 2) / sdense2_x + pow(gridman->ky(iy), 2) / sdense2_y < 1) {
if(!p.bandwidthLimiting() ||
pow(gridman->kx(ix), 2) / sdense2_x + pow(gridman->ky(iy), 2) / sdense2_y < 1) {
// integrate over energy
// @TODO: Replace with integration quad library
double integral = 0;
for(int fine_grid_i = 1; fine_grid_i <= fine_grid_num; fine_grid_i++) {
double energy = gridman->energyLossGrid()[eii] + (float)fine_grid_i / fine_grid_num * energy_delta;
integral += plasmonDoubleDifferentialCrossSection(energy, theta, atom_density, electron_density) * energy_delta / fine_grid_num;
}
_cache[eii](ix, iy) = (float)sqrt(integral);
std::function<double(double)> func = [&](double x) {
return plasmonDoubleDifferentialCrossSection(x, theta, atom_density, electron_density);
};
F.params = &func;
gsl_integration_qag(&F,
gridman->energyLossGrid()[eii],
gridman->energyLossGrid()[eii + 1],
0,
1e-7,
1000000,
GSL_INTEG_GAUSS61,
w,
&result,
&error);
_cache[eii](ix, iy) = (float) sqrt(result);
}
}
}
_cache[eii](0, 0) = _cache[eii](0, 1);
// IO::dumpWave(output::fmt("w%d_k.dat", eii), _cache[eii]);
// printf("wl %e\n", gridman->kx(1) * p.wavelength());
// printf("wl %e\n", gridman->kx(lx/2) * p.wavelength());
// printf("wl %e\n", p.wavelength());
_cache[eii].backwardFFT();
_cache[eii].normalize();
// IO::dumpWave(output::fmt("w%d_r.dat", eii), _cache[eii]);
// exit(0);
_cache[eii].backwardFFT();
}
gsl_integration_workspace_free(w);
_cache_populated = true;
}
const Wave & scatteringFunction(unsigned int energy_loss_bin_index) const {
const Wave &scatteringFunction(unsigned int energy_loss_bin_index) const {
return _cache[energy_loss_bin_index];
}
......@@ -129,7 +141,8 @@ namespace stemsalabim { namespace atomic {
* @param theta scattering angle in rad
* @return double differential cross section in angstroms^2/eV
*/
double plasmonDoubleDifferentialCrossSection(double E, double theta, double atom_density, double electron_density) {
double
plasmonDoubleDifferentialCrossSection(double E, double theta, double atom_density, double electron_density) {
// calculate fermi energy in eV
// this factor is pow(3*pi^2, 2./3.)
......@@ -152,14 +165,16 @@ namespace stemsalabim { namespace atomic {
double factor = 2 * _incident_energy * 1000 * (theta * theta + characteristic_angle * characteristic_angle);
double g = 3. / 5. * fermi_energy / _plasmon_energy;
// plasmon energy dispersion
double Ep = _plasmon_energy + g * factor;
// double differential cross section of plasmon scattering in angstroms^2/eV
double B = 1 / (atom_density * a0 * 2 * pi * pi);
B = sin(theta) / (atom_density * a0 * pi);
B *= 2 / factor;
B *= E * _plasmon_fwhm * _plasmon_energy * _plasmon_energy / (pow(E * E - _plasmon_energy * _plasmon_energy - 2*g*factor*_plasmon_energy, 2.) + pow(E * _plasmon_fwhm, 2.));
B *= E *
_plasmon_fwhm *
_plasmon_energy *
_plasmon_energy /
(pow(E * E - _plasmon_energy * _plasmon_energy - 2 * g * factor * _plasmon_energy, 2.) +
pow(E * _plasmon_fwhm, 2.));
return B;
}
......
......@@ -30,6 +30,7 @@
#include <string>
#include <random>
#include <cstring>
#include <functional>
#include <unistd.h>
#include <iostream>
#include <sstream>
......@@ -400,6 +401,12 @@ namespace stemsalabim {
return N;
}
inline double unwrap(double x, void *p) {
auto fp = static_cast<std::function<double(double)> *>(p);
return (*fp)(x);
}
}
}
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment