For paper1 marco routine:
			 Numbers of remnant mass shell masses, exists also in mesa_set!
172747
def final_bottom_envelope_set1(self):
		'''
			
		'''
		inim=[]
		remnm=[]
		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)-1],'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)-1],'iso_massf','H-1')

			for j in range(len(mass_dummy)):
        			if h_dummy[j] > 0.05:
                			bottom_of_envelope = mass_dummy[j]
                			break
			inim.append(m1p65_last.get("mini"))
			remnm.append(bottom_of_envelope)
		print "M_initial | M_remn/bottom of envelope"
		for i in range(len(inim)):
			print inim[i],"|",remnm[i]
280502
def primary_mass(mass1, mass2):
    """"""
    mass1, mass2, input_is_array = ensurearray(mass1, mass2)
    mp = copy.copy(mass1)
    mask = mass1 < mass2
    mp[mask] = mass2[mask]
    return formatreturn(mp, input_is_array)
923300
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']):
		'''
			  
		'''

		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()
233117
def generate(self, M, age=9.6, feh=0.0,
                 ichrone='mist', n=1e4, bands=None, **kwargs):
        
        ichrone = get_ichrone(ichrone, bands=bands)
        if np.size(M) > 1:
            n = np.size(M)
        else:
            n = int(n)
        M2 = M * self.q_fn(n, qmin=np.maximum(self.qmin,self.minmass/M))

        P = self.P_fn(n)
        ecc = self.ecc_fn(n,P)

        mass = np.ascontiguousarray(np.ones(n)*M)
        mass2 = np.ascontiguousarray(M2)
        age = np.ascontiguousarray(age)
        feh = np.ascontiguousarray(feh)
        pri = ichrone(mass, age, feh, return_df=True, bands=bands)
        sec = ichrone(mass2, age, feh, return_df=True, bands=bands)

        BinaryPopulation.__init__(self, primary=pri, secondary=sec,
                                  period=P, ecc=ecc, **kwargs)
        return self
1350487
def find_first_TP(self):

        

        star_mass         = self.get('star_mass')
        he_lumi           = self.get('log_LHe')
        h_lumi            = self.get('log_LH')
        mx2_bot           = self.get('mx2_bot')*star_mass
        try:
           h1_boundary_mass  = self.get('h1_boundary_mass')
           he4_boundary_mass = self.get('he4_boundary_mass')
        except:
           try:
                h1_boundary_mass  = self.get('he_core_mass')
                he4_boundary_mass = self.get('c_core_mass')
           except:
                pass


        TP_bot=np.array(self.get('conv_mx2_bot'))*np.array(self.get('star_mass'))
        TP_top=np.array(self.get('conv_mx2_top'))*np.array(self.get('star_mass'))
        lum_array=[]
        activate=False
        models=[]
        pdcz_size=[]
        for i in range(len(h1_boundary_mass)):
            if (h1_boundary_mass[i]-he4_boundary_mass[i] <0.2) and (he4_boundary_mass[i]>0.2):
                if (mx2_bot[i]>he4_boundary_mass[i]) and (he_lumi[i]>h_lumi[i]):
                        if TP_top[i]>he4_boundary_mass[i]:
                                pdcz_size.append(TP_top[i]-TP_bot[i])
                                activate=True
                                lum_array.append(he_lumi[i])
                                models.append(i)
                                #print(TP_bot[i],TP_top[i])
                if (activate == True) and (he_lumi[i]<h_lumi[i]):
                        #if fake tp
                        if max(pdcz_size)<1e-5:
                                active=False
                                lum_array=[]
                                models=[]
                                print('fake tp')
                        else:
                                break
        t0_model = models[np.argmax(lum_array)]
        return t0_model
822430
def _cont_norm_running_quantile_regions_mp(wl, fluxes, ivars, q, delta_lambda,
                                           ranges, n_proc=2, verbose=False):
    
    print("contnorm.py: continuum norm using running quantile")
    print("Taking spectra in %s chunks" % len(ranges))
    # nstars = fluxes.shape[0]
    nchunks = len(ranges)
    norm_fluxes = np.zeros(fluxes.shape)
    norm_ivars = np.zeros(ivars.shape)
    for i in xrange(nchunks):
        chunk = ranges[i, :]
        start = chunk[0]
        stop = chunk[1]
        if verbose:
            print('@Bo Zhang: Going to normalize Chunk [%d/%d], pixel:[%d, %d] ...'
                  % (i+1, nchunks, start, stop))
        output = _cont_norm_running_quantile_mp(
            wl[start:stop], fluxes[:, start:stop],
            ivars[:, start:stop], q, delta_lambda,
            n_proc=n_proc, verbose=verbose)
        norm_fluxes[:, start:stop] = output[0]
        norm_ivars[:, start:stop] = output[1]
    return norm_fluxes, norm_ivars
