From 6c6fe7ee52593f31932a9a155a98222899653b9d Mon Sep 17 00:00:00 2001
From: lbocklag <lars.bocklage@desy.de>
Date: Fri, 13 Dec 2024 12:54:53 +0100
Subject: [PATCH 1/2] added natural abundance to MoessbauerIsotope for precise
average mass calculation in material
---
docs/source/news/news.rst | 2 +
src/nexus/clib/classes/material.cpp | 14 ++++-
src/nexus/clib/classes/moessbauer.cpp | 21 ++++++-
src/nexus/clib/classes/moessbauer.h | 12 +++-
src/nexus/clib/classes/moessbauer.i | 7 ++-
src/nexus/clib/classes/sample.cpp | 4 +-
src/nexus/lib/moessbauer.py | 88 ++++++++++++++++++++-------
7 files changed, 115 insertions(+), 33 deletions(-)
diff --git a/docs/source/news/news.rst b/docs/source/news/news.rst
index a7ddd4dd..e2d92974 100644
--- a/docs/source/news/news.rst
+++ b/docs/source/news/news.rst
@@ -68,6 +68,8 @@ API changes:
It enables to use other than least-square prolems, like log-likelihood estimators.
It does not work with the ceres solver modules, which always solve least-square problems.
+#. :attr:`natural_abundance` attribute added to :class:`MoessbauerIsotope`.
+
changes:
diff --git a/src/nexus/clib/classes/material.cpp b/src/nexus/clib/classes/material.cpp
index 531b0f68..7d028c20 100644
--- a/src/nexus/clib/classes/material.cpp
+++ b/src/nexus/clib/classes/material.cpp
@@ -85,17 +85,25 @@ void Material::UpdateComposition()
for (const CompositionEntry& element : composition)
{
- double mass = element_maps::atomic_mass_map.find(element.first)->second;
+ double mass = element_maps::atomic_mass_map.find(element.first)->second; // averaged atomic mass over natural abundance
+ // correct mass for given abundance of resonant isotope
+ // X * M_resonant_isotope + Y M_other_isotopes = A (averaged atomic mass)
+ // X + Y = 1
+ // M_other_isotopes = (A - X * M_resonant_isotope) / (1 - X)
if (element.first == isotope->element)
- mass = abundance->value * isotope->mass + (1.0 - abundance->value) * mass;
+ {
+ const double mass_non_res_isotopes = (mass - isotope->natural_abundance * isotope->mass) / (1.0 - isotope->natural_abundance);
+
+ mass = abundance->value * isotope->mass + (1.0 - abundance->value) * mass_non_res_isotopes;
+ }
average_mole_mass += element.second / sum * mass;
}
// density in g/cm^3, convert to g/m^3 (1e6)
// average_mole_mass is in g/mole
- total_number_density = constants::kNa * (density->value * 1.0e6) / (average_mole_mass); // 1/m^3
+ total_number_density = constants::kNa * (density->value * 1.0e6) / average_mole_mass; // 1/m^3
for (const CompositionEntry& element : composition)
{
diff --git a/src/nexus/clib/classes/moessbauer.cpp b/src/nexus/clib/classes/moessbauer.cpp
index 7314bec8..ca396133 100644
--- a/src/nexus/clib/classes/moessbauer.cpp
+++ b/src/nexus/clib/classes/moessbauer.cpp
@@ -44,7 +44,8 @@ MoessbauerIsotope::MoessbauerIsotope(
Var* const gfactor_excited_,
Var* const quadrupole_ground_,
Var* const quadrupole_excited_,
- Var* const interference_term_
+ Var* const interference_term_,
+ const double natural_abundance_
) :
isotope(isotope_),
element(element_),
@@ -60,9 +61,23 @@ MoessbauerIsotope::MoessbauerIsotope(
gfactor_excited(gfactor_excited_),
quadrupole_ground(quadrupole_ground_),
quadrupole_excited(quadrupole_excited_),
- interference_term(interference_term_)
+ interference_term(interference_term_),
+ natural_abundance(natural_abundance_)
{
Update();
+
+ if (natural_abundance < 0.0)
+ {
+ errors::ErrorMessage("MoessbauerIsotope", "natural abundance smaller 0. Set to 0.");
+
+ natural_abundance = 0.0;
+ }
+ else if (natural_abundance > 1.0)
+ {
+ errors::ErrorMessage("MoessbauerIsotope", "natural abundance larger 1. Set to 1.");
+
+ natural_abundance = 1.0;
+ }
}
MoessbauerIsotope MoessbauerIsotope::Copy() const
@@ -172,4 +187,4 @@ void MoessbauerIsotope::PopulateFitHandler(OptimizerHandler* fit_handler)
fit_handler->PrepareFitVariable(quadrupole_excited);
fit_handler->PrepareFitVariable(interference_term);
-}
\ No newline at end of file
+}
diff --git a/src/nexus/clib/classes/moessbauer.h b/src/nexus/clib/classes/moessbauer.h
index ce52513c..235d2d7e 100644
--- a/src/nexus/clib/classes/moessbauer.h
+++ b/src/nexus/clib/classes/moessbauer.h
@@ -115,6 +115,10 @@ Args:
.. versionchanged:: 1.2.0
:class:`Var` input possible now.
+ natural_abundance (float): Natural abundance of the isotope.
+
+ .. versionadded:: 2.0.0
+
Attributes:
isotope (string): User identifier. Should be given in the format "mass number-element", e.g. "57-Fe".
@@ -214,6 +218,9 @@ Attributes:
It is given by :math:`\mu_g = g_g \mu_N I_g`, where :math:`I_g` is given in units of :math:`\hbar`.
magnetic_moment_excited (float): Magnetic moment of the excited state (eV/T).
It is given by :math:`\mu_e = g_e \mu_N I_e`, where :math:`I_e` is given in units of :math:`\hbar`.
+ natural_abundance (float): Natural abundance of the isotope.
+
+ .. versionadded:: 2.0.0
*/
struct MoessbauerIsotope {
public:
@@ -232,7 +239,8 @@ public:
Var* const gfactor_excited,
Var* const quadrupole_ground,
Var* const quadrupole_excited,
- Var* const interference_term
+ Var* const interference_term,
+ const double natural_abundance
);
/**
@@ -278,6 +286,8 @@ public:
mutable Var* interference_term = nullptr;
+ double natural_abundance = 0.0;
+
// derived parameters
int atomic_number = 0;
diff --git a/src/nexus/clib/classes/moessbauer.i b/src/nexus/clib/classes/moessbauer.i
index 54a134fd..2a5d5ff9 100644
--- a/src/nexus/clib/classes/moessbauer.i
+++ b/src/nexus/clib/classes/moessbauer.i
@@ -43,7 +43,8 @@
gfactor_excited = 0.0,
quadrupole_ground = 0.0,
quadrupole_excited = 0.0,
- interference_term = 0.0):
+ interference_term = 0.0,
+ natural_abundance = 0.0):
if isinstance(internal_conversion, (int, float)):
internal_conversion = Var(internal_conversion, 0.0, 1.0e9, False, "")
@@ -88,7 +89,8 @@
gfactor_excited,
quadrupole_ground,
quadrupole_excited,
- interference_term)
+ interference_term,
+ natural_abundance)
)
%}
@@ -180,6 +182,7 @@
output += " .magnetic_moment_ground (eV/T) = {}\n".format(self.magnetic_moment_ground)
output += " .magnetic_moment_excited (eV/T) = {}\n".format(self.magnetic_moment_excited)
output += " .nuclear_cross_section (converted to kbarn) = {}\n".format(self.nuclear_cross_section*1e25)
+ output += " .natural_abundance = {}\n".format(self.natural_abundance)
return output
}
}
diff --git a/src/nexus/clib/classes/sample.cpp b/src/nexus/clib/classes/sample.cpp
index ebe1d3ee..3b2d9699 100644
--- a/src/nexus/clib/classes/sample.cpp
+++ b/src/nexus/clib/classes/sample.cpp
@@ -59,9 +59,9 @@ namespace sample_definitions
static Composition kAirComposition = { std::make_pair("N", 1.562), std::make_pair("O", 0.42), std::make_pair("C", 0.0003), std::make_pair("Ar", 0.0094) };
static Var* const kDensity = new Var(0.0012041, 0.0, 0.0012041, false);
static Var* const kAbundance = new Var(0.0, 0.0, 0.0, false);
- static Var* const kMoessbuqerDummy = new Var(0.0, 0.0, 0.0, false);
+ static Var* const kMoessbauerDummy = new Var(0.0, 0.0, 0.0, false);
static Var* const kLambMoessbauer = new Var(0.0, 0.0, 0.0, false);
- static MoessbauerIsotope* const kNone = new MoessbauerIsotope("none", "", 0.0, 0.0, 0.0, kMoessbuqerDummy, Multipolarity::M1, kMoessbuqerDummy, 0.0, 0.0, kMoessbuqerDummy, kMoessbuqerDummy, kMoessbuqerDummy, kMoessbuqerDummy, kMoessbuqerDummy);
+ static MoessbauerIsotope* const kNone = new MoessbauerIsotope("none", "", 0.0, 0.0, 0.0, kMoessbauerDummy, Multipolarity::M1, kMoessbauerDummy, 0.0, 0.0, kMoessbauerDummy, kMoessbauerDummy, kMoessbauerDummy, kMoessbauerDummy, kMoessbauerDummy, 0.0);
static Material* const kAirMaterial = new Material(kAirComposition, kDensity, kNone, kAbundance, kLambMoessbauer, {}, "air");
static Var* const kThickness = new Var(0.0, 0.0, 0.0, false);
static Var* const kRoughness = new Var(0.0, 0.0, 0.0, false);
diff --git a/src/nexus/lib/moessbauer.py b/src/nexus/lib/moessbauer.py
index e06cc6e9..df8e9c59 100644
--- a/src/nexus/lib/moessbauer.py
+++ b/src/nexus/lib/moessbauer.py
@@ -49,7 +49,9 @@ none = clib.MoessbauerIsotope(
quadrupole_ground = 0.0,
quadrupole_excited = 0.0,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 0.0
)
r"""
none
@@ -79,7 +81,9 @@ K40 = clib.MoessbauerIsotope(
quadrupole_ground = -0.061,
quadrupole_excited = 1.0e-200, # unkwonwn
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 1.17e-4
)
r"""
K-40
@@ -110,7 +114,9 @@ Sc45 = clib.MoessbauerIsotope(
quadrupole_ground = -0.22,
quadrupole_excited = 0.318,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 1.0
)
r"""
Sc-45
@@ -142,7 +148,9 @@ Fe57 = clib.MoessbauerIsotope(
quadrupole_ground = 0.0,
quadrupole_excited = 0.187,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 0.02119
)
r"""
Fe-57
@@ -174,7 +182,9 @@ Ni61 = clib.MoessbauerIsotope(
quadrupole_ground = 0.162,
quadrupole_excited = -0.2,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 0.011399
)
r"""
Ni-61
@@ -205,7 +215,9 @@ Zn67 = clib.MoessbauerIsotope(
quadrupole_ground = 0.15,
quadrupole_excited = 0.0,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 0.0404
)
r"""
Zn-67
@@ -238,7 +250,9 @@ Ge73 = clib.MoessbauerIsotope(
quadrupole_ground = -0.173,
quadrupole_excited = -0.4,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 0.0776
)
r"""
Ge-73
@@ -268,7 +282,9 @@ Kr83 = clib.MoessbauerIsotope(
quadrupole_ground = 0.253,
quadrupole_excited = 0.495,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 0.115
)
r"""
Kr-83
@@ -299,7 +315,9 @@ Ru99 = clib.MoessbauerIsotope(
quadrupole_ground = 0.079,
quadrupole_excited = 0.231,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 0.1276
)
r"""
Ru-99
@@ -330,7 +348,9 @@ Sn119 = clib.MoessbauerIsotope(
quadrupole_ground = 0.0,
quadrupole_excited = -0.094,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 0.0859
)
r"""
Sn-119
@@ -360,7 +380,9 @@ Sb121 = clib.MoessbauerIsotope(
quadrupole_ground = -0.36,
quadrupole_excited = -0.48,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 0.5721
)
r"""
Sb-121
@@ -392,7 +414,9 @@ Te125 = clib.MoessbauerIsotope(
quadrupole_ground = 0.0,
quadrupole_excited = -0.31,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 0.0707
)
r"""
Te-125
@@ -425,7 +449,9 @@ I127 = clib.MoessbauerIsotope(
quadrupole_ground = -0.79,
quadrupole_excited = -0.71,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 1.0
)
r"""
I-127
@@ -457,7 +483,9 @@ Sm149 = clib.MoessbauerIsotope(
quadrupole_ground = 0.075,
quadrupole_excited = 1.01,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 0.1382
)
r"""
Sm-149
@@ -488,7 +516,9 @@ Eu151 = clib.MoessbauerIsotope(
quadrupole_ground = 0.903,
quadrupole_excited = 1.28,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 0.4781
)
r"""
Eu-151
@@ -518,7 +548,9 @@ Gd157 = clib.MoessbauerIsotope(
quadrupole_ground = 1.38,
quadrupole_excited = 2.46,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 0.1565
)
r"""
Gd-157
@@ -549,7 +581,9 @@ Dy161 = clib.MoessbauerIsotope(
quadrupole_ground = 2.507,
quadrupole_excited = 2.506,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 0.18889
)
r"""
Dy-161
@@ -580,7 +614,9 @@ Tm169 = clib.MoessbauerIsotope(
quadrupole_ground = 0.0,
quadrupole_excited = -1.3,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 1.0
)
r"""
Tm-169
@@ -609,7 +645,9 @@ Yb172 = clib.MoessbauerIsotope(
quadrupole_ground = 0.0,
quadrupole_excited = 2.16,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 0.21686
)
r"""
Yb-172
@@ -641,7 +679,9 @@ Yb174 = clib.MoessbauerIsotope(
quadrupole_ground = 0.0,
quadrupole_excited = 2.12,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 0.32025
)
r"""
Yb-174
@@ -673,7 +713,9 @@ Ta181 = clib.MoessbauerIsotope(
quadrupole_ground = 3.17,
quadrupole_excited = 3.71,
- interference_term = -0.08 # Tramell and Hannon, Phys.Rev. 180, 337 (1969)
+ interference_term = -0.08, # Tramell and Hannon, Phys.Rev. 180, 337 (1969)
+
+ natural_abundance = 0.9998799
)
r"""
Ta-181
@@ -703,7 +745,9 @@ Ir193 = clib.MoessbauerIsotope(
quadrupole_ground = 0.78,
quadrupole_excited = 0.0,
- interference_term = 0.0
+ interference_term = 0.0,
+
+ natural_abundance = 0.627
)
r"""
Ir-193
--
GitLab
From fd2e6e6c18b0185e9c078665ade6607c6b48f0ba Mon Sep 17 00:00:00 2001
From: lbocklag <lars.bocklage@desy.de>
Date: Fri, 13 Dec 2024 12:56:25 +0100
Subject: [PATCH 2/2] updated Swig file
---
src/nexus/clib/classes/moessbauer.i | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/src/nexus/clib/classes/moessbauer.i b/src/nexus/clib/classes/moessbauer.i
index 2a5d5ff9..a98d1403 100644
--- a/src/nexus/clib/classes/moessbauer.i
+++ b/src/nexus/clib/classes/moessbauer.i
@@ -107,7 +107,7 @@
if attr not in ("this", "_internal_conversion", "_mixing_ratio_E2M1", "_gfactor_ground", "_gfactor_excited", "_quadrupole_ground", "_quadrupole_excited", "_interference_term") and not hasattr(self, attr):
raise Exception("MoessbauerIsotope has no attribute {}".format(attr))
- if attr in ("internal_conversion", "mixing_ratio_E2M1", "gfactor_ground", "gfactor_excited", "quadrupole_ground", "quadrupole_excited", "interference_term") and isinstance(val, (int, float)):
+ if attr in ("internal_conversion", "mixing_ratio_E2M1", "gfactor_ground", "gfactor_excited", "quadrupole_ground", "quadrupole_excited", "interference_term", "natural_abundance") and isinstance(val, (int, float)):
setattr(getattr(self, attr), "value", val)
return
--
GitLab