Homec4science
Diffusion Gear (master) Browse

Gear/PyCool/tables_wiersma/scripts/plot_cooling_hdf5.py
d89b118970f3master

Last Change
Raw File

plot_cooling_hdf5.py
View Options

###########################################
#
# Read and plot data from CLOUDY tables
# using the hdf5 library
# Jeremie Despraz, 2012
#
# Usage: python2.7 plot_cooling_hdf5.py
#
###########################################
# to use IO utilities
from StringIO import StringIO
# to read and iterate over all cloudy output files
import os
# import the numpy library (for loadtxt)
from numpy import *
# import pylab (for matplotlib library)
from pylab import *
# import hdf library (for hdf5 handling)
import h5py
print 'python script starting'
# load data from all hdf5 files
dir = '../coolingtables/'
fnames = []
for file in os.listdir(dir):
if(file.endswith('.hdf5')):
fnames.append(file)
# plot only every 4th data file to avoid memory overflow
for k in range(len(fnames)/4):
path = dir + fnames[4*k]
table = h5py.File(path)
#print 'table keys:', table.keys()
#print 'table values:', table.values()
#print 'table items:', table.items()
# compute stuff
# net_cooling = Lambda / n_H^2
# electron_density_over_n_h = n_e/n_h
Lambda = table[u'Solar'][u'Net_cooling'][:,:] / table[u'Solar'][u'Electron_density_over_n_h'][:,:]
sumLambda = ones(len(Lambda))
for i in range(len(Lambda)):
sumLambda[i] = sum(Lambda[i,:])
# reproduce fig.1 form Wiersma, R. P. C., Schaye, J., and Smith, B. D. 2009, MNRAS, 393, 99
loglog(table[u'Solar'][u'Temperature_bins'],sumLambda)
hold(True)
#axis([1e4,1e8,1e-26,1e-20])
xlim([1e4,1e8])
xlabel('T [K]')
#legend((fnames[k],),loc='best')
title('Normalized, net cooling rate (pos for cooling, neg for heating) as a function\n of H density and T for solar abundances (Lambda/n_H^2 [erg s^-1 cm^3])\n')
#show()
legend(fnames[0:-1:4],loc='best') #(form 0 to end by steps of size 4)
show()
# compute stuff
path = dir + fnames[0]
table = h5py.File(path)
Lambda = table[u'Solar'][u'Net_cooling'][:,:] / table[u'Solar'][u'Electron_density_over_n_h'][:,:]
meanLambda = ones(len(Lambda))
for i in range(len(Lambda)):
meanLambda[i] = mean(Lambda[i,:])
# plot data
plot(table[u'Solar'][u'Temperature_bins'][:],table[u'Solar'][u'Net_cooling'][:,1:-1]) #last element = -1
xlabel('T [K]')
legend(('Lambda',),loc='best')
title('Normalized, net cooling rate (pos for cooling, neg for heating) as a function\n of H density and T for solar abundances (Lambda/n_H^2 [erg s^-1 cm^3])\n')
show()
loglog(table[u'Solar'][u'Temperature_bins'][:],Lambda[:,1:-1]) #logarithmic plot
xlim([1e4,1e8])
xlabel('T [K]')
legend(('log(Lambda)',),loc='best')
title('Normalized, net cooling rate (pos for cooling, neg for heating) as a function\n of H density and T for solar abundances (Lambda/n_H^2 [erg s^-1 cm^3])\n')
show()
loglog(table[u'Solar'][u'Temperature_bins'],meanLambda) #logarithmic plot
xlim([1e4,1e8])
xlabel('T [K]')
legend(('mean(Lambda)',),loc='best')
title('Normalized, net cooling rate (pos for cooling, neg for heating) as a function\n of H density and T for solar abundances (Lambda/n_H^2 [erg s^-1 cm^3])\n')
show()
print 'python script terminated'
c4science · Help