Homec4science
Diffusion Mtools (master) Browse

mtools/examples/gasNvsFe.py
b97ef7bb51aamaster

Last Change
Raw File

gasNvsFe.py
View Options

#!/usr/bin/env python
from optparse import OptionParser
from numpy import histogram as histo
import Ptools as pt
from pNbody import *
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = pt.add_postscript_options(parser)
parser = pt.add_ftype_options(parser)
parser = pt.add_reduc_options(parser)
parser = pt.add_center_options(parser)
parser = pt.add_select_options(parser)
parser = pt.add_cmd_options(parser)
parser = pt.add_display_options(parser)
parser = pt.add_info_options(parser)
parser = pt.add_limits_options(parser)
parser = pt.add_log_options(parser)
parser.add_option("--data",
action="store",
type="string",
dest="data",
default = None,
help="data")
(options, args) = parser.parse_args()
pt.check_files_number(args)
files = args
return files,options
def add_data_FeHistogram(data):
if data == 'scl':
Fe,f = pt.io.read_ascii('Sculptor-Fe-histogram.txt',[0,2])
pt.plot(Fe, f, 'r--', linewidth=1,c='r')
elif data == 'fnx':
Fe,f = pt.io.read_ascii('Fornax-Fe-histogram.txt',[0,2])
pt.plot(Fe, f, 'r--', linewidth=1,c='r')
elif data == 'car':
Fe,f = pt.io.read_ascii('Carina-Fe-histogram.txt',[0,2])
pt.plot(Fe, f, 'r--', linewidth=1,c='r')
elif data == 'sex':
Fe,f = pt.io.read_ascii('Sextans-Fe-histogram.txt',[0,2])
pt.plot(Fe, f, 'r--', linewidth=1,c='r')
elif data == 'leoII':
Fe,f = pt.io.read_ascii('LeoII-Fe-histogram.txt',[0,2])
pt.plot(Fe, f, 'r--', linewidth=1,c='r')
#######################################
# MakePlot
#######################################
def MakePlot(files,ps):
class Option:
def __init__(self):
pass
opt = Option()
opt.size = (5.12,5.12)
opt.cmd = "nb.atime=3000;nb.histocenter();nb=nb.selectc(nb.rxyz()<10)"
opt.xmin = -4.
opt.xmax = -0.5
opt.ymin = 0
opt.ymax = 0.15
opt.nxh = 0
opt.select ="gas"
opt.nopoints = False
opt.accumulate = False
opt.density = False
opt.log = None
opt.colorbar = False
# other default options
opt.ftype = "gadget"
opt.reduc = None
opt.center = None
opt.info = None
opt.display = None
opt.data = None
opt.ps = ps
# init plot
pt.InitPlot([],opt)
pt.pcolors
pt.labelfont=16
fig = pt.gcf()
fig.subplots_adjust(left =0.14)
fig.subplots_adjust(right =0.98)
fig.subplots_adjust(bottom=0.12)
fig.subplots_adjust(top =0.98)
fig.subplots_adjust(wspace=0.0)
fig.subplots_adjust(hspace=0.0)
# read files
for file in files:
#nb = Nbody(file,ftype=opt.ftype)
nb = file
# apply options
nb = pt.do_reduc_options(nb,opt)
nb = pt.do_select_options(nb,opt)
nb = pt.do_center_options(nb,opt)
nb = pt.do_cmd_options(nb,opt)
nb = pt.do_info_options(nb,opt)
nb = pt.do_display_options(nb,opt)
################
# get values
################
x = nb.Fe()
y = nb.Mg()-nb.Fe()
db = 0.1
bins = arange(-4,0,db)
#n,bins,patches = pt.hist(Fe,bins,normed=True,fc=pt.pcolors[file],lw=0.5)
#c = x > -3
#x = compress(c,x)
h = libutil.myhistogram(x,bins)
h = h.astype(float)/sum(h)
pt.bar(bins,h,width=0.9*(bins[1]-bins[0]),color='k',alpha=0.3)
x = bins
y = h
#############################################################
# new method to compute histogram
#############################################################
# ok, ne change quasiment rien, car n_stars nearly identical for all star particles
stellar_mass = 9.8e-08
N_i = stellar_mass*2.846012674 * 1.0e10
m_tau = (nb.mass/stellar_mass*(100.0**(-0.35) - 0.1**(-0.35)) + 0.1**(-0.35))**(-1.0/0.35)
n_stars = N_i*(m_tau**(-1.35) - 0.1**(-1.35))/(100.0**(-1.35) - 0.1**(-1.35))
'''
# normal histogram
n = len(Fe)
shape = (len(bins),)
mn = min(bins)
mx = max(bins)
mass = ones(n)
val = ones(n)
# scale between 0 and 1
Fes = (Fe-mn)/(mx-mn)
h = mapping.mkmap1d(Fes.astype(float32),mass.astype(float32),val.astype(float32),shape)
h = h/sum(h)
pt.plot(bins,h,'b')
'''
# weighted histogram
n = len(x)
shape = (len(bins),)
mn = min(bins)
mx = max(bins)
mass = n_stars
unity = ones(n)
# scale between 0 and 1
Fes = (x-mn)/(mx-mn)
h = mapping.mkmap1d(Fes.astype(float32),mass.astype(float32),unity.astype(float32),shape)
h = h/sum(h)
#pt.plot(bins,h,'g')
#############################################################
if file == files[0]:
xmin,xmax,ymin,ymax = pt.SetLimits(opt.xmin,opt.xmax,opt.ymin,opt.ymax,x,y,opt.log)
if file==files[0]:
add_data_FeHistogram(opt.data)
#if opt.data!=None:
# pt.legend([r'$\rm{%s}$'%opt.data],'upper right',shadow=True)
# plot axis
pt.SetAxis(xmin,xmax,ymin,ymax,opt.log)
pt.xlabel(r'$\left[ \rm{Fe}/\rm{H} \right]$',fontsize=pt.labelfont)
pt.ylabel(r'$\rm{Stellar\,\,Fraction}$',fontsize=pt.labelfont)
pt.grid(False)
pt.EndPlot(files,opt)
########################################################################
# MAIN
########################################################################
#if __name__ == '__main__':
# files,opt = parse_options()
# pt.InitPlot(files,opt)
# pt.pcolors
# MakePlot(files,opt)
# pt.EndPlot(files,opt)
c4science · Help