Page MenuHomec4science
Files F73636581

plot_cylindrical_model_miyamoto.py
No OneTemporary
Actions

File Metadata

Created
Tue, Jul 23, 08:24

plot_cylindrical_model_miyamoto.py
View Options

#!/usr/bin/env python
import sys
import Ptools as pt
from numpy import *
from pNbody import libmiyamoto
stats = pt.io.read_dmp(sys.argv[1])
n = stats['n']
M = stats['M']
G = stats['G']
hR = stats['hR']
hz = stats['hz']
a = stats['a']
b = stats['b']
Rmax = stats['Rmax']
Zmax = stats['Zmax']
eps = stats['eps']
nz = stats['nz']
nr = stats['nr']
rmin = stats['rmin']
rmax = stats['rmax']
zmin = stats['zmin']
zmax = stats['zmax']
###################
# plot
###################
r = stats['R']
nn = sum(stats['nn'],axis=1)
rho = stats['rho']
phi = stats['phi']
sigma_z = stats['sigma_z']
z = zeros(len(r))
rho = rho[:,nz/2] # valeur dans le plan
phi = phi[:,nz/2] # valeur dans le plan
sigma_z = sigma_z[:,nz/2]
# theorical values
M = 1.1
rho_th = libmiyamoto.Density(1,M,a,b,r,z)
phi_th = libmiyamoto.Potential(1,M,a,b,r,z)
sigma_z_th = libmiyamoto.Sigma_z(1,M,a,b,r,z)
#############
# number per bins
#############
pt.subplot(2,2,1)
pt.plot(r,nn)
pt.semilogx()
pt.semilogy()
pt.xlabel('Radius')
pt.ylabel('Number')
#############
# density
#############
pt.subplot(2,2,2)
pt.plot(r,rho,'b')
pt.plot(r,rho_th,'r--')
pt.semilogx()
#pt.semilogy()
pt.xlabel('Radius')
pt.ylabel('Density')
#############
# potential
#############
pt.subplot(2,2,3)
pt.plot(r,phi,'b')
pt.plot(r,phi_th,'r--')
pt.semilogx()
#pt.semilogy()
pt.xlabel('Radius')
pt.ylabel('Potential')
#############
# sigma
#############
pt.subplot(2,2,4)
pt.plot(r,sigma_z,'b')
pt.plot(r,sigma_z_th,'r--')
pt.semilogx()
#pt.semilogy()
pt.xlabel('Radius')
pt.ylabel('Sigma')
####################
# velocity curves
####################
fig = pt.figure()
pt.subplot(1,1,1)
pt.plot(r,stats['vc'],'k')
pt.plot(r,stats['vm'],'y')
pt.plot(r,stats['sr'],'r')
pt.plot(r,stats['sp'],'g')
pt.plot(r,stats['sz'],'b')
pt.xlabel('Radius')
pt.ylabel('Velocity')
pt.title('Rotation curve and velocity dispertions')
pt.legend(('vc','vm','n'))
####################
# freq.
####################
fig = pt.figure()
pt.subplot(1,1,1)
pt.plot(r,stats['kappa'],'r')
pt.plot(r,stats['omega'],'g')
#pt.plot(r,stats['nu'],'b')
pt.xlabel('Radius')
pt.ylabel('Frequencies')
pt.title('Frequences')
pt.legend(('kappa','omega','nu'))
pt.show()

Event Timeline

c4science · Help