#coding:sjis
'''
Created on 2011/03/18
@author: uchiyama
'''
"""
importの使い方・考え方は以下が分かりやすい。
http://1-byte.jp/2010/09/29/learning_python5/
そもそもこのシリーズは図解が多くて良さげなので後で読む。
"""
from math import cos,sin,pi, atan, exp, sqrt
from scipy import integrate
"""
普遍であるべき観測量
ここにグローバルな変数として書くのは正しいのだろうか。
"""
#postion of Sgr A* in the Galactic coordinate (degree)
L_GC=-0.056
B_GC=-0.046
# 1 pc in the unit of cm
#CM_PER_PC=3.09e18
#Position of the Earth from the plane (pc)
Z0=16.0
#Distance to Sgr A* (pc)
R0=8.5e3
#The rotation angle
BD1=15.0
BD2=1.0
"""
座標変換に必要な関数。
本当はscipyの座標変換使えばもっと賢くできる気がする。
"""
def sin_d(degree):
return sin(2*pi*degree/360.0);
def cos_d(degree):
return cos(2*pi*degree/360.0);
def dl(l):
return l-L_GC
def db_bulge(b):
return b-B_GC-atan(Z0/R0)
def db(b):
return b-B_GC
def r_f(s,l,b):
l0=dl(l)
b0=db(b)
(R0**2+(s*cos_d(b0))**2-2*R0*s*cos_d(l0)*cos_d(b0))**0.5
def z_f(s, l, b):
b0=db(b)
return s*sin_d(b0)
def x_f_bulge(s, l, b):
l0=dl(l)
b0=db(b)
return R0-s*cos_d(b0)*cos_d(l0)
def y_f_bulge(s, l, b):
l0=dl(l)
b0=db_bulge(b)
s*cos_d(b0)*sin_d(l0)
def z_f_bulge(s, l, b):
b0=db_bulge(b)
return s*sin_d(b0)+Z0
def X_f(x, y, z):
return x*cos_d(BD1)+y*sin_d(BD1)
def Y_f(x, y, z):
return -1*x*sin_d(BD1)*cos_d(BD2)+y*cos_d(BD1)*cos_d(BD2)+z*sin_d(BD2);
def Z_f(x, y, z):
return -1*x*sin_d(BD1)*sin_d(BD2)+y*cos_d(BD1)*sin_d(BD2)+z*cos_d(BD2);
"""
以下星質量モデル。
基本的にMuno+06に基づく…が数式に誤りたくさん+説明不足。
http://adsabs.harvard.edu/abs/2006ApJS..165..173M
Kent+91, Launhardt+02を元に情報を補足、整理。
http://adsabs.harvard.edu/abs/1991ApJ...378..131K
http://adsabs.harvard.edu/abs/2002A%26A...384..112L
更に詳しくは
http://www-utheal.phys.s.u-tokyo.ac.jp/~uchiyama/wordpress/wp-content/uploads/2011/05/SMM.pdf
"""
def NSC(x,l,b):
norm1=3.3e6 # M_solar pc^{-3}
norm2=norm1*8.98e7/3.3e6
n1=2.0
n2=3.0
R_c=0.22 # pc
R_lim=6.0
R_cutoff=200.0
r=r_f(x,l,b)
z=z_f(x,l,b)
"""
powや**は（defaultで）使えるのにsqrtはmathをimportしないと使えない。
意味がよく分からないが、そもそもsqrt不要と言われるとそれが正しい気もする。
"""
R=sqrt(pow(r,2)+pow(z,2))
"""
元々は3項演算子を使っていたがpythonの3項演算子に慣れずネストの仕方が分からないのでif文で書く。
"""
if R < R_lim :
return norm1/(1+pow((R/R_c),n1))
elif R < R_cutoff :
return norm2/(1+pow((R/R_c),n2))
else:
return 0
def NSD(x,l,b):
n1=-0.1
n2=-3.5
n3=-10.0
r_1=120.0
r_2=220.0
r_3=2000.0
r_d=1.0
z_d=45.0
norm1=300.0
norm2=norm1*(r_1/r_d)**(-n2+n1)
norm3=norm2*(r_2/r_d)**(-n3+n2)
r=r_f(x,l,b)+1e-30
z=z_f(x,l,b)
if r <=0 :
return 0.0
elif r < r_1 :
return norm1*(abs(r/r_d)**n1)*exp(-1.0*z/z_d)
elif r < r_2 :
return norm2*(abs(r/r_d)**n2)*exp(-1.0*z/z_d)
elif r < r_3 :
return norm3*(abs(r/r_d)**n3)*exp(-1.0*z/z_d)
else:
return 0.0
def GB(x,l,b):
"""
Unit of value returned by this function
M_solar pc^{-3}
"""
#norm Unit: M_solar pc^{-3}
norm=8.0
#a_x,a_y,a_z : pc
a_x=1100.0
a_y=360.0
a_z=220.0
c_para=3.2
c_perp=1.6
x_b=x_f_bulge(x, l, b);
y=y_f_bulge(x, l, b);
z=z_f_bulge(x, l, b);
"""
Pythonはべき乘を表す演算子**もあるけどpowもある。
そして小数乘もできるし、負の数を小数乗して複素数を
返すこともできる。
"""
R_perp=((abs(X_f(x_b,y,z))/a_x)**c_perp+(abs(Y_f(x_b,y,z))/a_y)**c_perp)**(1/c_perp)
Rs=(pow(abs(Z_f(x_b,y,z))/a_z,c_para)+pow(R_perp,c_para))**(1/c_para);
return norm*exp(-1*Rs);
def GD(x,l,b):
"""
Unit of this function
M_solar pc^{-3}
"""
#norm Unit: M_solar pc^{-3}
norm=5.0
#r_d, z_d : pc
r_d=2700.0;
z_d=200.0;
r=r_f(x, l, b)
z=z_f(x, l, b)
return norm*exp(-1.0*r/r_d)*exp(-1.0*(z/z_d))
def Stellar_density_f(x, l, b):
"""
Unit of this function
M_solar pc^{-3}
"""
return NSC(x,l,b)+NSD(x,l,b)+GB(x,l,b)+GD(x,l,b)
def SMM(l, b, s=16e3 , region='SMM'):
"""
Pythonで関数はオブジェクトなので下記の様なことができてしまう。
正直キモイ。ポインタを渡したくなる…。以下を参照。
http://python.g.hatena.ne.jp/muscovyduck/20080707/p2
"""
if region=='NSC':
func=NSC
elif region=='NSD':
func=NSD
elif region=='GD':
func=GD
elif region=='GB':
func=GB
elif region=='SMM':
func=SMM
else :
exit("3rd argument must be one of NSC/NSD/GD/GB/SMM.")
"""
SciPy(さいぱい。すしぱいと読んでいた…)を使った積分。
http://docs.scipy.org/doc/scipy/scipy-ref.pdf
の1.4.1を見る。
lambda式の意図は
http://atkonn.blogspot.com/2008/02/python-python27-lambda.html
をみるとなんとなく分かる。
"""
res=integrate.quad(lambda x: func(x,l,b), 0.0, s)[0]
"""
The unit of returned value is M_solar/pc^2/str
"""
return res/(4.0*pi);
