diff --git a/PyChem/PyChem/liblifetime.py b/PyChem/PyChem/liblifetime.py index a6b8f60..7d37abe 100644 --- a/PyChem/PyChem/liblifetime.py +++ b/PyChem/PyChem/liblifetime.py @@ -1,42 +1,83 @@ from numpy import * def WriteParams(f): """ Write lifetime parameters. For the moment, the values are simply hard coded. :param f: file :type f: the file handler :returns:none """ f.write("#### Livetime ####\n") f.write("\n") f.write("-40.1107251082866 5.50992173040633 0.782431795366473\n") f.write(" 141.929566656232 -15.88948185575660 -3.255779246324010\n") f.write("-261.365531429482 17.07350618651300 9.866058678313810\n") f.write("\n") -def WriteH5Params(ParentGroup): +def WriteH5Params(ParentGroup,model="Poirier"): """ Write lifetime parameters. For the moment, the values are simply hard coded. :param f: file :type f: the file handler :returns:none """ - coeff_z=array([ - [-40.1107251082866, 5.50992173040633, 0.782431795366473], - [ 141.929566656232, -15.88948185575660, -3.255779246324010], - [-261.365531429482, 17.07350618651300, 9.866058678313810] - ]) + if model=="Poirier": + + """" + This is the default model used in the Poirier Thesis + """ + + coeff_z=array([ + [-40.1107251082866, 5.50992173040633, 0.782431795366473], + [ 141.929566656232, -15.88948185575660, -3.255779246324010], + [-261.365531429482, 17.07350618651300, 9.866058678313810] + ]) + + elif model=="POPIII20170112": + + """This is a very simple model designed to extrapolate lifetime of very massive stars. + It is based on nothing else than a very simple extrapolation of the Poirier Model and + corresonds to + + tau = b * (m)**a # tau in Myr, m in Msol + + + u_lt0 is computed for units compatibility reasons with the Poirier model + """ + + + SEC_PER_MEGAYEAR=3.155e13 + UnitTime_in_s = 148849920000000.0 + UnitTime_in_Megayears0=UnitTime_in_s / SEC_PER_MEGAYEAR; + u_lt0 = -log10(UnitTime_in_Megayears0*1e6); + + a = -0.5 + b = 7e-4 * 1e10 / 1e6 + + coeff_z=array([ + [0, 0, 0], + [0, 0, a], + [0, 0, log10(b)-u_lt0] + ]) + + else: + print "no stellar lifetime model corresonds to %s"%(model) + print "!!! STOP !!!" + sys.exit() + + + print "Lifetime model using : %s"%(model) Group = ParentGroup.create_group("LiveTimes") Dataset = Group.create_dataset("coeff_z", data=coeff_z) Dataset.attrs['nx'] = coeff_z.shape[0] Dataset.attrs['ny'] = coeff_z.shape[1] diff --git a/PyChem/scripts/pychem_SSP_discret_evolution3 b/PyChem/scripts/pychem_SSP_discret_evolution3 index 34c291c..649de43 100755 --- a/PyChem/scripts/pychem_SSP_discret_evolution3 +++ b/PyChem/scripts/pychem_SSP_discret_evolution3 @@ -1,1469 +1,1483 @@ #!/usr/bin/env python from optparse import OptionParser from PyChem import chemistry from numpy import * import sys import string import Ptools as pt def parse_options(): usage = "usage: %prog [options] file" parser = OptionParser(usage=usage) parser = pt.add_limits_options(parser) parser = pt.add_log_options(parser) parser = pt.add_postscript_options(parser) parser.add_option("--x", action="store", dest="x", type="string", default = 'Fe', help="x value to plot", metavar=" STRING") parser.add_option("--y", action="store", dest="y", type="string", default = 'Mg', help="y value to plot", metavar=" STRING") parser.add_option("--X", action="store", dest="X", type="string", default = None, help="x value to plot", metavar=" STRING") parser.add_option("--Y", action="store", dest="Y", type="string", default = None, help="y value to plot", metavar=" STRING") parser.add_option("--dt", action="store", dest="dt", type="float", default = 0.1, help="dt", metavar=" FLOAT") parser.add_option("--tmin", action="store", dest="tmin", type="float", default = 1, help="tmin", metavar=" FLOAT") parser.add_option("--tmax", action="store", dest="tmax", type="float", default = 14000, help="tmax", metavar=" FLOAT") parser.add_option("--mstar", action="store", dest="mstar", type="float", default = 1e5, help="initial mass of the SSP in solar mass", metavar=" FLOAT") parser.add_option("--mgasfact", action="store", dest="mgasfact", type="float", default = 1.5, help="ratio between mgas and mstars", metavar=" FLOAT") parser.add_option("--tstar", action="store", dest="tstar", type="float", default = 0, help="formation time of the SSP", metavar=" FLOAT") parser.add_option("-o", action="store", dest="obs", type="string", default = 'Y', help="observable to plot", metavar=" STRING") parser.add_option("--timeunit", action="store", dest="timeunit", type="string", default = None, help="unit of time", metavar=" STRING") parser.add_option("--NumberOfTables", action="store", dest="NumberOfTables", type="int", default = 1, help="NumberOfTables", metavar=" INT") parser.add_option("--DefaultTable", action="store", dest="DefaultTable", type="int", default = 0, help="DefaultTable", metavar=" INT") parser.add_option("--seed", action="store", dest="seed", type="int", default = None, help="initial random seed", metavar=" INT") parser.add_option("-Z","--Z", action="store", dest="Z", type="float", default = 0.02, help="metallicity", metavar=" FLOAT") parser.add_option("--plot_continuous_only", action="store_true", dest="plot_continuous_only", default = False, help="plot continuous values only", metavar=" FLOAT") parser.add_option("--plot_discrete_only", action="store_true", dest="plot_discrete_only", default = False, help="plot discrete values only", metavar=" FLOAT") parser.add_option("--optimal_sampling", action="store_true", dest="optimal_sampling", default = False, help="use optimal sampling for the discrete case", metavar=" FLOAT") (options, args) = parser.parse_args() pt.check_files_number(args) files = args return files,options ######################################################################## # M A I N ######################################################################## SOLAR_MASS = 1.989e33 UnitLength_in_cm = 3.085e+21 UnitMass_in_g = 1.989e+43 UnitVelocity_in_cm_per_s = 20725573.785998672 UnitTime_in_s = 148849920000000.0 feedback_energy = 10**51 /UnitMass_in_g/(UnitVelocity_in_cm_per_s)**2 # parse options files,opt = parse_options() if opt.seed!=None: random.seed(opt.seed) # units opt.ux = 1.0 if opt.timeunit == None: opt.ux = 1.0 elif opt.timeunit == 'Myr': opt.ux = UnitTime_in_s/(3600*24*365)/1e6 elif opt.timeunit == 'Gyr': opt.ux = UnitTime_in_s/(3600*24*365)/1e9 # to plot opt.obs = string.split(opt.obs,',') if opt.X!=None: opt.X = string.split(opt.X,',') if opt.Y!=None: opt.Y = string.split(opt.Y,',') # init chimie chemistry.init_chimie(files[0],opt.NumberOfTables,opt.DefaultTable) elts = chemistry.get_elts_labels() nelts = chemistry.get_nelts() print elts # some parameters dt = opt.dt tmax = opt.tmax / opt.ux t_star = opt.tstar # formation time of the SSP M0 = opt.mstar *SOLAR_MASS/UnitMass_in_g # gas mass of the SSP (in code unit) Mstacont = M0 # star mass Mstaj = M0 Mstajcont = M0 Mgascont = opt.mgasfact*M0 # gas mass Mgasj = opt.mgasfact*M0 # gas mass Mgasjcont = opt.mgasfact*M0 # gas mass # stellar initial element abondance metals = zeros(chemistry.get_nelts(),float) metals[-1]= opt.Z #metals[0] = 0.001771 *1 #metals[1] = 0.00091245 *1 #metals[2] = 0.0 *1 #metals[3] = 0.02 *1 # metal abondance (=metallicity z) # gas initial element abondance metalgcont = zeros(chemistry.get_nelts(),float) metalgj = zeros(chemistry.get_nelts(),float) metalgjcont = zeros(chemistry.get_nelts(),float) # minimum live time of a star minlivetime = chemistry.star_lifetime(metals[-1],chemistry.get_Mmax()) maxlivetime = chemistry.star_lifetime(metals[-1],chemistry.get_Mmin()) # times times = arange(dt,tmax,dt).astype(float) nt = len(times) # output arrays TIMES = zeros(nt,float) # time NSNIas = zeros(nt,float) # number of SNIa NSNIIs = zeros(nt,float) # number of SNII NSNIaconts = zeros(nt,float) # number of SNIa (continuous value) NSNIIconts = zeros(nt,float) # number of SNII (continuous value) NDYINs = zeros(nt,float) # number of dying stars NDYINconts = zeros(nt,float) # number of dying stars Escont = zeros(nt,float) # injected energy Es = zeros(nt,float) # injected energy Mecont = zeros(nt,float) # ejected mass (continuous) Mej = zeros(nt,float) # ejected mass (computed from NSN) Mejcont = zeros(nt,float) # ejected mass (computed from NSN,continuous) Mstaconts = zeros(nt,float) # star mass (continuous) Mstajs = zeros(nt,float) # star mass Mstajconts = zeros(nt,float) # star mass Mgasconts = zeros(nt,float) # gas mass Mgasjs = zeros(nt,float) # gas mass Mgasjconts = zeros(nt,float) # gas mass # elements to follow # ejected dataecont = {} dataej = {} dataejcont = {} Meltconts = {} # ejected elt mass Meltjs = {} # ejected elt mass Meltjconts = {} # ejected elt mass for elt in elts: dataecont[elt] = zeros(nt,float) dataej[elt] = zeros(nt,float) dataejcont[elt] = zeros(nt,float) Meltconts[elt] = zeros(nt,float) # ejected elt mass Meltjs[elt] = zeros(nt,float) # ejected elt mass Meltjconts[elt] = zeros(nt,float) # ejected elt mass # gas datagcont = {} datagj = {} datagjcont = {} for elt in elts: datagcont[elt] = zeros(nt,float) datagj[elt] = zeros(nt,float) datagjcont[elt] = zeros(nt,float) # optimal Mecl = opt.mstar # in Msol m_max_star = chemistry.optimal_sampling_init_norm(Mecl) current_mass = m_max_star # in Msol M_diced = 0 ############################################################################################## # # main loop # ############################################################################################## for ti,time in enumerate(times): tpdt = time + dt if (tpdt - t_star >= minlivetime) and (tpdt - t_star < maxlivetime): m1 = chemistry.star_mass_from_age(metals[-1],tpdt - t_star) if (time - t_star >= minlivetime): m2 = chemistry.star_mass_from_age(metals[-1],time - t_star) else: m2 = chemistry.get_Mmax() if (m1>m2): m1=m2 raise "m1>m2 !!!" ############################ # optimal sampling stuffs if opt.optimal_sampling: NSNII_opt = 0 NSNIa_opt = 0 NDYIN_opt = 0 RSNIa=0 while ( current_mass*SOLAR_MASS/UnitMass_in_g > m1 ): if chemistry.optimal_sampling_stop_loop(current_mass): print "optimal_sampling_stop_loop !!!" sys.exit() #print "current mass = %g m1=%g m2=%g"%(current_mass,m1*1e10,m2*1e10) if ( (chemistry.get_SNII_Mmin() <= current_mass*SOLAR_MASS/UnitMass_in_g) and (current_mass*SOLAR_MASS/UnitMass_in_g <= chemistry.get_SNII_Mmax()) ): #print m1*1e10,m2*1e10,"one SNII" NSNII_opt = NSNII_opt+1 if ( (chemistry.get_DYIN_Mmin() <= current_mass*SOLAR_MASS/UnitMass_in_g) and (current_mass*SOLAR_MASS/UnitMass_in_g <= chemistry.get_DYIN_Mmax()) ): #print m1*1e10,m2*1e10,"one DYIN" NDYIN_opt = NDYIN_opt+1 RSNIa = RSNIa + chemistry.SNIa_single_rate(current_mass*SOLAR_MASS/UnitMass_in_g) current_mass = chemistry.optimal_sampling_get_next_mass(current_mass) # because the function SNIa_single_rate returns a continuous value, # we need to discretize it if RSNIa>0: fNSNIa = RSNIa-floor(RSNIa) NSNIa = floor(RSNIa) if random.random() < fNSNIa: NSNIa = NSNIa+1 NSNIa_opt = NSNIa # number of SNIa NSNIa = chemistry.SNIa_rate(m1,m2)*M0 # number of SNII NSNII = chemistry.SNII_rate(m1,m2)*M0 # number of dying stars NDYIN = chemistry.DYIN_rate(m1,m2)*M0 ######################### # discrete stuff NSNIIcont = NSNII NSNIacont = NSNIa NDYINcont = NDYIN # compute discrete values (SNIa) fNSNIa = NSNIa-floor(NSNIa) # fraction NSNIa = floor(NSNIa) # integer if random.random() < fNSNIa: NSNIa = NSNIa+1 # compute discrete values (SNII) fNSNII = NSNII-floor(NSNII) # fraction NSNII = floor(NSNII) # integer if random.random() < fNSNII: NSNII = NSNII+1 # compute discrete values (DYIN) fNDYIN = NDYIN-floor(NDYIN) # fraction NDYIN = floor(NDYIN) # integer if random.random() < fNDYIN: NDYIN = NDYIN+1 # #########################$ if opt.optimal_sampling: NSNII = NSNII_opt NSNIa = NSNIa_opt NDYIN = NDYIN_opt ################################################## # mass and yields ejection ################################################## # standard way (continuous explosion) (_cont) EjectedMass = chemistry.Total_mass_ejection_Z(m1,m2,M0,metals) TotalEjectedEltMasscont = EjectedMass[2:] TotalEjectedGasMasscont = EjectedMass[0] # new way (discrete explosion) (_j) EjectedMassj = chemistry.Total_single_mass_ejection_Z(0.5*(m1+m2),metals,float(NSNII),float(NSNIa),float(NDYIN)) TotalEjectedEltMassj = EjectedMassj[2:] TotalEjectedGasMassj = EjectedMassj[0] # new way (continuous explosion) (_jcont) EjectedMassjcont = chemistry.Total_single_mass_ejection_Z(0.5*(m1+m2),metals,NSNIIcont,NSNIacont,NDYINcont) TotalEjectedEltMassjcont = EjectedMassjcont[2:] TotalEjectedGasMassjcont = EjectedMassjcont[0] # other check #SNII_EjectedMass = chemistry.SNII_Total_single_mass_ejection( 0.5*(m1+m2) ,metals) * NSNII #SNIa_EjectedMass = chemistry.SNIa_Total_single_mass_ejection( 0.5*(m1+m2) ,metals) * NSNIa #TotalEjectedEltMassjcont = SNII_EjectedMass[2:] + SNIa_EjectedMass[2:] #TotalEjectedGasMassjcont = SNII_EjectedMass[0] + SNIa_EjectedMass[0] ''' # all dying stars (= SNII + lighter masses) DYIN_TotalEjectedGasMassj = chemistry.DYIN_Total_single_mass_ejection( 0.5*(m1+m2) ,metals)[0] * NSNIa # !!! need to add SNIa DYIN_TotalEjectedGasMassjcont = chemistry.DYIN_Total_single_mass_ejection( 0.5*(m1+m2) ,metals)[0] * NSNIacont # !!! need to add SNIa ''' ################################################## # energy ################################################## # energy injected TotalInjectedEnergycont = feedback_energy* (NSNIacont + NSNIIcont) TotalInjectedEnergy = feedback_energy* (NSNIa + NSNII) # mass fraction of ejected elements emetalcont = TotalEjectedEltMasscont /TotalEjectedGasMasscont emetalj = TotalEjectedEltMassj /TotalEjectedGasMassj emetaljcont= TotalEjectedEltMassjcont/TotalEjectedGasMassjcont # metal enrichment metalgcont = ((metalgcont *Mgascont) + TotalEjectedEltMasscont ) metalgj = ((metalgj *Mgasj) + TotalEjectedEltMassj ) metalgjcont = ((metalgjcont*Mgasjcont) + TotalEjectedEltMassjcont ) # Total Mass Mstacont = Mstacont - TotalEjectedGasMasscont Mstaj = Mstaj - TotalEjectedGasMassj Mstajcont = Mstajcont - TotalEjectedGasMassjcont # Gas Mass Mgascont = Mgascont + TotalEjectedGasMasscont Mgasj = Mgasj + TotalEjectedGasMassj Mgasjcont = Mgasjcont + TotalEjectedGasMassjcont # metal enrichment metalgcont = metalgcont /Mgascont metalgj = metalgj /Mgasj metalgjcont = metalgjcont/Mgasjcont # record output TIMES[ti] = time NSNIas[ti] = NSNIa NSNIIs[ti] = NSNII NSNIaconts[ti] = NSNIacont NSNIIconts[ti] = NSNIIcont NDYINs[ti] = NDYIN NDYINconts[ti] = NDYINcont Escont[ti] = TotalInjectedEnergycont Es[ti] = TotalInjectedEnergy Mecont[ti] = TotalEjectedGasMasscont Mej[ti] = TotalEjectedGasMassj Mejcont[ti]= TotalEjectedGasMassjcont Mstaconts[ti] = Mstacont # stellar mass Mstajs[ti] = Mstaj # stellar mass Mstajconts[ti] = Mstajcont # stellar mass Mgasconts[ti] = Mgascont Mgasjs[ti] = Mgasj Mgasjconts[ti] = Mgasjcont for j,elt in enumerate(elts): dataecont[elt][ti] = emetalcont[j] datagcont[elt][ti] = metalgcont[j] dataej[elt][ti] = emetalj[j] datagj[elt][ti] = metalgj[j] dataejcont[elt][ti] = emetaljcont[j] datagjcont[elt][ti] = metalgjcont[j] Meltconts[elt][ti] = TotalEjectedEltMasscont[j] Meltjs[elt][ti] = TotalEjectedEltMassj[j] Meltjconts[elt][ti] = TotalEjectedEltMassjcont[j] # remove time=0 values c = (TIMES!=0) TIMES = compress(c,TIMES) * opt.ux NSNIas = compress(c,NSNIas) NSNIIs = compress(c,NSNIIs) NSNIaconts = compress(c,NSNIaconts) NSNIIconts = compress(c,NSNIIconts) NDYINs = compress(c,NDYINs) NDYINconts = compress(c,NDYINconts) Escont = compress(c,Escont) Es = compress(c,Es) Mecont = compress(c,Mecont) Mej = compress(c,Mej) Mejcont = compress(c,Mejcont) Mstaconts = compress(c,Mstaconts) Mstajs = compress(c,Mstajs) Mstajconts= compress(c,Mstajconts) Mgasconts = compress(c,Mgasconts) Mgasjs = compress(c,Mgasjs) Mgasjconts= compress(c,Mgasjconts) for elt in elts: dataecont[elt] = compress(c,dataecont[elt]) datagcont[elt] = compress(c,datagcont[elt]) dataej[elt] = compress(c,dataej[elt]) datagj[elt] = compress(c,datagj[elt]) dataejcont[elt] = compress(c,dataejcont[elt]) datagjcont[elt] = compress(c,datagjcont[elt]) Meltconts[elt] = compress(c,Meltconts[elt]) Meltjs[elt] = compress(c,Meltjs[elt]) Meltjconts[elt] = compress(c,Meltjconts[elt]) # get Solar Abundances SolarAbun = chemistry.get_elts_SolarMassAbundances() Xsconts = {} Xsjs = {} Xsjconts= {} Xgconts = {} Xgjs = {} Xgjconts= {} for elt in elts: Xsconts[elt] = log10(dataecont[elt] / SolarAbun[elt] + 1.0e-10) Xgconts[elt] = log10(datagcont[elt] / SolarAbun[elt] + 1.0e-10) Xsjs[elt] = log10(dataej[elt] / SolarAbun[elt] + 1.0e-10) Xgjs[elt] = log10(datagj[elt] / SolarAbun[elt] + 1.0e-10) Xsjconts[elt] = log10(dataejcont[elt] / SolarAbun[elt] + 1.0e-10) Xgjconts[elt] = log10(datagjcont[elt] / SolarAbun[elt] + 1.0e-10) Xscont = Xsconts[opt.x] Yscont = Xsconts[opt.y] Xgcont = Xgconts[opt.x] Ygcont = Xgconts[opt.y] Xsj = Xsjs[opt.x] Ysj = Xsjs[opt.y] Xgj = Xgjs[opt.x] Ygj = Xgjs[opt.y] Xsjcont = Xsjconts[opt.x] Ysjcont = Xsjconts[opt.y] Xgjcont = Xgjconts[opt.x] Ygjcont = Xgjconts[opt.y] ################################################################################################################## # # plot # ################################################################################################################## # some plotting options opt.plot_continuous=True opt.plot_discrete=True opt.plot_semicontnuous=False if opt.plot_continuous_only: opt.plot_continuous=True opt.plot_discrete=False opt.plot_semicontnuous=False if opt.plot_discrete_only: opt.plot_continuous=False opt.plot_discrete=True opt.plot_semicontnuous=False pt.InitPlot([],opt) pt.pcolors fig = pt.gcf() #fig.subplots_adjust(left=0.02) #fig.subplots_adjust(right=0.99) #fig.subplots_adjust(bottom=0.1) #fig.subplots_adjust(top=0.99) #fig.subplots_adjust(wspace=0.0) fig.subplots_adjust(hspace=0.25) np = len(opt.obs) def PlotYields(i): xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit) ylabel = r"$\rm{[%s/H],[%s/H],[%s/%s]}$"%(opt.x,opt.y,opt.y,opt.x) pt.subplot(np,1,i) datas = [] if opt.plot_continuous: data = pt.DataPoints(TIMES,Xscont,tpe='line',label='Xs',color='b',linestyle='-') datas.append(data) data = pt.DataPoints(TIMES,Yscont,tpe='line',label='Ys',color='g',linestyle='-') datas.append(data) data = pt.DataPoints(TIMES,Yscont-Xscont,tpe='line',label='Ys-Xs',color='r',linestyle='-',linewidth=2) datas.append(data) if opt.plot_discrete: data = pt.DataPoints(TIMES,Xsj,tpe='line',label='Xs',color='b',linestyle='--') datas.append(data) data = pt.DataPoints(TIMES,Ysj,tpe='line',label='Ys',color='g',linestyle='--') datas.append(data) data = pt.DataPoints(TIMES,Ysj-Xsj,tpe='line',label='Ys-Xs',color='r',linestyle='--',linewidth=2) datas.append(data) if opt.plot_semicontnuous: data = pt.DataPoints(TIMES,Xsjcont,tpe='line',label='Xs',color='b',linestyle=':') datas.append(data) data = pt.DataPoints(TIMES,Ysjcont,tpe='line',label='Ys',color='g',linestyle=':') datas.append(data) data = pt.DataPoints(TIMES,Ysjcont-Xsjcont,tpe='line',label='Ys-Xs',color='r',linestyle=':',linewidth=2) datas.append(data) for d in datas: pt.plot(d.x,d.y,color=d.color,ls=d.linestyle,lw=d.linewidth) xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log) ymin = max(-2,ymin) # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotGasElts(i): xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit) ylabel = r"$\rm{gas\,[%s/H],[%s/H],[%s/%s]}$"%(opt.x,opt.y,opt.y,opt.x) pt.subplot(np,1,i) datas = [] data = pt.DataPoints(TIMES,Xgcont,tpe='line',label='Xs',color='b') datas.append(data) data = pt.DataPoints(TIMES,Ygcont,tpe='line',label='Ys',color='g') datas.append(data) data = pt.DataPoints(TIMES,Ygcont-Xgcont,tpe='line',label='Ys-Xs',color='r',linestyle='--') datas.append(data) data = pt.DataPoints(TIMES,Xgj,tpe='line',label='Xs',color='b',linestyle='--') datas.append(data) data = pt.DataPoints(TIMES,Ygj,tpe='line',label='Ys',color='g',linestyle='--') datas.append(data) data = pt.DataPoints(TIMES,Ygj-Xgj,tpe='line',label='Ys-Xs',color='r',linestyle='-') datas.append(data) data = pt.DataPoints(TIMES,Xgjcont,tpe='line',label='Xs',color='b') #datas.append(data) data = pt.DataPoints(TIMES,Ygjcont,tpe='line',label='Ys',color='g') #datas.append(data) data = pt.DataPoints(TIMES,Ygjcont-Xgjcont,tpe='line',label='Ys-Xs',color='r',linestyle='--') #datas.append(data) for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color,ls=d.linestyle) xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log) ymin = -2 # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotEnergy(i): xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit) ylabel = r"$\rm{Supernova\,Energy}$" pt.subplot(np,1,i) datas = [] data = pt.DataPoints(TIMES,Escont,tpe='line',label='Es',color='k') datas.append(data) data = pt.DataPoints(TIMES,Es,tpe='line',label='Es',color='r') datas.append(data) for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color) if opt.log==None: log='y' else: log = opt.log + 'y' xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log) # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotCumEnergy(i): xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit) ylabel = r"$\rm{Supernova\,Energy}$" pt.subplot(np,1,i) CumEscont = add.accumulate(Escont) CumEs = add.accumulate(Es) datas = [] data = pt.DataPoints(TIMES,CumEscont,tpe='line',label='Es',color='k') datas.append(data) data = pt.DataPoints(TIMES,CumEs,tpe='line',label='Es',color='r') datas.append(data) for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color) if opt.log==None: log='y' else: log = opt.log + 'y' xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log) # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotEjectedMass(i): xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit) ylabel = r"$\rm{Ejected\,Mass}$" pt.subplot(np,1,i) datas = [] data = pt.DataPoints(TIMES,Mecont,tpe='line',label='Es',color='r') datas.append(data) data = pt.DataPoints(TIMES,Mej,tpe='line',label='Es',color='k') datas.append(data) data = pt.DataPoints(TIMES,Mejcont,tpe='line',label='Es',color='g') #datas.append(data) for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color) if opt.log==None: log=opt.log else: log = opt.log xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log) # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotCumEjectedMass(i): xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit) ylabel = r"$\rm{Cumulative\,Ejected\,Mass}$" pt.subplot(np,1,i) cumcont = add.accumulate(Mecont) cumj = add.accumulate(Mej) cumjcont = add.accumulate(Mejcont) datas = [] data = pt.DataPoints(TIMES,cumcont,tpe='line',label='Es',color='r') datas.append(data) data = pt.DataPoints(TIMES,cumj,tpe='line',label='Es',color='k') datas.append(data) data = pt.DataPoints(TIMES,cumjcont,tpe='line',label='Es',color='g') #datas.append(data) for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color) if opt.log==None: log=opt.log else: log = opt.log xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log) # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotStellarMass(i): xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit) ylabel = r"$\rm{Stellar\,Mass}$" pt.subplot(np,1,i) + ''' datas = [] data = pt.DataPoints(TIMES,Mstaconts,tpe='line',label='Es',color='r') datas.append(data) data = pt.DataPoints(TIMES,Mstajs,tpe='line',label='Es',color='k') datas.append(data) data = pt.DataPoints(TIMES,Mstajconts,tpe='line',label='Es',color='g') #datas.append(data) + ''' + + + datas = [] + + if opt.plot_discrete: + data = pt.DataPoints(TIMES,Mstajs,tpe='line',label='Es',color='k') + datas.append(data) + + if opt.plot_continuous: + data = pt.DataPoints(TIMES,Mstaconts,tpe='line',label='Es',color='r') + datas.append(data) + for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color) if opt.log==None: log = opt.log else: log = opt.log opt.ymin = 0 xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log) # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotNDYIN(i): xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit) ylabel = r"$\rm{Number\,dying\,stars}$" pt.subplot(np,1,i) datas = [] data = pt.DataPoints(TIMES,NDYINs,tpe='line',label='Es',color='k') datas.append(data) data = pt.DataPoints(TIMES,NDYINconts,tpe='line',label='Es',color='r') datas.append(data) for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color) log=opt.log xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log) # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotCumNDYIN(i): xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit) ylabel = r"$\rm{Cum\,Number\,of\,dying\,stars}$" pt.subplot(np,1,i) cum = add.accumulate(NDYINs) cumc= add.accumulate(NDYINconts) datas = [] data = pt.DataPoints(TIMES,cum,tpe='line',label='Es',color='k') datas.append(data) data = pt.DataPoints(TIMES,cumc,tpe='line',label='Es',color='r') datas.append(data) for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color) log=opt.log xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log) # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotNSNII(i): xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit) ylabel = r"$\rm{Number\,of\,SNII}$" pt.subplot(np,1,i) datas = [] if opt.plot_discrete: data = pt.DataPoints(TIMES,NSNIIs,tpe='line',label='Es',color='k') datas.append(data) if opt.plot_continuous: data = pt.DataPoints(TIMES,NSNIIconts,tpe='line',label='Es',color='r') datas.append(data) for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color) log=opt.log xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log) # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotCumNSNII(i): xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit) ylabel = r"$\rm{Cum\,Number\,of\,SNII}$" pt.subplot(np,1,i) cum = add.accumulate(NSNIIs) cumc= add.accumulate(NSNIIconts) datas = [] if opt.plot_discrete: data = pt.DataPoints(TIMES,cum,tpe='line',label='Es',color='k') datas.append(data) if opt.plot_continuous: data = pt.DataPoints(TIMES,cumc,tpe='line',label='Es',color='r') datas.append(data) for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color) log=opt.log xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log) # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotNSNIa(i): xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit) ylabel = r"$\rm{Number\,of\,SNIa}$" pt.subplot(np,1,i) datas = [] if opt.plot_discrete: data = pt.DataPoints(TIMES,NSNIas,tpe='line',label='Es',color='k') datas.append(data) if opt.plot_continuous: data = pt.DataPoints(TIMES,NSNIaconts,tpe='line',label='Es',color='r') datas.append(data) for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color) log=opt.log xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log) # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotCumNSNIa(i): xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit) ylabel = r"$\rm{Cum\,Number\,of\,SNIa}$" pt.subplot(np,1,i) cum = add.accumulate(NSNIas) cumc= add.accumulate(NSNIaconts) datas = [] data = pt.DataPoints(TIMES,cum,tpe='line',label='Es',color='k') datas.append(data) data = pt.DataPoints(TIMES,cumc,tpe='line',label='Es',color='r') datas.append(data) for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color) log=opt.log xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log) # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotNSNIIDYIN(i): xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit) ylabel = r"$\rm{Number\,of\,SNII\,and\,dying}$" pt.subplot(np,1,i) datas = [] data = pt.DataPoints(TIMES,NSNIIs+NDYINs,tpe='line',label='Es',color='k') datas.append(data) data = pt.DataPoints(TIMES,NSNIIconts+NDYINconts,tpe='line',label='Es',color='r') datas.append(data) for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color) log=opt.log xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log) # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotFevsTime(i): xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit) ylabel = r"$\rm{Fe}$" pt.subplot(np,1,i) datas = [] data = pt.DataPoints(TIMES,Ygj,tpe='line',label='Es',color='k') datas.append(data) data = pt.DataPoints(TIMES,Ygjcont,tpe='line',label='Es',color='r') datas.append(data) for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color) log=opt.log xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log) # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotMassFevsTime(i): ''' !!! this is not very well written. Need to be improved ''' xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit) ylabel = r"$\rm{Fe}$" pt.subplot(np,1,i) Xgj = add.accumulate(Meltconts['Fe']) Ygj = add.accumulate(Meltconts['Mg']) #Xgj = log10(Meltconts['Mg']+1e-20) #Ygj = log10(Meltconts['Fe']+1e-20) pt.plot(TIMES,Xgj) pt.plot(TIMES,Ygj) datas = [] data = pt.DataPoints(TIMES,Ygj,tpe='line',label='Es',color='k') #datas.append(data) #data = pt.DataPoints(TIMES,Ygjcont,tpe='line',label='Es',color='r') #datas.append(data) for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color) log=None #xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log) xmin = 0 xmax = 14000 ymin = -20 ymax = -8 # plot #pt.SetAxis(xmin,xmax,ymin,ymax,log=log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotDiff(i): xlabel = r"$\rm{[%s/H]$"%(opt.x) ylabel = r"$\rm{[%s/%s]}$"%(opt.y,opt.x) pt.subplot(np,1,i) datas = [] if opt.plot_continuous: data = pt.DataPoints(Xgcont,Ygcont-Xgcont ,tpe='line',label='Xs',color='r') datas.append(data) if opt.plot_discrete: data = pt.DataPoints(Xgj,Ygj-Xgj ,tpe='line',label='Xs',color='k') datas.append(data) if opt.plot_semicontnuous: data = pt.DataPoints(Xgjcont,Ygjcont-Xgjcont,tpe='line',label='Xs',color='g') datas.append(data) for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color) log = opt.log #xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log) xmin = -4 xmax = 0 ymin = -2 ymax = 2 # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotDXvsDY(i): # X = X/H X1 = opt.X[0] X2 = opt.X[1] print X1,X2 Y1 = opt.Y[0] Y2 = opt.Y[1] print Y1,Y2 if X2=='H': xcont = Xgconts[X1] xj = Xgjs[X1] xjcont = Xgjconts[X1] xlabel = r"$\rm{[%s/%s]}$"%(X1,"H") else: xcont = Xgconts[X1] - Xgconts[X2] xj = Xgjs[X1] - Xgjs[X2] xjcont = Xgjconts[X1] - Xgjconts[X2] xlabel = r"$\rm{[%s/%s]}$"%(X1,X2) #Xgcont,Ygcont-Xgcont #Xgconts[opt.x] if Y2=='H': ycont = Xgconts[Y1] yj = Xgjs[Y1] yjcont = Xgjconts[Y1] ylabel = r"$\rm{[%s/%s]}$"%(Y1,"H") else: ycont = Xgconts[Y1] - Xgconts[Y2] yj = Xgjs[Y1] - Xgjs[Y2] yjcont = Xgjconts[Y1] - Xgjconts[Y2] ylabel = r"$\rm{[%s/%s]}$"%(Y1,Y2) pt.subplot(np,1,i) datas = [] if opt.plot_continuous: data = pt.DataPoints(xcont,ycont, tpe='line',label='Xs',color='r') datas.append(data) if opt.plot_discrete: data = pt.DataPoints(xj,yj, tpe='line',label='Xs',color='k') datas.append(data) if opt.plot_semicontnuous: data = pt.DataPoints(xjcont,yjcont,tpe='line',label='Xs',color='g') datas.append(data) for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color) log = None xmin = opt.xmin xmax = opt.xmax ymin = opt.ymin ymax = opt.ymax xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(xmin,xmax,ymin,ymax,datas,log) print xlabel print ylabel # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) def PlotStarMass(i): xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit) ylabel = r"$\rm{Mass}$" pt.subplot(np,1,i) datas = [] data = pt.DataPoints(TIMES,Mstas,tpe='line',label='Mstas',color='k') datas.append(data) for d in datas: if d.tpe=='line' or d.tpe=='both': pt.plot(d.x,d.y,color=d.color) log = opt.log xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log) # plot pt.SetAxis(xmin,xmax,ymin,ymax,log=log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) pdict = {'Y':PlotYields,'Yg':PlotGasElts,'ESN':PlotEnergy,'CumESN':PlotCumEnergy,'NSNII':PlotNSNII,'CumNSNII':PlotCumNSNII,'NSNIa':PlotNSNIa,'CumNSNIa':PlotCumNSNIa,'NDYIN':PlotNDYIN,'CumNDYIN':PlotCumNDYIN, 'D':PlotDiff,'MSTAR':PlotStarMass,'EjMass':PlotEjectedMass,'CumEjMass':PlotCumEjectedMass,'StellarMass':PlotStellarMass, 'NSNIIDYIN':PlotNSNIIDYIN,'Fe':PlotFevsTime,'DXvsDY':PlotDXvsDY,'MFe':PlotMassFevsTime} ip = 1 for obs in opt.obs: pdict[obs](ip) ip = ip + 1 pt.EndPlot(files,opt) diff --git a/PyChem/scripts/pychem_generate_hdf5_parameters b/PyChem/scripts/pychem_generate_hdf5_parameters index c585825..4ed1f39 100755 --- a/PyChem/scripts/pychem_generate_hdf5_parameters +++ b/PyChem/scripts/pychem_generate_hdf5_parameters @@ -1,196 +1,199 @@ #!/usr/bin/env python from optparse import OptionParser from PyChem import liblifetime from PyChem import libimf from PyChem import libSNIa from PyChem import libSNII from PyChem import libDYIN from PyChem import libmetalsSNII from PyChem import libmetalsSNIa from PyChem import libmetalsDYIN from PyChem import libSolarAbundances import sys import h5py import getpass import time def parse_options(): usage = "usage: %prog [options] file" parser = OptionParser(usage=usage) parser.add_option("-o","--outputfile", action="store", dest="outputfile", type="string", default = 'mychimie.dat', help="Output file", metavar=" NAME") parser.add_option("-p","--parameterfile", action="store", dest="parameterfile", type="string", default = 'chimieparam.py', help="Parameter file", metavar=" NAME") (options, args) = parser.parse_args() files = args return files,options ########################################################### # Main ########################################################### version = "v0.0" files,opt = parse_options() f = h5py.File(opt.outputfile,'w') # compute cmdline cmdline = "" for elt in sys.argv: cmdline = cmdline + " " + elt cmdline = cmdline[1:] try: nelts except NameError: nelts=None ###################### # header ###################### headerGroup = f.create_group("Header") # some attributes headerGroup.attrs["version"] = version headerGroup.attrs["author"] = getpass.getuser() headerGroup.attrs["date"] = time.ctime() headerGroup.attrs["cmdline"] = cmdline ListOfDYINZs=None ListOfDYINYieldsFiles=None ListOfDYINHeliumCoreFiles=None execfile(opt.parameterfile) ###################### # data ###################### dataGroup = f.create_group("Data") # Lifetime -liblifetime.WriteH5Params(dataGroup) +if not 'StellarLifeTimeModel' in locals(): + StellarLifeTimeModel="Poirier" # set default model + +liblifetime.WriteH5Params(dataGroup,model=StellarLifeTimeModel) # IMF libimf.WriteH5Params(dataGroup,Mmin,Mmax,m_s,a_s) # SNII print "SNII" libSNII.WriteH5Params(dataGroup,SNII_Mmin,SNII_Mmax) # DYIN print "DYIN" libDYIN.WriteH5Params(dataGroup,DYIN_Mmin,DYIN_Mmax) # SNIa print "SNIa" libSNIa.WriteH5Params(dataGroup,SNIa_Mpl,SNIa_Mpu,SNIa_a,SNIa_Mdl1,SNIa_Mdu1,SNIa_bb1,SNIa_Mdl2,SNIa_Mdu2,SNIa_bb2) # Metal print "SNII Metal" SNIIelts = libmetalsSNII.WriteH5Params(dataGroup,Mmin,Mmax,m_s,a_s,SNIIYieldsFile,SNIIHeliumCoreFile,elts,nelts,Mmin1SNII,Mmax1SNII,Mmin2SNII,Mmax2SNII,dMSNII,nSNII) # Metal print "SNIa Metal" SNIaelts = libmetalsSNIa.WriteH5Params(dataGroup,SNIaFile,MeanWDMass,elts,nelts) if (SNIIelts!=SNIaelts): raise "SNIIelts !=SNIaelts" elts = SNIaelts # Metal print "DYIN Metal" SNIIelts = libmetalsDYIN.WriteH5Params(dataGroup,Mmin,Mmax,m_s,a_s,DYINYieldsFile,DYINHeliumCoreFile,elts,nelts,Mmin1DYIN,Mmax1DYIN,Mmin2DYIN,Mmax2DYIN,dMDYIN,nDYIN) if ListOfDYINZs!=None and ListOfDYINYieldsFiles!=None and ListOfDYINHeliumCoreFiles!=None: print "DYIN Metal Z" SNIIelts = libmetalsDYIN.WriteH5ParamsWithZ(dataGroup,Mmin,Mmax,m_s,a_s,ListOfDYINYieldsFiles,ListOfDYINHeliumCoreFiles,ListOfDYINZs,elts,nelts,Mmin1DYIN,Mmax1DYIN,Mmin2DYIN,Mmax2DYIN,dMDYIN,nDYIN,ZminDYIN,ZmaxDYIN,nZDYIN) # SolarAbundances print "Solar Abundances" SolarMassAbundances = libSolarAbundances.Get(SolarAbundancesFile,elts) ###################### # data attributes nelts=len(elts) # add attribute to the data dataGroup.attrs["nelts"] = nelts dataGroup.attrs["elts"] = elts dataGroup.attrs["MeanWDMass"] = MeanWDMass dataGroup.attrs["SolarMassAbundances"] = SolarMassAbundances dataGroup.attrs["SNIIYieldsFile"] = str(SNIIYieldsFile) dataGroup.attrs["SNIIHeliumCoreFile"] = str(SNIIHeliumCoreFile) dataGroup.attrs["DYINYieldsFile"] = str(DYINYieldsFile) dataGroup.attrs["DYINHeliumCoreFile"] = str(DYINHeliumCoreFile) dataGroup.attrs["SNIaFile"] = str(SNIaFile) dataGroup.attrs["SolarAbundancesFile"] = str(SolarAbundancesFile) f.close() diff --git a/PyChem/scripts/pychem_lifetime b/PyChem/scripts/pychem_lifetime index 71c8d40..4f51ab5 100755 --- a/PyChem/scripts/pychem_lifetime +++ b/PyChem/scripts/pychem_lifetime @@ -1,107 +1,108 @@ #!/usr/bin/env python from optparse import OptionParser import Ptools as pt from pNbody import * from pNbody import units import string from scipy import optimize from PyChem import chemistry UnitLength_in_cm = 3.085e+21 UnitMass_in_g = 1.989e+43 UnitVelocity_in_cm_per_s = 20725573.785998672 UnitTime_in_s = 148849920000000.0 def parse_options(): usage = "usage: %prog [options] file" parser = OptionParser(usage=usage) parser = pt.add_postscript_options(parser) parser = pt.add_ftype_options(parser) parser = pt.add_reduc_options(parser) parser = pt.add_center_options(parser) parser = pt.add_select_options(parser) parser = pt.add_cmd_options(parser) parser = pt.add_display_options(parser) parser = pt.add_info_options(parser) parser = pt.add_limits_options(parser) parser = pt.add_log_options(parser) parser.add_option("--x", action="store", dest="x", type="string", default = 'r', help="x value to plot", metavar=" STRING") parser.add_option("--y", action="store", dest="y", type="string", default = 'T', help="y value to plot", metavar=" STRING") parser.add_option("--z", action="store", dest="z", type="string", default = None, help="z value to plot", metavar=" STRING") parser.add_option("--legend", action="store_true", dest="legend", default = False, help="add a legend") (options, args) = parser.parse_args() pt.check_files_number(args) files = args return files,options ######################################################################## # MAIN ######################################################################## files,opt = parse_options() file = files[0] chemistry.init_chimie(file) metals = zeros(chemistry.get_nelts(),float) - +print chemistry.get_Mmin() +print chemistry.get_Mmax() minlivetime = chemistry.star_lifetime(metals[-1],chemistry.get_Mmax()) maxlivetime = chemistry.star_lifetime(metals[-1],chemistry.get_Mmin()) print "minlivetime =%g"%(minlivetime) print "maxlivetime =%g"%(maxlivetime) diff --git a/PyChem/scripts/pychem_plot_yieldsAbundances_withZ b/PyChem/scripts/pychem_plot_yieldsAbundances_withZ new file mode 100755 index 0000000..21abcef --- /dev/null +++ b/PyChem/scripts/pychem_plot_yieldsAbundances_withZ @@ -0,0 +1,320 @@ +#!/usr/bin/env python + + + +from optparse import OptionParser +import Ptools as pt +from pNbody import * +from pNbody import units +import string + +from scipy import optimize + +from PyChem import chemistry + +UnitLength_in_cm = 3.085e+21 +UnitMass_in_g = 1.989e+43 +UnitVelocity_in_cm_per_s = 20725573.785998672 +UnitTime_in_s = 148849920000000.0 + + +def parse_options(): + + + usage = "usage: %prog [options] file" + parser = OptionParser(usage=usage) + + parser = pt.add_postscript_options(parser) + parser = pt.add_ftype_options(parser) + parser = pt.add_reduc_options(parser) + parser = pt.add_center_options(parser) + parser = pt.add_select_options(parser) + parser = pt.add_cmd_options(parser) + parser = pt.add_display_options(parser) + parser = pt.add_info_options(parser) + parser = pt.add_limits_options(parser) + parser = pt.add_log_options(parser) + + parser.add_option("--x", + action="store", + dest="x", + type="string", + default = 'r', + help="x value to plot", + metavar=" STRING") + + parser.add_option("--y", + action="store", + dest="y", + type="string", + default = 'T', + help="y value to plot", + metavar=" STRING") + + parser.add_option("--z", + action="store", + dest="z", + type="string", + default = None, + help="z value to plot", + metavar=" STRING") + + parser.add_option("--legend", + action="store_true", + dest="legend", + default = False, + help="add a legend") + + parser.add_option("--elts", + action="store", + type="string", + dest="elts", + default = None, + help="element name") + + parser.add_option("--NSNII", + action="store", + dest="NSNII", + type="float", + default = 1, + help="number of SNII", + metavar=" FLOAT") + + + parser.add_option("--NSNIa", + action="store", + dest="NSNIa", + type="float", + default = 1, + help="number of SNIa", + metavar=" FLOAT") + + + parser.add_option("--NDYIN", + action="store", + dest="NDYIN", + type="float", + default = 1, + help="number of NDYIN", + metavar=" FLOAT") + + parser.add_option("-Z","--Z", + action="store", + dest="Z", + type="float", + default = 0.02, + help="metallicity", + metavar=" FLOAT") + + (options, args) = parser.parse_args() + + + pt.check_files_number(args) + + files = args + + return files,options + + + + + +####################################### +# MakePlot +####################################### + + +def MakePlot(files,opt): + + + # some inits + datas = [] + linest = pt.LineStyles() + + if opt.elts !=None: + opt.elts = string.split(opt.elts,',') + + # read files + for file in files: + + + chemistry.init_chimie(file) + + nelts = chemistry.get_allnelts() + labels = chemistry.get_allelts_labels() + print labels + + + mmax = chemistry.get_Mmax() + mmin = chemistry.get_Mmin() + # mass range + ms = arange(mmin,mmax,1e-12).astype(float32) + + metals = zeros(chemistry.get_nelts(),float) + metals[-1]= opt.Z + + + SolarAbun = chemistry.get_elts_SolarMassAbundances() + + + ################ + # get values + ################ + + if opt.elts==None: + elts = labels + elts.remove('Ej') + elts.remove('Ejnp') + elts.remove('Metals') + + else: + elts = opt.elts + + + + ################################################ + # SNII (no Z dependence) + ################################################ + + colors = pt.Colors(n=len(elts)) + + for elt in elts: + nelt = labels.index(elt) + nFe = labels.index('Fe') + print "(SNII)",elt,nelt + + x = ms + y = zeros(len(x)) + for i,m in enumerate(ms): + + + + EjectedMass = chemistry.SNII_single_mass_ejection(m) + EjectedGasMass = EjectedMass[0] + EjectedEltMass = EjectedMass[nelt] + EjectedFeMass = EjectedMass[nFe] + + + if EjectedGasMass>0: + y[i] = log10( (EjectedEltMass/EjectedFeMass) / (SolarAbun[elt]/SolarAbun['Fe']) + 1.0e-100) + else: + y[i] = -100 + + + + + + # do some cleaning + c = y>-100 + x = compress(c,x) + y = compress(c,y) + + + x = x / (1.989e+33/UnitMass_in_g) # to solar mass + + + xlabel = r"$\rm{Star\,\,Mass [M_{\odot}]}$" + ylabel = r"$\rm{[X/Fe]}$" + #label = "%s (%s)"%(labels[nelt],file) + label = "%s"%(labels[nelt]) + data = pt.DataPoints(x,y,label=label,color=colors.get(),tpe='line',linestyle=linest.current()) + datas.append(data) + + linest.next() + + + # now, plot + for d in datas: + + if d.tpe=='points' or d.tpe=='both': + pt.scatter(d.x,d.y,c=d.color,s=5,linewidths=0,marker='o',vmin=opt.zmin,vmax=opt.zmax) + + if d.tpe=='line' or d.tpe=='both': + pt.plot(d.x,d.y,color=d.color,ls=d.linestyle) + + + pt.text(d.x[-1]+5,d.y[-1],d.label,color=d.color, horizontalalignment='left', verticalalignment='center') + + + if opt.legend: + pt.LegendFromDataPoints(datas,loc='lower right') + + + + + ################################################ + # SNIa (no Z dependence, no mass dependence) + ################################################ + + colors = pt.Colors(n=len(elts)) + + for elt in elts: + nelt = labels.index(elt) + nFe = labels.index('Fe') + label = elt + print "(SNIa)",elt,nelt + + + EjectedMass = chemistry.SNIa_single_mass_ejection(1) + + EjectedGasMass = EjectedMass[0] + EjectedEltMass = EjectedMass[nelt] + EjectedFeMass = EjectedMass[nFe] + + + + + if EjectedGasMass>0: + z = log10( (EjectedEltMass/EjectedFeMass) / (SolarAbun[elt]/SolarAbun['Fe']) + 1.0e-100) + + x = (0,3) + y = (z,z) + color = colors.get() + + pt.plot(x,y,color=color) + pt.text(x[0]-3,y[0],label,color=color, horizontalalignment='right', verticalalignment='center') + + + + data = pt.DataPoints(x,y) + datas.append(data) + + + + + + + # set limits and draw axis + xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log) + + # plot axis + pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log) + pt.xlabel(xlabel,fontsize=pt.labelfont) + pt.ylabel(ylabel,fontsize=pt.labelfont) + pt.grid(False) + + + +######################################################################## +# MAIN +######################################################################## + + + +if __name__ == '__main__': + files,opt = parse_options() + pt.InitPlot(files,opt) + pt.pcolors + MakePlot(files,opt) + pt.EndPlot(files,opt) + + + + + + + + + + + diff --git a/PyChem/shared/Readme b/PyChem/shared/Readme index ffa6cf8..b59d927 100644 --- a/PyChem/shared/Readme +++ b/PyChem/shared/Readme @@ -1,248 +1,255 @@ ############################################################ # used chemical parameter files ############################################################ * for dSph tests (dSph-1324-8new !!! Si->Ba : Fe,Mg,O,Si,Metals) ../scripts/pychem_generate_parameters -p chemistryparameters/chimieparam.yr.py -o chimie.yr.dat # hdf5 version ../scripts/pychem_generate_hdf5_parameters -p chemistryparameters/chimieparam.yr.py -o chimie.yr.h5 * for dSph with new code (compatible with paper simulations, only 4 elts, Fe,Mg,O,Metals) ../scripts/pychem_generate_parameters -p chemistryparameters/chimieparam.yr.4.py -o chimie.yr.4.dat * for dSph2 ("Fe","Mg","O","Ba","Metals") ../scripts/pychem_generate_parameters -p chemistryparameters/chimieparam.yr.2.py -o chimie.yr.2.dat * salpeter IMF ../scripts/pychem_generate_parameters -p chemistryparameters/chimieparam.salpeter.py -o chimie.salpeter.dat * 2 different IMF ../scripts/pychem_generate_parameters -p chemistryparameters/chimieparam.yr.py -o chimie.yr.dat.0 ../scripts/pychem_generate_parameters -p chemistryparameters/chimieparam.yr.py.1 -o chimie.yr.dat.1 * yields from Francois (Mon Apr 18 13:11:44 CEST 2011) ../scripts/pychem_generate_parameters -p chemistryparameters/chimieparam.francois.py -o chimie.francois.dat ./plot_SNII_yields.py --log xy --legend --xmin 0.05 --xmax 100 --ymin 1e-4 --ymax 1 --elts=Ej,Ejnp,Fe,Mg,O chimie.yr.dat chimie.francois.dat ./plot_SNII_interated_yields.py --log xy --legend --xmin 0.05 --xmax 100 --ymin 1e-4 --ymax 1 --elts=Ej,Ejnp,Fe,Mg,O chimie.yr.dat chimie.francois.dat * chimie.francois-met.dat = chimie.francois.dat with Metals from chimie.yr.dat * yields from Francois + SNIa from Francois ../scripts/pychem_generate_parameters -p chemistryparameters/chimieparam.francois-002.py -o chimie.francois-002.dat * yields from Woosley (Tue May 3 15:18:43 CEST 2011) ../scripts/pychem_generate_parameters -p chemistryparameters/chimieparam.woosley.py -o chimie.woosley.dat ( * Element Solar Abundances ( becomes ( Element Solar Mass Abundances (=Mx/MH) ( ( ../scripts/pychem_generate_parameters -p chemistryparameters/chimieparam.yr.py -o chimie.yr-new.dat * truncate SNII tables at 12.9 instead of 10 (Thu Apr 25 15:29:19 CEST 2013) ../scripts/pychem_generate_parameters -p chemistryparameters/chimieparam.yr-25042013.py -o chimie.yr-25042013.dat * AGB ../scripts/pychem_generate_parameters-31052013 -p chemistryparameters/chimieparam.yr-31052013.py -o chimie.yr-31052013.dat ../scripts/pychem_generate_hdf5_parameters-31052013 -p chemistryparameters/chimieparam.yr-31052013.py -o chimie.yr-31052013.h5 * with Cobalt ../scripts/pychem_generate_hdf5_parameters -p chemistryparameters/chimieparam-AGB+OMgSFeZnSrYBaEuCo-20062014.h5.py -o chemistry-AGB+OMgSFeZnSrYBaEuCo-20062014.h5 * with Ca, Ni, Mn, Cr ../scripts/pychem_generate_hdf5_parameters -p chemistryparameters/chimieparam-AGB+OMgSFeZnSrYBaEuCoCaNiMnCr-12022015.h5.py -o chemistry-AGB+OMgSFeZnSrYBaEuCoCaNiMnCr-12022015.h5 * for agora cp chemistryparameters/chimieparam-AGB+OMgSFeZnSrYBaEu-16072013.h5.py chemistryparameters/chimieparam-agora-24092015.h5.py cp tables/SolarAbundances/Grevesse98.txt tables/SolarAbundances/Grevesse98+agora.txt ../scripts/pychem_generate_hdf5_parameters -p chemistryparameters/chimieparam-agora-24092015.h5.py -o chemistry-agora-24092015.h5 * (Thu Dec 22 17:11:21 CET 2016) would like to re-create chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5 cp chemistryparameters/chimieparam-AGB+OMgSFeZnSrYBaEu-16072013.h5.py chemistryparameters/chimieparam-AGB+OMgSFeZnSrYBaEu_TOPH-16072013.h.py ../scripts/pychem_generate_hdf5_parameters -p chemistryparameters/chimieparam-AGB+OMgSFeZnSrYBaEu-16072013.h5.py -o qq.h5 # donne des diff. par rapport à chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5 # meme problème pour reproduir chemistry-agora-24092015.h5 ../scripts/pychem_generate_hdf5_parameters -p chemistryparameters/chimieparam-AGB+OMgSFeZnSrYBaEu_TOPH-16072013.h.py -o chemistry-AGB+OMgSFeZnSrYBaEu_TOPH-16072013.h5 + # change the stars mass (ensure the same num of SNII compared to a Kroupa IMF) + ../scripts/pychem_generate_hdf5_parameters -p chemistryparameters/chimieparam-AGB+OMgSFeZnSrYBaEu_TOPHv2-16072013.h.py -o chemistry-AGB+OMgSFeZnSrYBaEu_TOPHv2-16072013.h5 + + # checks + ../scripts/pychem_sample_IMF_and_properties chemistry-AGB+OMgSFeZnSrYBaEu_TOPHv2-16072013.h5 --mstar 152200 + ../scripts/pychem_plot_SNII_yields --log xy --xmin 0.05 --xmax 1000 --ymin 1e-10 --ymax 1 chemistry-AGB+OMgSFeZnSrYBaEu_TOPHv2-16072013.h5 - - + ../scripts/pychem_plot_lifetime --log xy --xmin 0.05 --xmax 1000 --ymin 1e-1 --ymax 1e10 chemistry-AGB+OMgSFeZnSrYBaEu_TOPHv2-16072013.h5 + + ############################################################ # open and read the chemistry parameter file ############################################################ ../scripts/pychem_read chimie.yr.dat ############################################################ # plot IMF ############################################################ ../scripts/pychem_plot_IMF --log xy chimie.yr.dat chimie.salpeter.dat --legend ############################################################ # sample IMF ############################################################ ../scripts/pychem_sample_IMF chimie.yr.dat --mstar 1e4 ../scripts/pychem_sample_IMF_and_properties chimie.yr.dat --mstar 1e4 ############################################################ # lifetime ############################################################ Plot the stars lifetime as a function of mass and metallicity ../scripts/pychem_plot_lifetime --legend --log xy --xmin 0.05 --xmax 100 --ymin 1e-1 --ymax 1e10 chimie.yr.dat ../scripts/pychem_lifetime chimie.yr.dat.1 ############################################################ # Helium core ############################################################ ./pychem_plot_HeliumCore.py tables/HeliumCore.dat tables/AGBs/OriginalData/karakas-HeliumCore-0.0200.txt ############################################################ # yields with Z ############################################################ # only for DYIN ../scripts/pychem_plot_DYIN_yields_withZ --log xy --xmin 0.05 --xmax 100 --ymin 1e-10 --ymax 1 -Z 1e-1 chemistry.simple.h5 # total (single) ejection ../scripts/pychem_plot_TOTAL_yields_withZ --log xy --legend --xmin 0.05 --xmax 100 --ymin 1e-4 --ymax 1 --NSNII=1 --NSNIa=0 --NDYIN=1 chemistry.simple.h5 ../scripts/pychem_plot_TOTAL_yields_withZ --log xy --legend --xmin 0.05 --xmax 100 --ymin 1e-4 --ymax 1 --NSNII=1 --NSNIa=0 --NDYIN=0 chemistry.simple.h5 ../scripts/pychem_plot_TOTAL_yields_withZ --log xy --legend --xmin 0.05 --xmax 100 --ymin 1e-4 --ymax 1 --NSNII=0 --NSNIa=1 --NDYIN=0 chemistry.simple.h5 ../scripts/pychem_plot_TOTAL_yields_withZ --log xy --legend --xmin 0.05 --xmax 100 --ymin 1e-4 --ymax 1 --NSNII=1 --NSNIa=1 --NDYIN=1 chemistry.simple.h5 # total (integrated) ejection ../scripts/pychem_plot_TOTAL_integrated_yields_withZ --log xy --xmin 0.05 --xmax 100 chemistry.simple.h5 ############################################################ # yields ############################################################ Plot the Mass fraction of ejected elements due to the explotion of one SNII of mass m ../scripts/pychem_plot_SNII_yields --log xy --legend --xmin 0.05 --xmax 100 --ymin 1e-4 --ymax 1 chimie.yr.dat ../scripts/pychem_plot_DYIN_yields --log xy --legend --xmin 0.05 --xmax 100 --ymin 1e-4 --ymax 1 chimie.yr.dat ../scripts/pychem_plot_SNII_yields --log xy --legend --xmin 0.05 --xmax 100 --ymin 1e-4 --ymax 1 --elts=Ej,Ejnp,Fe56,Mg24,O16,Metals chimie.yr.dat # compute in addition, Mg/Fe ratio ../scripts/pychem_plot_SNII_yields --log y --legend --xmin 8 --xmax 30 --ymin 1e-4 --ymax 1 --MgFe --elts=Mg,Fe chimie.yr.dat # total (single) ejection ../scripts/pychem_plot_TOTAL_yields --log xy --legend --xmin 0.05 --xmax 100 --ymin 1e-4 --ymax 1 --NSNII=1 --NSNIa=0 --NDYIN=1 chimie.yr.dat ../scripts/pychem_plot_TOTAL_yields --log xy --legend --xmin 0.05 --xmax 100 --ymin 1e-4 --ymax 1 --NSNII=1 --NSNIa=0 --NDYIN=0 chimie.yr.dat ../scripts/pychem_plot_TOTAL_yields --log xy --legend --xmin 0.05 --xmax 100 --ymin 1e-4 --ymax 1 --NSNII=0 --NSNIa=1 --NDYIN=0 chimie.yr.dat ../scripts/pychem_plot_TOTAL_yields --log xy --legend --xmin 0.05 --xmax 100 --ymin 1e-4 --ymax 1 --NSNII=1 --NSNIa=1 --NDYIN=1 chimie.yr.dat # total (integrated) ejection ../scripts/pychem_plot_TOTAL_integrated_yields --log xy --xmin 0.05 --xmax 100 chimie.yr.dat Plot the integrated yields (!! the zero depends on the zero of the function...) ../scripts/pychem_plot_SNII_interated_yields --log xy --legend --xmin 0.05 --xmax 100 --ymin 1e-4 --ymax 1 chimie.yr.dat Estimation of the energy lost by stars ../scripts/pychem_mass_ejection chimie.yr.dat.1 ############################################################ # chemical evolution ############################################################ # chemical enrichment due to a SSP ../scripts/pychem_SSP_evolution chimie.yr.dat --xmin 0 --xmax 30 --y Mg --timeunit Myr -o Y,ESN ../scripts/pychem_SSP_evolution chimie.yr.dat --xmin 0.001 --xmax 13 --y Mg --timeunit Gyr -o Y,ESN ../scripts/pychem_SSP_evolution chimie.yr.dat --xmin 0.001 --xmax 20 --y Mg --timeunit Gyr -o Y,ESN --log x ../scripts/pychem_SSP_evolution chimie.yr.dat --xmin 1 --xmax 20000 --y Mg --timeunit Myr -o Y,ESN --log x ../scripts/pychem_SSP_evolution chimie.yr.dat --xmin 1 --xmax 20000 --y Mg --timeunit Myr -o D # evolution of the stellar mass ../scripts/pychem_SSP_evolution chimie.yr.dat --xmin 1 --xmax 20000 --y Mg --timeunit Myr -o MSTAR --log x --ymin 0 --ymax 1e-5 --NumberOfTables 2 --DefaultTable 0 ../scripts/pychem_SSP_evolution chimie.yr.dat --xmin 1 --xmax 20000 --y Mg --timeunit Myr -o MSTAR --log x --ymin 0 --ymax 1e-5 --NumberOfTables 2 --DefaultTable 1 (il faudrait un pgm style SSP_evolution.py mais ou l'on plot ce qui nous interesse. ou alors ecrire les output de SSP_evolution.py et le plotter avec un autre pgm...) # using different tables ../scripts/pychem_SSP_evolution chimie.yr.dat --xmin 1 --xmax 20000 --y Mg --timeunit Myr -o Y,ESN --log x --NumberOfTables 2 ../scripts/pychem_SSP_evolution chimie.yr.dat --xmin 1 --xmax 20000 --y Mg --timeunit Myr -o Y,ESN --log x --NumberOfTables 2 --DefaultTable 1 # using the SN trick ../scripts/pychem_SSP_evolution chimie.yr.dat --mstar 1e7 --xmin 0 --xmax 50 --y Mg --timeunit Myr -o NSNII,CumNSNII --mstar 3e4 --dt .1 --discsn ../scripts/pychem_SSP_evolution chimie.yr.dat --mstar 1e7 --xmin 0 --xmax 3000 --y Mg --timeunit Myr -o NSNIa,CumNSNIa --mstar 3e4 --dt .1 --discsn ../scripts/pychem_SSP_discret_evolution2 chimie.yr.dat --mstar 1e7 --xmin 0 --discsn --x Fe --y Mg --timeunit Myr --dt .01 --xmax 3000 -o CumNSNII,CumNSNIa,CumNDYIN,CumEjMass,StellarMass,Y,Yg,D # check SN monte carlo ../scripts/pychem_SSP_discret_evolution2 chimie.yr.dat --mstar 105 --xmin 0 --discsn --x Fe --y Mg --timeunit Myr --dt .01 --xmax 3000 -o CumNSNII,CumNSNIa,CumNDYIN,CumEjMass diff --git a/src/chimie.c b/src/chimie.c index 19c631c..eb54794 100644 --- a/src/chimie.c +++ b/src/chimie.c @@ -1,8676 +1,8676 @@ #include #include #include #include #include #include #include "allvars.h" #include "proto.h" #ifdef CHIMIE #ifdef CHIMIE_FROM_HDF5 #include #include "hdf5io.h" #endif #ifdef PYCHEM #include #include #include #include #include /* **************************************************** these variables are already defined in Gadget (or not needed) **************************************************** */ #define TO_DOUBLE(a) ( (PyArrayObject*) PyArray_CastToType(a, PyArray_DescrFromType(NPY_DOUBLE) ,0) ) #endif /* PYCHEM */ const char *get_filename_ext(const char *filename) { const char *dot = strrchr(filename, '.'); if(!dot || dot == filename) return ""; return dot + 1; } /****************************************************************************************/ /* /* /* /* COMMON CHIMIE PART /* /* /* /****************************************************************************************/ #define MAXPTS 10 #define MAXDATASIZE 200 #define MAXDATAZSIZE 110 /******************************************************** hdf5 reading routines *********************************************************/ #ifdef CHIMIE_FROM_HDF5 #define HEADER_GRP "/Header" #define DATA_GRP "/Data" /******************************************************** Chemistry functions *********************************************************/ struct ChemistryHeaderStruct { char * version; char * author; char * date; }; struct ChemistryDataAttributesStruct { int nelts; char **elts; double *SolarMassAbundances; double MeanWDMass; char *SNIIYieldsFile; char *SNIIHeliumCoreFile; char *DYINYieldsFile; char *DYINHeliumCoreFile; char *SNIaFile; char *SolarAbundancesFile; }; struct ChemistryYieldsTableStruct { int nbins; double min; double step; char *label; double *data; }; struct ChemistryYieldsTable2DStruct { int nx; int ny; double x0; double y0; double dx; double dy; char *label; double **data; }; struct ChemistryLiveTimesStruct { int nx; int ny; double **coeff_z; }; struct ChemistryIMFStruct { int n; double *ms; double *as; char *bs; double Mmin; double Mmax; }; struct ChemistrySNIaStruct { double Mpl; double Mpu; double a; double Mdl1; double Mdu1; double bb1; double Mdl2; double Mdu2; double bb2; struct ElementsDataStruct Metals; }; struct ChemistrySNIIStruct { double Mmin; double Mmax; int npts; int nelts; char **elts; struct ChemistryYieldsTableStruct *table; }; struct ChemistryDYINStruct { double Mmin; double Mmax; int npts; int nelts; char **elts; struct ChemistryYieldsTableStruct *table; struct ChemistryYieldsTable2DStruct *tableZ; int Zflag; int nptZs; }; /*! Read a table containing yields * stored in the dataset named "name" */ struct ChemistryYieldsTableStruct readYieldsTable(hid_t group,char *name) { struct ChemistryYieldsTableStruct table; hid_t dset; herr_t status; table.data = readDatasetAsArrayDouble(group,name); /* read attributes */ dset = H5Dopen( group , name, H5P_DEFAULT); table.nbins = readAttributeAsInt(dset,"nbins"); table.min = readAttributeAsDouble(dset,"min"); table.step = readAttributeAsDouble(dset,"step"); table.label = readAttributeAsString(dset,"label"); //printYieldsTable(table); status = H5Dclose(dset); return table; } /*! Read a table containing 2d yields * stored in the dataset named "name" */ struct ChemistryYieldsTable2DStruct readYieldsTable2D(hid_t group,char *name) { struct ChemistryYieldsTable2DStruct table; hid_t dset; herr_t status; table.data = readDatasetAsArray2DDouble_v0(group,name); /* read attributes */ dset = H5Dopen( group , name, H5P_DEFAULT); table.nx = readAttributeAsInt(dset,"nx"); table.ny = readAttributeAsInt(dset,"ny"); table.x0 = readAttributeAsDouble(dset,"x0"); table.y0 = readAttributeAsDouble(dset,"y0"); table.dx = readAttributeAsDouble(dset,"dx"); table.dy = readAttributeAsDouble(dset,"dy"); table.label = readAttributeAsString(dset,"label"); //printYieldsTable(table); status = H5Dclose(dset); return table; } /*! Read the attributes linked to the data group */ int readDataAttributes(hid_t table,struct ChemistryDataAttributesStruct *Param) { hid_t group; herr_t status; group = H5Gopen(table, DATA_GRP,H5P_DEFAULT); Param->nelts = readAttributeAsInt(group,"nelts"); Param->elts = readAttributeAsArrayString(group,"elts"); Param->SolarMassAbundances = readAttributeAsArrayDouble(group,"SolarMassAbundances"); Param->MeanWDMass = readAttributeAsDouble(group,"MeanWDMass"); Param->elts = readAttributeAsArrayString(group,"elts"); Param->SNIIYieldsFile = readAttributeAsString(group,"SNIIYieldsFile"); Param->SNIIHeliumCoreFile = readAttributeAsString(group,"SNIIHeliumCoreFile"); Param->DYINYieldsFile = readAttributeAsString(group,"DYINYieldsFile"); Param->DYINHeliumCoreFile = readAttributeAsString(group,"DYINHeliumCoreFile"); Param->SNIaFile = readAttributeAsString(group,"SNIaFile"); Param->SolarAbundancesFile = readAttributeAsString(group,"SolarAbundancesFile"); status = H5Gclose (group); return 0; } /*! Print the attributes linked to the data group */ int printDataAttributes(struct ChemistryDataAttributesStruct Parameters) { int i; printf("\n"); printf("Data attribute content:\n\n"); printf("\t nelts = %d\n",Parameters.nelts); printf("\t elts = "); for (i=0; iversion = readAttributeAsString(group,"version"); Header->author = readAttributeAsString(group,"author"); Header->date = readAttributeAsString(group,"date"); status = H5Gclose (group); return 0; } /*! Print the attributes linked to the header group */ int printHeader(struct ChemistryHeaderStruct Header) { printf("\n"); printf("Header content:\n\n"); printf("\t version = %s\n",Header.version); printf("\t author = %s\n",Header.author); printf("\t date = %s\n",Header.date); printf("\n"); return 0; } /*! Print the content of a yield table */ int printYieldsTable(struct ChemistryYieldsTableStruct table,char * space) { int i; printf("%s%s:\n\n",space,table.label); printf("%s\t nbins = %d\n",space,table.nbins); printf("%s\t min = %g\n",space,table.min); printf("%s\t step = %g\n",space,table.step); printf("%s\t label = %s\n",space,table.label); printf("\n"); for(i=0;i<3;i++) printf("%s\t [%d] %g\n",space,i,table.data[i]); printf("%s\t ...\n",space); for(i=table.nbins-3;icoeff_z = readDatasetAsArray2DDouble_v0(subgroup,"coeff_z"); dset = H5Dopen( subgroup , "coeff_z", H5P_DEFAULT); livetimes->nx = readAttributeAsInt(dset,"nx"); livetimes->ny = readAttributeAsInt(dset,"ny"); status = H5Dclose (dset); status = H5Gclose (subgroup); status = H5Gclose (group); return 0; } /*! Print the content of an LiveTimes struct */ int printLiveTimes(struct ChemistryLiveTimesStruct livetimes) { int i,j; printf("Output from LiveTimes:\n\n"); printf("\t coeff_z = \n"); for (i=0; iMmin = readAttributeAsDouble(subgroup,"Mmin"); imf->Mmax = readAttributeAsDouble(subgroup,"Mmax"); imf->n = readAttributeAsInt(subgroup,"n"); imf->ms = readAttributeAsArrayDouble(subgroup,"ms"); imf->as = readAttributeAsArrayDouble(subgroup,"as"); status = H5Gclose (subgroup); status = H5Gclose (group); return 0; } /*! Print the content of an IMF */ int printIMF(struct ChemistryIMFStruct imf) { int i; printf("Output from IMF:\n\n"); printf("\t n = %d\n",imf.n); printf("\t Mmin = %g\n",imf.Mmin); printf("\t Mmax = %g\n",imf.Mmax); printf("\t as = "); for (i=0; iMpl = readAttributeAsDouble(subgroup,"Mpl"); snia->Mpu = readAttributeAsDouble(subgroup,"Mpu"); snia->a = readAttributeAsDouble(subgroup,"a"); snia->Mdl1 = readAttributeAsDouble(subgroup,"Mdl1"); snia->Mdu1 = readAttributeAsDouble(subgroup,"Mdu1"); snia->bb1 = readAttributeAsDouble(subgroup,"bb1"); snia->Mdl2 = readAttributeAsDouble(subgroup,"Mdl2"); snia->Mdu2 = readAttributeAsDouble(subgroup,"Mdu2"); snia->bb2 = readAttributeAsDouble(subgroup,"bb2"); snia->Metals = readGroupAsElementsData(subgroup,"Metals"); status = H5Gclose (subgroup); status = H5Gclose (group); return 0; } /*! Print the content of an SNIa struct */ int printSNIa(struct ChemistrySNIaStruct snia) { int i; printf("Output from SNIa:\n\n"); printf("\t Mpl = %g\n",snia.Mpl ); printf("\t Mpu = %g\n",snia.Mpu ); printf("\t a = %g\n",snia.a ); printf("\t Mdl1 = %g\n",snia.Mdl1); printf("\t Mdu1 = %g\n",snia.Mdu1); printf("\t bb1 = %g\n",snia.bb1 ); printf("\t Mdl2 = %g\n",snia.Mdl2); printf("\t Mdu2 = %g\n",snia.Mdu2); printf("\t bb2 = %g\n",snia.bb2 ); printf("\n"); printf("\t Metals:\n\n"); for (i=0; iMmin = readAttributeAsDouble(subgroup,"Mmin"); snii->Mmax = readAttributeAsDouble(subgroup,"Mmax"); snii->npts = readAttributeAsInt(subgroup,"npts"); snii->nelts = readAttributeAsInt(subgroup,"nelts"); snii->elts = readAttributeAsArrayString(subgroup,"elts"); snii->table = malloc( snii->nelts * sizeof(struct ChemistryYieldsTableStruct) ); /* read the yields */ for(i=0;inelts;i++) snii->table[i]=readYieldsTable(subgroup,snii->elts[i]); status = H5Gclose (subgroup); status = H5Gclose (group); return 0; } /*! Print the content of an SNII struct */ int printSNII(struct ChemistrySNIIStruct snii) { int i; printf("Output from SNII:\n\n"); printf("\t Mmin = %g\n",snii.Mmin); printf("\t Mmax = %g\n",snii.Mmax); printf("\t npts = %d\n",snii.npts); printf("\t nelts = %d\n",snii.nelts); printf("\t elts = "); for (i=0; iMmin = readAttributeAsDouble(subgroup,"Mmin"); dyin->Mmax = readAttributeAsDouble(subgroup,"Mmax"); dyin->npts = readAttributeAsInt(subgroup,"npts"); dyin->nelts = readAttributeAsInt(subgroup,"nelts"); dyin->elts = readAttributeAsArrayString(subgroup,"elts"); dyin->table = malloc( dyin->nelts * sizeof(struct ChemistryYieldsTableStruct) ); /* read the yields */ for(i=0;inelts;i++) dyin->table[i]=readYieldsTable(subgroup,dyin->elts[i]); dyin->Zflag = readAttributeAsInt(subgroup,"Zflag"); if (dyin->Zflag==1) { /* read data with Z dependences */ subsubgroup = H5Gopen(subgroup,"MetallicityDependent",H5P_DEFAULT); dyin->nptZs = readAttributeAsInt(subsubgroup,"nptZs"); dyin->tableZ = malloc( dyin->nelts * sizeof(struct ChemistryYieldsTable2DStruct) ); /* read the yields */ for(i=0;inelts;i++) dyin->tableZ[i]=readYieldsTable2D(subsubgroup,dyin->elts[i]); status = H5Gclose (subsubgroup); /* end read data with Z dependences */ } status = H5Gclose (subgroup); status = H5Gclose (group); return 0; } /*! Print the content of an DYIN struct */ int printDYIN(struct ChemistryDYINStruct dyin) { int i; printf("Output from DYIN:\n\n"); printf("\t Mmin = %g\n",dyin.Mmin); printf("\t Mmax = %g\n",dyin.Mmax); printf("\t npts = %d\n",dyin.npts); printf("\t nelts = %d\n",dyin.nelts); if (dyin.Zflag==1) printf("\t nptZs = %d\n",dyin.nptZs); printf("\t elts = "); for (i=0; i=All.ChimieNumberOfParameterFiles) { printf("\n set_table : i>= %d !!!\n\n",All.ChimieNumberOfParameterFiles); endrun(88809); } else { Cp = &Cps[i]; Elt = Elts[i]; MassFracSNII = MassFracSNIIs[i]; /* all this is useless, no ?*/ MassFracSNIa = MassFracSNIas[i]; MassFracDYIN = MassFracDYINs[i]; SingleMassFracSNII = SingleMassFracSNIIs[i]; SingleMassFracSNIa = SingleMassFracSNIas[i]; SingleMassFracDYIN = SingleMassFracDYINs[i]; EjectedMass = EjectedMasss[i]; SingleEjectedMass = SingleEjectedMasss[i]; } } /*! Read the chemistry table (hdf5 format) */ #ifdef CHIMIE_FROM_HDF5 void read_chimie_h5(char * filename,int it) { hid_t table; herr_t status; struct ChemistryHeaderStruct ChemistryHeader; struct ChemistryDataAttributesStruct ChemistryBasicParameters; struct ChemistryLiveTimesStruct ChemistryLiveTimes; struct ChemistryIMFStruct ChemistryIMF; struct ChemistrySNIaStruct ChemistrySNIa; struct ChemistrySNIIStruct ChemistrySNII; struct ChemistryDYINStruct ChemistryDYIN; table = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT); readHeader(table,&ChemistryHeader); if (verbose && ThisTask==0) printHeader(ChemistryHeader); readDataAttributes(table,&ChemistryBasicParameters); if (verbose && ThisTask==0) printDataAttributes(ChemistryBasicParameters); readLiveTimes(table,&ChemistryLiveTimes); if (verbose && ThisTask==0) printLiveTimes(ChemistryLiveTimes); readIMF(table,&ChemistryIMF); if (verbose && ThisTask==0) printIMF(ChemistryIMF); readSNIa(table,&ChemistrySNIa); if (verbose && ThisTask==0) printSNIa(ChemistrySNIa); readSNII(table,&ChemistrySNII); if (verbose && ThisTask==0) printSNII(ChemistrySNII); readDYIN(table,&ChemistryDYIN); if (verbose && ThisTask==0) printDYIN(ChemistryDYIN); status = H5Fclose(table); Cps[it].MassMinForEjecta = MAX_REAL_NUMBER; /* convert to chimie struct */ int i,j,k; /* Livetimes */ for(i=0;i<3;i++) for(j=0;j<3;j++) Cps[it].coeff_z[i][j] = ChemistryLiveTimes.coeff_z[i][j]; /* IMF Parameters */ Cps[it].Mmin = ChemistryIMF.Mmin; Cps[it].Mmax = ChemistryIMF.Mmax; Cps[it].n = ChemistryIMF.n; if (Cps[it].n>0) for (i=0;i MAXDATASIZE = %d !!!\n\n",Cps[it].nptsSNII,MAXDATASIZE); printf("\n Cps[it].nptsDYIN = %d > MAXDATASIZE = %d !!!\n\n",Cps[it].nptsDYIN,MAXDATASIZE); endrun(88800); } /* allocate memory */ MassFracSNIIs[it] = malloc((Cps[it].nelts+2) * sizeof(double)); /* really needed ? */ MassFracSNIas[it] = malloc((Cps[it].nelts+2) * sizeof(double)); MassFracDYINs[it] = malloc((Cps[it].nelts+2) * sizeof(double)); EjectedMasss[it] = malloc((Cps[it].nelts+2) * sizeof(double)); SingleMassFracSNIIs[it] = malloc((Cps[it].nelts+2) * sizeof(double)); SingleMassFracSNIas[it] = malloc((Cps[it].nelts+2) * sizeof(double)); SingleMassFracDYINs[it] = malloc((Cps[it].nelts+2) * sizeof(double)); SingleEjectedMasss[it] = malloc((Cps[it].nelts+2) * sizeof(double)); /**/ /* injected metals SNII */ /**/ for (i=0;i0) for (i=0;i MAXDATASIZE = %d !!!\n\n",Cps[it].nptsSNII,MAXDATASIZE); endrun(88800); } /* allocate memory */ MassFracSNIIs[it] = malloc((Cps[it].nelts+2) * sizeof(double)); /* really needed ? */ MassFracSNIas[it] = malloc((Cps[it].nelts+2) * sizeof(double)); MassFracDYINs[it] = malloc((Cps[it].nelts+2) * sizeof(double)); EjectedMasss[it] = malloc((Cps[it].nelts+2) * sizeof(double)); SingleMassFracSNIIs[it] = malloc((Cps[it].nelts+2) * sizeof(double)); SingleMassFracSNIas[it] = malloc((Cps[it].nelts+2) * sizeof(double)); SingleMassFracDYINs[it] = malloc((Cps[it].nelts+2) * sizeof(double)); SingleEjectedMasss[it] = malloc((Cps[it].nelts+2) * sizeof(double)); /* injected metals */ for (i=0;in; /* convert m in msol */ m = m*All.CMUtoMsol; if (n==0) return Cp->bs[0]* pow(m,Cp->as[0]); else { for (i=0;ims[i]) return Cp->bs[i]* pow(m,Cp->as[i]); return Cp->bs[n]* pow(m,Cp->as[n]); } } /*! This function returns the mass fraction between m1 and m2 * per mass unit, using the current IMF */ static double get_imf_M(double m1, double m2) { int i; int n; double p; double integral=0; double mmin,mmax; n = Cp->n; /* convert m in msol */ m1 = m1*All.CMUtoMsol; m2 = m2*All.CMUtoMsol; if (n==0) { p = Cp->as[0]+1; integral = (Cp->bs[0]/p) * ( pow(m2,p) - pow(m1,p) ); //printf("--> %g %g %g %g int=%g\n",m1,m2,pow(m2,p), pow(m1,p),integral); } else { integral = 0; /* first */ if (m1ms[0]) { mmin = m1; mmax = dmin(Cp->ms[0],m2); p = Cp->as[0] + 1; integral += (Cp->bs[0]/p) * ( pow(mmax,p) - pow(mmin,p) ); } /* last */ if (m2>Cp->ms[n-1]) { mmin = dmax(Cp->ms[n-1],m1); mmax = m2; p = Cp->as[n] + 1; integral += (Cp->bs[n]/p) * ( pow(mmax,p) - pow(mmin,p) ); } /* loop over other segments */ for (i=0;ims[i ],m1); mmax = dmin(Cp->ms[i+1],m2); if (mminas[i+1] + 1; integral += (Cp->bs[i+1]/p) * ( pow(mmax,p) - pow(mmin,p) ); } } } return integral; } /*! This function returns the number fraction between m1 and m2 * per mass unit, using the current IMF */ static double get_imf_N(double m1, double m2) { int i; int n; double p; double integral=0; double mmin,mmax; n = Cp->n; /* convert m in msol */ m1 = m1*All.CMUtoMsol; m2 = m2*All.CMUtoMsol; if (n==0) { p = Cp->as[0]; integral = (Cp->bs[0]/p) * ( pow(m2,p) - pow(m1,p) ); } else { integral = 0; /* first */ if (m1ms[0]) { mmin = m1; mmax = dmin(Cp->ms[0],m2); p = Cp->as[0]; integral += (Cp->bs[0]/p) * ( pow(mmax,p) - pow(mmin,p) ); } /* last */ if (m2>Cp->ms[n-1]) { mmin = dmax(Cp->ms[n-1],m1); mmax = m2; p = Cp->as[n]; integral += (Cp->bs[n]/p) * ( pow(mmax,p) - pow(mmin,p) ); } /* loop over other segments */ for (i=0;ims[i ],m1); mmax = dmin(Cp->ms[i+1],m2); if (mminas[i+1]; integral += (Cp->bs[i+1]/p) * ( pow(mmax,p) - pow(mmin,p) ); } } } /* convert into mass unit mass unit */ integral = integral *All.CMUtoMsol; return integral; } /*! Sample the imf using monte carlo approach */ static double imf_sampling() { int i; int n; double m; double f; double pmin,pmax; n = Cp->n; /* init random */ //srandom(irand); f = (double)random()/(double)RAND_MAX; if (n==0) { pmin = pow(Cp->Mmin,Cp->as[0]); pmax = pow(Cp->Mmax,Cp->as[0]); m = pow(f*(pmax - pmin) + pmin ,1./Cp->as[0]); return m* All.MsoltoCMU; } else { if (ffs[0]) { pmin = pow(Cp->Mmin ,Cp->as[0]); m = pow(Cp->imf_Ntot*Cp->as[0]/Cp->bs[0]* (f-0) + pmin ,1./Cp->as[0]); return m* All.MsoltoCMU; } for (i=0;ifs[i+1]) { pmin = pow(Cp->ms[i] ,Cp->as[i+1]); m = pow(Cp->imf_Ntot*Cp->as[i+1]/Cp->bs[i+1]* (f-Cp->fs[i]) + pmin ,1./Cp->as[i+1]); return m* All.MsoltoCMU; } } /* last portion */ pmin = pow(Cp->ms[n-1] ,Cp->as[n]); m = pow(Cp->imf_Ntot*Cp->as[n]/Cp->bs[n]* (f-Cp->fs[n-1]) + pmin ,1./Cp->as[n]); return m* All.MsoltoCMU; } } static double imf_sampling_from_random(double f) { int i; int n; double m; double pmin,pmax; n = Cp->n; if (n==0) { pmin = pow(Cp->Mmin,Cp->as[0]); pmax = pow(Cp->Mmax,Cp->as[0]); m = pow(f*(pmax - pmin) + pmin ,1./Cp->as[0]); return m* All.MsoltoCMU; } else { if (ffs[0]) { pmin = pow(Cp->Mmin ,Cp->as[0]); m = pow(Cp->imf_Ntot*Cp->as[0]/Cp->bs[0]* (f-0) + pmin ,1./Cp->as[0]); return m* All.MsoltoCMU; } for (i=0;ifs[i+1]) { pmin = pow(Cp->ms[i] ,Cp->as[i+1]); m = pow(Cp->imf_Ntot*Cp->as[i+1]/Cp->bs[i+1]* (f-Cp->fs[i]) + pmin ,1./Cp->as[i+1]); return m* All.MsoltoCMU; } } /* last portion */ pmin = pow(Cp->ms[n-1] ,Cp->as[n]); m = pow(Cp->imf_Ntot*Cp->as[n]/Cp->bs[n]* (f-Cp->fs[n-1]) + pmin ,1./Cp->as[n]); return m* All.MsoltoCMU; } } /*! This function initializes the imf parameters defined in the chemistry file */ void init_imf(void) { float integral = 0; float p; float cte; int i,n; double mmin,mmax; n = Cp->n; if (n==0) { p = Cp->as[0]+1; integral = integral + ( pow(Cp->Mmax,p)-pow(Cp->Mmin,p))/(p) ; Cp->bs[0] = 1./integral ; } else { cte = 1.0; if (Cp->Mmin < Cp->ms[0]) { p = Cp->as[0]+1; integral = integral + (pow(Cp->ms[0],p) - pow(Cp->Mmin,p))/p; } for (i=0;ims[i],( Cp->as[i] - Cp->as[i+1] )); p = Cp->as[i+1]+1; integral = integral + cte*(pow(Cp->ms[i+1],p) - pow(Cp->ms[i],p))/p; } if (Cp->Mmax > Cp->ms[-1]) { cte = cte* pow( Cp->ms[n-1] , ( Cp->as[n-1] - Cp->as[n] ) ); p = Cp->as[n]+1; integral = integral + cte*(pow(Cp->Mmax,p) - pow(Cp->ms[n-1],p))/p; } /* compute all b */ Cp->bs[0] = 1./integral; for (i=0;ibs[i+1] = Cp->bs[i] * pow( Cp->ms[i],( Cp->as[i] - Cp->as[i+1] )); } } //if (verbose && ThisTask==0) // { // printf("-- bs -- \n"); // for (i=0;ibs[i]); // printf("\n"); // } mmin = Cp->Mmin *All.MsoltoCMU; /* in CMU */ mmax = Cp->Mmax *All.MsoltoCMU; /* in CMU */ Cp->imf_Ntot = get_imf_N(mmin,mmax) *All.MsoltoCMU; /* in CMU */ /* init fs : mass fraction at ms */ if (n>0) { for (i=0;ims[i] *All.MsoltoCMU; /* in CMU */ Cp->fs[i] = All.MsoltoCMU*get_imf_N(mmin,mmax)/Cp->imf_Ntot; } } } #ifdef CHIMIE_OPTIMAL_SAMPLING /* integral of m\xi with the kroupa normalisation ( in k) input : mass in Msol */ double imf_int_mxi(double left,double right, double k) { return get_imf_M(left*All.MsoltoCMU,right*All.MsoltoCMU)/(Cp->bs[0])*k; } /* integral of \xi with the kroupa normalisation ( in k) input : mass in Msol */ double imf_int_xi(double left,double right, double k) { return get_imf_N(left*All.MsoltoCMU,right*All.MsoltoCMU)/ (Cp->bs[0]*All.CMUtoMsol) * k; } /* compute the normalisation and the maximum stellar mass i.e, Cp->k and Cp->m_max input : mass in Msol output : mass in Msol */ double optimal_init_norm(double Mecl) { double a,b,c; double mb; double Mmax; Mmax = Cp->Mmax; Cp->k = 1.; /* kroupa normalisation */ //we use 150Msun for the maximum possible stellar mass (TODO: put it as a parameter) a = Cp->Mmin; c = Mmax; b = (c+a)/2.; //solve the two equations to find m_max and k iteratively while(((c/b)-(a/b)) > 0.00001) { mb = imf_int_mxi(Cp->Mmin,b,Cp->k)/imf_int_xi(b,Mmax,Cp->k)+b; if(mb < Mecl) a = b; else c = b; b = (c+a)/2; } //store the results Cp->m_max = b; Cp->k = (Mecl-Cp->m_max) / imf_int_mxi(Cp->Mmin,Cp->m_max,Cp->k); return Cp->m_max; } /* compute the normalisation and the maximum stellar mass for one single stellar particle output in StP[m].OptIMF_CurrentMass, in Msol */ double init_optimal_imf_sampling(int i) { double a,b,c; double mb; double Mmax; double Mecl; int m; m = P[i].StPIdx; Mmax = Cp->Mmax; Mecl = StP[m].InitialMass*All.CMUtoMsol; /* to Msol */ StP[m].OptIMF_k = 1.; /* kroupa normalisation */ StP[m].OptIMF_N_WD = 0; //we use 150Msun for the maximum possible stellar mass (TODO: put it as a parameter) a = Cp->Mmin; c = Mmax; b = (c+a)/2.; //solve the two equations to find m_max and k iteratively while(((c/b)-(a/b)) > 0.00001) { mb = imf_int_mxi(Cp->Mmin,b,StP[m].OptIMF_k)/imf_int_xi(b,Mmax,StP[m].OptIMF_k)+b; if(mb < Mecl) a = b; else c = b; b = (c+a)/2; } //store the results StP[m].OptIMF_m_max = b; StP[m].OptIMF_k = (Mecl-StP[m].OptIMF_m_max) / imf_int_mxi(Cp->Mmin,StP[m].OptIMF_m_max,StP[m].OptIMF_k); StP[m].OptIMF_CurrentMass = StP[m].OptIMF_m_max; return StP[m].OptIMF_m_max; } /* compute the next mass given the previous input : mass in Msol output : mass in Msol */ double optimal_get_next_mass(double m){ double a,b,c,mb; const double err = 1.e-8; a = Cp->Mmin; c = m; b = (c+a)/2; mb = imf_int_mxi(b,m,Cp->k); while(fabs((mb-b)/mb) > err) { mb = imf_int_mxi(b,m,Cp->k); if(mb < b) c = b; else a = b; b = (c+a)/2; } return b; } /* comute the next mass given the previous, according to the imf of a given particle output : mass in Msol */ double optimal_get_next_mass_for_one_particle(int j){ double a,b,c,mb; const double err = 1.e-8; double m; m = StP[j].OptIMF_CurrentMass; a = Cp->Mmin; c = m; b = (c+a)/2; mb = imf_int_mxi(b,m,StP[j].OptIMF_k); while(fabs((mb-b)/mb) > err) { mb = imf_int_mxi(b,m,StP[j].OptIMF_k); if(mb < b) c = b; else a = b; b = (c+a)/2; } return b; } /* stop the imf loop input : mass in Msol */ int optimal_stop_loop(double m) { int bool; if(imf_int_mxi(Cp->Mmin,m,Cp->k) < Cp->Mmin) bool = 1; else bool = 0; return bool; } /* find m1 such as the stellar mass, according to the IMF, between m1 and m2 is mp. input : masses in Msol */ double optimal_get_m1_from_m2(double m2, double mp) { double a,b,c; double m1,mb; a = Cp->Mmin; c = m2; b = (c+a)/2.; //solve the two equations to find m_max and k iteratively while(((c/b)-(a/b)) > 0.00001) { mb = imf_int_mxi(b,m2,Cp->k); if(mb < mp) c = b; else a = b; b = (c+a)/2; } //store the results m1 = b; return m1; } float optimal_get_k_normalisation_factor() { return Cp->k; } #endif /* CHIMIE_OPTIMAL_SAMPLING */ /*! This function initializes the chemistry parameters */ void init_chimie(void) { int i,nf; double u_lt; double UnitLength_in_kpc; double UnitMass_in_Msol; char filename[500]; char ext[100]; /* check some flags */ #ifndef COSMICTIME if (All.ComovingIntegrationOn) { if(ThisTask == 0) printf("Code wasn't compiled with COSMICTIME support enabled!\n"); endrun(-88800); } #endif #ifdef CHIMIE_AGB if(ThisTask == 0) printf("AGB is enabled\n"); #else if(ThisTask == 0) printf("AGB is disabled\n"); #endif UnitLength_in_kpc = All.UnitLength_in_cm / KPC_IN_CM; UnitMass_in_Msol = All.UnitMass_in_g / SOLAR_MASS; //u_lt = -log10( 4.7287e11*sqrt(pow(UnitLength_in_kpc,3)/UnitMass_in_Msol)); /*Sat Dec 25 23:27:10 CET 2010 */ u_lt = -log10(All.UnitTime_in_Megayears*1e6); allocate_chimie(); for (nf=0;nfcoeff_z[2][2] = Cp->coeff_z[2][2] + u_lt; for (i=0;i<3;i++) Cp->coeff_z[1][i] = Cp->coeff_z[1][i]/2.0; /* init imf parameters */ init_imf(); /* init SNII parameters */ //if (Cp->n==0) // { // //Cp->SNII_cte[0] = Cp->bs[0]/Cp->as[0]; // Cp->SNII_cte = Cp->bs[0]/Cp->as[0]; // Cp->SNII_a = Cp->as[0]; // } //else // { // //for (i=0;in+1;i++) /* if multiple power law in the SNII mass range */ // // Cp->SNII_cte[i] = Cp->bs[i]/Cp->as[i]; // Cp->SNII_cte = Cp->bs[Cp->n]/Cp->as[Cp->n]; // Cp->SNII_a = Cp->as[Cp->n]; // } /* init DYIN parameters */ //if (Cp->n==0) // { // //Cp->DYIN_cte[0] = Cp->bs[0]/Cp->as[0]; // Cp->DYIN_cte = Cp->bs[0]/Cp->as[0]; // Cp->DYIN_a = Cp->as[0]; // } //else // { // //for (i=0;in+1;i++) /* if multiple power law in the DYIN mass range */ // // Cp->DYIN_cte[i] = Cp->bs[i]/Cp->as[i]; // Cp->DYIN_cte = Cp->bs[Cp->n]/Cp->as[Cp->n]; // Cp->DYIN_a = Cp->as[Cp->n]; // } /* init SNIa parameters */ Cp->SNIa_a1 = Cp->SNIa_a; Cp->SNIa_b1 = (Cp->SNIa_a1+1)/(pow(Cp->SNIa_Mdu1,Cp->SNIa_a1+1)-pow(Cp->SNIa_Mdl1,Cp->SNIa_a1+1)); Cp->SNIa_cte1 = Cp->SNIa_b1/Cp->SNIa_a1; Cp->SNIa_a2 = Cp->SNIa_a; Cp->SNIa_b2 = (Cp->SNIa_a2+1)/(pow(Cp->SNIa_Mdu2,Cp->SNIa_a2+1)-pow(Cp->SNIa_Mdl2,Cp->SNIa_a2+1)); Cp->SNIa_cte2 = Cp->SNIa_b2/Cp->SNIa_a2; /* init SNIa parameters */ if (Cp->n==0) { Cp->SNIa_cte = Cp->bs[0]/Cp->as[0]; Cp->SNIa_a = Cp->as[0]; } else { Cp->SNIa_cte = Cp->bs[Cp->n]/Cp->as[Cp->n]; Cp->SNIa_a = Cp->as[Cp->n]; } //for (i=0;inelts+2;i++) // Elt[i].MminSNII = log10(Elt[i].MminSNII); // //for (i=0;inelts+2;i++) // Elt[i].MminDYIN = log10(Elt[i].MminDYIN); /* output info */ //if (verbose && ThisTask==0) // { // //printf("-- SNII_cte -- \n"); // //for (i=0;in+1;i++) // // printf("%g ",Cp->SNII_cte[i]); // //printf("%g ",Cp->SNII_cte); // // printf("\n"); //} /* output info */ //if (verbose && ThisTask==0) // { // printf("-- DYIN_cte -- \n"); // //for (i=0;in+1;i++) // // printf("%g ",Cp->DYIN_cte[i]); // //printf("%g ",Cp->DYIN_cte); // // printf("\n"); // } /* check that the masses are higher than the last IMF elbow */ if (Cp->n>0) { if (Cp->SNIa_Mpl < Cp->ms[Cp->n-1]) { printf("\nSNIa_Mpl = %g < ms[n-1] = %g !!!\n\n",Cp->SNIa_Mpl,Cp->ms[Cp->n-1]); //printf("This is not supported by the current implementation !!!\n"); //endrun(88801); printf("We assume that there will be no SNIa !!!\n"); Cp->SNIa_Mpl = -1; Cp->SNIa_Mpu = -1; } if (Cp->SNIa_Mpu < Cp->ms[Cp->n-1]) { printf("\nSNIa_Mpu = %g < ms[n-1] = %g !!!\n\n",Cp->SNIa_Mpu,Cp->ms[Cp->n-1]); //printf("This is not supported by the current implementation !!!\n"); //endrun(88802); printf("We assume that there will be no SNIa !!!\n"); Cp->SNIa_Mpl = -1; Cp->SNIa_Mpu = -1; } if (Cp->SNII_Mmin < Cp->ms[Cp->n-1]) { printf("\nSNII_Mmin = %g < ms[n-1] = %g !!!\n\n",Cp->SNII_Mmin,Cp->ms[Cp->n-1]); printf("This is not supported by the current implementation !!!\n"); endrun(88803); } if (Cp->SNII_Mmax < Cp->ms[Cp->n-1]) { printf("\nSNII_Mmax = %g < ms[n-1] = %g !!!\n\n",Cp->SNII_Mmax,Cp->ms[Cp->n-1]); printf("This is not supported by the current implementation !!!\n"); endrun(88804); } if (Cp->DYIN_Mmin < Cp->ms[Cp->n-1]) { printf("\nDYIN_Mmin = %g < ms[n-1] = %g !!!\n\n",Cp->DYIN_Mmin,Cp->ms[Cp->n-1]); printf("This is not supported by the current implementation !!!\n"); endrun(88805); } if (Cp->DYIN_Mmax < Cp->ms[Cp->n-1]) { printf("\nDYIN_Mmax = %g < ms[n-1] = %g !!!\n\n",Cp->DYIN_Mmax,Cp->ms[Cp->n-1]); printf("This is not supported by the current implementation !!!\n"); endrun(88806); } } /* number of WD per solar mass */ Cp->SNIa_NWDN = Cp->SNIa_cte * ( pow(Cp->SNIa_Mpu,Cp->SNIa_a)-pow(Cp->SNIa_Mpl,Cp->SNIa_a) ); } /* init short cuts to some elements */ FE = get_ElementIndex("Fe"); METALS = get_ElementIndex("Metals"); /* write chime stat header */ #ifdef CHIMIE_STATS int k; fprintf(FdChimieStatsSNs,"# ThisTask TimeStep Time ID Mass NSNII NSNIa MassSN x y z "); fprintf(FdChimieStatsSNs,"EjectedGasMass "); for (k=0;knelts); for(i=2;inelts+2;i++) printf("%s ",&Elt[i].label); printf("\n"); } /* check number of elements */ if (NELEMENTS != Cp->nelts) { printf("(Taks=%d) NELEMENTS (=%d) != Cp->nelts (=%d) : please check !!!\n\n",ThisTask,NELEMENTS,Cp->nelts); endrun(88807); } /* check that iron is the first element */ //if ((strcmp("Fe",Elt[2].label))!=0) // { // printf("(Taks=%d) first element (=%s) is not %s !!!\n\n",ThisTask,Elt[2].label,FIRST_ELEMENT); // endrun(88808); // } } /*! This function print some info on the chimie parameters */ void info(int it) { int i,j; printf("\nTable %d\n\n",it); printf("%g %g %g\n", Cps[it].coeff_z[0][0],Cps[it].coeff_z[0][1],Cps[it].coeff_z[0][2]); printf("%g %g %g\n", Cps[it].coeff_z[1][0],Cps[it].coeff_z[1][1],Cps[it].coeff_z[1][2]); printf("%g %g %g\n", Cps[it].coeff_z[2][0],Cps[it].coeff_z[2][1],Cps[it].coeff_z[2][2]); printf("\n"); printf("\nIMF\n"); printf("%g %g\n",Cps[it].Mmin,Cps[it].Mmax); printf("%d\n",Cps[it].n); for (i=0;i %g %g\n",Elts[it][i].MminSNII,Elts[it][i].StepSNII); for (j=0;j %g %g\n",Elts[it][i].MminDYIN,Elts[it][i].StepDYIN); for (j=0;jnelts; } /*! Return the solar mass abundance of elt i */ float get_SolarMassAbundance(i) { return Elt[i+2].SolarMassAbundance; } /*! Return the label of element i */ char* get_Element(i) { return Elt[i+2].label; } int get_ElementIndex(char *elt) { /* example: get_SolarMassAbundance(get_ElementIndex("Fe")); */ int i=-1; for(i=0;inelts;i++) if( strcmp(Elt[i+2].label, elt) == 0 ) return i; printf("element %s is unknown\n",elt); endrun(777123); return i; } /*! Return the lifetime of a star of mass m and metallicity z */ double star_lifetime(double z,double m) { /* z is the mass fraction of metals, ie, the metallicity */ /* m is the stellar mass in code unit */ /* Return t in code time unit */ int i; double a,b,c; double coeff[3]; double logm,twologm,logm2,time; /* convert m in msol */ m = m*All.CMUtoMsol; - + for (i=0;i<3;i++) coeff[i] = ( Cp->coeff_z[i][0]*z+Cp->coeff_z[i][1] )*z+Cp->coeff_z[i][2]; a = coeff[0]; b = coeff[1]; c = coeff[2]; logm = log10(m); twologm = 2.0 * logm; logm2 = logm*logm; time = pow(10.,(a*logm2+b*twologm+c)); time = time * All.HubbleParam; /* correct from hubble as not done in coeff_z */ return time; } /*! Return the mass of a star having a livetime t and a metallicity z */ double star_mass_from_age(double z,double t) { /* z is the mass fraction of metals, ie, the metallicity */ /* t is the star life time (in code unit) */ /* return the stellar mass (in code unit) that has a lifetime equal to t */ /* this is the inverse of star_lifetime */ int i; double a,b,c; double coeff[3]; double m; t = t/All.HubbleParam; /* correct from hubble as not done in coeff_z */ for (i=0;i<3;i++) coeff[i] = ( Cp->coeff_z[i][0]*z+Cp->coeff_z[i][1] )*z+Cp->coeff_z[i][2]; a = coeff[0]; b = coeff[1]; c = coeff[2]; m = -(b+sqrt(b*b-a*(c-log10(t))))/a; m = pow(10,m); /* here, m is in solar mass */ m = m*All.MsoltoCMU; /* Msol to mass unit */ return m; } /****************************************************************************************/ /* /* Supernova rate : number of supernova per mass unit /* /****************************************************************************************/ double DYIN_rate(double m1,double m2) { /* compute the number of stars between m1 and m2 masses in code unit */ double RDYIN; double md,mu; RDYIN = 0.0; /* convert m in msol */ m1 = m1*All.CMUtoMsol; m2 = m2*All.CMUtoMsol; /* find md, mu */ md = dmax(m1,Cp->DYIN_Mmin); mu = dmin(m2,Cp->DYIN_Mmax); if (mu<=md) /* no dying stars in that mass range */ return 0.0; /* to code units */ md = md * All.MsoltoCMU; mu = mu * All.MsoltoCMU; /* compute number in the mass range */ RDYIN = get_imf_N(md,mu); return RDYIN; } double SNII_rate(double m1,double m2) { /* compute the number of SNII between m1 and m2 masses in code unit */ double RSNII; double md,mu; RSNII = 0.0; /* convert m in msol */ m1 = m1*All.CMUtoMsol; m2 = m2*All.CMUtoMsol; /* (1) find md, mu */ md = dmax(m1,Cp->SNII_Mmin); mu = dmin(m2,Cp->SNII_Mmax); if (mu<=md) /* no SNII in that mass range */ return 0.0; /* !!!!! here we should use get_imf_N !!!! */ /* to ensure the full imf */ //RSNII = Cp->SNII_cte * (pow(mu,Cp->SNII_a)-pow(md,Cp->SNII_a)); /* number per solar mass */ ///* convert in number per solar mass to number per mass unit */ //RSNII = RSNII *All.CMUtoMsol; /* to code units */ md = md * All.MsoltoCMU; mu = mu * All.MsoltoCMU; /* compute number in the mass range */ RSNII = get_imf_N(md,mu); return RSNII; } double SNIa_rate(double m1,double m2) { /* compute the number of SNIa between m1 and m2 masses in code unit */ double RSNIa; double md,mu; RSNIa = 0.0; if (Cp->SNIa_Mpl<0) return RSNIa; /* convert m in msol */ m1 = m1*All.CMUtoMsol; m2 = m2*All.CMUtoMsol; /* RG contribution */ md = dmax(m1,Cp->SNIa_Mdl1); mu = dmin(m2,Cp->SNIa_Mdu1); if (mdSNIa_bb1 * Cp->SNIa_cte1 * (pow(mu,Cp->SNIa_a1)-pow(md,Cp->SNIa_a1)); /* MS contribution */ md = dmax(m1,Cp->SNIa_Mdl2); mu = dmin(m2,Cp->SNIa_Mdu2); if (mdSNIa_bb2 * Cp->SNIa_cte2 * (pow(mu,Cp->SNIa_a2)-pow(md,Cp->SNIa_a2)); /* WD contribution */ md = dmax(m1,Cp->SNIa_Mpl); /* select stars that have finished their life -> WD */ mu = Cp->SNIa_Mpu; /* no upper bond */ if (mu<=md) /* no SNIa in that mass range */ return 0.0; RSNIa = RSNIa * Cp->SNIa_cte * (pow(mu,Cp->SNIa_a)-pow(md,Cp->SNIa_a)); /* number per solar mass */ /* convert in number per solar mass to number per mass unit */ RSNIa = RSNIa *All.CMUtoMsol; return RSNIa; } double SNIa_single_rate(double m) { /* compute the number of SNIa per mass unit (in CMU) if the mass of the star is m m is given in code units. NOTE : as single stars are follow a normal imf distribution and here the RG and MS follow another IMF, we need to introduce the normalisation factor */ double RSNIa; RSNIa = 0.0; if (Cp->SNIa_Mpl<0) return RSNIa; /* convert m in msol */ m = m*All.CMUtoMsol; /* the star mass is greater than the WD interval */ if (m <= Cp->SNIa_Mpl ) { /* according to its mass the star is a RG */ if ( ( m > Cp->SNIa_Mdl1 ) && ( m <= Cp->SNIa_Mdu1 ) ) RSNIa = Cp->SNIa_bb1 * Cp->SNIa_NWDN * Cp->SNIa_b1 * pow(m,Cp->SNIa_a1 ) / get_imf(m*All.MsoltoCMU); /* according to its mass the star is a MS */ if ( ( m > Cp->SNIa_Mdl2 ) && ( m <= Cp->SNIa_Mdu2 ) ) RSNIa = Cp->SNIa_bb2 * Cp->SNIa_NWDN * Cp->SNIa_b2 * pow(m,Cp->SNIa_a2 ) / get_imf(m*All.MsoltoCMU); } return RSNIa; } void DYIN_mass_ejection(double m1,double m2) { /* Compute the mass fraction and yields of dying stars with masses between m1 and m2. Store the result in the global variable`` MassFracDYIN``:: MassFracDYIN[0] = total gas MassFracDYIN[1] = helium core (i.e. alpha(m)) MassFracDYIN[i] = frac mass elt i. */ double l1,l2; int i1,i2,i1p,i2p,j; double f1,f2; double v1,v2; /* convert m in msol */ m1 = m1*All.CMUtoMsol; m2 = m2*All.CMUtoMsol; /* this was not in Poirier... */ m1 = dmax(m1,Cp->DYIN_Mmin); m2 = dmin(m2,Cp->DYIN_Mmax); if ( m2<=m1 ) { for (j=0;jnelts+2;j++) MassFracDYIN[j] = 0; return; } j = 0; l1 = ( log10(m1) - Elt[j].MminDYIN) / Elt[j].StepDYIN ; l2 = ( log10(m2) - Elt[j].MminDYIN) / Elt[j].StepDYIN ; if (l1 < 0.0) l1 = 0.0; if (l2 < 0.0) l2 = 0.0; i1 = (int)l1; i2 = (int)l2; i1p = i1 + 1; i2p = i2 + 1; f1 = l1 - i1; f2 = l2 - i2; /* check (yr) */ if (i1<0) i1=0; if (i2<0) i2=0; /* --------- TOTAL GAS ---------- */ j = 0; v1 = f1 * ( Elt[j].ArrayDYIN[i1p] - Elt[j].ArrayDYIN[i1] ) + Elt[j].ArrayDYIN[i1]; v2 = f2 * ( Elt[j].ArrayDYIN[i2p] - Elt[j].ArrayDYIN[i2] ) + Elt[j].ArrayDYIN[i2]; MassFracDYIN[j] = v2-v1; /* --------- He core therm ---------- */ j = 1; v1 = f1 * ( Elt[j].ArrayDYIN[i1p] - Elt[j].ArrayDYIN[i1] ) + Elt[j].ArrayDYIN[i1]; v2 = f2 * ( Elt[j].ArrayDYIN[i2p] - Elt[j].ArrayDYIN[i2] ) + Elt[j].ArrayDYIN[i2]; MassFracDYIN[j] = v2-v1; /* ---------------------------- */ /* --------- Metals ---------- */ /* ---------------------------- */ j = 2; l1 = ( log10(m1) - Elt[j].MminDYIN) / Elt[j].StepDYIN ; l2 = ( log10(m2) - Elt[j].MminDYIN) / Elt[j].StepDYIN ; if (l1 < 0.0) l1 = 0.0; if (l2 < 0.0) l2 = 0.0; i1 = (int)l1; i2 = (int)l2; i1p = i1 + 1; i2p = i2 + 1; f1 = l1 - i1; f2 = l2 - i2; /* check (yr) */ if (i1<0) i1=0; if (i2<0) i2=0; for (j=2;jnelts+2;j++) { v1 = f1 * ( Elt[j].ArrayDYIN[i1p] - Elt[j].ArrayDYIN[i1] ) + Elt[j].ArrayDYIN[i1]; v2 = f2 * ( Elt[j].ArrayDYIN[i2p] - Elt[j].ArrayDYIN[i2] ) + Elt[j].ArrayDYIN[i2]; MassFracDYIN[j] = v2-v1; } } void DYIN_mass_ejection_Z(double m1,double m2, double Z) { /* Compute the mass fraction and yields of dying stars with masses between m1 and m2. Store the result in the global variable`` MassFracDYIN``:: MassFracDYIN[0] = total gas MassFracDYIN[1] = helium core (i.e. alpha(m)) MassFracDYIN[i] = frac mass elt i. */ double l1,l2; int i1,i2,i1p,i2p,k; int j1,j1p; double f1,f2,g1; double v1,v2,v11,v12; /* convert m in msol */ m1 = m1*All.CMUtoMsol; m2 = m2*All.CMUtoMsol; /* this was not in Poirier... */ m1 = dmax(m1,Cp->DYIN_Mmin); m2 = dmin(m2,Cp->DYIN_Mmax); if ( m2<=m1 ) { for (k=0;knelts+2;k++) MassFracDYIN[k] = 0; return; } /* find i1,i1p, index in m1 and m2 */ k = 0; l1 = ( log10(m1) - Elt[k].MminDYIN) / Elt[k].StepDYIN ; l2 = ( log10(m2) - Elt[k].MminDYIN) / Elt[k].StepDYIN ; if (l1 < 0.0) l1 = 0.0; if (l2 < 0.0) l2 = 0.0; i1 = (int)l1; i2 = (int)l2; i1p = i1 + 1; i2p = i2 + 1; f1 = l1 - i1; f2 = l2 - i2; /* check (yr) */ if (i1<0) i1=0; if (i2<0) i2=0; /* find j1,j1p, index in Z */ l1 = ( log10(Z) - Elt[k].MminZDYIN) / Elt[k].StepZDYIN ; if (l1 < 0.0) l1 = 0.0; j1 = (int)l1; j1p = j1 + 1; g1 = l1 - j1; /* check (yr) */ if (j1<0) j1=0; /* --------- TOTAL GAS ---------- */ k = 0; /* mass 1 */ v11 = f1 * ( Elt[k].ArrayZDYIN[i1p][j1] - Elt[k].ArrayZDYIN[i1][j1] ) + Elt[k].ArrayZDYIN[i1][j1]; v12 = f1 * ( Elt[k].ArrayZDYIN[i1p][j1p] - Elt[k].ArrayZDYIN[i1][j1p] ) + Elt[k].ArrayZDYIN[i1][j1p]; v1 = g1 * ( v12 - v11 ) + v11; /* mass 2 */ v11 = f2 * ( Elt[k].ArrayZDYIN[i2p][j1] - Elt[k].ArrayZDYIN[i2][j1] ) + Elt[k].ArrayZDYIN[i2][j1]; v12 = f2 * ( Elt[k].ArrayZDYIN[i2p][j1p] - Elt[k].ArrayZDYIN[i2][j1p] ) + Elt[k].ArrayZDYIN[i2][j1p]; v2 = g1 * ( v12 - v11 ) + v11; MassFracDYIN[k] = v2-v1; /* --------- He core therm ---------- */ k = 1; /* mass 1 */ v11 = f1 * ( Elt[k].ArrayZDYIN[i1p][j1] - Elt[k].ArrayZDYIN[i1][j1] ) + Elt[k].ArrayZDYIN[i1][j1]; v12 = f1 * ( Elt[k].ArrayZDYIN[i1p][j1p] - Elt[k].ArrayZDYIN[i1][j1p] ) + Elt[k].ArrayZDYIN[i1][j1p]; v1 = g1 * ( v12 - v11 ) + v11; /* mass 2 */ v11 = f2 * ( Elt[k].ArrayZDYIN[i2p][j1] - Elt[k].ArrayZDYIN[i2][j1] ) + Elt[k].ArrayZDYIN[i2][j1]; v12 = f2 * ( Elt[k].ArrayZDYIN[i2p][j1p] - Elt[k].ArrayZDYIN[i2][j1p] ) + Elt[k].ArrayZDYIN[i2][j1p]; v2 = g1 * ( v12 - v11 ) + v11; MassFracDYIN[k] = v2-v1; /* ---------------------------- */ /* --------- Metals ---------- */ /* ---------------------------- */ k = 2; l1 = ( log10(m1) - Elt[k].MminDYIN) / Elt[k].StepDYIN ; l2 = ( log10(m2) - Elt[k].MminDYIN) / Elt[k].StepDYIN ; if (l1 < 0.0) l1 = 0.0; if (l2 < 0.0) l2 = 0.0; i1 = (int)l1; i2 = (int)l2; i1p = i1 + 1; i2p = i2 + 1; f1 = l1 - i1; f2 = l2 - i2; /* check (yr) */ if (i1<0) i1=0; if (i2<0) i2=0; /* find j1,j1p, index in Z */ l1 = ( log10(Z) - Elt[k].MminZDYIN) / Elt[k].StepZDYIN ; if (l1 < 0.0) l1 = 0.0; j1 = (int)l1; j1p = j1 + 1; g1 = l1 - j1; /* check (yr) */ if (j1<0) j1=0; for (k=2;knelts+2;k++) { /* mass 1 */ v11 = f1 * ( Elt[k].ArrayZDYIN[i1p][j1] - Elt[k].ArrayZDYIN[i1][j1] ) + Elt[k].ArrayZDYIN[i1][j1]; v12 = f1 * ( Elt[k].ArrayZDYIN[i1p][j1p] - Elt[k].ArrayZDYIN[i1][j1p] ) + Elt[k].ArrayZDYIN[i1][j1p]; v1 = g1 * ( v12 - v11 ) + v11; /* mass 2 */ v11 = f2 * ( Elt[k].ArrayZDYIN[i2p][j1] - Elt[k].ArrayZDYIN[i2][j1] ) + Elt[k].ArrayZDYIN[i2][j1]; v12 = f2 * ( Elt[k].ArrayZDYIN[i2p][j1p] - Elt[k].ArrayZDYIN[i2][j1p] ) + Elt[k].ArrayZDYIN[i2][j1p]; v2 = g1 * ( v12 - v11 ) + v11; MassFracDYIN[k] = v2-v1; } } void DYIN_single_mass_ejection(double m1) { /* Compute the mass fraction and yields of a dying stars of masse m1. Store the result in the global variable ``SingleMassFracDYIN``:: SingleMassFracDYIN[0] = total gas SingleMassFracDYIN[1] = helium core (i.e. alpha(m)) SingleMassFracDYIN[i] = frac mass elt i. */ double l1; int i1,i1p,j; double f1; double v1; /* convert m in msol */ m1 = m1*All.CMUtoMsol; /* this was not in Poirier... */ if ( (m1<=Cp->DYIN_Mmin) || (m1>=Cp->DYIN_Mmax) ) { for (j=0;jnelts+2;j++) SingleMassFracDYIN[j] = 0; return; } j = 0; l1 = ( log10(m1) - Elt[j].MminDYIN) / Elt[j].StepDYIN ; if (l1 < 0.0) l1 = 0.0; i1 = (int)l1; i1p = i1 + 1; f1 = l1 - i1; /* check (yr) */ if (i1<0) i1=0; /* --------- TOTAL GAS ---------- */ j = 0; v1 = f1 * ( Elt[j].MetalDYIN[i1p] - Elt[j].MetalDYIN[i1] ) + Elt[j].MetalDYIN[i1]; SingleMassFracDYIN[j] = v1; /* --------- He core therm ---------- */ j = 1; v1 = f1 * ( Elt[j].MetalDYIN[i1p] - Elt[j].MetalDYIN[i1] ) + Elt[j].MetalDYIN[i1]; SingleMassFracDYIN[j] = v1; /* ---------------------------- */ /* --------- Metals ---------- */ /* ---------------------------- */ j = 2; l1 = ( log10(m1) - Elt[j].MminDYIN) / Elt[j].StepDYIN ; if (l1 < 0.0) l1 = 0.0; i1 = (int)l1; i1p = i1 + 1; f1 = l1 - i1; /* check (yr) */ if (i1<0) i1=0; for (j=2;jnelts+2;j++) { v1 = f1 * ( Elt[j].MetalDYIN[i1p] - Elt[j].MetalDYIN[i1] ) + Elt[j].MetalDYIN[i1]; SingleMassFracDYIN[j] = v1; } } void DYIN_single_mass_ejection_Z(double m1,double Z) { /* Compute the mass fraction and yields of a dying stars of masse m1. Store the result in the global variable ``SingleMassFracDYIN``:: SingleMassFracDYIN[0] = total gas SingleMassFracDYIN[1] = helium core (i.e. alpha(m)) SingleMassFracDYIN[i] = frac mass elt i. */ double l1; int i1,i1p,k; int j1,j1p; double f1,g1; double v1,v11,v12; /* convert m in msol */ m1 = m1*All.CMUtoMsol; /* this was not in Poirier... */ if ( (m1<=Cp->DYIN_Mmin) || (m1>=Cp->DYIN_Mmax) ) { for (k=0;knelts+2;k++) SingleMassFracDYIN[k] = 0; return; } k = 0; /* find i1,i1p, index in m */ l1 = ( log10(m1) - Elt[k].MminDYIN) / Elt[k].StepDYIN ; if (l1 < 0.0) l1 = 0.0; i1 = (int)l1; i1p = i1 + 1; f1 = l1 - i1; /* check (yr) */ if (i1<0) i1=0; /* find j1,j1p, index in Z */ l1 = ( log10(Z) - Elt[k].MminZDYIN) / Elt[k].StepZDYIN ; if (l1 < 0.0) l1 = 0.0; j1 = (int)l1; j1p = j1 + 1; g1 = l1 - j1; /* check (yr) */ if (j1<0) j1=0; /* --------- TOTAL GAS ---------- */ k = 0; v11 = f1 * ( Elt[k].MetalZDYIN[i1p][j1] - Elt[k].MetalZDYIN[i1][j1] ) + Elt[k].MetalZDYIN[i1][j1]; v12 = f1 * ( Elt[k].MetalZDYIN[i1p][j1p] - Elt[k].MetalZDYIN[i1][j1p] ) + Elt[k].MetalZDYIN[i1p][j1]; v1 = g1 * ( v12 - v11 ) + v11; SingleMassFracDYIN[k] = v1; /* --------- He core therm ---------- */ k = 1; v11 = f1 * ( Elt[k].MetalZDYIN[i1p][j1] - Elt[k].MetalZDYIN[i1][j1] ) + Elt[k].MetalZDYIN[i1][j1]; v12 = f1 * ( Elt[k].MetalZDYIN[i1p][j1p] - Elt[k].MetalZDYIN[i1][j1p] ) + Elt[k].MetalZDYIN[i1p][j1]; v1 = g1 * ( v12 - v11 ) + v11; SingleMassFracDYIN[k] = v1; /* ---------------------------- */ /* --------- Metals ---------- */ /* ---------------------------- */ k = 2; /* find i1,i1p, index in m */ l1 = ( log10(m1) - Elt[k].MminDYIN) / Elt[k].StepDYIN ; if (l1 < 0.0) l1 = 0.0; i1 = (int)l1; i1p = i1 + 1; f1 = l1 - i1; /* check (yr) */ if (i1<0) i1=0; /* find j1,j1p, index in Z */ l1 = ( log10(Z) - Elt[k].MminZDYIN) / Elt[k].StepZDYIN ; if (l1 < 0.0) l1 = 0.0; j1 = (int)l1; j1p = j1 + 1; g1 = l1 - j1; /* check (yr) */ if (j1<0) j1=0; for (k=2;knelts+2;k++) { v11 = f1 * ( Elt[k].MetalZDYIN[i1p][j1] - Elt[k].MetalZDYIN[i1][j1] ) + Elt[k].MetalZDYIN[i1][j1]; v12 = f1 * ( Elt[k].MetalZDYIN[i1p][j1p] - Elt[k].MetalZDYIN[i1][j1p] ) + Elt[k].MetalZDYIN[i1p][j1]; v1 = g1 * ( v12 - v11 ) + v11; SingleMassFracDYIN[k] = v1; } } void SNII_mass_ejection(double m1,double m2) { /* .. warning:: here, we we do not limit the computation to SNII !!! Compute the mass fraction and yields of SNII stars with masses between m1 and m2. Store the result in the global variable ``MassFracSNII``:: MassFracSNII[0] = total gas MassFracSNII[1] = 1-helium core (i.e. non processed elts) MassFracSNII[i] = frac mass elt i. */ double l1,l2; int i1,i2,i1p,i2p,j; double f1,f2; double v1,v2; /* convert m in msol */ m1 = m1*All.CMUtoMsol; m2 = m2*All.CMUtoMsol; /* this was not in Poirier... */ m1 = dmax(m1,Cp->SNII_Mmin); m2 = dmin(m2,Cp->SNII_Mmax); if ( m2<=m1 ) { for (j=0;jnelts+2;j++) MassFracSNII[j] = 0; return; } j = 0; l1 = ( log10(m1) - Elt[j].MminSNII) / Elt[j].StepSNII ; l2 = ( log10(m2) - Elt[j].MminSNII) / Elt[j].StepSNII ; if (l1 < 0.0) l1 = 0.0; if (l2 < 0.0) l2 = 0.0; i1 = (int)l1; i2 = (int)l2; i1p = i1 + 1; i2p = i2 + 1; f1 = l1 - i1; f2 = l2 - i2; /* check (yr) */ if (i1<0) i1=0; if (i2<0) i2=0; /* --------- TOTAL GAS ---------- */ j = 0; v1 = f1 * ( Elt[j].ArraySNII[i1p] - Elt[j].ArraySNII[i1] ) + Elt[j].ArraySNII[i1]; v2 = f2 * ( Elt[j].ArraySNII[i2p] - Elt[j].ArraySNII[i2] ) + Elt[j].ArraySNII[i2]; MassFracSNII[j] = v2-v1; /* --------- He core therm ---------- */ j = 1; v1 = f1 * ( Elt[j].ArraySNII[i1p] - Elt[j].ArraySNII[i1] ) + Elt[j].ArraySNII[i1]; v2 = f2 * ( Elt[j].ArraySNII[i2p] - Elt[j].ArraySNII[i2] ) + Elt[j].ArraySNII[i2]; MassFracSNII[j] = v2-v1; /* ---------------------------- */ /* --------- Metals ---------- */ /* ---------------------------- */ j = 2; l1 = ( log10(m1) - Elt[j].MminSNII) / Elt[j].StepSNII ; l2 = ( log10(m2) - Elt[j].MminSNII) / Elt[j].StepSNII ; if (l1 < 0.0) l1 = 0.0; if (l2 < 0.0) l2 = 0.0; i1 = (int)l1; i2 = (int)l2; i1p = i1 + 1; i2p = i2 + 1; f1 = l1 - i1; f2 = l2 - i2; /* check (yr) */ if (i1<0) i1=0; if (i2<0) i2=0; for (j=2;jnelts+2;j++) { v1 = f1 * ( Elt[j].ArraySNII[i1p] - Elt[j].ArraySNII[i1] ) + Elt[j].ArraySNII[i1]; v2 = f2 * ( Elt[j].ArraySNII[i2p] - Elt[j].ArraySNII[i2] ) + Elt[j].ArraySNII[i2]; MassFracSNII[j] = v2-v1; } } // void SNII_mass_ejection_Z(double m1,double m2, double Z) // { // // /* // Compute the mass fraction and yields of SNII with masses between m1 and m2. // Store the result in the global variable`` MassFracSNII``:: // // MassFracSNII[0] = total gas // MassFracSNII[1] = helium core (i.e. alpha(m)) // MassFracSNII[i] = frac mass elt i. // // */ // // double l1,l2; // int i1,i2,i1p,i2p,k; // int j1,j1p; // double f1,f2,g1; // double v1,v2,v11,v12; // // /* convert m in msol */ // m1 = m1*All.CMUtoMsol; // m2 = m2*All.CMUtoMsol; // // /* this was not in Poirier... */ // m1 = dmax(m1,Cp->SNII_Mmin); // m2 = dmin(m2,Cp->SNII_Mmax); // // // if ( m2<=m1 ) // { // for (k=0;knelts+2;k++) // MassFracSNII[k] = 0; // return; // } // // // // // /* find i1,i1p, index in m1 and m2 */ // // k = 0; // l1 = ( log10(m1) - Elt[k].MminSNII) / Elt[k].StepSNII ; // l2 = ( log10(m2) - Elt[k].MminSNII) / Elt[k].StepSNII ; // // if (l1 < 0.0) l1 = 0.0; // if (l2 < 0.0) l2 = 0.0; // // i1 = (int)l1; // i2 = (int)l2; // // i1p = i1 + 1; // i2p = i2 + 1; // // f1 = l1 - i1; // f2 = l2 - i2; // // /* check (yr) */ // if (i1<0) i1=0; // if (i2<0) i2=0; // // // /* find j1,j1p, index in Z */ // // l1 = ( log10(Z) - Elt[k].MminZSNII) / Elt[k].StepZSNII ; // // if (l1 < 0.0) l1 = 0.0; // // j1 = (int)l1; // j1p = j1 + 1; // g1 = l1 - j1; // // /* check (yr) */ // if (j1<0) j1=0; // // // // /* --------- TOTAL GAS ---------- */ // k = 0; // // /* mass 1 */ // v11 = f1 * ( Elt[k].ArrayZSNII[i1p][j1] - Elt[k].ArrayZSNII[i1][j1] ) + Elt[k].ArrayZSNII[i1][j1]; // v12 = f1 * ( Elt[k].ArrayZSNII[i1p][j1p] - Elt[k].ArrayZSNII[i1][j1p] ) + Elt[k].ArrayZSNII[i1p][j1]; // v1 = g1 * ( v12 - v11 ) + v11; // // /* mass 2 */ // v11 = f2 * ( Elt[k].ArrayZSNII[i2p][j1] - Elt[k].ArrayZSNII[i2][j1] ) + Elt[k].ArrayZSNII[i2][j1]; // v12 = f2 * ( Elt[k].ArrayZSNII[i2p][j1p] - Elt[k].ArrayZSNII[i2][j1p] ) + Elt[k].ArrayZSNII[i2p][j1]; // v2 = g1 * ( v12 - v11 ) + v11; // // MassFracSNII[k] = v2-v1; // // /* --------- He core therm ---------- */ // k = 1; // // /* mass 1 */ // v11 = f1 * ( Elt[k].ArrayZSNII[i1p][j1] - Elt[k].ArrayZSNII[i1][j1] ) + Elt[k].ArrayZSNII[i1][j1]; // v12 = f1 * ( Elt[k].ArrayZSNII[i1p][j1p] - Elt[k].ArrayZSNII[i1][j1p] ) + Elt[k].ArrayZSNII[i1p][j1]; // v1 = g1 * ( v12 - v11 ) + v11; // // /* mass 2 */ // v11 = f2 * ( Elt[k].ArrayZSNII[i2p][j1] - Elt[k].ArrayZSNII[i2][j1] ) + Elt[k].ArrayZSNII[i2][j1]; // v12 = f2 * ( Elt[k].ArrayZSNII[i2p][j1p] - Elt[k].ArrayZSNII[i2][j1p] ) + Elt[k].ArrayZSNII[i2p][j1]; // v2 = g1 * ( v12 - v11 ) + v11; // // MassFracSNII[k] = v2-v1; // // /* ---------------------------- */ // /* --------- Metals ---------- */ // /* ---------------------------- */ // // k = 2; // // l1 = ( log10(m1) - Elt[k].MminSNII) / Elt[k].StepSNII ; // l2 = ( log10(m2) - Elt[k].MminSNII) / Elt[k].StepSNII ; // // if (l1 < 0.0) l1 = 0.0; // if (l2 < 0.0) l2 = 0.0; // // i1 = (int)l1; // i2 = (int)l2; // // i1p = i1 + 1; // i2p = i2 + 1; // // f1 = l1 - i1; // f2 = l2 - i2; // // /* check (yr) */ // if (i1<0) i1=0; // if (i2<0) i2=0; // // /* find j1,j1p, index in Z */ // // l1 = ( log10(Z) - Elt[k].MminZSNII) / Elt[k].StepZSNII ; // // if (l1 < 0.0) l1 = 0.0; // // j1 = (int)l1; // j1p = j1 + 1; // g1 = l1 - j1; // // /* check (yr) */ // if (j1<0) j1=0; // // for (k=2;knelts+2;k++) // { // // /* mass 1 */ // v11 = f1 * ( Elt[k].ArrayZSNII[i1p][j1] - Elt[k].ArrayZSNII[i1][j1] ) + Elt[k].ArrayZSNII[i1][j1]; // v12 = f1 * ( Elt[k].ArrayZSNII[i1p][j1p] - Elt[k].ArrayZSNII[i1][j1p] ) + Elt[k].ArrayZSNII[i1p][j1]; // v1 = g1 * ( v12 - v11 ) + v11; // // /* mass 2 */ // v11 = f2 * ( Elt[k].ArrayZSNII[i2p][j1] - Elt[k].ArrayZSNII[i2][j1] ) + Elt[k].ArrayZSNII[i2][j1]; // v12 = f2 * ( Elt[k].ArrayZSNII[i2p][j1p] - Elt[k].ArrayZSNII[i2][j1p] ) + Elt[k].ArrayZSNII[i2p][j1]; // v2 = g1 * ( v12 - v11 ) + v11; // // MassFracSNII[k] = v2-v1; // // } // // } void SNII_single_mass_ejection(double m1) { /* .. warning:: here, we we do not limit the computation to SNII !!! Compute the mass fraction and yields of a SNII stars of masse m1. Store the result in the global variable ``SingleMassFracSNII``:: SingleMassFracSNII[0] = total gas SingleMassFracSNII[1] = 1-helium core (i.e. non processed elts) SingleMassFracSNII[i] = frac mass elt i. */ double l1; int i1,i1p,j; double f1; double v1; /* convert m in msol */ m1 = m1*All.CMUtoMsol; /* this was not in Poirier... */ if ( (m1<=Cp->SNII_Mmin) || (m1>=Cp->SNII_Mmax) ) { for (j=0;jnelts+2;j++) SingleMassFracSNII[j] = 0; return; } j = 0; l1 = ( log10(m1) - Elt[j].MminSNII) / Elt[j].StepSNII ; if (l1 < 0.0) l1 = 0.0; i1 = (int)l1; i1p = i1 + 1; f1 = l1 - i1; /* check (yr) */ if (i1<0) i1=0; /* --------- TOTAL GAS ---------- */ j = 0; v1 = f1 * ( Elt[j].MetalSNII[i1p] - Elt[j].MetalSNII[i1] ) + Elt[j].MetalSNII[i1]; SingleMassFracSNII[j] = v1; /* --------- He core therm ---------- */ j = 1; v1 = f1 * ( Elt[j].MetalSNII[i1p] - Elt[j].MetalSNII[i1] ) + Elt[j].MetalSNII[i1]; SingleMassFracSNII[j] = v1; /* ---------------------------- */ /* --------- Metals ---------- */ /* ---------------------------- */ j = 2; l1 = ( log10(m1) - Elt[j].MminSNII) / Elt[j].StepSNII ; if (l1 < 0.0) l1 = 0.0; i1 = (int)l1; i1p = i1 + 1; f1 = l1 - i1; /* check (yr) */ if (i1<0) i1=0; for (j=2;jnelts+2;j++) { v1 = f1 * ( Elt[j].MetalSNII[i1p] - Elt[j].MetalSNII[i1] ) + Elt[j].MetalSNII[i1]; SingleMassFracSNII[j] = v1; } } // void SNII_single_mass_ejection_Z(double m1,double Z) // { // // // /* // Compute the mass fraction and yields of a SNII stars of masse m1. // Store the result in the global variable ``SingleMassFracSNII``:: // // SingleMassFracSNII[0] = total gas // SingleMassFracSNII[1] = helium core (i.e. alpha(m)) // SingleMassFracSNII[i] = frac mass elt i. // */ // // double l1; // int i1,i1p,k; // int j1,j1p; // double f1,g1; // double v1,v11,v12; // // /* convert m in msol */ // m1 = m1*All.CMUtoMsol; // // /* this was not in Poirier... */ // if ( (m1<=Cp->SNII_Mmin) || (m1>=Cp->SNII_Mmax) ) // { // for (k=0;knelts+2;k++) // SingleMassFracSNII[k] = 0; // return; // } // // // k = 0; // // // /* find i1,i1p, index in m */ // // l1 = ( log10(m1) - Elt[k].MminSNII) / Elt[k].StepSNII ; // // if (l1 < 0.0) l1 = 0.0; // // i1 = (int)l1; // i1p = i1 + 1; // f1 = l1 - i1; // // /* check (yr) */ // if (i1<0) i1=0; // // // /* find j1,j1p, index in Z */ // // l1 = ( log10(Z) - Elt[k].MminZSNII) / Elt[k].StepZSNII ; // // if (l1 < 0.0) l1 = 0.0; // // j1 = (int)l1; // j1p = j1 + 1; // g1 = l1 - j1; // // /* check (yr) */ // if (j1<0) j1=0; // // // // /* --------- TOTAL GAS ---------- */ // k = 0; // v11 = f1 * ( Elt[k].MetalZSNII[i1p][j1] - Elt[k].MetalZSNII[i1][j1] ) + Elt[k].MetalZSNII[i1][j1]; // v12 = f1 * ( Elt[k].MetalZSNII[i1p][j1p] - Elt[k].MetalZSNII[i1][j1p] ) + Elt[k].MetalZSNII[i1p][j1]; // v1 = g1 * ( v12 - v11 ) + v11; // // SingleMassFracSNII[k] = v1; // // /* --------- He core therm ---------- */ // k = 1; // v11 = f1 * ( Elt[k].MetalZSNII[i1p][j1] - Elt[k].MetalZSNII[i1][j1] ) + Elt[k].MetalZSNII[i1][j1]; // v12 = f1 * ( Elt[k].MetalZSNII[i1p][j1p] - Elt[k].MetalZSNII[i1][j1p] ) + Elt[k].MetalZSNII[i1p][j1]; // v1 = g1 * ( v12 - v11 ) + v11; // // SingleMassFracSNII[k] = v1; // // /* ---------------------------- */ // /* --------- Metals ---------- */ // /* ---------------------------- */ // // k = 2; // // /* find i1,i1p, index in m */ // // l1 = ( log10(m1) - Elt[k].MminSNII) / Elt[k].StepSNII ; // // if (l1 < 0.0) l1 = 0.0; // // i1 = (int)l1; // i1p = i1 + 1; // f1 = l1 - i1; // // /* check (yr) */ // if (i1<0) i1=0; // // // /* find j1,j1p, index in Z */ // // l1 = ( log10(Z) - Elt[k].MminZSNII) / Elt[k].StepZSNII ; // // if (l1 < 0.0) l1 = 0.0; // // j1 = (int)l1; // j1p = j1 + 1; // g1 = l1 - j1; // // /* check (yr) */ // if (j1<0) j1=0; // // for (k=2;knelts+2;k++) // { // // v11 = f1 * ( Elt[k].MetalZSNII[i1p][j1] - Elt[k].MetalZSNII[i1][j1] ) + Elt[k].MetalZSNII[i1][j1]; // v12 = f1 * ( Elt[k].MetalZSNII[i1p][j1p] - Elt[k].MetalZSNII[i1][j1p] ) + Elt[k].MetalZSNII[i1p][j1]; // v1 = g1 * ( v12 - v11 ) + v11; // // SingleMassFracSNII[k] = v1; // } // // } void SNIa_mass_ejection(double m1,double m2) { /* Compute the total mass and element mass per mass unit of SNIa stars with masses between m1 and m2. Store the result in the global variable ``MassFracSNIa``:: MassFracSNIa[0] = total gas MassFracSNIa[1] = unused MassFracSNIa[i] = frac mass elt i. */ int j; double NSNIa; /* number of SNIa per mass unit between time and time+dt */ NSNIa = SNIa_rate(m1,m2); /* ejected mass in gas per mass unit */ MassFracSNIa[0] = Cp->Mco * All.MsoltoCMU * NSNIa; /* ejected elements in gas per mass unit */ for (j=2;jnelts+2;j++) MassFracSNIa[j] = NSNIa* Elt[j].MSNIa * All.MsoltoCMU; /* unused */ MassFracSNIa[1]=-1; } void SNIa_single_mass_ejection(double m1) { /* Compute the total mass mass of element of a SNIa stars of masse m1. Store the result in the global variable ``SingleMassFracSNIa``:: SingleMassFracSNIa[0] = total gas SingleMassFracSNIa[1] = unused SingleMassFracSNIa[i] = frac mass elt i. */ int j; if (Cp->SNIa_Mpl<0) // no SNIa { /* total ejected gas mass */ SingleMassFracSNIa[0] = 0; /* ejected mass per element */ for (j=2;jnelts+2;j++) SingleMassFracSNIa[j] = 0; /* unused */ SingleMassFracSNIa[1] = 0; } else { /* total ejected gas mass */ SingleMassFracSNIa[0] = Cp->Mco * All.MsoltoCMU; /* ejected mass per element */ for (j=2;jnelts+2;j++) SingleMassFracSNIa[j] = Elt[j].MSNIa * All.MsoltoCMU; /* unused */ SingleMassFracSNIa[1] = -1; } } void Total_mass_ejection(double m1,double m2,double M0,double *z) { /* Sum the contribution in mass and yields of all stars in the mass range m1,m2. Store the result in the global variable EjectedMass:: EjectedMass[0] = total gas EjectedMass[1] = UNUSED EjectedMass[i+2] = frac mass elt i. FOR THE MOMENT:: - contrib of SNII (= all stars) - contrib of SNIa EjectedMass[0] = ejected Mass from SNII + Mco * number of SNIa EjectedMass[i] = (SNII elts created ) + (SNII elts existing) + (SNIa elts) */ int j; /* compute SNII mass ejection -> MassFracSNII */ SNII_mass_ejection(m1,m2); /* compute SNIa mass ejection -> MassFracSNIa */ /* not really a mass fraction */ SNIa_mass_ejection(m1,m2); /* compute DYIN mass ejection -> MassFracDYIN */ /* not really a mass fraction */ #ifdef CHIMIE_AGB DYIN_mass_ejection(m1,m2); #endif /* total ejected gas mass */ #ifdef CHIMIE_AGB EjectedMass[0] = M0 * ( MassFracDYIN[0] + MassFracSNII[0] + MassFracSNIa[0] ); #else EjectedMass[0] = M0 * ( MassFracSNII[0] + MassFracSNIa[0] ); #endif /* ejected mass per element */ for (j=2;jnelts+2;j++) #ifdef CHIMIE_AGB EjectedMass[j] = M0*( MassFracDYIN[j] +z[j-2]*MassFracDYIN[1] + MassFracSNII[j] +z[j-2]*MassFracSNII[1] + MassFracSNIa[j] ); #else EjectedMass[j] = M0*( MassFracSNII[j] +z[j-2]*MassFracSNII[1] + MassFracSNIa[j] ); #endif /* not used */ EjectedMass[1] = -1; } void Total_mass_ejection_Z(double m1,double m2,double M0,double *z) { /* Sum the contribution in mass and yields of all stars in the mass range m1,m2. Store the result in the global variable EjectedMass:: EjectedMass[0] = total gas EjectedMass[1] = UNUSED EjectedMass[i+2] = frac mass elt i. FOR THE MOMENT:: - contrib of SNII (= all stars) - contrib of SNIa EjectedMass[0] = ejected Mass from SNII + Mco * number of SNIa EjectedMass[i] = (SNII elts created ) + (SNII elts existing) + (SNIa elts) */ int j; float Z; /* metalicity */ Z = z[METALS]; /* compute SNII mass ejection -> MassFracSNII */ SNII_mass_ejection(m1,m2); /* compute SNIa mass ejection -> MassFracSNIa */ /* not really a mass fraction */ SNIa_mass_ejection(m1,m2); /* compute DYIN mass ejection -> MassFracDYIN */ /* not really a mass fraction */ #ifdef CHIMIE_AGB DYIN_mass_ejection_Z(m1,m2,Z); #endif /* total ejected gas mass */ #ifdef CHIMIE_AGB EjectedMass[0] = M0 * ( MassFracDYIN[0] + MassFracSNII[0] + MassFracSNIa[0] ); #else EjectedMass[0] = M0 * ( MassFracSNII[0] + MassFracSNIa[0] ); #endif /* ejected mass per element */ for (j=2;jnelts+2;j++) #ifdef CHIMIE_AGB EjectedMass[j] = M0*( MassFracDYIN[j] +z[j-2]*MassFracDYIN[1] + MassFracSNII[j] +z[j-2]*MassFracSNII[1] + MassFracSNIa[j] ); #else EjectedMass[j] = M0*( MassFracSNII[j] +z[j-2]*MassFracSNII[1] + MassFracSNIa[j] ); #endif /* not used */ EjectedMass[1] = -1; } void DYIN_Total_single_mass_ejection(double m1,double *z) { /* Mass and element ejected by a single dying stars of mass m1. This takes into account processed and non processed gas The results are stored in:: SingleEjectedMass[0] = gas mass SingleEjectedMass[1] = unsued SingleEjectedMass[i+2] = frac mass elt i */ int j; float M0; M0 = m1; /* compute dying stars mass ejection -> SingleMassFracDYIN */ DYIN_single_mass_ejection(m1); /* total ejected gas mass */ SingleEjectedMass[0] = M0 * SingleMassFracDYIN[0]; /* ejected mass per element */ for (j=2;jnelts+2;j++) SingleEjectedMass[j] = M0*(SingleMassFracDYIN[j] +z[j-2]*SingleMassFracDYIN[1]); /* not used */ SingleEjectedMass[1] = -1; } void SNII_Total_single_mass_ejection(double m1,double *z) { /* Mass and element ejected by a single SNII of mass m1. This takes into account processed and non processed gas The results are stored in:: SingleEjectedMass[0] = gas mass SingleEjectedMass[1] = unsued SingleEjectedMass[i+2] = frac mass elt i */ int j; float M0; M0 = m1; /* compute SNII mass ejection -> SingleMassFracSNII */ SNII_single_mass_ejection(m1); /* total ejected gas mass */ SingleEjectedMass[0] = M0 * SingleMassFracSNII[0]; /* ejected mass per element */ for (j=2;jnelts+2;j++) SingleEjectedMass[j] = M0*(SingleMassFracSNII[j] +z[j-2]*SingleMassFracSNII[1]); /* not used */ SingleEjectedMass[1] = -1; } void SNIa_Total_single_mass_ejection(double m1, double *z) { int j; /* Mass and element ejected by a single SNIa of mass m1. The results are stored in:: SingleEjectedMass[0] = gas mass SingleEjectedMass[1] = unsued SingleEjectedMass[i+2] = frac mass elt i */ /* compute SNIa mass ejection -> SingleMassFracSNIa */ SNIa_single_mass_ejection(m1); /* total ejected gas mass */ SingleEjectedMass[0] = SingleMassFracSNIa[0]; /* ejected mass per element */ for (j=2;jnelts+2;j++) SingleEjectedMass[j] = SingleMassFracSNIa[j]; } void Total_single_mass_ejection(double m1,double *z,double NSNII,double NSNIa,double NDYIN) { /* Sum the contribution in mass and yields of one star for mass m1. Store the result in the global variable EjectedMass:: SingleEjectedMass[0] = total gas SingleEjectedMass[1] = UNUSED SingleEjectedMass[i+2] = frac mass elt i. FOR THE MOMENT:: - contrib of SNII (= all stars) - contrib of SNIa SingleEjectedMass[0] = ejected Mass from SNII + Mco * number of SNIa SingleEjectedMass[i] = (SNII elts created ) + (SNII elts existing) + (SNIa elts) */ int j; float M0; M0 = m1; /* compute SNII mass ejection -> SingleMassFracSNII */ SNII_single_mass_ejection(m1); /* compute SNII mass ejection -> SingleMassFracSNIa */ SNIa_single_mass_ejection(m1); /* compute DYIN mass ejection -> SingleMassFracDYIN */ #ifdef CHIMIE_AGB DYIN_single_mass_ejection(m1); #endif /* total ejected gas mass */ #ifdef CHIMIE_AGB SingleEjectedMass[0] = M0 * ( SingleMassFracDYIN[0]*NDYIN + SingleMassFracSNII[0]*NSNII ) + SingleMassFracSNIa[0]*NSNIa; #else SingleEjectedMass[0] = M0 * ( SingleMassFracSNII[0]*NSNII ) + SingleMassFracSNIa[0]*NSNIa; #endif /* ejected mass per element */ for (j=2;jnelts+2;j++) #ifdef CHIMIE_AGB SingleEjectedMass[j] = M0*( SingleMassFracDYIN[j]*NDYIN +z[j-2]*SingleMassFracDYIN[1]*NDYIN + SingleMassFracSNII[j]*NSNII +z[j-2]*SingleMassFracSNII[1]*NSNII ) + SingleMassFracSNIa[j]*NSNIa; #else SingleEjectedMass[j] = M0*( SingleMassFracSNII[j]*NSNII +z[j-2]*SingleMassFracSNII[1]*NSNII ) + SingleMassFracSNIa[j]*NSNIa; #endif /* not used */ SingleEjectedMass[1] = -1; } void Total_single_mass_ejection_Z(double m1,double *z,double NSNII,double NSNIa,double NDYIN) { /* Sum the contribution in mass and yields of one star for mass m1. Store the result in the global variable EjectedMass:: SingleEjectedMass[0] = total gas SingleEjectedMass[1] = UNUSED SingleEjectedMass[i+2] = frac mass elt i. FOR THE MOMENT:: - contrib of SNII (= all stars) - contrib of SNIa SingleEjectedMass[0] = ejected Mass from SNII + Mco * number of SNIa SingleEjectedMass[i] = (SNII elts created ) + (SNII elts existing) + (SNIa elts) */ int j; float M0; float Z; /* metalicity */ M0 = m1; Z = z[METALS]; //z[METALS]=0; /* compute SNII mass ejection -> SingleMassFracSNII */ SNII_single_mass_ejection(m1); /* compute SNII mass ejection -> SingleMassFracSNIa */ SNIa_single_mass_ejection(m1); /* compute DYIN mass ejection -> SingleMassFracDYIN */ #ifdef CHIMIE_AGB DYIN_single_mass_ejection_Z(m1,Z); #endif /* total ejected gas mass */ #ifdef CHIMIE_AGB SingleEjectedMass[0] = M0 * ( SingleMassFracDYIN[0]*NDYIN + SingleMassFracSNII[0]*NSNII ) + SingleMassFracSNIa[0]*NSNIa; #else SingleEjectedMass[0] = M0 * ( SingleMassFracSNII[0]*NSNII ) + SingleMassFracSNIa[0]*NSNIa; #endif /* ejected mass per element */ for (j=2;jnelts+2;j++) #ifdef CHIMIE_AGB SingleEjectedMass[j] = M0*( SingleMassFracDYIN[j]*NDYIN +z[j-2]*SingleMassFracDYIN[1]*NDYIN + SingleMassFracSNII[j]*NSNII +z[j-2]*SingleMassFracSNII[1]*NSNII ) + SingleMassFracSNIa[j]*NSNIa; #else SingleEjectedMass[j] = M0*( SingleMassFracSNII[j]*NSNII +z[j-2]*SingleMassFracSNII[1]*NSNII ) + SingleMassFracSNIa[j]*NSNIa; #endif /* not used */ SingleEjectedMass[1] = -1; } /****************************************************************************************/ /* /* /* /* GADGET ONLY PART /* /* /* /****************************************************************************************/ static double hubble_a, atime, hubble_a2, fac_mu, fac_vsic_fix, a3inv, fac_egy; #ifdef FEEDBACK static double fac_pow; #endif #ifdef PERIODIC static double boxSize, boxHalf; #ifdef LONG_X static double boxSize_X, boxHalf_X; #else #define boxSize_X boxSize #define boxHalf_X boxHalf #endif #ifdef LONG_Y static double boxSize_Y, boxHalf_Y; #else #define boxSize_Y boxSize #define boxHalf_Y boxHalf #endif #ifdef LONG_Z static double boxSize_Z, boxHalf_Z; #else #define boxSize_Z boxSize #define boxHalf_Z boxHalf #endif #endif #if defined(CHIMIE_THERMAL_FEEDBACK) && defined(CHIMIE_COMPUTE_THERMAL_FEEDBACK_ENERGY) void chimie_compute_energy_int(int mode) { int i; double DeltaEgyInt; double Tot_DeltaEgyInt; DeltaEgyInt = 0; Tot_DeltaEgyInt = 0; if (mode==1) { LocalSysState.EnergyInt1 = 0; LocalSysState.EnergyInt2 = 0; } for(i = 0; i < N_gas; i++) { if (P[i].Type==0) { if (mode==1) #ifdef DENSITY_INDEPENDENT_SPH LocalSysState.EnergyInt1 += P[i].Mass * SphP[i].EntropyPred / (GAMMA_MINUS1) * pow(SphP[i].EgyWtDensity*a3inv, GAMMA_MINUS1); #else LocalSysState.EnergyInt1 += P[i].Mass * SphP[i].EntropyPred / (GAMMA_MINUS1) * pow(SphP[i].Density*a3inv, GAMMA_MINUS1); #endif else #ifdef DENSITY_INDEPENDENT_SPH LocalSysState.EnergyInt2 += P[i].Mass * SphP[i].EntropyPred / (GAMMA_MINUS1) * pow(SphP[i].EgyWtDensity*a3inv, GAMMA_MINUS1); #else LocalSysState.EnergyInt2 += P[i].Mass * SphP[i].EntropyPred / (GAMMA_MINUS1) * pow(SphP[i].Density*a3inv, GAMMA_MINUS1); #endif } } if (mode==2) { DeltaEgyInt = LocalSysState.EnergyInt2 - LocalSysState.EnergyInt1; MPI_Reduce(&DeltaEgyInt, &Tot_DeltaEgyInt, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); LocalSysState.EnergyThermalFeedback -= DeltaEgyInt; } } #endif #if defined(CHIMIE_KINETIC_FEEDBACK) && defined(CHIMIE_COMPUTE_KINETIC_FEEDBACK_ENERGY) void chimie_compute_energy_kin(int mode) { int i; double DeltaEgyKin; double Tot_DeltaEgyKin; DeltaEgyKin = 0; Tot_DeltaEgyKin = 0; if (mode==1) { LocalSysState.EnergyKin1 = 0; LocalSysState.EnergyKin2 = 0; } for(i = 0; i < N_gas; i++) { if (P[i].Type==0) { if (mode==1) LocalSysState.EnergyKin1 += 0.5 * P[i].Mass * (P[i].Vel[0]*P[i].Vel[0]+P[i].Vel[1]*P[i].Vel[1]+P[i].Vel[2]*P[i].Vel[2]); else LocalSysState.EnergyKin2 += 0.5 * P[i].Mass * (P[i].Vel[0]*P[i].Vel[0]+P[i].Vel[1]*P[i].Vel[1]+P[i].Vel[2]*P[i].Vel[2]); } } if (mode==2) { DeltaEgyKin = LocalSysState.EnergyKin2 - LocalSysState.EnergyKin1; MPI_Reduce(&DeltaEgyKin, &Tot_DeltaEgyKin, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); LocalSysState.EnergyKineticFeedback -= DeltaEgyKin; } } #endif #ifdef CHIMIE_THERMAL_FEEDBACK void chimie_apply_thermal_feedback(void) { int i; double EgySpec,NewEgySpec,DeltaEntropy; for(i = 0; i < N_gas; i++) { if (P[i].Type==0) { if (SphP[i].DeltaEgySpec > 0) { //printf("(%d) Step=%d i=%08d particle receive feedback\n",ThisTask,All.NumCurrentTiStep,i); /* spec energy at current step (allways compute energy budget based on predicted entropy) */ #ifdef DENSITY_INDEPENDENT_SPH EgySpec = SphP[i].EntropyPred / GAMMA_MINUS1 * pow(SphP[i].EgyWtDensity*a3inv, GAMMA_MINUS1); #else EgySpec = SphP[i].EntropyPred / GAMMA_MINUS1 * pow(SphP[i].Density *a3inv, GAMMA_MINUS1); #endif /* new egyspec */ NewEgySpec = EgySpec + SphP[i].DeltaEgySpec; LocalSysState.EnergyThermalFeedback -= SphP[i].DeltaEgySpec*P[i].Mass; /* new entropy */ #ifdef DENSITY_INDEPENDENT_SPH DeltaEntropy = GAMMA_MINUS1*NewEgySpec/pow(SphP[i].EgyWtDensity*a3inv, GAMMA_MINUS1) - SphP[i].EntropyPred; #else DeltaEntropy = GAMMA_MINUS1*NewEgySpec/pow(SphP[i].Density *a3inv, GAMMA_MINUS1) - SphP[i].EntropyPred; #endif SphP[i].EntropyPred += DeltaEntropy; SphP[i].Entropy += DeltaEntropy; #ifdef DENSITY_INDEPENDENT_SPH SphP[i].EntVarPred = pow(SphP[i].EntropyPred, 1/GAMMA); #endif /* set the adiabatic period for SNIa */ if (SphP[i].NumberOfSNIa>0) SphP[i].SNIaThermalTime = All.Time; /* set the adiabatic period for SNII */ if (SphP[i].NumberOfSNII>0) SphP[i].SNIIThermalTime = All.Time; /* reset variables */ SphP[i].DeltaEgySpec = 0; SphP[i].NumberOfSNIa = 0; SphP[i].NumberOfSNII = 0; } } } } #endif #ifdef CHIMIE_KINETIC_FEEDBACK void chimie_apply_wind(void) { /* apply wind */ int i; double e1,e2; double phi,costh,sinth,vx,vy,vz; for(i = 0; i < N_gas; i++) { if (P[i].Type==0) { if (SphP[i].WindFlag) { phi = get_ChimieKineticFeedback_random_number(P[i].ID)*PI*2.; costh = 1.-2.*get_ChimieKineticFeedback_random_number(P[i].ID+1); sinth = sqrt(1.-pow(costh,2)); vx = All.ChimieWindSpeed*sinth*cos(phi); vy = All.ChimieWindSpeed*sinth*sin(phi); vz = All.ChimieWindSpeed*costh; e1 = 0.5*P[i].Mass * ( SphP[i].VelPred[0]*SphP[i].VelPred[0] + SphP[i].VelPred[1]*SphP[i].VelPred[1] + SphP[i].VelPred[2]*SphP[i].VelPred[2]); P[i].Vel[0] += vx; P[i].Vel[1] += vy; P[i].Vel[2] += vz; SphP[i].VelPred[0] += vx; SphP[i].VelPred[1] += vy; SphP[i].VelPred[2] += vz; e2 = 0.5*P[i].Mass * ( SphP[i].VelPred[0]*SphP[i].VelPred[0] + SphP[i].VelPred[1]*SphP[i].VelPred[1] + SphP[i].VelPred[2]*SphP[i].VelPred[2]); LocalSysState.EnergyKineticFeedback -= e2-e1; SphP[i].WindFlag = 0; } } } } #endif /*! This function is the driver routine for the calculation of chemical evolution */ void chimie(void) { double t0, t1; t0 = second(); /* measure the time for the full chimie computation */ #ifdef CHIMIE_TIMEBET if((All.Time - All.ChimieTimeLastChimie) >= All.ChimieTimeBetChimie) { #endif if (ThisTask==0) printf("Start Chimie computation.\n"); if(All.ComovingIntegrationOn) { /* Factors for comoving integration of hydro */ hubble_a = All.Omega0 / (All.Time * All.Time * All.Time) + (1 - All.Omega0 - All.OmegaLambda) / (All.Time * All.Time) + All.OmegaLambda; hubble_a = All.Hubble * sqrt(hubble_a); hubble_a2 = All.Time * All.Time * hubble_a; fac_mu = pow(All.Time, 3 * (GAMMA - 1) / 2) / All.Time; fac_egy = pow(All.Time, 3 * (GAMMA - 1)); fac_vsic_fix = hubble_a * pow(All.Time, 3 * GAMMA_MINUS1); a3inv = 1 / (All.Time * All.Time * All.Time); atime = All.Time; #ifdef FEEDBACK fac_pow = fac_egy*atime*atime; #endif } else { hubble_a = hubble_a2 = atime = fac_mu = fac_vsic_fix = a3inv = fac_egy = 1.0; #ifdef FEEDBACK fac_pow = 1.0; #endif } stars_density(); /* compute density */ #ifdef CHIMIE_ONE_SN_ONLY if(All.ChimieOneSN==0) /* explode only if not one sn only*/ #endif do_chimie(); /* chimie */ #ifdef CHIMIE_TIMEBET All.ChimieTimeLastChimie += All.ChimieTimeBetChimie; All.ChimieTi_begstep = All.Ti_Current; if (ThisTask==0) printf("Chimie: next Chimie=%g \n",All.ChimieTimeLastChimie); #endif if (ThisTask==0) printf("Chimie computation done.\n"); t1 = second(); All.CPU_Chimie += timediff(t0, t1); #ifdef CHIMIE_TIMEBET } #endif } /*! This function is the driver routine for the calculation of chemical evolution */ void do_chimie(void) { long long ntot, ntotleft; int i, j, k, n, m, ngrp, maxfill, source, ndone; int *nbuffer, *noffset, *nsend_local, *nsend, *numlist, *ndonelist; int level, sendTask, recvTask, nexport, place; double tstart, tend, sumt, sumcomm; double timecomp = 0, timecommsumm = 0, timeimbalance = 0, sumimbalance; int flag_chimie; MPI_Status status; int do_it; int Ti0,Ti1,Ti2; double t1,t2,t01,t02; double tmin,tmax; double minlivetime,maxlivetime; double m1,m2,M0; double NSNIa,NSNII,NDYIN; double NSNIa_tot,NSNII_tot,NDYIN_tot,NSNIa_totlocal,NSNII_totlocal,NDYIN_totlocal; double EgySN,EgySNlocal; double EgySNThermal,EgySNKinetic; int Nchim,Nchimlocal; int Nwind,Nwindlocal; int Nflag,Nflaglocal; int Noldwind,Noldwindlocal; double metals[NELEMENTS]; double FeH; float MinRelMass=1e-3; #ifdef DETAILED_CPU_OUTPUT_IN_CHIMIE double *timecomplist; double *timecommsummlist; double *timeimbalancelist; #endif #ifdef PERIODIC boxSize = All.BoxSize; boxHalf = 0.5 * All.BoxSize; #ifdef LONG_X boxHalf_X = boxHalf * LONG_X; boxSize_X = boxSize * LONG_X; #endif #ifdef LONG_Y boxHalf_Y = boxHalf * LONG_Y; boxSize_Y = boxSize * LONG_Y; #endif #ifdef LONG_Z boxHalf_Z = boxHalf * LONG_Z; boxSize_Z = boxSize * LONG_Z; #endif #endif #ifdef COMPUTE_VELOCITY_DISPERSION double v1m,v2m; #endif /* `NumStUpdate' gives the number of particles on this processor that want a chimie computation */ for(n = 0, NumStUpdate = 0; n < N_gas+N_stars; n++) { #ifndef CHIMIE_TIMEBET if(P[n].Ti_endstep == All.Ti_Current) #else /*nothing*/ /* as this is performed every ChimieTimeBetChimie, do it over all particles */ #endif if(P[n].Type == ST) { m = P[n].StPIdx; #ifdef FOF_SFR if (StP[m].NStars != 0) #endif if ( (P[n].Mass/StP[m].InitialMass) > MinRelMass) NumStUpdate++; } if(P[n].Type == 0) SphP[n].dMass = 0.; } numlist = malloc(NTask * sizeof(int) * NTask); MPI_Allgather(&NumStUpdate, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD); for(i = 0, ntot = 0; i < NTask; i++) ntot += numlist[i]; free(numlist); noffset = malloc(sizeof(int) * NTask); /* offsets of bunches in common list */ nbuffer = malloc(sizeof(int) * NTask); nsend_local = malloc(sizeof(int) * NTask); nsend = malloc(sizeof(int) * NTask * NTask); ndonelist = malloc(sizeof(int) * NTask); i = 0; /* first gas particle, because stars may be hidden among gas particles */ ntotleft = ntot; /* particles left for all tasks together */ NSNIa_tot = 0; NSNII_tot = 0; NDYIN_tot = 0; NSNIa_totlocal = 0; NSNII_totlocal = 0; NDYIN_totlocal = 0; EgySN = 0; EgySNlocal =0; Nchimlocal = 0; Nchim = 0; Nwindlocal = 0; Nwind = 0; Noldwindlocal = 0; Noldwind = 0; Nflaglocal = 0; Nflag = 0; while(ntotleft > 0) { for(j = 0; j < NTask; j++) nsend_local[j] = 0; /* do local particles and prepare export list */ tstart = second(); for(nexport = 0, ndone = 0; i < N_gas+N_stars && nexport < All.BunchSizeChimie - NTask; i++) { /* only active particles and stars */ #ifndef CHIMIE_TIMEBET if((P[i].Ti_endstep == All.Ti_Current)&&(P[i].Type == ST)) #else if(P[i].Type == ST) /* as this is performed every ChimieTimeBetChimie, do it over all particles */ #endif { if(P[i].Type != ST) { printf("P[i].Type != ST, we better stop.\n"); printf("N_gas=%d (type=%d) i=%d (type=%d)\n",N_gas,P[N_gas].Type,i,P[i].Type); printf("Please, check that you do not use PEANOHILBERT\n"); endrun(777001); } m = P[i].StPIdx; #ifdef FOF_SFR if (StP[m].NStars != 0) #endif if ( (P[i].Mass/StP[m].InitialMass) > MinRelMass) { flag_chimie = 0; /******************************************/ /* do chimie */ /******************************************/ /*****************************************************/ /* look if a SN may have explode during the last step /*****************************************************/ /***********************************************/ /***********************************************/ /* set the right table base of the metallicity */ set_table(0); /***********************************************/ /***********************************************/ #ifdef FOF_SFR /* minimum live time for a given metallicity */ minlivetime = star_lifetime(StP[m].Metal[NELEMENTS-1],StP[m].MassMax); /* maximum live time for a given metallicity */ maxlivetime = star_lifetime(StP[m].Metal[NELEMENTS-1],StP[m].MassMin); #else /* minimum live time for a given metallicity */ minlivetime = star_lifetime(StP[m].Metal[NELEMENTS-1],Cp->Mmax*All.MsoltoCMU); /* maximum live time for a given metallicity */ maxlivetime = star_lifetime(StP[m].Metal[NELEMENTS-1],Cp->Mmin*All.MsoltoCMU); #endif if (All.ComovingIntegrationOn) { /* FormationTime on the time line */ Ti0 = log(StP[m].FormationTime/All.TimeBegin) / All.Timebase_interval; /* Beginning of time step on the time line */ #ifndef CHIMIE_TIMEBET Ti1 = P[i].Ti_begstep; #else Ti1 = All.ChimieTi_begstep; #endif /* End of time step on the time line */ Ti2 = All.Ti_Current; #ifdef COSMICTIME t01 = get_cosmictime_difference(Ti0,Ti1); t02 = get_cosmictime_difference(Ti0,Ti2); #endif } else { #ifndef CHIMIE_TIMEBET t1 = All.TimeBegin + (P[i].Ti_begstep * All.Timebase_interval); #else t1 = All.TimeBegin + (All.ChimieTi_begstep * All.Timebase_interval); #endif t2 = All.TimeBegin + (All.Ti_Current * All.Timebase_interval); t01 = t1-StP[m].FormationTime; t02 = t2-StP[m].FormationTime; } /* now treat all cases */ do_it=1; #ifdef CHIMIE_ONE_SN_ONLY if (All.Time<0.1) do_it=0; else m1=m2=1; /* fix to 1 in order numerical problems with lifetime */ #else /* beginning of interval */ if (t01>=minlivetime) if (t01>=maxlivetime) do_it=0; /* nothing to do */ else m2 = star_mass_from_age(StP[m].Metal[NELEMENTS-1],t01); else m2 = Cp->Mmax*All.MsoltoCMU; /* end of interval */ if (t02<=maxlivetime) if (t02<=minlivetime) do_it=0; /* nothing to do */ else m1 = star_mass_from_age(StP[m].Metal[NELEMENTS-1],t02); else m1 = Cp->Mmin*All.MsoltoCMU; #endif //printf("Time=%g t01=%g t02=%g id=%d minlivetime=%g maxlivetime=%g \n",All.Time,t01,t02,P[i].ID,minlivetime,maxlivetime); /* if some of the stars in the SSP explode between t1 and t2 */ if (do_it) { #ifdef CHIMIE_SN_ON_SINGLE_PART printf("(%d) CHIMIE_SN_ON_SINGLE_PART : %d -> %d dist=%g",ThisTask,P[i].ID,StP[m].IDMinDist,StP[m].MinDist); #endif Nchimlocal++; StP[m].Flag = 1; /* mark it as active */ if (m1>m2) { printf("m1=%g (%g Msol) > m2=%g (%g Msol) !!!\n\n",m1,m1*All.CMUtoMsol,m2,m2*All.CMUtoMsol); endrun(777002); } M0 = StP[m].InitialMass; for (k=0;k0) { current_mass = StP[m].MassMax; /* start by the most massive star represented by the stellar particle */ printf("(%d) FOF_SFR : INDIVIDUAL STARS id=%8d %8.4g %8.4g %8d current_mass=%g MassSSP=%g\n",ThisTask,P[i].ID,StP[m].MassMin*All.CMUtoMsol,StP[m].MassMax*All.CMUtoMsol,StP[m].NStars,current_mass*All.CMUtoMsol,P[i].Mass*All.CMUtoMsol,StP[m].MassSSP*All.CMUtoMsol); fprintf(FdFOF_Chimie,"%15g : id=%8d INDIVIDUAL STARS %8.4g %8.4g %8d current_mass=%g MassSSP=%g %g\n",All.Time,P[i].ID,StP[m].MassMin*All.CMUtoMsol,StP[m].MassMax*All.CMUtoMsol,StP[m].NStars,current_mass*All.CMUtoMsol,P[i].Mass*All.CMUtoMsol,StP[m].MassSSP*All.CMUtoMsol); if(current_mass > m2) { printf(" current_mass = %g > m2 = %g, this seems odd (id=%d)!!!\n ",current_mass*All.CMUtoMsol,m2*All.CMUtoMsol,P[i].ID); fprintf(FdFOF_Chimie,"%15g : id=%8d !!! current_mass = %g > m2 = %g, this seems odd !!!\n ",All.Time,P[i].ID,current_mass*All.CMUtoMsol,m2*All.CMUtoMsol); endrun(888029); } while ( (current_mass > m1) && (StP[m].NStars>0)) { /* stop the loop if the mass goes beyond the minimum stellar mass responsible of ejetas */ if (current_mass < Cp->MassMinForEjecta*All.MsoltoCMU) { StP[m].NStars=0; printf("(%d) FOF_SFR : id=%8d current_mass=%g ---> mass below MassMinForEjecta=%g\n",ThisTask,P[i].ID,current_mass*All.CMUtoMsol,Cp->MassMinForEjecta); fprintf(FdFOF_Chimie,"%15g : id=%8d current_mass=%g ---> mass below MassMinForEjecta=%g\n",All.Time,P[i].ID,current_mass*All.CMUtoMsol,Cp->MassMinForEjecta); break; } /**************** SNII ****************/ printf("(%d) FOF_SFR : id=%8d current_mass=%g m1=%g\n",ThisTask,P[i].ID,current_mass*All.CMUtoMsol,m1*All.CMUtoMsol); fprintf(FdFOF_Chimie,"%15g : id=%8d current_mass=%g m1=%g\n",All.Time,P[i].ID,current_mass*All.CMUtoMsol,m1*All.CMUtoMsol); if ( (Cp->SNII_Mmin*All.MsoltoCMU <= current_mass) && (current_mass <= Cp->SNII_Mmax*All.MsoltoCMU) ) { NSNII++; //StP[m].NStars--; printf("(%d) FOF_SFR : id=%8d one SNII\n",ThisTask,P[i].ID); fprintf(FdFOF_Chimie,"%15g : id=%8d one SNII\n",All.Time,P[i].ID); SNII_Total_single_mass_ejection(current_mass,metals); StP[m].TotalEjectedGasMass += SingleEjectedMass[0]; /* gas mass */ for (k=0;kDYIN_Mmin*All.MsoltoCMU <= current_mass) && (current_mass <= Cp->DYIN_Mmax*All.MsoltoCMU) ) { NDYIN++; //StP[m].NStars--; printf("(%d) FOF_SFR : id=%8d one DYIN\n",ThisTask,P[i].ID); fprintf(FdFOF_Chimie,"%15g : id=%8d one DYIN\n",All.Time,P[i].ID); DYIN_Total_single_mass_ejection(current_mass,metals); StP[m].TotalEjectedGasMass += SingleEjectedMass[0]; /* gas mass */ for (k=0;k0) fprintf(FdFOF_Chimie,"%15g : id=%8d NSNIa=%d, m1=%g m2=%g\n",All.Time,P[i].ID,NSNIa,m1/All.MsoltoCMU,m2/All.MsoltoCMU); } /************************************************************************** second case: we are dealing with a particle containing a portion of IMF **************************************************************************/ else { m2 = dmin(m2,StP[m].MassMax); m1 = dmax(m1,StP[m].MassMin); /* number of SNIa */ NSNIa = SNIa_rate(m1,m2)*StP[m].MassSSP; /* number of SNII */ NSNII = SNII_rate(m1,m2)*StP[m].MassSSP; /* number of DYIN */ #ifdef CHIMIE_AGB NDYIN = DYIN_rate(m1,m2)*StP[m].MassSSP; #endif printf("(%d) FOF_SFR : IMF PORTION id=%8d Mass=%g MassSSP=%g NSNII=%g NDYIN=%g NSNIa=%g\n",ThisTask,P[i].ID,P[i].Mass*All.CMUtoMsol,StP[m].MassSSP*All.CMUtoMsol,NSNII,NDYIN,NSNIa); fprintf(FdFOF_Chimie,"%15g : id=%8d IMF PORTION Mass=%g MassSSP=%g NSNII=%g NDYIN=%g NSNIa=%g\n",All.Time,P[i].ID,P[i].Mass*All.CMUtoMsol,StP[m].MassSSP*All.CMUtoMsol,NSNII,NDYIN,NSNIa); #ifdef CHIMIE_MC_SUPERNOVAE double fNSNIa,fNSNII,fNDYIN; /* discretize SNIa */ fNSNIa = NSNIa-floor(NSNIa); NSNIa = floor(NSNIa); if (get_Chimie_random_number(P[i].ID) < fNSNIa) NSNIa = NSNIa+1; /* discretize SNII */ fNSNII = NSNII-floor(NSNII); NSNII = floor(NSNII); //if (get_Chimie_random_number(P[i].ID) < fNSNII) // NSNII = NSNII+1; /* discretize DYIN */ #ifdef CHIMIE_AGB fNDYIN = NDYIN-floor(NDYIN); NDYIN = floor(NDYIN); //if (get_Chimie_random_number(P[i].ID) < fNDYIN) // NDYIN = NDYIN+1; #endif /* CHIMIE_AGB */ #endif /* compute ejectas */ Total_single_mass_ejection(0.5*(m1+m2),metals,NSNII,NSNIa,NDYIN); StP[m].TotalEjectedGasMass = SingleEjectedMass[0]; /* gas mass */ for (k=0;k m2) { printf(" current_mass = %g > m2 = %g, this seems odd (id=%d)!!!\n ",current_mass*All.CMUtoMsol,m2*All.CMUtoMsol,P[i].ID); endrun(888027); } while (current_mass > m1) { /**************** SNII ****************/ if ( (Cp->SNII_Mmin*All.MsoltoCMU <= current_mass) && (current_mass <= Cp->SNII_Mmax*All.MsoltoCMU) ) { NSNII++; SNII_Total_single_mass_ejection(current_mass,metals); StP[m].TotalEjectedGasMass += SingleEjectedMass[0]; /* gas mass */ for (k=0;kDYIN_Mmin*All.MsoltoCMU <= current_mass) && (current_mass <= Cp->DYIN_Mmax*All.MsoltoCMU) ) { NDYIN++; DYIN_Total_single_mass_ejection(current_mass,metals); StP[m].TotalEjectedGasMass += SingleEjectedMass[0]; /* gas mass */ for (k=0;k0) { #ifdef CHIMIE_STATS if (StP[m].TotalEjectedGasMass>0) { fprintf(FdChimieStatsSNs,"%4d %10d %8.4f %8d %10.3g %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f ",ThisTask,All.NumCurrentTiStep,All.Time, P[i].ID,P[i].Mass,NSNII,NSNIa,0.5*(m1+m2)*All.CMUtoMsol,P[i].Pos[0],P[i].Pos[1],P[i].Pos[2]); fprintf(FdChimieStatsSNs,"%10.3g ",StP[m].TotalEjectedGasMass); for (k=0;k0) flag_chimie=1; /* compute injected energy */ StP[m].TotalEjectedEgySpec = All.ChimieSupernovaEnergy* (NSNIa + NSNII) /StP[m].TotalEjectedGasMass; StP[m].NumberOfSNIa = NSNIa; StP[m].NumberOfSNII = NSNII; EgySNlocal += All.ChimieSupernovaEnergy* (NSNIa + NSNII); NSNIa_totlocal += NSNIa; NSNII_totlocal += NSNII; NDYIN_totlocal += NDYIN; /* correct mass particle */ P[i].Mass = P[i].Mass-StP[m].TotalEjectedGasMass; if(P[i].Mass<0) endrun(777023); } else { printf("CHIMIE : mass wants to be less than zero...\n"); printf("CHIMIE : ID=%d Mass=%g StP[m].TotalEjectedGasMass=%g\n",P[i].ID,P[i].Mass*All.CMUtoMsol,StP[m].TotalEjectedGasMass*All.CMUtoMsol); } } /******************************************/ /* end do chimie */ /******************************************/ ndone++; if (flag_chimie) { for(j = 0; j < NTask; j++) Exportflag[j] = 0; chimie_evaluate(i, 0); for(j = 0; j < NTask; j++) { if(Exportflag[j]) { for(k = 0; k < 3; k++) { ChimieDataIn[nexport].Pos[k] = P[i].Pos[k]; ChimieDataIn[nexport].Vel[k] = P[i].Vel[k]; } ChimieDataIn[nexport].ID = P[i].ID; ChimieDataIn[nexport].Timestep = P[i].Ti_endstep - P[i].Ti_begstep; ChimieDataIn[nexport].Hsml = StP[m].Hsml; ChimieDataIn[nexport].Density = StP[m].Density; ChimieDataIn[nexport].Y = StP[m].Y; #ifdef CHIMIE_KINETIC_FEEDBACK ChimieDataIn[nexport].NgbMass = StP[m].NgbMass; #endif ChimieDataIn[nexport].TotalEjectedGasMass = StP[m].TotalEjectedGasMass; for(k = 0; k < NELEMENTS; k++) ChimieDataIn[nexport].TotalEjectedEltMass[k] = StP[m].TotalEjectedEltMass[k]; ChimieDataIn[nexport].TotalEjectedEgySpec = StP[m].TotalEjectedEgySpec; ChimieDataIn[nexport].NumberOfSNIa = StP[m].NumberOfSNIa; ChimieDataIn[nexport].NumberOfSNII = StP[m].NumberOfSNII; #ifdef WITH_ID_IN_HYDRA ChimieDataIn[nexport].ID = P[i].ID; #endif #if CHIMIE_EJECTA_RADIUS == 2 ChimieDataIn[nexport].Pressure = StP[m].Pressure; #endif ChimieDataIn[nexport].Index = i; ChimieDataIn[nexport].Task = j; nexport++; nsend_local[j]++; } } } } } } tend = second(); timecomp += timediff(tstart, tend); qsort(ChimieDataIn, nexport, sizeof(struct chimiedata_in), chimie_compare_key); for(j = 1, noffset[0] = 0; j < NTask; j++) noffset[j] = noffset[j - 1] + nsend_local[j - 1]; tstart = second(); MPI_Allgather(nsend_local, NTask, MPI_INT, nsend, NTask, MPI_INT, MPI_COMM_WORLD); tend = second(); timeimbalance += timediff(tstart, tend); /* now do the particles that need to be exported */ for(level = 1; level < (1 << PTask); level++) { tstart = second(); for(j = 0; j < NTask; j++) nbuffer[j] = 0; for(ngrp = level; ngrp < (1 << PTask); ngrp++) { maxfill = 0; for(j = 0; j < NTask; j++) { if((j ^ ngrp) < NTask) if(maxfill < nbuffer[j] + nsend[(j ^ ngrp) * NTask + j]) maxfill = nbuffer[j] + nsend[(j ^ ngrp) * NTask + j]; } if(maxfill >= All.BunchSizeChimie) break; sendTask = ThisTask; recvTask = ThisTask ^ ngrp; if(recvTask < NTask) { if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0) { /* get the particles */ MPI_Sendrecv(&ChimieDataIn[noffset[recvTask]], nsend_local[recvTask] * sizeof(struct chimiedata_in), MPI_BYTE, recvTask, TAG_CHIMIE_A, &ChimieDataGet[nbuffer[ThisTask]], nsend[recvTask * NTask + ThisTask] * sizeof(struct chimiedata_in), MPI_BYTE, recvTask, TAG_CHIMIE_A, MPI_COMM_WORLD, &status); } } for(j = 0; j < NTask; j++) if((j ^ ngrp) < NTask) nbuffer[j] += nsend[(j ^ ngrp) * NTask + j]; } tend = second(); timecommsumm += timediff(tstart, tend); /* now do the imported particles */ tstart = second(); for(j = 0; j < nbuffer[ThisTask]; j++) chimie_evaluate(j, 1); tend = second(); timecomp += timediff(tstart, tend); /* do a block to measure imbalance */ tstart = second(); MPI_Barrier(MPI_COMM_WORLD); tend = second(); timeimbalance += timediff(tstart, tend); /* get the result */ tstart = second(); for(j = 0; j < NTask; j++) nbuffer[j] = 0; for(ngrp = level; ngrp < (1 << PTask); ngrp++) { maxfill = 0; for(j = 0; j < NTask; j++) { if((j ^ ngrp) < NTask) if(maxfill < nbuffer[j] + nsend[(j ^ ngrp) * NTask + j]) maxfill = nbuffer[j] + nsend[(j ^ ngrp) * NTask + j]; } if(maxfill >= All.BunchSizeChimie) break; sendTask = ThisTask; recvTask = ThisTask ^ ngrp; if(recvTask < NTask) { if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0) { /* send the results */ MPI_Sendrecv(&ChimieDataResult[nbuffer[ThisTask]], nsend[recvTask * NTask + ThisTask] * sizeof(struct chimiedata_out), MPI_BYTE, recvTask, TAG_CHIMIE_B, &ChimieDataPartialResult[noffset[recvTask]], nsend_local[recvTask] * sizeof(struct chimiedata_out), MPI_BYTE, recvTask, TAG_CHIMIE_B, MPI_COMM_WORLD, &status); /* add the result to the particles */ for(j = 0; j < nsend_local[recvTask]; j++) { source = j + noffset[recvTask]; place = ChimieDataIn[source].Index; // for(k = 0; k < 3; k++) // SphP[place].HydroAccel[k] += HydroDataPartialResult[source].Acc[k]; // // SphP[place].DtEntropy += HydroDataPartialResult[source].DtEntropy; //#ifdef FEEDBACK // SphP[place].DtEgySpecFeedback += HydroDataPartialResult[source].DtEgySpecFeedback; //#endif // if(SphP[place].MaxSignalVel < HydroDataPartialResult[source].MaxSignalVel) // SphP[place].MaxSignalVel = HydroDataPartialResult[source].MaxSignalVel; //#ifdef COMPUTE_VELOCITY_DISPERSION // for(k = 0; k < VELOCITY_DISPERSION_SIZE; k++) // SphP[place].VelocityDispersion[k] += HydroDataPartialResult[source].VelocityDispersion[k]; //#endif } } } for(j = 0; j < NTask; j++) if((j ^ ngrp) < NTask) nbuffer[j] += nsend[(j ^ ngrp) * NTask + j]; } tend = second(); timecommsumm += timediff(tstart, tend); level = ngrp - 1; } MPI_Allgather(&ndone, 1, MPI_INT, ndonelist, 1, MPI_INT, MPI_COMM_WORLD); for(j = 0; j < NTask; j++) ntotleft -= ndonelist[j]; } free(ndonelist); free(nsend); free(nsend_local); free(nbuffer); free(noffset); /* do final operations on results */ tstart = second(); for(i = 0; i < N_gas; i++) { if (P[i].Type==0) { P[i].Mass += SphP[i].dMass; SphP[i].dMass = 0.; } } tend = second(); timecomp += timediff(tstart, tend); /* collect some timing information */ MPI_Reduce(&timecomp, &sumt, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); MPI_Reduce(&timecommsumm, &sumcomm, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); MPI_Reduce(&timeimbalance, &sumimbalance, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); if(ThisTask == 0) { All.CPU_ChimieCompWalk += sumt / NTask; All.CPU_ChimieCommSumm += sumcomm / NTask; All.CPU_ChimieImbalance += sumimbalance / NTask; } #ifdef DETAILED_CPU_OUTPUT_IN_CHIMIE numlist = malloc(sizeof(int) * NTask); timecomplist = malloc(sizeof(double) * NTask); timecommsummlist = malloc(sizeof(double) * NTask); timeimbalancelist = malloc(sizeof(double) * NTask); MPI_Gather(&NumStUpdate, 1, MPI_INT, numlist, 1, MPI_INT, 0, MPI_COMM_WORLD); MPI_Gather(&timecomp, 1, MPI_DOUBLE, timecomplist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Gather(&timecommsumm, 1, MPI_DOUBLE, timecommsummlist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Gather(&timeimbalance, 1, MPI_DOUBLE, timeimbalancelist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); if(ThisTask == 0) { fprintf(FdTimings, "\n chimie\n\n"); fprintf(FdTimings, "Nupdate "); for (i=0;i= (All.Time-All.ChimieWindTime)) Nwindlocal++; //else // if (SphP[i].WindTime > All.TimeBegin-2*All.ChimieWindTime) // Noldwindlocal++; if (SphP[i].WindFlag) Nflaglocal++; } } MPI_Reduce(&Nwindlocal, &Nwind, 1, MPI_INT , MPI_SUM, 0, MPI_COMM_WORLD); MPI_Reduce(&Noldwindlocal, &Noldwind, 1, MPI_INT , MPI_SUM, 0, MPI_COMM_WORLD); MPI_Allreduce(&Nflaglocal, &Nflag, 1, MPI_INT , MPI_SUM, MPI_COMM_WORLD); #else EgySNKinetic = 0; #endif /* write some info */ if (ThisTask==0) { fprintf(FdChimie, "%15g %10d %15g %15g %15g %15g %15g %10d %10d %10d\n",All.Time,Nchim,NSNIa_tot,NSNII_tot,EgySN,EgySNThermal,EgySNKinetic,Nwind,Noldwind,Nflag); fflush(FdChimie); } /* this is no longer used */ // if (Nflag>0) // { // SetMinTimeStepForActives=1; // if (ThisTask==0) // fprintf(FdLog,"%g : !!! set min timestep for active particles !!!\n",All.Time); // } #ifdef CHIMIE_ONE_SN_ONLY if (EgySN>0) All.ChimieOneSN=1; MPI_Bcast(&All.ChimieOneSN, 1, MPI_INT, 0, MPI_COMM_WORLD); #endif } /*! This function is the 'core' of the Chemie computation. A target * particle is specified which may either be local, or reside in the * communication buffer. */ void chimie_evaluate(int target, int mode) { int j, n, startnode, numngb,numngb_inbox,k; FLOAT *pos,*vel; //FLOAT *vel; //FLOAT mass; double h, h2; double acc[3]; double dx, dy, dz; double wk, r, r2, u=0; double hinv=1, hinv3; int target_stp; double density; double Y; #ifdef CHIMIE_KINETIC_FEEDBACK double ngbmass; double p; #endif double aij; double ejectedGasMass; double ejectedEltMass[NELEMENTS]; double ejectedEgySpec; double NumberOfSNIa; double NumberOfSNII; double mass_k; double NewMass; double fv,vi2,vj2; double EgySpec,NewEgySpec; double DeltaEntropy; double DeltaVel[3]; #ifndef LONGIDS unsigned int id; /*!< particle identifier */ #else unsigned long long id; /*!< particle identifier */ #endif #if CHIMIE_EJECTA_RADIUS == 2 double rho,E51,n0,P04,pressure,P0; #endif if(mode == 0) { pos = P[target].Pos; vel = P[target].Vel; id = P[target].ID; target_stp = P[target].StPIdx; h = StP[target_stp].Hsml; density = StP[target_stp].Density; Y = StP[target_stp].Y; #ifdef CHIMIE_KINETIC_FEEDBACK ngbmass = StP[target_stp].NgbMass; #endif ejectedGasMass = StP[target_stp].TotalEjectedGasMass; for(k=0;k boxHalf_X) dx -= boxSize_X; if(dx < -boxHalf_X) dx += boxSize_X; if(dy > boxHalf_Y) dy -= boxSize_Y; if(dy < -boxHalf_Y) dy += boxSize_Y; if(dz > boxHalf_Z) dz -= boxSize_Z; if(dz < -boxHalf_Z) dz += boxSize_Z; #endif r2 = dx * dx + dy * dy + dz * dz; if(r2 < h2) { numngb++; r = sqrt(r2); u = r * hinv; if(u < 0.5) { wk = hinv3 * (KERNEL_COEFF_1 + KERNEL_COEFF_2 * (u - 1) * u * u); } else { wk = hinv3 * KERNEL_COEFF_5 * (1.0 - u) * (1.0 - u) * (1.0 - u); } /* normalisation using mass */ //aij = P[j].Mass*wk/density; /* normalisation using volume */ /* !!! si on utilise, il faut stoquer une nouvelle variable : OldDensity, car density est modifié plus bas... */ //aij = P[j].Mass/SphP[j].Density*wk/volume; /* !!! we shoud better use an OldDensity here, as .Density may slighly change after */ #if CHIMIE_ELTS_WEIGHTING == 0 aij = P[j].Mass * wk; // 0: mass + kernel weighting (default, Revaz & Jablonka 2012) #elif CHIMIE_ELTS_WEIGHTING == 1 aij = P[j].Mass/SphP[j].Density * wk; // 1: volume + kernel weighting #elif CHIMIE_ELTS_WEIGHTING == 2 aij = wk; // 2: number + kernel weighting #elif CHIMIE_ELTS_WEIGHTING == 3 aij = P[j].Mass; // 3: mass weighting #elif CHIMIE_ELTS_WEIGHTING == 4 aij = P[j].Mass/SphP[j].Density; // 4: volume weighting #elif CHIMIE_ELTS_WEIGHTING == 5 aij = 1; // 5: number weighting #endif aij = aij/Y; // normalisation #ifdef CHIMIE_STATS /* before the interaction */ fprintf(FdChimieStatsGas,"%4d %10d %8.4f %8d %8d %10.3g %10.3g %10.3g %10.3g ",ThisTask,All.NumCurrentTiStep,All.Time, P[j].ID,id,aij,P[j].Pos[0],P[j].Pos[1],P[j].Pos[2]); fprintf(FdChimieStatsGas,"%10.3g ",P[j].Mass); #ifdef CHIMIE_SMOOTH_METALS for (k=0;k= 0); /* Now collect the result at the right place */ if(mode == 0) { // for(k = 0; k < 3; k++) // SphP[target].HydroAccel[k] = acc[k]; // SphP[target].DtEntropy = dtEntropy; //#ifdef FEEDBACK // SphP[target].DtEgySpecFeedback = dtEgySpecFeedback; //#endif // SphP[target].MaxSignalVel = maxSignalVel; //#ifdef COMPUTE_VELOCITY_DISPERSION // for(k = 0; k < VELOCITY_DISPERSION_SIZE; k++) // SphP[target].VelocityDispersion[k] = VelocityDispersion[k]; //#endif } else { // for(k = 0; k < 3; k++) // HydroDataResult[target].Acc[k] = acc[k]; // HydroDataResult[target].DtEntropy = dtEntropy; //#ifdef FEEDBACK // HydroDataResult[target].DtEgySpecFeedback = dtEgySpecFeedback; //#endif // HydroDataResult[target].MaxSignalVel = maxSignalVel; //#ifdef COMPUTE_VELOCITY_DISPERSION // for(k = 0; k < VELOCITY_DISPERSION_SIZE; k++) // HydroDataResult[target].VelocityDispersion[k] = VelocityDispersion[k]; //#endif } } /*! This is a comparison kernel for a sort routine, which is used to group * particles that are going to be exported to the same CPU. */ int chimie_compare_key(const void *a, const void *b) { if(((struct chimiedata_in *) a)->Task < (((struct chimiedata_in *) b)->Task)) return -1; if(((struct chimiedata_in *) a)->Task > (((struct chimiedata_in *) b)->Task)) return +1; return 0; } /****************************************************************************************/ /* /* /* /* PYTHON INTERFACE /* /* /* /****************************************************************************************/ #ifdef PYCHEM static PyObject * chemistry_CodeUnits_to_SolarMass_Factor(PyObject *self, PyObject *args) { return Py_BuildValue("d",All.CMUtoMsol); } static PyObject * chemistry_SolarMass_to_CodeUnits_Factor(PyObject *self, PyObject *args) { return Py_BuildValue("d",All.MsoltoCMU); } static PyObject * chemistry_SetVerbosityOn(PyObject *self, PyObject *args) { verbose=1; return Py_BuildValue("i",0); } static PyObject * chemistry_SetVerbosityOff(PyObject *self, PyObject *args) { verbose=0; return Py_BuildValue("i",0); } static PyObject * chemistry_InitDefaultParameters(void) { /* list of Gadget parameters */ /* System of units */ All.UnitLength_in_cm = 3.085e+21; /* 1.0 kpc */ All.UnitMass_in_g = 1.989e+43; /* 1.0e10 solar masses */ All.UnitVelocity_in_cm_per_s = 20725573.785998672; /* 207 km/sec */ All.GravityConstantInternal = 0; All.HubbleParam = 1; /* other usefull constants */ All.UnitTime_in_s = All.UnitLength_in_cm / All.UnitVelocity_in_cm_per_s; All.UnitTime_in_Megayears=All.UnitTime_in_s / SEC_PER_MEGAYEAR; All.CMUtoMsol = All.UnitMass_in_g/SOLAR_MASS/All.HubbleParam; /* convertion factor from Code Mass Unit to Solar Mass */ All.MsoltoCMU = 1/All.CMUtoMsol; /* convertion factor from Solar Mass to Code Mass Unit */ return Py_BuildValue("i",1); } static PyObject * SetParameters(PyObject *dict) { PyObject *key; PyObject *value; int ivalue; float fvalue; double dvalue; /* check that it is a PyDictObject */ if(!PyDict_Check(dict)) { PyErr_SetString(PyExc_AttributeError, "argument is not a dictionary."); return NULL; } if (PyDict_Size(dict)==0) return Py_BuildValue("i",0); Py_ssize_t pos=0; while(PyDict_Next(dict,&pos,&key,&value)) { if(PyString_Check(key)) { /* System of units */ if(strcmp(PyString_AsString(key), "UnitLength_in_cm")==0) { if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value)) All.UnitLength_in_cm = PyFloat_AsDouble(value); } if(strcmp(PyString_AsString(key), "UnitMass_in_g")==0) { if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value)) All.UnitMass_in_g = PyFloat_AsDouble(value); } if(strcmp(PyString_AsString(key), "UnitVelocity_in_cm_per_s")==0) { if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value)) All.UnitVelocity_in_cm_per_s = PyFloat_AsDouble(value); } if(strcmp(PyString_AsString(key), "HubbleParam")==0) { if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value)) All.HubbleParam = PyFloat_AsDouble(value); } if(strcmp(PyString_AsString(key), "GravityConstantInternal")==0) { if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value)) All.GravityConstantInternal = PyFloat_AsDouble(value); } } } /* other usefull constants */ All.UnitTime_in_s = All.UnitLength_in_cm / All.UnitVelocity_in_cm_per_s; All.UnitTime_in_Megayears=All.UnitTime_in_s / SEC_PER_MEGAYEAR; All.CMUtoMsol = All.UnitMass_in_g/SOLAR_MASS/All.HubbleParam; /* convertion factor from Code Mass Unit to Solar Mass */ All.MsoltoCMU = 1/All.CMUtoMsol; /* convertion factor from Solar Mass to Code Mass Unit */ return Py_BuildValue("i",1); } static PyObject * chemistry_SetParameters(PyObject *self, PyObject *args) { PyObject *dict; /* here, we can have either arguments or dict directly */ if(PyDict_Check(args)) { dict = args; } else { if (! PyArg_ParseTuple(args, "O",&dict)) return NULL; } SetParameters(dict); return Py_BuildValue("i",1); } static PyObject * chemistry_GetParameters(void) { PyObject *dict; PyObject *key; PyObject *value; dict = PyDict_New(); /* System of units */ key = PyString_FromString("UnitLength_in_cm"); value = PyFloat_FromDouble(All.UnitLength_in_cm); PyDict_SetItem(dict,key,value); key = PyString_FromString("UnitMass_in_g"); value = PyFloat_FromDouble(All.UnitMass_in_g); PyDict_SetItem(dict,key,value); key = PyString_FromString("UnitVelocity_in_cm_per_s"); value = PyFloat_FromDouble(All.UnitVelocity_in_cm_per_s); PyDict_SetItem(dict,key,value); key = PyString_FromString("HubbleParam"); value = PyFloat_FromDouble(All.HubbleParam); PyDict_SetItem(dict,key,value); key = PyString_FromString("GravityConstantInternal"); value = PyFloat_FromDouble(All.GravityConstantInternal); PyDict_SetItem(dict,key,value); return Py_BuildValue("O",dict); } /*********************************/ /* */ /*********************************/ static PyObject * chemistry_init_chimie(PyObject *self, PyObject *args, PyObject *kwds) { int NumberOfTables=1; int DefaultTable=0; PyObject *paramsDict=NULL; paramsDict= PyDict_New(); //PyObject *filename; //if (! PyArg_ParseTuple(args, "Oii",&filename,&NumberOfTables,&DefaultTable)) // { // PyErr_SetString(PyExc_ValueError,"init_chimie, error in parsing."); // return NULL; // } static char *kwlist[] = {"filename","NumberOfTables","DefaultTable","params", NULL}; PyObject *filename=PyString_FromString("chimie.yr.dat"); /* this fails with python2.6, I do not know why ??? */ if (! PyArg_ParseTupleAndKeywords(args, kwds, "|OiiO",kwlist,&filename,&NumberOfTables,&DefaultTable,¶msDict)) { PyErr_SetString(PyExc_ValueError,"init_chimie, error in parsing arguments."); return NULL; } if (!PyString_Check(filename)) { PyErr_SetString(PyExc_ValueError,"Argument must be a string."); return NULL; } /* copy filename */ All.ChimieParameterFile = PyString_AsString(filename); /* set number of tables */ All.ChimieNumberOfParameterFiles = NumberOfTables; /* check if the file exists */ if(!(fopen(All.ChimieParameterFile, "r"))) { PyErr_SetString(PyExc_ValueError,"The parameter file does not exists."); return NULL; } /* use default parameters */ chemistry_InitDefaultParameters(); /* check if units are given */ /* check that it is a PyDictObject */ if(!PyDict_Check(paramsDict)) { PyErr_SetString(PyExc_AttributeError, "argument is not a dictionary."); return NULL; } else { SetParameters(paramsDict); } init_chimie(); /* by default, set the first one */ set_table(DefaultTable); return Py_BuildValue("O",Py_None); } /*********************************/ /* */ /*********************************/ static PyObject * chemistry_info(PyObject *self, PyObject *args, PyObject *kwds) { int DefaultTable=0; static char *kwlist[] = {"DefaultTable", NULL}; /* this fails with python2.6, I do not know why ??? */ if (! PyArg_ParseTupleAndKeywords(args, kwds, "|i",kwlist,&DefaultTable)) { PyErr_SetString(PyExc_ValueError,"init_chimie, error in parsing arguments."); return NULL; } info(0); return Py_BuildValue("O",Py_None); } /*********************************/ /* */ /*********************************/ static PyObject * chemistry_set_table(PyObject *self, PyObject *args, PyObject *kwds) { int i; if (! PyArg_ParseTuple(args, "i",&i)) return PyString_FromString("error"); /* set the table */ set_table(i); return Py_BuildValue("d",0); } /*********************************/ /* */ /*********************************/ static PyObject * chemistry_get_imf(self, args) PyObject *self; PyObject *args; { PyArrayObject *m,*imf; int i; if (! PyArg_ParseTuple(args, "O",&m)) return PyString_FromString("error"); m = TO_DOUBLE(m); /* create an output */ imf = (PyArrayObject *) PyArray_SimpleNew(m->nd,m->dimensions,PyArray_DOUBLE); //printf("--> %g\n",Cp->bs[0]); //for (i=0;in;i++) // printf("%g %g\n",Cp->ms[i],Cp->as[i]); for(i = 0; i < m->dimensions[0]; i++) { *(double *)(imf->data + i*(imf->strides[0])) = get_imf(*(double *)(m->data + i*(m->strides[0]))); } return PyArray_Return(imf); } /*********************************/ /* */ /*********************************/ static PyObject * chemistry_get_imf_M(self, args) PyObject *self; PyObject *args; { PyArrayObject *m1,*m2,*imf; int i; if (! PyArg_ParseTuple(args, "OO",&m1,&m2)) return PyString_FromString("error"); m1 = TO_DOUBLE(m1); m2 = TO_DOUBLE(m2); /* create an output */ imf = (PyArrayObject *) PyArray_SimpleNew(m1->nd,m1->dimensions,PyArray_DOUBLE); for(i = 0; i < imf->dimensions[0]; i++) { *(double *)(imf->data + i*(imf->strides[0])) = get_imf_M( *(double *)(m1->data + i*(m1->strides[0])), *(double *)(m2->data + i*(m2->strides[0])) ); } return PyArray_Return(imf); } /*********************************/ /* */ /*********************************/ static PyObject * chemistry_get_imf_N(self, args) PyObject *self; PyObject *args; { PyArrayObject *m1,*m2,*imf; int i; if (! PyArg_ParseTuple(args, "OO",&m1,&m2)) return PyString_FromString("error"); m1 = TO_DOUBLE(m1); m2 = TO_DOUBLE(m2); /* create an output */ imf = (PyArrayObject *) PyArray_SimpleNew(m1->nd,m1->dimensions,PyArray_DOUBLE); for(i = 0; i < imf->dimensions[0]; i++) { *(double *)(imf->data + i*(imf->strides[0])) = get_imf_N( *(double *)(m1->data + i*(m1->strides[0])), *(double *)(m2->data + i*(m2->strides[0])) ); } return PyArray_Return(imf); } /*********************************/ /* */ /*********************************/ static PyObject * chemistry_star_lifetime(self, args) PyObject *self; PyObject *args; { /* z is the mass fraction of metals, ie, the metallicity */ /* m is the star mass in code unit */ /* Return t in time unit */ double time,z,m; if (!PyArg_ParseTuple(args, "dd", &z, &m)) return NULL; time = star_lifetime(z,m); return Py_BuildValue("d",time); } static PyObject * chemistry_star_mass_from_age(self, args) PyObject *self; PyObject *args; { /* t : life time (in code unit) */ /* return the stellar mass (in code unit) that has a lifetime equal to t */ double time,z,m; if (!PyArg_ParseTuple(args, "dd", &z, &time)) return NULL; m = star_mass_from_age(z,time); return Py_BuildValue("d",m); } static PyObject * chemistry_DYIN_rate(self, args) PyObject *self; PyObject *args; { double m1,m2; double RDYIN; /* parse arguments */ if (!PyArg_ParseTuple(args, "dd", &m1,&m2)) return NULL; RDYIN = DYIN_rate(m1,m2); return Py_BuildValue("d",RDYIN); } static PyObject * chemistry_SNII_rate(self, args) PyObject *self; PyObject *args; { double m1,m2; double RSNII; /* parse arguments */ if (!PyArg_ParseTuple(args, "dd", &m1,&m2)) return NULL; RSNII = SNII_rate(m1,m2); return Py_BuildValue("d",RSNII); } static PyObject * chemistry_SNIa_rate(self, args) PyObject *self; PyObject *args; { double m1,m2; double RSNIa; /* parse arguments */ if (!PyArg_ParseTuple(args, "dd", &m1,&m2)) return NULL; RSNIa = SNIa_rate(m1,m2); return Py_BuildValue("d",RSNIa); } static PyObject * chemistry_SNIa_single_rate(self, args) PyObject *self; PyObject *args; { double m; double RSNIa; /* parse arguments */ if (!PyArg_ParseTuple(args, "d", &m)) return NULL; RSNIa = SNIa_single_rate(m); return Py_BuildValue("d",RSNIa); } static PyObject * chemistry_DYIN_mass_ejection(self, args) PyObject *self; PyObject *args; { double m1,m2; PyArrayObject *ArrMassDYIN; npy_intp ld[1]; int i; /* parse arguments */ if (!PyArg_ParseTuple(args, "dd", &m1,&m2)) return NULL; /* create output array */ ld[0]= Cp->nelts+2; ArrMassDYIN = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* compute dying stars ejection */ DYIN_mass_ejection(m1,m2); /* import values */ for (i=0;inelts+2;i++) *(double *)(ArrMassDYIN->data + (i)*(ArrMassDYIN->strides[0])) = MassFracDYIN[i]; /* convert in array */ return PyArray_Return(ArrMassDYIN); } static PyObject * chemistry_DYIN_mass_ejection_Z(self, args) PyObject *self; PyObject *args; { double m1,m2,Z; PyArrayObject *ArrMassDYIN; npy_intp ld[1]; int i; /* parse arguments */ if (!PyArg_ParseTuple(args, "ddd", &m1,&m2,&Z)) return NULL; /* create output array */ ld[0]= Cp->nelts+2; ArrMassDYIN = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* compute dying stars ejection */ DYIN_mass_ejection_Z(m1,m2, Z); /* import values */ for (i=0;inelts+2;i++) *(double *)(ArrMassDYIN->data + (i)*(ArrMassDYIN->strides[0])) = MassFracDYIN[i]; /* convert in array */ return PyArray_Return(ArrMassDYIN); } static PyObject * chemistry_DYIN_single_mass_ejection(self, args) PyObject *self; PyObject *args; { double m1; PyArrayObject *ArrMassDYIN; npy_intp ld[1]; int i; /* parse arguments */ if (!PyArg_ParseTuple(args, "d", &m1)) return NULL; /* create output array */ ld[0]= Cp->nelts+2; ArrMassDYIN = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* compute SN ejection */ DYIN_single_mass_ejection(m1); /* import values */ for (i=0;inelts+2;i++) *(double *)(ArrMassDYIN->data + (i)*(ArrMassDYIN->strides[0])) = SingleMassFracDYIN[i]; /* convert in array */ return PyArray_Return(ArrMassDYIN); } static PyObject * chemistry_DYIN_single_mass_ejection_Z(self, args) PyObject *self; PyObject *args; { double m1,Z; PyArrayObject *ArrMassDYIN; npy_intp ld[1]; int i; /* parse arguments */ if (!PyArg_ParseTuple(args, "dd", &m1, &Z)) return NULL; /* create output array */ ld[0]= Cp->nelts+2; ArrMassDYIN = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* compute SN ejection */ DYIN_single_mass_ejection_Z(m1,Z); /* import values */ for (i=0;inelts+2;i++) *(double *)(ArrMassDYIN->data + (i)*(ArrMassDYIN->strides[0])) = SingleMassFracDYIN[i]; /* convert in array */ return PyArray_Return(ArrMassDYIN); } static PyObject * chemistry_SNII_mass_ejection(self, args) PyObject *self; PyObject *args; { double m1,m2; PyArrayObject *ArrMassSNII; npy_intp ld[1]; int i; /* parse arguments */ if (!PyArg_ParseTuple(args, "dd", &m1,&m2)) return NULL; /* create output array */ ld[0]= Cp->nelts+2; ArrMassSNII = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* compute SN ejection */ SNII_mass_ejection(m1,m2); /* import values */ for (i=0;inelts+2;i++) *(double *)(ArrMassSNII->data + (i)*(ArrMassSNII->strides[0])) = MassFracSNII[i]; /* convert in array */ return PyArray_Return(ArrMassSNII); } static PyObject * chemistry_SNII_mass_ejection_Z(self, args) PyObject *self; PyObject *args; { double m1,m2,Z; PyArrayObject *ArrMassSNII; npy_intp ld[1]; int i; /* parse arguments */ if (!PyArg_ParseTuple(args, "ddd", &m1,&m2,&Z)) return NULL; /* create output array */ ld[0]= Cp->nelts+2; ArrMassSNII = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* compute SNII stars ejection */ SNII_mass_ejection_Z(m1,m2, Z); /* import values */ for (i=0;inelts+2;i++) *(double *)(ArrMassSNII->data + (i)*(ArrMassSNII->strides[0])) = MassFracSNII[i]; /* convert in array */ return PyArray_Return(ArrMassSNII); } static PyObject * chemistry_SNII_single_mass_ejection(self, args) PyObject *self; PyObject *args; { double m1; PyArrayObject *ArrMassSNII; npy_intp ld[1]; int i; /* parse arguments */ if (!PyArg_ParseTuple(args, "d", &m1)) return NULL; /* create output array */ ld[0]= Cp->nelts+2; ArrMassSNII = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* compute SN ejection */ SNII_single_mass_ejection(m1); /* import values */ for (i=0;inelts+2;i++) *(double *)(ArrMassSNII->data + (i)*(ArrMassSNII->strides[0])) = SingleMassFracSNII[i]; /* convert in array */ return PyArray_Return(ArrMassSNII); } static PyObject * chemistry_SNII_single_mass_ejection_Z(self, args) PyObject *self; PyObject *args; { double m1,Z; PyArrayObject *ArrMassSNII; npy_intp ld[1]; int i; /* parse arguments */ if (!PyArg_ParseTuple(args, "dd", &m1, &Z)) return NULL; /* create output array */ ld[0]= Cp->nelts+2; ArrMassSNII = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* compute SN ejection */ SNII_single_mass_ejection_Z(m1,Z); /* import values */ for (i=0;inelts+2;i++) *(double *)(ArrMassSNII->data + (i)*(ArrMassSNII->strides[0])) = SingleMassFracSNII[i]; /* convert in array */ return PyArray_Return(ArrMassSNII); } static PyObject * chemistry_SNIa_mass_ejection(self, args) PyObject *self; PyObject *args; { double m1,m2; PyArrayObject *ArrMassSNIa; npy_intp ld[1]; int i; /* parse arguments */ if (!PyArg_ParseTuple(args, "dd", &m1,&m2)) return NULL; /* create output array */ ld[0]= Cp->nelts+2; ArrMassSNIa = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* compute SN ejection */ SNIa_mass_ejection(m1,m2); /* import values */ for (i=0;inelts+2;i++) *(double *)(ArrMassSNIa->data + (i)*(ArrMassSNIa->strides[0])) = MassFracSNIa[i]; /* convert in array */ return PyArray_Return(ArrMassSNIa); } static PyObject * chemistry_SNIa_single_mass_ejection(self, args) PyObject *self; PyObject *args; { double m1; PyArrayObject *ArrMassSNIa; npy_intp ld[1]; int i; /* parse arguments */ if (!PyArg_ParseTuple(args, "d", &m1)) return NULL; /* create output array */ ld[0]= Cp->nelts+2; ArrMassSNIa = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* compute SN ejection */ SNIa_single_mass_ejection(m1); /* import values */ for (i=0;inelts+2;i++) *(double *)(ArrMassSNIa->data + (i)*(ArrMassSNIa->strides[0])) = SingleMassFracSNIa[i]; /* convert in array */ return PyArray_Return(ArrMassSNIa); } static PyObject * chemistry_Total_mass_ejection(self, args) PyObject *self; PyObject *args; { double m1,m2,M; PyArrayObject *zs; PyArrayObject *EMass; npy_intp ld[1]; int i; double *z; /* parse arguments */ if (!PyArg_ParseTuple(args, "dddO", &m1,&m2,&M,&zs)) return NULL; /* create output array */ ld[0]= Cp->nelts+2; EMass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* allocate memory for the metallicity array */ z = malloc((Cp->nelts) * sizeof(double)); /* export values */ for (i=0;inelts;i++) z[i]= *(double *)(zs->data + (i)*(zs->strides[0])); /* compute SN ejection */ Total_mass_ejection(m1,m2,M,z); /* import values */ for (i=0;inelts+2;i++) *(double *)(EMass->data + (i)*(EMass->strides[0])) = EjectedMass[i]; free(z); /* convert in array */ return PyArray_Return(EMass); } static PyObject * chemistry_Total_mass_ejection_Z(self, args) PyObject *self; PyObject *args; { double m1,m2,M; PyArrayObject *zs; PyArrayObject *EMass; npy_intp ld[1]; int i; double *z; /* parse arguments */ if (!PyArg_ParseTuple(args, "dddO", &m1,&m2,&M,&zs)) return NULL; /* create output array */ ld[0]= Cp->nelts+2; EMass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* allocate memory for the metallicity array */ z = malloc((Cp->nelts) * sizeof(double)); /* export values */ for (i=0;inelts;i++) z[i]= *(double *)(zs->data + (i)*(zs->strides[0])); /* compute SN ejection */ Total_mass_ejection_Z(m1,m2,M,z); /* import values */ for (i=0;inelts+2;i++) *(double *)(EMass->data + (i)*(EMass->strides[0])) = EjectedMass[i]; free(z); /* convert in array */ return PyArray_Return(EMass); } static PyObject * chemistry_DYIN_Total_single_mass_ejection(self, args) PyObject *self; PyObject *args; { double m1; PyArrayObject *zs; PyArrayObject *EMass; npy_intp ld[1]; int i; double *z; /* parse arguments */ if (!PyArg_ParseTuple(args, "dO", &m1,&zs)) return NULL; /* create output array */ ld[0]= Cp->nelts+2; EMass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* allocate memory for the metallicity array */ z = malloc((Cp->nelts) * sizeof(double)); /* export values */ for (i=0;inelts;i++) z[i]= *(double *)(zs->data + (i)*(zs->strides[0])); /* compute dying stars ejection */ DYIN_Total_single_mass_ejection(m1,z); /* import values */ for (i=0;inelts+2;i++) *(double *)(EMass->data + (i)*(EMass->strides[0])) = SingleEjectedMass[i]; free(z); /* convert in array */ return PyArray_Return(EMass); } static PyObject * chemistry_SNII_Total_single_mass_ejection(self, args) PyObject *self; PyObject *args; { double m1; PyArrayObject *zs; PyArrayObject *EMass; npy_intp ld[1]; int i; double *z; /* parse arguments */ if (!PyArg_ParseTuple(args, "dO", &m1,&zs)) return NULL; /* create output array */ ld[0]= Cp->nelts+2; EMass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* allocate memory for the metallicity array */ z = malloc((Cp->nelts) * sizeof(double)); /* export values */ for (i=0;inelts;i++) z[i]= *(double *)(zs->data + (i)*(zs->strides[0])); /* compute SN ejection */ SNII_Total_single_mass_ejection(m1,z); /* import values */ for (i=0;inelts+2;i++) *(double *)(EMass->data + (i)*(EMass->strides[0])) = SingleEjectedMass[i]; free(z); /* convert in array */ return PyArray_Return(EMass); } static PyObject * chemistry_SNIa_Total_single_mass_ejection(self, args) PyObject *self; PyObject *args; { double m1; PyArrayObject *zs; PyArrayObject *EMass; npy_intp ld[1]; int i; double *z; /* parse arguments */ if (!PyArg_ParseTuple(args, "dO", &m1,&zs)) return NULL; /* create output array */ ld[0]= Cp->nelts+2; EMass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* allocate memory for the metallicity array */ z = malloc((Cp->nelts) * sizeof(double)); /* export values */ for (i=0;inelts;i++) z[i]= *(double *)(zs->data + (i)*(zs->strides[0])); /* compute SN ejection */ SNIa_Total_single_mass_ejection(m1,z); /* import values */ for (i=0;inelts+2;i++) *(double *)(EMass->data + (i)*(EMass->strides[0])) = SingleEjectedMass[i]; free(z); /* convert in array */ return PyArray_Return(EMass); } static PyObject * chemistry_Total_single_mass_ejection(self, args) PyObject *self; PyObject *args; { double m1; double NSNII,NSNIa,NDYIN; PyArrayObject *zs; PyArrayObject *EMass; npy_intp ld[1]; int i; double *z; /* parse arguments */ if (!PyArg_ParseTuple(args, "dOddd", &m1,&zs,&NSNII,&NSNIa,&NDYIN)) return NULL; /* create output array */ ld[0]= Cp->nelts+2; EMass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* allocate memory for the metallicity array */ z = malloc((Cp->nelts) * sizeof(double)); /* export values */ for (i=0;inelts;i++) z[i]= *(double *)(zs->data + (i)*(zs->strides[0])); /* compute SN ejection */ Total_single_mass_ejection(m1,z,NSNII,NSNIa,NDYIN); /* import values */ for (i=0;inelts+2;i++) *(double *)(EMass->data + (i)*(EMass->strides[0])) = SingleEjectedMass[i]; free(z); /* convert in array */ return PyArray_Return(EMass); } static PyObject * chemistry_Total_single_mass_ejection_Z(self, args) PyObject *self; PyObject *args; { double m1; double NSNII,NSNIa,NDYIN; PyArrayObject *zs; PyArrayObject *EMass; npy_intp ld[1]; int i; double *z; /* parse arguments */ if (!PyArg_ParseTuple(args, "dOddd", &m1,&zs,&NSNII,&NSNIa,&NDYIN)) return NULL; /* create output array */ ld[0]= Cp->nelts+2; EMass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* allocate memory for the metallicity array */ z = malloc((Cp->nelts) * sizeof(double)); /* export values */ for (i=0;inelts;i++) z[i]= *(double *)(zs->data + (i)*(zs->strides[0])); /* compute SN ejection */ Total_single_mass_ejection_Z(m1,z,NSNII,NSNIa,NDYIN); /* import values */ for (i=0;inelts+2;i++) *(double *)(EMass->data + (i)*(EMass->strides[0])) = SingleEjectedMass[i]; free(z); /* convert in array */ return PyArray_Return(EMass); } /*********************************/ /* */ /*********************************/ static PyObject * chemistry_cooling_function(self, args) PyObject *self; PyObject *args; { /* on gives : u_energy metal = metal(i,2) parameters t_const,zmin,zmax,slz,tmin,tmax,slt,FeHSolar,cooling_data_max */ PyArrayObject *cooling_data; double u_energy,metal; double t_const,zmin,zmax,slz,tmin,tmax,slt,FeHSolar,cooling_data_max; double cooling,u_cutoff,T,Z; double rt, rz, ft, fz, v1, v2, v; int it,iz,itp,izp; /* parse arguments */ if (!PyArg_ParseTuple(args, "ddOddddddddd", &u_energy, &metal, &cooling_data,&t_const,&zmin,&zmax,&slz,&tmin,&tmax,&slt,&FeHSolar,&cooling_data_max)) return NULL; u_cutoff=(100)/t_const; cooling = 0.0; if (u_energy > u_cutoff) { T = log10( t_const*u_energy ); Z = log10( metal/FeHSolar + 1.e-10 ); if (Z>zmax) { /*print *,'Warning: Z>Zmax for',i*/ Z=zmax; } if (Z < zmin) { rt = (T-tmin)/slt; it = (int)rt; if (it < cooling_data_max ) it = (int)rt; else it = cooling_data_max; itp = it+1; ft = rt - it; fz = ( 10. + Z )/( 10. + zmin); //v1 = ft*(cooling_data( 1, itp)-cooling_data( 1,it) ) + cooling_data( 1,it ); v1 = ft * (*(double *) (cooling_data->data + 1*(cooling_data->strides[0]) + itp*cooling_data->strides[1]) - *(double *) (cooling_data->data + 1*(cooling_data->strides[0]) + it *cooling_data->strides[1])) + *(double *) (cooling_data->data + 1*(cooling_data->strides[0]) + it *cooling_data->strides[1]); //v2 = ft*(cooling_data( 0,itp )-cooling_data( 0, it ) ) + cooling_data( 0, it ); v2 = ft * (*(double *) (cooling_data->data + 0*(cooling_data->strides[0]) + itp*cooling_data->strides[1]) - *(double *) (cooling_data->data + 0*(cooling_data->strides[0]) + it *cooling_data->strides[1])) + *(double *) (cooling_data->data + 0*(cooling_data->strides[0]) + it *cooling_data->strides[1]); v = v2 + fz*(v1-v2); } else { rt = (T-tmin)/slt; rz = (Z-zmin)/slz+1.0; if (it < cooling_data_max ) it = (int)rt; else it = cooling_data_max; iz = (int)rz; itp = it+1; izp = iz+1; ft = rt - it; fz = rz - iz; //v1 = ft*(cooling_data( izp, itp)-cooling_data(izp,it)) + cooling_data( izp, it ); v1 = ft * (*(double *) (cooling_data->data + izp*(cooling_data->strides[0]) + itp*cooling_data->strides[1]) - *(double *) (cooling_data->data + izp*(cooling_data->strides[0]) + it *cooling_data->strides[1])) + *(double *) (cooling_data->data + izp*(cooling_data->strides[0]) + it *cooling_data->strides[1]); //v2 = ft*(cooling_data( iz, itp )-cooling_data(iz,it )) + cooling_data( iz, it ); v2 = ft * (*(double *) (cooling_data->data + iz *(cooling_data->strides[0]) + itp*cooling_data->strides[1]) - *(double *) (cooling_data->data + iz *(cooling_data->strides[0]) + it *cooling_data->strides[1])) + *(double *) (cooling_data->data + iz *(cooling_data->strides[0]) + it *cooling_data->strides[1]); v = v2 + fz*(v1-v2); } cooling = pow(10,v); } return Py_BuildValue("d",cooling); } /*********************************/ /* */ /*********************************/ static PyObject * chemistry_get_Mmax(self, args) PyObject *self; PyObject *args; { return Py_BuildValue("d",(double)Cp->Mmax * All.MsoltoCMU); } static PyObject * chemistry_get_Mmin(self, args) PyObject *self; PyObject *args; { return Py_BuildValue("d",(double)Cp->Mmin * All.MsoltoCMU); } static PyObject * chemistry_get_Mco(self, args) PyObject *self; PyObject *args; { return Py_BuildValue("d",(double)Cp->Mco * All.MsoltoCMU); } static PyObject * chemistry_get_SNIa_Mpl(self, args) PyObject *self; PyObject *args; { return Py_BuildValue("d",(double)Cp->SNIa_Mpl * All.MsoltoCMU); } static PyObject * chemistry_get_SNIa_Mpu(self, args) PyObject *self; PyObject *args; { return Py_BuildValue("d",(double)Cp->SNIa_Mpu * All.MsoltoCMU); } static PyObject * chemistry_get_SNII_Mmin(self, args) PyObject *self; PyObject *args; { return Py_BuildValue("d",(double)Cp->SNII_Mmin * All.MsoltoCMU); } static PyObject * chemistry_get_SNII_Mmax(self, args) PyObject *self; PyObject *args; { return Py_BuildValue("d",(double)Cp->SNII_Mmax * All.MsoltoCMU); } static PyObject * chemistry_get_DYIN_Mmin(self, args) PyObject *self; PyObject *args; { return Py_BuildValue("d",(double)Cp->DYIN_Mmin * All.MsoltoCMU); } static PyObject * chemistry_get_DYIN_Mmax(self, args) PyObject *self; PyObject *args; { return Py_BuildValue("d",(double)Cp->DYIN_Mmax * All.MsoltoCMU); } static PyObject * chemistry_get_imf_Ntot(self, args) PyObject *self; PyObject *args; { return Py_BuildValue("d",(double)Cp->imf_Ntot*All.CMUtoMsol); /* in code mass unit */ } static PyObject * chemistry_get_as(self, args) PyObject *self; PyObject *args; { PyArrayObject *as; npy_intp ld[1]; int i; /* create output array */ ld[0]= Cp->n+1; as = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* import values */ for (i=0;in+1;i++) *(double *)(as->data + (i)*(as->strides[0])) = Cp->as[i]; return PyArray_Return(as); } static PyObject * chemistry_get_bs(self, args) PyObject *self; PyObject *args; { PyArrayObject *bs; npy_intp ld[1]; int i; /* create output array */ ld[0]= Cp->n+1; bs = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* import values */ for (i=0;in+1;i++) *(double *)(bs->data + (i)*(bs->strides[0])) = Cp->bs[i]; return PyArray_Return(bs); } static PyObject * chemistry_get_fs(self, args) PyObject *self; PyObject *args; { PyArrayObject *fs; npy_intp ld[1]; int i; /* create output array */ ld[0]= Cp->n; fs = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* import values */ for (i=0;in;i++) *(double *)(fs->data + (i)*(fs->strides[0])) = Cp->fs[i]; return PyArray_Return(fs); } static PyObject * chemistry_get_allnelts(self, args) PyObject *self; PyObject *args; { return Py_BuildValue("i",(int)Cp->nelts+2); } static PyObject * chemistry_get_nelts(self, args) PyObject *self; PyObject *args; { return Py_BuildValue("i",(int)Cp->nelts); } static PyObject * chemistry_get_allelts_labels(self, args) PyObject *self; PyObject *args; { int i; PyObject *LabelList,*LabelString; LabelList = PyList_New((Py_ssize_t)Cp->nelts+2); for(i=0;inelts+2;i++) { LabelString = PyString_FromString(Elt[i].label); PyList_SetItem(LabelList, (Py_ssize_t)i,LabelString); } return Py_BuildValue("O",LabelList); } static PyObject * chemistry_get_elts_labels(self, args) PyObject *self; PyObject *args; { int i; PyObject *LabelList,*LabelString; LabelList = PyList_New((Py_ssize_t)Cp->nelts); for(i=2;inelts+2;i++) { LabelString = PyString_FromString(Elt[i].label); PyList_SetItem(LabelList, (Py_ssize_t)i-2,LabelString); } return Py_BuildValue("O",LabelList); } /* static PyObject * chemistry_get_elts_SolarMassAbundances(self, args) PyObject *self; PyObject *args; { int i; npy_intp ld[1]; PyArrayObject *AbList; ld[0] = Cp->nelts; AbList = (PyArrayObject *) PyArray_SimpleNew(1,ld,NPY_FLOAT); for(i=2;inelts+2;i++) *(float*)(AbList->data + (i-2)*(AbList->strides[0])) = (float) Elt[i].SolarMassAbundance; return PyArray_Return(AbList); } */ static PyObject * chemistry_get_elts_SolarMassAbundances(self, args) PyObject *self; PyObject *args; { int i; PyObject *AbDict,*LabelString,*AbVal; AbDict = PyDict_New(); for(i=2;inelts+2;i++) { AbVal = PyFloat_FromDouble(Elt[i].SolarMassAbundance); LabelString = PyString_FromString(Elt[i].label); PyDict_SetItem(AbDict,LabelString, AbVal); } return Py_BuildValue("O",AbDict); } static PyObject * chemistry_get_MSNIa(self, args) PyObject *self; PyObject *args; { PyArrayObject *MSNIa; npy_intp ld[1]; int i; /* create output array */ ld[0]= Cp->nelts; MSNIa = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* import values */ for (i=0;inelts;i++) *(double *)(MSNIa->data + (i)*(MSNIa->strides[0])) = Elt[i+2].MSNIa*All.MsoltoCMU; return PyArray_Return(MSNIa); } static PyObject * chemistry_get_MassFracSNII(self, args) PyObject *self; PyObject *args; { PyArrayObject *MassFrac; npy_intp ld[1]; int i; /* create output array */ ld[0]= Cp->nelts+2; MassFrac = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* import values */ for (i=0;inelts+2;i++) *(double *)(MassFrac->data + (i)*(MassFrac->strides[0])) = MassFracSNII[i]; return PyArray_Return(MassFrac); } static PyObject * chemistry_get_SingleMassFracSNII(self, args) PyObject *self; PyObject *args; { PyArrayObject *MassFrac; npy_intp ld[1]; int i; /* create output array */ ld[0]= Cp->nelts+2; MassFrac = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); /* import values */ for (i=0;inelts+2;i++) *(double *)(MassFrac->data + (i)*(MassFrac->strides[0])) = SingleMassFracSNII[i]; return PyArray_Return(MassFrac); } static PyObject * chemistry_imf_sampling(self, args) PyObject *self; PyObject *args; { PyArrayObject *ms; npy_intp ld[1]; int i; int n,seed; /* parse arguments */ if (!PyArg_ParseTuple(args, "ii", &n,&seed)) return NULL; /* create output array */ ld[0]= n; ms = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); srandom(seed); /* import values */ for (i=0;idata + (i)*(ms->strides[0])) = imf_sampling(); return PyArray_Return(ms); } static PyObject * chemistry_imf_init_seed(self, args) PyObject *self; PyObject *args; { int n,seed; /* parse arguments */ if (!PyArg_ParseTuple(args, "i",&seed)) return NULL; srandom(seed); return Py_BuildValue("i",1); } static PyObject * chemistry_imf_sampling_single(self, args) PyObject *self; PyObject *args; { return Py_BuildValue("f",imf_sampling()); } static PyObject * chemistry_imf_sampling_single_from_random(self, args) PyObject *self; PyObject *args; { double f; if (!PyArg_ParseTuple(args, "f",&f)) return NULL; return Py_BuildValue("f",imf_sampling_from_random(f)); } /*********************************/ /* */ /*********************************/ static PyObject * chemistry_SNII_rate_P(self, args) PyObject *self; PyObject *args; { PyArrayObject *ConstSN,*Msn; double m1,m2,md; double powSN1,powSN2; double RSNII; RSNII = 0.0; /* parse arguments */ if (!PyArg_ParseTuple(args, "ddddOO", &m1,&m2,&powSN1,&powSN2,&ConstSN,&Msn)) return NULL; if ( m1 < *(double *) (Msn->data + 2*(Msn->strides[0]) + 1*(Msn->strides[1])) ) md = *(double *) (Msn->data + 2*(Msn->strides[0]) + 1*(Msn->strides[1])); else md = m1; if (md >= m2) RSNII = 0; else RSNII = *(double *) (ConstSN->data + 2*ConstSN->strides[0]) *(pow(m2,powSN1)-pow(md,powSN1)); return Py_BuildValue("d",RSNII); } /*********************************/ /* */ /*********************************/ static PyObject * chemistry_SNIa_rate_P(self, args) PyObject *self; PyObject *args; { PyArrayObject *ConstSN,*Msn; double m1,m2,md,mu; double powSN1,powSN2; double RSNIa; double rate; int i; /* parse arguments */ if (!PyArg_ParseTuple(args, "ddddOO", &m1,&m2,&powSN1,&powSN2,&ConstSN,&Msn)) return NULL; RSNIa = 0.0; for (i=0;i<2;i++) { if ( m1 < *(double *) (Msn->data + i*(Msn->strides[0]) + 0*(Msn->strides[1])) ) md = *(double *) (Msn->data + i*(Msn->strides[0]) + 0*(Msn->strides[1])); else md = m1; if ( m2 > *(double *) (Msn->data + i*(Msn->strides[0]) + 1*(Msn->strides[1])) ) mu = *(double *) (Msn->data + i*(Msn->strides[0]) + 1*(Msn->strides[1])); else mu = m2; if (mddata + i*ConstSN->strides[0])*(pow(mu,powSN2)-pow(md,powSN2)); } if ( m1 < *(double *) (Msn->data + 2*(Msn->strides[0]) + 0*(Msn->strides[1])) ) md = *(double *) (Msn->data + 2*(Msn->strides[0]) + 0*(Msn->strides[1])); else md = m1; mu = *(double *) (Msn->data + 2*(Msn->strides[0]) + 1*(Msn->strides[1])); if (md >= mu) RSNIa = 0.0; else RSNIa = RSNIa*(pow(mu,powSN1)-pow(md,powSN1)); if (RSNIa<0) RSNIa = 0; return Py_BuildValue("d",RSNIa); } /*********************************/ /* */ /*********************************/ static PyObject * chemistry_SNII_mass_ejection_P(self, args) PyObject *self; PyObject *args; { PyArrayObject *ArrayOrigin,*ArrayStep,*ChemArray; double m1,m2; int NbElement; double l1,l2; int i1,i2,i1p,i2p,j; double f1,f2; double v1,v2; PyArrayObject *ArrMassSNII; npy_intp ld[1]; /* parse arguments */ if (!PyArg_ParseTuple(args, "ddOOOi", &m1,&m2,&ArrayOrigin,&ArrayStep,&ChemArray,&NbElement)) return NULL; /* create output array */ ld[0]= NbElement+2; ArrMassSNII = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE); l1 = ( log10(m1) - *(double *)(ArrayOrigin->data + 0*(ArrayOrigin->strides[0])) ) / *(double *)(ArrayStep->data + 0*(ArrayStep->strides[0])) ; l2 = ( log10(m2) - *(double *)(ArrayOrigin->data + 0*(ArrayOrigin->strides[0])) ) / *(double *)(ArrayStep->data + 0*(ArrayStep->strides[0])) ; if (l1 < 0.0) l1 = 0.0; if (l2 < 0.0) l2 = 0.0; i1 = (int)l1; i2 = (int)l2; i1p = i1 + 1; i2p = i2 + 1; f1 = l1 - i1; f2 = l2 - i2; /* check (yr) */ if (i1<1) i1=1; if (i2<1) i2=1; /* --------- TOTAL GAS ---------- */ j = NbElement; v1=f1* (*(double *)(ChemArray->data + (i1p)*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])) - *(double *)(ChemArray->data + (i1 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]))) + *(double *)(ChemArray->data + (i1 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])); v2=f2* (*(double *)(ChemArray->data + (i2p)*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])) - *(double *)(ChemArray->data + (i2 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]))) + *(double *)(ChemArray->data + (i2 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])); *(double *)(ArrMassSNII->data + (j)*(ArrMassSNII->strides[0])) = v2-v1; /* --------- He core therm ---------- */ j = NbElement+1; v1=f1* (*(double *)(ChemArray->data + (i1p)*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])) - *(double *)(ChemArray->data + (i1 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]))) + *(double *)(ChemArray->data + (i1 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])); v2=f2* (*(double *)(ChemArray->data + (i2p)*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])) - *(double *)(ChemArray->data + (i2 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]))) + *(double *)(ChemArray->data + (i2 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])); *(double *)(ArrMassSNII->data + (j)*(ArrMassSNII->strides[0])) = v2-v1; /* --------- Metals ---------- */ l1 = ( log10(m1) - *(double *)(ArrayOrigin->data + 1*(ArrayOrigin->strides[0])) ) / *(double *)(ArrayStep->data + 1*(ArrayStep->strides[0])) ; l2 = ( log10(m2) - *(double *)(ArrayOrigin->data + 1*(ArrayOrigin->strides[0])) ) / *(double *)(ArrayStep->data + 1*(ArrayStep->strides[0])) ; if (l1 < 0.0) l1 = 0.0; if (l2 < 0.0) l2 = 0.0; i1 = (int)l1; i2 = (int)l2; i1p = i1 + 1; i2p = i2 + 1; f1 = l1 - i1; f2 = l2 - i2; /* check (yr) */ if (i1<1) i1=1; if (i2<1) i2=1; for (j=0;jdata + (i1p)*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])) - *(double *)(ChemArray->data + (i1 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]))) + *(double *)(ChemArray->data + (i1 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])); v2=f2* (*(double *)(ChemArray->data + (i2p)*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])) - *(double *)(ChemArray->data + (i2 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1]))) + *(double *)(ChemArray->data + (i2 )*(ChemArray->strides[0]) + (j)*(ChemArray->strides[1])); *(double *)(ArrMassSNII->data + (j)*(ArrMassSNII->strides[0])) = v2-v1; } return PyArray_Return(ArrMassSNII); } #ifdef CHIMIE_OPTIMAL_SAMPLING static PyObject * chemistry_optimal_init_norm(self, args) PyObject *self; PyObject *args; { double m_max,Mecl; if (!PyArg_ParseTuple(args, "d", &Mecl)) return NULL; m_max = optimal_init_norm(Mecl); return Py_BuildValue("d",m_max); } static PyObject * chemistry_optimal_get_next_mass(self, args) PyObject *self; PyObject *args; { double m; double m_next; /* parse arguments */ if (!PyArg_ParseTuple(args, "d", &m)) return NULL; m_next = optimal_get_next_mass(m); return Py_BuildValue("d",m_next); } static PyObject * chemistry_stop_loop(self, args) PyObject *self; PyObject *args; { double m; int bool; /* parse arguments */ if (!PyArg_ParseTuple(args, "d", &m)) return NULL; bool = optimal_stop_loop(m); return Py_BuildValue("i",bool); } static PyObject * chemistry_optimal_get_m1_from_m2(self, args) PyObject *self; PyObject *args; { double m1,m2,mp; /* parse arguments */ if (!PyArg_ParseTuple(args, "dd", &m2,&mp)) return NULL; m1 = optimal_get_m1_from_m2(m2,mp); return Py_BuildValue("d",m1); } #endif /* definition of the method table */ static PyMethodDef chemistryMethods[] = { {"CodeUnits_to_SolarMass_Factor", chemistry_CodeUnits_to_SolarMass_Factor, METH_VARARGS, "convertion factor : CodeUnits -> SolarMass"}, {"SolarMass_to_CodeUnits_Factor", chemistry_SolarMass_to_CodeUnits_Factor, METH_VARARGS, "convertion factor : SolarMass -> CodeUnits"}, {"SetVerbosityOn", (PyCFunction)chemistry_SetVerbosityOn, METH_VARARGS, "Set verbosity to on"}, {"SetVerbosityOff", (PyCFunction)chemistry_SetVerbosityOff, METH_VARARGS, "Set verbosity to off"}, {"InitDefaultParameters", (PyCFunction)chemistry_InitDefaultParameters, METH_VARARGS, "Init default parameters"}, {"SetParameters", (PyCFunction)chemistry_SetParameters, METH_VARARGS, "Set gadget parameters"}, {"GetParameters", (PyCFunction)chemistry_GetParameters, METH_VARARGS, "get some gadget parameters"}, {"init_chimie", chemistry_init_chimie, METH_VARARGS| METH_KEYWORDS, "Init chimie."}, {"info", chemistry_info, METH_VARARGS| METH_KEYWORDS, "Get info on tables."}, {"set_table", chemistry_set_table, METH_VARARGS, "Set the chimie table."}, {"get_imf", chemistry_get_imf, METH_VARARGS, "Compute corresponding imf value."}, {"get_imf_M", chemistry_get_imf_M, METH_VARARGS, "Compute the mass fraction between m1 and m2."}, {"get_imf_N", chemistry_get_imf_N, METH_VARARGS, "Compute the fraction number between m1 and m2."}, {"star_lifetime", chemistry_star_lifetime, METH_VARARGS, "Compute star life time."}, {"star_mass_from_age", chemistry_star_mass_from_age, METH_VARARGS, "Return the stellar mass that has a lifetime equal to t."}, {"DYIN_rate", chemistry_DYIN_rate, METH_VARARGS, "Return the number of dying stars per unit mass with masses between m1 and m2."}, {"SNII_rate", chemistry_SNII_rate, METH_VARARGS, "Return the number of SNII per unit mass with masses between m1 and m2."}, {"SNIa_rate", chemistry_SNIa_rate, METH_VARARGS, "Return the number of SNIa per unit mass with masses between m1 and m2."}, {"SNIa_single_rate", chemistry_SNIa_single_rate, METH_VARARGS, "Return the number of SNIa per unit mass for a star of mass m."}, {"DYIN_mass_ejection", chemistry_DYIN_mass_ejection, METH_VARARGS, "Mass fraction of ejected elements, per unit mass due to the explotion of dying stars with masses between m1 and m2."}, {"DYIN_mass_ejection_Z", chemistry_DYIN_mass_ejection_Z, METH_VARARGS, "Mass fraction of ejected elements, per unit mass due to the explotion of dying stars with masses between m1 and m2 (metallicity dependency)."}, {"DYIN_single_mass_ejection", chemistry_DYIN_single_mass_ejection, METH_VARARGS, "Mass fraction of ejected elements due to the explotion of one dying star of mass m."}, {"DYIN_single_mass_ejection_Z", chemistry_DYIN_single_mass_ejection_Z, METH_VARARGS, "Mass fraction of ejected elements due to the explotion of one dying star of mass m (metallicity dependency)."}, {"SNII_mass_ejection", chemistry_SNII_mass_ejection, METH_VARARGS, "Mass fraction of ejected elements, per unit mass due to the explotion of SNII with masses between m1 and m2."}, //{"SNII_mass_ejection_Z", chemistry_SNII_mass_ejection_Z, METH_VARARGS, // "Mass fraction of ejected elements, per unit mass due to the explotion of SNII with masses between m1 and m2 (metallicity dependency)."}, {"SNII_single_mass_ejection", chemistry_SNII_single_mass_ejection, METH_VARARGS, "Mass fraction of ejected elements due to the explotion of one SNII of mass m."}, //{"SNII_single_mass_ejection_Z", chemistry_SNII_single_mass_ejection_Z, METH_VARARGS, // "Mass fraction of ejected elements due to the explotion of one SNII star of mass m (metallicity dependency)."}, {"SNIa_mass_ejection", chemistry_SNIa_mass_ejection, METH_VARARGS, "Mass fraction of ejected elements, per unit mass due to the explotion of SNIa with masses between m1 and m2."}, {"SNIa_single_mass_ejection", chemistry_SNIa_single_mass_ejection, METH_VARARGS, "Mass fraction of ejected elements due to the explotion of one SNIa of mass m."}, {"Total_mass_ejection", chemistry_Total_mass_ejection, METH_VARARGS, "Mass fraction of ejected elements, per unit mass due to the explotion of SNIa and SNII with masses between m1 and m2."}, {"Total_mass_ejection_Z", chemistry_Total_mass_ejection_Z, METH_VARARGS, "Mass fraction of ejected elements, per unit mass due to the explotion of SNIa and SNII with masses between m1 and m2 (metallicity dependency)."}, {"DYIN_Total_single_mass_ejection", chemistry_DYIN_Total_single_mass_ejection, METH_VARARGS, "Mass fraction of ejected elements (including processed and non processed elements) due to the explotion of one dying star of mass m."}, {"SNII_Total_single_mass_ejection", chemistry_SNII_Total_single_mass_ejection, METH_VARARGS, "Mass fraction of ejected elements (including processed and non processed elements) due to the explotion of one SNII of mass m."}, {"SNIa_Total_single_mass_ejection", chemistry_SNIa_Total_single_mass_ejection, METH_VARARGS, "Mass fraction of ejected elements (including processed and non processed elements) due to the explotion of one SNIa of mass m."}, {"Total_single_mass_ejection", chemistry_Total_single_mass_ejection, METH_VARARGS, "Mass fraction of ejected elements, per unit mass due to the explotion of one dying star of mass m1."}, {"Total_single_mass_ejection_Z", chemistry_Total_single_mass_ejection_Z, METH_VARARGS, "Mass fraction of ejected elements, per unit mass due to the explotion of one dying star of mass m1 (metallicity dependency)."}, {"get_Mmax", chemistry_get_Mmax, METH_VARARGS, "Get max star mass of the IMF, in code unit."}, {"get_Mmin", chemistry_get_Mmin, METH_VARARGS, "Get min star mass of the IMF, in code unit."}, {"get_Mco", chemistry_get_Mco, METH_VARARGS, "Get mean WD mass, in code unit."}, {"get_SNIa_Mpl", chemistry_get_SNIa_Mpl, METH_VARARGS, "Get min mass of SNIa, in code unit."}, {"get_SNIa_Mpu", chemistry_get_SNIa_Mpu, METH_VARARGS, "Get max mass of SNIa, in code unit."}, {"get_SNII_Mmin", chemistry_get_SNII_Mmin, METH_VARARGS, "Get min mass of SNII, in code unit."}, {"get_SNII_Mmax", chemistry_get_SNII_Mmax, METH_VARARGS, "Get max mass of SNII, in code unit."}, {"get_DYIN_Mmin", chemistry_get_DYIN_Mmin, METH_VARARGS, "Get min mass of DYIN, in code unit."}, {"get_DYIN_Mmax", chemistry_get_DYIN_Mmax, METH_VARARGS, "Get max mass of DYIN, in code unit."}, {"get_imf_Ntot", chemistry_get_imf_Ntot, METH_VARARGS, "Get number of stars in the imf, per unit mass."}, {"get_as", chemistry_get_as, METH_VARARGS, "Get power coefficients."}, {"get_bs", chemistry_get_bs, METH_VARARGS, "Get normalisation coefficients."}, {"get_fs", chemistry_get_fs, METH_VARARGS, "Get fs, mass fraction at ms."}, {"get_allnelts", chemistry_get_allnelts, METH_VARARGS, "Get the number of element considered, including ejected mass (Ej) and non processed ejected mass (Ejnp).."}, {"get_nelts", chemistry_get_nelts, METH_VARARGS, "Get the number of element considered."}, {"get_allelts_labels", chemistry_get_allelts_labels, METH_VARARGS, "Get the labels of elements, including ejected mass (Ej) and non processed ejected mass (Ejnp)."}, {"get_elts_labels", chemistry_get_elts_labels, METH_VARARGS, "Get the labels of elements."}, {"get_elts_SolarMassAbundances", chemistry_get_elts_SolarMassAbundances, METH_VARARGS, "Get the solar mass abundance of elements."}, {"get_MassFracSNII", chemistry_get_MassFracSNII, METH_VARARGS, "Get the mass fraction per element ejected by a set of SNII."}, {"get_SingleMassFracSNII", chemistry_get_SingleMassFracSNII, METH_VARARGS, "Get the mass fraction per element ejected by a SNII."}, {"get_MSNIa", chemistry_get_MSNIa, METH_VARARGS, "Get the mass per element ejected by a SNIa."}, {"cooling_function", chemistry_cooling_function, METH_VARARGS, "Compute cooling."}, {"imf_sampling", chemistry_imf_sampling, METH_VARARGS, "Sample imf with n points."}, {"imf_sampling_single", chemistry_imf_sampling_single, METH_VARARGS, "Sample imf with a single point."}, {"imf_init_seed", chemistry_imf_init_seed, METH_VARARGS, "Init the random seed."}, {"imf_sampling_single_from_random", chemistry_imf_sampling_single_from_random, METH_VARARGS, "Sample imf with a single point from a given random number."}, /* old poirier */ {"SNIa_rate_P", chemistry_SNIa_rate_P, METH_VARARGS, "Return the number of SNIa per unit mass and time. (Poirier version)"}, {"SNII_rate_P", chemistry_SNII_rate_P, METH_VARARGS, "Return the number of SNII per unit mass and time. (Poirier version)"}, {"SNII_mass_ejection_P", chemistry_SNII_mass_ejection_P, METH_VARARGS, "Mass ejection due to SNII per unit mass and time. (Poirier version)"}, #ifdef CHIMIE_OPTIMAL_SAMPLING {"optimal_sampling_init_norm", chemistry_optimal_init_norm, METH_VARARGS, "init normalistation"}, {"optimal_sampling_get_next_mass", chemistry_optimal_get_next_mass, METH_VARARGS, "next star mass"}, {"optimal_sampling_stop_loop", chemistry_stop_loop, METH_VARARGS, "return 1 if the loop for optimal sampling must be stopped "}, {"optimal_sampling_get_m1_from_m2", chemistry_optimal_get_m1_from_m2, METH_VARARGS, "for a givent mass m2, return m1, such that the mass in the IMF between m1 and m2 is mp."}, #endif {NULL, NULL, 0, NULL} /* Sentinel */ }; void initchemistry(void) { (void) Py_InitModule("chemistry", chemistryMethods); import_array(); } #endif /* PYCHEM */ #endif /* CHIMIE */