<jupyter_start><jupyter_text>Temporal Smearing $\tau$Here we compute the $\tau$ in post-processing, for models that cannot provide $\tau$ immediately, as it depends on distance to source and lensing material.We use equations derived by Macquart & Koay 2013. The following plot shows the lense position for which the effective lense distance $D_{eff}$ is maximized.<jupyter_code>zs = np.arange(1.,10.,1)
for z_s, color in zip(zs, rainbow(np.linspace(0,1,zs.size))):
    z_L = np.arange(0.1,z_s,0.1)
    d = Deff( z_L=z_L, z_s=z_s )
    plt.plot(z_L,d, color=color )
Colorbar(zs, label=r"$z_{source}$")
plt.xlabel(r"$z_{lense}$")
plt.ylabel(r"$D_{eff} = D_L D_{LS} / D_S$ / Gpc")<jupyter_output><empty_output><jupyter_text>IGMHere, we translate the likelihood of SM in the IGM, $P(SM_{IGM}|z_s)$, to likelihood of scattering time in the IGM, $P(\tau(SM_{IGM},z_s)|z_s)$.The likelihoods are of same shape, so the range of bins in $P$ has to be shifted and $P$ renormalized accordingly in order to arrive at $P(\tau)$.The plasma of the IGM is elongated along the whole LoS.We thus choose the effective $D_{\rm eff}$  such that it maximizes the observed effects of scattering.<jupyter_code>from scipy.optimize import fmin_slsqp as  fmin

def Deff_max( z_s, accuracy=1e-5 ):
    ### returns optimal redshift and Deff for lense between observer and source at redshift z_s
    Deff_min = lambda z_L : -Deff( z_s=z_s, z_L=z_L )
    ### find z_L that maximizes D_eff(z_s,z_L) for given z_s
    z_L, = fmin( Deff_min, z_s/2, iprint=0, acc=accuracy)
    return z_L, Deff( z_s=z_s, z_L=z_L )

model = 'primordial'  ## other IGM models only differ in RM
fig, axs = plt.subplots( 1, 2, figsize=(8,3))
## for each redshift
for z_s, color in zip(redshift_bins, rainbow( np.linspace(0,1,len(redshift_bins)))):
    ### read likelihood of SM: P(SM_IGM|z_s)
    P, SM = GetLikelihood_IGM( redshift=z_s, model=model, nside=4, measure='SM', absolute=False)
    PlotLikelihood( x=SM, P=P, measure='SM', color=color, ax=axs[0] )
    ### find optimal lense position
    z_L, D_eff = Deff_max( z_s )
    ### translate x-axis   SM -> tau
    x_tau = ScatteringTime( SM=SM, redshift=z_L, D_eff=D_eff)
    ### renormalize
    P *= np.diff(SM)/np.diff(x_tau)
    PlotLikelihood( x=x_tau, P=P, measure='tau', color=color, ax=axs[1] )
    Write2h5( filename=likelihood_file_IGM, datas=[P,x_tau], 
              keys=[ KeyIGM( redshift=z_s, model=model, nside=4, measure='tau', axis=axis) for axis in ['P','x'] ] 
            )

Colorbar(redshift_bins)
plt.tight_layout()
plt.show()<jupyter_output>/data/PreFRBLE/venv_PreFRBLE/lib/python3.4/site-packages/yt/utilities/cosmology.py:623: RuntimeWarning: invalid value encountered in log10
  zbins = np.logspace(np.log10(a + 1), np.log10(b + 1), bins) - 1<jupyter_text>Host GalaxyHere, we translate the likelihood of SM from host galaxies at redshift $z$, $P(SM_{Host}|z)$, to likelihood of scattering time, $P(\tau(SM_{Host},r_{gal})|z)$.For this we use the average size of ellispoids representing the galaxy, i. e. 2.7*half_mass radius, marking the distance where surface density falls below 1%.Since the host galaxies provide negligible $\tau$, choosing the average size for galaxies of all morphologies and masses is a justified simplification.<jupyter_code>model = 'Rodrigues18'

from PreFRBLE.file_system import Rodrigues_file_rgal


d = np.genfromtxt(Rodrigues_file_rgal, names=True)
r_gals = d['r_gal'] # kpc  half mass radius
r_gals *= 2.7e-6  # Gpc, required for D_eff, *2.7 to meet assumed size of galactic ellipsoid

for z, r_gal in zip( redshift_bins, r_gals ):
    ## read likelihood SM
    P, SM = GetLikelihood_Host( model=model, measure='SM', redshift=z )
    ## transform x-values  SM -> tau
    tau = ScatteringTime( SM=SM, redshift=z, D_eff=r_gal )
    ## renormalize
    P *= np.diff(SM)/np.diff(tau)
    
    ## write to file
    keys = [ KeyHost( model=model, redshift=z, measure='tau', axis=axis) for axis in ['P','x'] ]
    Write2h5( filename=likelihood_file_galaxy, keys=keys, datas=[P,tau] )

