Page MenuHomec4science
Files F84900346

gadget.py
No OneTemporary
Actions

File Metadata

Created
Wed, Sep 25, 11:47

gadget.py
View Options

####################################################################################################################################
#
# GADGET CLASS
#
####################################################################################################################################
class Nbody_gadget(NbodyDefault):
def __init__(self,p_name=None,pos=None,vel=None,mass=None,num=None,tpe=None,ftype=None,status='old',byteorder=sys.byteorder,pio='no',local=False,log=None):
NbodyDefault.__init__(self,p_name=p_name,pos=pos,vel=vel,mass=mass,num=num,tpe=tpe,ftype=ftype,status=status,byteorder=byteorder,pio=pio,local=local,log=log)
# initialize SSP
from pNbody.SSP import libvazdekis
# vazdekis_kb_mu1.3.txt : krupa 01 revisited
self.LObj = libvazdekis.VazdekisLuminosities(os.path.join(OPTDIR,'SSP','vazdekis_kb_mu1.3.txt'))
self.LObj.ExtrapolateMatrix(order=1,s=0)
self.LObj.CreateInterpolator()
self.LObj.Extrapolate2DMatrix()
def InitSpec(self):
self.getComovingIntegration()
def get_read_fcts(self):
return [self.read_particles]
def get_write_fcts(self):
return [self.write_particles]
def get_mxntpe(self):
return 6
def get_default_spec_vars(self):
'''
return specific variables default values for the class
'''
return {'massarr' :array([0,0,self.nbody,0,0,0]),
'atime' :0.,
'redshift' :0.,
'flag_sfr' :0,
'flag_feedback' :0,
'nall' :array([0,0,self.nbody,0,0,0]),
'flag_cooling' :0,
'num_files' :1,
'boxsize' :0.,
'omega0' :0.,
'omegalambda' :0.,
'hubbleparam' :0.,
'flag_age' :0.,
'hubbleparam' :0.,
'flag_metals' :0.,
'nallhw' :array([0,0,0,0,0,0]),
'flag_entr_ic' :0,
'flag_chimie_extraheader':0,
'critical_energy_spec' :0.,
'empty' :48*'',
'comovingintegration' :None
}
def get_massarr_and_nzero(self):
"""
return massarr and nzero
!!! when used in //, if only a proc has a star particle,
!!! nzero is set to 1 for all cpu, while massarr has a length of zero !!!
"""
if self.has_var('massarr') and self.has_var('nzero'):
if self.massarr !=None and self.nzero!=None:
print "warning : get_massarr_and_nzero : here we use massarr and nzero",self.massarr,self.nzero
return self.massarr,self.nzero
massarr = zeros(len(self.npart),float)
nzero = 0
# for each particle type, see if masses are equal
for i in range(len(self.npart)):
first_elt = sum((arange(len(self.npart)) < i) * self.npart)
last_elt = first_elt + self.npart[i]
if first_elt!=last_elt:
c = (self.mass[first_elt]==self.mass[first_elt:last_elt]).astype(int)
if sum(c)==len(c):
massarr[i] = self.mass[first_elt]
else:
nzero = nzero+len(c)
return massarr.tolist(),nzero
def getComovingIntegration(self):
"""
set true if the file has been runned using
the comoving integration scheme
"""
flag = True
# this is not very good, however, there is not other choice...
if self.omega0==0 and self.omegalambda==0:
flag = False
self.comovingintegration = flag
print "comovingintegration=",self.comovingintegration
def isComovingIntegrationOn(self):
"""
return true if the file has been runned using
the comoving integration scheme
"""
flag = True
# this is not very good, however, there is not other choice...
if self.omega0==0 and self.omegalambda==0:
flag = False
self.comovingintegration = flag
return flag
def setComovingIntegrationOn(self):
self.comovingintegration = True
def setComovingIntegrationOff(self):
self.comovingintegration = False
def read_particles(self,f):
'''
read gadget file
'''
##########################################
# read the header and send it to each proc
##########################################
tpl = (24,48,float,float,int32,int32,24,int32,int32,float,float,float,float,int32,int32,24,int32,int32,float,48)
header = io.ReadBlock(f,tpl,byteorder=self.byteorder,pio=self.pio)
npart,massarr,atime,redshift,flag_sfr,flag_feedback,nall,flag_cooling,num_files,boxsize,omega0,omegalambda,hubbleparam,flag_age,flag_metals,nallhw,flag_entr_ic,flag_chimie_extraheader,critical_energy_spec,empty = header
if fabs(flag_chimie_extraheader)>1: # if the header is empty, it may be indetermined
flag_chimie_extraheader=0
npart = fromstring(npart,int32)
massarr = fromstring(massarr,float)
nall = fromstring(nall,int32)
nallhw = fromstring(nallhw,int32)
if sys.byteorder != self.byteorder:
npart.byteswap(True)
massarr.byteswap(True)
nall.byteswap(True)
nallhw.byteswap(True)
if flag_metals: # gas metal properties
NELEMENTS = flag_metals
self.NELEMENTS = NELEMENTS
else:
NELEMENTS = 0
self.NELEMENTS = NELEMENTS
##########################################
# computes nzero
##########################################
# count number of particles that have non constant masses and then
# have masses storded further
if self.pio == 'no':
npart_tot = npart
npart_all = libutil.get_npart_all(npart,mpi.mpi_NTask())
npart = npart_all[mpi.mpi_ThisTask()] # local
npart_read = npart_tot
nbody_read = sum(npart_read)
npart_m_read = npart_tot * (massarr==0)
ngas = npart[0]
ngas_read = npart_tot[0]
nstars = npart[1]
nstars_read = npart_tot[1]
# compute nzero
nzero=0
mass = array([])
for i in range(len(npart_tot)):
if massarr[i] == 0:
nzero = nzero+npart_tot[i]
else:
mass = concatenate((mass,ones(npart[i])*massarr[i]))
else:
npart_tot = mpi.mpi_allreduce(npart)
npart_all = None # each proc read for himself
npart = npart # local
npart_read = None # each proc read for himself
nbody_read = sum(npart)
npart_m_read = None # each proc read for himself
ngas = npart[0]
ngas_read = ngas
nstars = npart[1]
nstars_read = nstars
# compute nzero
nzero=0
mass = array([])
for i in range(len(npart)):
if massarr[i] == 0:
nzero = nzero+npart[i]
else:
mass = concatenate((mass,ones(npart[i])*massarr[i]))
nbody = sum(npart)
nbody_tot = sum(npart_tot)
if npart_m_read !=None:
if sum(npart_m_read) != nzero:
raise "sum(npart_m) (%d) != nzero (%d)"%(sum(npart_m_read),nzero),npart_m_read
##########################################
# optionnally read extra header
##########################################
if flag_chimie_extraheader==1:
if mpi.mpi_IsMaster(): print "reading chimie extra-header..."
tpl = (int32, int(NELEMENTS)*4, 256-4-int(NELEMENTS)*4)
nelts,ChimieSolarAbundances,labels = io.ReadBlock(f,tpl,byteorder=self.byteorder,pio=self.pio)
nelts = int(nelts)
self.ChimieNelements = nelts
ChimieElements = string.split(labels,',')[:nelts]
self.ChimieElements = ChimieElements
ChimieSolarAbundances = fromstring(ChimieSolarAbundances,float32)
self.ChimieSolarAbundances = {}
for i,elt in enumerate(self.ChimieElements):
self.ChimieSolarAbundances[elt] = ChimieSolarAbundances[i]
else:
self.ChimieNelements = int(NELEMENTS)
self.ChimieElements = ['Fe','Mg','O','Metals']
self.ChimieSolarAbundances = {}
self.ChimieSolarAbundances['Fe'] = 0.001771
self.ChimieSolarAbundances['Mg'] = 0.00091245
self.ChimieSolarAbundances['O'] = 0.0108169
self.ChimieSolarAbundances['Metals'] = 0.02
##########################################
# read and send particles attribute
##########################################
if mpi.mpi_IsMaster(): print "reading pos..."
pos = io.ReadDataBlock(f,float32,shape=(nbody_read,3),byteorder=self.byteorder,pio=self.pio,npart=npart_read)
if mpi.mpi_IsMaster(): print "reading vel..."
vel = io.ReadDataBlock(f,float32,shape=(nbody_read,3),byteorder=self.byteorder,pio=self.pio,npart=npart_read)
if mpi.mpi_IsMaster(): print "reading num..."
num = io.ReadDataBlock(f,int32 ,shape=(nbody_read,) ,byteorder=self.byteorder,pio=self.pio,npart=npart_read)
##########################################
# read mass if needed
##########################################
if nzero != 0:
if mpi.mpi_IsMaster(): print "reading mass..."
massnzero = io.ReadDataBlock(f,float32,shape=(nzero,),byteorder=self.byteorder,pio=self.pio,npart=npart_m_read)
#mass = concatenate((mass,massnzero))
#if nzero==nbody: # this is maybe needed when pio='yes' Tue Feb 1 21:56:48 CET 2011
if nzero==nbody_tot:
mass = massnzero
else:
mass = array([])
i1 = 0
for i in range(len(npart)):
if npart[i]!=0: # particles belong to the class
if massarr[i] != 0:
mass = concatenate((mass,ones(npart[i])*massarr[i]))
else:
i2 = i1+npart[i]
mass = concatenate((mass,massnzero[i1:i2]))
i1 = i2
if i2!=len(massnzero):
raise "i2=",i2,"""!=len(massnzero)"""
if len(mass)!=nbody_tot:
raise "len(mass)=",len(mass),"!=nbody_tot"
'''
if massarr[i] == 0:
if npart[i]!=0:
if len(massnzero)!=npart[i]:
raise "this case is not taken into account, sorry !"
mass = concatenate((mass,massnzero))
else:
mass = concatenate((mass,ones(npart[i])*massarr[i]))
'''
# extentions
u = None
rho = None
rsp = None
opt = None
opt2= None
erd = None
dte = None
pot = None
tstar = None
minit = None
idp = None
metals = None
if not io.end_of_file(f,pio=self.pio,MPI=MPI):
if mpi.mpi_IsMaster(): print "reading u..."
u = io.ReadDataBlock(f,float32,shape=(ngas_read,),byteorder=self.byteorder,pio=self.pio,npart=None)
u = concatenate((u,zeros(nbody-ngas).astype(float32)))
if not io.end_of_file(f,pio=self.pio,MPI=MPI):
if mpi.mpi_IsMaster(): print "reading rho..."
rho = io.ReadDataBlock(f,float32,shape=(ngas_read,),byteorder=self.byteorder,pio=self.pio,npart=None)
rho = concatenate((rho,zeros(nbody-ngas).astype(float32)))
if not io.end_of_file(f,pio=self.pio,MPI=MPI):
if mpi.mpi_IsMaster(): print "reading rsp..."
rsp = io.ReadDataBlock(f,float32,shape=(ngas_read,),byteorder=self.byteorder,pio=self.pio,npart=None)
rsp = concatenate((rsp,zeros(nbody-ngas).astype(float32)))
# here it is the end of the minimal output
if flag_metals: # gas metal properties
if not io.end_of_file(f,pio=self.pio,MPI=MPI):
if mpi.mpi_IsMaster(): print "reading metals..."
metals = io.ReadDataBlock(f,float32,shape=(ngas_read,NELEMENTS),byteorder=self.byteorder,pio=self.pio,npart=None)
metals = concatenate((metals,zeros((nbody-ngas,NELEMENTS)).astype(float32)))
if flag_age: # stellar properties
if not io.end_of_file(f,pio=self.pio,MPI=MPI):
if mpi.mpi_IsMaster(): print "reading tstar..."
tstar = io.ReadDataBlock(f,float32,shape=(nstars_read,),byteorder=self.byteorder,pio=self.pio,npart=None)
tstar = concatenate((-1*ones(ngas).astype(float32),tstar,-1*ones(nbody-ngas-nstars).astype(float32)))
if not io.end_of_file(f,pio=self.pio,MPI=MPI):
if mpi.mpi_IsMaster(): print "reading minit..."
minit = io.ReadDataBlock(f,float32,shape=(nstars_read,),byteorder=self.byteorder,pio=self.pio,npart=None)
minit = concatenate((0*ones(ngas).astype(float32),minit,0*ones(nbody-ngas-nstars).astype(float32)))
if not io.end_of_file(f,pio=self.pio,MPI=MPI):
if mpi.mpi_IsMaster(): print "reading idp..."
idp = io.ReadDataBlock(f,int32,shape=(nstars_read,),byteorder=self.byteorder,pio=self.pio,npart=None)
idp = concatenate((-1*ones(ngas).astype(float32),idp,-1*ones(nbody-ngas-nstars).astype(float32)))
if not io.end_of_file(f,pio=self.pio,MPI=MPI):
if mpi.mpi_IsMaster(): print "reading rho_stars..."
rho_stars = io.ReadDataBlock(f,float32,shape=(nstars_read,),byteorder=self.byteorder,pio=self.pio,npart=None)
rho = concatenate((rho[:ngas],rho_stars,zeros(nbody-ngas-nstars).astype(float32)))
if not io.end_of_file(f,pio=self.pio,MPI=MPI):
if mpi.mpi_IsMaster(): print "reading rsp_stars..."
rsp_stars = io.ReadDataBlock(f,float32,shape=(nstars_read,),byteorder=self.byteorder,pio=self.pio,npart=None)
rsp = concatenate((rsp[:ngas],rsp_stars,zeros(nbody-ngas-nstars).astype(float32)))
if flag_metals: # stars metal properties
if not io.end_of_file(f,pio=self.pio,MPI=MPI):
metals_stars = io.ReadDataBlock(f,float32,shape=(nstars_read,NELEMENTS),byteorder=self.byteorder,pio=self.pio,npart=None)
metals = concatenate(( metals[:ngas,:], metals_stars, zeros((nbody-ngas-nstars,NELEMENTS)).astype(float32) ))
if not io.end_of_file(f,pio=self.pio,MPI=MPI):
opt = io.ReadDataBlock(f,float32,shape=(ngas_read,),byteorder=self.byteorder,pio=self.pio,npart=None)
if (len(opt)==ngas):
opt = concatenate((opt,zeros(nbody-ngas).astype(float32)))
if not io.end_of_file(f,pio=self.pio,MPI=MPI):
opt2 = io.ReadDataBlock(f,float32,shape=(ngas_read,),byteorder=self.byteorder,pio=self.pio,npart=None)
if (len(opt2)==ngas):
opt2 = concatenate((opt2,zeros(nbody-ngas).astype(float32)))
if not io.end_of_file(f,pio=self.pio,MPI=MPI):
erd = io.ReadDataBlock(f,float32,shape=(ngas_read,),byteorder=self.byteorder,pio=self.pio,npart=None)
if (len(erd)==ngas):
erd = concatenate((erd,zeros(nbody-ngas).astype(float32)))
if not io.end_of_file(f,pio=self.pio,MPI=MPI):
dte = io.ReadDataBlock(f,float32,shape=(ngas_read,),byteorder=self.byteorder,pio=self.pio,npart=None)
dte = concatenate((dte,zeros(nbody-ngas).astype(float32)))
if not io.end_of_file(f,pio=self.pio,MPI=MPI):
pot = io.ReadDataBlock(f,float32,shape=(ngas_read,),byteorder=self.byteorder,pio=self.pio,npart=None)
# make global
self.npart = npart
self.massarr = massarr
self.atime = atime
self.redshift = redshift
self.flag_sfr = flag_sfr
self.flag_feedback = flag_feedback
self.nall = nall
self.flag_cooling = flag_cooling
self.num_files = num_files
self.boxsize = boxsize
self.omega0 = omega0
self.omegalambda = omegalambda
self.hubbleparam = hubbleparam
self.flag_age = flag_age
self.flag_metals = flag_metals
self.nallhw = nallhw
self.flag_entr_ic = flag_entr_ic
self.flag_chimie_extraheader = flag_chimie_extraheader
self.critical_energy_spec = critical_energy_spec
self.empty = empty
self.nbody = nbody
self.pos = pos
self.vel = vel
self.mass = mass
self.num = num
self.tpe = array([],int32)
for i in range(len(npart)):
self.tpe = concatenate( (self.tpe,ones(npart[i])*i) )
self.nzero = nzero
self.u = u
self.rho = rho
self.tstar = tstar
self.minit = minit
self.idp = idp
self.metals = metals
self.opt = opt
self.opt2= opt2
self.erd = erd
self.dte = dte
self.pot = pot
self.rsp = rsp
if type(self.massarr) == ndarray:
self.massarr = self.massarr.tolist()
if type(self.nall) == ndarray:
self.nall = self.nall.tolist()
if type(self.nallhw) == ndarray:
self.nallhw = self.nallhw.tolist()
def write_particles(self,f):
'''
specific format for particle file
'''
# here, we must let the user decide if we creates
# the mass block or not, event if all particles have the same mass
massarr,nzero = self.get_massarr_and_nzero()
if self.pio == 'yes':
'''
here, we have to compute also
mass for each proc
'''
npart = self.npart
nall = self.npart_tot
num_files = mpi.NTask
npart_write = None
npart_m_write = None
else:
npart = self.npart
nall = self.npart_tot
npart_all = array(mpi.mpi_allgather(npart))
num_files = 1
npart_write = self.npart
npart_m_write = array(self.npart) * (array(self.massarr)==0)
# compute the global massarr and global nzero
nzero_tot = mpi.mpi_sum(nzero)
massarr_all = array(mpi.mpi_allgather(massarr))
massarr_tot = zeros(len(npart),float)
for i in range(len(npart)):
# keep only values where there are particles
massarr_all_red = compress(npart_all[:,i]!=0,massarr_all[:,i])
if len(massarr_all_red)>0:
if (massarr_all_red == massarr_all_red).all():
massarr_tot[i] = massarr_all[0,i]
else: # not equal
raise "this case is not implemented"
massarr_tot[i] = 0.0
nzero_tot = nzero_tot + sum(npart_write[:,i])
# now, re-compute nzero for the current node
massarr = massarr_tot
nzero = 0
for i in range(len(npart)):
if massarr[i] == 0:
nzero = nzero+npart[i]
# now that we have the right massarr and nzero,
# we can compute massnzero for each node
nzero_all = zeros((mpi.NTask,len(self.npart)))
if nzero != 0:
ni = 0
massnzero=array([],float32)
for i in range(len(self.npart)):
if npart[i]!=0 and massarr[i]==0.:
massnzero = concatenate((massnzero,self.mass[ni:ni+self.npart[i]]))
nzero_all[mpi.ThisTask,i] = nzero_all[mpi.ThisTask,i] + npart[i]
ni = ni + self.npart[i]
nzero_all = mpi.mpi_allreduce(nzero_all)
if self.pio == 'yes':
if nzero!=0 and len(massnzero)==0: # !!! because zere is a bug see warning in get_massarr_and_nzero
nzero = 0
else:
npart = self.npart_tot
nzero = mpi.mpi_allreduce(nzero) # to ensure that all nodes
# will do -> write mass if needed
# header
if sys.byteorder == self.byteorder:
npart = array(npart,int32).tostring()
massarr = array(massarr,float).tostring()
else:
npart = array(npart,int32).byteswap().tostring()
massarr = array(massarr,float).byteswap().tostring()
atime = self.atime
redshift = self.redshift
flag_sfr = self.flag_sfr
flag_feedback = self.flag_feedback
if sys.byteorder == self.byteorder:
nall = array(nall,int32).tostring()
else:
nall = array(nall,int32).byteswap().tostring()
flag_cooling = self.flag_cooling
num_files = num_files
boxsize = self.boxsize
omega0 = self.omega0
omegalambda = self.omegalambda
hubbleparam = self.hubbleparam
flag_age = self.flag_age
flag_metals = self.flag_metals
if sys.byteorder == self.byteorder:
nallhw = array(self.nallhw,float).tostring()
else:
nallhw = array(self.nallhw,float).byteswap().tostring()
flag_entr_ic = self.flag_entr_ic
flag_chimie_extraheader = self.flag_chimie_extraheader
critical_energy_spec = self.critical_energy_spec
empty = self.empty
# header
tpl = ((npart,24),(massarr,48),(atime,float),(redshift,float),(flag_sfr,int32),(flag_feedback,int32),(nall,24),(flag_cooling,int32),(num_files,int32),(boxsize,float),(omega0,float),(omegalambda,float),(hubbleparam,float),(flag_age,int32),(flag_metals,int32),(nallhw,24),(flag_entr_ic,int32),(flag_chimie_extraheader,int32),(critical_energy_spec,float),(empty,48))
io.WriteBlock(f,tpl,byteorder=self.byteorder)
# extra header
if self.flag_chimie_extraheader:
print "writing chimie extra-header..."
SolarAbundances = zeros(self.ChimieNelements,float32)
labels = ""
for i,elt in enumerate(self.ChimieElements):
SolarAbundances[i] = self.ChimieSolarAbundances[elt]
labels = labels + "%s,"%elt
labels_len = (256-4-self.ChimieNelements*4)
labels = labels + (labels_len-len(labels)) *" "
tpl = ((self.ChimieNelements,int32), (SolarAbundances,float32), (labels,len(labels)))
io.WriteBlock(f,tpl,byteorder=self.byteorder)
# positions
io.WriteArray(f,self.pos.astype(float32),byteorder=self.byteorder,pio=self.pio,npart=npart_write)
# velocities
io.WriteArray(f,self.vel.astype(float32),byteorder=self.byteorder,pio=self.pio,npart=npart_write)
# id
io.WriteArray(f,self.num.astype(int32),byteorder=self.byteorder,pio=self.pio,npart=npart_write)
print npart_write
# write mass if needed
if nzero != 0:
io.WriteArray(f,massnzero.astype(float32),byteorder=self.byteorder,pio=self.pio,npart=npart_m_write)
# write extension
if self.has_array('u'):
if self.u != None:
io.WriteArray(f,self.u[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None)
print "write u"
if self.has_array('rho'):
if self.rho != None:
io.WriteArray(f,self.rho[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None)
print "write rho"
if self.has_array('rsp'):
if self.rsp != None:
io.WriteArray(f,self.rsp[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None)
print "write rsp"
# this is the end of the minimal output
if self.flag_metals:
if self.has_array('metals'):
io.WriteArray(f,self.metals[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None)
print "write metals"
if self.flag_age:
if self.has_array('tstar'):
io.WriteArray(f,self.tstar[self.npart[0]:self.npart[0]+self.npart[1]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None)
print "write tstar"
if self.has_array('minit'):
io.WriteArray(f,self.minit[self.npart[0]:self.npart[0]+self.npart[1]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None)
print "write minit"
if self.has_array('idp'):
print
io.WriteArray(f,self.idp[self.npart[0]:self.npart[0]+self.npart[1]].astype(int32),byteorder=self.byteorder,pio=self.pio,npart=None)
print "write idp"
if self.has_array('rho'):
data = self.rho[self.npart[0]:self.npart[0]+self.npart[1]].astype(float32)
if len(data)>0:
io.WriteArray(f,data,byteorder=self.byteorder,pio=self.pio,npart=None)
print "write rho (stars)"
if self.has_array('rsp'):
data = self.rsp[self.npart[0]:self.npart[0]+self.npart[1]].astype(float32)
if len(data)>0:
io.WriteArray(f,data,byteorder=self.byteorder,pio=self.pio,npart=None)
print "write rsp (stars)"
if self.flag_metals:
if self.has_array('metals'):
data = self.metals[self.npart[0]:self.npart[0]+self.npart[1]].astype(float32)
if len(data)>0:
io.WriteArray(f,data,byteorder=self.byteorder,pio=self.pio,npart=None)
print "write metals (stars)"
if self.has_array('opt'):
if self.opt != None:
io.WriteArray(f,self.opt[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None)
if self.has_array('opt2'):
if self.opt2 != None:
io.WriteArray(f,self.opt2[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None)
if self.has_array('erd'):
if self.erd != None:
io.WriteArray(f,self.erd[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None)
if self.has_array('dte'):
if self.dte != None:
io.WriteArray(f,self.dte[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None)
if self.has_array('pot'):
if self.pot != None:
io.WriteArray(f,self.pot[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None)
def spec_info(self):
"""
Write spec info
"""
infolist = []
infolist.append("")
#infolist.append("nzero : %s"%self.nzero)
#infolist.append("npart : %s"%self.npart)
#infolist.append("massarr : %s"%self.massarr)
infolist.append("atime : %s"%self.atime)
infolist.append("redshift : %s"%self.redshift)
infolist.append("flag_sfr : %s"%self.flag_sfr)
infolist.append("flag_feedback : %s"%self.flag_feedback)
infolist.append("nall : %s"%self.nall)
infolist.append("flag_cooling : %s"%self.flag_cooling)
infolist.append("num_files : %s"%self.num_files)
infolist.append("boxsize : %s"%self.boxsize)
infolist.append("omega0 : %s"%self.omega0)
infolist.append("omegalambda : %s"%self.omegalambda)
infolist.append("hubbleparam : %s"%self.hubbleparam)
infolist.append("flag_age : %s"%self.flag_age)
infolist.append("flag_metals : %s"%self.flag_metals)
infolist.append("nallhw : %s"%self.nallhw)
infolist.append("flag_entr_ic : %s"%self.flag_entr_ic)
infolist.append("critical_energy_spec: %s"%self.critical_energy_spec)
infolist.append("")
if self.has_array('u'):
infolist.append("len u : %s"%len(self.u))
infolist.append("u[0] : %s"%self.u[0])
infolist.append("u[-1] : %s"%self.u[-1])
if self.has_array('rho'):
infolist.append("len rho : %s"%len(self.rho))
infolist.append("rho[0] : %s"%self.rho[0])
infolist.append("rho[-1] : %s"%self.rho[-1])
if self.has_array('rsp'):
infolist.append("len rsp : %s"%len(self.rsp))
infolist.append("rsp[0] : %s"%self.rsp[0])
infolist.append("rsp[-1] : %s"%self.rsp[-1])
if self.has_array('opt'):
infolist.append("len opt : %s"%len(self.opt))
infolist.append("opt[0] : %s"%self.opt[0])
infolist.append("opt[-1] : %s"%self.opt[-1])
if self.has_array('opt2'):
infolist.append("len opt2 : %s"%len(self.opt2))
infolist.append("opt2[0] : %s"%self.opt2[0])
infolist.append("opt2[-1] : %s"%self.opt2[-1])
if self.has_array('erd'):
infolist.append("len erd : %s"%len(self.erd))
infolist.append("erd[0] : %s"%self.erd[0])
infolist.append("erd[-1] : %s"%self.erd[-1])
if self.has_array('dte'):
infolist.append("len dte : %s"%len(self.dte))
infolist.append("dte[0] : %s"%self.dte[0])
infolist.append("dte[-1] : %s"%self.dte[-1])
if self.has_array('tstar'):
infolist.append("len tstar : %s"%len(self.tstar))
infolist.append("tstar[0] : %s"%self.tstar[0])
infolist.append("tstar[-1] : %s"%self.tstar[-1])
if self.has_array('idp'):
infolist.append("len idp : %s"%len(self.idp))
infolist.append("idp[0] : %s"%self.idp[0])
infolist.append("idp[-1] : %s"%self.idp[-1])
return infolist
#def select(self,tpe='gas'):
def select(self,*arg,**kw):
"""
Return an N-body object that contain only particles of a
certain type, defined by in gadget:
gas : gas particles
halo : halo particles
disk : disk particles
bulge : bulge particles
stars : stars particles
bndry : bndry particles
sph : gas with u > u_c
sticky : gas with u < u_c
"""
index = {'gas':0,'halo':1,'disk':2,'bulge':3,'stars':4,'bndry':5,'stars1':1,'halo1':2}
# this allows to write nb.select(('gas','disk'))
if len(arg)==1:
if type(arg[0])==types.TupleType:
arg = arg[0]
tpes = arg
# create the selection vector
c = zeros(self.nbody)
for tpe in tpes:
if type(tpe) == types.StringType:
if (tpe=='sph'):
c = c+(self.u>self.critical_energy_spec)*(self.tpe==0)
elif (tpe=='sticky'):
c = c+(self.u<self.critical_energy_spec)*(self.tpe==0)
elif (tpe=='diskbulge'):
c = c+(self.tpe==2)+(self.tpe==3)
elif (tpe=='wg'):
c = c+(self.u>0)*(self.tpe==0)
elif (tpe=='cg'):
c = c+(self.u<0)*(self.tpe==0)
elif (tpe=='all'):
return self
elif not index.has_key(tpe):
print "unknown type, do nothing %s"%(tpe)
return self
else:
i = index[tpe]
c = c+(self.tpe==i)
elif type(tpe) == types.IntType:
c = c+(self.tpe==tpe)
return self.selectc(c)
'''
elif type(tpe) == types.StringType:
if (tpe=='sph'):
nb = self.select('gas')
return nb.selectc((self.u>self.critical_energy_spec))
if (tpe=='sticky'):
nb = self.select('gas')
return nb.selectc((nb.u<=nb.critical_energy_spec))
if not index.has_key(tpe):
print "unknown type %s"%(tpe)
return self
i = index[tpe]
else:
i = tpe
if self.npart[i]==0:
#print "no particle of type %s"%(tpe)
return self.selectc(zeros(self.nbody))
n1 = sum(self.npart[:i])
n2 = n1 + self.npart[i]-1
return self.sub(n1,n2)
'''
def subdis(self,mode='dd',val=None):
"""
Equivalent of select
"""
return self.select(mode)
def Z(self):
"""
total metallicity
"""
return log10(self.metals[:,self.NELEMENTS-1] / self.ChimieSolarAbundances['Metals'] + 1.0e-20)
def Fe(self):
"""
metallicity Fe
"""
return log10(self.metals[:,0] / self.ChimieSolarAbundances['Fe'] + 1.0e-20)
def Mg(self):
"""
magnesium
"""
return log10(self.metals[:,1] / self.ChimieSolarAbundances['Mg'] + 1.0e-20)
def O(self):
"""
Oxygen
"""
return log10(self.metals[:,2] / self.ChimieSolarAbundances['O'] + 1.0e-20)
def Ba(self):
"""
Barium
"""
return log10(self.metals[:,3] / self.ChimieSolarAbundances['Ba'] + 1.0e-20)
def MgFe(self):
eps = 1e-20
MgFe = log10((self.metals[:,1]+eps)/(self.metals[:,0]+eps) / self.ChimieSolarAbundances['Mg'] * self.ChimieSolarAbundances['Fe'])
return MgFe
def BaFe(self):
eps = 1e-20
BaFe = log10((self.metals[:,3]+eps)/(self.metals[:,0]+eps) / self.ChimieSolarAbundances['Ba'] * self.ChimieSolarAbundances['Fe'])
return BaFe
def luminosity_spec(self,tnow=None):
"""
compute specific luminosity, using metalicity
u_mass = 1.e10
LuminosityElement*u_mass/1.0d6, ' x 10^6 Lsun'
"""
u_mass = 1.e10
u_time = 4.7287e6
SolarAbun_Fe = 0.001771
if self.tstar==None:
return array([],float32)
if tnow==None:
tnow = self.atime
Ages = (tnow-self.tstar)*u_time*1.0e-9 # ages in Gyr
Zs = self.Z()
# compute luminosities using LObj
L = self.LObj.Luminosities(Zs,Ages)
return L
def old_luminosity_spec(self,tnow=None):
"""
compute specific luminosity, using metalicity
u_mass = 1.e10
LuminosityElement*u_mass/1.0d6, ' x 10^6 Lsun'
"""
u_mass = 1.e10
u_time = 4.7287e6
SolarAbun_Fe = 0.001771
if self.tstar==None:
return 0
if tnow==None:
tnow = self.atime
age = (tnow-self.tstar)*u_time*1.0e-9
Fe = log10(self.metals[:,0] / SolarAbun_Fe + 1.0e-20)
met_bins = array([-1.68, -1.28,-0.68])
age_bins = array([0.007,0.02,0.07,0.10,0.11,0.13,0.14,0.16,0.18,0.20,0.22,0.25,0.28,
0.32,0.35,0.40,0.45,0.50,0.56,0.63,0.71,0.79,0.89,1.00,1.12,1.26,1.41,1.58,1.78,2.00,
2.24,2.51,2.82,3.16,3.55,3.98,4.47,5.01,5.62,6.31,7.08,7.94,8.91,10.00,11.22,12.59,14.13,15.85])
L0 = array([1./(0.003*2.0),1./(0.007*2.0),1./(0.054*2.0),1./0.172,1./0.122,1./0.126,1./0.132,1./0.138,1./0.145,
1./0.153,1./0.166,1./0.179,1./0.196,1./0.218,1./0.232,1./0.248,1./0.267,1./0.292,1./0.313,1./0.339,
1./0.371,1./0.397,1./0.429,1./0.464,1./0.471,1./0.461,1./0.597,1./0.640,1./0.681,1./0.746,1./0.813,
1./0.899,1./0.975,1./1.085,1./1.202,1./1.328,1./1.464,1./1.606,1./1.753,1./1.929,1./2.100,1./2.282,
1./2.485,1./2.762,1./2.987,1./3.230,1./3.504,1./3.794,1./4.100])
L1 = array([1./(0.003 * 2.0),1./(0.007 * 2.0),1./(0.054 * 2.0),1./0.154,1./0.126,1./0.134,1./0.141,1./0.149,
1./0.162,1./0.177,1./0.195,1./0.212,1./0.227,1./0.245,1./0.263,1./0.278,1./0.296,1./0.324,1./0.348,
1./0.376,1./0.412,1./0.444,1./0.481,1./0.516,1./0.550,1./0.543,1./0.590,1./0.688,1./0.747,1./0.819,
1./0.896,1./0.985,1./1.077,1./1.197,1./1.317,1./1.444,1./1.573,1./1.735,1./1.888,1./2.066,1./2.304,
1./2.502,1./2.729,1./3.008,1./3.256,1./3.564,1./3.740,1./4.212,1./4.491])
L2 = array([1./(0.003 * 2.0),1./(0.007 * 2.0),1./(0.054 * 2.0),1./0.163,1./0.169,1./0.159,1./0.172,1./0.184,
1./0.198,1./0.214,1./0.230,1./0.253,1./0.270,1./0.291,1./0.318,1./0.345,1./0.371,1./0.401,1./0.432,
1./0.464,1./0.501,1./0.540,1./0.583,1./0.633,1./0.684,1./0.710,1./0.694,1./0.870,1./0.989,1./1.092,
1./1.183,1./1.296,1./1.410,1./1.534,1./1.665,1./1.809,1./2.023,1./2.200,1./2.477,1./2.692,1./2.927,
1./3.214,1./3.490,1./3.783,1./4.087,1./4.411,1./4.742,1./5.115,1./5.515])
L3 = array([1./(0.003 * 2.0),1./(0.007 * 2.0),1./(0.054 * 2.0),1./0.166,1./0.177,1./0.179,1./0.192,1./0.206,
1./0.221,1./0.241,1./0.259,1./0.282,1./0.308,1./0.331,1./0.357,1./0.386,1./0.415,1./0.448,1./0.485,
1./0.524,1./0.562,1./0.614,1./0.671,1./0.725,1./0.785,1./0.770,1./0.937,1./1.068,1./1.189,1./1.318,
1./1.448,1./1.575,1./1.736,1./1.918,1./2.066,1./2.236,1./2.439,1./2.728,1./2.977,1./3.331,1./3.722,
1./4.063,1./4.387,1./4.717,1./5.047,1./5.474,1./5.972,1./6.636,1./7.018])
L = array([L0,L1,L2,L3])
#i = searchsorted(met_bins,Fe) # metallicity index
#j = searchsorted(age_bins,age) # age index
i = met_bins.searchsorted(Fe)
j = age_bins.searchsorted(age)
#return take(L,[i,j],[0,1])*self.mass
return L[i,j]
def luminosity(self,tnow=None):
return self.luminosity_spec(tnow)*self.mass
def Age(self):
'''in Gyrs (for treeasph units)'''
u_time = 4.7287e6
age = (self.atime-self.tstar)*u_time*1.0e-9
return age
#################################################################
# physical values
#################################################################
def Rho(self):
if self.comovingintegration:
print "Rho : converting to physical units (a=%5.3f h=%5.3f)"%(self.atime,self.hubbleparam)
return self.rho/self.atime**3*self.hubbleparam**2
else:
print "Rho : converting to physical units (h=%5.3f)"%(self.hubbleparam)
return self.rho*self.hubbleparam**2
def T(self):
gamma = self.unitsparameters.get('gamma')
xi = self.unitsparameters.get('xi')
ionisation = self.unitsparameters.get('ionisation')
mu = thermodyn.MeanWeight(xi,ionisation)
mh = ctes.PROTONMASS.into(self.localsystem_of_units)
k = ctes.BOLTZMANN.into(self.localsystem_of_units)
thermopars = {"k":k,"mh":mh,"mu":mu,"gamma":gamma}
T = where((self.u>0),thermodyn.Tru(self.Rho(),self.u,thermopars),0)
return T
def Tcool(self):
from pNbody import cooling
if self.metals == None:
FeH = zeros(self.nbody).astype(float32)
else:
FeH = self.metals[:,0]
#l = cooling.get_lambda_from_Density_EnergyInt_FeH(self.rho,self.u,FeH)
#dudt = l/self.rho
#tcool = self.u/dudt
tcool = cooling.get_cooling_time_from_Density_EnergyInt_FeH(self.rho,self.u,FeH)
return tcool.astype(float32)

Event Timeline

c4science · Help