Homec4science
Diffusion Gear (master) Browse

Gear/PyChem/shared/tables/MetalSNIa/.svn/text-base/libmetalsSNIa.py.svn-base
73ca2cd44116master

Last Change
Raw File

libmetalsSNIa.py.svn-base
View Options

from numpy import *
import sys,types,string
#####################################################
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,SNIaFile,MeanWDMass,elts,nelts):
# read ejecta
eltbs,datai = read_MetalFile(SNIaFile,[0,1,2,3],skipheader=False,cchar='#')
if elts==None:
elts=eltbs
data = datai
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(1)
for key in datai.keys():
if key[:len(elt)] == elt:
dataiso = dataiso+datai[key]
print "add",key,"to",elt
dataif[elt] = dataiso
data = dataif
# check length of elts
if len(elts)!=nelts:
raise "nelts(=%d) != len(elts)(=%d) !!!"%(nelts,len(elts))
f.write("\n")
f.write("# Mean WD Mass\n")
f.write("\n")
f.write("%f\n"%MeanWDMass)
f.write("\n")
f.write("# SNIa Metal Ejection\n")
f.write("\n")
f.write("# Number of elements\n")
f.write("%d\n"%nelts)
f.write("\n")
f.write("# Ej\n")
f.write("0.0\n")
f.write("\n")
f.write("# Ejnp\n")
f.write("0.0\n")
f.write("\n")
for elt in elts:
f.write("# %s\n"%elt)
f.write("%e\n"%data[elt])
f.write("\n")
return elts
c4science · Help