Commit fbbe2c3f authored by Jan Oliver Oelerich's avatar Jan Oliver Oelerich

[untested] converted everything from Angstroms to nm

parent 1b026e05
......@@ -64,9 +64,9 @@ void stemsalabim::Crystal::initFromXYZ(const std::string &field_file_content) {
algorithms::trim(tokens[4]);
algorithms::trim(tokens[8]);
_size_x = std::stod(tokens[0]) * 10;
_size_y = std::stod(tokens[4]) * 10;
_size_z = std::stod(tokens[8]) * 10;
_size_x = std::stod(tokens[0]);
_size_y = std::stod(tokens[4]);
_size_z = std::stod(tokens[8]);
} else {
output::error("Wrong syntax of extended XYZ file line 2.");
}
......@@ -77,9 +77,9 @@ void stemsalabim::Crystal::initFromXYZ(const std::string &field_file_content) {
algorithms::trim(tokens[1]);
algorithms::trim(tokens[2]);
_size_x = std::stod(tokens[0]) * 10;
_size_y = std::stod(tokens[1]) * 10;
_size_z = std::stod(tokens[2]) * 10;
_size_x = std::stod(tokens[0]);
_size_y = std::stod(tokens[1]);
_size_z = std::stod(tokens[2]);
}
......@@ -111,10 +111,10 @@ void stemsalabim::Crystal::initFromXYZ(const std::string &field_file_content) {
}
// x y z are columns 1,2,3 and the mean square displacement is column 4. Convert to angstrom
x = std::stod(tokens[1]) * 10;
y = std::stod(tokens[2]) * 10;
z = std::stod(tokens[3]) * 10;
msd = std::stod(tokens[4]) * 100;
x = std::stod(tokens[1]);
y = std::stod(tokens[2]);
z = std::stod(tokens[3]);
msd = std::stod(tokens[4]);
// set the site's occupant
if(with_slices) {
......@@ -145,45 +145,6 @@ void stemsalabim::Crystal::initFromXYZ(const std::string &field_file_content) {
}
void stemsalabim::Crystal::readFieldFile(const std::string &field_file_content) {
// some useful variable declarations
istringstream iss(field_file_content);
string line = "";
vector<string> tokens;
double x, y, z, field;
// ignore first and second lines.
getline(iss, line);
getline(iss, line);
// the rest of the lines is atoms with positions
while(getline(iss, line)) {
if(line.length() == 0)
continue;
tokens = algorithms::split(line, ' ');
if(tokens.size() != 4) {
output::error("Field file needs to have 4 columns!\n");
}
algorithms::trim(tokens[0]);
algorithms::trim(tokens[1]);
algorithms::trim(tokens[2]);
algorithms::trim(tokens[3]);
// x y z are columns 0,1,2 and the field is column 3. Convert to angstrom
x = std::stod(tokens[0]) * 10;
y = std::stod(tokens[1]) * 10;
z = std::stod(tokens[2]) * 10;
field = std::stod(tokens[3]);
// set the site's occupant
_fields.push_back(make_tuple(x, y, z, field));
}
}
void stemsalabim::Crystal::generateSlices() {
Params &p = Params::getInstance();
......
......@@ -58,8 +58,6 @@ namespace stemsalabim {
*/
void initFromXYZ(const std::string &field_file_content);
void readFieldFile(const std::string &field_file_content);
/*!
* Returns a reference to the vector of pointers to all Atom objects.
* @return atoms
......
......@@ -79,8 +79,7 @@ void IO::initNcFile(const shared_ptr<GridManager> &gridman, const shared_ptr<Cry
g.defineVar<float>("slice_potentials", vd({frameDim, slicesDim, gridxDim, gridyDim}));
const char angles_labels[3][6] = {"alpha", "beta", "gamma"};
vector<float> syslength = {(float) crystal->sizeX() / 10.0f, (float) crystal->sizeY() / 10.0f,
(float) crystal->sizeZ() / 10.0f};
vector<float> syslength = {(float) crystal->sizeX(), (float) crystal->sizeY(), (float) crystal->sizeZ()};
spatialVar.put<char>("xyz");
cellSpatialVar.put<char>("abc");
cellAngularVar.put<char>(angles_labels);
......@@ -272,32 +271,32 @@ void IO::writeParams(const shared_ptr<GridManager> &gridman) {
{
auto g = pg.defineGroup("probe");
g.att("c5", p.probeC5() / 10);
g.att("cs", p.probePhericalAberrationCoeff() / 10);
g.att("astigmatism_ca", p.probeAstigmatismCa() / 10);
g.att("defocus", p.meanDefocus() / 10);
g.att("fwhm_defoci", p.fwhmDefocus() / 10);
g.att("c5", p.probeC5());
g.att("cs", p.probePhericalAberrationCoeff());
g.att("astigmatism_ca", p.probeAstigmatismCa());
g.att("defocus", p.meanDefocus());
g.att("fwhm_defoci", p.fwhmDefocus());
g.att("num_defoci", p.numberOfDefoci());
g.att("astigmatism_ca", p.probeAstigmatismCa() / 10);
g.att("astigmatism_ca", p.probeAstigmatismCa());
g.att("astigmatism_angle", p.probeAstigmatismAngle());
g.att("min_apert", p.probeMinAperture());
g.att("max_apert", p.probeMaxAperture());
g.att("beam_energy", p.beamEnergy());
g.att("scan_density", p.probeDensity() * 10);
g.att("scan_density", p.probeDensity());
}
{
auto g = pg.defineGroup("specimen");
g.att("max_potential_radius", p.maxPotentialRadius() / 10);
g.att("max_potential_radius", p.maxPotentialRadius());
g.att("crystal_file", p.crystalFilename());
}
{
auto g = pg.defineGroup("grating");
g.att("density", p.samplingDensity() * 10);
g.att("density", p.samplingDensity());
g.att("nx", gridman->samplingX());
g.att("ny", gridman->samplingY());
g.att("slice_thickness", p.sliceThickness() / 10);
g.att("slice_thickness", p.sliceThickness());
}
{
......@@ -388,9 +387,9 @@ void IO::writeCrystal(const shared_ptr<FPConfManager> &fpman, const shared_ptr<C
size_t n = crystal->numberOfAtoms();
vector<float> cell_lengths(3);
cell_lengths[0] = (float) crystal->sizeX() / 10.0f;
cell_lengths[1] = (float) crystal->sizeY() / 10.0f;
cell_lengths[2] = (float) crystal->sizeZ() / 10.0f;
cell_lengths[0] = (float) crystal->sizeX();
cell_lengths[1] = (float) crystal->sizeY();
cell_lengths[2] = (float) crystal->sizeZ();
vector<short> elements_data(n);
vector<float> radii_data(n);
......@@ -407,18 +406,18 @@ void IO::writeCrystal(const shared_ptr<FPConfManager> &fpman, const shared_ptr<C
elements_data[i] = (short) distance(crystal->elements().begin(),
crystal->elements().find(atom->getElement()));
radii_data[i] = atom->getElement()->atomicRadius() / 10;
msd_data[i] = (float) (atom->getMSD() / 100.);
radii_data[i] = atom->getElement()->atomicRadius();
msd_data[i] = (float) (atom->getMSD());
double z = atom->getZ() + fpman->dz(i);
slice_data[i] = crystal->sliceIndex(z);
coords_data[3 * i + 0] = (float) (atom->getX() + fpman->dx(i)) / 10;
coords_data[3 * i + 1] = (float) (atom->getY() + fpman->dy(i)) / 10;
coords_data[3 * i + 2] = (float) (atom->getZ() + fpman->dz(i)) / 10;
coords_data[3 * i + 0] = (float) (atom->getX() + fpman->dx(i));
coords_data[3 * i + 1] = (float) (atom->getY() + fpman->dy(i));
coords_data[3 * i + 2] = (float) (atom->getZ() + fpman->dz(i));
lattice_coords_data[3 * i + 0] = (float) atom->getX() / 10;
lattice_coords_data[3 * i + 1] = (float) atom->getY() / 10;
lattice_coords_data[3 * i + 2] = (float) atom->getZ() / 10;
lattice_coords_data[3 * i + 0] = (float) atom->getX();
lattice_coords_data[3 * i + 1] = (float) atom->getY();
lattice_coords_data[3 * i + 2] = (float) atom->getZ();
i++;
}
......@@ -741,63 +740,27 @@ void IO::writeGrids(const shared_ptr<GridManager> &gridman) {
if(p.adf()) {
auto g = f.group("adf");
vector<double> probe_x_grid = gridman->adfXGrid();
vector<double> probe_y_grid = gridman->adfYGrid();
vector<double> slice_coords = gridman->adfSliceCoords();
vector<double> detector_grid = gridman->adfDetectorGrid();
for(size_t j = 0; j < probe_x_grid.size(); ++j)
probe_x_grid[j] /= 10;
for(size_t j = 0; j < probe_y_grid.size(); ++j)
probe_y_grid[j] /= 10;
for(size_t j = 0; j < slice_coords.size(); ++j)
slice_coords[j] /= 10;
g.var("adf_probe_x_grid").put(vs({0}), vs({probe_x_grid.size()}), probe_x_grid);
g.var("adf_probe_y_grid").put(vs({0}), vs({probe_y_grid.size()}), probe_y_grid);
g.var("adf_detector_grid").put(vs({0}), vs({detector_grid.size()}), detector_grid);
g.var("adf_slice_coords").put(vs({0}), vs({slice_coords.size()}), slice_coords);
g.var("adf_probe_x_grid").put(vs({0}), vs({gridman->adfXGrid().size()}), gridman->adfXGrid());
g.var("adf_probe_y_grid").put(vs({0}), vs({gridman->adfYGrid().size()}), gridman->adfYGrid());
g.var("adf_detector_grid").put(vs({0}), vs({gridman->adfDetectorGrid().size()}), gridman->adfDetectorGrid());
g.var("adf_slice_coords").put(vs({0}), vs({gridman->adfSliceCoords().size()}), gridman->adfSliceCoords());
}
if(p.cbed()) {
auto g = f.group("cbed");
vector<double> cbed_k_grid_x = gridman->outputKXGrid();
vector<double> cbed_k_grid_y = gridman->outputKYGrid();
vector<double> cbed_x_grid = gridman->cbedXGrid();
vector<double> cbed_y_grid = gridman->cbedYGrid();
vector<double> slice_cbed_coords = gridman->cbedSliceCoords();
for(size_t j = 0; j < cbed_x_grid.size(); ++j)
cbed_x_grid[j] /= 10;
for(size_t j = 0; j < cbed_y_grid.size(); ++j)
cbed_y_grid[j] /= 10;
for(size_t j = 0; j < slice_cbed_coords.size(); ++j)
slice_cbed_coords[j] /= 10;
g.var("cbed_x_grid").put(vs({0}), vs({cbed_k_grid_x.size()}), cbed_k_grid_x);
g.var("cbed_y_grid").put(vs({0}), vs({cbed_k_grid_y.size()}), cbed_k_grid_y);
g.var("cbed_probe_x_grid").put(vs({0}), vs({cbed_x_grid.size()}), cbed_x_grid);
g.var("cbed_probe_y_grid").put(vs({0}), vs({cbed_y_grid.size()}), cbed_y_grid);
g.var("cbed_slice_coords").put(vs({0}), vs({slice_cbed_coords.size()}), slice_cbed_coords);
g.var("cbed_x_grid").put(vs({0}), vs({gridman->outputKXGrid().size()}), gridman->outputKXGrid());
g.var("cbed_y_grid").put(vs({0}), vs({gridman->outputKYGrid().size()}), gridman->outputKYGrid());
g.var("cbed_probe_x_grid").put(vs({0}), vs({gridman->cbedXGrid().size()}), gridman->cbedXGrid());
g.var("cbed_probe_y_grid").put(vs({0}), vs({gridman->cbedYGrid().size()}), gridman->cbedYGrid());
g.var("cbed_slice_coords").put(vs({0}), vs({gridman->cbedSliceCoords().size()}), gridman->cbedSliceCoords());
}
{
auto g = f.group("params");
vector<double> defocus_grid = gridman->defoci();
vector<double> defocus_weights = gridman->defocusWeights();
for(size_t j = 0; j < defocus_grid.size(); ++j)
defocus_grid[j] /= 10;
g.var("defocus").put(vs({0}), vs({defocus_grid.size()}), defocus_grid);
g.var("defocus_weights").put(vs({0}), vs({defocus_weights.size()}), defocus_weights);
g.var("defocus").put(vs({0}), vs({gridman->defoci().size()}), gridman->defoci());
g.var("defocus_weights").put(vs({0}), vs({gridman->defocusWeights().size()}), gridman->defocusWeights());
}
} catch(NcException &e) {
......
......@@ -141,34 +141,6 @@ void Params::readParamsFromString(const string &prms) {
_adf_average_configurations = true;
}
// convert lengths to angstrom and densities to /angstrom
if(_cfg.exists("grating.density"))
_sample_density /= 10;
if(_cfg.exists("probe.scan_density"))
_probe_scan_density /= 10;
if(_cfg.exists("grating.slice_thickness"))
_slice_thickness *= 10;
if(_cfg.exists("specimen.max_potential_radius"))
_max_potential_radius *= 10;
if(_cfg.exists("probe.c5"))
_probe_c5_coeff *= 10;
if(_cfg.exists("probe.cs"))
_probe_spherical_aberration_coeff *= 10;
if(cliParamGiven("defocus") || _cfg.exists("probe.defocus"))
_probe_defocus *= 10;
if(_cfg.exists("probe.fwhm_defoci"))
_probe_fwhm_defoci *= 10;
if(_cfg.exists("probe.astigmatism_ca"))
_probe_astigmatism_ca *= 10;
if(_tmp_dir.empty()) {
// use the path of the output file directory for tmp
int found = _output_file.find_last_of("/\\");
......@@ -396,15 +368,4 @@ void Params::readParamsFromNCFile() {
}
// convert lengths to angstrom and densities to /angstrom
_sample_density /= 10;
_probe_scan_density /= 10;
_slice_thickness *= 10;
_max_potential_radius *= 10;
_probe_c5_coeff *= 10;
_probe_spherical_aberration_coeff *= 10;
_probe_defocus *= 10;
_probe_fwhm_defoci *= 10;
_probe_astigmatism_ca *= 10;
}
......@@ -529,8 +529,8 @@ namespace stemsalabim {
/// electron rest mass in keV. \f$m_0 c^2\f$ in keV
constexpr static double ELECTRON_REST_MASS = 510.99906;
/// product of planck's constant and speed of light in keV * Angstroms.
constexpr static double PLANCK_CONSTANT_SPEED_LIGHT = 12.39841875;
/// product of planck's constant and speed of light in keV * nm.
constexpr static double PLANCK_CONSTANT_SPEED_LIGHT = 1.239841875;
/// pi. Approximately.
constexpr static double pi = 3.14159265358979323846;
......
......@@ -172,8 +172,8 @@ void Simulation::run() {
for(shared_ptr<Slice> &slice: _crystal->slices())
printMaster(output::fmt("Slice %d [%.3fnm - %.3fnm] contains %d atoms.\n",
slice->id(),
slice->z() / 10,
(slice->z() + slice->thickness()) / 10,
slice->z(),
(slice->z() + slice->thickness()),
slice->numberOfAtoms()));
}
......@@ -189,7 +189,7 @@ void Simulation::run() {
// if a new defocus was started, do some output stuff.
if(!st.iconf()) {
printMaster(st, output::fmt("Starting Defocus %.2fnm.\n", st.defocus() / 10));
printMaster(st, output::fmt("Starting Defocus %.2fnm.\n", st.defocus()));
}
printMaster(st, "Starting configuration. \n");
......@@ -206,8 +206,8 @@ void Simulation::run() {
printMaster(st,
output::fmt("Slice %d [%.3fnm - %.3fnm] contains %d atoms.\n",
slice->id(),
slice->z() / 10,
(slice->z() + slice->thickness()) / 10,
slice->z(),
(slice->z() + slice->thickness()),
slice->numberOfAtoms()));
}
......
......@@ -246,22 +246,15 @@ namespace stemsalabim { namespace atomic {
val.get("symbol", "").asString(),
val.get("atomic_number", 0).asUInt()));
e->setLatticeStructure(val.get("lattice_structure", "").asString());
e->setAtomicWeight(val.get("atomic_weight", 0.f).asFloat());
try {
// convert pm to angstrom
e->setAtomicRadius(val.get("atomic_radius pm", 0.f).asFloat() / 100);
e->setAtomicRadius(val.get("atomic_radius", 0.f).asFloat());
} catch(const Json::LogicError &err) {
e->setAtomicRadius(0.0);
}
try {
e->setLatticeConstant(val.get("lattice_constant ang", 0.f).asFloat());
} catch(const Json::LogicError &err) {
e->setLatticeConstant(0.0);
}
std::vector<double> scattering_params_dt(12);
scattering_params_dt[0] = val["scattering_dt"].get("a1", 0).asDouble();
......
......@@ -35,66 +35,6 @@ namespace stemsalabim { namespace atomic {
* @tparam float precision type, float or double
*/
/*!
* regular modified cylindrical Bessel function order 0
* @param x argument
* @return regular modified cylindrical Bessel function order 0 of x
*/
static double bessi0(double x) {
double ax, ans;
double y;
if((ax = fabs(x)) < 3.75) {
y = x / 3.75, y = y * y;
ans = 1.0 +
y *
(3.5156229 +
y * (3.0899424 + y * (1.2067492 + y * (0.2659732 + y * (0.360768e-1 + y * 0.45813e-2)))));
} else {
y = 3.75 / ax;
ans = (exp(ax) / sqrt(ax)) *
(0.39894228 +
y *
(0.1328592e-1 +
y *
(0.225319e-2 +
y *
(-0.157565e-2 +
y *
(0.916281e-2 +
y * (-0.2057706e-1 + y * (0.2635537e-1 + y * (-0.1647633e-1 + y * 0.392377e-2))))))));
}
return ans;
}
/*!
* irregular modified cylindrical Bessel function order 0
* @param x argument
* @return irregular modified cylindrical Bessel function order 0 of x
*/
static double bessk0(double x) {
double y, ans;
if(x <= 2.0) {
y = x * x / 4.0;
ans = (-log(x / 2.0) * bessi0(x)) +
(-0.57721566 +
y *
(0.42278420 +
y * (0.23069756 + y * (0.3488590e-1 + y * (0.262698e-2 + y * (0.10750e-3 + y * 0.74e-5))))));
} else {
y = 2.0 / x;
ans = (exp(-x) / sqrt(x)) *
(1.25331414 +
y *
(-0.7832358e-1 +
y *
(0.2189568e-1 + y * (-0.1062446e-1 + y * (0.587872e-2 + y * (-0.251540e-2 + y * 0.53208e-3))))));
}
return ans;
}
class Scattering {
public:
......@@ -167,9 +107,9 @@ namespace stemsalabim { namespace atomic {
double r2 = r * r;
double f = pi2 * r2;
double sum = 2 *
(sp[0] * bessk0(2 * pi * r * sqrt(sp[1])) +
sp[2] * bessk0(2 * pi * r * sqrt(sp[3])) +
sp[4] * bessk0(2 * pi * r * sqrt(sp[5]))) +
(sp[0] * std::cyl_bessel_k(0, 2 * pi * r * sqrt(sp[1])) +
sp[2] * std::cyl_bessel_k(0, 2 * pi * r * sqrt(sp[3])) +
sp[4] * std::cyl_bessel_k(0, 2 * pi * r * sqrt(sp[5]))) +
sp[6] / sp[7] * exp(-f / sp[7]) +
sp[8] / sp[9] * exp(-f / sp[9]) +
sp[10] / sp[11] * exp(-f / sp[11]);
......@@ -207,7 +147,7 @@ namespace stemsalabim { namespace atomic {
std::map<std::string, std::vector<double>> _cache;
/// this is \f$2\pi^2 a_0 e\f$, with \f$e = 2R_ya_0\f$ (\f$R_y = 13.6 eV\f$)
constexpr static double PREFACTOR_ATOMPOT = 150.4120584;
constexpr static double PREFACTOR_ATOMPOT = 1.504120584;
/// pi
constexpr static double pi = 3.14159265358979323846;
};
......
This diff is collapsed.
......@@ -160,8 +160,8 @@ int main(int argc, const char **argv) {
Simulation::printMaster(st,
output::fmt("Slice %d [%.3fnm - %.3fnm] contains %d atoms.\n",
slice->id(),
slice->z() / 10,
(slice->z() + slice->thickness()) / 10,
slice->z(),
(slice->z() + slice->thickness()),
slice->numberOfAtoms()));
_io->writeCrystal(_fpman, _crystal);
......
......@@ -20,7 +20,7 @@
# compilation
add_executable(stemsalabim_test main.cpp mpi.cpp queue.cpp utilities.cpp fft.cpp buffer.cpp)
add_executable(stemsalabim_test main.cpp mpi.cpp queue.cpp utilities.cpp fft.cpp buffer.cpp math.cpp)
# linking ################################################
......
/*
* STEMsalabim: Magical STEM image simulations
*
* Copyright (c) 2016-2018 Jan Oliver Oelerich <jan.oliver.oelerich@physik.uni-marburg.de>
* Copyright (c) 2016-2018 Structure and Technology Research Laboratory, Philipps-Universität Marburg, Germany
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
*/
#include "catch.hpp"
#include <complex>
#include <cmath>
static double bessi0(double x) {
double ax, ans;
double y;
if((ax = fabs(x)) < 3.75) {
y = x / 3.75, y = y * y;
ans = 1.0 +
y *
(3.5156229 +
y * (3.0899424 + y * (1.2067492 + y * (0.2659732 + y * (0.360768e-1 + y * 0.45813e-2)))));
} else {
y = 3.75 / ax;
ans = (exp(ax) / sqrt(ax)) *
(0.39894228 +
y *
(0.1328592e-1 +
y *
(0.225319e-2 +
y *
(-0.157565e-2 +
y *
(0.916281e-2 +
y * (-0.2057706e-1 + y * (0.2635537e-1 + y * (-0.1647633e-1 + y * 0.392377e-2))))))));
}
return ans;
}
/*!
* irregular modified cylindrical Bessel function order 0
* @param x argument
* @return irregular modified cylindrical Bessel function order 0 of x
*/
static double bessk0(double x) {
double y, ans;
if(x <= 2.0) {
y = x * x / 4.0;
ans = (-log(x / 2.0) * bessi0(x)) +
(-0.57721566 +
y *
(0.42278420 +
y * (0.23069756 + y * (0.3488590e-1 + y * (0.262698e-2 + y * (0.10750e-3 + y * 0.74e-5))))));
} else {
y = 2.0 / x;
ans = (exp(-x) / sqrt(x)) *
(1.25331414 +
y *
(-0.7832358e-1 +
y *
(0.2189568e-1 + y * (-0.1062446e-1 + y * (0.587872e-2 + y * (-0.251540e-2 + y * 0.53208e-3))))));
}
return ans;
}
TEST_CASE("Bessel functions", "[math]") {
REQUIRE(std::cyl_bessel_i(0, 2.0) == Approx(bessi0(2.0)));
REQUIRE(std::cyl_bessel_k(0, 2.0) == Approx(bessk0(2.0)));
}
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