Homec4science
Diffusion Gear (master) Browse

Gear/PyChem/shared/tables/MetalSNII/.svn/text-base/libmetalsSNII.py.svn-base
73ca2cd44116master

Last Change
Raw File

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()
c4science · Help