Homec4science
Diffusion pNbody (master) Browse

pNbody/examples/plots/rotation_curve/plot_rotation_curve.py
179d38abe82dmaster

Last Change
Raw File

plot_rotation_curve.py
View Options

#!/usr/bin/env python
'''
Plot value computed in a spherical grid, as a function of the radius
Yves Revaz
Wed Jan 9 09:29:46 CET 2008
'''
import Ptools as pt
from pNbody import *
from pNbody import myNumeric
from pNbody import libdisk
from pNbody import libgrid
from pNbody import profiles
from pNbody import cosmo
import string
import sys
import os
from numpy import *
from optparse import OptionParser
from scipy.optimize import leastsq
from scipy.optimize import bisection
#from mpfit 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 = pt.add_legend_options(parser)
parser.add_option("--Rmax",
action="store",
dest="Rmax",
type="float",
default = 50.,
help="max radius of bins",
metavar=" FLOAT")
parser.add_option("--nR",
action="store",
dest="nR",
type="int",
default = 32,
help="number of bins in r",
metavar=" INT")
parser.add_option("--nt",
action="store",
dest="nt",
type="int",
default = 32,
help="number of bins in r",
metavar=" INT")
parser.add_option("--ErrTolTheta",
action="store",
dest="ErrTolTheta",
type="float",
default = 0.5,
help="Error tolerance theta",
metavar=" FLOAT")
parser.add_option("--fct",
action="store",
dest="fct",
type="string",
default = None,
help="transformation function ex : 'r**3' should be of the form M(r)",
metavar=" STR")
parser.add_option("--fctm",
action="store",
dest="fctm",
type="string",
default = None,
help="inverse transformation function ex : 'r**(1/3.)'",
metavar=" STR")
parser.add_option("--AdaptativeSoftenning",
action="store_true",
dest="AdaptativeSoftenning",
default = False,
help="AdaptativeSoftenning")
parser.add_option("--eps",
action="store",
dest="eps",
type="float",
default = 0.28,
help="smoothing length",
metavar=" FLOAT")
(options, args) = parser.parse_args()
if len(args) == 0:
print "you must specify a filename"
sys.exit(0)
files = args
return files,options
#######################################
# units
#######################################
# gear system of units
params = {}
params['UnitVelocity_in_cm_per_s'] = 20725573.785998672
params['UnitMass_in_g'] = 1.989e+43
params['UnitLength_in_cm'] = 3.085e+21
gear_units = units.Set_SystemUnits_From_Params(params)
out_units = units.UnitSystem('local',[units.Unit_km,units.Unit_Msol,units.Unit_s,units.Unit_K])
tokms = gear_units.convertionFactorTo(out_units.UnitVelocity)
#######################################
# MakePlot
#######################################
def MakePlot(files,opt):
# some inits
colors = pt.Colors(n=len(files))
datas = []
velocity_factor = 207.
# read files
for i,file in enumerate(files):
nb = Nbody(file,ftype=opt.ftype)
# 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)
#######################
# set the grid scaling
#######################
if opt.fct!=None:
print "g =lambda r:(%s)"%opt.fct
exec("g =lambda r:(%s)"%opt.fct)
else:
g = None
if opt.fctm!=None:
print "gm =lambda r:(%s)"%opt.fctm
exec("gm =lambda r:(%s)"%opt.fctm)
else:
gm = None
##############################
# compute the tree
##############################
nb.getTree(force_computation=True,ErrTolTheta=opt.ErrTolTheta)
##############################
# create cylindrical rt grid
##############################
Grt = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.Rmax,nr=opt.nR,nt=opt.nt,z=0,g=g,gm=gm)
# radial acceleration (selection)
Accx,Accy,Accz = Grt.get_AccelerationMap(nb,eps=opt.eps,UseTree=True,AdaptativeSoftenning=opt.AdaptativeSoftenning)
Ar = sqrt(Accx**2+Accy**2)
Ar = sum(Ar,1)/opt.nt
dPhi = Ar
# radius
R,t = Grt.get_rt()
# velocity
Vc = libdisk.Vcirc(R,dPhi) * tokms
datas.append(pt.DataPoints(R, Vc, color=colors.next(),label=i))
# plot
for d in datas:
pt.plot(d.x,d.y,color=d.color)
# set limits and draw axis
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
xlabel = r"$\rm{Radius}\,[\rm{kpc}]$"
ylabel = r"$\rm{Velocity}\,[\rm{km/s}]$"
# plot axis
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
pt.grid(False)
if opt.legend:
pt.LegendFromDataPoints(datas)
########################################################################
# MAIN
########################################################################
if __name__ == '__main__':
files,opt = parse_options()
pt.InitPlot(files,opt)
pt.pcolors
MakePlot(files,opt)
pt.EndPlot(files,opt)
c4science · Help