Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F106411039
libmetalsSNII.py.svn-base
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Mar 25, 13:09
Size
8 KB
Mime Type
text/x-python
Expires
Thu, Mar 27, 13:09 (2 d)
Engine
blob
Format
Raw Data
Handle
25179323
Attached To
rGEAR Gear
libmetalsSNII.py.svn-base
View Options
#!/usr/bin/env python
import Ptools as pt
from pNbody import myNumeric
from numpy import *
from scipy.integrate import cumtrapz
import string,types,sys
import libimf
#####################################################
def read_MetalFile(file,columns=None,lines=None,dtype=float,skipheader=False,cchar='#'):
#####################################################
"""[X,Y,Z] = READ('FILE',[1,4,13],lines=[10,1000])
Read columns 1,4 and 13 from 'FILE' from line 10 to 1000
into array X,Y and Z
file is either fd or name file
"""
def RemoveComments(l):
if l[0]==cchar:
return None
else:
return l
def toNumList(l):
return map(dtype,l)
if type(file) != types.FileType:
f = open(file,'r')
else:
f = file
# read header while there is one
while 1:
fpos = f.tell()
header = f.readline()
if header[0] != cchar:
f.seek(fpos)
header = None
break
else:
if skipheader:
header = None
else:
# create dict from header
header = string.strip(header[2:])
elts = string.split(header)
break
# now, read the file content
lines = f.readlines()
# remove trailing
lines = map(string.strip, lines)
# remove comments
#lines = map(RemoveComments, lines)
# split
lines = map(string.split, lines)
# convert into float
lines = map(toNumList, lines)
# convert into array
lines = array(map(array, lines))
# transpose
lines = transpose(lines)
if header != None:
iobs = {}
i = 0
for elt in elts:
iobs[elt]=i
i = i + 1
vals = {}
for key in iobs.keys():
vals[key] = lines[iobs[key]]
return elts,vals
# return
if columns == None:
return lines
else:
return lines.take(axis=0,indices=columns)
def WriteParams(f,Mmin,Mmax,m_s,a_s,MetalFiles,HeliumCoreFile,elts,nelts,Mmin1,Mmax1,Mmin2,Mmax2,dM,n):
Mi2,Resi,Hcoi = ReadAndInterpolateHeliumCore(HeliumCoreFile,Mmin2,Mmax2,dM)
datai = {}
eltbs = []
for MetalFile in MetalFiles:
Mi1_f,datai_f,elts_f = ReadAndInterpolate(MetalFile,Mmin1,Mmax1,dM)
Mi1 = Mi1_f
for elt in elts_f:
eltbs.append(elt)
datai[elt] = datai_f[elt]
# check if we have an entry for Metals
try:
eltbs.index('Metals')
compute_metals=False
except ValueError:
compute_metals=True
# compute metals
# here, we sum the contribution of all elements
if compute_metals:
print "compute Metals :"
metals = zeros(len(Mi1))
for elt in eltbs:
print "add ",elt
metals = metals + datai[elt]
print "--------"
eltbs.append('Metals')
datai['Metals'] = metals
if elts==None:
elts=eltbs
else:
nelts = len(elts)
# create a new data and elts
dataif = {}
for elt in elts:
# check if exists
if datai.has_key(elt):
dataif[elt] = datai[elt]
else:
# try to find isotopes
dataiso = zeros(len(Mi1))
for key in datai.keys():
if key[:len(elt)] == elt:
dataiso = dataiso+datai[key]
print "add",key,"to",elt
dataif[elt] = dataiso
datai = dataif
## check length of elts
if len(elts)!=nelts: # nelts-1 because metals is not present
raise "nelts(=%d) != len(elts)(=%d) !!!"%(nelts,len(elts))
# output steps
lStep1 = (log10(Mmax*1.25) - log10(Mmin1) )/float(n-1)
lStep2 = (log10(Mmax*1.25) - log10(Mmin2) )/float(n-1)
Ms1 = 10**(arange(n)*lStep1 + log10(Mmin1))
Ms2 = 10**(arange(n)*lStep2 + log10(Mmin2))
datas,Ress, Hcos = ComputeEjectats(Mi1,datai,elts, Mi2,Resi,Hcoi, Ms1,Ms2)
Idatas,IRess,IHcos = ComputeIntegrals(Mi1,datai,elts, Mi2,Resi,Hcoi ,Mmin,Mmax,m_s,a_s, Ms1,Ms2)
label = "#### Metal Parameters ####\n"
eltList = [ ("Ej",Ress,Mmin2,lStep2),("Ejnp",Hcos,Mmin2,lStep2) ]
for elt in elts:
eltList.append( (elt,datas[elt],Mmin1,lStep1) )
WriteOutput(f,label,eltList, n,nelts)
f.write("\n")
label = "#### Metal Ejecta Parameters ####\n"
eltList = [ ("IEj",IRess,Mmin2,lStep2),("IEjnp",IHcos,Mmin2,lStep2) ]
for elt in elts:
eltList.append( ("I%s"%elt,Idatas[elt],Mmin1,lStep1) )
WriteOutput(f,label,eltList, n,nelts)
return elts
def WriteOutput(f,label,eltsList,n,nelts):
f.write(label)
f.write("\n# Number of points / Number of elements\n")
f.write("%d %d\n"%(n,nelts))
for elt,data,Mmin,lStep in eltsList:
f.write("\n# %s\n"%elt)
f.write("\n%g %g\n"%(Mmin,lStep))
for i in range(len(data)):
f.write("%e\n"%data[i])
def ReadAndInterpolate(MetalFile,Mmin1,Mmax1,dM):
'''
Read data
and compute linear interpolation, using
Mmin1,dM
Mmin2,dM
'''
# read ejecta
elts,data = read_MetalFile(MetalFile,cchar='#')
# devide by mass ...
for elt in elts[1:]:
data[elt] = data[elt]/data['M']
# ... and remove mass from dict
M1 = data['M']
del data['M']
elts = elts[1:]
# linear interpolation
Mi1 = arange(Mmin1,Mmax1+dM,dM)
# interpolate elements
datai = {}
for elt in elts:
datai[elt] = myNumeric.lininterp1d(Mi1.astype(float32),M1.astype(float32),data[elt].astype(float32))
datai[elt] = where(datai[elt]<=0,0,datai[elt])
return Mi1,datai,elts
def ReadAndInterpolateHeliumCore(HeliumCoreFile,Mmin2,Mmax2,dM):
'''
Read data
and compute linear interpolation, using
Mmin1,dM
Mmin2,dM
return *i
'''
# read residual and helium core
M2,Res,Hco = pt.io.read_ascii(HeliumCoreFile,[0,1,2],skipheader=True,cchar='#')
Res = Res/M2
Hco = Hco/M2
# linear interpolation
Mi2 = arange(Mmin2,Mmax2+dM,dM)
# interpolate
# interpolate residual and helium core
Resi = myNumeric.lininterp1d(Mi2.astype(float32),M2.astype(float32),Res.astype(float32))
Hcoi = myNumeric.lininterp1d(Mi2.astype(float32),M2.astype(float32),Hco.astype(float32))
Resi = where(Resi<=0,0,Resi)
Hcoi = where(Hcoi<=0,0,Hcoi)
return Mi2,Resi,Hcoi
def ComputeEjectats(Mi1,datai,elts, Mi2,Resi,Hcoi, Ms1,Ms2):
'''
From the interpolated values *i,
return interpolation of ejectats
'''
###############################################
# output
###############################################
Ress = myNumeric.lininterp1d(Ms2.astype(float32),Mi2.astype(float32),(1-Resi).astype(float32))
Hcos = myNumeric.lininterp1d(Ms2.astype(float32),Mi2.astype(float32),(1-Hcoi).astype(float32))
datas = {}
for elt in elts:
datas[elt] = myNumeric.lininterp1d(Ms1.astype(float32),Mi1.astype(float32),datai[elt].astype(float32))
return datas,Ress,Hcos
def ComputeIntegrals(Mi1,datai,elts, Mi2,Resi,Hcoi, Mmin,Mmax,m_s,a_s, Ms1,Ms2):
'''
From the interpolated values *i,
integrate and return interpolation
'''
# imf
b_s = libimf.ComputeConstants(m_s,a_s,Mmin,Mmax)
imf1 = libimf.imfv(Mi1,array(m_s),array(a_s),array(b_s),Mmin,Mmax)
imf2 = libimf.imfv(Mi2,array(m_s),array(a_s),array(b_s),Mmin,Mmax)
Idata = {}
for elt in elts:
Idata[elt] = cumtrapz(imf1 * datai[elt],Mi1)
IRes = cumtrapz(imf2 * (1-Resi),Mi2)
IHco = cumtrapz(imf2 * (1-Hcoi),Mi2)
###############################################
# output
###############################################
Idatas = {}
for elt in elts:
Idatas[elt] = myNumeric.lininterp1d(Ms1.astype(float32),Mi1.astype(float32),Idata[elt].astype(float32))
IRess = myNumeric.lininterp1d(Ms2.astype(float32),Mi2.astype(float32),IRes.astype(float32))
IHcos = myNumeric.lininterp1d(Ms2.astype(float32),Mi2.astype(float32),IHco.astype(float32))
return Idatas,IRess,IHcos
def Plot():
###############################################
# plot
###############################################
# plot 1
pt.subplot(211)
pt.plot(M2,Res,'ko')
pt.plot(Mi2,Resi,'k-')
pt.plot(M2,Hco,'ro')
pt.plot(Mi2,Hcoi,'r-')
pt.plot(M1,Fe,'yo')
pt.plot(Mi1,Fei,'y-')
pt.plot(M1,Mg,'go')
pt.plot(Mi1,Mgi,'g-')
pt.plot(M1,Ox,'bo')
pt.plot(Mi1,Oxi,'b-')
pt.plot(M1,Me,'co')
pt.plot(Mi1,Mei,'c-')
pt.semilogx()
pt.semilogy()
pt.axis([Mmin,Mmax1,1e-4,1e0])
# plot 2
pt.subplot(212)
pt.plot(Mi1[1:],IFe*500,'y')
pt.plot(Mi1[1:],IMg*500,'g')
pt.plot(Mi1[1:],IOx*10,'b')
pt.plot(Mi1[1:],IMe*10,'c')
pt.plot(Mi2[1:],IRes,'k')
pt.plot(Mi2[1:],IHco,'r')
'''
pt.plot(Ms1,IFes*500,'o')
pt.plot(Ms1,IMgs*500,'o')
pt.plot(Ms1,IOxs*10 ,'o')
pt.plot(Ms1,IMes*10 ,'o')
pt.plot(Ms2,IRess,'o')
pt.plot(Ms2,IHcos,'o')
'''
pt.semilogx()
pt.axis([Mmin,Mmax1,0.,4*b])
pt.show()
Event Timeline
Log In to Comment