These Python scripts can be used to reproduce the figures in Bovy & Tremaine (2012): On the local dark matter density.
Figure 1 is produced by:
python plot_dvcapprox.pywhile Figure 2 is produced by:
python plot_sigmaz.pyThis code requires galpy, monoAbundanceMW, as well as the Numpy, Scipy, Matplotlib stack.
| #####Routines related to the dvc/dr approximation | |
| import os, os.path | |
| import cPickle as pickle | |
| import numpy | |
| from scipy import integrate | |
| from galpy.potential import DoubleExponentialDiskPotential, \ | |
| LogarithmicHaloPotential, NFWPotential, HernquistPotential, \ | |
| MiyamotoNagaiPotential | |
| from galpy.util import bovy_plot | |
| _R0= 2.2*3.5 #8. #kpc | |
| def a2b2approx(pot,dR=0.00001,z=0.): | |
| """ | |
| NAME: | |
| a2b2approx | |
| PURPOSE: | |
| approximate the A2-B2 by its value in the plane | |
| INPUT: | |
| pot - galpy potential | |
| OUTPUT: | |
| value | |
| HISTORY: | |
| 2012-04-30 - Written - Bovy (IAS) | |
| """ | |
| #Approximation: FR/R+dFR/dR | |
| zeroRforce= pot.Rforce(1.,z) | |
| return (zeroRforce+(pot.Rforce(1.+dR,z)-zeroRforce)/dR)/2. | |
| #return (zeroRforce-pot.R2deriv(1.,0))/2. | |
| def plot_dvcapprox(): | |
| nzs= 21 | |
| zs= numpy.linspace(1.,4.,nzs)/_R0 | |
| #First calculate Double-Exp disk with hR=3., hz=300. | |
| dp= DoubleExponentialDiskPotential(normalize=1.,hr=3.3/_R0,hz=.437/_R0,tol=0.001) | |
| kz= numpy.zeros_like(zs) | |
| for ii in range(nzs): kz[ii]= dp.zforce(1.,zs[ii]) | |
| approx= -kz/2./numpy.pi-a2b2approx(dp)/numpy.pi*zs | |
| approx2= -(1./3.3*_R0-1./(3.3/_R0)**2.)*(zs**2./2.-0.437/_R0*zs+(0.437/_R0)**2.-(0.437/_R0)**2.*numpy.exp(-zs/0.437*_R0)) | |
| sz= dp.dens(1.,0.)*2.*.437/_R0*(1.-numpy.exp(-zs/0.437*_R0)) | |
| print a2b2approx(dp), zs, (sz-approx)/sz, (sz-approx2)/sz | |
| bovy_plot.bovy_print() | |
| bovy_plot.bovy_plot(zs*_R0,(sz-approx)/sz,'r-',loglog=False, | |
| xrange=[0.1,5.], | |
| yrange=[0.000,.3], | |
| xlabel=r'$Z\ [\mathrm{kpc}]$', | |
| ylabel=r'$(\Sigma^{\mathrm{true}} - \Sigma^{\mathrm{approx}}) / \Sigma^{\mathrm{true}}$') | |
| bovy_plot.bovy_plot(zs*_R0,approx2,'r--',overplot=True,lw=1.) | |
| bovy_plot.bovy_plot(zs*_R0,(sz-approx/(1.-approx2))/sz, | |
| 'r--',overplot=True,lw=2.) | |
| bovy_plot.bovy_plot([0.,5.],[0.,0.],'-',color='0.65',overplot=True) | |
| #Another de disk, hR= 4.5,hz=300., KG | |
| dp= DoubleExponentialDiskPotential(normalize=1.,hr=4.5/_R0,hz=.4/_R0,tol=0.001) | |
| for ii in range(nzs): kz[ii]= dp.zforce(1.,zs[ii]) | |
| approx= -kz/2./numpy.pi-a2b2approx(dp)/numpy.pi*zs | |
| approx2= approx/(1.+(1./4.5*_R0-1./(4.5/_R0)**2.)*zs**2./4.) | |
| sz= dp.dens(1.,0.)*2.*.4/_R0*(1.-numpy.exp(-zs/0.4*_R0)) | |
| print "Here", a2b2approx(dp), zs, (sz-approx)/sz, (sz-approx2)/sz | |
| bovy_plot.bovy_plot(zs*_R0,(sz-approx)/sz,'c-',overplot=True) | |
| bovy_plot.bovy_plot(zs*_R0,(sz-approx2)/sz,'c-',overplot=True,lw=2.) | |
| #Logarithmic halo | |
| dp= LogarithmicHaloPotential(normalize=1.) | |
| for ii in range(nzs): kz[ii]= dp.zforce(1.,zs[ii]) | |
| approx= -kz/2./numpy.pi-a2b2approx(dp)/numpy.pi*zs | |
| print a2b2approx(dp) | |
| for ii in range(nzs): sz[ii]= 2.*integrate.quad((lambda x: dp.dens(1.,x)),0.,zs[ii])[0] | |
| bovy_plot.bovy_plot(zs*_R0,(sz-approx)/sz,'r-',overplot=True,lw=2.) | |
| #Logarithmic halo, oblate | |
| dp= LogarithmicHaloPotential(normalize=1.,q=0.8) | |
| for ii in range(nzs): kz[ii]= dp.zforce(1.,zs[ii]) | |
| approx= -kz/2./numpy.pi-a2b2approx(dp)/numpy.pi*zs | |
| print a2b2approx(dp) | |
| for ii in range(nzs): sz[ii]= 2.*integrate.quad((lambda x: dp.dens(1.,x)),0.,zs[ii])[0] | |
| bovy_plot.bovy_plot(zs*_R0,(sz-approx)/sz,'g-',overplot=True) | |
| #Logarithmic halo, prolate (?) | |
| dp= LogarithmicHaloPotential(normalize=1.,q=1.2) | |
| for ii in range(nzs): kz[ii]= dp.zforce(1.,zs[ii]) | |
| approx= -kz/2./numpy.pi-a2b2approx(dp)/numpy.pi*zs | |
| for ii in range(nzs): sz[ii]= 2.*integrate.quad((lambda x: dp.dens(1.,x)),0.,zs[ii])[0] | |
| bovy_plot.bovy_plot(zs*_R0,(sz-approx)/sz,'y-',overplot=True) | |
| #NFW | |
| dp= NFWPotential(a=4.5,normalize=1.) | |
| for ii in range(nzs): kz[ii]= dp.zforce(1.,zs[ii]) | |
| approx= -kz/2./numpy.pi-a2b2approx(dp)/numpy.pi*zs | |
| approx2= approx/(1.+(1./6.*_R0-1./(6./_R0)**2.)*zs**2./4.) | |
| print a2b2approx(dp) | |
| for ii in range(nzs): sz[ii]= 2.*integrate.quad((lambda x: dp.dens(1.,x)),0.,zs[ii])[0] | |
| print sz, approx, approx2, (sz-approx2)/sz | |
| bovy_plot.bovy_plot(zs*_R0,(sz-approx)/sz,'y-',overplot=True) | |
| bovy_plot.bovy_plot(zs*_R0,(sz-approx2)/sz,'y-',overplot=True,lw=2.) | |
| #Miyamoto-Nagai | |
| dp= MiyamotoNagaiPotential(a=0.5,b=.4/_R0,normalize=1.) | |
| for ii in range(nzs): kz[ii]= dp.zforce(1.,zs[ii]) | |
| approx= -kz/2./numpy.pi-a2b2approx(dp)/numpy.pi*zs | |
| print a2b2approx(dp) | |
| #sz= dp.dens(1.,0.)*2.*.4/_R0*(1.-numpy.exp(-zs/0.4*_R0)) | |
| for ii in range(nzs): sz[ii]= 2.*integrate.quad((lambda x: dp.dens(1.,x)),0.,zs[ii])[0] | |
| print sz, approx | |
| bovy_plot.bovy_plot(zs*_R0,(sz-approx)/sz,'m-',overplot=True) | |
| #Double-exp disk + NFW | |
| #dp= DoubleExponentialDiskPotential(normalize=.65,hr=3.5/_R0,hz=.4/_R0) | |
| dp= MiyamotoNagaiPotential(a=0.45,b=.4/_R0,normalize=.65) | |
| np= NFWPotential(a=4.5,normalize=.35) | |
| for ii in range(nzs): kz[ii]= dp.zforce(1.,zs[ii])+np.zforce(1.,zs[ii]) | |
| approx= -kz/2./numpy.pi-a2b2approx(dp)/numpy.pi*zs-a2b2approx(np)/numpy.pi*zs | |
| print a2b2approx(dp),a2b2approx(np), a2b2approx(dp)+a2b2approx(np) | |
| #sz= dp.dens(1.,0.)*2.*.4/_R0*(1.-numpy.exp(-zs/0.4*_R0)) | |
| for ii in range(nzs): sz[ii]= 2.*integrate.quad((lambda x: dp.dens(1.,x)),0.,zs[ii])[0] | |
| for ii in range(nzs): sz[ii]+= 2.*integrate.quad((lambda x: np.dens(1.,x)),0.,zs[ii])[0] | |
| print sz, approx | |
| bovy_plot.bovy_plot(zs*_R0,(sz-approx)/sz,'b-',overplot=True) | |
| bovy_plot.bovy_end_print('test.png') | |
| if __name__ == '__main__': | |
| plot_dvcapprox() |
| import os, os.path | |
| import cPickle as pickle | |
| import numpy | |
| from scipy import integrate | |
| from galpy.potential import DoubleExponentialDiskPotential, \ | |
| LogarithmicHaloPotential, NFWPotential, HernquistPotential, \ | |
| MiyamotoNagaiPotential | |
| from galpy.util import bovy_plot | |
| from matplotlib import pyplot | |
| import monoAbundanceMW | |
| import dvcapprox | |
| _R0= 8. #kpc | |
| def plot_dvcapprox(): | |
| nzs= 21 | |
| zs= numpy.linspace(1.5,4.,nzs)/_R0 | |
| #First calculate Double-Exp disk with hR=3.4, hz=392 | |
| dp= DoubleExponentialDiskPotential(normalize=1., | |
| hr=3.4/_R0,hz=.392/_R0,tol=0.001) | |
| kz= numpy.zeros_like(zs) | |
| for ii in range(nzs): kz[ii]= dp.zforce(1.,zs[ii]) | |
| approx= -kz/2./numpy.pi-dvcapprox.a2b2approx(dp)/numpy.pi*zs | |
| sz= dp.dens(1.,0.)*2.*.392/_R0*(1.-numpy.exp(-zs/0.392*_R0)) | |
| bovy_plot.bovy_print(fig_height=3.) | |
| line1= bovy_plot.bovy_plot(zs*_R0,(sz-approx)/sz,'-.',color='k', | |
| xrange=[0.,5.], | |
| yrange=[0.000,.21], | |
| xlabel=r'$Z\ [\mathrm{kpc}]$', | |
| ylabel=r'$(\Sigma^{\mathrm{true}} - \Sigma^{\mathrm{approx}}) / \Sigma^{\mathrm{true}}$') | |
| x1= (sz-approx)/sz | |
| #NFW | |
| dp= NFWPotential(a=2.8,normalize=1.) #Scale radius 22.25 kpc from Xue et al. | |
| for ii in range(nzs): kz[ii]= dp.zforce(1.,zs[ii]) | |
| approx= -kz/2./numpy.pi-dvcapprox.a2b2approx(dp)/numpy.pi*zs | |
| for ii in range(nzs): sz[ii]= 2.*integrate.quad((lambda x: dp.dens(1.,x)),0.,zs[ii])[0] | |
| line2= bovy_plot.bovy_plot(zs*_R0,(sz-approx)/sz,'--',overplot=True,color='k') | |
| x2= (sz-approx)/sz | |
| #Save to pickle | |
| savefile= open('approx.sav','wb') | |
| pickle.dump(zs,savefile) | |
| pickle.dump(x1,savefile) | |
| pickle.dump(x2,savefile) | |
| savefile.close() | |
| #Combination | |
| dp= DoubleExponentialDiskPotential(normalize=.85, | |
| hr=3.4/_R0,hz=.392/_R0,tol=0.001) | |
| np= NFWPotential(a=2.8,normalize=.15) | |
| for ii in range(nzs): kz[ii]= dp.zforce(1.,zs[ii])+np.zforce(1.,zs[ii]) | |
| approx= -kz/2./numpy.pi-dvcapprox.a2b2approx(dp)/numpy.pi*zs\ | |
| -dvcapprox.a2b2approx(np)/numpy.pi*zs | |
| print -dvcapprox.a2b2approx(dp)-dvcapprox.a2b2approx(np) | |
| sz= dp.dens(1.,0.)*2.*.392/_R0*(1.-numpy.exp(-zs/0.392*_R0)) | |
| for ii in range(nzs): sz[ii]+= 2.*integrate.quad((lambda x: np.dens(1.,x)),0.,zs[ii])[0] | |
| line3= bovy_plot.bovy_plot(zs*_R0,(sz-approx)/sz,'-',overplot=True,color='k') | |
| #Legend | |
| pyplot.legend((line1[0],line2[0],line3[0]), | |
| (r'$\mathrm{exponential\ disk\ w/}\ h_{R} = 3.4\ \mathrm{kpc}, h_{Z} = 392\ \mathrm{pc}$', | |
| r'$\mathrm{NFW\ w/\ scale\ radius}\ R_c = 22.25\ \mathrm{kpc}$', | |
| r'$\mathrm{combination\ w/\ 85\%\ disk}$'+'\n'+r'$\mathrm{contribution\ to}\ F_R$'), | |
| loc='upper left',#bbox_to_anchor=(.91,.375), | |
| # numpoints=100, | |
| prop={'size':10}, | |
| frameon=False) | |
| bovy_plot.bovy_end_print('../tex/approx.ps') | |
| return None | |
| #Mix of disks | |
| fehs= monoAbundanceMW.fehs() | |
| afes= monoAbundanceMW.afes() | |
| norm= 0. | |
| for ii in range(len(fehs)): | |
| norm+= monoAbundanceMW.abundanceDist(fehs[ii],afes[ii]) | |
| disks= [] | |
| for ii in range(len(fehs)): | |
| print fehs[ii], afes[ii], monoAbundanceMW.hr(fehs[ii],afes[ii]),\ | |
| monoAbundanceMW.hz(fehs[ii],afes[ii]) | |
| if monoAbundanceMW.hr(fehs[ii],afes[ii]) == 0.: continue | |
| if monoAbundanceMW.hr(fehs[ii],afes[ii]) < 2.: thishr= 2./_R0 | |
| else: thishr= monoAbundanceMW.hr(fehs[ii],afes[ii])/_R0 | |
| disks.append(DoubleExponentialDiskPotential(normalize=1.*monoAbundanceMW.abundanceDist(fehs[ii],afes[ii])/norm, | |
| hr=thishr, | |
| hz=monoAbundanceMW.hz(fehs[ii],afes[ii])/_R0)) | |
| #Calculate | |
| kz= numpy.zeros_like(zs) | |
| approx= numpy.zeros_like(zs) | |
| for dp in disks: | |
| for ii in range(nzs): | |
| kz[ii]+= dp.zforce(1.,zs[ii]) | |
| approx[ii]+= -kz[ii]/2./numpy.pi-dvcapprox.a2b2approx(dp)/numpy.pi*zs | |
| sz= numpy.zeros_like(zs) | |
| for dp in disks: | |
| sz+= dp.dens(1.,0.)*2.*dp._hz*(1.-numpy.exp(-zs/dp._hz)) | |
| bovy_plot.bovy_plot(zs*_R0,(sz-approx)/sz,'--',overplot=True,color='r') | |
| bovy_plot.bovy_end_print('/home/bovy/Desktop/test.png') | |
| if __name__ == '__main__': | |
| plot_dvcapprox() |
| #Plot the *correct* surface-mass density profile | |
| import math | |
| import cPickle as pickle | |
| import numpy | |
| from scipy import integrate, interpolate | |
| from galpy.util import bovy_plot | |
| import monoAbundanceMW | |
| from matplotlib import pyplot | |
| _R0= 8. #kpc | |
| _hR= 3.8 #kpc | |
| _hz= 0.9 #kpc | |
| _TWOPIG= 2.*numpy.pi*4.302 | |
| _NMC= 10001 | |
| _OUTPUTRESULT= False | |
| def sigmaz(z,hr=_hR,tsigmaUW=None,hrthick=_hR,varhz=False): | |
| k3= 1./hr+1./hrthick-1./_R0 | |
| tsigmaW= sigmaW(z) | |
| dsigmaWdz= 2.7 | |
| if tsigmaUW is None: tsigmaUW= sigmaUW(z) | |
| if varhz: | |
| hz= monoAbundanceMW.meanhz(z*1000.)/1000. | |
| else: | |
| hz= _hz | |
| return 1./_TWOPIG*(tsigmaW**2./hz - 2.*tsigmaW*dsigmaWdz + k3*tsigmaUW) | |
| def mcsigmaz(z,hr=_hR,tsigmaUW=None,hrthick=_hR,varhz=False): | |
| k3= 1./hr+1./hrthick-1./_R0 | |
| #Draw coefficients | |
| wa= numpy.random.normal(size=_NMC)*0.8+40.6 | |
| wb= numpy.random.normal(size=_NMC)*0.3+2.7 | |
| uwa= numpy.random.normal(size=_NMC)*100.+1522. | |
| uwb= numpy.random.normal(size=_NMC)*30.+366. | |
| samples= numpy.zeros(_NMC) | |
| if varhz: | |
| hz= monoAbundanceMW.meanhz(z*1000.)/1000. | |
| else: | |
| hz= _hz | |
| for ii in range(_NMC): | |
| tsigmaW= sigmaW(z,a=wa[ii],b=wb[ii]) | |
| dsigmaWdz= wb[ii] | |
| tsigmaUW= sigmaUW(z,a=uwa[ii],b=uwb[ii]) | |
| samples[ii]= 1./_TWOPIG*(tsigmaW**2./hz - 2.*tsigmaW*dsigmaWdz + k3*tsigmaUW) | |
| return samples | |
| def mbsigmaz(z,hr=_hR,sigmaUW=None): | |
| k1= 3./_R0/hr-2./hr**2. | |
| k2= -1/_R0/hr | |
| k3= 3./hr-2./_R0 | |
| #k3= 2./hr-1./_R0 | |
| sigmaW= 40.6+2.7*(z-2.5) | |
| dsigmaWdz= 2.7 | |
| if sigmaUW is None: sigmaUW= 1522.+366*(z-2.5) | |
| planarpart= numpy.zeros(len(z)) | |
| for ii in range(len(z)): planarpart[ii]= integrate.quad((lambda x: k1*sigmaU(x)**2+k2*sigmaV(x)**2.),0.,z[ii])[0] | |
| return 1./_TWOPIG*(planarpart+sigmaW**2./_hz - 2.*sigmaW*dsigmaWdz + k3*sigmaUW) | |
| def mcmbsigmaz(z,hr=_hR,tsigmaUW=None): | |
| k1= 3./_R0/hr-2./hr**2. | |
| k2= -1/_R0/hr | |
| k3= 3./hr-2./_R0 | |
| #Draw coefficients | |
| ua= numpy.random.normal(size=_NMC)*3.2+82.9 | |
| ub= numpy.random.normal(size=_NMC)*1.1+6.3 | |
| va= numpy.random.normal(size=_NMC)*3.1+62.2 | |
| vb= numpy.random.normal(size=_NMC)*1.0+4.1 | |
| wa= numpy.random.normal(size=_NMC)*0.8+40.6 | |
| wb= numpy.random.normal(size=_NMC)*0.3+2.7 | |
| uwa= numpy.random.normal(size=_NMC)*100.+1522. | |
| uwb= numpy.random.normal(size=_NMC)*30.+366. | |
| samples= numpy.zeros(_NMC) | |
| for ii in range(_NMC): | |
| tsigmaW= sigmaW(z,a=wa[ii],b=wb[ii]) | |
| dsigmaWdz= wb[ii] | |
| tsigmaUW= sigmaUW(z,a=uwa[ii],b=uwb[ii]) | |
| planarpart= integrate.quad((lambda x: k1*sigmaU(x,a=ua[ii],b=ub[ii])**2+k2*sigmaV(x,a=va[ii],b=vb[ii])**2.),0.,z)[0] | |
| samples[ii]= 1./_TWOPIG*(planarpart+tsigmaW**2./_hz - 2.*tsigmaW*dsigmaWdz + k3*tsigmaUW) | |
| return samples | |
| def sigmaW(z,a=40.6,b=2.7): | |
| return a+b*(z-2.5) | |
| def sigmaU(z,a=82.9,b=6.3): | |
| return a+b*(z-2.5) | |
| def sigmaV(z,a=62.2,b=4.1): | |
| return a+b*(z-2.5) | |
| def sigmaUW(z,a=1522.,b=366.): | |
| return a+b*(z-2.5) | |
| def meanV(z): | |
| return 220.-25.-30.*z | |
| def rhoOM(z): | |
| return 20.6*10.**-3.*(8.01**2./(8.01**2.+_R0**2.+z**2.)) | |
| def sigmaVIS(z): | |
| #Needs to be normalized later | |
| return 2000.*integrate.quad(rhoVIS,0.,z)[0] | |
| def rhoVIS(z): | |
| #Thin + thick + halo from Juric+ 2008 | |
| return numpy.exp(-z/0.3)+.12*numpy.exp(-z/0.9)+0.005*(_R0/numpy.sqrt(_R0**2.+(z/0.64)**2.))**2.8 | |
| def sigmaOM(z): | |
| return 2000.*integrate.quad(rhoOM,0.,z)[0] | |
| def rhoDM(z,R,alpha,beta,gamma,Rc,rhoo,q): | |
| r= numpy.sqrt(R**2.+(z/q)**2.) | |
| return rhoo*(r/_R0)**-alpha*((1.+(r/Rc)**beta)/(1.+(_R0/Rc)**beta))**-gamma | |
| def sigmaDM(z,R=_R0,alpha=1.,beta=1.,gamma=2.,Rc=10.8,rhoo=8.4*10.**-3.,q=1.): | |
| return 2000.*integrate.quad(rhoDM,0.,z,args=(R,alpha,beta,gamma,Rc,rhoo,q))[0] | |
| def plot_sigmaz(): | |
| zs= numpy.linspace(1.5,4.,1001) | |
| sigmazs= sigmaz(zs) | |
| #Also generate samples | |
| sampleszs= numpy.array([1.8,2.8,3.8]) | |
| dev= numpy.zeros((2,len(sampleszs))) | |
| sigmazatsamples= sigmaz(sampleszs) | |
| slopes= numpy.zeros(_NMC) | |
| for ii in range(len(sampleszs)): | |
| sigmazsamples= mcsigmaz(sampleszs[ii]) | |
| if ii == 0.: | |
| slopes-= sigmazsamples | |
| elif ii == (len(sampleszs)-1): | |
| slopes+= sigmazsamples | |
| dev[0,ii]= sigmazatsamples[ii]-sorted(sigmazsamples)[int(numpy.round(0.16*_NMC))] | |
| dev[1,ii]= sorted(sigmazsamples)[int(numpy.round(0.84*_NMC))]-sigmazatsamples[ii] | |
| print 'slope', (sigmazs[-1]-sigmazs[0])/(zs[-1]-zs[0]), '+/-',\ | |
| numpy.std(slopes)/(sampleszs[-1]-sampleszs[0]) | |
| mbsigmazs= mbsigmaz(zs) | |
| #Also generate samples | |
| mbdev= numpy.zeros((2,len(sampleszs))) | |
| mbsigmazatsamples= mbsigmaz(sampleszs) | |
| for ii in range(len(sampleszs)): | |
| sigmazsamples= mcmbsigmaz(sampleszs[ii]) | |
| mbdev[0,ii]= mbsigmazatsamples[ii]-sorted(sigmazsamples)[int(numpy.round(0.16*_NMC))] | |
| mbdev[1,ii]= sorted(sigmazsamples)[int(numpy.round(0.84*_NMC))]-mbsigmazatsamples[ii] | |
| #sigmazsuwzero= sigmaz(zs,sigmaUW=0.) | |
| sigmazs2kpc= sigmaz(zs,hrthick=2.) | |
| #Also generate samples | |
| dev2kpc= numpy.zeros((2,len(sampleszs))) | |
| sigmaz2kpcatsamples= sigmaz(sampleszs,hrthick=2.) | |
| slopes= numpy.zeros(_NMC) | |
| for ii in range(len(sampleszs)): | |
| sigmazsamples= mcsigmaz(sampleszs[ii],hrthick=2.) | |
| if ii == 0.: | |
| slopes-= sigmazsamples | |
| elif ii == (len(sampleszs)-1): | |
| slopes+= sigmazsamples | |
| dev2kpc[0,ii]= sigmaz2kpcatsamples[ii]-sorted(sigmazsamples)[int(numpy.round(0.16*_NMC))] | |
| dev2kpc[1,ii]= sorted(sigmazsamples)[int(numpy.round(0.84*_NMC))]-sigmaz2kpcatsamples[ii] | |
| print 'slope', (sigmazs2kpc[-1]-sigmazs2kpc[0])/(zs[-1]-zs[0]), '+/-',\ | |
| numpy.std(slopes)/(sampleszs[-1]-sampleszs[0]) | |
| sigmazsvarhz= sigmaz(zs,varhz=True,hrthick=2.) | |
| #Also generate samples | |
| devvarhz= numpy.zeros((2,len(sampleszs))) | |
| sigmazvarhzatsamples= sigmaz(sampleszs,hrthick=2.,varhz=True) | |
| slopes= numpy.zeros(_NMC) | |
| for ii in range(len(sampleszs)): | |
| sigmazsamples= mcsigmaz(sampleszs[ii],hrthick=2.,varhz=True) | |
| if ii == 0.: | |
| slopes-= sigmazsamples | |
| elif ii == (len(sampleszs)-1): | |
| slopes+= sigmazsamples | |
| dev[0,ii]= sigmazvarhzatsamples[ii]-sorted(sigmazsamples)[int(numpy.round(0.16*_NMC))] | |
| dev[1,ii]= sorted(sigmazsamples)[int(numpy.round(0.84*_NMC))]-sigmazvarhzatsamples[ii] | |
| print 'slope', (sigmazsvarhz[-1]-sigmazsvarhz[0])/(zs[-1]-zs[0]), '+/-',\ | |
| numpy.std(slopes)/(sampleszs[-1]-sampleszs[0]) | |
| #Calculate expected | |
| tzs= numpy.linspace(0.,5.,1001) | |
| vis= numpy.zeros(len(tzs)) | |
| for ii in range(len(tzs)): vis[ii]= sigmaVIS(tzs[ii]) | |
| vis/= sigmaVIS(1.1)/40. #Normalization at z=1.1 kpc | |
| vis+= 13. #ISM | |
| om= numpy.zeros(len(tzs)) | |
| for ii in range(len(tzs)): om[ii]= sigmaOM(tzs[ii])+vis[ii] | |
| shm= numpy.zeros(len(tzs)) | |
| for ii in range(len(tzs)): shm[ii]= sigmaDM(tzs[ii])+vis[ii] | |
| n97= numpy.zeros(len(tzs)) | |
| for ii in range(len(tzs)): n97[ii]= sigmaDM(tzs[ii],Rc=20.,rhoo=6.1*10.**-3.)+vis[ii] | |
| min= numpy.zeros(len(tzs)) | |
| for ii in range(len(tzs)): min[ii]= sigmaDM(tzs[ii],alpha=0.,beta=2.,gamma=1.,Rc=5.,rhoo=5.3*10.**-3.)+vis[ii] | |
| bovy_plot.bovy_print() | |
| bovy_plot.bovy_plot(zs,sigmazs,'k-', | |
| xrange=[0.,5.], | |
| yrange=[20.,160.], | |
| xlabel=r'$Z\ [\mathrm{kpc}]$', | |
| ylabel=r'$\Sigma(Z)\ [M_\odot\ \mathrm{pc}^{-2}]$', | |
| zorder=2) | |
| #Errors | |
| pyplot.errorbar(sampleszs,sigmazatsamples, | |
| yerr=dev, | |
| fmt=',',marker='None',ecolor='k',lw=0.5) | |
| bovy_plot.bovy_plot(zs,mbsigmazs,'k-',overplot=True) | |
| #Errors | |
| pyplot.errorbar(sampleszs,mbsigmazatsamples, | |
| yerr=mbdev, | |
| fmt=',',marker='None',ecolor='k',lw=0.5) | |
| bovy_plot.bovy_plot(zs,sigmazs2kpc,'k--',overplot=True) | |
| #Errors | |
| pyplot.errorbar(sampleszs,sigmaz2kpcatsamples, | |
| yerr=dev2kpc, | |
| fmt=',',marker='None',ecolor='k',lw=0.5) | |
| #bovy_plot.bovy_plot(zs,sigmazsvarhz,'k-.',overplot=True) | |
| #bovy_plot.bovy_plot(zs,sigmazsuwzero,'k-.',overplot=True,zorder=1) | |
| bovy_plot.bovy_plot(tzs,vis, | |
| '-',color='0.7',lw=2.,overplot=True,zorder=0) | |
| bovy_plot.bovy_plot(tzs,shm,'-',color='0.7',lw=2.,overplot=True,zorder=0) | |
| bovy_plot.bovy_plot(tzs,n97,'-',color='0.7',lw=2.,overplot=True,zorder=0) | |
| bovy_plot.bovy_plot(tzs,min,'-',color='0.7',lw=2.,overplot=True,zorder=0) | |
| bovy_plot.bovy_plot(tzs,om,'-',color='0.7',lw=2.,overplot=True,zorder=0) | |
| #Load approximation and apply | |
| savefile= open('approx.sav','rb') | |
| xzs= pickle.load(savefile) | |
| x1= pickle.load(savefile) | |
| x2= pickle.load(savefile) | |
| savefile.close() | |
| #Interpolate | |
| x1Interp= interpolate.InterpolatedUnivariateSpline(xzs,x1) | |
| x2Interp= interpolate.InterpolatedUnivariateSpline(xzs,x2) | |
| ul= sigmazs/(1.-x1Interp(zs/_R0)) | |
| ll= sigmazs/(1.-x2Interp(zs/_R0)) | |
| print 'slope', (ul[-1]-ul[0])/(zs[-1]-zs[0]) | |
| print 'slope', (ll[-1]-ll[0])/(zs[-1]-zs[0]) | |
| #Plot as fill_between | |
| pyplot.fill_between(zs,ul,ll,color='0.85',zorder=-10) | |
| #Labels | |
| bovy_plot.bovy_text(4.1,147.,r'$\mathrm{OM}$',rotation=30.) | |
| bovy_plot.bovy_text(4.1,130.,r'$\mathrm{SHM}$',rotation=25.) | |
| bovy_plot.bovy_text(4.1,110.,r'$\mathrm{N97}$',rotation=17.) | |
| bovy_plot.bovy_text(4.1,95.,r'$\mathrm{MIN}$',rotation=15.) | |
| bovy_plot.bovy_text(4.1,60.,r'$\mathrm{VIS}$') | |
| #Lines | |
| bovy_plot.bovy_text(1.4,35.5,r'$\mathrm{incorrect}\ \frac{\partial \bar{V}}{\partial R} = 0$'+'\n'+r'$\mathrm{approximation}$') | |
| bovy_plot.bovy_plot([1.6,1.8],[48.,mbsigmaz(numpy.array([1.8]))],'-',color='0.8',overplot=True) | |
| bovy_plot.bovy_text(2.7,68.,r'$\mathrm{correct}\ \frac{\partial V_c}{\partial R} = 0$'+'\n'+r'$\mathrm{approximation}$') | |
| bovy_plot.bovy_plot([2.5,2.65],[sigmaz(2.5),79.],'-',color='0.8',overplot=True) | |
| bovy_plot.bovy_text(1.4,125.,r'$\mathrm{correct}$'+'\n'+r'$+ h_{R,\mathrm{thick}} = 2\ \mathrm{kpc}$') | |
| bovy_plot.bovy_plot([1.8,2.4],[123.,sigmaz(2.4,2.)],'-',color='0.8',overplot=True) | |
| bovy_plot.bovy_text(.4,99.,r'$\mathrm{effect\ of}\ \frac{\partial V_c}{\partial R} = 0$'+'\n'+r'$\mathrm{approximation}$') | |
| bovy_plot.bovy_plot([1.7,2.5],[97.,sigmaz(2.5)/(1.-x1Interp(2.5/_R0))], | |
| '-',color='0.8',overplot=True) | |
| #TITLE | |
| #bovy_plot.bovy_text('$\mathrm{Dark\ matter\ near\ the\ Sun\ is\ A\!\!-\!\!OK}$', | |
| #title=True) | |
| #bovy_plot.bovy_end_print('/home/bovy/Desktop/correct_monibidin12.png') | |
| bovy_plot.bovy_end_print('test.png') | |
| #bovy_plot.bovy_end_print('../tex/correct_monibidin12.ps') | |
| if _OUTPUTRESULT: | |
| output= numpy.zeros((len(zs),4)) | |
| output[:,0]= zs | |
| output[:,1]= sigmazs2kpc | |
| output[:,2]= ll/sigmazs*sigmazs2kpc | |
| output[:,3]= ul/sigmazs*sigmazs2kpc | |
| numpy.savetxt('bovytremaine_mb12_Sigmaz.dat',output,delimiter='|') | |
| if __name__ == '__main__': | |
| numpy.random.seed(1) | |
| plot_sigmaz() |