def prodfac_explosion(self,mp,sefiles,sefiles_hout,isotopes,ini_abu=""):
'''
Marcos routine adapted, calc SN yields
'''
import utils as u
initial_abu=iniabu(ini_abu) #,iniabufile='../../frames/mppnp/USEEPP/"
cycle=sefiles_hout.se.cycles[-1]
first_cycle=sefiles_hout.se.cycles[0]
#mass_cut=1.60
###identify mcut at mass cell where rho >1 at first cycle,from core to surface
rho_first_cycle=sefiles_hout.get(first_cycle,"rho")
mass_first_cycle=sefiles_hout.get(first_cycle,"mass")
print rho_first_cycle
print cycle
print sefiles_hout.get("mini")
print sefiles_hout.get("zini")
idx_mcut=int(np.where(rho_first_cycle!=1)[0][0])
mass_cut=mass_first_cycle[idx_mcut]
# I take the max mass coord
mass_limit=mass_first_cycle[-1]
# I calculate average_mass_frac_decay
sefiles_hout.average_iso_abund_marco([mass_cut,mass_limit],cycle,True,2)
E_i_pre_15=array(mp.average_mass_frac_decay)*(mass_limit-mass_cut)
#follwoing needed for isotope identification,must be stable
# define back_ind presn
back_ind_pre=u.back_ind
yields=[]
iniabus=[]
prod_facs=[]
for i in range(len(isotopes)):
other_name_scheme=isotopes[i].split("-")[0].upper()+(5-len(isotopes[i])+1)*" "+isotopes[i].split("-")[1]
if other_name_scheme not in u.stable:
print "error: Chosen isotope is not stable: ",other_name_scheme
return 1
idx_iso=back_ind_pre[other_name_scheme]
yield_decayed=E_i_pre_15[idx_iso]
yields.append(yield_decayed)
iniabu_1=initial_abu.abu[ initial_abu.names.index(isotopes[i].split("-")[0].lower()+(5-len(isotopes[i])+1)*" "+isotopes[i].split("-")[1])]
print "iniabu for exp",iniabu_1
iniabus.append(iniabu_1)
total_yield= yield_decayed
total_prod_factor= total_yield/ ((mass_limit-mass_cut)*iniabu_1)
prod_facs.append(total_prod_factor)
print total_yield, total_prod_factor
ini_star_mass=sefiles.get("mini")
iniabus=np.array(iniabus)*(ini_star_mass)#-end_star_mass)
return np.array(prod_facs),isotopes,yields,iniabus,mass_cut,first_cycle
评论列表
文章目录
