Newer
Older
Gigg, Martyn Anthony
committed
//------------------------------------------------------------------------------
// Includes
//------------------------------------------------------------------------------
#include "MantidKernel/Atom.h"
#include "MantidKernel/StringTokenizer.h"
#include <boost/lexical_cast.hpp>
#include <boost/make_shared.hpp>
#include <boost/regex.hpp>
namespace Mantid {
namespace Kernel {
typedef Mantid::Kernel::StringTokenizer tokenizer;
typedef std::pair<std::string, std::string> str_pair;
using PhysicalConstants::Atom;
using PhysicalConstants::getAtom;
using PhysicalConstants::NeutronAtom;
/**
* Construct an "empty" material. Everything returns zero
*/
Material::Material()
: m_name(), m_chemicalFormula(), m_atomTotal(0.0),
m_numberDensity(0.0), m_temperature(0.0), m_pressure(0.0) {}
/**
* Construct a material object
* @param name :: The name of the material
* @param formula :: The chemical formula
* @param numberDensity :: Density in A^-3
* @param temperature :: The temperature in Kelvin (Default = 300K)
* @param pressure :: Pressure in kPa (Default: 101.325 kPa)
*/
Material::Material(const std::string &name,
const ChemicalFormula &formula,
const double numberDensity, const double temperature,
const double pressure)
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
: m_name(name), m_atomTotal(0.0), m_numberDensity(numberDensity),
m_temperature(temperature), m_pressure(pressure) {
m_chemicalFormula.assign(formula.begin(), formula.end());
this->countAtoms();
}
/**
* Construct a material object
* @param name :: The name of the material
* @param formula :: The chemical formula
* @param numberDensity :: Density in A^-3
* @param temperature :: The temperature in Kelvin (Default = 300K)
* @param pressure :: Pressure in kPa (Default: 101.325 kPa)
*/
Material::Material(const std::string &name,
const PhysicalConstants::NeutronAtom &atom,
const double numberDensity, const double temperature,
const double pressure)
: m_name(name), m_chemicalFormula(), m_atomTotal(1.0), m_numberDensity(numberDensity),
m_temperature(temperature), m_pressure(pressure) {
if (atom.a_number == 0) { // user specified atom
FormulaUnit unit {boost::make_shared<Atom>(atom), 1.};
m_chemicalFormula.push_back(unit);
} else if (atom.z_number > 0) { // single isotope
FormulaUnit unit {boost::make_shared<Atom>(getAtom(atom.z_number, atom.a_number)), 1.};
m_chemicalFormula.push_back(unit);
} else { // isotopic average
FormulaUnit unit {boost::make_shared<Atom>(atom), 1.};
m_chemicalFormula.push_back(unit);
}
}
// update the total atom count
void Material::countAtoms() {
m_atomTotal = std::accumulate(std::begin(m_chemicalFormula),
std::end(m_chemicalFormula),
0.,
[](double subtotal, const FormulaUnit &right) {
return subtotal + right.multiplicity;
});
}
/**
* Returns the name
* @returns A string containing the name of the material
*/
const std::string &Material::name() const { return m_name; }
const Material::ChemicalFormula &Material::chemicalFormula() const { return m_chemicalFormula; }
/**
* Get the number density
* @returns The number density of the material in A^-3
*/
double Material::numberDensity() const { return m_numberDensity; }
/**
* Get the temperature
* @returns The temperature of the material in Kelvin
*/
double Material::temperature() const { return m_temperature; }
/**
* Get the pressure
* @returns The pressure of the material
*/
double Material::pressure() const { return m_pressure; }
/**
* Get the coherent scattering cross section for a given wavelength.
* CURRENTLY this simply returns the value for the underlying element
* @param lambda :: The wavelength to evaluate the cross section
* @returns The value of the coherent scattering cross section at
* the given wavelength
*/
double Material::cohScatterXSection(const double lambda) const {
UNUSED_ARG(lambda);
const double weightedTotal = std::accumulate(std::begin(m_chemicalFormula),
std::end(m_chemicalFormula),
0.,
[](double subtotal, const FormulaUnit &right) {
return subtotal + right.atom->neutron.coh_scatt_xs * right.multiplicity;
}) / m_atomTotal;
if (! std::isnormal(weightedTotal)) {
return 0.;
} else {
return weightedTotal;
}
}
/**
* Get the incoherent scattering cross section for a given wavelength
* CURRENTLY this simply returns the value for the underlying element
* @param lambda :: The wavelength to evaluate the cross section
* @returns The value of the coherent scattering cross section at
* the given wavelength
*/
double Material::incohScatterXSection(const double lambda) const {
UNUSED_ARG(lambda);
const double weightedTotal = std::accumulate(std::begin(m_chemicalFormula),
std::end(m_chemicalFormula),
0.,
[](double subtotal, const FormulaUnit &right) {
return subtotal + right.atom->neutron.inc_scatt_xs * right.multiplicity;
}) / m_atomTotal;
if (! std::isnormal(weightedTotal)) {
return 0.;
} else {
return weightedTotal;
}
}
/**
* Get the total scattering cross section for a given wavelength
* CURRENTLY this simply returns the value for sum of the incoherent
* and coherent scattering cross sections.
* @param lambda :: The wavelength to evaluate the cross section
* @returns The value of the total scattering cross section at
* the given wavelength
*/
double Material::totalScatterXSection(const double lambda) const {
UNUSED_ARG(lambda);
const double weightedTotal = std::accumulate(std::begin(m_chemicalFormula),
std::end(m_chemicalFormula),
0.,
[](double subtotal, const FormulaUnit &right) {
return subtotal + right.atom->neutron.tot_scatt_xs * right.multiplicity;
}) / m_atomTotal;
if (! std::isnormal(weightedTotal)) {
return 0.;
} else {
return weightedTotal;
}
}
/**
* Get the absorption cross section for a given wavelength.
* CURRENTLY This assumes a linear dependence on the wavelength with the
* reference
* wavelength = NeutronAtom::ReferenceLambda angstroms.
* @param lambda :: The wavelength to evaluate the cross section
* @returns The value of the absoprtioncross section at
* the given wavelength
*/
double Material::absorbXSection(const double lambda) const {
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
const double weightedTotal = std::accumulate(std::begin(m_chemicalFormula),
std::end(m_chemicalFormula),
0.,
[](double subtotal, const FormulaUnit &right) {
return subtotal + right.atom->neutron.abs_scatt_xs * right.multiplicity;
}) / m_atomTotal;
if (! std::isnormal(weightedTotal)) {
return 0.;
} else {
return weightedTotal * (lambda / NeutronAtom::ReferenceLambda);
}
}
double Material::cohScatterLength(const double lambda) const {
UNUSED_ARG(lambda);
const double weightedTotal = std::accumulate(std::begin(m_chemicalFormula),
std::end(m_chemicalFormula),
0.,
[](double subtotal, const FormulaUnit &right) {
return subtotal + right.atom->neutron.coh_scatt_length * right.multiplicity;
}) / m_atomTotal;
if (! std::isnormal(weightedTotal)) {
return 0.;
} else {
return weightedTotal;
}
}
double Material::incohScatterLength(const double lambda) const {
UNUSED_ARG(lambda);
const double weightedTotal = std::accumulate(std::begin(m_chemicalFormula),
std::end(m_chemicalFormula),
0.,
[](double subtotal, const FormulaUnit &right) {
return subtotal + right.atom->neutron.inc_scatt_length * right.multiplicity;
}) / m_atomTotal;
if (! std::isnormal(weightedTotal)) {
return 0.;
} else {
return weightedTotal;
}
}
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
double
Material::cohScatterLengthReal(const double lambda) const {
UNUSED_ARG(lambda);
const double weightedTotal = std::accumulate(std::begin(m_chemicalFormula),
std::end(m_chemicalFormula),
0.,
[](double subtotal, const FormulaUnit &right) {
return subtotal + right.atom->neutron.coh_scatt_length_real * right.multiplicity;
}) / m_atomTotal;
if (! std::isnormal(weightedTotal)) {
return 0.;
} else {
return weightedTotal;
}
}
double Material::incohScatterLengthReal(
const double lambda) const {
UNUSED_ARG(lambda);
const double weightedTotal = std::accumulate(std::begin(m_chemicalFormula),
std::end(m_chemicalFormula),
0.,
[](double subtotal, const FormulaUnit &right) {
return subtotal + right.atom->neutron.inc_scatt_length_real * right.multiplicity;
}) / m_atomTotal;
if (! std::isnormal(weightedTotal)) {
return 0.;
} else {
return weightedTotal;
}
}
double Material::totalScatterLength(const double lambda) const {
UNUSED_ARG(lambda);
const double weightedTotal = std::accumulate(std::begin(m_chemicalFormula),
std::end(m_chemicalFormula),
0.,
[](double subtotal, const FormulaUnit &right) {
return subtotal + right.atom->neutron.tot_scatt_length * right.multiplicity;
}) / m_atomTotal;
if (! std::isnormal(weightedTotal)) {
return 0.;
} else {
return weightedTotal;
}
}
double Material::cohScatterLengthRealSqrd(const double lambda) const{
UNUSED_ARG(lambda);
const double weightedTotal = std::accumulate(std::begin(m_chemicalFormula),
std::end(m_chemicalFormula),
0.,
[](double subtotal, const FormulaUnit &right) {
const double value = right.atom->neutron.coh_scatt_length_real;
return subtotal + value * value * right.multiplicity;
}) / m_atomTotal;
if (! std::isnormal(weightedTotal)) {
return 0.;
} else {
return weightedTotal;
}
}
double Material::incohScatterLengthRealSqrd(
const double lambda) const{
UNUSED_ARG(lambda);
const double weightedTotal = std::accumulate(std::begin(m_chemicalFormula),
std::end(m_chemicalFormula),
0.,
[](double subtotal, const FormulaUnit &right) {
const double value = right.atom->neutron.inc_scatt_length_real;
return subtotal + value * value * right.multiplicity;
}) / m_atomTotal;
if (! std::isnormal(weightedTotal)) {
return 0.;
} else {
return weightedTotal;
}
}
double
Material::cohScatterLengthSqrd(const double lambda) const{
UNUSED_ARG(lambda);
const double weightedTotal = std::accumulate(std::begin(m_chemicalFormula),
std::end(m_chemicalFormula),
0.,
[](double subtotal, const FormulaUnit &right) {
const double value = right.atom->neutron.coh_scatt_length;
return subtotal + value * value * right.multiplicity;
}) / m_atomTotal;
if (! std::isnormal(weightedTotal)) {
return 0.;
} else {
return weightedTotal;
}
}
double Material::incohScatterLengthSqrd(
const double lambda) const{
UNUSED_ARG(lambda);
const double weightedTotal = std::accumulate(std::begin(m_chemicalFormula),
std::end(m_chemicalFormula),
0.,
[](double subtotal, const FormulaUnit &right) {
const double value = right.atom->neutron.inc_scatt_length;
return subtotal + value * value * right.multiplicity;
}) / m_atomTotal;
if (! std::isnormal(weightedTotal)) {
return 0.;
} else {
return weightedTotal;
}
}
double Material::totalScatterLengthSqrd(
const double lambda) const{
UNUSED_ARG(lambda);
const double weightedTotal = std::accumulate(std::begin(m_chemicalFormula),
std::end(m_chemicalFormula),
0.,
[](double subtotal, const FormulaUnit &right) {
const double value = right.atom->neutron.tot_scatt_length;
return subtotal + value * value * right.multiplicity;
}) / m_atomTotal;
if (! std::isnormal(weightedTotal)) {
return 0.;
} else {
return weightedTotal;
}
}
/** Save the object to an open NeXus file.
* @param file :: open NeXus file
* @param group :: name of the group to create
*/
void Material::saveNexus(::NeXus::File *file, const std::string &group) const {
file->makeGroup(group, "NXdata", 1);
file->putAttr("version", 2);
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
// determine how the information will be stored
std::string style = "formula"; // default is a chemical formula
if (m_chemicalFormula.size() == 0) {
style = "empty";
} else if (m_chemicalFormula.size() == 1) {
if (m_chemicalFormula[0].atom->symbol == "user") {
style = "userdefined";
}
}
file->putAttr("formulaStyle", style);
// write the actual information out
if (style == "formula") {
std::stringstream formula;
for ( const auto &formulaUnit : m_chemicalFormula) {
if (formulaUnit.atom->a_number != 0) {
formula << "(";
}
formula << formulaUnit.atom->symbol;
if (formulaUnit.atom->a_number != 0) {
formula << formulaUnit.atom->a_number << ")";
}
formula << formulaUnit.multiplicity << " ";
}
file->writeData("chemical_formula", formula.str());
} else if (style == "userdefined") {
file->writeData("coh_scatt_length_real", m_chemicalFormula[0].atom->neutron.coh_scatt_length_real);
file->writeData("coh_scatt_length_img", m_chemicalFormula[0].atom->neutron.coh_scatt_length_img);
file->writeData("inc_scatt_length_real", m_chemicalFormula[0].atom->neutron.inc_scatt_length_real);
file->writeData("inc_scatt_length_img", m_chemicalFormula[0].atom->neutron.inc_scatt_length_img);
file->writeData("coh_scatt_xs", m_chemicalFormula[0].atom->neutron.coh_scatt_xs);
file->writeData("inc_scatt_xs", m_chemicalFormula[0].atom->neutron.inc_scatt_xs);
file->writeData("tot_scatt_xs", m_chemicalFormula[0].atom->neutron.tot_scatt_xs);
file->writeData("abs_scatt_xs", m_chemicalFormula[0].atom->neutron.abs_scatt_xs);
file->writeData("tot_scatt_length", m_chemicalFormula[0].atom->neutron.tot_scatt_length);
file->writeData("coh_scatt_length", m_chemicalFormula[0].atom->neutron.coh_scatt_length);
file->writeData("inc_scatt_length", m_chemicalFormula[0].atom->neutron.inc_scatt_length);
}
file->writeData("number_density", m_numberDensity);
file->writeData("temperature", m_temperature);
file->writeData("pressure", m_pressure);
file->closeGroup();
}
/** Load the object from an open NeXus file.
* @param file :: open NeXus file
* @param group :: name of the group to open
*/
void Material::loadNexus(::NeXus::File *file, const std::string &group) {
file->openGroup(group, "NXdata");
file->getAttr("name", m_name);
int version;
file->getAttr("version", version);
if (version == 1) {
// Find the element
uint16_t element_Z, element_A;
file->readData("element_Z", element_Z);
file->readData("element_A", element_A);
try {
m_chemicalFormula.clear();
if (element_Z > 0) {
FormulaUnit formulaUnit {boost::make_shared<Atom>(Mantid::PhysicalConstants::getAtom(element_Z, element_A)), 1.};
m_chemicalFormula.push_back(formulaUnit);
} else {
FormulaUnit formulaUnit {boost::make_shared<Atom>(Mantid::PhysicalConstants::getNeutronAtom(element_Z, element_A)), 1.};
m_chemicalFormula.push_back(formulaUnit);
}
} catch (std::runtime_error &) { /* ignore and use the default */
}
} else if (version == 2) {
std::string style;
file->getAttr("formulaStyle", style);
if (style == "formula") {
std::string formula;
file->readData("chemical_formula", formula);
this->m_chemicalFormula = Material::parseChemicalFormula(formula);
this->countAtoms();
} else if (style == "userdefined") {
NeutronAtom neutron;
file->readData("coh_scatt_length_real", neutron.coh_scatt_length_real);
file->readData("coh_scatt_length_img", neutron.coh_scatt_length_img);
file->readData("inc_scatt_length_real", neutron.inc_scatt_length_real);
file->readData("inc_scatt_length_img", neutron.inc_scatt_length_img);
file->readData("coh_scatt_xs", neutron.coh_scatt_xs);
file->readData("inc_scatt_xs", neutron.inc_scatt_xs);
file->readData("tot_scatt_xs", neutron.tot_scatt_xs);
file->readData("abs_scatt_xs", neutron.abs_scatt_xs);
file->readData("tot_scatt_length", neutron.tot_scatt_length);
file->readData("coh_scatt_length", neutron.coh_scatt_length);
file->readData("inc_scatt_length", neutron.inc_scatt_length);
FormulaUnit formulaUnit {boost::make_shared<Atom>(neutron), 1.};
m_chemicalFormula.push_back(formulaUnit);
}
// the other option is empty which does not need to be addressed
}else {
throw std::runtime_error("Only know how to read version 1 or 2 for Material");
file->readData("number_density", m_numberDensity);
file->readData("temperature", m_temperature);
file->readData("pressure", m_pressure);
file->closeGroup();
}
namespace { // anonymous namespace to hide the function
str_pair
getAtomName(const std::string &text) // TODO change to get number after letters
{
// one character doesn't need
if (text.size() <= 1)
return std::make_pair(text, "");
// check the second character
const char *s;
s = text.c_str();
if ((s[1] >= '0' && s[1] <= '9') || s[1] == '.')
return std::make_pair(text.substr(0, 1), text.substr(1));
else
return std::make_pair(text.substr(0, 2), text.substr(2));
}
}
Material::ChemicalFormula
Material::parseChemicalFormula(const std::string chemicalSymbol) {
Material::ChemicalFormula CF;
tokenizer tokens(chemicalSymbol, " -",
Mantid::Kernel::StringTokenizer::TOK_IGNORE_EMPTY);
for (const auto &atom : tokens) {
std::string name;
float numberAtoms = 1;
uint16_t aNumber = 0;
// split out the isotope bit
if (atom.find('(') != std::string::npos) {
size_t end = atom.find(')');
if (end == std::string::npos) {
std::stringstream msg;
msg << "Failed to parse isotope \"" << atom << "\"";
throw std::runtime_error(msg.str());
}
// get the number of atoms
std::string numberAtomsStr = atom.substr(end + 1);
if (!numberAtomsStr.empty())
numberAtoms = boost::lexical_cast<float>(numberAtomsStr);
// split up the atom and isotope number
name = atom.substr(1, end - 1);
str_pair temp = getAtomName(name);
name = temp.first;
aNumber = boost::lexical_cast<uint16_t>(temp.second);
} else // for non-isotopes
str_pair temp = getAtomName(atom);
name = temp.first;
if (!temp.second.empty())
numberAtoms = boost::lexical_cast<float>(temp.second);
Material::FormulaUnit formulaUnit {boost::make_shared<PhysicalConstants::Atom>(getAtom(name, aNumber)),
static_cast<double>(numberAtoms)};
CF.push_back(formulaUnit);
} catch (boost::bad_lexical_cast &e) {
std::stringstream msg;
msg << "While trying to parse atom \"" << atom
<< "\" encountered bad_lexical_cast: " << e.what();
throw std::runtime_error(msg.str());