Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F111679218
plot_P.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
Sun, May 4, 09:06
Size
5 KB
Mime Type
text/x-python
Expires
Tue, May 6, 09:06 (2 d)
Engine
blob
Format
Raw Data
Handle
25967743
Attached To
rGEAR Gear
plot_P.py.svn-base
View Options
#!/usr/bin/env python
import string
from numpy import *
import sys
import pylab as pl
import libtreeasph
import Ptools as pt
unit_length = 1.0
unit_mass = 1.e10
u_lt = -log10( 4.7287e11*sqrt(unit_length**3/unit_mass))
NbElement=4
feedback_energy = 1.1682615e13*(unit_length/unit_mass**2) # Energy released by each supernova by solar mass (10^(51)erg (in TREEASPH units))
# UnitLength_in_cm = 3.085e+21
# UnitMass_in_g = 1.989e+43
# UnitTime_in_s = 148849920000000.0
# feedback_energy = 10**51 /UnitMass_in_g * (UnitTime_in_s/UnitLength_in_cm)**2
def read_chemistry(file):
f = open(file)
######################
# Livetime
######################
f.readline()
f.readline()
coeff_z1 = array(string.split(f.readline()))
coeff_z2 = array(string.split(f.readline()))
coeff_z3 = array(string.split(f.readline()))
coeff_z = concatenate( (coeff_z1,coeff_z2,coeff_z3))
coeff_z.shape = (3,3)
coeff_z = coeff_z.astype(float)
# Conversion into program time unit
coeff_z[3-1,3-1] = coeff_z[3-1,3-1] + u_lt
# Transform log(timelive) polynomial expression log(t) = a(z) log(m)^2 + b(z) log(m) + c(z)
# into log(t) = a(z) log(m)^2 + 2b(z) log(m) + c(z)
coeff_z[2-1,:] = coeff_z[2-1,:] / 2.0
######################
# Mass interval
######################
f.readline()
f.readline()
f.readline()
Mmin = float(f.readline())
Mmax = float(f.readline())
LogMmax2 = log10(Mmax)
TwoLogMmax = 2.0 * LogMmax2
LogMmax2 = LogMmax2**2
######################
# IMFs Power
######################
f.readline()
f.readline()
f.readline()
PowSN1 = float(f.readline())
PowSN2 = float(f.readline())
######################
# Mass Parameters for SNe Rates
######################
f.readline()
f.readline()
f.readline()
Msn1 = array(string.split(f.readline()))
Msn2 = array(string.split(f.readline()))
Msn3 = array(string.split(f.readline()))
Msn = concatenate( (Msn1,Msn2,Msn3))
Msn.shape = (3,2)
Msn = Msn.astype(float)
f.readline()
Mco = float(f.readline()) / unit_mass
######################
# b Parameters for SNIa Rates
######################
f.readline()
f.readline()
f.readline()
b1 = float(f.readline())
b2 = float(f.readline())
ConstSN = array([0,0,0],float)
ConstSN[3-1] = (1.0+PowSN1)/(Mmax**(1.0+PowSN1)-Mmin**(1.0+PowSN1))
ConstSN[1-1] = (1.0+PowSN2)/(Msn[1-1,2-1]**(1.0+PowSN2)-Msn[1-1,1-1]**(1.0+PowSN2))
ConstSN[2-1] = (1.0+PowSN2)/(Msn[2-1,2-1]**(1.0+PowSN2)-Msn[2-1,1-1]**(1.0+PowSN2))
ConstSN[1-1] = b1*ConstSN[1-1]*ConstSN[3-1]
ConstSN[2-1] = b2*ConstSN[2-1]*ConstSN[3-1]
######################
# Array size and step
######################
f.readline()
f.readline()
f.readline()
ChemArraySize = int(f.readline())
######################
# >>>Matter<<<< ejection and Helium core
######################
ArrayOrigin = array([0,0],float)
ArrayStep = array([0,0],float)
for p in [0,1]:
f.readline()
f.readline()
f.readline()
line = array(string.split(f.readline()))
ArrayOrigin[p] = float(line[0])
ArrayStep[p] = float(line[1])
ArrayOrigin = log10(ArrayOrigin)
######################
# Matter ejection for IMF power x=-1.35
# Integral (( 1 - Helium core )*m^x) for IMF power x=-1.35
######################
'''
ChemArray[:,0] metal
ChemArray[:,1] Fe
ChemArray[:,2] Mg
ChemArray[:,3] Ox
ChemArray[:,4] Matter ejection
ChemArray[:,5] Integral (( 1 - Helium core )*m^x)
'''
ChemArray = zeros((110,NbElement+2),float)
MIa = zeros(NbElement,float)
for p in [0,1]:
f.readline()
f.readline()
f.readline()
for j in xrange(ChemArraySize):
ChemArray[j,NbElement+p] = float(f.readline()) * ConstSN[3-1]
#########################
# Number of elements ####
#########################
f.readline()
f.readline()
f.readline()
nelt = int(f.readline())
if nelt != NbElement:
raise "nelt != NbElement"
#########################
# Global metals ejection for IMF power x=-1.35
#########################
for p in xrange(NbElement):
f.readline()
f.readline()
f.readline()
for j in xrange(ChemArraySize):
ChemArray[j,p] = float(f.readline()) * ConstSN[3-1]
f.readline()
f.readline()
f.readline()
MIa[p] = float(f.readline())
MIa[p] = MIa[p] / unit_mass
ConstSN[1-1] = ( ConstSN[1-1]/PowSN2/PowSN1 ) * unit_mass
ConstSN[2-1] = ( ConstSN[2-1]/PowSN2/PowSN1 ) * unit_mass
ConstSN[3-1] = ( ConstSN[3-1]/PowSN1 ) * unit_mass
return coeff_z,Mmin,Mmax,PowSN1,PowSN2,ConstSN,Msn,ArrayOrigin,ArrayStep,ChemArray,Mco,MIa
coeff_z,Mmin,Mmax,PowSN1,PowSN2,ConstSN,Msn,ArrayOrigin,ArrayStep,ChemArray,Mco,MIa = read_chemistry('chimie.dat')
'''
ChemArray[:,0] metal
ChemArray[:,1] Fe
ChemArray[:,2] Mg
ChemArray[:,3] Ox
ChemArray[:,4] Matter ejection
ChemArray[:,5] Integral (( 1 - Helium core )*m^x)
'''
#b = 0.13466316207
#ChemArray = ChemArray/b
n = len(ChemArray)
n = 106
IFe = ChemArray[0:n,1]
IMg = ChemArray[0:n:,2]
IOx = ChemArray[0:n:,3]
IMe = ChemArray[0:n:,0]
IRes = ChemArray[0:n:,4]
IHco = ChemArray[0:n:,5]
M1 = 10**( ArrayOrigin[0]+ arange(n)*ArrayStep[0] )
M2 = 10**( ArrayOrigin[1]+ arange(n)*ArrayStep[1] )
print ArrayStep[0]
pt.plot(M1,IRes)
pt.plot(M1,IHco)
pt.plot(M2,IFe*500)
pt.plot(M2,IMg*500)
pt.plot(M2,IOx*10)
pt.plot(M2,IMe*10)
pt.semilogx()
pt.show()
Event Timeline
Log In to Comment