## check results
PlotLikelihoodEvolution( measure='SM', scenario={'Host':[model]} )
plt.ylim(1e-3,1)
plt.show()
PlotLikelihoodEvolution( measure='tau', scenario={'Host':[model]} )
plt.ylim(1e-3,1)
plt.show()<jupyter_output><empty_output><jupyter_text>Intervening GalaxiesHere, we translate the likelihood of SM from galaxies intervening LoS at $z_L$, $P(SM_{Inter}|z_L)$, to likelihood of scattering time for source at $z_s$, $P(\tau(SM_{Inter},z_s)|z_s)$.For the intervening galaxies, $D_{\rm eff}$ is determined by the redshift of the intervening galaxy.The contribution of galaxies at different redshift has to be weighed by the expected amount of galaxies intersecting LoS at this redshift $\pi_{\rm inter}$.This is given by the integrand of the average amount of galaxies along each LoS $$N_{\rm Inter} = \int \pi~n_{\rm gal}~r^2_{\rm gal} \frac{d_H(z)}{1+z}\text{d}z = \int \pi_{\rm inter}(z) \text{d}z = \sum n_{\rm inter}.$$<jupyter_code>### plot amount of galaxies expected to intervene LoS at redshift z, n_inter(z),
### as well as expectet number of galaxies in LoS to source at redshift z, N_inter(z)

model = 'Rodrigues18'

#n_inter = nInter( redshift_bins[-1], r=0.01/hubble_constant, n=0.02*hubble_constant**3 )
n_inter = nInter( redshift_bins[-1], model=model )
N_inter = np.cumsum( n_inter )

plt.plot( redshift_bins, N_inter, label=r"$N_{Inter}$" )
plt.plot( redshift_bins, n_inter, linestyle=':', label=r"$n_{Inter}$" )
plt.xlabel('redshift')
plt.ylabel(r"$N_{Inter}$")
plt.yscale( 'log')
plt.legend()
## compute P(t|z) and write to file

model = 'Rodrigues18'
fig, axs = plt.subplots( 1, 2, figsize=(8,3))

weights = nInter(redshift_bins[-1], model=model )  ## intersection probability for each redshift bin
redshift_L = redshift_bins - np.diff(redshift_range)/2  ## idealized possible intersection redshift, central redshift of bin

## for each source redshift
for i_s, (z_s, color) in enumerate( zip( redshift_bins, rainbow( np.linspace(0,1,len(redshift_bins)) ) ) ):
    Ps, xs = [], []
    ## for each possible redshift of intersection < z_s
    for i_L , (z_L, color_L) in enumerate( zip( redshift_L[:i_s+1], rainbow( np.linspace(0,1,len(redshift_bins)) ) ) ):
        ## read raw likelihood of SM from intervening galaxies at that redshift: P(SM_inter|z_L)
        P, SM = GetLikelihood_inter( redshift=redshift_bins[i_L], model=model, measure='SM')
        if i_s == len(redshift_bins) - 1:
            PlotLikelihood( x=SM, P=P, measure='SM', color=color_L, ax=axs[0] )
        ## calculate effective lense distance
        D_eff = Deff( z_s=z_s, z_L=z_L )
        ### translate x-axis
        x_tau = ScatteringTime( SM=SM, redshift=z_L, D_eff=D_eff )
        ## renormalize
        P *= np.diff(SM)/np.diff(x_tau)
        Ps.append(P)
        xs.append(x_tau)
    ## integrate results for all intervening redshifts, weighing and renormalizing with differential and cumulative intersection probability, respectively
    ### result is renormalized to 1, can be read out renormalized to N_inter with GetLikelihood_Inter( ..., N_inter=True )
    P, x = LikelihoodsAdd( Ps, xs, weights=weights[:i_s+1], renormalize=1)
    PlotLikelihood( x=x, P=P, measure='tau', color=color, ax=axs[1] )
    Write2h5( filename=likelihood_file_galaxy, datas=[P,x], 
             keys=[ KeyInter( redshift=z_s, model=model.split('/')[0], measure='tau', axis=axis) for axis in ['P','x'] ] 
            )

Colorbar(redshift_bins)
plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>validity of equationsHere we show, that we can always use equations in Marcquardt & Koay 2013 for the case $r_{\rm diff} > l_0$<jupyter_code>lambda_0 = 0.23 # m
#lambda_0 = 0.50 # m
l_0 = 1.496e11  # m  == 1 AU
kpc2m = 3.086e19 # m
fig, ax = plt.subplots()

SM_unit = 3.241e-8 # kpc m^-20/3 to 10^12 m^-17/3

## range of possible SM_eff 
SM_eff = 10.**np.linspace(10,30,100)/kpc2m ## kpc m^-20/3

#r_lo = ( 8e9 * lambda_0**-1 * (SM_unit*SM_eff)**-0.5 )**1.2 ## r < l_0, 
r_lo = 8e9 * lambda_0**-1 * (SM_unit*SM_eff)**-0.5  ## r < l_0, 
r_hi = 3.7e9 * lambda_0**-1.2 * (SM_unit*SM_eff)**-0.6      ## r > l_0

plt.loglog( SM_eff, r_lo/l_0, label=r"$r_{\rm diff} < l_0$", lw=3 )
plt.plot( SM_eff, r_hi/l_0, label=r"$r_{\rm diff} > l_0$" , lw=3)
plt.axhline(y=l_0/l_0, label=r"$l_0$", lw=3)
plt.ylabel('$r$ / AU',fontsize=16)
plt.xlabel('SM / (kpc m$^{-20/3}$)',fontsize=16)
plt.tick_params(axis='both', which='major', labelsize=16)

plt.legend(fontsize=16)
#AllSidesTicks(ax)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14479690975067305, 'fraction_numerical': 0.023293924850755005, 'mean_word_length': 4.548051948051948, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '74012', 'license': 'MIT', 'n_tokens_mistral': 3378, 'n_tokens_neox': 3171, 'n_words': 857, 'path': 'notebooks/model/temporal_smearing.ipynb', 'repo_name': 'FRBs/PreFRBLE', 'size': '8543'}