diff --git a/Core/Material/MaterialBySLDImpl.cpp b/Core/Material/MaterialBySLDImpl.cpp index d5b72f48d326ad92ad79a22f0c45dfc95eb022f8..03e2209a3c52ff9ca99d29f8993813ac8b8781b7 100644 --- a/Core/Material/MaterialBySLDImpl.cpp +++ b/Core/Material/MaterialBySLDImpl.cpp @@ -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() << "}"; } diff --git a/Core/Material/MaterialBySLDImpl.h b/Core/Material/MaterialBySLDImpl.h index 5e43f4d3b36053bc97be9fa47df5a7b536ac5683..07c81512390257dda458ff3018ccb8d1d66dbf1c 100644 --- a/Core/Material/MaterialBySLDImpl.h +++ b/Core/Material/MaterialBySLDImpl.h @@ -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_ */ diff --git a/Core/Material/MaterialFactoryFuncs.cpp b/Core/Material/MaterialFactoryFuncs.cpp index 2597127cfa08b9b7d2749fcd9fb9fe5885ea0c1b..02f90ab174a897a2ac1b14656c5b4d92720f0a63 100644 --- a/Core/Material/MaterialFactoryFuncs.cpp +++ b/Core/Material/MaterialFactoryFuncs.cpp @@ -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."); } diff --git a/Core/Material/MaterialFactoryFuncs.h b/Core/Material/MaterialFactoryFuncs.h index 4406b996decc7a782f1315a5d224c99d6bedec50..56e6a3b090f419ea485bfee28c972b55e6be3a51 100644 --- a/Core/Material/MaterialFactoryFuncs.h +++ b/Core/Material/MaterialFactoryFuncs.h @@ -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_ */ diff --git a/Tests/UnitTests/Core/Other/MaterialTest.h b/Tests/UnitTests/Core/Other/MaterialTest.h index b1f29a57d15280367503553d724719ba3d68d13d..19a115148ccf519af07ee5aabc866ac06597bfb3 100644 --- a/Tests/UnitTests/Core/Other/MaterialTest.h +++ b/Tests/UnitTests/Core/Other/MaterialTest.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); }