Nuclide/
mmodel.rs

1/*
2    Mass model computing the binding energy
3    
4    Translated from Duflo-Zuker's 10-parameter Fortran 90 code
5    
6    In: total nucleons, number of protons
7    Out: Binding energy in MeV
8*/
9
10pub(crate) fn mass_model(a: usize, z: usize) -> f64 {
11    
12  const COEF : [f64;10] = [
13        0.7043,17.7418,16.2562,37.5562,53.9017,
14        0.4711,2.1307, 0.0210,40.5356,6.0632
15        ];
16        
17  let mut y : [f64;2] = [0f64,0f64];
18  let mut oei : [f64;2] = [0f64,0f64];
19  let mut dei : [f64;2] = [0f64,0f64];
20  let mut qx : [f64;2] = [0f64,0f64];
21  let mut dx : [f64;2] = [0f64,0f64];    
22  let mut pp : [f64;2] = [0f64,0f64];     
23  let mut n2 : [f64;2] = [0f64,0f64];                           
24  let mut op : [f64;2] = [0f64,0f64]; 
25  let mut onp :  [[[f64;9];2];2] = [[[0f64;9];2];2];
26  
27  let n = a-z;
28  let nuclei : [usize;2] = [n,z]; 
29  
30  
31  let a = n+z;
32  let t = n.abs_diff(z);
33  let r = (a as f64).cbrt();
34  let rc = r * (1.0-0.25*(t as f64/a as f64).powi(2));
35  let ra = rc.powi(2)/r;
36  let z2 = (z*(z-1)) as f64;
37  let initial_term = (-z2 + 0.76 * z2.powf(2.0/3.0))/rc;
38  for deformed in 0..2{
39  
40   let mut ju = 0f64;
41   
42    if deformed == 1{
43      ju = 4f64
44    }
45    
46    let mut term : [f64;10] = [0f64;10];
47    let mut noc : [[f64;18];2] = [[0f64;18];2];
48    let mut os : [f64;2] = [0f64,0f64]; 
49                        
50    term[0] = initial_term;
51    
52    for idx in 0..2{
53      n2[idx] = (2*(nuclei[idx]/2)) as f64;
54      let mut ncum = 0i64;
55      let mut i = 0i64;
56      let mut idd = 0i64;
57
58        loop{
59       
60         i+=1;
61         if i%2==1{
62            idd = i+1;
63         }
64         else if i%2==0{
65           idd = (i*(i-2))/4;
66         
67         }
68         ncum+=idd;
69       
70            if ncum < nuclei[idx] as i64{
71               noc[idx][(i-1) as usize] = idd as f64;
72            }
73            else if ncum >= nuclei[idx] as i64{
74           
75               break;
76            }
77         }
78         
79       let imax = i+1;
80       let mut ip = (i-1)/2;
81       let ipm = i/2;
82       pp[idx] = ip as f64;
83       let moc : f64 = nuclei[idx] as f64 - ncum as f64 + idd as f64;
84       noc[idx][(i-1) as usize] = moc-ju;
85       noc[idx][i as usize] = ju;
86                           
87       if i&1 == 1{
88         oei[idx] = (moc as i64+ (ip as i64*(ip as i64-1i64))) as f64;
89                
90         dei[idx] = (ip*(ip+1)+2) as f64;
91       }
92       else if i&1 == 0{
93         oei[idx] = ((moc as i64)-(ju as i64)) as f64;
94         dei[idx] = ((ip+1)*(ip+2)+2) as f64;
95       }
96         
97     
98       qx[idx] = oei[idx]*(dei[idx]-oei[idx]-ju)/dei[idx];
99       
100       dx[idx] = qx[idx]*(2.0*oei[idx]-dei[idx]);
101       if deformed == 1{
102         qx[idx] = qx[idx]/dei[idx].sqrt();
103       }
104       
105       if deformed == 1{
106         for i in 1..(imax+1){
107           let ip2 = (i-1)/2;
108           onp[idx][0][ip2 as usize]=0f64;
109           onp[idx][1][ip2 as usize]=0f64;
110         }
111     
112       }
113       
114      
115       for i in 1..(imax+1){
116         ip = (i-1)/2;
117         let fact = (((ip+1)*(ip+2)) as f64).sqrt();
118
119         onp[idx][0][ip as usize]+=noc[idx][(i-1) as usize]/fact;
120   
121         let mut vm = -1.0;
122         if i&1==1{
123           vm=0.5* ip as f64;
124         }
125         onp[idx][1][ip as usize]+=noc[idx][(i-1) as usize]*vm;
126       }
127
128      op[idx] = 0f64;
129      os[idx] = 0f64; 
130       for ip in 0..(ipm+1){
131         let den = (((ip+1)*(ip+2)) as f64).powf(3.0/2.0);
132         let one = onp[idx][1][ip as usize];
133         let zero = onp[idx][0][ip as usize];
134       
135         op[idx] = op[idx]+zero;
136      
137         os[idx] = os[idx] + one* (1.0+zero)* ((ip*ip) as f64/den) 
138                   + one*(1.0-zero)*((4i64*ip as i64-5i64) as f64/den);
139            
140       }
141    
142       op[idx] = op[idx]*op[idx];
143       } // end loop
144      
145       term[1] = op[0]+op[1];
146       
147       term[2] = -term[1]/ra;
148       term[1] += os[0]+os[1];
149       term[3] = -(t as f64)*(t as f64+2.0)/r.powi(2);
150       term[4] = -term[3]/ra;
151       
152       if deformed == 0{
153          term[5] = dx[0]+dx[1];
154          term[6] = -term[5]/ra;
155          let px = pp[0].sqrt()+pp[1].sqrt();
156          term[7] = qx[0] * qx[1]*(2f64.powf(px));
157       }
158      else if deformed == 1{
159        term[8] = qx[0]*qx[1];
160      }
161      term[4] = t as f64*(1.0-t as f64)/(a as f64*ra.powi(3)) + term[4];
162      
163    
164      
165      if(n2[0] != nuclei[0] as f64) && (n2[1] != nuclei[1] as f64){ term[9]= t as f64/a as f64;}                        
166      if n > z{                                                          
167        if(n2[0] == nuclei[0] as f64)&&(n2[1] != nuclei[1] as f64){term[9]= 1.0-t as f64/a as f64;}                    
168        if(n2[0] != nuclei[0] as f64)&&(n2[1]==nuclei[1] as f64){term[9]= 1.0;}
169        }                        
170      else {                                                                     
171        if(n2[0]==nuclei[0] as f64)&&(n2[1] != nuclei[1] as f64){term[9]= 1.0;}                        
172        if(n2[0]!=nuclei[0] as f64)&&(n2[1]==nuclei[1] as f64){term[9]= 1.0-t as f64/a as f64;}                    
173      }                                                                     
174      if(n2[1]==nuclei[1] as f64) && (n2[0] == nuclei[0] as f64){term[9]= 2.0-t as f64/a as f64;}            
175      
176     
177      for i in term[1..].iter_mut(){
178         *i=*i/ra;
179      }
180      for i in 0..10{
181         y[deformed]+=term[i]*COEF[i];     
182       }
183    }// end loop
184    
185    let de = y[1]-y[0];
186    if de > 0f64 || z < 50{
187      return y[0] 
188    }  
189    y[1]
190  }
Nuclide/mmodel.rs

Nuclide/
mmodel.rs
