Page MenuHomec4science
Files F115522132

pychem_chemivol
No OneTemporary
Actions

File Metadata

Created
Sun, Jun 1, 13:12

pychem_chemivol
View Options

#!/usr/bin/env python
from optparse import OptionParser
from PyChem import chemistry
from numpy import *
import sys
import string
import Ptools as pt
from pNbody import units
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("--dt",
action="store",
dest="dt",
type="float",
default = 1.0,
help="dt",
metavar=" FLOAT")
parser.add_option("--tmax",
action="store",
dest="tmax",
type="float",
default = 14000.0,
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("--mgas",
action="store",
dest="mgas",
type="float",
default = 1e5,
help="initial mass of the gas reservoir",
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="outputfile",
type="string",
default = None,
help="output file",
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")
(options, args) = parser.parse_args()
pt.check_files_number(args)
files = args
return files,options
########################################################################
# some objects
########################################################################
class Star:
def __init__(self,tstar,mass,metals,name=None):
self.tstar = tstar
self.mass = mass
self.metals = metals
self.M0 = mass
self.name = name
def compute_ejectats(self,time,dt):
minlivetime = chemistry.star_lifetime(self.metals[-1],chemistry.get_Mmax())
maxlivetime = chemistry.star_lifetime(self.metals[-1],chemistry.get_Mmin())
ejectats = {}
ejectats["TotalEjectedEltMass"] = zeros(nelts)
ejectats["TotalEjectedGasMass"] = 0
ejectats["TotalInjectedEnergy"] = 0
ejectats["TotalSNII"] = 0
ejectats["TotalSNIa"] = 0
ejectats["TotalDYIN"] = 0
tpdt = time + dt
if (tpdt - self.tstar >= minlivetime) and (tpdt - self.tstar < maxlivetime):
m1 = chemistry.star_mass_from_age(self.metals[-1],tpdt - self.tstar)
if (time - self.tstar >= minlivetime):
m2 = chemistry.star_mass_from_age(self.metals[-1],time - self.tstar)
else:
m2 = chemistry.get_Mmax()
if (m1>m2):
m1=m2
raise "m1>m2 !!!"
# number of SNIa
NSNIa = chemistry.SNIa_rate(m1,m2)*self.M0
# number of SNII
NSNII = chemistry.SNII_rate(m1,m2)*self.M0
# number of dying stars
NDYIN = chemistry.DYIN_rate(m1,m2)*self.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
#
#########################
# standard way (continuous explosion) (_cont)
EjectedMass = chemistry.Total_mass_ejection_Z(m1,m2,self.M0,self.metals)
TotalEjectedEltMasscont = EjectedMass[2:]
TotalEjectedGasMasscont = EjectedMass[0]
# new way (discrete explosion) (_j)
EjectedMassj = chemistry.Total_single_mass_ejection_Z(0.5*(m1+m2),self.metals,float(NSNII),float(NSNIa),float(NDYIN))
TotalEjectedEltMassj = EjectedMassj[2:]
TotalEjectedGasMassj = EjectedMassj[0]
# energy injected
TotalInjectedEnergy = feedback_energy* (NSNIa + NSNII)
# remove mass
self.mass = self.mass - TotalEjectedGasMasscont
# outputs
ejectats["TotalEjectedEltMass"] = TotalEjectedEltMasscont
ejectats["TotalEjectedGasMass"] = TotalEjectedGasMasscont
ejectats["TotalInjectedEnergy"] = TotalInjectedEnergy
ejectats["TotalSNII"] = NSNII
ejectats["TotalSNIa"] = NSNIa
ejectats["TotalDYIN"] = NDYIN
return ejectats
########################################################################
# some functinos
########################################################################
def SFR(t):
"SFR in mass solar per year"
t = t*UnitsConvTime_gear2Myrs # in Myr
if t < 2000:
return 0.0005;
else:
return 0.;
########################################################################
# M A I N
########################################################################
##################################
# units
##################################
# gear system of units
params = {}
params['UnitVelocity_in_cm_per_s'] = 20725573.785998672
params['UnitMass_in_g'] = 1.989e+43
params['UnitLength_in_cm'] = 3.085678e+21
gear_units = units.Set_SystemUnits_From_Params(params)
# units sfr (Msol/yr)
sfr_units = units.UnitSystem('local',[units.Unit_cm,units.Unit_Msol,units.Unit_yr,units.Unit_K])
# units time (Myr)
time_units = units.UnitSystem('local',[units.Unit_cm,units.Unit_Msol,units.Unit_Myr,units.Unit_K])
# some conversion factors
UnitsConvEnergy_cgs2gear = 1.0/gear_units.convertionFactorTo(units.cgs.UnitEnergy)
UnitsConvSfr_sfr2gear = 1.0/( gear_units.convertionFactorTo(sfr_units.UnitMass) / gear_units.convertionFactorTo(sfr_units.UnitTime) )
UnitsConvMass_Msol2gear = 1.0/gear_units.convertionFactorTo(sfr_units.UnitMass)
UnitsConvMass_gear2Msol = 1/UnitsConvMass_Msol2gear
UnitsConvTime_Myrs2gear = 1.0/gear_units.convertionFactorTo(time_units.UnitTime)
UnitsConvTime_gear2Myrs = 1/UnitsConvTime_Myrs2gear
# other
feedback_energy = 10**51 * UnitsConvEnergy_cgs2gear
# parse options
files,opt = parse_options()
if opt.seed!=None:
random.seed(opt.seed)
if opt.xmax==None:
opt.xmax = 14000
# to plot
# init chimie
chemistry.init_chimie(files[0],opt.NumberOfTables,opt.DefaultTable)
elts = chemistry.get_elts_labels()
nelts = chemistry.get_nelts()
# some parameters
dt = opt.dt *UnitsConvTime_Myrs2gear
tmax = opt.tmax *UnitsConvTime_Myrs2gear
mstar = opt.mstar *UnitsConvMass_Msol2gear # mass of a SSP
Mgas = opt.mgas *UnitsConvMass_Msol2gear # mass of the gas reservoir
# stellar initial element abondance for the gas
SolarAbun = chemistry.get_elts_SolarMassAbundances()
Zratio = opt.Z/SolarAbun['Metals']
metals = zeros(chemistry.get_nelts(),float)
for i,elt in enumerate(elts):
metals[i] = Zratio*SolarAbun[elt]
for i,elt in enumerate(elts):
print "%6s X=%g (solar value =%g )"%(elt,metals[i],SolarAbun[elt])
# times
times = arange(dt,tmax,dt).astype(float)
Stars = []
n = len(times)
# outputs
OUT = {}
OUT["Times"] = zeros(n)
OUT["Sfr"] = zeros(n,float)
OUT["Sfr_th"] = zeros(n,float)
OUT["TotStarMass"] = zeros(n,float)
OUT["NStars"] = zeros(n,float)
OUT["TotGasMass"] = zeros(n,float)
OUT["TotMetals"] = zeros((n,nelts),float)
OUT["TotSNII"] = zeros(n,int)
OUT["TotSNIa"] = zeros(n,int)
OUT["TotDYIN"] = zeros(n,int)
##############################################################################################
#
# main loop
#
##############################################################################################
for ti,time in enumerate(times):
###################
# star formation
###################
# N is the number of IMF of mass mstar formed during the timestep
dMdt = SFR(time)*UnitsConvSfr_sfr2gear
dM_star = dMdt * dt
N = dM_star/mstar
# compute discrete values
fN = N-floor(N) # fraction
N = floor(N) # integer
if random.random() < fN:
N = N+1
Mnewstars = N*mstar
###################
# create the Stars
###################
for i in xrange(N):
Stars.append(Star(time,mstar,metals))
NStars = len(Stars)
#######################
# do chemical evolution
#######################
# record current amount of metals
Mx = metals * Mgas
TotalEjectedEltMass = zeros(nelts)
TotalEjectedGasMass = 0
TotalInjectedEnergy = 0
TotalSNII = 0
TotalSNIa = 0
TotalDYIN = 0
for star in Stars:
ejectats = star.compute_ejectats(time,dt)
totalEjectedEltMass = ejectats["TotalEjectedEltMass"]
totalEjectedGasMass = ejectats["TotalEjectedGasMass"]
totalInjectedEnergy = ejectats["TotalInjectedEnergy"]
totalSNII = ejectats["TotalSNII"]
totalSNIa = ejectats["TotalSNIa"]
totalDYIN = ejectats["TotalDYIN"]
TotalEjectedEltMass += totalEjectedEltMass
TotalEjectedGasMass += totalEjectedGasMass
TotalInjectedEnergy += totalInjectedEnergy
TotalSNII += totalSNII
TotalSNIa += totalSNIa
TotalDYIN += totalDYIN
Mgas = Mgas + TotalEjectedGasMass - Mnewstars
if Mgas<0:
print Mgas
raise "Mgas < 0"
# compute new gas metallicity
Mx = Mx + TotalEjectedEltMass
metals = Mx/Mgas
######################################
# some statistics
######################################
Mstars_tot = 0
for star in Stars:
Mstars_tot+=star.mass
OUT["Times"][ti] = time
OUT["TotStarMass"][ti] = Mstars_tot
OUT["TotGasMass"][ti] = Mgas
OUT["NStars"][ti] = NStars
OUT["Sfr"][ti] = Mnewstars/dt
OUT["Sfr_th"][ti] = dMdt
OUT["TotSNII"][ti] = TotalSNII
OUT["TotSNIa"][ti] = TotalSNIa
OUT["TotDYIN"][ti] = TotalDYIN
# some elements
OUT["TotMetals"][ti] = metals
######################################
# print some info
######################################
print "Step = %08d Time = %g NStar = %06d Mgas = %8.3g "%(ti, time*UnitsConvTime_gear2Myrs, NStars, Mgas*UnitsConvMass_gear2Msol )
if opt.outputfile!=None:
import pickle
f = open(opt.outputfile,'w')
pickle.dump(opt,f)
pickle.dump(gear_units,f)
pickle.dump(files[0],f)
pickle.dump(OUT,f)
f.close()
######################################
# crate a pNbody-object
######################################
from pNbody import *
pos = ones((NStars,3),float32)
nb = Nbody(p_name='snap.dat',ftype='gadget',pos=pos,status='new')
nb.set_tpe(1)
N=NStars
nb.tstar = zeros(N,float32)
nb.minit = zeros(N,float32)
nb.idp = zeros(N,int32)
nb.metals= zeros((N,nelts),float32)
nb.mass = zeros((N),float32)
for i in xrange(NStars):
nb.tstar[i] = Stars[i].tstar
nb.minit[i] = Stars[i].M0
nb.idp[i] = 0
nb.metals[i] = Stars[i].metals
nb.mass[i] = Stars[i].mass
nb.u = zeros(nb.nbody).astype(float32)
nb.rsp = ones(nb.nbody).astype(float32)
nb.rho = ones(nb.nbody).astype(float32)
# set the chimie extra-header
nb.flag_chimie_extraheader = 1
nb.ChimieNelements = int(chemistry.get_nelts())
nb.ChimieElements = chemistry.get_elts_labels()
nb.ChimieSolarMassAbundances = chemistry.get_elts_SolarMassAbundances()
NELEMENTS = nb.ChimieNelements
nb.flag_metals=NELEMENTS
nb.write()
##################################################################################################################
#
# plot
#
##################################################################################################################
sys.exit()
pt.InitPlot([],opt)
pt.pcolors
fig = pt.gcf()
#pt.plot(Times,TotStarMass)
#pt.plot(Times,TotGasMass)
pt.plot(Times,Sfr)
#pt.plot(Times,FeH)
pt.show()

Event Timeline

c4science · Help