Commit 88c481f5 authored by Jan Oliver Oelerich's avatar Jan Oliver Oelerich

Normalized plasmon scattering.

parent 8d4c40c6
......@@ -67,7 +67,8 @@ cbed: {
plasmon_scattering: {
enabled = false; # enable plasmon scattering
simple_mode = true; # No energy resolution, only E = 0 and E > 0
simple_mode = true; # No energy resolution, only E = 0 and 0 < E < max_energy
plural_scattering = true; # Include plural scattering. Makes simulation pretty slow if simple_mode = false
max_energy = 10; # max energy of the energy grid considered for plasmon energy transfer in eV
energy_grid_density = 10; # density of the energy grid in 1/eV
mean_free_path = 120; # mean free path of a plasmon in nm
......
......@@ -368,6 +368,7 @@ void IO::writeParams(const shared_ptr<GridManager> &gridman, const string &conf_
auto g = pg.defineGroup("plasmon_scattering");
g.att("enabled", p.plasmonsEnabled());
g.att("simple_mode", p.plasmonsSimple());
g.att("plural_scattering", p.plasmonsPluralScattering());
g.att("max_energy", p.plasmonsMaxEnergy());
g.att("energy_grid_density", p.plasmonsEnergyGridDensity());
g.att("mean_free_path", p.plasmonMeanFreePath());
......
......@@ -138,6 +138,7 @@ void Params::readParamsFromString(const string &prms) {
if(_cfg.exists("plasmon_scattering")) {
_cfg.lookupValue("plasmon_scattering.enabled", _plasmons_enabled);
_cfg.lookupValue("plasmon_scattering.simple_mode", _plasmons_simple);
_cfg.lookupValue("plasmon_scattering.plural_scattering", _plasmons_plural);
_cfg.lookupValue("plasmon_scattering.max_energy", _plasmons_max_energy);
_cfg.lookupValue("plasmon_scattering.energy_grid_density", _plasmons_energy_grid_density);
_cfg.lookupValue("plasmon_scattering.mean_free_path", _plasmons_mean_free_path);
......@@ -257,6 +258,7 @@ void Params::broadcast(int rank) {
mpi_env.broadcast(_plasmons_enabled, rank);
mpi_env.broadcast(_plasmons_simple, rank);
mpi_env.broadcast(_plasmons_plural, rank);
mpi_env.broadcast(_plasmons_max_energy, rank);
mpi_env.broadcast(_plasmons_energy_grid_density, rank);
mpi_env.broadcast(_plasmons_mean_free_path, rank);
......@@ -395,6 +397,7 @@ void Params::readParamsFromNCFile(const std::string & path) {
_plasmons_enabled = gg.is("enabled");
_plasmons_simple = gg.is("simple_mode");
_plasmons_plural = gg.is("plural_scattering");
_plasmons_max_energy = gg.att<float>("max_energy");
_plasmons_energy_grid_density = gg.att<float>("energy_grid_density");
_plasmons_mean_free_path = gg.att<float>("mean_free_path");
......
......@@ -512,7 +512,6 @@ namespace stemsalabim {
/*!
* Return true if plasmon scattering is enabled.
* @return true if enabled
*/
bool plasmonsEnabled() const {
return _plasmons_enabled;
......@@ -520,15 +519,20 @@ namespace stemsalabim {
/*!
* Return true when the simple plasmon mode is activated
* @return true if enabled
*/
bool plasmonsSimple() const {
return _plasmons_simple;
}
/*!
* Return true when plural plasmon scattering is enabled
*/
bool plasmonsPluralScattering() const {
return _plasmons_plural;
}
/*!
* Return the maximal energy for the plasmon energy grid
* @return plasmon maximal energy
*/
float plasmonsMaxEnergy() const {
return _plasmons_max_energy;
......@@ -536,7 +540,6 @@ namespace stemsalabim {
/*!
* Return the density in 1/eV for the plasmon energy grid.
* @return density in 1/eV
*/
float plasmonsEnergyGridDensity() const {
return _plasmons_energy_grid_density;
......@@ -544,7 +547,6 @@ namespace stemsalabim {
/*!
* Plasmon mean free path in nm
* @return mean free path in nm
*/
float plasmonMeanFreePath() const {
return _plasmons_mean_free_path;
......@@ -552,7 +554,6 @@ namespace stemsalabim {
/*!
* Plasmon energy in eV
* @return plasmon energy in eV
*/
float plasmonEnergy() const {
return _plasmons_energy;
......@@ -560,7 +561,6 @@ namespace stemsalabim {
/*!
* Plasmon energy full with half maximum in eV
* @return plasmon energy FWHM in eV
*/
float plasmonFWHM() const {
return _plasmons_fwhm;
......@@ -649,6 +649,7 @@ namespace stemsalabim {
// plasmon scattering
bool _plasmons_enabled{false};
bool _plasmons_simple{false};
bool _plasmons_plural{true};
float _plasmons_max_energy{10};
float _plasmons_energy_grid_density{10};
float _plasmons_mean_free_path{110};
......
......@@ -322,18 +322,21 @@ void Simulation::calculatePixel(ScanPoint &point, const double defocus) {
if(p.plasmonsEnabled()) {
// scattering
for(unsigned int energy_i = 1; energy_i < _gridman->energyLossGrid().size(); energy_i++) {
for(unsigned long energy_i = _gridman->energyLossGrid().size()-1ul; energy_i > 0; energy_i--) {
waves[energy_i] *= sqrt(plasmons.beta0());
// plural scattering includes single scattering with energy_j = 0
for(unsigned int energy_j = 0; energy_j < energy_i; energy_j++) {
for(unsigned int energy_j = 0; energy_j < (p.plasmonsPluralScattering() ? energy_i : 1u); energy_j++) {
tmpwave = waves[energy_j];
tmpwave *= plasmons.scatteringFunction(energy_i - energy_j - 1);
tmpwave *= sqrt(plasmons.beta1());
waves[energy_i] += tmpwave;
}
}
waves[0] *= sqrt(plasmons.beta0());
}
for(auto &c: waves)
......@@ -347,10 +350,21 @@ void Simulation::calculatePixel(ScanPoint &point, const double defocus) {
for(auto &c: waves)
propagate(c, false);
//wave.forwardFFT();
// we are now in fourier space, so let's normalize the inelastic scattered wave functions
if(p.plasmonsEnabled()) {
if(!p.plasmonsEnabled() && p.normalizeAlways())
double total_intensity = 0;
double desired_intensity = 1-waves[0].integratedIntensity();
for(unsigned long ei = 1; ei < _gridman->energyLossGrid().size(); ei++)
total_intensity += waves[ei].integratedIntensity();
for(unsigned long ei = 1; ei < _gridman->energyLossGrid().size(); ei++)
waves[ei] *= sqrt(desired_intensity / total_intensity);
} else if(p.normalizeAlways()) {
waves[0].normalize();
}
if(point.adf && _gridman->adfStoreSlice(slice->id())) {
......
......@@ -100,7 +100,7 @@ namespace stemsalabim { namespace atomic {
}
}
}
_cache[eii](0, 0) = _cache[eii](0, 1);
//_cache[eii](0, 0) = _cache[eii](0, 1);
_cache[eii].normalize();
_cache[eii].backwardFFT();
}
......
......@@ -406,10 +406,21 @@ namespace stemsalabim {
/*!
* Get the square sum of all data.
*/
double integratedIntensity() const {
double integratedIntensity() {
bool transformed = false;
if(!_is_kspace) {
forwardFFT();
transformed = true;
}
double sum = 0.0;
for(unsigned int i = 0; i < _lx * _ly; i++)
sum += pow(std::abs(_data[i]), 2);
if(transformed)
backwardFFT();
return sum;
}
......
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