The file src/32-dark/make_figures.py of the validation package of the paper 32-dark, rendered with an anchor per line; the file itself is beside this page.
1#!/usr/bin/env python32# =====================================================================3# make_figures.py -- consolidated, de-duplicated dark-sector figures.4#5# Replaces the redundant 3-panel cons2_rar_btfr (RAR|BTFR|rotcurve) and6# the overlapping cons3_interpolation (RAR|angle). The RAR appeared7# twice and the BTFR panel was "just a slope". New, non-redundant set:8#9# fig_rar.pdf (phenomenology): RAR (with the rational-alt knee10# and both asymptotes) | exponential-disk rotation curve11# fig_mechanism.pdf (first principles): meridian first-passage collapse12# (walk -> e^{-a sqrt(2s)}) | chart angle alpha(x)13#14# The BTFR (v^4=GM a0) is a one-line corollary and stays in the text.15# =====================================================================16import numpy as np17import matplotlib; matplotlib.use('Agg'); import matplotlib.pyplot as plt1819plt.rcParams.update({'font.size':10,'axes.titlesize':10,'axes.labelsize':10,20 'legend.fontsize':8,'xtick.labelsize':8.5,'ytick.labelsize':8.5})2122c=2.998e8; G=6.674e-11; Msun=1.989e30; kpc=3.086e19; Mpc=1000*kpc23a0=c*(67.4*1000/Mpc)/(2*np.pi) # H0 = 67.4 km/s/Mpc, the rate the octant lemma entails (14-entr)24nu =lambda x:1.0/(1-np.exp(-np.sqrt(x))) # derived = McGaugh RAR fit25nu_rat =lambda x:0.5*(1+np.sqrt(1+4.0/x)) # 'simple' rational alt (same limits)26alpha =lambda x:np.arcsin(np.sqrt(np.exp(-np.sqrt(x))))2728# ---------------------------------------------------------------------29# FIGURE 1 -- galactic phenomenology: RAR | rotation curve30# ---------------------------------------------------------------------31fig,ax=plt.subplots(1,2,figsize=(9.6,4.0))3233gb=np.logspace(-13,-8,400) # baryonic acceleration, m/s^234x=gb/a035ax[0].loglog(gb, gb*nu(x), 'b', lw=1.8, label=r'two-chart $g_{\rm obs}=g_{\rm bar}/(1-e^{-\sqrt{x}})$')36ax[0].loglog(gb, gb, 'k--', lw=0.9, label='Newtonian (slope 1)')37ax[0].loglog(gb, np.sqrt(gb*a0), 'g-.', lw=1.0, label=r'deep $\sqrt{g_{\rm bar}a_0}$ (slope $\frac{1}{2}$)')38ax[0].loglog(gb, gb*nu_rat(x), color='m', ls=':', lw=1.3, label=r'rational $\mu$ (alt. knee)')39ax[0].axvline(a0,color='gray',ls=':',lw=0.8); ax[0].text(a0*1.15,1.5e-13,r'$a_0$',color='gray',fontsize=9)40ax[0].set_xlabel(r'$g_{\rm bar}\ \mathrm{[m\,s^{-2}]}$'); ax[0].set_ylabel(r'$g_{\rm obs}\ \mathrm{[m\,s^{-2}]}$')41ax[0].set_title('Radial acceleration relation'); ax[0].legend(loc='upper left'); ax[0].set_xlim(1e-13,1e-8)4243from scipy.special import i0,i1,k0,k144Mdisk,Rd=5e10*Msun,3*kpc45r=np.linspace(0.1,30,300)*kpc; xd=r/(2*Rd)46vb2=(G*Mdisk/Rd)*xd**2*(i0(xd)*k0(xd)-i1(xd)*k1(xd))*247gbar=vb2/r; vobs=np.sqrt(r*gbar*nu(gbar/a0))48ax[1].plot(r/kpc, np.sqrt(vb2)/1000, 'k--', lw=1.1, label='baryonic')49ax[1].plot(r/kpc, vobs/1000, 'b', lw=1.8, label='registered (two-chart)')50ax[1].set_xlabel('$r$ [kpc]'); ax[1].set_ylabel('$v$ [km s$^{-1}$]')51ax[1].set_title(r'Rotation curve (exponential disk, $5\times10^{10}M_\odot$)')52ax[1].legend(loc='lower right'); ax[1].set_ylim(0,None)53fig.tight_layout(); fig.savefig('../figures/fig_rar.pdf'); print("fig_rar.pdf saved")5455# ---------------------------------------------------------------------56# FIGURE 2 -- the interpolation from first principles:57# first-passage collapse | chart angle58# ---------------------------------------------------------------------59fig,ax=plt.subplots(1,2,figsize=(9.6,4.0))6061def escape_prob(a,q,trials=20000,seed=0):62 rng=np.random.default_rng(seed); esc=063 for _ in range(trials):64 pos=065 while True:66 if rng.random()<q: break67 pos+=1 if rng.random()<0.5 else -168 if pos>=a: esc+=1; break69 return esc/trials70pts=[]71for a in (3,5,8,10):72 for q in (0.01,0.02,0.04):73 pts.append((a*np.sqrt(2*q), escape_prob(a,q,seed=a*100+int(q*1000))))74pts=np.array(pts); u=np.linspace(0,3,100)75ax[0].plot(pts[:,0],pts[:,1],'o',ms=5,color='steelblue',label='killed-walk simulation')76ax[0].plot(u,np.exp(-u),'r-',lw=1.6,label=r'$e^{-A}$, $A=a\,\sqrt{2s}=\sqrt{x}$')77ax[0].set_xlabel(r'$A=a\,\sqrt{2s}$'); ax[0].set_ylabel('escape (non-registration) prob.')78ax[0].set_title('Meridian first-passage law (Prop. 3)'); ax[0].legend(loc='upper right')7980xx=np.logspace(-3,3,400)81ax[1].semilogx(xx, np.degrees(alpha(xx)), 'r', lw=1.8)82ax[1].axhline(45,color='gray',ls=':',lw=0.7); ax[1].axhline(37,color='gray',ls=':',lw=0.5)83ax[1].axvline(1,color='gray',ls=':',lw=0.7); ax[1].text(1.3,5,r'$a_0$',color='gray',fontsize=9)84ax[1].text(1e-3*1.4,40,r'spectral',fontsize=8,color='dimgray'); ax[1].text(1e2,6,r'coordinate',fontsize=8,color='dimgray')85ax[1].set_xlabel(r'$x=g_{\rm bar}/a_0$'); ax[1].set_ylabel(r'chart angle $\alpha$ [deg]')86ax[1].set_title(r'Chart angle: $\sin^2\alpha=e^{-\sqrt{x}}$ ($0^\circ$ space, $90^\circ$ spectrum)')87ax[1].set_ylim(0,92)88fig.tight_layout(); fig.savefig('../figures/fig_mechanism.pdf'); print("fig_mechanism.pdf saved")89