def prod_fac_massgrid(self,weighted=True,log=True,runs=[],isotopes=[],elements=[],cycles=[],plot_set_diagram=True,color=['r','b','g','k'],marker_type=['o','p','s','D'],line_style=['--','-','-.',':'],markersize=[6,6,6,6],line_width=[14,14,14,14],title='',withlabel=True,withlabelelement=False,label='',plot_lines=True,exp_only=False,pre_exp=False,delay=True,exp_dir='',yields_output=False): ''' Plots behaviour star mass dependent yields of different isotopes - specify dirs at the beginning Beware of different z, runs : If array is empty function uses all available directories. !! If len of elements longer than zero, than isotopes will be ignored and elements used!! ''' runs=[] HDF5_surf=[] HDF5_out=[] HDF5_restart=[] for i in range(len(self.run_dirs_name)): HDF5_surf.append(self.runs_H5_surf[i]) HDF5_out.append(self.runs_H5_out[i]) HDF5_restart.append(self.runs_H5_restart[i]) legend_k=0 sefiles=[] sefiles_hout=[] sefiles_restart=[] for i in range(len(HDF5_surf)): sefiles.append(se(HDF5_surf[i])) sefiles_hout.append(se(HDF5_out[i])) sefiles_restart.append(se(HDF5_restart[i])) isotopes11=sefiles[0].se.isotopes if len(isotopes)>0: if isotopes[0]=='allstable': isotopes=get_stable(isotopes11,get_elements=False) color=len(isotopes)*[color[0]] marker_type=len(isotopes)*[marker_type[0]] line_style=len(isotopes)*[line_style[0]] markersize=len(isotopes)*[markersize[0]] line_width=len(isotopes)*[line_width[0]] print isotopes if len(elements)>0: if elements[0]=='allstable': elements=get_stable(isotopes11,get_elements=True) color=len(elements)*[color[0]] marker_type=len(elements)*[marker_type[0]] line_style=len(elements)*[line_style[0]] markersize=len(elements)*[markersize[0]] line_width=len(elements)*[line_width[0]] print elements ####dealign with explosion pinput #test to take right dir: # if delay: exp_type='delay' else: exp_type='rapid' slist = os.listdir(exp_dir) expr = re.compile(exp_type) slist=(filter(expr.search,slist)) exp_runs_H5_restart=[] sefiles_exp=[] #to dinstiungish between pre exp and exp sources if pre_exp==False: slist = os.listdir(exp_dir) expr = re.compile(exp_type) slist=(filter(expr.search,slist)) for element in slist: run_path=exp_dir+'/'+element if not os.path.isdir(run_path): continue if os.path.isdir(run_path+"/H5_restart"): sefiles1 = os.listdir(run_path+"/H5_restart") if (filter(expr.search,sefiles1)) <1: print "Warning: No hdf5 restart files found in "+run_path+"/H5_restart" else: exp_runs_H5_restart.append(run_path+"/H5_restart") sefiles_exp.append(se(run_path+"/H5_restart")) else: exp_runs_H5_restart=HDF5_restart for k in range(len(HDF5_restart)): sefiles_exp.append(se(HDF5_restart[k])) z_index_files=[] z_values=[] j=-1 for i in range(len(HDF5_surf)): j+=1 star_z=sefiles[i].get("zini") if star_z not in z_values: z_values.append(star_z) z_index_files.append([]) z_index_files[ z_values.index(star_z)].append(i) max_yield=[] color_iso=-1 yields_1=[] t=0 if len(elements)>0: isotopes=elements legend_k=0 for w in range(len(z_index_files)): star_mass_array=[] yields=[] legend_k+=1 #iso_yields=[] #iniabu_yields_folded=[] production_factor=[] for i in range(len(isotopes)): #iso_yields.append(np.zeros(len( z_index_files[w] ))) #iniabu_yields_folded.append(np.zeros(len( z_index_files[w] ))) if exp_only==False: production_factor.append(np.zeros(len( z_index_files[w] ))) else: production_factor.append([]) ttt=0 for k in z_index_files[w]: if type(sefiles[k].get("mini")) == np.ndarray: star_mass=sefiles[k].get("mini")[0] star_z=sefiles[k].get("zini")[0] else: star_mass=sefiles[k].get("mini") star_z=sefiles[k].get("zini") print 'mini',star_mass #star_mass=sefiles[k].get("mini")[0] #star_z=sefiles[k].get("zini")[0] if cycles[k][1]==-1: print 'cycles' #print sefiles_restart[k].se.cycles endcycle=int(sefiles_restart[k].se.cycles[-2]) else: endcycle=cycles[k][1] cycs=range(cycles[k][0],endcycle,cycles[k][2]) if endcycle not in cycs: cycs=cycs+[endcycle] star_mass11=sefiles[k].get(cycs,'mass') if not len(star_mass11)==len(cycs): cycs1=[int(cycs[0])] for cyc in cycs[1:]: if cyc <= max(cycs1): continue a=sefiles[k].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=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 print cycs[-5:] if exp_only==False: star_mass_array.append(star_mass) prod_factor,isotopes_prod_fac,yields,iniabu,mass_frac_ini,remn_mass = self.weighted_yields(sefiles[k],sefiles_restart[k],isotopes,elements,cycs) #print 'yield output###################################:' #print 'wind: ',yields,prod_factor,iniabu for tt in range(len(prod_factor)): if yields_output==False: prod_factor[tt]=yields[tt]/iniabu[tt] else: prod_factor[tt]=yields[tt] else: prod_factor=[] if star_mass<=8: continue for pp in range(len(isotopes)): production_factor[pp].append([0]) star_mass_array.append(star_mass) mass=star_mass metallicity=star_z if mass>8: #in case no ppn_exp run dir is available, skip explosion contribution #this is because mass cut is then larger then actual mass of star...no mass lost for t in range(len(exp_runs_H5_restart)): if type(sefiles_exp[t].get("mini")) == np.ndarray: mass1=sefiles_exp[t].get("mini")[0] metallicity1=sefiles_exp[t].get("zini")[0] else: mass1=sefiles_exp[t].get("mini") metallicity1=sefiles_exp[t].get("zini") ###for currupted files corrupt=False if type(mass1) == list: corrupt=True if type(metallicity1) ==list: corrupt=True if ((float(mass1) < 1) or (float(mass1) >70 )) or type(mass1) == type(None): corrupt=True if ((float(metallicity1) < 0.00001) or (float(metallicity1) >1. )) or type(metallicity1) == type(None): corrupt=True if corrupt == True: metallicity1=0.02 rundir11=exp_runs_H5_restart[t].split('/')[-2] mass1=float(rundir11.split('.')[0][1:]) #print 'corruption, assinge new values',mass1,metallicity1 if mass1 == mass and metallicity1 == metallicity: #sefiles_re2=sefiles_exp[t] #calculate exp part import nugridse as mp reload(mp) #Have to find the right cycle in restart file #.se.cycles does not work when reading all files #therefore identify file with last cycle and read it: refiles=sefiles_exp[t].se.files tosort=[] refiles1=[] for i in range(len(refiles)): cycle_file=refiles[i][-18:-11] if cycle_file.isdigit(): tosort.append(int(cycle_file)) refiles1.append(refiles[i]) idx_sorted=sorted(range(len(tosort)), key=lambda k: tosort[k]) idx_file = idx_sorted[-1] if pre_exp==True: sefiles_re_cycle=mp.se(exp_runs_H5_restart[t]) else: sefiles_re_cycle=mp.se(exp_runs_H5_restart[t],refiles1[idx_file]) if len(sefiles_re_cycle.se.cycles)==0: #in the case there is not cycle in the last file...set1.2 M2- sefiles_re_cycle=mp.se(exp_runs_H5_restart[t],refiles1[idx_sorted[-2]]) if pre_exp==True: cycleend=cycs[-1] else: cycleend=sefiles_re_cycle.se.cycles[-1] #print 'identifies as ',exp_runs_H5_restart[t],'with file',refiles1[idx_file] isotopes_exp,yields_exp,iniabu_exp,mass_cut_exp,first_cycle_exp= self.weighted_yields_explosion(mp,sefiles_re_cycle,sefiles[k],isotopes=isotopes,elements=elements,cycleend=cycleend,delay=delay) if exp_only==False: #add pre-exp yields for tt in range(len(prod_factor)): #print 'star above 8M' #print isotopes_prod_fac[tt],yields[tt],yields_exp[tt] #the iniabu from the explosion is not added #because it is already in the iniabu from wind if yields_output==False: prod_factor[tt]=(yields[tt]+yields_exp[tt])/(iniabu[tt] + iniabu_exp[tt])#(iniabu[tt]) else: prod_factor[tt]=yields[tt]+yields_exp[tt] #print 'yields wind:',yields[tt] #print 'prodfac wind:',(yields[tt]/iniabu[tt]) #print 'wind+exp yields : ',yields[tt]+yields_exp[tt] else: for tt in range(len(isotopes_exp)): if yields_output==False: prod_factor.append((yields_exp[tt])/(iniabu_exp[tt])) else: prod_factor.append(yields_exp[tt]) remn_mass=mass_cut_exp #print 'exp yield',yields_exp #print 'prodfac exp:',(yields_exp[tt]/iniabu_exp[tt]) #print 'wind+exp prodfac:',prod_factor[tt] break #print 'iso exp',len(isotopes_exp) #prod_factor=np.log10(prod_factor) print len(isotopes),len(prod_factor),len(production_factor) for i in range(len(isotopes)): production_factor[i][ttt]=prod_factor[i] ttt+=1 ###plotting #if plot_set_diagram==True: #if plot_set_diagram==True: # fig_2=plt.figure(isotopes[h]) # #fig_2=plt.figure() # ax_2 = fig_2.add_subplot(1,1,1) # plt.rcParams.update({'font.size': 16}) # plt.rc('xtick', labelsize=16) # plt.rc('ytick', labelsize=16) # if log==True: # ax_2.set_yscale('log') color_iso=0 legend_k=0 for h in range(len(isotopes)): #yield_1=iso_yields[i] #mass_1=star_mass_array yield_1=[] mass_1=[] iniabu_folded=[] prod_fac_sorted=[] indices_array=sorted(range(len(star_mass_array)),key=lambda x:star_mass_array[x]) for i in indices_array: # yield_1.append(iso_yields[h][i]) mass_1.append(star_mass_array[i]) # iniabu_folded.append(iniabu_yields_folded[h][i]) prod_fac_sorted.append(production_factor[h][i]) ####Plotting prodfactor #plt.figure(fig_2.number) fig_2=plt.figure(isotopes[h]) ax_2 = fig_2.add_subplot(1,1,1) plt.rcParams.update({'font.size': 16}) plt.rc('xtick', labelsize=16) plt.rc('ytick', labelsize=16) if log==True: ax_2.set_yscale('log') ax_2.legend() ax_2.set_xlabel("M/M$_{\odot}$",fontsize=16) ax_2.minorticks_on() if yields_output==False: ax_2.set_ylabel("Production factor",fontsize=16) else: ax_2.set_ylabel("Yields/M$_{\odot}$",fontsize=16) ax_2.set_title(title) if len(label)>0: label=", "+label if len(elements)==0: plot_quantity=isotopes[h] plot_quantity="$^{"+plot_quantity.split("-")[1]+"}$"+plot_quantity.split("-")[0] else: plot_quantity=elements[h] if withlabel==True: if withlabelelement==True: plt.plot(mass_1,prod_fac_sorted,marker=marker_type[legend_k],color=color[color_iso],markersize=markersize[legend_k],mfc=color[color_iso],linewidth=line_width[legend_k],linestyle=line_style[color_iso],label=plot_quantity+" , Z="+str(star_z)+label ) else: plt.plot(mass_1,prod_fac_sorted,marker=marker_type[legend_k],color=color[color_iso],markersize=markersize[legend_k],mfc=color[color_iso],linewidth=line_width[legend_k],linestyle=line_style[color_iso],label="Z="+str(star_z)+label ) plt.legend(loc=2) print 'legend created' else: plt.plot(mass_1,prod_fac_sorted,marker=marker_type[legend_k],color=color[color_iso],markersize=markersize[legend_k],mfc=color[color_iso],linewidth=line_width[legend_k],linestyle=line_style[color_iso]) m_max=max(star_mass_array)+2 if plot_lines==True: plt.plot([0,m_max],[1,1],"k--",linewidth=3) plt.plot([0,m_max],[2,2],"k--",linewidth=1) plt.plot([0,m_max],[0.5,0.5],"k--",linewidth=1) color_iso+=1 legend_k+=1 ##### x_imf=[] y_imf=[] m_max=max(star_mass_array)+2. #if plot_set_diagram==True: # plt.figure(fig_2.number) ## ###plot figure 1 # plt.legend() ## plt.xlabel("M/M$_{\odot}$",fontsize=20) # plt.minorticks_on() # plt.ylabel("Overproduction factor",fontsize=20) # plt.title(title)