libmetalsSNII.py
73ca2cd44116master

libmetalsSNII.py
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='#'):
			#####################################################
			"""
			Read a file containing metal information:


			:param fd: file name
			:type fd: string or file handler

			:param columns: columns to read
			:type columns: python list

			:param lines: lines to read
			:type lines: python list

			:param dtype: output dtype
			:type dtype: dtype

			:param skipheader: skip header
			:type skipheader: bool

			:param cchar: comment charater
			:type cchar: string

			:returns: data


			``[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``


			"""


			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):
			"""Write the SNIa parameters in a file:


			:param f: file
			:type f: file handler

			:param Mmin: min SNII mass
			:type Mmin: float

			:param Mmax: max SNII mass
			:type Mmax: float

			:param m_s: IMF masses
			:type m_s: ndarray

			:param a_s: IMF slopes
			:type a_s: ndarray

			:param MetalFiles: file name of the file containing yields information
			:type MetalFiles: string

			:param HeliumCoreFile: file name of the file containing helium core information
			:type HeliumCoreFile: string

			:param elts: list of elements
			:type elts: python list

			:param nelts: number of elements
			:type nelts: int

			:param Mmin1: min star mass for yields interpolation
			:type Mmin1: float

			:param Mmax1: min SNII mass for yields interpolation
			:type Mmax1: float

			:param Mmin2: min star mass for helium core interpolation
			:type Mmin2: float

			:param Mmax2: max star mass for helium core interpolation
			:type Mmax2: float

			:param dM: mass steps
			:type dM: float

			:param n: ???
			:type n: int


			:returns: elts
			"""


			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
			print "nelts(=%d) != len(elts)(=%d) !!!"%(nelts,len(elts))
			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):
			"""
			Write a block of yields
			"""


			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 yields file and compute a linear interpolation of the values in
			the range Mmin1,Mmax1:

			:param MetalFiles: file name of the file containing yields information
			:type MetalFiles: string

			:param Mmin1: min star mass for the interpolation
			:type Mmin1: float

			:param Mmax1: max star mass for the interpolation
			:type Mmax1: float

			:param dM: mass steps
			:type dM: float

			:returns Mi2,Resi,Hcoi:

			.. warning:: values in solar mass are converted in mass fraction


			"""




			# read ejecta
			elts,data = read_MetalFile(MetalFile,cchar='#')

			# devide by mass in order to get a mass fraction
			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 helium core file and compute a linear interpolation of the values in
			the range Mmin2,Mmax2:

			:param HeliumCoreFile: file name of the file containing helium core information
			:type HeliumCoreFile: string

			:param Mmin2: min star mass for the interpolation
			:type Mmin2: float

			:param Mmax2: max star mass for the interpolation
			:type Mmax2: float

			:param dM: mass steps
			:type dM: float

			:returns Mi2,Resi,Hcoi:


			'''

			# 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 the 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():
			"""
			Do some plotting
			"""

			###############################################
			# 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()

Gear/PyChem/PyChem/libmetalsSNII.py73ca2cd44116master

libmetalsSNII.pyView Options

Gear/PyChem/PyChem/libmetalsSNII.py
73ca2cd44116master

libmetalsSNII.py
View Options