docs: def plot_surf_composition(self,hsfe_hsls = True,lsfe_hsls = False,rbsr_hsls = False,rbzr_hsls = False,sry_srzr = False,rbfe_sfe = False,iso_ratio= False,isotopes_x = ['Mg-25','Mg-24'],isotopes_y = ['Mg-26','Mg-24'],sparsity = 5000,marker1=['o'],symbols1=['k-'],color=['k']):
  		'''
  			Set only one of the vriables True to plot quantitie4s.
  			In case you set iso_ratio = True:
  				isotooes_x and isototopes_y are used.  
  		'''
  
  		import nugridse as mp
  		import matplotlib.pyplot as pl
  		import sys
  
  		plotting=True
  
  		### SET THIS IF PLOTTING = TRUE ####################
  		###>>>>!!!! ONLY ONE OPTION CAN BE TRUE !!!!<<<< ###
  		#test=[hsfe_hsls,lsfe_hsls,
  		#hsfe_hsls = False
  		#lsfe_hsls = False
  		#rbsr_hsls = False
  		#rbzr_hsls = False
  		#sry_srzr = False
  		#rbfe_sfe = False
  		#iso_ratio= True
  		###################################
  
                  for i in range(len(self.runs_H5_surf)):
                         	m3z2m2=se(self.runs_H5_surf[i])
  			run_mass=[m3z2m2]
  
  			markers=[marker1[i]]#['ks']#,'r^','gh','bd','mv','cH','kD']
  			symbol=[symbols1[i]]#['k-']#,'r--','g-','b:','m-','c--','k-']
  			color1=color[i]
  			ls_element = ['Sr-84','Sr-86','Sr-87','Sr-88','Y-89','Zr-90','Zr-91','Zr-92','Zr-94','Zr-96']
  			hs_element = ['Ba-130','Ba-132','Ba-134','Ba-135','Ba-136','Ba-137','Ba-138','La-138','La-139','Nd-142','Nd-143','Nd-144','Nd-145','Nd-146','Nd-148','Nd-150','Sm-144','Sm-147','Sm-148','Sm-149','Sm-150','Sm-152','Sm-154']
  			# notice that Pb (3rd s-process index) is not included here.
  			s_element =  ['Sr-84','Sr-86','Sr-87','Sr-88','Y-89','Zr-90','Zr-91','Zr-92','Zr-94','Zr-96','Ba-130','Ba-132','Ba-134','Ba-135','Ba-136','Ba-137','Ba-138','La-138','La-139','Nd-142','Nd-143','Nd-144','Nd-145','Nd-146','Nd-148','Nd-150','Sm-144','Sm-147','Sm-148','Sm-149','Sm-150','Sm-152','Sm-154']	
  
  			ls_element = ['Sr-84','Sr-86','Sr-87','Sr-88','Y-89','Zr-90','Zr-91','Zr-92','Zr-94','Zr-96']
  			hs_element = ['Ba-130','Ba-132','Ba-134','Ba-135','Ba-136','Ba-137','Ba-138','La-138','La-139','Nd-142','Nd-143','Nd-144','Nd-145','Nd-146','Nd-148','Nd-150','Sm-144','Sm-147','Sm-148','Sm-149','Sm-150','Sm-152','Sm-154']
  			# notice that Pb (3rd s-process index) is not included here.
  			s_element =  ['Sr-84','Sr-86','Sr-87','Sr-88','Y-89','Zr-90','Zr-91','Zr-92','Zr-94','Zr-96','Ba-130','Ba-132','Ba-134','Ba-135','Ba-136','Ba-137','Ba-138','La-138','La-139','Nd-142','Nd-143','Nd-144','Nd-145','Nd-146','Nd-148','Nd-150','Sm-144','Sm-147','Sm-148','Sm-149','Sm-150','Sm-152','Sm-154']	
  
  			s_ini = []
  			ls_ini = []
  			hs_ini = []
  			pb_ini = []
  			fe_ini = []
  			rb_ini = []
  			sr_ini = []
  			zr_ini = []
  			y_ini = []
  			zra_ini = []
  			zrb_ini = []
  			gda_ini = []
  			gdb_ini = []
  
  			if plotting:
  
  				if iso_ratio:
  
  					isotopes_x = ['Mg-25','Mg-24']
  					isotopes_y = ['Mg-26','Mg-24']
  
  					#isotopes_x = ['Zr-96','Zr-94']
  					#isotopes_y = ['Gd-152','Gd-154']
  					#isotopes_y = ['Zr-90','Zr-94']
  					#isotopes_y = ['Zr-91','Zr-94']
  					#isotopes_y = ['Zr-92','Zr-94']
  
  					#isotopes_x = ['Ba-135','Ba-136']
  					#isotopes_y = ['Ba-138','Ba-136']
  					#isotopes_y = ['Ba-134','Ba-136']
  					#isotopes_y = ['Ba-137','Ba-136']
  
  					#isotopes_x = ['Sr-84','Sr-86']
  					#isotopes_y = ['Sr-87','Sr-86']
  					#isotopes_y = ['Sr-88','Sr-86']
  
  					# get initial value
  					initial_isotopic_ratio_x = float(run_mass[0].se.get(min(run_mass[0].se.cycles),'iso_massf',isotopes_x[0]))/float(run_mass[0].se.get(min(run_mass[0].se.cycles),'iso_massf',isotopes_x[1]))
  					initial_isotopic_ratio_y = float(run_mass[0].se.get(min(run_mass[0].se.cycles),'iso_massf',isotopes_y[0]))/float(run_mass[0].se.get(min(run_mass[0].se.cycles),'iso_massf',isotopes_y[1]))
  
  					# sparcity for cycles I am looking at.
  					#sparsity = 5000
  
  					I_want_iso_ratio = False
  					I_want_delta = True
  
  					isotopic_ratio_x  = []
  					isotopic_ratio_y = []
  					isotopic_ratio_x_tps = []
  					isotopic_ratio_y_tps = []
  					isotopic_ratio_x_tps_co = []
  					isotopic_ratio_y_tps_co = []
  
  					k = 0
  
  					for i in run_mass:
  						dum_isotopic_ratio_x = []
  						dum_isotopic_ratio_y = []
  						dumdum_isotopic_ratio_x = []
  						dumdum_isotopic_ratio_y = []
  						dum_isotopic_ratio_x_tps = []
  						dum_isotopic_ratio_y_tps = []
  						dum_isotopic_ratio_x_tps_co = []
  						dum_isotopic_ratio_y_tps_co = []
  						co_ratio=[]
  						for j in i.se.cycles[0::sparsity]:
  							print 'j= ', j
  							dum_isotopic_ratio_x.append(float(i.se.get(j,'iso_massf',isotopes_x[0]))/float(i.se.get(j,'iso_massf',isotopes_x[1])))
  							dum_isotopic_ratio_y.append(float(i.se.get(j,'iso_massf',isotopes_y[0]))/float(i.se.get(j,'iso_massf',isotopes_y[1])))
  							co_ratio.append((float(i.se.get(j,'iso_massf','C-12')*4.))/(float(i.se.get(j,'iso_massf','O-16'))*3))
  
  							if (len(co_ratio)>1):
  								if (co_ratio[len(co_ratio)-1]>(co_ratio[len(co_ratio)-2]+0.02)):
  									dum_isotopic_ratio_x_tps.append(float(i.se.get(j,'iso_massf',isotopes_x[0]))/float(i.se.get(j,'iso_massf',isotopes_x[1])))
  									dum_isotopic_ratio_y_tps.append(float(i.se.get(j,'iso_massf',isotopes_y[0]))/float(i.se.get(j,'iso_massf',isotopes_y[1])))
  
  									if (co_ratio[len(co_ratio)-1]>1.):
  										dum_isotopic_ratio_x_tps_co.append(float(i.se.get(j,'iso_massf',isotopes_x[0]))/float(i.se.get(j,'iso_massf',isotopes_x[1])))
  										dum_isotopic_ratio_y_tps_co.append(float(i.se.get(j,'iso_massf',isotopes_y[0]))/float(i.se.get(j,'iso_massf',isotopes_y[1])))
  
  							if I_want_delta:
  								dumdum_isotopic_ratio_x=(np.array(dum_isotopic_ratio_x)/initial_isotopic_ratio_x-1.)*1000.
  								dumdum_isotopic_ratio_y=(np.array(dum_isotopic_ratio_y)/initial_isotopic_ratio_y-1.)*1000.
  								dumdum_isotopic_ratio_x_tps=(np.array(dum_isotopic_ratio_x_tps)/initial_isotopic_ratio_x-1.)*1000.
  								dumdum_isotopic_ratio_y_tps=(np.array(dum_isotopic_ratio_y_tps)/initial_isotopic_ratio_y-1.)*1000.
  								dumdum_isotopic_ratio_x_tps_co=(np.array(dum_isotopic_ratio_x_tps_co)/initial_isotopic_ratio_x-1.)*1000.
  								dumdum_isotopic_ratio_y_tps_co=(np.array(dum_isotopic_ratio_y_tps_co)/initial_isotopic_ratio_y-1.)*1000.
  
  						if I_want_iso_ratio:
  							isotopic_ratio_x.append(dum_isotopic_ratio_x)
  							isotopic_ratio_y.append(dum_isotopic_ratio_y)
  						else:
  							isotopic_ratio_x.append(dumdum_isotopic_ratio_x)
  							isotopic_ratio_y.append(dumdum_isotopic_ratio_y)
  							isotopic_ratio_x_tps.append(dumdum_isotopic_ratio_x_tps)
  							isotopic_ratio_y_tps.append(dumdum_isotopic_ratio_y_tps)
  							isotopic_ratio_x_tps_co.append(dumdum_isotopic_ratio_x_tps_co)
  							isotopic_ratio_y_tps_co.append(dumdum_isotopic_ratio_y_tps_co)
  
  						k = k+1
  
  					markersss=['ko','b^','rh','gd','mv','cH','kD']
  					symbol=['k-','b--','r-.','g:','m-','c--','k-']
  
  					mass_label =[]
  					for i in run_mass:
  						mini=float(i.se.get('mini'))
  						zini=float(i.se.get('zini'))
  						label=str(mini)+'$M_{\odot}$, Z='+str(zini)
  						mass_label.append(label)
  
  					metallicity_label=['M3z2m2','M3z2m2_nomol','M3z1m2','3 Msun, Z=0.02, PI13']
  					metallicity_label=[label]
  
  					array_to_plot_x = isotopic_ratio_x
  					array_to_plot_y = isotopic_ratio_y
  					array_to_plot_x_tps = isotopic_ratio_x_tps
  					array_to_plot_y_tps = isotopic_ratio_y_tps
  					array_to_plot_x_tps_co = isotopic_ratio_x_tps_co
  					array_to_plot_y_tps_co = isotopic_ratio_y_tps_co
  
  					if I_want_delta:
  						pl.axhline(y=0.,linewidth=2, color='k')
  						pl.axvline(x=0.,linewidth=2, color='k')
  					else:
  						#	pl.axhline(y=initial_isotopic_ratio_y/initial_isotopic_ratio_y,linewidth=2, color='k')
  						pl.axhline(y=initial_isotopic_ratio_y,linewidth=2, color='k')
  						#	pl.axvline(x=initial_isotopic_ratio_x/initial_isotopic_ratio_x,linewidth=2, color='k')
  						pl.axvline(x=initial_isotopic_ratio_x,linewidth=2, color='k')
  
  					for k in range(0,len(array_to_plot_x)):
  						if I_want_delta:
  							pl.plot(array_to_plot_x[k],array_to_plot_y[k],symbol[k],markersize=10.,linewidth=4.,label=str(metallicity_label[k]))
  							pl.plot(array_to_plot_x_tps[k],array_to_plot_y_tps[k],markersss[k],markersize=10.,linewidth=4.)
  							pl.plot(array_to_plot_x_tps_co[k],array_to_plot_y_tps_co[k],markersss[k],markersize=20.,linewidth=4.)
  
  						else:
  							#		pl.loglog(array_to_plot_x[k],array_to_plot_y[k],symbol[k],markersize=10.,linewidth=2.,label=str(mass_label[k])+'Msun ')
  							pl.loglog(array_to_plot_x[k],array_to_plot_y[k],symbol[k],markersize=10.,linewidth=2.,label=str(mass_label[k])+'Msun ')
  
  					pl.legend(numpoints=1,loc=2,prop={'size':20})
  					size=20
  					pl.xticks(size=size)
  					pl.yticks(size=size)
  					#pl.axis([0.08,0.5,0.1,1.])
  					#pl.grid()
  
  					if I_want_delta:
  						#pl.xlabel(r'$\delta$($^{25}$Mg/$^{24}$Mg)', fontsize=30)
  						#pl.ylabel(r'$\delta$($^{26}$Mg/$^{24}$Mg)', fontsize=30)
  
  						pl.xlabel(r'$\delta$($^{96}$Zr/$^{94}$Zr)', fontsize=20)
  						pl.ylabel(r'$\delta$($^{152}$Gd/$^{154}$Gd)', fontsize=20)
  						#pl.ylabel(r'$\delta$($^{90}$Zr/$^{94}$Zr)', fontsize=30)
  						#pl.ylabel(r'$\delta$($^{91}$Zr/$^{94}$Zr)', fontsize=30)
  						#pl.ylabel(r'$\delta$($^{92}$Zr/$^{94}$Zr)', fontsize=30)
  					
  						#pl.xlabel('d($^{135}$Ba/$^{136}$Ba)', fontsize=20)
  						#pl.ylabel('d($^{138}$Ba/$^{136}$Ba)', fontsize=20)
  						#pl.ylabel('d($^{134}$Ba/$^{136}$Ba)', fontsize=20)
  						#pl.ylabel('d($^{137}$Ba/$^{136}$Ba)', fontsize=20)
  
  						#pl.xlabel(r'$\delta$($^{84}$Sr/$^{86}$Sr)', fontsize=20)
  						#pl.ylabel(r'$\delta$($^{87}$Sr/$^{86}$Sr)', fontsize=20)
  						#pl.ylabel(r'$\delta$($^{88}$Sr/$^{86}$Sr)', fontsize=20)
  					else:
  						pl.xlabel('$^{25}$Mg/$^{24}$Mg', fontsize=20)
  						pl.ylabel('$^{26}$Mg/$^{24}$Mg', fontsize=20)
  
  						#pl.xlabel('$^{96}$Zr/$^{94}$Zr', fontsize=30)
  						#pl.ylabel('$^{152}$Gd/$^{154}$Gd', fontsize=30)
  						#pl.ylabel('$^{90}$Zr/$^{94}$Zr', fontsize=30)
  						#pl.ylabel('$^{91}$Zr/$^{94}$Zr', fontsize=30)
  						#pl.ylabel('$^{92}$Zr/$^{94}$Zr', fontsize=30)
  					
  						#pl.xlabel('d($^{135}$Ba/$^{136}$Ba)', fontsize=20)
  						#pl.ylabel('d($^{138}$Ba/$^{136}$Ba)', fontsize=20)
  						#pl.ylabel('d($^{134}$Ba/$^{136}$Ba)', fontsize=20)
  						#pl.ylabel('d($^{137}$Ba/$^{136}$Ba)', fontsize=20)
  
  						#pl.xlabel('$^{84}$Sr/$^{86}$Sr', fontsize=20)
  						#pl.ylabel('$^{87}$Sr/$^{86}$Sr', fontsize=20)
  						#pl.ylabel('$^{88}$Sr/$^{86}$Sr', fontsize=20)
  						
  					sys.exit()             
  
  			for i in run_mass:
  				dum_s_ini = 0.
  				dum_ls_ini = 0.
  				dum_hs_ini = 0.
  				dum_pb_ini = 0.
  				dum_fe_ini = 0.
  				dum_rb_ini = 0.
  				dum_sr_ini = 0.
  				dum_zr_ini = 0.	
  				dum_pb_ini = float(i.se.get(min(i.se.cycles),'iso_massf','Pb-204'))+float(i.se.get(min(i.se.cycles),'iso_massf','Pb-206'))+float(i.se.get(min(i.se.cycles),'iso_massf','Pb-207'))+float(i.se.get(min(i.se.cycles),'iso_massf','Pb-208'))
  				#
  				dum_fe_ini = float(i.se.get(min(i.se.cycles),'iso_massf','Fe-54'))+float(i.se.get(min(i.se.cycles),'iso_massf','Fe-56'))+float(i.se.get(min(i.se.cycles),'iso_massf','Fe-57'))+float(i.se.get(min(i.se.cycles),'iso_massf','Fe-58'))	
  				dum_rb_ini = float(i.se.get(min(i.se.cycles),'iso_massf','Rb-85'))+float(i.se.get(min(i.se.cycles),'iso_massf','Rb-87'))
  				dum_sr_ini = float(i.se.get(min(i.se.cycles),'iso_massf','Sr-84'))+float(i.se.get(min(i.se.cycles),'iso_massf','Sr-86'))+float(i.se.get(min(i.se.cycles),'iso_massf','Sr-87'))+float(i.se.get(min(i.se.cycles),'iso_massf','Sr-88'))
  				dum_zr_ini = float(i.se.get(min(i.se.cycles),'iso_massf','Zr-90'))+float(i.se.get(min(i.se.cycles),'iso_massf','Zr-91'))+float(i.se.get(min(i.se.cycles),'iso_massf','Zr-92'))+float(i.se.get(min(i.se.cycles),'iso_massf','Zr-94'))+float(i.se.get(min(i.se.cycles),'iso_massf','Zr-96'))
  				dum_y_ini = float(i.se.get(min(i.se.cycles),'iso_massf','Y-89'))
  				dum_zra_ini = float(i.se.get(min(i.se.cycles),'iso_massf','Zr-96'))
  				dum_gda_ini = float(i.se.get(min(i.se.cycles),'iso_massf','Gd-152'))
  				dum_zrb_ini = float(i.se.get(min(i.se.cycles),'iso_massf','Zr-94'))
  				dum_gdb_ini = float(i.se.get(min(i.se.cycles),'iso_massf','Gd-154'))
  
  
  				for j in ls_element:
  					dum_ls_ini = dum_ls_ini + float(i.se.get(min(i.se.cycles),'iso_massf',j))
  				for j in hs_element:
  					dum_hs_ini = dum_hs_ini + float(i.se.get(min(i.se.cycles),'iso_massf',j))
  				for j in s_element:
  					dum_s_ini = dum_s_ini + float(i.se.get(min(i.se.cycles),'iso_massf',j))
  				fe_ini.append(dum_fe_ini)
  				pb_ini.append(dum_pb_ini)
  				ls_ini.append(dum_ls_ini)    #  /float(len(ls_element)))
  				hs_ini.append(dum_hs_ini)     #/float(len(hs_element)))
  				s_ini.append(dum_s_ini)  #/float(len(s_element)))				
  				rb_ini.append(dum_rb_ini)
  				sr_ini.append(dum_sr_ini)
  				zr_ini.append(dum_zr_ini)
  				y_ini.append(dum_y_ini)
  				zra_ini.append(dum_zra_ini)
  				gda_ini.append(dum_gda_ini)
  				zrb_ini.append(dum_zrb_ini)
  				gdb_ini.append(dum_gdb_ini)
  
  
  
  			# sparcity for cycles I am looking at.
  			sparsity = 5000
  			s_fe  = []
  			ls_fe = []
  			hs_fe = []
  			pb_fe = []
  			hs_ls = []
  			rb_sr = []
  			rb_zr = []
  			rb_fe = []
  			sr_y = []
  			sr_zr = []
  			sr_zr_tps=[]
  			sr_y_tps=[]
  			sr_zr_tps_co=[]
  			sr_y_tps_co=[]
  			s_fe_tps  = []
  			ls_fe_tps = []
  			hs_fe_tps = []
  			pb_fe_tps = []
  			hs_ls_tps = []
  			rb_sr_tps = []
  			rb_zr_tps = []
  			rb_fe_tps = []
  			s_fe_tps_co  = []
  			ls_fe_tps_co = []
  			hs_fe_tps_co = []
  			pb_fe_tps_co = []
  			hs_ls_tps_co = []
  			rb_sr_tps_co = []
  			rb_zr_tps_co = []
  			rb_fe_tps_co = []
  			gd=[]
  			zr=[]
  			gd_tps=[]
  			zr_tps=[]
  			gd_tps_co=[]
  			zr_tps_co=[]
  			co=[]
  			k = 0
  
  			for i in run_mass:
  
  				jjjj=0
  				dum_s_fe = []
  				dum_ls_fe = []
  				dum_hs_fe = []
  				dum_pb_fe = []
  				dum_hs_ls = []
  				dum_rb_sr = []
  				dum_rb_zr = []
  				dum_rb_fe = []
  				dum_sr_y = []
  				dum_sr_zr = []
  				dum_sr_zr_tps=[]
  				dum_sr_y_tps=[]
  				dum_sr_zr_tps_co=[]
  				dum_sr_y_tps_co=[]
  				dum_s_fe_tps = []
  				dum_ls_fe_tps = []
  				dum_hs_fe_tps = []
  				dum_pb_fe_tps = []
  				dum_hs_ls_tps = []
  				dum_rb_sr_tps = []
  				dum_rb_fe_tps = []
  				dum_rb_zr_tps = []
  				dum_sr_y_tps = []
  				dum_sr_zr_tps = []
  				dum_s_fe_tps_co = []
  				dum_ls_fe_tps_co = []
  				dum_hs_fe_tps_co = []
  				dum_pb_fe_tps_co = []
  				dum_hs_ls_tps_co = []
  				dum_rb_sr_tps_co = []
  				dum_rb_zr_tps_co = []
  				dum_rb_fe_tps_co = []
  				dum_rb_sr_tps_co = []
  				dum_sr_y_tps_co = []
  				dum_sr_zr_tps_co = []
  				dum_gd=[]
  				dum_zrr=[]
  				dum_gd_tps=[]
  				dum_zr_tps=[]
  				dum_gd_tps_co=[]
  				dum_zr_tps_co=[]
  				dum_co=[]
  
  				for j in i.se.cycles[0::sparsity]:
  					dum_s = 0.
  					dum_ls = 0.
  					dum_hs = 0.
  					dum_fe = 0.
  					dum_pb = 0.
  					dum_rb = 0.
  					dum_sr = 0.
  					dum_zr = 0.
  					dum_y = 0.
  					dum_c = 0.
  					dum_o = 0.
  					dum_gda = 0.
  					dum_gdb = 0.
  					dum_zra = 0.
  					dum_zrb = 0.
  
  					print j
  					dum_fe = float(i.se.get(j,'iso_massf','Fe-54'))+float(i.se.get(j,'iso_massf','Fe-56'))+float(i.se.get(j,'iso_massf','Fe-57'))+float(i.se.get(j,'iso_massf','Fe-58')) 
  					dum_pb = float(i.se.get(j,'iso_massf','Pb-204'))+float(i.se.get(j,'iso_massf','Pb-206'))+float(i.se.get(j,'iso_massf','Pb-207'))+float(i.se.get(j,'iso_massf','Pb-208'))
  					dum_rb = float(i.se.get(j,'iso_massf','Rb-85'))+float(i.se.get(j,'iso_massf','Rb-87'))
  					dum_sr = float(i.se.get(j,'iso_massf','Sr-84'))+float(i.se.get(j,'iso_massf','Sr-86'))+float(i.se.get(j,'iso_massf','Sr-87'))+float(i.se.get(j,'iso_massf','Sr-88'))
  					dum_zr = float(i.se.get(j,'iso_massf','Zr-90'))+float(i.se.get(j,'iso_massf','Zr-91'))+float(i.se.get(j,'iso_massf','Zr-92'))+float(i.se.get(j,'iso_massf','Zr-94'))+float(i.se.get(j,'iso_massf','Zr-96'))
  					dum_y = float(i.se.get(j,'iso_massf','Y-89'))
  					dum_c = float(i.se.get(j,'iso_massf','C-12'))
  					dum_o = float(i.se.get(j,'iso_massf','O-16'))
  					dum_gda = float(i.se.get(j,'iso_massf','Gd-152'))
  					dum_gdb = float(i.se.get(j,'iso_massf','Gd-154'))
  					dum_zra = float(i.se.get(j,'iso_massf','Zr-96'))
  					dum_zrb = float(i.se.get(j,'iso_massf','Zr-94'))
  
  					dum_c=(float((i.se.get((int(j)),'iso_massf','C-12')))+float((i.se.get(j,'iso_massf','C-13'))))
  					dum_o=(float((i.se.get((int(j)),'iso_massf','O-16')))+float((i.se.get(j,'iso_massf','O-17'))))
  
  					dum_c12=(float((i.se.get((int(j)),'iso_massf','C-12'))))
  					dum_c13=(float((i.se.get((int(j)),'iso_massf','C-13'))))
  
  					dum_co.append(((dum_c/dum_o)*(16./12.)))
  
  					for jj in ls_element:
  						dum_ls = dum_ls + float(i.se.get(j,'iso_massf',jj)) #/float(len(ls_element)))
  					for jj in hs_element:
  						dum_hs = dum_hs + float(i.se.get(j,'iso_massf',jj))  #/float(len(hs_element)))
  					for jj in s_element:
  						dum_s = dum_s + float(i.se.get(j,'iso_massf',jj)) #/float(len(s_element)))
  					
  					dum_s_fe.append(log10((dum_s/dum_fe)/(s_ini[k]/fe_ini[k])))
  					dum_ls_fe.append(log10((dum_ls/dum_fe)/(ls_ini[k]/fe_ini[k])))
  					dum_hs_fe.append(log10((dum_hs/dum_fe)/(hs_ini[k]/fe_ini[k])))
  					dum_pb_fe.append(log10((dum_pb/dum_fe)/(pb_ini[k]/fe_ini[k])))
  					dum_hs_ls.append(log10((dum_hs/dum_ls)/(hs_ini[k]/ls_ini[k])))
  					dum_rb_sr.append(log10((dum_rb/dum_sr)/(rb_ini[k]/sr_ini[k])))
  					dum_rb_zr.append(log10((dum_rb/dum_zr)/(rb_ini[k]/zr_ini[k])))
  					dum_sr_y.append(log10((dum_sr/dum_y)/(sr_ini[k]/y_ini[k])))
  					dum_sr_zr.append(log10((dum_sr/dum_zr)/(sr_ini[k]/zr_ini[k])))
  					dum_gd.append((dum_gda/dum_gdb)/(gda_ini[k]/gdb_ini[k]))
  					dum_zrr.append((dum_zra/dum_zrb)/(zra_ini[k]/zrb_ini[k]))
  			#                dum_rb_fe.append(log10((dum_rb/dum_fe)/(rb_ini[k]/fe_ini[k]))-0.2) ## correction for Lambert 1995 data
  					dum_rb_fe.append(log10((dum_rb/dum_fe)/(rb_ini[k]/fe_ini[k]))-0.068) ## correction for Zamora 2009 data
  
  					if (len(dum_co)>1):
  						if (dum_co[len(dum_co)-1]>(dum_co[len(dum_co)-2]+0.02)):
  							dum_sr_zr_tps.append(log10((dum_sr/dum_zr)/(sr_ini[k]/zr_ini[k])))
  							dum_sr_y_tps.append(log10((dum_sr/dum_y)/(sr_ini[k]/y_ini[k])))
  							dum_s_fe_tps.append(log10((dum_s/dum_fe)/(s_ini[k]/fe_ini[k])))
  							dum_ls_fe_tps.append(log10((dum_ls/dum_fe)/(ls_ini[k]/fe_ini[k])))
  							dum_hs_fe_tps.append(log10((dum_hs/dum_fe)/(hs_ini[k]/fe_ini[k])))
  							dum_pb_fe_tps.append(log10((dum_pb/dum_fe)/(pb_ini[k]/fe_ini[k])))
  							dum_hs_ls_tps.append(log10((dum_hs/dum_ls)/(hs_ini[k]/ls_ini[k])))
  							dum_rb_sr_tps.append(log10((dum_rb/dum_sr)/(rb_ini[k]/sr_ini[k])))
  							dum_rb_zr_tps.append(log10((dum_rb/dum_zr)/(rb_ini[k]/zr_ini[k])))
  			#                                dum_rb_fe_tps.append(log10((dum_rb/dum_fe)/(rb_ini[k]/fe_ini[k]))-0.2) ## correction for Lambert 1995 data
  							dum_rb_fe_tps.append(log10((dum_rb/dum_fe)/(rb_ini[k]/fe_ini[k]))-0.068) ## correction for Zamora 2009 data
  							dum_gd_tps.append((((dum_gda/dum_gdb)/(gda_ini[k]/gdb_ini[k]))-1)*1000.)
  							dum_zr_tps.append((((dum_zra/dum_zrb)/(zra_ini[k]/zrb_ini[k]))-1)*1000.)
  
  							if (dum_co[len(dum_co)-1]>1.):
  								dum_sr_zr_tps_co.append(log10((dum_sr/dum_zr)/(sr_ini[k]/zr_ini[k])))
  								dum_sr_y_tps_co.append(log10((dum_sr/dum_y)/(sr_ini[k]/y_ini[k])))
  								dum_s_fe_tps_co.append(log10((dum_s/dum_fe)/(s_ini[k]/fe_ini[k])))
  								dum_ls_fe_tps_co.append(log10((dum_ls/dum_fe)/(ls_ini[k]/fe_ini[k])))
  								dum_hs_fe_tps_co.append(log10((dum_hs/dum_fe)/(hs_ini[k]/fe_ini[k])))
  								dum_pb_fe_tps_co.append(log10((dum_pb/dum_fe)/(pb_ini[k]/fe_ini[k])))
  								dum_hs_ls_tps_co.append(log10((dum_hs/dum_ls)/(hs_ini[k]/ls_ini[k])))
  								dum_rb_sr_tps_co.append(log10((dum_rb/dum_sr)/(rb_ini[k]/sr_ini[k])))
  								dum_rb_zr_tps_co.append(log10((dum_rb/dum_zr)/(rb_ini[k]/zr_ini[k])))
  			###                                        dum_rb_fe_tps_co.append(log10((dum_rb/dum_fe)/(rb_ini[k]/fe_ini[k]))-0.2)
  								dum_rb_zr_tps_co.append(log10((dum_rb/dum_zr)/(rb_ini[k]/zr_ini[k])))
  			#                                        dum_rb_fe_tps_co.append(log10((dum_rb/dum_fe)/(rb_ini[k]/fe_ini[k]))-0.2) ## correction for Lambert 1995 data
  								dum_rb_fe_tps_co.append(log10((dum_rb/dum_fe)/(rb_ini[k]/fe_ini[k]))-0.068) ## correction for Zamora 2009 data
  								dum_gd_tps_co.append((((dum_gda/dum_gdb)/(gda_ini[k]/gdb_ini[k]))-1)*1000.)
  								dum_zr_tps_co.append((((dum_zra/dum_zrb)/(zra_ini[k]/zrb_ini[k]))-1)*1000.)
  							jjjj=jjjj+1
  					
  					#
  				s_fe.append(dum_s_fe)
  				ls_fe.append(dum_ls_fe)
  				hs_fe.append(dum_hs_fe)
  				pb_fe.append(dum_pb_fe)
  				hs_ls.append(dum_hs_ls)
  				rb_sr.append(dum_rb_sr)
  				rb_zr.append(dum_rb_zr)
  				rb_fe.append(dum_rb_fe)
  				sr_y.append(dum_sr_y)
  				sr_zr.append(dum_sr_zr)
  				sr_zr_tps.append(dum_sr_zr_tps)
  				sr_y_tps.append(dum_sr_y_tps)
  				sr_zr_tps_co.append(dum_sr_zr_tps_co)
  				sr_y_tps_co.append(dum_sr_y_tps_co)
  				s_fe_tps.append(dum_s_fe_tps)
  				ls_fe_tps.append(dum_ls_fe_tps)
  				hs_fe_tps.append(dum_hs_fe_tps)
  				pb_fe_tps.append(dum_pb_fe_tps)
  				hs_ls_tps.append(dum_hs_ls_tps)
  				rb_sr_tps.append(dum_rb_sr_tps)
  				rb_zr_tps.append(dum_rb_zr_tps)
  				rb_fe_tps.append(dum_rb_fe_tps)
  				s_fe_tps_co.append(dum_s_fe_tps_co)
  				ls_fe_tps_co.append(dum_ls_fe_tps_co)
  				hs_fe_tps_co.append(dum_hs_fe_tps_co)
  				pb_fe_tps_co.append(dum_pb_fe_tps_co)
  				hs_ls_tps_co.append(dum_hs_ls_tps_co)
  				rb_sr_tps_co.append(dum_rb_sr_tps_co)
  				rb_zr_tps_co.append(dum_rb_zr_tps_co)
  				rb_fe_tps_co.append(dum_rb_fe_tps_co)
  				gd.append(dum_gd)
  				zr.append(dum_zrr)
  				gd_tps.append(dum_gd_tps)
  				zr_tps.append(dum_zr_tps)
  				gd_tps_co.append(dum_gd_tps_co)
  				zr_tps_co.append(dum_zr_tps_co)
  				co.append(dum_co)
  				k = k+1
  
  			if plotting:
  
  				mass_label =[]
  				for i in run_mass:
  					mini=float(i.se.get('mini'))
  					zini=float(i.se.get('zini'))
  					label=str(mini)+'$M_{\odot}$, Z='+str(zini)
  					mass_label.append(label)
  
  				metallicity_label=['3 Msun, Z=0.02','3 Msun, Z=0.02, no mol diff','3 Msun, Z=0.01','M3 set1.2']
  
  				fig = plt.figure(0)            # Figure object
  				ax = fig.add_subplot(1,1,1)     # Axes object: one row, one column, first plot (one plot!)
  
  				mpl.rcParams['xtick.major.size'] = 20
  				#mpl.rcParams['xtick.major.width'] = 4
  				mpl.rcParams['xtick.minor.size'] = 10
  				#mpl.rcParams['xtick.minor.width'] = 2
  				
  				mpl.rcParams['ytick.major.size'] = 20
  				#mpl.rcParams['ytick.major.width'] = 4
  				mpl.rcParams['ytick.minor.size'] = 10
  				#mpl.rcParams['ytick.minor.width'] = 2
  				
  				if hsfe_hsls:
  					array_to_plot_x = hs_ls
  					array_to_plot_x_tps = hs_ls_tps
  					array_to_plot_x_co = hs_ls_tps_co
  					
  					array_to_plot_y = hs_fe
  					array_to_plot_y_tps = hs_fe_tps
  					array_to_plot_y_co = hs_fe_tps_co
  
  					pl.axhline(y=0,linewidth=2, color='k')
  					pl.axvline(x=0,linewidth=2, color='k')
  					for k in range(0,len(array_to_plot_x)):
  						#if k > 0:
  						pl.plot(array_to_plot_x[k],array_to_plot_y[k],symbol[k],markersize=10.,linewidth=3.,color=color1)
  						#pl.plot(array_to_plot_x[k],array_to_plot_y[k],symbol[k],markersize=10.,linewidth=3.,label=str(metallicity_label[k]))
  						#if k >= 1:
  						pl.plot(array_to_plot_x_tps[k],array_to_plot_y_tps[k],markers[k],markersize=12.,color=color1)
  						pl.plot(array_to_plot_x_co[k],array_to_plot_y_co[k],markers[k],markersize=20.,color=color1,label=str(mass_label[k]))
  
  					pl.legend(numpoints=1,loc=2,prop={'size':20})
  
  					pl.xlabel('[hs/ls]', fontsize=20)
  					pl.ylabel('[hs/Fe]', fontsize=20)
  					y_min = -0.1
  					y_max = 1.0
  					pl.ylim(y_min,y_max)
  					pl.xlim(-0.1,0.5)
  
  					size=20
  					pl.xticks(size=size)
  					pl.yticks(size=size)
  
  					ax = pl.gca()
  
  					for line in ax.xaxis.get_ticklines():
  						line.set_markersize(20)
  						line.set_markeredgewidth(3)
  
  					for line in ax.yaxis.get_ticklines():
  						line.set_markersize(20)
  						line.set_markeredgewidth(3)
  
  				if lsfe_hsls:
  
  					array_to_plot_x = hs_ls
  					array_to_plot_x_tps = hs_ls_tps
  					array_to_plot_x_co = hs_ls_tps_co
  					
  					array_to_plot_y = ls_fe
  					array_to_plot_y_tps = ls_fe_tps
  					array_to_plot_y_co = ls_fe_tps_co
  
  					pl.axhline(y=0,linewidth=2, color='k')
  					pl.axvline(x=0,linewidth=2, color='k')
  					for k in range(0,len(array_to_plot_x)):
  						#                if k > 0:
  						pl.plot(array_to_plot_x[k],array_to_plot_y[k],symbol[k],markersize=10.,linewidth=4.,label=str(mass_label[k]))
  						pl.plot(array_to_plot_x_tps[k],array_to_plot_y_tps[k],markers[k],markersize=12.)
  						pl.plot(array_to_plot_x_co[k],array_to_plot_y_co[k],markers[k],markersize=20.)
  
  					pl.legend(numpoints=1,loc='upper right',prop={'size':20})
  
  					pl.xlabel('[hs/ls]', fontsize=20)
  					pl.ylabel('[ls/Fe]', fontsize=20)
  					y_min = -0.1
  					y_max = 0.8
  					pl.ylim(y_min,y_max)
  					pl.xlim(-0.2,0.5)
  
  					size=20
  					pl.xticks(size=size)
  					pl.yticks(size=size)
  
  					ax = pl.gca()
  
  					for line in ax.xaxis.get_ticklines():
  						line.set_markersize(20)
  						line.set_markeredgewidth(3)
  
  					for line in ax.yaxis.get_ticklines():
  						line.set_markersize(20)
  						line.set_markeredgewidth(3)
  
  				if rbsr_hsls:
  					array_to_plot_x = hs_ls
  					array_to_plot_y = rb_sr
  
  					pl.axhline(y=0,linewidth=2, color='k')
  					pl.axvline(x=0,linewidth=2, color='k')
  					for k in range(0,len(array_to_plot_x)):
  						#                if k > 0:
  						if k >= 1:
  							pl.plot(array_to_plot_x[k],array_to_plot_y[k],symbol[k],markersize=10.,linewidth=4.,label=str(mass_label[k]))
  						else:
  							pl.plot(array_to_plot_x[k],array_to_plot_y[k],symbol[k],markersize=10.,linewidth=4.,label=str(mass_label[k]))
  
  					pl.legend(numpoints=1,loc='upper right',prop={'size':20})
  
  					pl.xlabel('[hs/ls]', fontsize=20)
  					pl.ylabel('[Rb/Sr]', fontsize=20)
  					y_min = -0.4
  					y_max = 1.5
  					x = ax.get_position()
  
  					pl.ylim(y_min,y_max)
  					pl.xlim(-0.7,0.4)
  
  					size=30
  					pl.xticks(size=size)
  					pl.yticks(size=size)
  
  					ax = pl.gca()
  
  					for line in ax.xaxis.get_ticklines():
  						line.set_markersize(20)
  						line.set_markeredgewidth(3)
  					for line in ax.yaxis.get_ticklines():
  						line.set_markersize(20)
  						line.set_markeredgewidth(3)
  
  				if rbzr_hsls:
  					array_to_plot_x = hs_ls
  					array_to_plot_y = rb_zr
  
  					pl.axhline(y=0,linewidth=2, color='k')
  					pl.axvline(x=0,linewidth=2, color='k')
  					for k in range(0,len(array_to_plot_x)):
  			#                if k > 0:
  						if k >= 1:
  							pl.plot(array_to_plot_x[k],array_to_plot_y[k],symbol[k],markersize=10.,linewidth=2.,label=str(mass_label[k]))
  							pl.plot(hs_ls_tps[k],rb_zr_tps[k],markers[k],markersize=12.)
  							pl.plot(hs_ls_tps_co[k],rb_zr_tps_co[k],markers[k],markersize=20.)
  						else:
  							pl.plot(array_to_plot_x[k],array_to_plot_y[k],symbol[k],markersize=10.,linewidth=2.,label=str(mass_label[k]))
  							pl.plot(hs_ls_tps[k],rb_zr_tps[k],markers[k],markersize=12.)
  							pl.plot(hs_ls_tps_co[k],rb_zr_tps_co[k],markers[k],markersize=20.)
  
  					pl.legend(numpoints=1,loc='upper right',prop={'size':20})
  
  					pl.xlabel('[hs/ls]', fontsize=20)
  					pl.ylabel('[Rb/Zr]', fontsize=20)
  					y_min = -0.4
  					y_max = 0.4
  					pl.ylim(y_min,y_max)
  					pl.xlim(-0.5,0.5)
  
  					size=20
  					pl.xticks(size=size)
  					pl.yticks(size=size)
  
  					ax = pl.gca()
  
  					for line in ax.xaxis.get_ticklines():
  						line.set_markersize(20)
  						line.set_markeredgewidth(3)
  
  					for line in ax.yaxis.get_ticklines():
  						line.set_markersize(20)
  						line.set_markeredgewidth(3)
  
  				if rbfe_sfe:
  
  					array_to_plot_x = s_fe
  					array_to_plot_y = rb_fe
  
  					pl.axhline(y=0,linewidth=3, color='k')
  					pl.axvline(x=0,linewidth=3, color='k')
  					for k in range(0,len(array_to_plot_x)):
  			#                if k > 0:
  						if k >= 1:
  			#                        pl.plot(array_to_plot_x[k],array_to_plot_y[k],symbol[k],markersize=10.,linewidth=2.,label=str(mass_label[k])+'Msun '+str(metallicity_label[k]))
  							pl.plot(array_to_plot_x[k],array_to_plot_y[k],symbol[k],markersize=10.,linewidth=2.,label=str(metallicity_label[k]))
  							pl.plot(s_fe_tps[k],rb_fe_tps[k],markers[k],markersize=12.)
  							pl.plot(s_fe_tps_co[k],rb_fe_tps_co[k],markers[k],markersize=20.)
  						else:
  			#                        pl.plot(array_to_plot_x[k],array_to_plot_y[k],symbol[k],markersize=10.,linewidth=2.,label=str(mass_label[k])+'Msun '+str(metallicity_label[k]))
  							pl.plot(array_to_plot_x[k],array_to_plot_y[k],symbol[k],markersize=10.,linewidth=2.,label=str(metallicity_label[k]))
  							pl.plot(s_fe_tps[k],rb_fe_tps[k],markers[k],markersize=12.)
  							pl.plot(s_fe_tps_co[k],rb_fe_tps_co[k],markers[k],markersize=20.)
  
  					pl.legend(numpoints=1,loc=4,prop={'size':20})
  
  					pl.xlabel('[s/Fe]', fontsize=20)
  					pl.ylabel('[Rb/Fe]', fontsize=20)
  					y_min = -0.4
  					y_max = 0.4
  					pl.ylim(y_min,y_max)
  					pl.xlim(-0.1,2.0)
  
  					size=20
  					pl.xticks(size=size)
  					pl.yticks(size=size)
  
  					ax = pl.gca()
  
  					for line in ax.xaxis.get_ticklines():
  						line.set_markersize(25)
  						line.set_markeredgewidth(3)
  
  					for line in ax.yaxis.get_ticklines():
  						line.set_markersize(25)
  						line.set_markeredgewidth(3)
  
  				if sry_srzr:
  					array_to_plot_x = sr_zr
  					array_to_plot_y = sr_y
  
  					pl.axhline(y=0,linewidth=2, color='k')
  					pl.axvline(x=0,linewidth=2, color='k')
  					for k in range(0,len(array_to_plot_x)):
  						if k >= 1:
  							#                if k > 0:
  							pl.plot(array_to_plot_x[k],array_to_plot_y[k],symbol[k],markersize=10.,linewidth=3.5,label=str(metallicity_label[k]))
  							pl.plot(sr_zr_tps[k],sr_y_tps[k],markers[k],markersize=12.)
  							pl.plot(sr_zr_tps_co[k],sr_y_tps_co[k],markers[k],markersize=20.)
  						else:
  							pl.plot(array_to_plot_x[k],array_to_plot_y[k],symbol[k],markersize=10.,linewidth=3.5,label=str(metallicity_label[k]))
  							pl.plot(sr_zr_tps[k],sr_y_tps[k],markers[k],markersize=12.,label='C/O < 1')
  							pl.plot(sr_zr_tps_co[k],sr_y_tps_co[k],markers[k],markersize=20.,label='C/O > 1')
  
  					pl.xlabel('[Sr/Zr]', fontsize=20)
  					pl.ylabel('[Sr/Y]', fontsize=20)
  					y_min = -0.1
  					y_max = 0.1
  					pl.ylim(y_min,y_max)
  					pl.xlim(-0.3,0.3)
  					
  					size=20
  					pl.xticks(size=size)
  					pl.yticks(size=size)
  
  					ax = pl.gca()
  
  					for line in ax.xaxis.get_ticklines():
  						line.set_markersize(25)
  						line.set_markeredgewidth(3)
  
  					for line in ax.yaxis.get_ticklines():
  						line.set_markersize(25)
  						line.set_markeredgewidth(3)
  
  					box = ax.get_position()
  
  					ax.set_position([box.x0, box.y0, box.width * 0.9, box.height])
  
  					pl.legend(loc=2, borderaxespad=1.0,prop={'size':25})
  
  			box = ax.get_position()