Nuclide/nstruct/
isotope.rs

1use crate::traits::{ChemElement,Isotope};
2use crate::nstruct::core::Nuclide;
3
4use crate::constant::*;
5use crate::mmodel::mass_model;
6use crate::decay::DecayMode;
7use crate::nuclidedata::half_life::HALF_LIFE;
8use crate::nuclidedata::decay_chain::DECAY_CHAIN;
9use crate::Particle;
10use crate::rng::rand;
11
12use crate::decay::is_mode_recorded;
13use crate::decay::dhelper::decay_mode_idx;
14use crate::decay::dhelper::decay_select;
15use crate::decay::decayimpl::decayindex;
16
17
18
19
20
21impl Isotope for Nuclide {
22
23    
24    fn mass_deficit(&self) -> f64 {
25        let (p,n) = self.proton_neutron();
26        ( PROTONMASS* p as f64
27            + ELECTRONMASS * p as f64
28            + NEUTRONMASS * n as f64)
29            - self.am()
30    }
31    
32    fn binding_energy(&self) -> f64 {
33        let (z, a) = (
34            self.proton_neutron().0,
35            self.proton_neutron().0 + self.proton_neutron().1,
36        );
37        mass_model(a, z)
38    }
39    
40   
41   fn branching_ratio<T: DecayMode>(&self) -> f64{
42      let idx = self.nuclide_index() * 6 + 5; 
43      match decayindex::<T>(idx){
44        Some(x) => x as f64/FLOAT_64,
45        None => f64::NAN,
46      }
47    }
48    
49    fn half_life<T: DecayMode>(&self) -> f64 {
50       let branch = self.branching_ratio::<T>();
51       HALF_LIFE[self.nuclide_index()]/branch
52       
53    }
54    
55    fn mean_lifetime<T: DecayMode>(&self) -> f64 {
56        //reciprocal of ln(2) average lifespan of a particle
57        self.half_life::<T>() * std::f64::consts::LOG2_E
58    }
59
60    /// Approximation of decay constant
61    fn decay_constant<T: DecayMode>(&self) -> f64 {
62        self.mean_lifetime::<T>().recip()
63    }
64    
65    
66    /// Returns probability of decay 
67    fn decay_probability<T: DecayMode>(&self, time: f64) -> f64{
68           1.0 - (-self.decay_constant::<T>() * time).exp()
69    }
70    
71    fn daughter_theoretical<T: DecayMode>(&self) -> Option<Self>{
72    	T::decay_result(self)
73    }
74
75    //lowest probability is 1/u64::MAX
76    ///Returns true if the nuclide would have decayed in the time given. The nuclide remains unchanged
77    fn decay_time<T: DecayMode>(&self, time: f64) -> bool {
78        let prob =
79            ((1.0 - (-self.decay_constant::<T>() * time).exp()) * FLOAT_64) as u64;
80
81        prob > rand()
82    }
83
84    /// Returns the probable decay modes  as a string
85    fn decay_string(&self) -> String {
86        let mut unu_alea =
87            ((DECAY_CHAIN[self.nuclide_index() * 6] as f64 / FLOAT_64) * 100.0)
88                .to_string();
89        unu_alea.truncate(4);
90        unu_alea.push_str("% ");
91        let mut doua_alea = ((DECAY_CHAIN[self.nuclide_index() * 6 + 1] as f64
92            / FLOAT_64)
93            * 100.0)
94            .to_string();
95        doua_alea.truncate(4);
96        doua_alea.push_str("% ");
97        let mut trei_alea = ((DECAY_CHAIN[self.nuclide_index() * 6 + 2] as f64
98            / FLOAT_64)
99            * 100.0)
100            .to_string();
101        trei_alea.truncate(4);
102        trei_alea.push_str("% ");
103        let mut patru_alea = ((DECAY_CHAIN[self.nuclide_index() * 6 + 3] as f64
104            / FLOAT_64)
105            * 100.0)
106            .to_string();
107        patru_alea.truncate(4);
108        patru_alea.push_str("% ");
109        let mut cinci_alea = ((DECAY_CHAIN[self.nuclide_index() * 6 + 4] as f64
110            / FLOAT_64)
111            * 100.0)
112            .to_string();
113        cinci_alea.truncate(4);
114        cinci_alea.push_str("% ");
115
116        let mut decay_str = vec![];
117
118        let decay_vector = DECAY_CHAIN[self.nuclide_index() * 6 + 5].to_be_bytes();
119
120        for i in decay_vector[..5].iter() {
121            match i {
122                1 => decay_str.push("α; "),
123                2 => decay_str.push("p; "),
124                3 => decay_str.push("2p; "),
125                4 => decay_str.push("n; "),
126                5 => decay_str.push("2n; "),
127                6 => decay_str.push("EC; "),
128                7 => decay_str.push("2EC; "),
129                8 => decay_str.push("β− + p; "),
130                9 => decay_str.push("β+; "),
131                10 => decay_str.push("2β+; "),
132                11 => decay_str.push("β−; "),
133                12 => decay_str.push("2β−; "),
134                13 => decay_str.push("β− + n; "),
135                14 => decay_str.push("β− + 2n; "),
136                15 => decay_str.push("β− + 3n; "),
137                16 => decay_str.push("β+ + p; "),
138                17 => decay_str.push("β+ + 2p; "),
139                18 => decay_str.push("β+ + 3p; "),
140                19 => decay_str.push("β− + α; "),
141                20 => decay_str.push("β+ + α; "),
142                21 => decay_str.push("β− + d; "),
143                22 => decay_str.push("β− + t; "),
144                23 => decay_str.push("SF; "),
145                24 => decay_str.push("β− + SF; "),
146                25 => decay_str.push("β+ + SF; "),
147                26 => decay_str.push("C-14; "),
148                27 => decay_str.push("Ne-20; "),
149                28 => decay_str.push("Ne-24; "),
150                29 => decay_str.push("Ne-20 + Ne-24; "),
151                30 => decay_str.push("Si-32; "),
152                31 => decay_str.push("Si-34; "),
153
154                _ => decay_str.push("Null"),
155            }
156        }
157        let mut decayvec = vec![];
158        decayvec.push(unu_alea);
159        decayvec.push(decay_str[0].to_string());
160        decayvec.push(doua_alea);
161        decayvec.push(decay_str[1].to_string());
162        decayvec.push(trei_alea);
163        decayvec.push(decay_str[2].to_string());
164        decayvec.push(patru_alea);
165        decayvec.push(decay_str[3].to_string());
166        decayvec.push(cinci_alea);
167        decayvec.push(decay_str[4].to_string());
168
169        if decayvec[0] == "0% " {
170            "Stable".to_string()
171        } else {
172            match decayvec.iter().position(|r| r == "Null") {
173                Some(x) => decayvec.truncate(x - 1),
174                None => decayvec.truncate(10),
175            }
176            decayvec.join("")
177        }
178    }
179    
180    fn daughter_energetic<T: DecayMode>(&mut self) -> (f64,Vec<Particle>){
181        T::decay(self)
182    }
183    
184    
185    fn daughter<T: DecayMode>(&self) -> Option<Self>{
186        if is_mode_recorded::<T>(self){
187          return T::decay_result(self)
188        }
189        None
190    }
191    
192    /// Q-value (total energy) of a nuclear decay, regardless of whether it is observed
193    /// # NAN
194    /// Returns NAN if this decay mode results in a nonexistent nuclide
195    fn decay_q<T: DecayMode>(&self) -> f64{
196        T::q_value(self)
197    } 
198    
199    /**
200    Returns the name and isotope number of the nuclide
201
202       ```
203       use ::Nuclide::{Nuclide, Isotope,::decay::TotalDecay};
204       let mut uranium = "U-238".parse::<Nuclide>().unwrap();
205
206      // total particle energy of the nuclide and vector of decay particles
207       let (particle_energy,particle_vector) = uranium.decay::<TotalDecay>(5E+20);
208
209        // final nuclide in the U-238 decay series
210        assert_eq!(uranium.to_string(), "Pb-206");
211       ```
212    */
213    fn decay<T: DecayMode>(&mut self, mut time: f64) -> (f64, Vec<Particle>) {
214        let mut total_energy = 0f64;
215        let mut particlevec = vec![];
216        
217        if T::decay_index() == 254u8{
218       
219        loop {
220        
221        // Obtain index of decay branch
222        let (b,b_idx) = decay_mode_idx(self);
223       
224        let idx = self.nuclide_index() * 6usize + b_idx as usize; 
225
226        let ratio = (DECAY_CHAIN[idx] as f64)/FLOAT_64;
227        
228        let mean_time = (HALF_LIFE[self.nuclide_index()]/ratio) * std::f64::consts::LOG2_E; 
229
230        if mean_time >= time || mean_time == f64::INFINITY{
231          return (total_energy,particlevec)
232        }
233        
234        let (energy, p_vector) = decay_select(self,b);
235        particlevec.extend_from_slice(&p_vector[..]);
236        total_energy+=energy;
237        time-=mean_time;
238        
239        }
240        
241        } // end if
242      
243        loop{
244         
245             if !is_mode_recorded::<T>(self){
246                return (total_energy,particlevec)
247             }
248      
249             let mean_time = self.mean_lifetime::<T>();
250             
251             if mean_time >= time || mean_time.is_infinite() || mean_time.is_nan(){
252                return (total_energy,particlevec)
253             }
254             
255             let (energy, p_vector) = T::decay(self);
256          
257             particlevec.extend_from_slice(&p_vector[..]);
258             total_energy+=energy;
259             
260             time-=mean_time;
261          } 
262    }
263    
264 }
Nuclide/nstruct/isotope.rs

Nuclide/nstruct/
isotope.rs
