162,164c162
< 				if len(multi_dir)>0:
< 					element_1=element.split('/')[-1]	
< 						
---
> 								
169,170c167
< 					element_1=element
< 				run_dirs_name.append(element_1)
---
> 				run_dirs_name.append(element)
193d189
< 
280,282c276,278
< 			for j in range(len(mass_dummy)):
<         			if h_dummy[j] > 0.05:
<                 			bottom_of_envelope = mass_dummy[j]
---
> 			for i in range(len(mass_dummy)):
>         			if h_dummy[i] > 0.05:
>                 			bottom_of_envelope = mass_dummy[i]
290,362d285
<         def final_bottom_envelope(self):
<                 '''
<                         For paper1 extension:
<                          Numbers of remnant mass shell masses, exists also in mesa_set!
<                 '''
<                 inim=[]
<                 remnm=[]
< 		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')
< 			#print "use run ",self.extra_label
< 			for j in range(len(mass_dummy)):
<                                 if h_dummy[j] > 1e-1:
<                                         bottom_of_envelope = mass_dummy[j]
<                                         break
< 			#check core
< 			#if c_dummy[0]>o_dummy[0]:
< 			#	##c/o
< 			#	c_core_center.append("T") 
< 		#		o_mg_ne_core=False
< 			#	for j in range(len(mass_dummy)):
< 			#		##find onemg core
< 			#		if c_dummy[j]<o_dummy[j] and mass_dummy[j]<0.5:
< 			#			o_mg_ne_core=True
< 			#			small_co_core.append(mass_dummy[j])
< 			#			break
< 			#		#then
< 			#	for j in range(len(mass_dummy)):
< 			#		#in any case ccore exists with this boundary
< 			#		if c_dummy[i]<1e-1:
< 			#			c_core.append(mass_dummy[j])
< 			#			if o_mg_ne_core[-1] ==False:
< 			#				o_core.append(0)
< 			#				small_co_core.append(0)
< 			#			break
< 			#	#in case of o core..
< 			#	for j in range(len(mass_dummy)):
< 			#		if (c_dummy[j]>o_dummy[j] and mass_dummy[j]>0.5) and o_mg_ne_core[-1] ==True:
< 			#			o_core.append(mass_dummy[j])
< 			#else:
<                         #        o_core.append(0)
< 			#	c_core_center.append("F")
<                         #        o_mg_ne_core=True
< 			#	small_co_core.append(0)
< 			#	for j in range(len(mass_dummy)):
<                         #                ##find onemg core
<                         #                if c_dummy[j]>o_dummy[j]:
< 			#			o_core[-1]=mass_dummy[j]
< 			#			break
< 		#		for j in range(len(mass_dummy)):
< 			#		if c_dummy[j]<1e-1:
<                          #               	c_core.append(mass_dummy[j])
< 					
< 
<                         inim.append(m1p65_last.get("mini")[0])
<                         remnm.append(bottom_of_envelope)
<                 print "M_initial | M_remn/bottom of envelope"
<                 for i in range(len(inim)):
<                         print inim[i],"|",remnm[i]
< 
< 
< 
< 
< 
< 
469d391
< 				print "mini get done"
479c401
< 				prod_factor,isotopes_prod_fac,yields,iniabu,remn_mass =self.weighted_yields(sefiles[k],sefiles_hout[k],isotopes,elements,cycles[k][0],endcycle,cycles[k][2])
---
> 				prod_factor,isotopes_prod_fac,yields,iniabu =self.weighted_yields(sefiles[k],sefiles_hout[k],isotopes,elements,cycles[k][0],endcycle,cycles[k][2])
549,554c471
< 				if len(elements)==0:
< 					plot_quantity=isotopes[h]
< 					plot_quantity="$^{"+plot_quantity.split("-")[1]+"}$"+plot_quantity.split("-")[0]
< 				else:
< 					plot_quantity=elements[h]
<                                 plt.plot(mass_1,prod_fac_sorted,marker=marker_type[legend_k],color=color[color_iso],markersize=10,mfc=color[color_iso],linewidth=line_width[legend_k],linestyle=line_style[color_iso],label=plot_quantity+" , Z="+str(star_z)+label  )				
---
>                                 plt.plot(mass_1,prod_fac_sorted,marker=marker_type[legend_k],color=color[color_iso],markersize=10,mfc=color[color_iso],linewidth=line_width[legend_k],linestyle=line_style[color_iso],label=isotopes[h]+" , Z="+str(star_z)+label  )				
561c478
< 					plt.ylabel("Log. production factor",fontsize=20)
---
> 					plt.ylabel("Production factor",fontsize=20)
747c664
< 					prod_factor,isotopes_prod_fac,yields,iniabu,remn_mass =self.weighted_yields(sefiles[k],sefiles_hout[k],isotopes,startcycle,endcycle,cycles[k][2])
---
> 					prod_factor,isotopes_prod_fac,yields,iniabu =self.weighted_yields(sefiles[k],sefiles_hout[k],isotopes,startcycle,endcycle,cycles[k][2])
1092c1009
< 					prod_factor,isotopes_prod_fac,yields,iniabu,remn_mass =self.weighted_yields(sefiles[k],sefiles_hout[k],isotopes,startcycle,endcycle,cycles[k][2])
---
> 					prod_factor,isotopes_prod_fac,yields,iniabu =self.weighted_yields(sefiles[k],sefiles_hout[k],isotopes,startcycle,endcycle,cycles[k][2])
1619c1536
< 	def set_triple_isotope_plot(self,runs=[],xiso=['Zr',96,'Zr',94],yiso=['Zr',92,'Zr',94],graintype=['sic'],C_star_only=True,dcycle=100,USEEPP_path='USEEPP',extra_label=[],sens_runs=True):
---
> 	def set_triple_isotope_plot(self,runs=[],xiso=['Zr',96,'Zr',94],yiso=['Zr',92,'Zr',94],graintype=['sic'],C_star_only=True,dcycle=100,USEEPP_path='USEEPP',sens_runs=True):
1633a1551
>                 extra_label=[]
1637,1638c1555
< 			if len(extra_label)==0:
<                         	extra_label=self.extra_label
---
>                         extra_label=self.extra_label
1657c1574
<                         	legend=str(mass)+"M$_{\odot}$, "+str(z)+"Z$_{\odot}$, "+extra_label[i]
---
>                         	legend=str(mass)+"M$_{\odot}$, "+str(z)+"Z$_{\odot}$" #+", "+extra_label[i]
1671,2658d1587
< 	def write_GCE_input_fxt2(self,input_file='maltov1.new4',input_file2='isotopic_table_set1.1_yields_exp.txt',output_file='maltov1.new5',exp_type='delay',label_masses='set1.2'):
< 	
< 	
<                 import re
<                 f1=open(input_file)
<                 lines=f1.readlines()
<                 f1.close()
<                 isotopes=[]
< 		isotopes_out=[]
<                 iso_line_idx=[]
<                 i=0
< 		first_change=True
< 
< 		##########read input file##############
<                 #GCE isotopes
<                 for line in lines:
<                         if line[0] != ' ' and line[0:3] != 'rem':
<                                 if line[0].isalpha():
<                                         a= line[:6].split()[0]
< 					isotopes_out.append(a)
<                                         match = re.match(r"([a-z]+)([0-9]+)",a, re.I)
<                                         items = match.groups()
<                                         isotopes.append(items[0].capitalize()+'-'+items[1])
<                                         iso_line_idx.append(i)
<                         if line[0:3] == 'rem':
< 				isotopes_out.append('rem')
<                                 isotopes.append('rem')
<                                 iso_line_idx.append(i)
<                         i+=1
< 
< 		bb_massfrac=[0]*len(isotopes)
< 		type1a_yields=[0]*len(isotopes)
< 		nova_yields=[0]*len(isotopes)
< 		metallicities=[0]*len(isotopes)
< 		masses=[0]*len(isotopes)
< 		yields=[0]*len(isotopes)
< 		label=[0]*len(isotopes)
< 
< 		i=-1
< 		j=-1
< 		while (True):
< 			i+=1
<                         if i>len(lines)-1:
<                                 break
< 			line=lines[i]
< 			if i in iso_line_idx:
< 				j+=1
< 			#read data
<                         if 'big bang' in line:
< 				
<                                 bb_massfrac[j]=float(line[:13])
<                                 continue
<                         #type 1a input ignore
<                         if 'number of type1a' in line:
<                                 number=int(line[0])
< 				type1a_yields[j]=[0]*number
< 				k=0
< 				for sn in range(i+1,i+number+1):
< 					type1a_yields[j][k]=(float(lines[sn][:13]))				
< 					k+=1
< 				i=i+number	
< 				continue
<                         #nova input ignore
<                         if 'number of nova' in line:
<                                 number=int(line[0])
< 				k=0
< 				nova_yields[j]=[0]*number
<                                 for nova in range(i+1,i+number+1):
<                                          nova_yields[j][k]=(float(lines[nova][:13])) 
<                                 	 k+=1
< 				i=i+number   
<                                 continue
< 			if 'number of metallicity values' in line:
< 				number_metallicities=int(line[0])
< 				metallicities[j]=[0]*(number_metallicities)
< 				masses[j]=[0]*(number_metallicities)
< 				yields[j]=[0]*(number_metallicities)
< 				label[j]=[0]*(number_metallicities)	
< 				z_counter=-1
< 				continue
< 			if '* metallicity' in line:
< 				z_counter+=1
< 				metallicities[j][z_counter]=float(line[:13])				
< 				continue
< 			if '* number of masses @ this metallicity' in line:
< 				number_masses=int(line[:2])
< 				masses[j][z_counter]=[0]*(number_masses)
< 				yields[j][z_counter]=[0]*(number_masses)
< 				label[j][z_counter]=[0]*(number_masses)
< 				k=0
< 				for mass in range(i+1,i+number_masses+1):
< 					masses[j][z_counter][k]=float(lines[mass][:13])
< 					yields[j][z_counter][k]=float(lines[mass][14:24])
< 					#check for label
< 					if len(lines[mass])>27:
< 						label1=lines[mass][27:].strip()
< 									
< 					label[j][z_counter][k] = label1
< 					k+=1
< 				i=i+number_masses
<                                 continue
<                 #print bb_massfrac
<                 #print type1a_yields
<                 #print nova_yields
<                 #print metallicities
<                 #print masses
<                 #print yields
< 		##print label	
< 		
< 		##############read new simulation input###########
< 
< 		f1=open(input_file2)
< 		lines_txt=f1.readlines()
< 		f1.close()
< 		isotopes_1=[]
< 		for i in range(len(lines_txt)):
< 			if lines_txt[i][0]=='H':
< 				continue
< 			if 'specie' in lines_txt[i]:
< 				continue
< 			name1=lines_txt[i].split('&')[1].strip()
< 			name1= name1.replace(' ','')
< 			match = re.match(r"([a-z]+)([0-9]+)",name1, re.I)
< 			items = match.groups()
< 			isotopes_1.append(items[0].capitalize()+'-'+items[1])
< 
< 		#print isotopes_1
< 
<                 #check if GCE isotopes are all available in simulation input
<                 #sefiles=se(self.runs_H5_surf[0])
<                 isotopes_2=[]
<                 no_iso=[]
<                 available_iso=[]
<                 for j in range(len(isotopes_1)):
<                         if is_stable(isotopes_1[j])=='t':
<                                 isotopes_2.append(isotopes_1[j])
<                 for j in range(len(isotopes)):
<                         #if isotope is not available stop       
<                         if isotopes[j] not in isotopes_2:
<                                 print 'Isotope ',isotopes[j],'is not available'
<                                 no_iso.append(j)
<                                 #print 'STOP'
<                         else:
<                                 available_iso.append(isotopes[j])
<                 print 'All isotopes for GCE available'
< 		
< 		label_done=False
< 		for i in range(len(lines_txt)):
< 			#header
< 			if 'specie' in lines_txt[i]:
< 				lins=lines_txt[i].split('&')[2:-1]
< 				print lins
< 				label1=[]
< 				idx_runs=[]
< 				for k in range(len(lins)):
< 					if exp_type in lins[k]:
< 						label1.append(lins[k])
< 						idx_runs.append(k)
< 				#print idx_runs
< 				label_done=True
< 				cols=[0]*len(label1)
< 				for k in range(len(cols)):
< 					cols[k]=[]
< 				continue
< 			if label_done==False:
< 				continue
< 
< 			#if 'hline' in lines_txt[i] and i<15:
< 			#	continue
< 			#if 'hline' in lines_txt[i]:
< 			#	break
< 
< 			#print lines_txt
< 			lins=lines_txt[i].split('&')[2:-1]
< 			j=0
< 			#print cols[j]
< 			for k in idx_runs:
< 				cols[j].append(float(lins[k].strip('\n\t').replace("\\", "")))
<                 		j+=1
< 		#print cols
< 		
< 
<                 #labelout,cols,remn_masses_1 = self.write_prod_fact_stellarwinds(cycles=cycles,isotopes=available_iso,ascii_1=False,latex=False,GCE_input_fxt=True)
< 
< 
< 		#from paper1
<                 if label_masses == 'set1.1':
<                         remn_masses=[12.349,14.120,14.232]
< 			metallicity_new=0.01
<                 if label_masses == 'set1.2':
<                         metallicity_new=0.02
< 			remn_masses=[12.132,13.974,13.738,12.495,13.428]
<                 if label_masses == 'set1.3a':
<                         metallicity_new=0.006
< 			remn_masses=[12.349,14.120,14.232]
<                 if label_masses== 'set1.5a':
<                         metallicity_new=0.0001
< 			remn_masses=[12.349,14.120,14.232]
< 
<                 #add values for missing isotopes, value=99
< 		print 'testteest'
< 		print no_iso
<                 for i in no_iso:
< 			print no_iso
<                         for j in range(len(cols)):
< 				if i == no_iso[-1]:
< 				# here come the right remnant masses for massive stars
< 				#not just pre-collapse masses
< 					print 'adding normal value for rem'
< 					cols[j].insert(i,remn_masses[j])
< 				else:
<                                 	cols[j].insert(i,0)
< 
<                 #metallicity1='{:.4E}'.format(float(metallicity))
<                 masses_new1=[]
< 		labelout=label1
<                 for i in range(len(labelout)):
<                         mass=labelout[i].split('Msun')[0]
<                         #mass='{:.4E}'.format(float(mass))
<                         masses_new1.append(float(mass))
<                 #print 'Available input masses',masses
< 
< 	
< 		yields_new=[0]*len(cols[0])
< 		masses_new=[0]*len(cols[0])
< 		label_new=[0]*len(cols[0])
< 		
< 		#for every mass
< 		for k in range(len(yields_new)):	
< 			yields_new[k]=[]
< 			masses_new[k]=masses_new1
< 			label_new[k]=[label_masses]*len(masses_new1)
< 			for w in range(len(cols)):
< 				yields_new[k].append(cols[w][k])
< 
< 		#print metallicity_new
<                 #print yields_new
< 		#print label_new
< 		#print masses   
< 
< 		##### Merge data###############################
< 
< 
< 		
< 
< 
< 		metallicities_out=len(isotopes)*[0]
< 		for t in range(len(isotopes)):
< 			pos_idx=0
< 			#assuming that if this value is still in it, fresh file from fxt without any new z
< 			#print metallicities,float('1.9000E-04')	
< 			if float('1.9000E-04') in metallicities[0]:
< 				metallicities_out[t]=metallicities[t][:2]+[metallicity_new]
< 				pos_idx=2
< 			else:
< 				metal_sum=[]
< 				st=0	
< 				for metal in metallicities[t]:
< 					#metallicities_out[t]=metallicities[t]
< 					if (metal>metallicity_new) and (metallicity_new not in metal_sum):
< 						metal_sum.append(metallicity_new)
< 						metal_sum.append(metallicities[t][st])
< 					else:
< 						metal_sum.append(metallicities[t][st])
< 					st+=1
< 				#in the case the added Z is larger than all available Z
< 				if (metallicity_new not in metal_sum):
< 					if len(metal_sum) == len(metallicities[t]):
< 						metal_sum.append(metallicity_new)
< 				metallicities_out[t]=metal_sum
< 				pos_idx=metallicities_out[t].index(metallicity_new)		
< 		
< 		masses_out=len(isotopes)*[0]
< 		yields_out=len(isotopes)*[0]
< 		label_out=len(isotopes)*[0]
< 		print metallicities_out[0]
< 		for t in range(len(isotopes)):
< 			masses_out[t]=len(metallicities_out)*[0]
< 			yields_out[t]=len(metallicities_out)*[0]
< 			label_out[t]=len(metallicities_out)*[0]
< 			w=0	
< 			for wt in range(len(metallicities_out[0])):
< 				#if case that it is a fresh file 
< 				if float('1.9000E-04') in metallicities[0] and wt <2:
< 					masses_out[t][wt]=masses[t][w]
< 					yields_out[t][wt]=yields[t][w]
< 					label_out[t][wt]=label[t][w]
< 					w+=1				
< 				else:
< 					#if we are a t the position to add new metallicity
< 					#or in primary case of this function to add yield values of exp
< 					#print isotopes[t]
< 					if wt==pos_idx:
< 						#print 'add new z or primary add yields to old Z'
< 					#	print t
< 						#number of to added masses different from already existing data
< 						yield_sum=[]
< 						for idx1 in range(len(masses[t][w])):
< 							if masses[t][w][idx1] in masses_new[t]:
< 								idx2=masses_new[t].index(masses[t][w][idx1])
< 								if not isotopes[t]=='rem':
< 									yield_sum.append( float(yields[t][w][idx1] + yields_new[t][idx2] ))
< 								else:
< 									yield_sum.append( float(yields_new[t][idx2] ))
< 							else:
< 									yield_sum.append(float(yields[t][w][idx1]))
< 						masses_out[t][wt]=masses[t][w]
< 						#print yield_sum
< 						yields_out[t][wt]=yield_sum
< 						label_out[t][wt]=label[t][w]
< 						#print wt,w
< 					else:
< 						#print metallicities_out
< 						#print 'else',wt,w
< 						#print masses_out[t],masses[t]
< 						masses_out[t][wt]=masses[t][w]
< 						yields_out[t][wt]=yields[t][w]
< 						label_out[t][wt]=label[t][w]
< 					w+=1
< 
< 
< 		###### Prepare output in right format######################
< 		#isotopes_out defined already above
< 		bb_massfrac_out=bb_massfrac
<                 type1a_yields_out=type1a_yields
<                 nova_yields_out=nova_yields
<                 #metallicities_out=metallicities
<                 #masses_out=masses
<                 #yields_out=yields
<                 #label_out=label 		
< 		newline=''
< 		for i in range(len(isotopes)):
< 			newline+= (isotopes_out[i].ljust(29)+'* symbol name'+'\n') 	
< 			newline+= ('  '+'{:.4E}'.format(bb_massfrac_out[i]).ljust(27)+'* big bang mass fraction'+'\n')
< 			newline+= (str(len(type1a_yields_out[i])).ljust(29)+'* number of type1a'+'\n')
< 			ia_label=['w7 tny86','sub chandra 1','sub chandra 2','sub chandra 7','sub chandra 8','nse 1a']
< 			for k in range(len(type1a_yields_out[i])):
< 				newline+= ('  '+'{:.4E}'.format(type1a_yields_out[i][k]).ljust(29)+ia_label[k]+'\n')
< 			newline+= (str(len(nova_yields_out[i])).ljust(29)+'* number of nova'+'\n')
< 			nova_label='* nova yield 1.0 1.25 & 1.35 msun'
< 			for k in range(len(nova_yields_out[i])):
< 				if k==0:
< 					newline+= ('  '+'{:.4E}'.format(nova_yields_out[i][k]).ljust(27)+nova_label+'\n')
< 				else:
< 					newline+= ('  '+'{:.4E}'.format(nova_yields_out[i][k]).ljust(27)+'\n')
< 			newline+= (str(len(metallicities_out[i])).ljust(29)+'* number of metallicity values'+'\n')
< 			for k in range(len(metallicities_out[i])):
< 				newline+= ('  '+'{:.4E}'.format(metallicities_out[i][k]).ljust(27)+'* metallicity'+'\n') 
< 				newline+= (str(len(masses_out[i][k])).ljust(29)+'* number of masses @ this metallicity'+'\n')
< 				for w in range(len(masses_out[i][k])):
< 					#the following is needed for the right label setting
< 					if w==0:
< 						newline+= ('  '+'{:.4E}'.format(masses_out[i][k][w])+'  '+'{:.4E}'.format(yields_out[i][k][w])+(7*' ')+label_out[i][k][w]+'\n')
< 						label_old=label_out[i][k][w]	
< 					elif label_out[i][k][w] != label_old:
< 						newline+= ('  '+'{:.4E}'.format(masses_out[i][k][w])+'  '+'{:.4E}'.format(yields_out[i][k][w])+(7*' ')+label_out[i][k][w]+'\n')	
< 						label_old=label_out[i][k][w]
< 					else:
< 						newline+= ('  '+'{:.4E}'.format(masses_out[i][k][w])+'  '+'{:.4E}'.format(yields_out[i][k][w])+'\n')
< 			
< 
< 
< 		###write out data
< 		
< 		#to fix bug where fxt GCe cannot deal with 99 yields
< 		newline=newline.replace('9.9000E+01','0.0000E+00')
< 		f2=open(output_file,'w')
< 		f2.write(newline)
< 		f2.close()
< 
< 		
< 
< 	def write_GCE_input_fxt1(self,input_file='maltov1.orig',output_file='maltov1.new',first_change=True,label_masses='set1.2',cycles=20*[[0,-1,500]]):
< 	
< 	
<                 import re
<                 f1=open(input_file)
<                 lines=f1.readlines()
<                 f1.close()
<                 isotopes=[]
< 		isotopes_out=[]
<                 iso_line_idx=[]
<                 i=0
< 
< 
< 		##########read input file##############
<                 #GCE isotopes
<                 for line in lines:
<                         if line[0] != ' ' and line[0:3] != 'rem':
<                                 if line[0].isalpha():
<                                         a= line[:6].split()[0]
< 					isotopes_out.append(a)
<                                         match = re.match(r"([a-z]+)([0-9]+)",a, re.I)
<                                         items = match.groups()
<                                         isotopes.append(items[0].capitalize()+'-'+items[1])
<                                         iso_line_idx.append(i)
<                         if line[0:3] == 'rem':
< 				isotopes_out.append('rem')
<                                 isotopes.append('rem')
<                                 iso_line_idx.append(i)
<                         i+=1
< 
< 		bb_massfrac=[0]*len(isotopes)
< 		type1a_yields=[0]*len(isotopes)
< 		nova_yields=[0]*len(isotopes)
< 		metallicities=[0]*len(isotopes)
< 		masses=[0]*len(isotopes)
< 		yields=[0]*len(isotopes)
< 		label=[0]*len(isotopes)
< 
< 		i=-1
< 		j=-1
< 		while (True):
< 			i+=1
<                         if i>len(lines)-1:
<                                 break
< 			line=lines[i]
< 			if i in iso_line_idx:
< 				j+=1
< 			#read data
<                         if 'big bang' in line:
< 				
<                                 bb_massfrac[j]=float(line[:13])
<                                 continue
<                         #type 1a input ignore
<                         if 'number of type1a' in line:
<                                 number=int(line[0])
< 				type1a_yields[j]=[0]*number
< 				k=0
< 				for sn in range(i+1,i+number+1):
< 					type1a_yields[j][k]=(float(lines[sn][:13]))				
< 					k+=1
< 				i=i+number	
< 				continue
<                         #nova input ignore
<                         if 'number of nova' in line:
<                                 number=int(line[0])
< 				k=0
< 				nova_yields[j]=[0]*number
<                                 for nova in range(i+1,i+number+1):
<                                          nova_yields[j][k]=(float(lines[nova][:13])) 
<                                 	 k+=1
< 				i=i+number   
<                                 continue
< 			if 'number of metallicity values' in line:
< 				number_metallicities=int(line[0])
< 				metallicities[j]=[0]*(number_metallicities)
< 				masses[j]=[0]*(number_metallicities)
< 				yields[j]=[0]*(number_metallicities)
< 				label[j]=[0]*(number_metallicities)	
< 				z_counter=-1
< 				continue
< 			if '* metallicity' in line:
< 				z_counter+=1
< 				metallicities[j][z_counter]=float(line[:13])				
< 				continue
< 			if '* number of masses @ this metallicity' in line:
< 				number_masses=int(line[:2])
< 				masses[j][z_counter]=[0]*(number_masses)
< 				yields[j][z_counter]=[0]*(number_masses)
< 				label[j][z_counter]=[0]*(number_masses)
< 				k=0
< 				for mass in range(i+1,i+number_masses+1):
< 					masses[j][z_counter][k]=float(lines[mass][:13])
< 					yields[j][z_counter][k]=float(lines[mass][14:24])
< 					#check for label
< 					if len(lines[mass])>27:
< 						label1=lines[mass][27:].strip()
< 									
< 					label[j][z_counter][k] = label1
< 					k+=1
< 				i=i+number_masses
<                                 continue
<                 #print bb_massfrac
<                 #print type1a_yields
<                 #print nova_yields
<                 #print metallicities
<                 #print masses
<                 #print yields
< 		#print label	
< 		
< 		##############read new simulation input###########
< 
<                 #check if GCE isotopes are all available in simulation input
<                 sefiles=se(self.runs_H5_surf[0])
<                 isotopes_1=sefiles.se.isotopes
<                 isotopes_2=[]
<                 no_iso=[]
<                 available_iso=[]
<                 for j in range(len(isotopes_1)):
<                         if is_stable(isotopes_1[j])=='t':
<                                 isotopes_2.append(isotopes_1[j])
<                 for j in range(len(isotopes)):
<                         #if isotope is not available stop       
<                         if isotopes[j] not in isotopes_2:
<                                 print 'Isotope ',isotopes[j],'is not available'
<                                 no_iso.append(j)
<                                 #print 'STOP'
<                         else:
<                                 available_iso.append(isotopes[j])
<                 print 'All isotopes for GCE available'
<                 #cycles=20*[[0,-1,500]]
<                 labelout,cols,remn_masses_1 = self.write_prod_fact_stellarwinds(cycles=cycles,isotopes=available_iso,ascii_1=False,latex=False,GCE_input_fxt=True)
< 
<                 #to include remnants
<                 remn_masses=[]
<                 #remn_masses_1=20*[3]
<                 for i in range(len(remn_masses_1)):
<                         #remnant_1='{:.4E}'.format(float(remn_masses_1[i]))
<                         remn_masses.append(remn_masses_1[i])
< 
<                 #add values for missing isotopes, value=99
< 		print 'testteest'
< 		print no_iso
<                 for i in no_iso:
< 			print no_iso
<                         for j in range(len(cols)):
< 				if i == no_iso[-1]:
< 				# to add the remnant masses
< 					print 'adding normal value for rem'
< 					cols[j].insert(i,remn_masses[j])
< 				else:
<                                 	cols[j].insert(i,0)
< 
<                 #because the set runs have the same Z:
<                 #but first item in label is 'specie', so skip this
< 		#trick here, measure z only from first input file
< 		#massive stars are therefore not checked and treated as other z
<                 metallicity_new=float(labelout[1].split('=')[1])
<                 #metallicity1='{:.4E}'.format(float(metallicity))
<                 masses_new1=[]
<                 for i in range(1,len(labelout)):
<                         mass=labelout[i].split('Msun')[0]
<                         #mass='{:.4E}'.format(float(mass))
<                         masses_new1.append(float(mass))
<                 #print 'Available input masses',masses
< 
< 	
< 		yields_new=[0]*len(cols[0])
< 		masses_new=[0]*len(cols[0])
< 		label_new=[0]*len(cols[0])
< 		
< 		#for every mass
< 		for k in range(len(yields_new)):	
< 			yields_new[k]=[]
< 			masses_new[k]=masses_new1
< 			label_new[k]=[label_masses]*len(masses_new1)
< 			for w in range(len(cols)):
< 				yields_new[k].append(cols[w][k])
< 
< 		#print metallicity_new
<                 #print yields_new
< 		#print label_new
< 		#print masses   
< 
< 		##### Merge data###############################
< 		metallicities_out=len(isotopes)*[0]
< 		for t in range(len(isotopes)):
< 			pos_idx=0
< 			#assuming that if this value is still in it, fresh file from fxt without any new z
< 			#print metallicities,float('1.9000E-04')	
< 			if float('1.9000E-04') in metallicities[0]:
< 				metallicities_out[t]=metallicities[t][:2]+[metallicity_new]
< 				pos_idx=2
< 			else:
< 				metal_sum=[]
< 				st=0	
< 				for metal in metallicities[t]:
< 					#metallicities_out[t]=metallicities[t]
< 					if (metal>metallicity_new) and (metallicity_new not in metal_sum):
< 						metal_sum.append(metallicity_new)
< 						metal_sum.append(metallicities[t][st])
< 					else:
< 						metal_sum.append(metallicities[t][st])
< 					st+=1
< 				#in the case the added Z is larger than all available Z
< 				if len(metal_sum) == len(metallicities[t]):
< 					metal_sum.append(metallicity_new)
< 				metallicities_out[t]=metal_sum
< 				pos_idx=metallicities_out[t].index(metallicity_new)		
< 		
< 		masses_out=len(isotopes)*[0]
< 		yields_out=len(isotopes)*[0]
< 		label_out=len(isotopes)*[0]
< 		print metallicities_out[0]
< 		for t in range(len(isotopes)):
< 			masses_out[t]=len(metallicities_out)*[0]
< 			yields_out[t]=len(metallicities_out)*[0]
< 			label_out[t]=len(metallicities_out)*[0]
< 			w=0	
< 			for wt in range(len(metallicities_out[0])):
< 				#if case that it is a fresh file 
< 				if float('1.9000E-04') in metallicities[0] and wt <2:
< 					masses_out[t][wt]=masses[t][w]
< 					yields_out[t][wt]=yields[t][w]
< 					label_out[t][wt]=label[t][w]
< 					w+=1				
< 				else:
< 					#if we are a t the position to add new metallicity
< 					if wt==pos_idx:
< 						print 'add new z'
< 						masses_out[t][wt]=masses_new[t]
< 						yields_out[t][wt]=yields_new[t]
< 						label_out[t][wt]=label_new[t]
< 						print wt,w
< 					else:
< 						print metallicities_out
< 						print 'else',wt,w
< 						print masses_out[t],masses[t]
< 						masses_out[t][wt]=masses[t][w]
< 						yields_out[t][wt]=yields[t][w]
< 						label_out[t][wt]=label[t][w]
< 						w+=1
< 
< 
< 		###### Prepare output in right format######################
< 		#isotopes_out defined already above
< 		bb_massfrac_out=bb_massfrac
<                 type1a_yields_out=type1a_yields
<                 nova_yields_out=nova_yields
<                 #metallicities_out=metallicities
<                 #masses_out=masses
<                 #yields_out=yields
<                 #label_out=label 		
< 		newline=''
< 		for i in range(len(isotopes)):
< 			newline+= (isotopes_out[i].ljust(29)+'* symbol name'+'\n') 	
< 			newline+= ('  '+'{:.4E}'.format(bb_massfrac_out[i]).ljust(27)+'* big bang mass fraction'+'\n')
< 			newline+= (str(len(type1a_yields_out[i])).ljust(29)+'* number of type1a'+'\n')
< 			ia_label=['w7 tny86','sub chandra 1','sub chandra 2','sub chandra 7','sub chandra 8','nse 1a']
< 			for k in range(len(type1a_yields_out[i])):
< 				newline+= ('  '+'{:.4E}'.format(type1a_yields_out[i][k]).ljust(29)+ia_label[k]+'\n')
< 			newline+= (str(len(nova_yields_out[i])).ljust(29)+'* number of nova'+'\n')
< 			nova_label='* nova yield 1.0 1.25 & 1.35 msun'
< 			for k in range(len(nova_yields_out[i])):
< 				if k==0:
< 					newline+= ('  '+'{:.4E}'.format(nova_yields_out[i][k]).ljust(27)+nova_label+'\n')
< 				else:
< 					newline+= ('  '+'{:.4E}'.format(nova_yields_out[i][k]).ljust(27)+'\n')
< 			newline+= (str(len(metallicities_out[i])).ljust(29)+'* number of metallicity values'+'\n')
< 			for k in range(len(metallicities_out[i])):
< 				newline+= ('  '+'{:.4E}'.format(metallicities_out[i][k]).ljust(27)+'* metallicity'+'\n') 
< 				newline+= (str(len(masses_out[i][k])).ljust(29)+'* number of masses @ this metallicity'+'\n')
< 				for w in range(len(masses_out[i][k])):
< 					#the following is needed for the right label setting
< 					if w==0:
< 						newline+= ('  '+'{:.4E}'.format(masses_out[i][k][w])+'  '+'{:.4E}'.format(yields_out[i][k][w])+(7*' ')+label_out[i][k][w]+'\n')
< 						label_old=label_out[i][k][w]	
< 					elif label_out[i][k][w] != label_old:
< 						newline+= ('  '+'{:.4E}'.format(masses_out[i][k][w])+'  '+'{:.4E}'.format(yields_out[i][k][w])+(7*' ')+label_out[i][k][w]+'\n')	
< 						label_old=label_out[i][k][w]
< 					else:
< 						newline+= ('  '+'{:.4E}'.format(masses_out[i][k][w])+'  '+'{:.4E}'.format(yields_out[i][k][w])+'\n')
< 			
< 
< 
< 		###write out data
< 
< 		f2=open(output_file,'w')
< 		f2.write(newline)
< 		f2.close()
< 
< 
< 
< 
<         def write_GCE_input_fxt(self,input_file='maltov1.orig',output_file='maltov1.new',first_change=True,label_masses='set1.2',cycles=20*[[0,-1,500]]):
< 
<                 '''
<                         Reads input_file of input for GCE code of fxt and replaces
<                         data with available data from the simulations. Also creates
<                         a file *output_file*_differences to keep track of the replacements.
<                 '''
< 
< 		import re
< 		f1=open(input_file)
< 		lines=f1.readlines()
< 		f1.close()
< 		isotopes=[]
< 		iso_line_idx=[]
< 		i=0
< 		#GCE isotopes
< 		for line in lines:
< 			if line[0] != ' ' and line[0:3] != 'rem':
< 				if line[0].isalpha():
< 					a= line[:6].split()[0]
< 					match = re.match(r"([a-z]+)([0-9]+)",a, re.I)
< 					items = match.groups()
< 					isotopes.append(items[0].capitalize()+'-'+items[1])
< 					iso_line_idx.append(i)
< 			if line[0:3] == 'rem':
<                 		isotopes.append('rem')
<                 		iso_line_idx.append(i)
< 			i+=1
< 		
< 		print isotopes
< 		print iso_line_idx
< 		#isotopes=['H-1', 'H-2', 'He-3', 'He-4', 'Li-6', 'Li-7']
< 
< 		#check if GCE isotopes are all available in simulation input
< 		sefiles=se(self.runs_H5_surf[0])
< 		isotopes_1=sefiles.se.isotopes
< 		isotopes_2=[]
< 		no_iso=[]
< 		available_iso=[]
< 		for j in range(len(isotopes_1)):
< 			if is_stable(isotopes_1[j])=='t':
< 				isotopes_2.append(isotopes_1[j])
< 		for j in range(len(isotopes)):
< 			#if isotope is not available stop	
< 			if isotopes[j] not in isotopes_2:
< 				print 'Isotope ',isotopes[j],'is not available'
< 				no_iso.append(j)
< 				#print 'STOP'
< 			else:
< 				available_iso.append(isotopes[j])
< 		print 'All isotopes for GCE available'
< 		#cycles=20*[[0,-1,500]]
< 		label,cols,remn_masses_1 = self.write_prod_fact_stellarwinds(cycles=cycles,isotopes=available_iso,ascii_1=False,latex=False,GCE_input_fxt=True)
< 		
< 
< 		#to include remnants
< 		remn_masses=[]
< 		#remn_masses_1=20*[3]
< 		for i in range(len(remn_masses_1)):
< 		        remnant_1='{:.4E}'.format(float(remn_masses_1[i]))
<         		remn_masses.append(remnant_1)
< 
< 
< 		
< 		#add values for missing isotopes, value=99
< 		for i in no_iso:
< 			for j in range(len(cols)):
< 				cols[j].insert(i,0)		
< 		print isotopes
< 		print cols
< 		#because the set runs have the same Z:
< 		#but first item in label is 'specie', so skip this
< 		metallicity=label[-1].split('=')[1]
< 		metallicity='{:.4E}'.format(float(metallicity))
< 		masses=[]
< 		for i in range(1,len(label)):
< 			mass=label[i].split('Msun')[0]
< 			mass='{:.4E}'.format(float(mass))
< 			masses.append(mass)
< 		print 'Available input masses',masses
< 
< 		#Write them into output_file
< 
< 		print 'Writing in file ',output_file
< 				
< 		i=-1
< 		j=-1
< 		new_lines=''
< 		new_data_added=False	
< 		while (True):
< 			i+=1
< 			if i>len(lines)-1:
< 				break
< 			#if 200<i:
< 			#	break
< 			line=lines[i]
< 			#for each isotope
< 			if i in iso_line_idx:
< 				#if new data was not added for last isotope
< 				#jump back to last metallicity line and add new metallicity
< 				if new_data_added==False and i>iso_line_idx[0]:
< 					i=last_metallicity_line
< 					need_to_add_z=True
< 					continue
< 				j+=1
< 				#if j>1:
< 				 #      break
< 				#if models do not have isotope needed, take one from fxt
< 				if j in no_iso and iso_line_idx[-1]>i:
< 					print 'Isotope data for isotope',isotopes[j],'taken from fxt'
< 					for old_input_line in range(i,149+i):
< 						new_lines+=lines[old_input_line]
< 					print 'jump to next isotope'
< 					print 'j:',j
< 					i=iso_line_idx[j+1]-1
< 					#break
< 					continue
< 				print 'isotope',isotopes[j]
< 				print 'set new data false'
< 				new_data_added=False
< 				need_to_add_z=False
< 
< 				to_skip_1=0
< 				to_skip_2=0
< 				to_skip_3=0
< 				label_masses_done=False
< 				new_lines+=line
< 				continue
< 			if 'big bang' in line:
< 				new_lines+=line
< 				continue
< 			#type 1a input ignore
< 			if 'type1a' in line:
< 				to_skip_1=int(line[0])
< 
< 				new_lines+=line
< 				continue
< 			if to_skip_1>0:
< 				to_skip_1-=1
< 				new_lines+=line
< 				continue
< 			#nova input ignore
< 			if 'number of nova' in line:
< 				to_skip_2=int(line[0])
< 				new_lines+=line
< 				continue
< 			if to_skip_2>0:
< 				to_skip_2-=1
< 				new_lines+=line
< 				continue
< 			#input from models
< 			#how many metallicities? 2 standard values from fxt + 
< 			if 'number of metallicity values' in line:
< 				#z_value=0
< 				#check if z is already available:
< 				following_z_array=[]
< 				#the following if statement because iso_line_idx does not contain
< 				#iso after 'remn', no further index j
< 				if iso_line_idx[-1]<i:
< 					max1=len(lines)
< 				else:
< 					max1=iso_line_idx[j+1]
< 				for t in range(i,max1):
< 					if '* metallicity' in lines[t]: 
< 						following_z_array.append(lines[t][2:12])
< 				#if z is already there
< 				print metallicity,following_z_array
< 				if metallicity in following_z_array:
< 					new_lines+=line
< 					print 'Z found'
< 					z_values=0
< 				#else add new z	
< 				else:
< 					print 'add +1 to metal number for element'
< 					new_metal_number=int(line[:5]) +1
< 					new_lines+=str(new_metal_number)+line[len(str(new_metal_number)):]
< 					z_values=int(line[0])
< 
<                                 #if first_change == True:
<                                 z_values=3
< 				continue
< 		       #go over metallicities
< 			if '* metallicity' in line:
< 				last_metallicity_line=i
< 				z_values-=1
< 				skip_masses=0
< 				add_own_masses=False
< 				if z_values>0:
< 					new_lines+=line
< 				else:
< 					#check if available Z is larger then Z from input simulation
< 					print float(metallicity),float(line[2:12]),'teststest'
< 					print line[2:12].strip()
< 					if (float(metallicity)< float(line[2:12]) and (new_data_added==False)) or need_to_add_z==True:
< 						       #add own Z
< 							print 'adding new metallicity'
< 							new_lines+='  '+metallicity+line[12:]
< 							new_data_added=True
< 							line_replaced=i
< 					#if z is already there, compare strings
< 					elif metallicity == line[2:12].strip() and (new_data_added==False):
< 						print 'location of z found'
< 						add_own_masses=True
< 						new_lines+=line
< 						new_data_added=True
< 					else:
< 						#add massive stars, M>=8
< 						new_lines+=line
< 						if first_change==False:
< 							skip_masses=11
< 						else:
< 							skip_masses=11
< 						#in remnant case
< 						print 'skipping masses'
< 				continue
< 			#go over masses
< 			if 'number of masses @ this metallicity' in line:
< 				#for available data
< 				masses_available=0
< 				if z_values>0:
< 					to_skip_3=int(line[0:5])
< 					new_lines+=line
< 					continue
< 				#for new input
< 				else:
< 					if skip_masses>0:
< 						once_done=False
< 						masses_available=int(line[:2])
< 						#when reaching the last 'isotope' remn:
< 						print masses_available,skip_masses
< 						if iso_line_idx[-1]<i:
< 							new_lines+=line
< 						else:
< 							if skip_masses==15:
< 								new_lines+=line
< 							else:
< 								if masses_available<=15:
< 									new_lines+=line
< 									once_done=True
< 								else:
< 									new_lines+=str( masses_available-skip_masses+1 )+line[len(str(masses_available-skip_masses+1)):]
< 						continue
< 					else:
< 						#if z is already there
< 						if add_own_masses==True:
< 							new_sum_masses=len(masses)+int(line[:2])
< 							new_lines+=str(new_sum_masses)+'  '+line[len(str(new_sum_masses))+2:]
< 							mass_counter=-int(line[:2])
< 						else:
< 							print 'add new Z'
< 							new_lines+=str(len(masses))+'  '+line[len(str(len(masses)))+2:]
< 							mass_counter=0
< 						continue
< 			#to skip due to available data
< 			if to_skip_3>0:
< 				to_skip_3-=1
< 				new_lines+=line
< 				continue
< 			if masses_available>0:
< 				#print 'this is executed too many times'
< 				masses_available-=1
< 				skip_masses-=1
< 				if iso_line_idx[-1]<i:
< 					new_lines+=line
< 					continue
< 				if once_done==True:
< 					new_lines+=line
< 					continue
< 				if skip_masses<1:
< 					new_lines+=line
< 				if first_change==False:
< 					new_lines+=line
< 				continue
< 			print 'this line is not executed...'
< 			#add new masses
< 			if mass_counter < (len(masses)):
< 				#if z is already there and available masses needs to be skipped to add new later on
<                                 if add_own_masses==True and mass_counter<0:
< 					mass_counter+=1
< 					new_lines+=line
< 					i-=1
< 					continue
< 				#if create new z or add new masses to available z
< 				#when reaching the last 'isotope' remn:
< 				if iso_line_idx[-1]<i:
< 					yield_1=remn_masses[mass_counter]
< 				#for all real isotopes:
< 				else:
< 					yield_1 = cols[mass_counter][j]
< 
< 				yield_line='  '+masses[mass_counter]+'  '+'{:.4E}'.format(float(yield_1))
< 				if label_masses_done==False:
< 					new_lines+=yield_line+7*' '+label_masses+'\n'
< 					label_masses_done=True
< 				else:
< 					new_lines+=yield_line+'\n'
< 				mass_counter+=1
< 				continue
< 			# jump to line with last metallicity which was replaced by new data
< 			elif mass_counter == (len(masses)):
< 				print 'jump to original'
< 				if need_to_add_z==False:
< 					i=line_replaced-1
< 					mass_counter+=1
< 				continue
< 			#all input for isotope j done, jump to line where next isotope starts
< 
< 			print 'jump to next isotope'
< 			print 'j:',j
< 			#i=iso_line_idx[j+1]-1
< 			#for rem_line in range(11771,11920):
< 			#       new_lines+=lines[rem_line]
< 			#       print new_lines
<         		break
< 
< #write output file
< 		f2=open(output_file,'w')
< 		f2.write(new_lines)
< 		f2.close()		
< 
< 
< 
< 
2660c1589
< 	def write_prod_fact_stellarwinds(self,cycles=[],isotopes=["C-12","N-14","O-16","Ne-22"],elements=[],yields_output=True,headers=['this is a test','to write strings in the first col'],file_name="isotopic_table_set1.2_prodfac.txt",ascii_1=False,latex=False,GCE_input_fxt=True):
---
> 	def write_prod_fact_stellarwinds(self,cyclerange=[[5000,8000,1000],[5000,8000,1000]],isotopes=["C-12","N-14","O-16","Ne-22"],weighted=True, label="fig:testlabel",caption="",table_header=[],file_1="test",table_size="normal",ascii_1=True,latex=False):
2667d1595
< 		case_label_and_first_column = ['specie']
2670d1597
< 		remn_masses=[]
2673d1599
< 			sefiles_hout=se(self.runs_H5_out[i])
2676a1603,1613
> 			for w in range(len(isotopes)):
> 				indx_iso.append(sefiles.se.isotopes.index(isotopes[w]))
> 
> 			if len(cyclerange)==0:
> 				cyclestart=sefiles.se.cycles[0]
> 				cycleend=sefiles.se.cycles[-1]
> 				cyclesparse=500
> 			else:
> 				cyclestart=cyclerange[i][0]
> 				cycleend=cyclerange[i][1]
> 				cyclesparse=cyclerange[i][2]
2678,2712c1615,1623
< 			metallicity=sefiles.get("zini")[0]
< 			case_label_and_first_column.append(str(mass)+" Msun, Z="+str(metallicity))
< 			##if all stable ones:
< 			if isotopes[0]=="all":
< 				isotopes_1=sefiles.se.isotopes
< 				isotopes=[]
< 				for j in range(len(isotopes_1)):
< 					if is_stable(isotopes_1[j])=='t':
< 						isotopes.append(isotopes_1[j])			
< 			print isotopes		
< 			if len(elements)>0:
< 				if elements[0]=="all":
< 					elements=sefiles.se.elements
< 					elements=get_stable(sefiles.se.isotopes,get_elements=True)[:12]		
< 					#exclude neutrons
< 					print elements
<                         if cycles[i][1]==-1:
<                         	endcycle=int(sefiles.se.cycles[-1]) #+ cycles[k][2]  #1000
<                                 print endcycle
<                                	endcycle= int( round( endcycle,-3) -2000 )
<                                 print "endcycle",endcycle
<                         else:
<                                 endcycle=cycles[i][1]
< 			prod_factor,isotopes_prod_fac,yields,iniabu,remn_mass =self.weighted_yields(sefiles,sefiles_hout,isotopes,elements,cycles[i][0],endcycle,cycles[i][2])
< 			#print 'isotopes',isotopes_prod_fac
< 			#print 'prod factor:',prod_factor
< 			#print 'yield:',yields
< 			remn_masses.append(remn_mass)	
< 			if yields_output==False:	
< 				yields_all.append(list(prod_factor))
< 			else:
< 				yields_all.append(list(yields))
< 		cols=yields_all
< 		print cols
< 		print case_label_and_first_column
---
> 		    	yield_folded,yield_unfolded=self.weighted_yields(sefiles=sefiles,cyclestart=cyclestart,cycleend=cycleend,sparse=cyclesparse,star_mass=mass)
> 			for i in indx_iso:
> 				iso_names.append(sefiles.se.isotopes[i])
> 				if weighted==True:
> 					yields_1.append(yield_folded[i])
> 				else:
> 					yields_1.append(yield_unfolded[i])	
> 			yields_all.append(yields_1)
> 		cols=[iso_names] + yields_all
2715,2717c1626
< 			self.write_latex_table_set1(isotopes_prod_fac,cols,headers,file_name,case_label_and_first_column)
< 
< 			#self.write_latex_table(header_file="Weighted stellar winds",cols=cols,label=label,caption=caption,table_header=table_header,latexfile=file_1,table_size=table_size)
---
> 			self.write_latex_table(header_file="Weighted stellar winds",cols=cols,label=label,caption=caption,table_header=table_header,latexfile=file_1,table_size=table_size)
2721,2755d1629
< 		#for GCE model of fxt	
< 		if GCE_input_fxt==True:
< 			return case_label_and_first_column,cols,remn_masses	
< 	
< 
< 
< 	def write_latex_table_set1(self,isotopes_for_table,data,headers=['this is a test','to write strings in the first col'],file_name_table_species_wind = 'isotopic_table_set1.2_prodfac.txt',case_label_and_first_column = ['specie','1.65 Msun','2 Msun','3 Msun','5 Msun','15 Msun','20 Msun','25 Msun',]):
< 
< 		'''
< 			Table write method from Paper1 scripts (from Marco)
< 
< 			isotopes_for_table - isotopes for table (first column), preferable stable
< 			data - 2d array e.g. [tab_w_1p65,tab_w_2,tab_w_3,tab_w_5,tab_w_15,tab_w_20,tab_w_25]
< 
< 
< 		'''
< 
< 
< 		# headers and format
< 		form_str='%7.3E'
< 		sep_string=' & '
< 
< 		all_data=[]
< 		all_data.append(isotopes_for_table)
< 		for i in range(len(data)):
<         		datal=list(form_str%data[i][j] for j in range(len(data[i])))
<         		all_data.append(datal)
< 
< 
< 		### attempt to add the trailing '\\' to a latex table
< 		final_col=len(isotopes_for_table)*[r'\\ ']
< 		case_label_and_first_column.append(r'\\ ')
< 		all_data.append(final_col)
< 		att.write(file_name_table_species_wind,headers,case_label_and_first_column,all_data,sep=sep_string)
< 
2844c1718
< 			row_type	3Msun Z=...	4Msun Z=...
---
> 
2883,2884c1757
< 				a="%.3E" % a
<                                 line+= "	"+ a		#str(round(float(str(a).split("e")[0]),3))+"E"+str(a).split("e")[1]
---
>                                 line+= "	"+str(round(float(str(a).split("e")[0]),3))+"E"+str(a).split("e")[1]
3364,3366c2237,2239
<      		   	This function is under construction....will be done after my thesis.
< 			Idea: time dependent prod factors and yields..
< 			
---
> 				Uses H5_surf
> 	  	    This function returns the wind yields and ejected masses for stable species.  The command
>      		   is
3457,3459c2330,2332
< 	  	    This function returns the production factors, isotope names,yields and initial
< 		    abundances of isotopes or elements
< 		    If elements>0 than ignore isotope input. 
---
> 				Uses H5_surf
> 	  	    This function returns the wind yields and ejected masses for stable species.  The command
>      		   is
3467,3468c2340
< 		print "start,etc:",cyclestart,cycleend,sparse
< 		print "mini get done, second part"	
---
> 		print "start,etc:",cyclestart,cycleend,sparse	
3469a2342
> 		
3497c2370
< 				c12=self.get_elem_abu(sefiles,cycs,elements[i])
---
> 				c12=sefiles.get(cycs,'elem_massf',elements[i])
3503,3505c2376,2377
< 				iniabu_1=self.get_elem_abu(sefiles,0,elements[i])
< 				envelope_yield=envelope_mass*self.get_elem_abu(sefiles,cycleend,elements[i])
< 
---
> 				iniabu_1=sefiles.get(0,'elem_massf',elements[i])
> 				envelope_yield=envelope_mass*sefiles.get(cycleend,'elem_massf',elements[i])
3546,3547c2418
< 		remn_mass=h_free_core
< 		return np.array(prod_facs),isotopes,yields,iniabus,remn_mass
---
> 		return np.array(prod_facs),isotopes,yields,iniabus
3779,3966c2650,2651
< 
< 
< 	def surface_plot(self,runs=[],extra_label=[],steps=1000,timestep=[],color=["r","b","k","g","r","b","k","g","r"],line_style=["-","--","-.","-.","--","-"],setting='iso_massf_decay',iniabupath='/home/christian/PPN/forum.astro.keele.ac.uk/frames/mppnp/USEEPP/iniab2.0E-02GN93.ppn'):
< 
< 		'''
< 			Plots surface abus
< 
< 		'''	
< 		import utils as u
<                 sefiles=[]
<                 legend=[]
<                 HDF5_surf=[]
<                 if len(runs) ==0:
<                         runs=self.run_dirs_name
<                         if len(extra_label)==0:
<                                 extra_label=self.extra_label
< 		for i in range(len(self.runs_H5_surf)):
< 			sefiles.append(se(self.runs_H5_surf[i]))
< 
< 		hs=["Ba","La","Nd","Sm"]
< 		ls=["Sr","Y-","Zr"]
< 		#get initial abu	
< 		ba_solar=0
< 		la_solar=0
< 		nd_solar=0
< 		sm_solar=0
< 		sr_solar=0
< 		y_solar=0
< 		zr_solar=0
< 		fe_solar=0
< 		pb_solar=0
< 		a=u.iniabu(iniabupath)
< 		for iso in a.names:
< 			iso_namescheme=iso.capitalize()[:2].split()[0]+"-"+iso[3:].split()[-1]	
< 			if "Ba-" in iso_namescheme:
< 				ba_solar+=a.habu[iso]
< 			if "La-" in iso_namescheme:
< 				la_solar+=a.habu[iso]
< 			if "Nd-" in iso_namescheme:
< 				nd_solar+=a.habu[iso]
< 			if "Sm-" in iso_namescheme:
< 				sm_solar+=a.habu[iso]
< 			if "Sr-" in iso_namescheme:
< 				sr_solar+=a.habu[iso]
< 			if "Y-" in iso_namescheme:
< 				y_solar+=a.habu[iso]
< 			if "Zr-" in iso_namescheme:
< 				zr_solar+=a.habu[iso]
< 			if "Fe-" in iso_namescheme:
< 				fe_solar+=a.habu[iso]
< 			if "Pb-" in iso_namescheme:
< 				pb_solar+=a.habu[iso]
< 
< 		for i in range(len(sefiles)):
< 			print "open",sefiles[i]
< 
< 			ls_fe_abu=[]
< 			hs_ls_abu=[]
< 			pb_hs_abu=[]
< 			hs_fe_abu=[]
< 			star_mass=sefiles[i].get("mini")[0]
<                         star_z=sefiles[i].get("zini")[0]
< 			legend=str(star_mass)+"$M_{\odot}$, Z="+str(star_z)
< 			cycles=[]
< 			if timestep>0:
< 				cycles.append(sefiles[i].se.cycles[0])
< 				timesteps=timestep
< 				for j in range(0,len(sefiles[i].se.cycles),1):
< 					if sefiles[i].get('star_age')[j]>timesteps:
< 						cycles.append(sefiles[i].se.cycles[j])				
< 						timesteps+=timestep
< 						print 'take cycle',cycles[j]
< 			else:
< 				for j in range(0,len(sefiles[i].se.cycles),steps):
< 					cycles.append(sefiles[i].se.cycles[j])					
< 			
< 
< 			ba_abu=np.zeros(len(cycles))
< 			la_abu=np.zeros(len(cycles))
< 			nd_abu=np.zeros(len(cycles))
<                         sm_abu=np.zeros(len(cycles))
<                         sr_abu=np.zeros(len(cycles))
<                         y_abu=np.zeros(len(cycles))
<                         zr_abu=np.zeros(len(cycles))
<                         fe_abu=np.zeros(len(cycles))
< 			pb_abu=np.zeros(len(cycles))	
< 			
< 			k=-1
< 			for cyc in cycles:
< 				k+=1
< 				print "read cycle",cyc
< 				if "iso" in setting:
< 					for iso_namescheme in sefiles[i].se.isotopes:
< 						if "Ba-" in iso_namescheme:
< 							print iso_namescheme
< 							ba_abu[k]+=sefiles[i].get(cyc,setting,iso_namescheme)
< 						if "La-" in iso_namescheme:
< 							la_abu[k]+=sefiles[i].get(cyc,setting,iso_namescheme)
< 							print iso_namescheme
< 						if "Nd-" in iso_namescheme:
< 							nd_abu[k]+=sefiles[i].get(cyc,setting,iso_namescheme)
< 							print iso_namescheme
< 						if "Sm-" in iso_namescheme:
< 							sm_abu[k]+=sefiles[i].get(cyc,setting,iso_namescheme)
< 							print iso_namescheme
< 						if "Sr-" in iso_namescheme:
< 							sr_abu[k]+=sefiles[i].get(cyc,setting,iso_namescheme)
< 							print iso_namescheme
< 						if "Y-" in iso_namescheme:
< 							y_abu[k]+=sefiles[i].get(cyc,setting,iso_namescheme)
< 							print iso_namescheme
< 						if "Zr-" in iso_namescheme:
< 							zr_abu[k]+=sefiles[i].get(cyc,setting,iso_namescheme)
< 							print iso_namescheme
< 						if "Fe-" in iso_namescheme:					
< 							fe_abu[k]+=sefiles[i].get(cyc,setting,iso_namescheme)
< 							print iso_namescheme
< 						if "Pb-" in iso_namescheme:
< 							pb_abu[k]+=sefiles[i].get(cyc,setting,iso_namescheme)
< 				else:
< 					for iso_namescheme in sefiles[i].se.elements:
<                                                 if "Ba" in iso_namescheme:
<                                                         ba_abu[k]+=sefiles[i].get(cyc,setting,iso_namescheme)
<                                                 if "La" in iso_namescheme:
<                                                         la_abu[k]+=sefiles[i].get(cyc,setting,iso_namescheme)
<                                                 if "Nd" in iso_namescheme:
<                                                         nd_abu[k]+=sefiles[i].get(cyc,setting,iso_namescheme)
<                                                 if "Sm" in iso_namescheme:
<                                                         sm_abu[k]+=sefiles[i].get(cyc,setting,iso_namescheme)
<                                                 if "Sr" in iso_namescheme:
<                                                         sr_abu[k]+=sefiles[i].get(cyc,setting,iso_namescheme)
<                                                 if "Y" == iso_namescheme:
<                                                         y_abu[k]+=sefiles[i].get(cyc,setting,iso_namescheme)
<                                                 if "Zr" in iso_namescheme:
<                                                         zr_abu[k]+=sefiles[i].get(cyc,setting,iso_namescheme)
<                                                 if "Fe" in iso_namescheme:
< 							fe_abu[k]+=sefiles[i].get(cyc,setting,iso_namescheme)
< 						#pb needs to be added here		
< 	
< 				#print fe_abu[k],sr_abu[k],zr_abu[k],y_abu[k]
< 				#print sr_solar,zr_solar,y_solar,fe_solar
< 				print "PB:",pb_solar,pb_abu[k]
< 				ls_fe_abu.append( np.log10((fe_solar/fe_abu[k])*((sr_abu[k] * zr_abu[k] * y_abu[k]/(sr_solar * zr_solar * y_solar))**(1./3.) )))
< 				#print ((fe_solar/fe_abu[k])*((sr_abu[k] * zr_abu[k] * y_abu[k]/(sr_solar * zr_solar * y_solar))**(1./3.) ))
< 				ls_ba=np.log10((ba_solar/ba_abu[k])*((sr_abu[k] * zr_abu[k] * y_abu[k]/(sr_solar * zr_solar * y_solar))**(1./3.) ))
< 				ls_la=np.log10((la_solar/la_abu[k])*((sr_abu[k] * zr_abu[k] * y_abu[k]/(sr_solar * zr_solar * y_solar))**(1./3.) ))
< 				ls_nd=np.log10((nd_solar/nd_abu[k])*((sr_abu[k] * zr_abu[k] * y_abu[k]/(sr_solar * zr_solar * y_solar))**(1./3.) ))
< 				ls_sm=np.log10((sm_solar/sm_abu[k])*((sr_abu[k] * zr_abu[k] * y_abu[k]/(sr_solar * zr_solar * y_solar))**(1./3.) ))
< 				hs_fe=np.log10((fe_solar/fe_abu[k])*((ba_abu[k] * la_abu[k] * nd_abu[k] * sm_abu[k]/( ba_solar * la_solar * nd_solar*sm_solar))**(1./4.) ))	
< 				hs_pb=np.log10((pb_solar/pb_abu[k])*((ba_abu[k] * la_abu[k] * nd_abu[k] * sm_abu[k]/( ba_solar * la_solar * nd_solar*sm_solar))**(1./4.) ))
< 				hs_fe_abu.append(hs_fe)
< 				hs_ls_abu.append(hs_fe - ls_fe_abu[-1])   #( -1*(ls_ba) + -1*(ls_la) + -1*(ls_nd) + -1*(ls_sm))/4.)
< 				pb_hs_abu.append( -1* hs_pb )
< 				#ls_fe_abu.append( np.log10( 1./3. *( sr_abu[k]/sr_solar + y_abu[k]/y_solar + zr_abu[k]/zr_solar)/(fe_abu[k]/fe_solar))) 
< 				#hs_ls_abu.append( np.log10( 1./4. *( ba_abu[k]/ba_solar + la_abu[k]/la_solar + nd_abu[k]/nd_solar + sm_abu[k]/sm_solar) / ( 1./3. *( sr_abu[k]/sr_solar + y_abu[k]/y_solar + zr_abu[k]/zr_solar)))) 
< 				
< 
< 
< 			print hs_ls_abu
< 			print "------------------"
< 			print ls_fe_abu
< 			print "plot",legend
< 			fig1=plt.figure(333)
< 			plt.rcParams.update({'font.size': 16})
<                 	plt.rc('xtick', labelsize=16)
<                 	plt.rc('ytick', labelsize=16)
<                 	plt.plot([0,0],[-10,10],linewidth=2,color="K")
<                 	plt.plot([-10,10],[0,0],linewidth=2,color="K")
< 			ax=fig1.add_subplot(1,1,1)
< 			ax.plot(ls_fe_abu,hs_fe_abu,marker='<',label=legend)
< 			plt.xlabel("[ls/Fe]")
< 			plt.ylabel("[hs/Fe]")
< 			plt.legend()
< 			
< 			fig2=plt.figure(555)
<                         plt.rcParams.update({'font.size': 16})
<                         plt.rc('xtick', labelsize=16)
<                         plt.rc('ytick', labelsize=16)
<                         plt.plot([0,0],[-10,10],linewidth=2,color="K")
<                         plt.plot([-10,10],[0,0],linewidth=2,color="K")
< 			ax2=fig2.add_subplot(1,1,1)
< 			ax2.plot(ls_fe_abu,pb_hs_abu,marker='<',label=legend)
< 			plt.xlabel("[ls/Fe]")
< 			plt.ylabel("[Pb/hs]")
< 			plt.legend()
< 
< 
< 
---
>  
> 	
4345,4367d3029
< 	def get_elem_abu(self,sefiles,cycles,element):
< 		'''
< 			Get abundance of element without using elem_massf since
< 			its not working correctly...
< 			reads also hdf5 out files
< 		'''
< 		isotopes=sefiles.se.isotopes
< 		if type(cycles)==int:
< 			len_cycles=1
< 		else:
< 			len_cycles=len(cycles)
< 		type_hdf5= type(sefiles.get(sefiles.se.cycles[0],"H-1"))
< 		print type_hdf5
< 		if type_hdf5 == np.float64:
< 		 	abu=np.zeros(len_cycles)
< 			for i in range(len(isotopes)):
< 				if element in isotopes[i].split("-")[0]: 	
< 					abu+=sefiles.get(cycles,isotopes[i])
< 		## hdf5 out
< 		#else:
< 			#abu
< 
< 		return abu	
4444c3106
< 					print "Warning: h5 files are not available in "+run_path
---
> 					print "Warning: h5 files are not available in "+run_path+"/"+h5_dir
4703,4726c3365
< 
< 	def lifetime(self,label=""):
< 		'''
< 			Calculate stellar lifetime till first TP 
< 			dependent of initial mass
< 		'''
<                 plt.rcParams.update({'font.size': 20})
<                 plt.rc('xtick', labelsize=20)
<                 plt.rc('ytick', labelsize=20)
<                 t0_model=self.set_find_first_TP()
<                 m=self.run_historydata
<                 i=0
< 		age=[]
< 		mass=[]
<                 for case in m:
< 			##lifetime till first TP
< 			age.append(np.log10(case.get("star_age")[t0_model[i]]))
< 			mass.append(case.get("star_mass")[0])
< 			i+=1
< 		plt.plot(mass,age,"*",markersize=10,linestyle="-",label=label)
< 		plt.xlabel('star mass $[M_{\odot}]$',fontsize=20)
< 		plt.ylabel('Logarithmic stellar lifetime',fontsize=20)	
< 		
< 		
---
> 	
4734c3373
< 		print "M_ini: | M_h1 | M_he4  | M_final (withoutloss!)"
---
> 		print "M_ini: | M_h1 | M_he4 | M_final (withoutloss!)"
4739,4744c3378,3379
< 			#if case.header_attr["initial_mass"]>=5: 
< 			#	c12bndy=case.get('c12_boundary_mass')
< 			#else:
< 			#	c12bndy=0
< 			print star_mass[0],"|",h1bndy[-1],"|",he4bndy[-1],"|","|",star_mass[-1]
< 
---
> 			#.get('c12_boundary_mass')
> 			print star_mass[0],"|",h1bndy[-1],"|",he4bndy[-1],"|",star_mass[-1]
4746c3381
< 	def TPAGB_core_growth(self,fig, label="",color='k',marker_type='<',linestyle='-'):
---
> 	def TPAGB_core_growth(self, label="",color='k',marker_type='<',linestyle='-'):
4758c3393
< 		plt.figure(fig)
---
> 		plt.figure()
4864,4868c3499
< 	def set_plot_CO(self,startfirstTP=True,xtime=True,label=[]):
< 
< 	
< 		color=['r','b','k','g']
< 
---
> 	def set_plot_CO(self):
4870,4882d3500
<                 plt.rcParams.update({'font.size': 20})
<                 plt.rc('xtick', labelsize=20)
<                 plt.rc('ytick', labelsize=20)
< 
<                 m=self.run_historydata
<                 i=0
<                 if startfirstTP==True:
<                         t0_model=self.set_find_first_TP()
< 
<                 else:
<                         t0_model=len(self.run_historydata)*[0]
< 
< 
4885,4886d3502
< 		if len(label)==0:
< 			label=self.run_label[i]
4888,4890c3504,3507
< 			model_number=case.get('model_number')[t0_model[i]:]
<         		C=case.get('surface_c12')[t0_model[i]:]
<         		O=case.get('surface_o16')[t0_model[i]:]
---
> 	       		star_age=case.get('star_age')
> 			model_number=case.get('model_number')
>         		C=case.get('surface_c12')
>         		O=case.get('surface_o16')
4892,4897c3509
< 			age=case.get("star_age")[t0_model[i]:]
<                         t0_age=age[0]
<                         age=age-t0_age
< 			if xtime==True:
< 				model_number=age
<         		plot(model_number,CO,'-',color=color[i],label=label[i])
---
>         		plot(model_number[noffset:],CO[noffset:],symbs[i],label=self.run_label[i])
4900,4903c3512
< 		if xtime==True:
< 			plt.xlabel("Star age [yr]")
< 		else:
< 			xlabel('model number')
---
> 		xlabel('model number')
4960,4968c3569,3570
< 			plt.figure(55)
< 			model=case.get("model_number")
< 			vdiv=case.get("v_div_csound_surf")
< 			plt.plot(model,vdiv,label=self.run_label[i])
< 			plt.legend()
< 			plt.figure(i)
< 			plt.plot(model,vdiv,label=self.run_label[i])
< 			plt.legend()
<         		i += 1
---
> 			case.plot('model_number','v_div_csound_surf',legend=self.run_label[i],shape=symbs[i])
>         	i += 1
4971a3574,3575
>     		if xlim_mod_min >= 0:
>         		xlim(xlim_mod_min,xlim_mod_max)
4974,4991d3577
<         def set_plot_timesteps(self):
<                 m=self.run_historydata
<                 figure(55)
<                 i=0
<                 for case in m:
<                         plt.figure(55)
<                         model=case.get("model_number")
<                         vdiv=case.get("log_dt")
<                         plt.plot(model,vdiv,label=self.run_label[i])
<                         plt.legend()
<                         plt.figure(i)
<                         plt.plot(model,vdiv,label=self.run_label[i])
<                         plt.legend()
<                         i += 1
<                 legend(loc=2)
<                 xlabel('model number')
<                 ylabel('Log(dt)')
< 
4996,4999d3581
<                 plt.rcParams.update({'font.size': 20})
<                 plt.rc('xtick', labelsize=20)
<                 plt.rc('ytick', labelsize=20)
<                 
5001c3583
< 		if xaxis=="time":
---
>                 if xaxis=="time":
5029,5030c3611,3612
< 			xlabel('M/M$_{\odot}$',fontsize=22)
<     		ylabel('log($\|\dot{M}\|$)',fontsize=22)
---
> 			xlabel('"M/M$_{\odot}$"')
>     		ylabel('log_abs_mdot')
5045,5103c3627
< 
< 	def set_plot_kip_special(self,startfirstTP=True,xtime=True,label=[],color=[]):
< 		'''
< 			Kippenhahn which plots only h and he free bndry,
< 			label and color can be chosen.
< 			if label>0 then color must be set too!
< 			color=["r","b","g","k"]
< 			
< 
< 		'''
< 		plt.rcParams.update({'font.size': 24})
< 		plt.rc('xtick', labelsize=24)
< 		plt.rc('ytick', labelsize=24)
< 
<                 m=self.run_historydata
<                 i=0
<                 if startfirstTP==True:
<                         t0_model=self.set_find_first_TP()
< 
<                 else:
<                 	t0_model=len(self.run_historydata)*[0]
< 
< #                if xtime==True:
<  #                       for case in m:
< #				
<  #                               case.kippenhahn(i,'time',t0_model=t0_model[i],c12_bm=False)
<    #                             title(self.run_label[i])
<   #                              i += 1
<     #            else:
< 		#color=["r","b","g","k"]
< 		for case in m:
< 			h1_boundary_mass  = case.get('h1_boundary_mass')[t0_model[i]:]
< 			he4_boundary_mass = case.get('he4_boundary_mass')[t0_model[i]:]
< 			star_mass         = case.get('star_mass')[t0_model[i]:]
< 			mx1_bot           = case.get('mx1_bot')[t0_model[i]:]*star_mass
< 			model = case.get("model_number")[t0_model[i]:]
< 			age=case.get("star_age")[t0_model[i]:]
< 			t0_age=age[0]
< 			age=age-t0_age
< 			if xtime==True:
< 				model=age
< 			if len(label)>0:
< 				plt.plot(model,h1_boundary_mass,color=color[i],label="$^1$H bndry, "+label[i])
< 				plt.plot(model,he4_boundary_mass,"--",color=color[i],label="$^4$He bndry, "+ label[i])
< 				#plt.plot(model,star_mass,color=color[i],label="Total mass")
< 			else:
< 				plt.plot(model,h1_boundary_mass,label="h1_boundary_mass")
< 				plt.plot(model,he4_boundary_mass,"--",label="he4_boundary_mass")
< 			#plt.plot(model,mx1_bot,label="convective boundary")
< 			#title(self.run_label[i])
< 			i += 1
< 		if xtime==True:
< 			plt.xlabel('stellar age',size=28)
< 			if startfirstTP==True:
< 				plt.xlabel('t - t$_0$ $\mathrm{[yr]}$',size=28)
< 		else:	
< 			plt.xlabel('model number',size=28)
< 		plt.ylabel("M/M$_{\odot}$",size=28)		
< 		plt.legend()	
---
> 	
5238,5259c3762
< 
< 
< def get_stable(specie_list,get_elements=True):
< 
< 	'''
< 		Input isotope list or element list, get stable list back
< 		if get_elements=True return elements of all stables in isotope_list
< 
< 	'''
< 
< 	stable_list=[]
< 	for specie in specie_list:
< 		if is_stable(specie) == 't':
< 			if get_elements==True and  len(specie.split("-"))==2:
< 				if specie.split("-")[0] not in stable_list:
< 					stable_list.append(specie.split("-")[0])	
< 			else:
< 				stable_list.append(specie)
< 	
< 	return stable_list
< 
< 
---
> 			
5261,5266c3764
< 
< 	'''
< 		Input either isotope e.g. "C-12" or element "C".
< 		if stable returns "t"
< 		Recent correction in this list: Na-23 was set f
< 	'''                
---
>                 
5273c3771
< 	"Na-23":"t","Mg-23":"f","Mg-24":"t","Mg-25":"t","Mg-26":"t",
---
> 	"Na-23":"f","Mg-23":"f","Mg-24":"t","Mg-25":"t","Mg-26":"t",
5486,5500c3984,3989
< 
< 	##input is isotope
< 	if len(isotope.split("-"))==2:
< 		isotope="".join(isotope.split())
< 		for i in range(len(isotopes.keys())):
<                 	#print "|"+isotope+stable_list[i][0]+"|"
< 			if isotope == isotopes.keys()[i]:
< 				return isotopes[isotope]
< 	#if element
< 	if len(isotope.split("-"))==1:
< 		for i in range(len(isotopes)):
< 			if isotope == isotopes.keys()[i].split("-")[0]  and isotopes[isotopes.keys[i]] == "t":
< 				return "t"
< 		return "f"
< 
---
> 	
> 	isotope="".join(isotope.split())
> 	for i in range(len(isotopes.keys())):
>                 #print "|"+isotope+stable_list[i][0]+"|"
> 		if isotope == isotopes.keys()[i]:
> 			return isotopes[isotope]