1363334
def electron_shell_str(shell, shellidx=None):
    '''
    '''
    am = shell['angular_momentum']
    amchar = lut.amint_to_char(am)
    amchar = amchar.upper()

    shellidx_str = ''
    if shellidx is not None:
        shellidx_str = 'Index {} '.format(shellidx)

    exponents = shell['exponents']
    coefficients = shell['coefficients']
    ncol = len(coefficients) + 1

    point_places = [8 * i + 15 * (i - 1) for i in range(1, ncol + 1)]
    s = "Shell: {}Region: {}: AM: {}\n".format(shellidx_str, shell['region'], amchar)
    s += "Function: {}\n".format(shell['function_type'])
    s += write_matrix([exponents, *coefficients], point_places)
    return s
828954
def mass_3d(self, R, Rs, rho0, r_trunc):
        

        x = R * Rs ** -1

        func = (r_trunc ** 2 * (-2 * x * (1 + r_trunc ** 2) + 4 * (1 + x) * r_trunc * np.arctan(x / r_trunc) -
                                2 * (1 + x) * (-1 + r_trunc ** 2) * np.log(Rs) + 2 * (1 + x) * (-1 + r_trunc ** 2) * np.log(Rs * (1 + x)) +
                                2 * (1 + x) * (-1 + r_trunc ** 2) * np.log(Rs * r_trunc) -
                                (1 + x) * (-1 + r_trunc ** 2) * np.log(Rs ** 2 * (x ** 2 + r_trunc ** 2)))) / (2. * (1 + x) * (1 + r_trunc ** 2) ** 2)

        m_3d = 4*np.pi*Rs ** 3 * rho0 * func
        return m_3d
230692
def mass_3d(self, R, Rs, rho0, r_core):
        
        b = r_core * Rs ** -1
        x = R * Rs ** -1

        M_0 = 4 * np.pi * Rs**3 * rho0

        return M_0 * (x * (1+x) ** -1 * (-1+b) ** -1 + (-1+b) ** -2 *
                      ((2*b-1)*np.log(1/(1+x)) + b **2 * np.log(x / b + 1)))
1030258
def _create_waveforms(self):
        

        # find eccentricity and semi major axis over time until e=0.
        e_vals, a_vals, t_vals = self._t_of_e(a0=self.a0, f0=self.f0,
                                              t_start=self.t_start, ef=None,
                                              t_obs=self.t_obs)

        f_mrg = 0.02/(self.m1 + self.m2)
        a_mrg = ((self.m1+self.m2)/f_mrg**2)**(1/3)

        # limit highest frequency to ISCO even though this is not innermost orbit for eccentric
        # binaries
        # find where binary goes farther than observation time or merger frequency limit.
        a_ind_start = np.asarray([np.where(a_vals[i] > a_mrg[i])[0][0] for i in range(len(a_vals))])
        t_ind_start = np.asarray([np.where(t_vals[i] < self.t_obs[i])[0][0]
                                 for i in range(len(t_vals))])

        ind_start = (a_ind_start*(a_ind_start >= t_ind_start)
                     + t_ind_start*(a_ind_start < t_ind_start))

        self.ef = np.asarray([e_vals[i][ind] for i, ind in enumerate(ind_start)])

        # higher resolution over the eccentricities seen during observation
        self.e_vals, self.a_vals, self.t_vals = self._t_of_e(a0=a_vals[:, -1],
                                                             ef=self.ef,
                                                             t_obs=self.t_obs)

        self.freqs_orb = np.sqrt((self.m1[:, np.newaxis]+self.m2[:, np.newaxis])/self.a_vals**3)

        # tile for efficient calculation across modes.
        for attr in ['e_vals', 'a_vals', 't_vals', 'freqs_orb']:
            arr = getattr(self, attr)
            new_arr = (np.flip(
                       np.tile(arr, self.n_max).reshape(len(arr)*self.n_max, len(arr[0])), -1))
            setattr(self, attr, new_arr)

        for attr in ['m1', 'm2', 'z', 'dist']:
            arr = getattr(self, attr)
            new_arr = np.repeat(arr, self.n_max)[:, np.newaxis]
            setattr(self, attr, new_arr)

        # setup modes
        self.n = np.tile(np.arange(1, self.n_max + 1), self.length)[:, np.newaxis]

        self._hcn_func()

        # reshape hc
        self.hc = self.hc.reshape(self.length, self.n_max, self.hc.shape[-1])
        self.freqs = np.reshape(self.n*self.freqs_orb/(1+self.z)
                                * ct.c,
                                (self.length, self.n_max, self.freqs_orb.shape[-1]))

        self.hc, self.freqs = np.squeeze(self.hc), np.squeeze(self.freqs)
        return
