<jupyter_start><jupyter_text>Mass model of the Galaxy<jupyter_code>%matplotlib inline
from matplotlib import rcParams
rcParams["savefig.dpi"] = 100
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gammainc<jupyter_output><empty_output><jupyter_text>Observational constraints[Piffl et al. 2014](http://adsabs.harvard.edu/abs/2014A%26A...562A..91P) Mvir = (1.2-1.3)e12 Msun[Deason et al. 2012](http://adsabs.harvard.edu/abs/2012MNRAS.425.2840D) [at 150 kpc] M = (5-10)e11 Msun[Gnedin et al. 2010](http://adsabs.harvard.edu/abs/2010ApJ...720L.108G) [at 80 kpc] M = (6.9 +3.0-1.2)e11 Msun[Deason et al. 2012](http://adsabs.harvard.edu/abs/2012MNRAS.424L..44D) [at 50 kpc] M = (4.2 ± 0.4)e11 Msun[McMillan 2016](http://arxiv.org/abs/1608.00971) M_disk,stellar = 5.4e10, Mvir = 1.3e12 Msun[Bovy & Rix 2013](http://adsabs.harvard.edu/abs/2013ApJ...779..115B) [at 4-9 kpc] M_disk = 5.3e10, M_disk,stellar = 4.6e10[Nakanishi & Sofue](http://adsabs.harvard.edu/abs/2016PASJ...68....5N) M_gas = 8e9 Msun[Bland-Hawthorn & Gerhard 2016](http://arxiv.org/abs/1602.07702) M_NSD = (1.4 ± 0.6)e9 Msun, r_h,NSD = 90 pc,M_bulge = (1.4-1.7)e10 Msun, M_thin = (4 ± 1)e10, r_thin = 2.6 ± 0.5 kpc, M_thick = (8 ± 3)e9 Msun, r_thick = 2 ± 0.2 kpc, R_0 = 8.2 ± 0.1 kpc, V_0 = 238 ± 15 km/s, V_0/R_0 = 29.0 ± 1.8 km/s/kpc, M(8.2kpc) = 1.08e11 Msun[Launhardt et all. 2002](http://adsabs.harvard.edu/abs/2002A%26A...384..112L) [at 120 pc] NSD mass = (8 ± 2)e8 Msun, scale radius = 90 pc[Feldmeier et al. 2014](http://adsabs.harvard.edu/abs/2014A%26A...570A...2F) [at 10 pc] total mass = (3 ± 1)e7 Msun[Chatzopoulos et al. 2015](http://adsabs.harvard.edu/abs/2015MNRAS.447..948C) [at 1 and 4 pc] R_0 = 8.33 ± 0.11 kpc, M_BH = (4.23 ± 0.14)e6 Msun, M_NSC = (1.8 ± 0.3)e7 Msun, r_h,NSC = 4.2 ± 0.4 pc, M_NSC(1pc) = 0.89e6 Msun<jupyter_code>robs =   np.array([ 0.001, 0.004, 0.01,  0.12, 8.2,     50.,    80.,    150.  ])
Mobs =   np.array([ 5.1e6, 1.3e7, 2.6e7, 8.e8, 1.08e11, 4.2e11, 6.9e11, 9.0e11 ])
Mobs_l = np.array([ 4.6e6, 1.1e7, 1.6e7, 6.e8, 9.37e10, 3.8e11, 5.0e11, 5.0e11 ])
Mobs_u = np.array([ 5.6e6, 1.5e7, 3.6e7, 1.e9, 1.24e11, 4.6e11, 9.9e11, 1.1e12 ])<jupyter_output><empty_output><jupyter_text>Nuclear star cluster mass distribution from [Chatzopoulos et al. 2015](http://adsabs.harvard.edu/abs/2015MNRAS.447..948C)<jupyter_code>def M_dehnen( x, gam ):
    return np.power(x/(1.+x), 3.-gam)

def Mass_NSC( r ):
    mfrac1 = 1./106.45       # fraction of mass in first component
    mfrac = [ mfrac1, 1.-mfrac1 ]
    rh = 0.0042              # half-mass radius of the nuclear star cluster in kpc
    gam = [ 0.51, 0.07 ]     # inner logarithmic slope
    ascale = [ 99., 2376. ]  # scale length in arcsec
    arcsec = 4.e-5           # 1 arcsec in kpc at the distance of the Galactic Center
    asc = np.array(ascale)*arcsec
    part = [ frac*M_dehnen(r/a, g) for (a,g,frac) in zip(asc,gam,mfrac) ]
    parth = [ frac*M_dehnen(rh/a, g) for (a,g,frac) in zip(asc,gam,mfrac) ]
    fracm = np.minimum( np.sum(part)/np.sum(parth)/2., 1. )
    return Mnsc*fracm<jupyter_output><empty_output><jupyter_text>Galactic mass components: nuclear star cluster, bulge, disk, and dark matter halo<jupyter_code>def NSC():
    M = 1.8e7   # mass of the nuclear star cluster in Msun
    return M

def Bulge():
    M = 1.4e10  # mass of stellar bulge/bar in Msun (in G05 was 1e10)
    a = 0.4     # scale length of stellar bulge in kpc (in G05 was 0.6)
    return M, a

def Disk():
    M = 5.6e10  # mass of stellar and gaseous disk in Msun (in G05 was 4e10)
    a = 2.6     # scale length of stellar disk in kpc (in G05 was 5)
    b = 0.3     # scale height of stellar disk in kpc
    return M, a, b

def Halo():
    M = 1.2e12  # mass of dark matter halo
    rs = 20.    # halo scale radius, in kpc
    xm = 2.2    # scaled radius of maximum circular velocity
    return M, rs, xm

def SMBH():
    M = 4.2e6   # mass of central black hole
    return M

Mnsc = NSC()
Mbulge, abulge = Bulge()
Mdisk, adisk, bdisk = Disk()
Mhalo, rs, xm = Halo()
MBH = SMBH()

Mvir = Mhalo + Mdisk + Mbulge + Mnsc + MBH

kms2 = 4.30e-6                         # conversion from GMsun/kpc to (km/s)^2
Rvir = 56.*np.power(Mvir/1.e10, 1./3.) # virial radius in kpc, for delta0=340
c = Rvir/rs                            # halo concentration parameter
Mh = Mhalo/(np.log(1.+c)-c/(1.+c))
print('M_vir = %.2e Msun  R_vir = %.1f kpc  c_vir = %.1f'%(Mvir, Rvir, c))

#print 'M_NSC = %.2e Msun'%(Mass_NSC(0.01))
# from galaxy_mass_model import galaxy_mass_model

# gm = galaxy_mass_model()

# MBH = gm.M_BH
# print(gm.M_BH)
plt.figure(figsize=(8,6))
plt.xlim(-3.5, 2.5)
plt.ylim(6.4, 12.4)
rcParams['lines.linewidth'] = 1.5
rcParams['xtick.major.size'] = 6
rcParams['ytick.major.size'] = 6
rcParams['xtick.labelsize'] = 14
rcParams['ytick.labelsize'] = 14
plt.xlabel(r'$\log{\,r}\ (\mathrm{kpc})$', fontsize=18)
plt.ylabel(r'$\log{\,M}\ (M_\odot)$', fontsize=18)

lgr = np.arange(-3.5, 2.6, 0.05)
r = 10.**lgr

# best model
Mnsc_g = np.array([ Mass_NSC(rr) for rr in r ])
#Mbulge_g = Mbulge*r**2/(r + abulge)**2
#Mdisk_g = Mdisk*r**3/(r**2 + (adisk+np.sqrt(0.**2+bdisk**2))**2)**1.5
#Mhalo_g = Mh*(np.log(1.+r/rs) - r/rs/(1.+r/rs))
Mbulge_g = np.array([ Mass_Bulge(rr) for rr in r ])
Mdisk_g = np.array([ Mass_Disk(rr, 0.) for rr in r ])
Mhalo_g = np.array([ Mass_Halo(rr) for rr in r ])
Mtot = MBH + Mnsc_g + Mbulge_g + Mdisk_g + Mhalo_g

#plt.plot(lgr, np.log10(Mnsc_g), 'k--')
#plt.plot(lgr, np.log10(Mbulge_g), 'k:')
#plt.plot(lgr, np.log10(Mdisk_g), 'k-.')
plt.plot(lgr, np.log10(Mtot), 'k-')

#plt.text(1.2, 7.65, 'nuclear cluster', fontsize=12)
#plt.text(1.9, 9.85, 'bulge', fontsize=12)
#plt.text(1.9, 10.45, 'disk', fontsize=12)
#plt.text(1.9, 11.4, 'halo', fontsize=12)

# Sersic fit, used in Gnedin, Ostriker & Tremaine 2014
nser = 4.     # Sersic index (in G14 was 2.2)
aser = 4.     # effective radius, in kpc
bn = 2.*nser-1./3.+0.0098765/nser+0.0018/nser**2
Mser = 5.e10*gammainc(2*nser, bn*(r/aser)**(1./nser))
#plt.plot(lgr, np.log10(Mser + Mhalo_g + MBH), 'g-')
    
# Gnedin 2005 model
Mbulge_g5 = 1e10*r**2/(r + 0.6)**2
Mdisk_g5 = 4e10*r**3/(r**2 + (5.+0.3)**2)**1.5
Mhalo_g5 = Mh/1.2*(np.log(1.+r/rs) - r/rs/(1.+r/rs))
Mtot_g5 = MBH + Mbulge_g5 + Mdisk_g5 + Mhalo_g5
#plt.plot(lgr, np.log10(Mtot_g5), 'g-')
#plt.text(1., 7.7, 'Gnedin+05', color='g', fontsize=12)

# Kenyon 2008 model, updated in Brown et al. 2014
Mbulge_k = 3.76e9*r**2/(r + 0.1)**2
Mdisk_k = 6e10*r**3/(r**2 + (2.75+bdisk)**2)**1.5
Mtot_k = MBH + Mbulge_k + Mdisk_k + Mhalo_g/1.2
#plt.plot(lgr, np.log10(Mtot_k), 'b-')
#plt.text(1., 8.3, 'Kenyon+08', color='b', fontsize=12)

# observational points
plt.scatter(np.log10(robs), np.log10(Mobs), s=20, marker='s', color='k')

yerr1 = np.log10(Mobs) - np.log10(Mobs_l)
yerr2 = np.log10(Mobs_u) - np.log10(Mobs)
plt.errorbar(np.log10(robs), np.log10(Mobs), yerr=[yerr1,yerr2], ecolor='k', capthick=0, linestyle='None')

plt.show()
#plt.savefig('galactic_mass_compare.png')<jupyter_output><empty_output><jupyter_text>Escape velocity curve<jupyter_code># best model
pot = -Mbulge/(r+abulge) -Mnsc/r -MBH/r -Mdisk/np.sqrt(0**2+(adisk+np.sqrt(r**2+bdisk**2))**2) -Mh/r*np.log(1.+r/rs)
Vesc = np.sqrt(-2.*pot*kms2)

# Kenyon 2008 model
pot_k = -3.76e9/(r+0.1) -MBH/r -6e10/np.sqrt(0**2+(2.75+np.sqrt(r**2+bdisk**2))**2) -Mh/1.2/r*np.log(1.+r/rs)
Vesc_k = np.sqrt(-2.*pot_k*kms2)

plt.figure(figsize=(6,4))
plt.xlim(-3, 2.5)
plt.ylim(0, 1000)
rcParams['lines.linewidth'] = 1.5
rcParams['xtick.major.size'] = 6
rcParams['ytick.major.size'] = 6
rcParams['xtick.labelsize'] = 12
rcParams['ytick.labelsize'] = 12
plt.xlabel(r'$\log{\,r}\ (\mathrm{kpc})$', fontsize=18)
plt.ylabel(r'$V_{esc}\ (\mathrm{km\, s}^{-1})$', fontsize=18)

plt.plot(lgr, Vesc, 'k-')
plt.plot(lgr, Vesc_k, 'b-')

plt.show()
plt.figure(figsize=(6,4))
plt.xlim(0, 100)
plt.ylim(200, 800)
plt.xlabel(r'$r\ (\mathrm{kpc})$', fontsize=18)
plt.ylabel(r'$V_{esc}\ (\mathrm{km\, s}^{-1})$', fontsize=18)

plt.plot(r, Vesc, 'k-')
plt.plot(r, Vesc_k, 'b-')

plt.show()
for lev in [ -3, -2, -1, 0., 1., 2. ]:
    l = np.fabs(lgr-lev) < 0.001
    print r[l], Vesc[l], Vesc_k[l]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20250061289531748, 'fraction_numerical': 0.10149546457465065, 'mean_word_length': 3.550474065811489, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '507527', 'license': 'MIT', 'n_tokens_mistral': 4398, 'n_tokens_neox': 3831, 'n_words': 812, 'path': 'notebooks/oleg_galactic_model.ipynb', 'repo_name': 'adrn/uncluster', 'size': '8158'}