sciforge_lib/physics/nuclear/
nuclide.rs

1use std::fmt;
2
3pub const PROTON_MASS: f64 = crate::constants::PROTON_MASS_AMU;
4pub const NEUTRON_MASS: f64 = crate::constants::NEUTRON_MASS_AMU;
5pub const ELECTRON_MASS: f64 = crate::constants::ELECTRON_MASS_AMU;
6pub const AMU_TO_MEV: f64 = crate::constants::AMU_TO_MEV;
7pub const BOLTZMANN_KEV: f64 = crate::constants::KELVIN_TO_KEV;
8
9#[derive(Clone, Debug)]
10pub struct Nuclide {
11    pub z: u32,
12    pub a: u32,
13    pub mass_excess_mev: f64,
14    pub binding_energy_per_nucleon: f64,
15    pub half_life_seconds: Option<f64>,
16    pub name: &'static str,
17}
18
19impl Nuclide {
20    pub fn n(&self) -> u32 {
21        self.a - self.z
22    }
23
24    pub fn atomic_mass_amu(&self) -> f64 {
25        self.z as f64 * PROTON_MASS + self.n() as f64 * NEUTRON_MASS
26            - self.binding_energy() / AMU_TO_MEV
27    }
28
29    pub fn binding_energy(&self) -> f64 {
30        self.binding_energy_per_nucleon * self.a as f64
31    }
32
33    pub fn is_stable(&self) -> bool {
34        self.half_life_seconds.is_none()
35    }
36
37    pub fn decay_constant(&self) -> Option<f64> {
38        self.half_life_seconds.map(|t| core::f64::consts::LN_2 / t)
39    }
40
41    pub fn activity_bq(&self, num_atoms: f64) -> Option<f64> {
42        self.decay_constant().map(|l| l * num_atoms)
43    }
44}
45
46impl fmt::Display for Nuclide {
47    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
48        write!(
49            f,
50            "{}(Z={}, A={}, B/A={:.3} MeV)",
51            self.name, self.z, self.a, self.binding_energy_per_nucleon
52        )
53    }
54}
55
56pub fn semi_empirical_mass(z: u32, a: u32) -> f64 {
57    let a_f = a as f64;
58    let z_f = z as f64;
59    let n_f = (a - z) as f64;
60
61    let a_v = 15.56;
62    let a_s = 17.23;
63    let a_c = 0.7;
64    let a_a = 23.285;
65    let a_p = 12.0;
66
67    let volume = a_v * a_f;
68    let surface = -a_s * a_f.powf(2.0 / 3.0);
69    let coulomb = -a_c * z_f * (z_f - 1.0) / a_f.powf(1.0 / 3.0);
70    let asymmetry = -a_a * (n_f - z_f).powi(2) / a_f;
71    let pairing = if !a.is_multiple_of(2) {
72        0.0
73    } else if z.is_multiple_of(2) {
74        a_p / a_f.powf(0.5)
75    } else {
76        -a_p / a_f.powf(0.5)
77    };
78
79    volume + surface + coulomb + asymmetry + pairing
80}
81
82pub fn binding_energy_per_nucleon_semf(z: u32, a: u32) -> f64 {
83    semi_empirical_mass(z, a) / a as f64
84}
85
86pub fn nuclear_radius_fm(a: u32) -> f64 {
87    1.2 * (a as f64).powf(1.0 / 3.0)
88}
89
90pub fn nuclear_density_kg_m3() -> f64 {
91    2.3e17
92}
93
94pub fn separation_energy_proton(z: u32, a: u32) -> f64 {
95    if z == 0 || a <= 1 {
96        return 0.0;
97    }
98    semi_empirical_mass(z - 1, a - 1) + 7.289 - semi_empirical_mass(z, a)
99}
100
101pub fn separation_energy_neutron(z: u32, a: u32) -> f64 {
102    if a <= 1 {
103        return 0.0;
104    }
105    semi_empirical_mass(z, a - 1) + 8.071 - semi_empirical_mass(z, a)
106}
107
108pub fn separation_energy_alpha(z: u32, a: u32) -> f64 {
109    if z < 2 || a < 4 {
110        return 0.0;
111    }
112    semi_empirical_mass(z - 2, a - 4) + 28.296 - semi_empirical_mass(z, a)
113}
114
115pub fn valley_of_stability_z(a: u32) -> f64 {
116    let a_f = a as f64;
117    a_f / (2.0 + 0.0155 * a_f.powf(2.0 / 3.0))
118}
119
120pub fn neutron_excess(z: u32, a: u32) -> i32 {
121    a as i32 - 2 * z as i32
122}
123
124pub fn isospin_z(z: u32, a: u32) -> f64 {
125    (a as f64 - 2.0 * z as f64) / 2.0
126}
127
128pub fn liquid_drop_fission_parameter(z: u32, a: u32) -> f64 {
129    let z_f = z as f64;
130    let a_f = a as f64;
131    z_f * z_f / (50.88 * a_f)
132}
133
134pub fn fission_barrier_estimate_mev(z: u32, a: u32) -> f64 {
135    let x = liquid_drop_fission_parameter(z, a);
136    if x >= 1.0 {
137        return 0.0;
138    }
139    0.38 * (0.75 - x) * semi_empirical_mass(z, a).abs()
140}
141
142pub fn nuclear_skin_thickness_fm(z: u32, a: u32) -> f64 {
143    let n = a - z;
144    0.9 * (n as f64 - z as f64) / a as f64
145}
146
147pub fn weizsacker_mass_excess_mev(z: u32, a: u32) -> f64 {
148    let z_f = z as f64;
149    let n_f = (a - z) as f64;
150    z_f * 7.289 + n_f * 8.071 - semi_empirical_mass(z, a)
151}
152
153pub fn proton_drip_line_a(z: u32) -> u32 {
154    let mut best_a = 2 * z;
155    let mut best_sep = separation_energy_proton(z, best_a);
156    for a in (z + 1)..=(3 * z) {
157        let sep = separation_energy_proton(z, a);
158        if sep > best_sep {
159            best_sep = sep;
160            best_a = a;
161        }
162    }
163    best_a
164}
165
166pub fn neutron_drip_line_a(z: u32) -> u32 {
167    let mut a = 2 * z;
168    loop {
169        let sep = separation_energy_neutron(z, a);
170        if sep < 0.0 || a > 4 * z + 50 {
171            return a;
172        }
173        a += 1;
174    }
175}
176
177pub fn magic_number_nearest(n: u32) -> u32 {
178    let magic = [2, 8, 20, 28, 50, 82, 126, 184];
179    let mut best = magic[0];
180    let mut best_dist = n.abs_diff(best);
181    for &m in &magic[1..] {
182        let d = n.abs_diff(m);
183        if d < best_dist {
184            best_dist = d;
185            best = m;
186        }
187    }
188    best
189}
190
191pub fn is_doubly_magic(z: u32, n: u32) -> bool {
192    let magic = [2, 8, 20, 28, 50, 82, 126];
193    magic.contains(&z) && magic.contains(&n)
194}
sciforge_lib/physics/nuclear/nuclide.rs

sciforge_lib/physics/nuclear/
nuclide.rs
