def get_path_from_sys(i):
slash = 0
path_list = []
while (slash < 7):
path_list.append(arg_list[i])
path = ''.join(path_list)
if path[-1:] != '/':
path_list.append(' ')
slash = path.count("/")
i+=1
return path, i
python类append()的实例源码
def get_datetime_from_sys(i):
time = []
time.append(arg_list[i])
time.append(' ')
time.append(arg_list[i+1])
datetime = ''.join(time)
i +=2
return datetime,i
# Initialize argument list from system command input
def get_dss_filename(i):
dss_ext = False
filename_list = []
if arg_list[i][-3:].lower() == 'dss':
filename = arg_list[i]
filename = ''.join(filename)
dss_ext = True
i+=1
while dss_ext == False:
filename_list.append(arg_list[i])
filename = ' '.join(filename_list)
if arg_list[i][-3:].lower() == 'dss':
dss_ext = True
i+=1
return filename, i
def get_path_from_sys(i):
slash = 0
path_list = []
while (slash < 7):
path_list.append(arg_list[i])
path = ''.join(path_list)
if path[-1:] != '/':
path_list.append(' ')
slash = path.count("/")
i+=1
return path, i
def get_datetime_from_sys(i):
time = []
time.append(arg_list[i])
time.append(' ')
time.append(arg_list[i+1])
datetime = ''.join(time)
i +=2
return datetime,i
def get_data_list(i):
data_count = int(arg_list[i])
i+=1
count = 0
data_list = []
while count < data_count:
data_list.append(float(arg_list[i]))
count+=1
i+=1
return data_list
# Initialize argument list from system command input
def format_time(timespan, precision=3):
"""Formats the timespan in a human readable form"""
if timespan >= 60.0:
# we have more than a minute, format that in a human readable form
parts = [("d", 60 * 60 * 24), ("h", 60 * 60), ("min", 60), ("s", 1)]
time = []
leftover = timespan
for suffix, length in parts:
value = int(leftover / length)
if value > 0:
leftover = leftover % length
time.append('{0}{1}'.format(str(value), suffix))
if leftover < 1:
break
return " ".join(time)
# Unfortunately the unicode 'micro' symbol can cause problems in
# certain terminals.
# See bug: https://bugs.launchpad.net/ipython/+bug/348466
# Try to prevent crashes by being more secure than it needs to
# E.g. eclipse is able to print a µ, but has no sys.stdout.encoding set.
units = ["s", "ms", 'us', "ns"] # the save value
if hasattr(sys.stdout, 'encoding') and sys.stdout.encoding:
try:
'\xb5'.encode(sys.stdout.encoding)
units = ["s", "ms", '\xb5s', "ns"]
except Exception:
pass
scaling = [1, 1e3, 1e6, 1e9]
if timespan > 0.0:
order = min(-int(math.floor(math.log10(timespan)) // 3), 3)
else:
order = 3
return "{1:.{0}g} {2}".format(precision, timespan * scaling[order],
units[order])
def _format_time(timespan, precision=3):
"""Formats the timespan in a human readable form"""
import math
if timespan >= 60.0:
# we have more than a minute, format that in a human readable form
# Idea from http://snipplr.com/view/5713/
parts = [("d", 60*60*24),("h", 60*60),("min", 60), ("s", 1)]
time = []
leftover = timespan
for suffix, length in parts:
value = int(leftover / length)
if value > 0:
leftover = leftover % length
time.append(u'%s%s' % (str(value), suffix))
if leftover < 1:
break
return " ".join(time)
# Unfortunately the unicode 'micro' symbol can cause problems in
# certain terminals.
# See bug: https://bugs.launchpad.net/ipython/+bug/348466
# Try to prevent crashes by being more secure than it needs to
# E.g. eclipse is able to print a µ, but has no sys.stdout.encoding set.
units = [u"s", u"ms",u'us',"ns"] # the save value
if hasattr(sys.stdout, 'encoding') and sys.stdout.encoding:
try:
u'\xb5'.encode(sys.stdout.encoding)
units = [u"s", u"ms",u'\xb5s',"ns"]
except:
pass
scaling = [1, 1e3, 1e6, 1e9]
if timespan > 0.0:
order = min(-int(math.floor(math.log10(timespan)) // 3), 3)
else:
order = 3
return u"%.*g %s" % (precision, timespan * scaling[order], units[order])
def _format_time(timespan, precision=3):
"""Formats the timespan in a human readable form"""
if timespan >= 60.0:
# we have more than a minute, format that in a human readable form
# Idea from http://snipplr.com/view/5713/
parts = [("d", 60*60*24),("h", 60*60),("min", 60), ("s", 1)]
time = []
leftover = timespan
for suffix, length in parts:
value = int(leftover / length)
if value > 0:
leftover = leftover % length
time.append(u'%s%s' % (str(value), suffix))
if leftover < 1:
break
return " ".join(time)
# Unfortunately the unicode 'micro' symbol can cause problems in
# certain terminals.
# See bug: https://bugs.launchpad.net/ipython/+bug/348466
# Try to prevent crashes by being more secure than it needs to
# E.g. eclipse is able to print a µ, but has no sys.stdout.encoding set.
units = [u"s", u"ms",u'us',"ns"] # the save value
if hasattr(sys.stdout, 'encoding') and sys.stdout.encoding:
try:
u'\xb5'.encode(sys.stdout.encoding)
units = [u"s", u"ms",u'\xb5s',"ns"]
except:
pass
scaling = [1, 1e3, 1e6, 1e9]
if timespan > 0.0:
order = min(-int(math.floor(math.log10(timespan)) // 3), 3)
else:
order = 3
return u"%.*g %s" % (precision, timespan * scaling[order], units[order])
def remnant_lifetime_agb(self):
'''
For paper1 extension:
bottom_envelope
Numbers of remnant mass shell masses, exists also in mesa_set + star age!
'''
inim=[]
remnm=[]
time11=[]
tottime=[]
c_core=[]
o_core=[]
small_co_core=[]
c_core_center=[]
for i in range(len(self.runs_H5_surf)):
m1p65_last=se(self.runs_H5_out[i])
mass_dummy=m1p65_last.se.get(m1p65_last.se.cycles[len(m1p65_last.se.cycles)-2],'mass')
top_of_envelope=mass_dummy[len(mass_dummy)-1]
h_dummy=m1p65_last.se.get(m1p65_last.se.cycles[len(m1p65_last.se.cycles)-2],'iso_massf','H-1')
c_dummy=m1p65_last.se.get(m1p65_last.se.cycles[len(m1p65_last.se.cycles)-2],'iso_massf','C-12')
o_dummy=m1p65_last.se.get(m1p65_last.se.cycles[len(m1p65_last.se.cycles)-2],'iso_massf','O-16')
for j in range(len(mass_dummy)):
if h_dummy[j] > 1e-1:
bottom_of_envelope = mass_dummy[j]
break
inim.append(m1p65_last.get("mini"))
remnm.append(bottom_of_envelope)
###Calculate the lifetime (MS)
sefiles=m1p65_last
cycs=[]
for k in range(5,len(sefiles.se.cycles),5):
cycs.append(int(sefiles.se.cycles[k]))
w=0
for cyc in cycs:
c12_center=sefiles.get(cyc,'C-12')[0]
#c12_center=c12[w][0]
w+=1
if c12_center>1e-1:
time1=(sefiles.get(cyc,'age')*sefiles.get('age_unit'))/31557600.
time11.append(time1)
break
tottime.append(sefiles.get(int(sefiles.se.cycles[-1]),'age')/31557600.)
print "M_initial | M_remn/bottom of envelope | total lifetime"
for i in range(len(inim)):
print inim[i],"|",'{:.3E}'.format(remnm[i]),"|",'{:.3E}'.format(tottime[i])
def set_plot_profile_decay(self,cycles=20*[-1],mass_range=20*[[0,0]],ylim=20*[[0,0]],isotopes=[],linestyle=[],save_dir=''):
'''
Plots HRDs
end_model - array, control how far in models a run is plottet, if -1 till end
symbs_1 - set symbols of runs
'''
if len(linestyle)==0:
linestyle=200*['-']
import nugridse as mp
import utils as u
#print self.runs_H5_restart
for i in range(len(self.runs_H5_restart)):
sefiles=mp.se(self.runs_H5_restart[i])
cycle=cycles[i]
if cycle==-1:
cycle=int(sefiles.se.cycles[-1])
if mass_range[i][0] ==0 and mass_range[i][1]==0:
mass_range[i][1]=sefiles.get(cycle,'mass')[-1]
sefiles.read_iso_abund_marco(mass_range[i],cycle)
u.stable_specie()
sefiles.decay(sefiles.mass_frac)
idx_species=[]
for k in range(len(isotopes)):
other_name_scheme=isotopes[k].split("-")[0].upper()+(5-len(isotopes[k])+1)*" "+isotopes[k].split("-")[1]
#other_name_scheme=other_name_scheme.capitalize()
idx_specie=u.back_ind[other_name_scheme]
idx_species.append(idx_specie)
mass_abu_array=[]
for idx_specie in idx_species:
mass_abu_array.append([])
for idx_mass in range(len(mp.decayed_multi_d)):
mass_abu_array[-1].append(mp.decayed_multi_d[idx_mass][idx_specie])
#plotting
plt.figure(self.run_dirs_name[i])
#print len(mp.used_masses),len(mass_abu_array[0])
#print mass_abu_array[0]
for k in range(len(isotopes)):
plt.plot(mp.used_masses,mass_abu_array[k],linestyle=linestyle[k],label=isotopes[k])
plt.legend()
plt.yscale('log')
#print sefiles.get(cycle,'mass')[-1]
plt.xlabel('M/Msun')
plt.ylabel('$X_i$')
plt.xlim(mass_range[i][0],mass_range[i][1])
if (ylim[i][0]>0 or ylim[i][1]>0) or (ylim[i][0]>0 and ylim[i][1]>0):
plt.ylim(ylim[i][0],ylim[i][1])
if len(save_dir)>0:
star_mass=sefiles.get("mini")
star_z=sefiles.get("zini")
plt.savefig(save_dir+'/'+self.run_dirs_name[i]+'_decay_profiles.png')
def set_pocket(self,runs=[],cycle=[45532,47566],mass_cell=[0.641981,0.641981],isotopes=[],x_charge=False,mass_number_range=[],decayed=False,iso_label=False,title="PDCZ",colors=["red","blue","green","black"],yax_log=True,filename_norm="iniab2.0E-02GN93.ppn",elem_norm=''):
'''
plots isotopic composition for different runs, each with specific cycle number cycle[i] and specific mass cell mass_cell[i].
Initial abundace file for normalization can be chosen as filename_norm, but MUST be in run directory
decayed - plot only decayed isotopes
mass_number_range - array of min mass number and max mass number, isotopes inbetween will be plottet, if set, isotope array will be ignored
iso_label - set labels of isotopes True or False
title -
colors -
filename_norm - not yet included, normalization with initial abundance file
....
e.g.
setse.set_pocket(cycle=[32840,29390],mass_cell=[0.6300,0.6076],isotopes=["He-4","C-12","O-16","Ne-22"],mass_number_range=[80,140],decayed=True)
'''
import nugridse as mp
sefiles=[]
legend=[]
HDF5_out=[]
extra_label=[]
if len(runs) ==0:
HDF5_out=self.runs_H5_out
runs=self.run_dirs_name
extra_label=self.extra_label
else:
for i in range(len(self.run_dirs_name)):
if self.run_dirs_name[i] in runs:
HDF5_out.append(self.runs_H5_out[i])
extra_label.append(self.extra_label[i])
for i in range(len(HDF5_out)):
reload(mp)
file_norm=HDF5_out[i][:-6]+filename_norm
sefiles=se(HDF5_out[i])
mass=sefiles.get("mini")
z=sefiles.get("zini")
legend=str(mass)+"M$_{\odot}$ Z= "+str(z)+", "+extra_label[i]
if len(isotopes)==0:
self.plot_abu_atmasscoord(mp,sefiles,cycle[i],mass_cell[i],["He-4","C-12","O-16","Ne-22"],x_charge,mass_number_range,decayed,file_norm,elem_norm,yax_log,label=iso_label,legend=legend,color=colors[i],title=title)
#plt.figure(111);self.pocket_composition(mp,sefiles=sefiles,cycle=cycle[i], massbot=massrange[i][0],masstop=massrange[i][1],isotopes=["He-4","C-12","O-16","Ne-22"],label=True,legend=legend,color=colors[i],title=title)
#plt.figure(222);self.pocket_composition(mp,sefiles=sefiles,cycle=cycle[i], massbot=massrange[i][0],masstop=massrange[i][1],isotopes=['Fe-56','Co-60','Ni-61','Cu-65','Zn-67','Ga-71','Ge-73','As-75','Se-77','Br-82','Kr-84','Rb-87','Sr-88','Y-89','Zr-92','Nb-94','Zr-96','Mo-96','Ba-137','Ba-138','La-139','Ce-140','Nd-142','Sm-144','Pb-206'],label=True,legend=legend,color=colors[i],title=title)
else:
#plt.figure(111);self.pocket_composition(mp,sefiles=sefiles,cycle=cycle[i], massbot=massrange[i][0],masstop=massrange[i][1],isotopes=isotopes,label=True,legend=legend,color=colors[i],title=title)
self.plot_abu_atmasscoord(mp,sefiles,cycle[i],mass_cell[i],isotopes,x_charge,mass_number_range,decayed,file_norm,elem_norm,yax_log,label=iso_label,legend=legend,color=colors[i],title=title)
def set_surface_abundances(self,runs=[],cycle=[-1,-1],isotopes=[],x_charge=True,mass_number_range=[],decayed=False,iso_label=False,title="",colors=["red","blue","green","black"],yax_log=True,filename_norm="iniab2.0E-02GN93.ppn",elem_norm=''):
'''
plots isotopic composition for different runs, each with specific cycle number cycle[i] and specific mass cell mass_cell[i].
Initial abundace file for normalization can be chosen as filename_norm, but MUST be in run directory
decayed - plot only decayed isotopes
mass_number_range - array of min mass number and max mass number, isotopes inbetween will be plottet, if set, isotope array will be ignored
iso_label - set labels of isotopes True or False
title -
colors -
filename_norm - not yet included, normalization with initial abundance file
....
e.g.
setse.set_pocket(cycle=[32840,29390],mass_cell=[0.6300,0.6076],isotopes=["He-4","C-12","O-16","Ne-22"],mass_number_range=[80,140],decayed=True)
'''
import nugridse as mp
sefiles=[]
legend=[]
HDF5_out=[]
extra_label=[]
if len(runs) ==0:
HDF5_out=self.runs_H5_out
runs=self.run_dirs_name
extra_label=self.extra_label
else:
for i in range(len(self.run_dirs_name)):
if self.run_dirs_name[i] in runs:
HDF5_out.append(self.runs_H5_out[i])
extra_label.append(self.extra_label[i])
for i in range(len(HDF5_out)):
reload(mp)
file_norm=HDF5_out[i][:-6]+filename_norm
sefiles=se(HDF5_out[i])
if cycle[i] ==-1:
cycle_1=sefiles.se.cycles[-1]
else:
cycle_1=cycle[i]
masscell=sefiles.get(cycle_1,"mass")[0]
mass=sefiles.get("mini")
z=sefiles.get("zini")
legend=str(mass)+"M$_{\odot}$ Z= "+str(z)+", "+extra_label[i]
if len(isotopes)==0:
self.plot_abu_atmasscoord(mp,sefiles,cycle_1,masscell,isotopes,x_charge,mass_number_range,decayed,file_norm,elem_norm,yax_log,label=iso_label,legend=legend,color=colors[i],title=title)
#plt.figure(111);self.pocket_composition(mp,sefiles=sefiles,cycle=cycle[i], massbot=massrange[i][0],masstop=massrange[i][1],isotopes=["He-4","C-12","O-16","Ne-22"],label=True,legend=legend,color=colors[i],title=title)
#plt.figure(222);self.pocket_composition(mp,sefiles=sefiles,cycle=cycle[i], massbot=massrange[i][0],masstop=massrange[i][1],isotopes=['Fe-56','Co-60','Ni-61','Cu-65','Zn-67','Ga-71','Ge-73','As-75','Se-77','Br-82','Kr-84','Rb-87','Sr-88','Y-89','Zr-92','Nb-94','Zr-96','Mo-96','Ba-137','Ba-138','La-139','Ce-140','Nd-142','Sm-144','Pb-206'],label=True,legend=legend,color=colors[i],title=title)
else:
#plt.figure(111);self.pocket_composition(mp,sefiles=sefiles,cycle=cycle[i], massbot=massrange[i][0],masstop=massrange[i][1],isotopes=isotopes,label=True,legend=legend,color=colors[i],title=title)
self.plot_abu_atmasscoord(mp,sefiles,cycle_1,masscell,isotopes,x_charge,mass_number_range,decayed,file_norm,elem_norm,yax_log,label=iso_label,legend=legend,color=colors[i],title=title)
def write_gce_input_wind_ejection_rate(self,file_name="isotopic_table.txt",final_models=[]):
'''
Calculates total mass ejection rate from stellar wind phase.
'''
e_kin_array=[]
mz=[]
m_final=[]
mz=[]
ini_m=[]
p=-1
for i in range(len(self.runs_H5_out)):
p+=1
sefiles=se(self.runs_H5_out[i])
sefiles_restart=se(self.runs_H5_restart[i])
mass=sefiles.get("mini")
metallicity=sefiles.get("zini")
mz1='(M='+str(round(mass,2))+',Z='+str(metallicity)+')'
endcycle=int(sefiles.se.cycles[-1])
if len(final_models)>0:
cycs=self.get_last_cycle([0,final_models[p]+500,500],sefiles,sefiles_restart)
else:
cycs=self.get_last_cycle([0,-1,500],sefiles,sefiles_restart)
endcycle=cycs[-1]
h_mass_frac=sefiles.get(endcycle,'H-1')
mass_array=sefiles.get(endcycle,'mass')
#earlier: 1.e-1 but now as in set1 scripts 0.05
h_free_core=mass_array[np.where(h_mass_frac<5.e-2)[0][-1]]
m_final.append(h_free_core)
mz.append(mz1)
ini_m.append(mass)
f1=open(file_name,'r')
lines=f1.readlines()
f1.close()
i=-1
line1=''
while (True):
i+=1
if i>len(lines)-1:
break
line=lines[i]
line1+=lines[i]
for k in range(len(mz)):
if mz[k] in lines[i]:
## check ahead for H Mfinal:
#for h in range(i,i+10):
# if 'H Lifetime:' in lines[h]:
# lifetime=float(lines[h].split(':')[1])
# break
# if h==i+9:
# print 'ERROR: no lifetime in table found!!'
# return
ejection=(ini_m[k]-m_final[k]) #/lifetime
line1+=('H Wind ejection: '+'{:.3E}'.format(ejection)+'\n')
break
f1=open(file_name,'w')
f1.write(line1)
f1.close()
def write_gce_input_lifetimes(self,file_name="isotopic_table.txt",final_cycles=[]):
'''
Calculates and writes lifetime in file_name
X_carbon_center > 0.4
'''
###Get lifetimes###
print '##############GETTING THE LIFETIMES##############'
time=[]
model=[]
mz=[]
for i in range(len(self.runs_H5_surf)):
sefiles=se(self.runs_H5_surf[i])
mass=sefiles.get("mini")
metallicity=sefiles.get("zini")
mz1='(M='+str(round(mass,2))+',Z='+str(metallicity)+')'
cycs=[]
if final_cycles[i]>-1:
final_cycle=final_cycles[i]
else:
final_cycle=int(sefiles.se.cycles[-1])
time1=(sefiles.get(final_cycle,'age')*sefiles.get('age_unit'))/31557600.
time.append(time1)
mz.append(mz1)
f1=open(file_name,'r')
lines=f1.readlines()
f1.close()
i=-1
line1=''
print 'Adding for ',mz,'lifetime : ',time
while (True):
i+=1
if i>len(lines)-1:
break
line=lines[i]
line1+=lines[i]
for k in range(len(mz)):
if mz[k] in lines[i]:
line1+=('H Lifetime: '+'{:.3E}'.format(time[k])+'\n')
break
f1=open(file_name,'w')
f1.write(line1)
f1.close()
def get_last_cycle(self,cycles,sefiles,sefiles_restart):
'''
Because there are a lot of issues with cycles we need this routine :(
'''
if cycles[1]==-1:
#if star_mass>9:
endcycle=int(sefiles_restart.se.cycles[-2])
#else:
# endcycle=int(sefiles[k].se.cycles[-2])
else:
endcycle=cycles[1]
cycs=range(cycles[0],endcycle,cycles[2])
if endcycle not in cycs:
cycs=cycs+[endcycle]
print 'first, last cycle: ',cycles[0],endcycle
star_mass11=sefiles.get(cycs,'mass')
#do a check for missing cycles
if not len(star_mass11)==len(cycs):
cycs1=[int(cycs[0])]
for cyc in cycs[1:]:
if cyc <= max(cycs1):
continue
a=sefiles.get(cyc,'mass')
cyctest=int(cyc)
#if cycle 0 or a non-existing cycle, or cycle was alrady added,
#go one cycle further
while ((type(a)==list) or (a==0)):
if cyctest<int(cycs[-1]):
print int(cycs[-1])
print 'found bad cycle',cyctest
cyctest+=1
a=sefiles.get(cyctest,'mass')
elif cyctest==int(cycs[-1]):
cyctest=int(sefiles_restart.se.cycles[-3])
a=sefiles.get(cyctest,'mass')
h=-4
while ((type(a)==list) or (a==0)):
cyctest=int(sefiles_restart.se.cycles[h])
a=sefiles.get(cyctest,'mass')
h-=1
for kk in cycs1:
if kk>=cyctest:
cycs1=cycs1[:kk]
break
else:
break
print 'Use cycle ',cyctest
if cyctest>max(cycs1):
cycs1.append(cyctest)
print len(cycs),len(cycs1)
cycs=cycs1
return cycs
def weighted_yields_pre_explosion(self,mp,sefiles,sefiles_hout,endcycle,mass_cut,isotopes):
'''
Marcos routine adapted, calc pre SN yields
explosion cycle is endcycle, at which pre SN yield is taken from
'''
import utils as u
cycle=endcycle
first_cycle=sefiles.se.cycles[0]
#mass_cut=1.60
mass_last_cycle=sefiles_hout.get(cycle,"mass")
idx_mcut=int(np.where(mass_last_cycle>mass_cut)[0][0])
# I take the max mass coord
mass_limit=mass_last_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=sefiles.get(0,isotopes[i])
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
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
def multi_DUP(self,t0_model=[],h_core_mass=False,plot_fig=True,linestyle=[],marker=[],color=[]):
'''
Plot (!) DUP parameter lambda versus mass.
There is a algorithm finding the first TP (and set t0_model) for reach run dir
which worked fine for set1. But there could be a problem with SAGB stars
so use with caution! You can use a manual input (the array t0_model) to
skip this step.
It also returns peak lum and corresponding
model number of each TP.
h_core_mass - If True: plot dependence from h free core , else star mass
t0_model - Array of first TPs of models.
examples of set1:
t0_model paramter:
z1e-2:1.65-5: [13690,3120,3163,5306]
z2e-2:1.65-5: [16033,6214,3388,5368]
dirs=["M1.65Z2.0e-02/LOGS","M2.00Z2.0e-02/LOGS","M3.00Z2.0e-02/LOGS","M5.00Z2.0e-02/LOGS"]
dirs=["M1.65Z1.0e-02/LOGS","M2.00Z1.0e-02/LOGS","M3.00Z1.0e-02/LOGS","M5.00Z1.0e-02/LOGS"]
note: M3_1e-2 shows decrease in core mass in the end
set.multi_DUP(t0_model=[13690,3120,3163,5306])
set.multi_DUP(t0_model=[16033,6214,3388,5368])
'''
if len(t0_model)==0:
t0_model=self.set_find_first_TP()
dirs=self.run_LOGS
historydata=self.run_historydata
marker_type=marker
#line_style=['--','-','-.',':']
peak_lum_model_array=[]
h1_mass_min_DUP_model_array=[]
for i in range(len(dirs)):
#color,marker_type)
peak_lum_model,h1_mass_min_DUP_model = historydata[i].find_TP_attributes(0,t0_model[i],color[i],marker_type[i],h_core_mass,no_fig=plot_fig)
peak_lum_model_array.append(peak_lum_model)
h1_mass_min_DUP_model.append(h1_mass_min_DUP_model)
return peak_lum_model_array,h1_mass_min_DUP_model_array
###the following methods allow are parts of Falks vis3.py file, thanks Falk
def plot_TPDCZ(self,fig=13,marker=[]):
'''Plots TPDCZ vs TP'''
minis=[]
for hh in range(len(self.run_historydata)):
historydata=self.run_historydata[hh]
sefiles_out=se(self.runs_H5_out[hh])
sefiles=sefiles_out
mini=float(sefiles_out.se.get('mini'))
zini=float(sefiles_out.se.get('zini'))
label=str(mini)+'$M_{\odot}$, Z='+str(zini)
minis.append(mini)
TPstart,TPmods,TP_max_env,TPend,min_m_TP,max_m_TP,DUPmods,DUPm_min_h,c13_pocket_min,c13_pocket_max=self.TPAGB_properties(historydata,sefiles_out)
###Temperature
pulse_peak_T=[]
pulse_number=[]
plt.figure(fig)
for k in range(len(TPmods)):
model=int(float(TPmods[k]))
T_profile=sefiles_out.get(model,'temperature')*sefiles_out.get('temperature_unit')
mass=sefiles_out.get(model,'mass')
#plt.figure(2)
print 'testse'
print TPmods[k]
#plt.plot(mass,T_profile,label='cycle '+str(TPmods[k]),marker=marker[k],markevery=100)
#plt.figure(11)
#models=sefiles_out.se.cycles
#models=map(int,models)
#print 'model taken: ',models[min(range(len(models)), key=lambda i: abs(models[i]-model))]
idx_mass=min(range(len(mass)), key=lambda i: abs(mass[i]-min_m_TP[k]))
pulse_peak_T.append(max(T_profile[idx_mass:]))
print k+1
print '{:.3E}'.format(T_profile[idx_mass])
pulse_number.append(k+1)
plt.plot(pulse_number,pulse_peak_T,marker=marker[hh],linestyle='--',label=label)
plt.legend()
plt.ylabel('$T_{max,PDCZ}$')
plt.xlabel('TP number')
