Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88762967
plot_both.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, Oct 20, 13:50
Size
6 KB
Mime Type
text/x-python
Expires
Tue, Oct 22, 13:50 (2 d)
Engine
blob
Format
Raw Data
Handle
21817468
Attached To
rGEAR Gear
plot_both.py.svn-base
View Options
#!/usr/bin/env python
import string
from numpy import *
import sys
import pylab as pl
import libtreeasph
import libmetalsSNII
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
##################################
# poirier
##################################
coeff_z,Mmin,Mmax,PowSN1,PowSN2,ConstSN,Msn,ArrayOrigin,ArrayStep,ChemArray,Mco,MIa = read_chemistry('chimie.dat')
#b = 0.13466316207
#ChemArray = ChemArray/b
n = len(ChemArray)
n = 106
IFe = ChemArray[0:n,1] # Fe
IMg = ChemArray[0:n:,2] # Mg
IOx = ChemArray[0:n:,3] # Ox
IMe = ChemArray[0:n:,0] # metal
IRes = ChemArray[0:n:,4] # Matter ejection
IHco = ChemArray[0:n:,5] # Integral (( 1 - Helium core )*m^x)
M1 = 10**( ArrayOrigin[0]+ arange(n)*ArrayStep[0] )
M2 = 10**( ArrayOrigin[1]+ arange(n)*ArrayStep[1] )
##################################
# yves
##################################
ChAr,PaAr = libmetalsSNII.Read('chemyves.dat')
n = len(ChAr[0,:])
j = 0
MyRes = 10**( log10(PaAr[j,0])+ arange(n)* PaAr[j,1] )
IyRes = ChAr[j]
j = 1
MyHco = 10**( log10(PaAr[j,0])+ arange(n)* PaAr[j,1] )
IyHco = ChAr[j]
j = 2
MyFe = 10**( log10(PaAr[j,0])+ arange(n)* PaAr[j,1] )
IyFe = ChAr[j]
j = 3
MyMg = 10**( log10(PaAr[j,0])+ arange(n)* PaAr[j,1] )
IyMg = ChAr[j]
j = 4
MyOx = 10**( log10(PaAr[j,0])+ arange(n)* PaAr[j,1] )
IyOx = ChAr[j]
j = 5
MyMe = 10**( log10(PaAr[j,0])+ arange(n)* PaAr[j,1] )
IyMe = ChAr[j]
##################################
# plot
##################################
pt.plot(M1,IRes,'k')
pt.plot(M1,IHco,'r')
pt.plot(M2,IFe*500,'y')
pt.plot(M2,IMg*500,'g')
pt.plot(M2,IOx*10,'b')
pt.plot(M2,IMe*10,'c')
pt.plot(MyRes,IyRes,'k--')
pt.plot(MyHco,IyHco,'r--')
pt.plot(MyFe,IyFe*500,'y--')
pt.plot(MyMg,IyMg*500,'g--')
pt.plot(MyOx,IyOx*10,'b--')
pt.plot(MyMe,IyMe*10,'c--')
pt.semilogx()
pt.show()
Event Timeline
Log In to Comment