Script to calculate nuclear Q-values from AME2012 (Atomic Mass Evaluation) data.

calcq.py

import numpy as np                                                                                                                                                                                                                                                            

def get_ame_masses(path):                                                                                                                                                                                                                                                     
    masses = np.zeros((400,400))                                                                                                                                                                                                                                              
    for line in open(path):                                                                                                                                                                                                                                                   
        if line[0]=='#':                                                                                                                                                                                                                                                      
            continue                                                                                                                                                                                                                                                          
        line = line.split()                                                                                                                                                                                                                                                   
        #print line                                                                                                                                                                                                                                                           
        if '-' in line[5]:                                                                                                                                                                                                                                                    
            continue                                                                                                                                                                                                                                                          
        if line[5][-1]=='#':                                                                                                                                                                                                                                                  
            line[5] = line[5][:-1]                                                                                                                                                                                                                                            
        A = int(line[3])                                                                                                                                                                                                                                                      
        Z = int(line[2])                                                                                                                                                                                                                                                      
        mass = 0.931494028*1.0e-6*float(line[5])                                                                                                                                                                                                                              
        masses[A,Z]=mass                                                                                                                                                                                                                                                      
    return masses                                                                                                                                                                                                                                                             
                                                                                                                                                                                                                                                                              
def get_ame_q(A,Z):                                                                                                                                                                                                                                                           
    m1 = mass_table[A,Z]                                                                                                                                                                                                                                                      
    m2 = mass_table[A,Z-1]                                                                                                                                                                                                                                                    
    if m1 == 0.0 or m2 == 0.0:                                                                                                                                                                                                                                                
        return False                                                                                                                                                                                                                                                          
    else:                                                                                                                                                                                                                                                                     
        return (m1 - m2 - 0.000511)*1000                                                                                                                                                                                                                                      

if __name__=="__main__":                                                                                                                                                                                                                                                      
    mass_table = get_ame_masses('./mass.mas12')                                                                                                                                                                                                                               
    print get_ame_q(20,10)