Page MenuHomec4science
Files F102081559

plotElt.py
No OneTemporary
Actions

File Metadata

Created
Sun, Feb 16, 22:21

plotElt.py
View Options

#!/usr/bin/env python
import Ptools as pt
from pNbody import myNumeric
from numpy import *
from scipy.integrate import cumtrapz
from scipy.interpolate import interp2d
import types
import string
import pickle
#####################################################
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)
files = ["karakas-0.0001.txt",
"karakas-0.0040.txt",
"karakas-0.0080.txt",
"karakas-0.0200.txt"]
Zs = array([0.0001,0.0040,0.0080,0.0200])
#files = ["karakas-0.0001.txt"]
f=open("/tmp/ab.dmp",'r')
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d, Axes3D #<-- Note the capitalization!
fig = plt.figure()
ax = Axes3D(fig) #<-- Note the difference from your original code...
for i,file in enumerate(files):
# read ejecta
elts,data = read_MetalFile(file,cchar='#')
# devide by mass in order to get a mass fraction
for elt in elts[1:]:
data[elt] = data[elt]/data['M']
#pt.scatter(data['M'],data['Na'])
#pt.plot(data['M'],data['Na'],label=file)
ax.scatter( ones(len(data['Na']))*log10(Zs[i]), log10(data['M']) ,data['Na'])
# linear interp of the 4 files
Mmin =pickle.load(f)
lStep1 =pickle.load(f)
n =pickle.load(f)
datas =pickle.load(f)
Ms = 10**( arange(n)*lStep1 + log10(Mmin) )
#pt.plot(Ms,datas)
# get the final matrix
nZ =pickle.load(f)
lZmin =pickle.load(f)
dZ =pickle.load(f)
datasZsf =pickle.load(f)
#for i in xrange(nZ):
i=5
#pt.scatter( Ms, datasZsf[:,i],marker='x')
#pt.loglog()
#pt.axis([0.05 ,100 ,1e-7,1e-4 ])
#pt.legend()
#pt.show()
X = arange(nZ)*dZ + lZmin
Y = log10(Ms)
X, Y = meshgrid(X, Y)
Z = datasZsf
#cset = ax.contour(X, Y, Z, 16, extend3d=True)
ax.plot_wireframe(X, Y, Z, rstride=10, cstride=10)
plt.show()

Event Timeline

c4science · Help