789568
def _cont_norm_running_quantile_mp(wl, fluxes, ivars, q, delta_lambda,
                                   n_proc=2, verbose=False):
    
    nStar = fluxes.shape[0]

    # start mp.Pool
    mp_results = []
    pool = mp.Pool(processes=n_proc)
    for i in xrange(nStar):
        mp_results.append(pool.apply_async(\
            _find_cont_running_quantile,
            (wl, fluxes[i, :].reshape((1, -1)), ivars[i, :].reshape((1, -1)),
             q, delta_lambda), {'verbose': False}))
        if verbose:
            print('@Bo Zhang: continuum normalizing star [%d/%d] ...'\
                  % (i + 1, nStar))
    # close mp.Pool
    pool.close()
    pool.join()

    # reshape results --> cont
    cont = np.zeros_like(fluxes)
    for i in xrange(nStar):
        cont[i, :] = mp_results[i].get() #.flatten()
    norm_fluxes = np.ones(fluxes.shape)
    norm_fluxes[cont!=0] = fluxes[cont!=0] / cont[cont!=0]
    norm_ivars = cont**2 * ivars

    print('@Bo Zhang: continuum normalization finished!')
    return norm_fluxes, norm_ivars
724426
def Fierz_to_Bern_chrom(C, dd, parameters):
    """"""
    e = sqrt(4 * pi * parameters['alpha_e'])
    gs = sqrt(4 * pi * parameters['alpha_s'])
    if dd == 'sb' or dd == 'db':
        mq = parameters['m_b']
    elif dd == 'ds':
        mq = parameters['m_s']
    else:
        KeyError("Not sure what to do with quark mass for flavour {}".format(dd))
    return {
        '7gamma' + dd : gs**2 / e / mq * C['F7gamma' + dd ],
        '8g' + dd : gs / mq * C['F8g' + dd ],
        '7pgamma' + dd : gs**2 / e /mq * C['F7pgamma' + dd],
        '8pg' + dd : gs / mq * C['F8pg' + dd]
            }
672117
def mass_2d(self, r, rho0, Ra, Rs):
        
        Ra, Rs = self._sort_ra_rs(Ra, Rs)
        sigma0 = self.rho2sigma(rho0, Ra, Rs)
        m_2d = 2 * np.pi * sigma0 * Ra * Rs / (Rs - Ra) * (np.sqrt(Ra ** 2 + r ** 2) - Ra - np.sqrt(Rs ** 2 + r ** 2) + Rs)
        return m_2d
782098
def mass_2d(self, R, Rs, rho0):
        
        x = R/Rs
        gx = self.g_(x)
        m_2d = 4*rho0*Rs*R**2*gx/x**2 * np.pi
        return m_2d
839179
def meval(self, f, *args, **kw):
        
        # !! PART OF ORIGINAL DOCSTRING INCOMPATIBLE WITH CLASS INTERFACE !!
        # Example
        # -------
        # We set up multiple stationary models, one for a reference (ground)
        # state, and two for biased states, and group them in a
        # MultiStationaryModel.
        # >>> from pyemma.thermo import StationaryModel, MEMM
        # >>> m_1 = StationaryModel(f=[1.0, 0], label='biased 1')
        # >>> m_2 = StationaryModel(f=[2.0, 0], label='biased 2')
        # >>> m_mult = MEMM([m_1, m_2], [0, 0], label='unbiased')
        # Compute the stationary distribution for the two biased models
        # >>> m_mult.meval('stationary_distribution')
        # [array([ 0.73105858,  0.26894142]), array([ 0.88079708,  0.11920292])]
        # We set up multiple Markov state models for different temperatures
        # and group them in a MultiStationaryModel.
        # >>> import numpy as np
        # >>> from pyemma.msm import MSM
        # >>> from pyemma.thermo import MEMM
        # >>> b = 20  # transition barrier in kJ / mol
        # >>> temps = np.arange(300, 500, 25)  # temperatures 300 to 500 K
        # >>> p_trans = [np.exp(- b / kT) for kT in 0.00831*temps ]
        # >>> # build MSMs for different temperatures
        # >>> msms = [MSM(P=np.array([[1.0-p, p], [p, 1.0-p]])) for p in p_trans]
        # >>> # build Multi-MSM
        # >>> msm_mult = MEMM(pi=msms[0].stationary_distribution, label='300 K', models=msms)
        # Compute the timescales and see how they decay with temperature
        # Greetings to Arrhenius.
        # >>> np.hstack(msm_mult.meval('timescales'))
        # array([ 1523.83827932,   821.88040004,   484.06386176,   305.87880068,
        #          204.64109413,   143.49286817,   104.62539128,    78.83331598])
        # !! END OF INCOMPATIBLE PART !!
        return [_call_member(M, f, *args, **kw) for M in self.models]
951970
def ngram(n, iter_tokens):
    
    z = len(iter_tokens)
    return (iter_tokens[i:i+n] for i in range(z-n+1))
942246
def _array_star(args):
    
    fn, cls, genelist, kwargs = args
    return _array(fn, cls, genelist, **kwargs)