Skip to content
Snippets Groups Projects
Commit c4c3925c authored by Yurov, Dmitry's avatar Yurov, Dmitry
Browse files

Changes to python-user interface for MaterialBySLD

Redmine: #1946
Also related changes to unit tests
parent 87588b9a
No related branches found
No related tags found
No related merge requests found
......@@ -17,23 +17,22 @@
namespace
{
// Returns SLD-like complex value, where real part is SLD and imaginary one is
// wavelength-independent absorptive term
inline complex_t getSLD(double sld, double abs_term)
inline double getWlPrefactor(double wavelength)
{
return complex_t(sld, -abs_term / 2.0);
return wavelength * wavelength / M_PI;
}
inline double getWlPrefactor(double wavelength)
inline complex_t getSLD(double sld_real, double sld_imag)
{
return wavelength * wavelength / M_PI;
return complex_t(sld_real, -sld_imag);
}
}
MaterialBySLDImpl::MaterialBySLDImpl(const std::string& name, double sld,
double abs_term,
kvector_t magnetization)
: MagneticMaterialImpl(name, magnetization), m_sld(sld), m_abs_term(abs_term)
MaterialBySLDImpl::MaterialBySLDImpl(const std::string& name, double sld_real, double sld_imag,
kvector_t magnetization)
: MagneticMaterialImpl(name, magnetization)
, m_sld_real(sld_real)
, m_sld_imag(sld_imag)
{}
MaterialBySLDImpl* MaterialBySLDImpl::clone() const
......@@ -48,23 +47,23 @@ complex_t MaterialBySLDImpl::refractiveIndex(double wavelength) const
complex_t MaterialBySLDImpl::refractiveIndex2(double wavelength) const
{
return 1.0 - getWlPrefactor(wavelength) * getSLD(m_sld, m_abs_term);
return 1.0 - getWlPrefactor(wavelength) * getSLD(m_sld_real, m_sld_imag);
}
complex_t MaterialBySLDImpl::materialData() const
{
return complex_t(m_sld, m_abs_term);
return complex_t(m_sld_real, m_sld_imag);
}
complex_t MaterialBySLDImpl::scalarSubtrSLD(const WavevectorInfo& wavevectors) const
{
double wavelength = wavevectors.getWavelength();
return 1.0 / getWlPrefactor(wavelength) - getSLD(m_sld, m_abs_term);
return 1.0 / getWlPrefactor(wavelength) - getSLD(m_sld_real, m_sld_imag);
}
void MaterialBySLDImpl::print(std::ostream& ostr) const
{
ostr << "MaterialBySLD:" << getName() << "<" << this << ">{ "
<< "sld=" << m_sld << ", absorp_term=" << m_abs_term
<< "sld_real=" << m_sld_real << ", sld_imag = " << m_sld_imag
<< ", B=" << magnetization() << "}";
}
......@@ -24,8 +24,9 @@
class BA_CORE_API_ MaterialBySLDImpl : public MagneticMaterialImpl
{
public:
friend BA_CORE_API_ Material MaterialBySLD(const std::string& name, double sld, double abs_term,
kvector_t magnetization);
friend BA_CORE_API_ Material createMaterialBySLDInNativeUnits(const std::string& name,
double sld_real, double sld_imag,
kvector_t magnetization);
virtual ~MaterialBySLDImpl() = default;
......@@ -54,13 +55,11 @@ public:
void print(std::ostream &ostr) const override;
private:
MaterialBySLDImpl(const std::string& name, double sld, double abs_term,
kvector_t magnetization);
MaterialBySLDImpl(const std::string& name, double sld_real, double sld_imag,
kvector_t magnetization);
double m_sld; //!< product of number density and coherent scattering length
//!< (scattering length density)
double m_abs_term; //!< product of number density and
//!< absorption cross-section normalized to wavelength
double m_sld_real; //!< complex-valued scattering length density
double m_sld_imag; //!< imaginary part of scattering length density (negative by default)
};
#endif /* MATERIALBYSLDIMPL_H_ */
......@@ -47,17 +47,24 @@ Material HomogeneousMaterial()
return HomogeneousMaterial("vacuum", 0.0, 0.0, kvector_t{});
}
Material MaterialBySLD(const std::string& name, double sld, double abs_term,
Material MaterialBySLD()
{
return MaterialBySLD("vacuum", 0.0, 0.0, kvector_t{});
}
Material MaterialBySLD(const std::string& name, double sld_real, double sld_imag,
kvector_t magnetization)
{
std::unique_ptr<MaterialBySLDImpl> mat_impl(
new MaterialBySLDImpl(name, sld, abs_term, magnetization));
return Material(std::move(mat_impl));
// A^{-2} = 100 nm^{-2}
return createMaterialBySLDInNativeUnits(name, sld_real * 100.0, sld_imag * 100.0, magnetization);
}
Material MaterialBySLD()
Material createMaterialBySLDInNativeUnits(const std::string& name, double sld_real, double sld_imag,
kvector_t magnetization)
{
return MaterialBySLD("vacuum", 0.0, 0.0, kvector_t{});
std::unique_ptr<MaterialBySLDImpl> mat_impl(
new MaterialBySLDImpl(name, sld_real, sld_imag, magnetization));
return Material(std::move(mat_impl));
}
Material createAveragedMaterial(const Material& layer_mat,
......@@ -94,7 +101,8 @@ Material createAveragedMaterial(const Material& layer_mat,
complex_t (*avrData)(const Material&)
= [](const Material& mat) { return mat.materialData(); };
const complex_t avr_mat_data = averageData<complex_t>(layer_mat, regions, avrData);
return MaterialBySLD(avr_mat_name, avr_mat_data.real(), avr_mat_data.imag(), mag_avr);
return createMaterialBySLDInNativeUnits(avr_mat_name, avr_mat_data.real(), avr_mat_data.imag(),
mag_avr);
} else
throw std::runtime_error("Error in CalculateAverageMaterial: unknown material type.");
}
......
......@@ -30,9 +30,12 @@ BA_CORE_API_ Material HomogeneousMaterial(const std::string& name, double delta,
//! @ingroup materials
//! Constructs a material with _name_, _refractive_index_ and _magnetization_ (in A/m). Alternatively,
//! \f$\delta\f$ and \f$\beta\f$ for refractive index \f$n = 1 - \delta + i \beta\f$ can be passed directly.
//! With no parameters given, constructs default (vacuum) material with \f$n = 1\f$ and zero magnetization.
//! Constructs a material with _name_, _refractive_index_ and _magnetization_ (in A/m).
//! Alternatively,
//! \f$\delta\f$ and \f$\beta\f$ for refractive index \f$n = 1 - \delta + i \beta\f$ can be passed
//! directly.
//! With no parameters given, constructs default (vacuum) material with \f$n = 1\f$ and zero
//! magnetization.
BA_CORE_API_ Material HomogeneousMaterial(const std::string& name, complex_t refractive_index,
kvector_t magnetization = kvector_t());
......@@ -42,25 +45,39 @@ BA_CORE_API_ Material MaterialBySLD();
//! @ingroup materials
//! Constructs a wavelength-independent material with given sld and absorptive term
//! Absorptive term is wavelength-independent (normalized to a wavelength)
//! and can be considered as inverse of imaginary part of complex scattering length density:
//! \f$ SLD = (N b_{coh}, N \sigma_{abs}(\lambda_0) / \lambda_0) \f$
//! Here \f$ N \f$ - material number density,
//! \f$ b_{coh} \f$ - bound coherent scattering length
//! \f$ \sigma_{abs}(\lambda_0) \f$ - absorption cross-section at \f$ \lambda_0 \f$ wavelength.
//! With no parameters given, constructs default (vacuum) material with zero sld and zero magnetization.
//! Constructs a wavelength-independent material with a given complex-valued
//! scattering lenght density (SLD).
//! SLD values for a wide variety of materials can be found on
//! https://sld-calculator.appspot.com/
//! and
//! https://www.ncnr.nist.gov/resources/activation/
//! By convention, SLD imaginary part is treated as negative by default, which corresponds to
//! attenuation of the signal.
//! With no parameters given, MaterialBySLD constructs default (vacuum) material with zero sld
//! and zero magnetization.
//! @param name: material name
//! @param sld: scattering length density, \f$ nm^{-2} \f$
//! @param abs_term: wavelength-independent absorptive term, \f$ nm^{-2} \f$
//! @param sld_real: real part of the scattering length density, inverse square angstroms
//! @param sld_imag: imaginary part of the scattering length density, inverse square angstroms
//! @param magnetization: magnetization (in A/m)
BA_CORE_API_ Material MaterialBySLD(const std::string& name, double sld, double abs_term,
BA_CORE_API_ Material MaterialBySLD(const std::string& name, double sld_real, double sld_imag,
kvector_t magnetization = kvector_t());
#ifndef SWIG
//! @ingroup materials
//! Constructs a wavelength-independent material with a given complex-valued
//! scattering lenght density (SLD). SLD units are \f$ nm^{-2} \f$.
BA_CORE_API_ Material createMaterialBySLDInNativeUnits(const std::string& name, double sld_real,
double sld_imag, kvector_t magnetization);
//! @ingroup materials
//! Creates averaged material. Square refractive index of returned material is arithmetic mean over
//! _regions_ and _layer_mat_. Magnetization (if present) is averaged linearly.
BA_CORE_API_ Material createAveragedMaterial(const Material& layer_mat, const std::vector<HomogeneousRegion>& regions);
BA_CORE_API_ Material createAveragedMaterial(const Material& layer_mat,
const std::vector<HomogeneousRegion>& regions);
#endif //SWIG
#endif /* MATERIALFACTORYFUNCS_H_ */
......@@ -33,8 +33,8 @@ TEST_F(MaterialTest, MaterialConstruction)
EXPECT_EQ(material_data, material2.materialData());
EXPECT_EQ(magnetism, material2.magnetization());
Material material3
= MaterialBySLD("MagMaterial", material_data.real(), material_data.imag(), magnetism);
Material material3 = createMaterialBySLDInNativeUnits("MagMaterial", material_data.real(),
material_data.imag(), magnetism);
EXPECT_EQ("MagMaterial", material3.getName());
EXPECT_EQ(material_data, material3.materialData());
EXPECT_EQ(magnetism, material3.magnetization());
......@@ -53,8 +53,6 @@ TEST_F(MaterialTest, MaterialConstruction)
EXPECT_EQ(default_magnetism, material5.magnetization());
Material material6 = MaterialBySLD("Material", material_data.real(), material_data.imag());
EXPECT_EQ("Material", material6.getName());
EXPECT_EQ(material_data, material6.materialData());
EXPECT_EQ(default_magnetism, material6.magnetization());
}
......@@ -107,7 +105,7 @@ TEST_F(MaterialTest, DefaultMaterials)
TEST_F(MaterialTest, ComputationTest)
{
// Reference data for Fe taken from
// https://sld-calculator.appspot.com/save
// https://sld-calculator.appspot.com
// http://www.ati.ac.at/~neutropt/scattering/table.html
const double bc = 9.45e-6; // nm, bound coherent scattering length
const double abs_cs = 2.56e-10; // nm^2, absorption cross-section for 2200 m/s neutrons
......@@ -116,20 +114,31 @@ TEST_F(MaterialTest, ComputationTest)
const double avog_number = 6.022e+23; // mol^{-1}, Avogadro number
const double density = 7.874e-21; // g/nm^3, Fe material density
const double number_density = avog_number * density / mol_mass; // 1/nm^3, Fe number density
const double sld = number_density * bc;
const double abs_term = number_density * abs_cs / basic_wavelength;
const double sld_real = number_density * bc;
const double sld_imag = number_density * abs_cs / ( 2.0 * basic_wavelength);
const complex_t sld_ref(8.0241e-04, // nm^{-2}, reference data
6.0448e-8); // taken from https://sld-calculator.appspot.com/
// checking input data and reference values are the same
EXPECT_NEAR(2.0 * (sld_ref.real() - sld_real) / (sld_ref.real() + sld_real), 0.0, 1e-3);
EXPECT_NEAR(2.0 * (sld_ref.imag() - sld_imag) / (sld_ref.imag() + sld_imag), 0.0, 1e-3);
const complex_t sld_ref(8.0241e-04,
-6.0448e-8); // nm^{-2}, reference data, wavelength-independent
const double wl_factor_2200 = basic_wavelength * basic_wavelength / M_PI;
const double wl_factor_1100 = 4.0 * basic_wavelength * basic_wavelength / M_PI;
Material material = MaterialBySLD("Fe", sld, abs_term);
const complex_t sld_res_2200
// MaterialBySLD accepts sld in AA^{-2}
Material material = MaterialBySLD("Fe", sld_real / 100, sld_imag / 100);
complex_t sld_res_2200
= (1.0 - material.refractiveIndex2(basic_wavelength)) / wl_factor_2200;
const complex_t sld_res_1100
complex_t sld_res_1100
= (1.0 - material.refractiveIndex2(2.0 * basic_wavelength)) / wl_factor_1100;
// change the sign of imaginary part according to internal convention
sld_res_2200 = std::conj(sld_res_2200);
sld_res_1100 = std::conj(sld_res_1100);
EXPECT_NEAR(0.0, (sld_ref.real() - sld_res_2200.real()) / sld_ref.real(), 1e-3);
EXPECT_NEAR(0.0, (sld_ref.imag() - sld_res_2200.imag()) / sld_ref.imag(), 1e-3);
EXPECT_NEAR(0.0, (sld_ref.real() - sld_res_1100.real()) / sld_ref.real(), 1e-3);
......@@ -168,7 +177,7 @@ TEST_F(MaterialTest, AveragedMaterialTest)
EXPECT_EQ(material_avr2.magnetization(), kvector_t(0.5, 0.0, 0.0));
EXPECT_TRUE(material_avr2.typeID() == MATERIAL_TYPES::RefractiveMaterial);
const Material material3 = MaterialBySLD("Material3", 0.5, 0.5, magnetization);
const Material material3 = createMaterialBySLDInNativeUnits("Material3", 0.5, 0.5, magnetization);
EXPECT_THROW(createAveragedMaterial(material3, regions), std::runtime_error);
const Material material4 = HomogeneousMaterial();
......@@ -234,25 +243,10 @@ TEST_F(MaterialTest, MaterialMove)
kvector_t magnetism = kvector_t(3.0, 4.0, 5.0);
Material material = HomogeneousMaterial("MagMaterial", refIndex, magnetism);
Material move(std::move(material));
EXPECT_EQ("MagMaterial", move.getName());
EXPECT_EQ(material_data, move.materialData());
EXPECT_EQ(magnetism, move.magnetization());
EXPECT_TRUE(material.isEmpty());
}
TEST_F(MaterialTest, MaterialAssignment)
{
Material material = MaterialBySLD("Material", 1.0, 1.0);
Material material_ref = MaterialBySLD("Material", 1.0, 1.0);
material = material;
EXPECT_TRUE(material == material_ref);
material = std::move(material);
EXPECT_TRUE(material == material_ref);
Material material2 = std::move(material);
Material moved_material(std::move(material));
EXPECT_EQ("MagMaterial", moved_material.getName());
EXPECT_EQ(material_data, moved_material.materialData());
EXPECT_EQ(magnetism, moved_material.magnetization());
EXPECT_TRUE(material.isEmpty());
EXPECT_THROW(material2 = material, Exceptions::NullPointerException);
EXPECT_THROW(Material material2 = material, Exceptions::NullPointerException);
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment