From 689e16b51192aff5d1777545acc1424ca495fc2f Mon Sep 17 00:00:00 2001 From: Daniel Gilman Date: Thu, 29 Aug 2024 08:43:33 -1000 Subject: [PATCH 1/3] implement new galacticus truncation routine --- pyHalo/Halos/HaloModels/TNFW.py | 24 +- .../transfer_function_density_profile.py | 834 ++++++++++++++++++ pyHalo/Halos/halo_base.py | 2 + pyHalo/Halos/tidal_truncation.py | 26 +- .../test_halos/test_galacticus_truncation.py | 27 + 5 files changed, 898 insertions(+), 15 deletions(-) create mode 100644 pyHalo/Halos/galacticus_truncation/transfer_function_density_profile.py diff --git a/pyHalo/Halos/HaloModels/TNFW.py b/pyHalo/Halos/HaloModels/TNFW.py index a748fb2..bafd719 100644 --- a/pyHalo/Halos/HaloModels/TNFW.py +++ b/pyHalo/Halos/HaloModels/TNFW.py @@ -1,7 +1,7 @@ from pyHalo.Halos.halo_base import Halo import numpy as np from pyHalo.Halos.tnfw_halo_util import tnfw_mass_fraction - +from scipy.integrate import quad class TNFWFieldHalo(Halo): @@ -99,6 +99,7 @@ def profile_args(self): if not hasattr(self, '_profile_args'): truncation_radius_kpc = self._truncation_class.truncation_radius_halo(self) + self._profile_args = (self.c, truncation_radius_kpc) return self._profile_args @@ -113,10 +114,17 @@ def bound_mass(self): Computes the mass inside the virial radius (with truncation effects included) :return: the mass inside r = c * r_s """ - if hasattr(self, '_kwargs_lenstronomy'): - tau = self._kwargs_lenstronomy[0]['r_trunc'] / self._kwargs_lenstronomy[0]['Rs'] - else: - params_physical = self.params_physical - tau = params_physical['r_trunc_kpc'] / params_physical['rs'] - f = tnfw_mass_fraction(tau, self.c) - return f * self.mass + def _integrand(r, rs, rhos, rt): + x = r/rs + tau = rt/rs + return 4 * tau **2 * np.pi * r ** 2 * rhos / x / (1+x)**2 / (x**2 + tau**2) + + params_physical = self.params_physical + rhos = params_physical['rhos'] + rs = params_physical['rs'] + rt = params_physical['r_trunc_kpc'] + r200 = rs * self.c + #tau = params_physical['r_trunc_kpc'] / params_physical['rs'] + #f = tnfw_mass_fraction(tau, self.c) * self._rescale_norm + bound_mass = quad(_integrand, 0, r200, args=(rs, rhos, rt))[0] + return bound_mass diff --git a/pyHalo/Halos/galacticus_truncation/transfer_function_density_profile.py b/pyHalo/Halos/galacticus_truncation/transfer_function_density_profile.py new file mode 100644 index 0000000..b65cfd7 --- /dev/null +++ b/pyHalo/Halos/galacticus_truncation/transfer_function_density_profile.py @@ -0,0 +1,834 @@ +import warnings +import numpy as np +from scipy.special import hyp2f1 +from scipy.integrate import quad +from scipy.interpolate import interp1d + + +# Implement the transfer function of subhalo density profile. + +def _r_te_Du_2024(x, alpha=1.0, beta=3.0, gamma=1.0, delta=2.0): + """ + Return the effective tidal (truncation) radius of a subhalo given the bound mass fraction. + :param x: bound mass fraction + :param alpha: alpha parameter in the generalized NFW profile + :param beta: beta parameter in the generalized NFW profile + :param gamma: gamma parameter in the generalized NFW profile + :param delta: power index of the density truncation + :return: the effective tidal (truncation) radius in units of infall virial radius + """ + error_status = 0 + + if (delta == 2.0): + if (alpha == 1.0): + if (beta == 3.0): + if (gamma == 0.0): + A = -9.98429803e-01 + B = 5.50290877e-04 + C = 1.23471540e+00 + elif (gamma == 0.5): + A = -0.19210388 + B = 0.42676044 + C = 1.46250084 + elif (gamma == 1.0): + A = 0.68492777 + B = 0.66438857 + C = 2.07766512 + elif (gamma == 1.5): + A = 0.9839032 + B = 0.76882264 + C = 2.51841765 + else: + error_status = 1 + elif (beta == 4.0): + if (gamma == 0.0): + A = 39.42152817 + B = 0.58827911 + C = 4.05363252 + elif (gamma == 0.5): + A = 2.06550731 + B = 0.38713431 + C = 3.63592324 + elif (gamma == 1.0): + A = 13.98955016 + B = 0.72217703 + C = 5.00880454 + elif (gamma == 1.5): + A = 83.39811168 + B = 1.21616656 + C = 6.62016776 + else: + error_status = 1 + else: + error_status = 1 + elif (alpha == 2.0): + if (beta == 3.0): + if (gamma == 0.0): + A = -1.02415147e+00 + B = 1.64345425e-03 + C = 2.41297704e+00 + else: + error_status = 1 + elif (beta == 4.0): + if (gamma == 0.0): + A = 3.27271231e-06 + B = -1.94602678e+00 + C = 5.84716369e+00 + else: + error_status = 1 + # throw a warning message if the x is too small, where the + # fitting function might not work well. + if (np.min(x) < 0.04): + warnings.warn( + 'for [alpha=2, beta=4, gamma=0], the fitting function might not work reliably when the bound mass fraction is below 0.04') + else: + error_status = 1 + else: + error_status = 1 + elif (delta == 3.0): + if (alpha == 1.0): + if (beta == 3.0): + if (gamma == 0.0): + A = -0.93094144 + B = 0.04703346 + C = 0.76842857 + elif (gamma == 0.5): + A = -0.2290658 + B = 0.41234666 + C = 1.39936438 + elif (gamma == 1.0): + A = 0.90933265 + B = 0.63676731 + C = 2.18508243 + elif (gamma == 1.5): + A = 0.83533898 + B = 0.73399035 + C = 2.43215333 + else: + error_status = 1 + elif (beta == 4.0): + if (gamma == 0.0): + A = 41.3266337 + B = 0.60818104 + C = 4.07045865 + elif (gamma == 0.5): + A = 0.42123475 + B = -0.38164340 + C = 3.62931001 + elif (gamma == 1.0): + A = 15.91886121 + B = 0.71939954 + C = 4.98186081 + elif (gamma == 1.5): + A = 104.0534991 + B = 1.24696989 + C = 6.62957129 + else: + error_status = 1 + else: + error_status = 1 + elif (alpha == 2.0): + if (beta == 3.0): + if (gamma == 0.0): + A = -9.64566747e-01 + B = -2.37053865e-03 + C = 2.51094575e+00 + else: + error_status = 1 + elif (beta == 4.0): + if (gamma == 0.0): + A = 5.79321402e-07 + B = -1.62255246e+00 + C = 5.47069835e+00 + else: + error_status = 1 + # throw a warning message if the x is too small, where the + # fitting function might not work well. + if (np.min(x) < 0.04): + warnings.warn( + 'for [alpha=2, beta=4, gamma=0], the fitting function might not work reliably when the bound mass fraction is below 0.04') + else: + error_status = 1 + else: + error_status = 1 + else: + raise Exception('[delta] value not supported yet, set [delta] = 2 or 3') + + if (error_status): + raise Exception('unsupported combination of [alpha, beta, gamma]') + + return (1.0 + A) * x ** B / (1.0 + A * x ** (2.0 * B)) / np.exp(C * (1.0 - x)) + + +############################################################################### + +def _f_t_Du_2024(x, alpha=1.0, beta=3.0, gamma=1.0, delta=2.0): + """ + Return the density normalization of a subhalo given the bound mass fraction. + :param x: bound mass fraction + :param alpha: alpha parameter in the generalized NFW profile + :param beta: beta parameter in the generalized NFW profile + :param gamma: gamma parameter in the generalized NFW profile + :param delta: power index of the density truncation + :return: the density normalization + """ + error_status = 0 + + if (delta == 2.0): + if (alpha == 1.0): + if (beta == 3.0): + if (gamma == 0.0): + D = 1.48242247 + E = 0.63639230 + elif (gamma == 0.5): + D = 1.17122717 + E = 0.45000174 + elif (gamma == 1.0): + D = 0.75826635 + E = 0.23376409 + elif (gamma == 1.5): + D = -9.95103536e-01 + E = 1.31454694e-04 + else: + error_status = 1 + elif (beta == 4.0): + if (gamma == 0.0): + D = 1.18519061 + E = 0.73729078 + elif (gamma == 0.5): + D = 0.75235790 + E = 0.47771706 + elif (gamma == 1.0): + D = 0.36838147 + E = 0.22307797 + elif (gamma == 1.5): + D = 0.28500656 + E = 0.10528364 + else: + error_status = 1 + else: + error_status = 1 + elif (alpha == 2.0): + if (beta == 3.0): + if (gamma == 0.0): + D = 1.37068893 + E = 0.83149579 + else: + error_status = 1 + elif (beta == 4.0): + if (gamma == 0.0): + D = 0.83836120 + E = 1.23071929 + else: + error_status = 1 + else: + error_status = 1 + else: + error_status = 1 + elif (delta == 3.0): + if (alpha == 1.0): + if (beta == 3.0): + if (gamma == 0.0): + D = 1.40209850 + E = 0.63252747 + elif (gamma == 0.5): + D = 1.08663071 + E = 0.45228951 + elif (gamma == 1.0): + D = 1.43637351 + E = -0.24907289 + elif (gamma == 1.5): + D = 0.08092773 + E = 0.08490512 + else: + error_status = 1 + elif (beta == 4.0): + if (gamma == 0.0): + D = 1.08789526 + E = 0.72799030 + elif (gamma == 0.5): + D = 0.68794863 + E = 0.48619566 + elif (gamma == 1.0): + D = 0.33591171 + E = 0.25084764 + elif (gamma == 1.5): + D = 0.15135479 + E = 0.12838997 + else: + error_status = 1 + else: + error_status = 1 + elif (alpha == 2.0): + if (beta == 3.0): + if (gamma == 0.0): + D = 1.31446961 + E = 0.82402378 + else: + error_status = 1 + elif (beta == 4.0): + if (gamma == 0.0): + D = 0.77798730 + E = 1.21267992 + else: + error_status = 1 + else: + error_status = 1 + else: + error_status = 1 + else: + raise Exception('[delta] value not supported yet, set [delta] = 2 or 3') + + if (error_status): + raise Exception('unsupported combination of [alpha, beta, gamma]') + + # For x > 1, set f_t = 1. + f_t = np.where(x <= 1.0, (1.0 + D) * x ** E / (1.0 + D * x ** (2.0 * E)), 1.0) + + # f_t should not be larger than 1. + f_t = np.where(f_t <= 1.0, f_t, 1.0) + + return f_t + + +############################################################################### +# +# Note that the fittting functions above are derived in terms of the bound mass fraction +# +# x = M_{bound} / M_{bound,0}, +# +# and a fixed halo concentration of c_vir = 20.6. +# For a different mass definition or halo concentration, the tracks of r_te(x) and f_t(x) +# might be different. So we redefine r_te and f_t in terms of the mass ratio +# +# y = M_{bound} / M_{mx,0}, +# +# where M_{mx,0} is the enclosed mass within the infall Rmax (the radius where the circular +# velocity reaches its maximum). +# +def _M_enclosed_generalized_NFW_dimensionless(x, alpha=1.0, beta=3.0, gamma=1.0): + """ + Return the dimensionless enclosed mass within radius x = r / r_s for the generalized NFW profile. + :param x: dimensionless radius + :param alpha: alpha parameter in the generalized NFW profile + :param beta: beta parameter in the generalized NFW profile + :param gamma: gamma parameter in the generalized NFW profile + :return: the dimensionless enclosed mass + """ + status = 1 + + # For special cases, use simple analytic expressions if available. + if (alpha == 1.0): + if (beta == 3.0): + if (gamma == 0.0): + mass = 4.0 \ + * np.pi \ + * ( \ + +np.log(+1.0 + x) \ + - x \ + * (+2.0 + 3.0 * x) \ + / 2.0 \ + / (+1.0 + x) ** 2 \ + ) + elif (gamma == 0.5): + mass = +( \ + -8.0 \ + * np.pi \ + * np.sqrt(x) \ + * (3.0 + 4.0 * x) \ + ) \ + / ( \ + +3.0 \ + * (1.0 + x) ** 1.5 \ + ) \ + + 8.0 \ + * np.pi \ + * np.arcsinh(np.sqrt(x)) + elif (gamma == 1.0): + mass = +4.0 \ + * np.pi \ + * ( \ + +np.log(+1.0 + x) \ + - x \ + / (+1.0 + x) \ + ) + elif (gamma == 1.5): + mass = +8.0 \ + * np.pi \ + * ( \ + - np.sqrt(x / (1.0 + x)) \ + + np.arcsinh(np.sqrt(x)) \ + ) + else: + status = 0 + elif (beta == 4.0): + mass = 4.0 * np.pi / (3.0 - gamma) * (x / (1.0 + x)) ** (3.0 - gamma) + else: + status = 0 + elif (alpha == 2.0): + if (beta == 3.0): + if (gamma == 0.0): + mass = -4.0 \ + * np.pi \ + * ( \ + +x \ + / np.sqrt(1.0 + x ** 2) \ + + np.log(np.sqrt(1.0 + x ** 2) - x) \ + ) + elif (gamma == 1.0): + mass = 2.0 * np.pi * np.log(1.0 + x ** 2) + else: + status = 0 + elif (beta == 4.0): + if (gamma == 0.0): + mass = +2.0 \ + * np.pi \ + * ( \ + +np.arctan(x) \ + - x / (1.0 + x ** 2) + ) + elif (gamma == 1.0): + mass = 4.0 * np.pi * (1.0 - 1.0 / np.sqrt(1.0 + x ** 2)) + else: + status = 0 + else: + status = 0 + else: + status = 0 + + if (status == 0): + # Use the general expression. + mass = +4.0 \ + * np.pi \ + * x ** (3.0 - gamma) \ + * hyp2f1((3.0 - gamma) / alpha, (beta - gamma) / alpha, 1.0 + (3.0 - gamma) / alpha, -x ** alpha) \ + / (3.0 - gamma) + + return mass + + +def _R_max_generalized_NFW_dimensionless(alpha=1.0, beta=3.0, gamma=1.0): + """ + Return the dimensionless Rmax for the generalized NFW profile. + :param x: dimensionless radius + :param alpha: alpha parameter in the generalized NFW profile + :param beta: beta parameter in the generalized NFW profile + :param gamma: gamma parameter in the generalized NFW profile + :return: the dimensionless Rmax + """ + status = 1 + + # Use expressions for special cases if available. + if (alpha == 1.0): + if (beta == 3.0): + if (gamma < 2.0): + if (gamma == 0.0): + Rmax = 4.424700646872270 + elif (gamma == 0.5): + Rmax = 3.289276561384110 + elif (gamma == 1.0): + Rmax = 2.162581587064612 + elif (gamma == 1.5): + Rmax = 1.054966571869124 + else: + status = 0 + else: + raise Exception('Rmax does not exist') + elif (beta == 4.0): + if (gamma < 2.0): + Rmax = 2.0 - gamma + else: + raise Exception('Rmax does not exist') + else: + status = 0 + elif (alpha == 2.0): + if (beta == 3.0): + if (gamma < 2.0): + if (gamma == 0.0): + Rmx = 2.919847688299723 + else: + status = 0 + else: + raise Exception('Rmax does not exist') + elif (beta == 4.0): + if (gamma < 2.0): + if (gamma == 0.0): + Rmx = 1.825255642519694 + else: + status = 0 + else: + raise Exception('Rmax does not exist') + else: + status = 0 + else: + status = 0 + + if (status == 0): + # Use numerical solver. + def root(x): + return (1.0 + x ** alpha) ** ((gamma - beta) / alpha) \ + - hyp2f1((3.0 - gamma) / alpha, (beta - gamma) / alpha, 1.0 + (3.0 - gamma) / alpha, -x ** alpha) \ + / (3.0 - gamma) + + Rmax = newton(root, 1.0) + + return Rmax + + +def _M_total_generalized_NFW_Truncated_dimensionless(tau, alpha=1.0, beta=3.0, gamma=1.0, delta=2.0): + """ + Return the dimensionless total mass of a subhalo with truncated generalized NFW profile: + rho(x) = rho_org(x) / (1.0 + (x / tau)**2), + where rho_org(x) is the original density profile (without truncation), x = r / r_s, and tau = r_te / r_s. + :param tau: dimensionless truncation radius + :param alpha: alpha parameter in the generalized NFW profile + :param beta: beta parameter in the generalized NFW profile + :param gamma: gamma parameter in the generalized NFW profile + :param delta: power index of the density truncation + :return: the dimensionless enclosed mass + """ + status = 1 + + if (delta == 2.0): + if (alpha == 1.0): + if (beta == 3.0): + if (gamma == 0.0): + mass = 4.0 * np.pi \ + * tau ** 2 / (tau ** 2 + 1.0) ** 3 \ + * ( \ + 0.5 * (-1.0 + (np.pi - tau) * tau) * (3.0 * tau ** 2 - 1.0) \ + + tau ** 2 * (tau ** 2 - 3.0) * np.log(tau) \ + ) + elif (gamma == 0.5): + atan_tau = np.arctan(tau) + mass = 4.0 * np.pi \ + * tau ** 2 / (tau ** 2 + 1.0) ** (9.0 / 4.0) \ + / 6.0 \ + * ( \ + 3.0 * np.sqrt(2.0) * np.pi * np.sqrt(tau) * (tau ** 2 + 1.0) \ + * (np.sin(2.5 * atan_tau) - np.cos(2.5 * atan_tau)) \ + + 2.0 * np.sqrt(2.0) \ + * (np.sqrt(1.0 + np.sqrt(2.0 * tau) + tau) - np.sqrt(1.0 - np.sqrt(2.0 * tau) + tau)) \ + * np.sqrt(tau * (tau ** 2 + 1.0) * (+1.0 - tau + np.sqrt(tau ** 2 + 1.0))) \ + * np.cos(0.25 * (np.pi + 2.0 * atan_tau)) \ + - 2.0 * np.sqrt(2.0) \ + * (np.sqrt(1.0 + np.sqrt(2.0 * tau) + tau) + np.sqrt(1.0 - np.sqrt(2.0 * tau) + tau)) \ + * np.sqrt(tau * (tau ** 2 + 1.0) * (-1.0 + tau + np.sqrt(tau ** 2 + 1.0))) \ + * np.sin(0.25 * (np.pi + 2.0 * atan_tau)) \ + + ( \ + +( \ + +2.0 * np.sqrt(2.0) + 4.0 * np.sqrt(tau) + 3.0 * np.sqrt(2.0) * tau \ + - 2.0 * np.sqrt(2.0) * tau ** 2 - 4.0 * tau ** 2.5 - np.sqrt( + 2.0) * tau ** 3 \ + ) * np.sqrt(1.0 - np.sqrt(2.0 * tau) + tau) \ + + ( \ + +2.0 * np.sqrt(2.0) - 4.0 * np.sqrt(tau) + 3.0 * np.sqrt(2.0) * tau \ + - 2.0 * np.sqrt(2.0) * tau ** 2 + 4.0 * tau ** 2.5 - np.sqrt( + 2.0) * tau ** 3 \ + ) * np.sqrt(1.0 + np.sqrt(2.0 * tau) + tau) \ + ) \ + * np.sqrt((-1.0 + tau + np.sqrt(tau ** 2 + 1.0)) / tau) \ + * np.cos(0.25 * (np.pi + 2.0 * atan_tau)) \ + + ( \ + +( \ + +2.0 * np.sqrt(2.0) + 4.0 * np.sqrt(tau) + 3.0 * np.sqrt(2.0) * tau \ + - 2.0 * np.sqrt(2.0) * tau ** 2 - 4.0 * tau ** 2.5 - np.sqrt( + 2.0) * tau ** 3 \ + ) * np.sqrt(1.0 - np.sqrt(2.0 * tau) + tau) \ + - ( \ + +2.0 * np.sqrt(2.0) - 4.0 * np.sqrt(tau) + 3.0 * np.sqrt(2.0) * tau \ + - 2.0 * np.sqrt(2.0) * tau ** 2 + 4.0 * tau ** 2.5 - np.sqrt( + 2.0) * tau ** 3 \ + ) * np.sqrt(1.0 + np.sqrt(2.0 * tau) + tau) \ + ) \ + * np.sqrt((+1.0 - tau + np.sqrt(tau ** 2 + 1.0)) / tau) \ + * np.sin(0.25 * (np.pi + 2.0 * atan_tau)) \ + - 12.0 \ + * ( \ + +(tau ** 2 - 1.0) * np.arcsin( + np.sqrt(0.5 * (1.0 + tau - np.sqrt(tau ** 2 + 1.0)))) \ + + tau \ + * np.log(tau + np.sqrt(tau ** 2 + 1.0) + np.sqrt( + 2.0 * tau * (tau + np.sqrt(tau ** 2 + 1.0)))) \ + ) \ + * np.sqrt(tau) \ + * np.cos(0.25 * (np.pi + 2.0 * atan_tau)) \ + - 12.0 \ + * ( \ + +2.0 * tau * np.arcsin(np.sqrt(0.5 * (1.0 + tau - np.sqrt(tau ** 2 + 1.0)))) \ + - 0.5 * (tau ** 2 - 1.0) \ + * np.log(tau + np.sqrt(tau ** 2 + 1.0) + np.sqrt( + 2.0 * tau * (tau + np.sqrt(tau ** 2 + 1.0)))) \ + ) \ + * np.sqrt(tau) \ + * np.sin(0.25 * (np.pi + 2.0 * atan_tau)) \ + ) + elif (gamma == 1.0): + mass = 4.0 * np.pi \ + * tau ** 2 / (tau ** 2 + 1.0) ** 2 \ + * ( \ + (tau ** 2 - 1.0) * np.log(tau) \ + + tau * np.pi \ + - (tau ** 2 + 1.0) + ) + elif (gamma == 1.5): + atan_tau = np.arctan(tau) + mass = 4.0 * np.pi \ + * tau ** 1.5 / (tau ** 2 + 1.0) ** (5.0 / 4.0) \ + / 2.0 \ + * ( \ + +np.sqrt(2.0) * np.pi * np.sqrt(tau ** 2 + 1.0) \ + * (np.cos(1.5 * atan_tau) + np.sin(1.5 * atan_tau)) \ + - np.sqrt(2.0) \ + * ( \ + +np.sqrt(1.0 + np.sqrt(2.0 * tau) + tau) \ + + np.sqrt(1.0 - np.sqrt(2.0 * tau) + tau) \ + ) \ + * np.sqrt((tau ** 2 + 1.0) * (-1.0 + tau + np.sqrt(tau ** 2 + 1.0))) \ + * np.cos(0.25 * (np.pi + 2.0 * atan_tau)) \ + - np.sqrt(2.0) \ + * ( \ + +np.sqrt(1.0 + np.sqrt(2.0 * tau) + tau) \ + - np.sqrt(1.0 - np.sqrt(2.0 * tau) + tau) \ + ) \ + * np.sqrt((tau ** 2 + 1.0) * (+1.0 - tau + np.sqrt(tau ** 2 + 1.0))) \ + * np.sin(0.25 * (np.pi + 2.0 * atan_tau)) \ + - 2.0 \ + * ( \ + +2.0 * tau * np.arcsin(np.sqrt(0.5 * (1.0 + tau - np.sqrt(tau ** 2 + 1.0)))) \ + + np.log(tau + np.sqrt(tau ** 2 + 1.0) + np.sqrt( + 2.0 * tau * (tau + np.sqrt(tau ** 2 + 1.0)))) \ + ) \ + * np.cos(0.25 * (np.pi + 2.0 * atan_tau)) \ + - 2.0 \ + * ( \ + +2.0 * np.arcsin(np.sqrt(0.5 * (1.0 + tau - np.sqrt(tau ** 2 + 1.0)))) \ + - tau \ + * np.log(tau + np.sqrt(tau ** 2 + 1.0) + np.sqrt( + 2.0 * tau * (tau + np.sqrt(tau ** 2 + 1.0)))) \ + ) \ + * np.sin(0.25 * (np.pi + 2.0 * atan_tau)) \ + ) + else: + status = 0 + else: + status = 0 + elif (alpha == 2.0): + if (beta == 4.0): + if (gamma == 0.0): + mass = np.pi ** 2 * tau ** 2 / (tau + 1.0) ** 2 + else: + status = 0 + else: + status = 0 + else: + status = 0 + elif (delta == 3.0): + if (alpha == 1.0 and beta == 3.0 and gamma == 1.0): + mass = 4.0 * np.pi \ + * tau ** 2 / (tau ** 3 - 1.0) ** 2 \ + * ( \ + -tau * (tau ** 3 - 1.0) \ + + 2.0 * np.sqrt(3.0) * np.pi / 9.0 \ + * (tau - 1.0) ** 2 \ + * (2.0 * tau + 1.0) \ + + tau \ + * (tau ** 3 + 2.0) \ + * np.log(tau) \ + ) + else: + status = 0 + else: + status = 0 + + if (status == 0): + # Use numerical integration. + def integral(x): + return 4.0 * np.pi * x ** 2 \ + / x ** gamma \ + / (1.0 + x ** alpha) ** ((beta - gamma) / alpha) \ + / (1.0 + (x / tau) ** delta) + + mass = quad(integral, 0.0, np.inf) + + return mass + + +# Mass table used to compute the radius enclosing a certain mass for truncated generalized NFW profile. +N_per_decade = 10 +x_min = 1.0e-6 +x_max = 1.0e+4 +N_tab = int(np.log10(x_max / x_min) * N_per_decade) +x_tab = np.geomspace(x_min, x_max, N_tab) +M_tab = np.zeros(N_tab) + +Mass_table_Initialized = False +alpha_used = None +beta_used = None +gamma_used = None +delta_used = None +M_min = np.inf +M_max = 0.0 + + +def _Truncation_Radius(M, f_t, alpha=1.0, beta=3.0, gamma=1.0, delta=2.0): + """ + Return the truncation radius of a subhalo given the bound mass and normalization of density profile. + The density profile is assumed to be + rho(x) = rho_org(x) * f_t / (1.0 + (x / tau)**2), + where rho_org(x) is the original density profile (without truncation), x = r / r_s, and tau = r_te / r_s. + :param M: dimensionless mass of the subhalo + :param f_t: normalization of the density profile + :param alpha: alpha parameter in the generalized NFW profile + :param beta: beta parameter in the generalized NFW profile + :param gamma: gamma parameter in the generalized NFW profile + :param delta: power index of the density truncation + :return: the dimensionless truncation radius + """ + + global N_per_decade, x_min, x_max, N_tab, x_tab, M_tab + global Mass_table_Initialized, alpha_used, beta_used, gamma_used, delta_used, M_min, M_max + + M_unnormalized = M / f_t + + retabulate = True + if (Mass_table_Initialized and alpha_used == alpha and beta_used == beta and gamma_used == gamma \ + and delta_used == delta and np.min(M_unnormalized) >= M_min and np.max(M_unnormalized) <= M_max): + retabulate = False + + if (retabulate): + M_min = _M_total_generalized_NFW_Truncated_dimensionless(x_min, alpha, beta, gamma, delta) + M_max = _M_total_generalized_NFW_Truncated_dimensionless(x_max, alpha, beta, gamma, delta) + while (np.min(M_unnormalized) < M_min): + x_min = x_min / 2.0 + M_min = _M_total_generalized_NFW_Truncated_dimensionless(x_min, alpha, beta, gamma, delta) + + while (np.max(M_unnormalized) > M_max): + x_max = x_max * 2.0 + M_max = _M_total_generalized_NFW_Truncated_dimensionless(x_max, alpha, beta, gamma, delta) + + N_tab = int(np.log10(x_max / x_min) * N_per_decade) + x_tab = np.geomspace(x_min, x_max, N_tab) + M_tab = _M_total_generalized_NFW_Truncated_dimensionless(x_tab, alpha, beta, gamma, delta) + + Mass_table_Initialized = True + alpha_used = alpha + beta_used = beta + gamma_used = gamma + delta_used = delta + + x_M_interp_loglog = interp1d(np.log(M_tab), np.log(x_tab), kind='cubic') + + return np.exp(x_M_interp_loglog(np.log(M_unnormalized))) + + +def Convert_to_reference_model(alpha=1.0, beta=3.0, gamma=1.0): + """ + Compute the conversion factors from parameter y = M_{bound} / M_{mx,0} to parameter + x = M_{bound} / M_{bound,0} in the reference model. + :param alpha: alpha parameter in the generalized NFW profile + :param beta: beta parameter in the generalized NFW profile + :param gamma: gamma parameter in the generalized NFW profile + :return: conversion factors [Mmx_ref / Mbound_ref_0, Rmax_ref / Rvir_ref] + """ + error_status = 0 + + if (alpha == 1.0): + if (beta == 3.0): + if (gamma == 0.0): + Mbound_ref_0 = 1.1969934e9 # in units of Msun + elif (gamma == 0.5): + Mbound_ref_0 = 1.1793903e9 # in units of Msun + elif (gamma == 1.0): + Mbound_ref_0 = 1.1605002e9 # in units of Msun + elif (gamma == 1.5): + Mbound_ref_0 = 1.1395515e9 # in units of Msun + else: + error_status = 1 + elif (beta == 4.0): + if (gamma == 0.0): + Mbound_ref_0 = 1.0426536e9 # in units of Msun + elif (gamma == 0.5): + Mbound_ref_0 = 1.0353760e9 # in units of Msun + elif (gamma == 1.0): + Mbound_ref_0 = 1.0281669e9 # in units of Msun + elif (gamma == 1.5): + Mbound_ref_0 = 1.0210254e9 # in units of Msun + else: + error_status = 1 + else: + error_status = 1 + elif (alpha == 2.0): + if (beta == 3.0): + if (gamma == 0.0): + Mbound_ref_0 = 1.1343262e9 # in units of Msun + else: + error_status = 1 + elif (beta == 4.0): + if (gamma == 0.0): + Mbound_ref_0 = 1.0194281e9 # in units of Msun + else: + error_status = 1 + else: + error_status = 1 + else: + error_status = 1 + + if (error_status): + raise Exception('unsupported combination of [alpha, beta, gamma]') + + Mvir_ref = 1.0e9 # in units of Msun, note that Mbound_ref_0 is slightly larger than Mvir_ref + Rvir_ref = 0.02632439488393568 # in units of Mpc + c_ref = 20.587820037799197 + Rmax_dimensionless = _R_max_generalized_NFW_dimensionless(alpha, beta, gamma) + Rmax_ref = Rmax_dimensionless / c_ref * Rvir_ref + Mmx_ref = _M_enclosed_generalized_NFW_dimensionless(Rmax_dimensionless, alpha, beta, gamma) \ + / _M_enclosed_generalized_NFW_dimensionless(c_ref, alpha, beta, gamma) \ + * Mvir_ref + + return Mmx_ref / Mbound_ref_0, Rmax_ref / Rvir_ref + + +def compute_r_te_and_f_t(Mbound, MvirInfall, RvirInfall, cInfall, alpha=1.0, beta=3.0, gamma=1.0, delta=2.0): + """ + Return the effective tidal (truncation) radius and the density normalization of a subhalo given the + bound mass and infall halo properties. + :param Mbound: bound mass + :param MvirInfall: virial mass at infall + :param RvirInfall: virial radius at infall + :param cInfall: halo concentration at infall + :param alpha: alpha parameter in the generalized NFW profile + :param beta: beta parameter in the generalized NFW profile + :param gamma: gamma parameter in the generalized NFW profile + :param delta: power index of the density truncation + :return: the effective tidal (truncation) radius + """ + + Rs = RvirInfall / cInfall + Rmax_dimensionless = _R_max_generalized_NFW_dimensionless(alpha, beta, gamma) + + # Characteristic mass. + M0 = MvirInfall / _M_enclosed_generalized_NFW_dimensionless(cInfall, alpha, beta, gamma) + + M_mx = M0 * _M_enclosed_generalized_NFW_dimensionless(Rmax_dimensionless, alpha, beta, gamma) + + y = Mbound / M_mx + + # Convert y to parameter x = MBoud / MvirInfall in the reference model. + y_scale, R_scale = Convert_to_reference_model(alpha, beta, gamma) + + # r_te = r_te_Du_2024(y * y_scale, alpha, beta, gamma, delta) / R_scale * Rmax + ## Make sure that the truncation radius is smaller than the virial radius at infall. + # r_te = np.where(r_te < RvirInfall, r_te, RvirInfall) + + f_t = _f_t_Du_2024(y * y_scale, alpha, beta, gamma, delta) + + # Instead of directly using the fitting function for r_te, we first compute f_t + # using the fitting function and then compute r_te assuming that the total mass + # of the truncated halo equals to the bound mass. r_te computed in this way is + # in good agreement with the fitting function, but is more reliable where the + # fitting function is not calibrated to simulations. + Mbound_dimensionless = Mbound / M0 + r_te = Rs * _Truncation_Radius(Mbound_dimensionless, f_t, alpha, beta, gamma) + + return r_te, f_t diff --git a/pyHalo/Halos/halo_base.py b/pyHalo/Halos/halo_base.py index 67de2ad..320e763 100644 --- a/pyHalo/Halos/halo_base.py +++ b/pyHalo/Halos/halo_base.py @@ -66,6 +66,8 @@ def rescale_normalization(self, factor): delattr(self, '_params_physical') if hasattr(self, '_kwargs_lenstronomy'): delattr(self, '_kwargs_lenstronomy') + if hasattr(self, '_nfw_params'): + delattr(self, '_nfw_params') @property @abstractmethod diff --git a/pyHalo/Halos/tidal_truncation.py b/pyHalo/Halos/tidal_truncation.py index 7770d1b..0689deb 100644 --- a/pyHalo/Halos/tidal_truncation.py +++ b/pyHalo/Halos/tidal_truncation.py @@ -4,6 +4,8 @@ from colossus.lss import peaks from colossus.halo import splashback from pyHalo.Halos.galacticus_truncation.interp_mass_loss import InterpGalacticus, InterpGalacticusKeeley24 +from pyHalo.Halos.galacticus_truncation.transfer_function_density_profile import compute_r_te_and_f_t + class TruncationSplashBack(object): @@ -349,8 +351,19 @@ def truncation_radius_halo(self, halo): :return: the truncation radius in physical kpc """ - return self.truncation_radius(halo.mass, halo.c, - halo.time_since_infall, halo.z) + halo_mass = halo.mass + infall_concentration = halo.c + time_since_infall = halo.time_since_infall + + log10c = np.log10(infall_concentration) + log10mbound_over_minfall = self._mass_loss_interp(log10c, time_since_infall) + m_bound = halo_mass * 10 ** log10mbound_over_minfall + _, rs, r200 = self._lens_cosmo.NFW_params_physical(halo_mass, + infall_concentration, + halo.z) + r_te, f_t = compute_r_te_and_f_t(m_bound, halo_mass, r200, infall_concentration) + halo.rescale_normalization(f_t) + return r_te def truncation_radius(self, halo_mass, infall_concentration, time_since_infall, z_eval_angles): @@ -365,13 +378,12 @@ def truncation_radius(self, halo_mass, infall_concentration, """ log10c = np.log10(infall_concentration) log10mbound_over_minfall = self._mass_loss_interp(log10c, time_since_infall) - point = self._make_params_in_bounds_tau_evaluate((log10c, log10mbound_over_minfall)) - log10tau = float(self._tau_mf_interpolation(point)) - tau = 10 ** log10tau - _, rs, _ = self._lens_cosmo.NFW_params_physical(halo_mass, + m_bound = halo_mass * 10 ** log10mbound_over_minfall + _, rs, r200 = self._lens_cosmo.NFW_params_physical(halo_mass, infall_concentration, z_eval_angles) - return tau * rs + r_te, _ = compute_r_te_and_f_t(m_bound, halo_mass, r200, infall_concentration) + return r_te diff --git a/tests/test_halos/test_galacticus_truncation.py b/tests/test_halos/test_galacticus_truncation.py index 97f2929..2991c24 100644 --- a/tests/test_halos/test_galacticus_truncation.py +++ b/tests/test_halos/test_galacticus_truncation.py @@ -1,6 +1,8 @@ import pytest import numpy.testing as npt from pyHalo.Halos.tidal_truncation import TruncationGalacticus, TruncationGalacticusKeeley24 +from pyHalo.Halos.HaloModels.TNFW import TNFWSubhalo +from pyHalo.Halos.concentration import ConcentrationDiemerJoyce import numpy as np from pyHalo.Halos.lens_cosmo import LensCosmo @@ -101,11 +103,36 @@ def __init__(self): self.c = 16.0 self.time_since_infall = 5.0 self.z = 0.5 + def rescale_normalization(self, rescale): + return halo = DummyHalo() rt_kpc = self.tg.truncation_radius_halo(halo) npt.assert_equal(True, np.isfinite(rt_kpc)) npt.assert_equal(True, isinstance(rt_kpc, float)) + def test_tg(self): + lens_cosmo = LensCosmo(0.5, 2.0,None) + truncation_class = TruncationGalacticus(lens_cosmo) + concentration_class = ConcentrationDiemerJoyce(lens_cosmo.cosmo, scatter=False) + mass = 10 ** 8 + x = 0 + y = 0 + r3d = None + z = 0.5 + sub_flag = True + halo_args = {} + unique_tag = 1 + tnfw_profile_1 = TNFWSubhalo(mass, x, y, r3d, z, sub_flag, + lens_cosmo, halo_args, truncation_class, concentration_class, unique_tag) + tnfw_profile_1._c = 3.0 + tnfw_profile_1._time_since_infall = 6.0 + + m_bound = tnfw_profile_1.bound_mass + rescale_norm = tnfw_profile_1._rescale_norm + kwargs, _ = tnfw_profile_1.lenstronomy_params + npt.assert_equal(rescale_norm<1, True) + npt.assert_equal(m_bound < mass, True) + if __name__ == '__main__': pytest.main() From 3214bb294ed836ab52d02dae66499e2733e5eed1 Mon Sep 17 00:00:00 2001 From: Daniel Gilman Date: Thu, 5 Sep 2024 01:04:20 -0500 Subject: [PATCH 2/3] implement host halo concentration dependence --- .../custom_mass_concentration_relations.ipynb | 49 +++++++++++++++++++ pyHalo/Halos/HaloModels/TNFW.py | 21 +++----- .../galacticus_truncation/interp_mass_loss.py | 19 ++++--- .../galacticus_truncation/johnsonSUparams.py | 5 +- pyHalo/Halos/tidal_truncation.py | 12 ++--- pyHalo/PresetModels/cdm.py | 8 +-- pyHalo/PresetModels/wdm.py | 17 ++++--- .../test_halos/test_galacticus_truncation.py | 4 +- tests/test_halos/test_tnfw_halo.py | 2 +- tests/test_halos/test_truncation.py | 2 +- 10 files changed, 96 insertions(+), 43 deletions(-) diff --git a/example_notebooks/custom_mass_concentration_relations.ipynb b/example_notebooks/custom_mass_concentration_relations.ipynb index 0eea95c..0454801 100644 --- a/example_notebooks/custom_mass_concentration_relations.ipynb +++ b/example_notebooks/custom_mass_concentration_relations.ipynb @@ -292,6 +292,55 @@ "ax2.set_ylim(0.9, 1.4)\n" ] }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "from pyHalo.Halos.lens_cosmo import InfallDistributionGalacticus2024 \n", + "from pyHalo.Cosmology.cosmology import Cosmology" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGwCAYAAAC3qV8qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAApXElEQVR4nO3df1TVdZ7H8dcF5EcIl9CVyx1BqPUXZuqIIeo008oRk6m1nMpZUprc3DWoEDVxCyszUacyKZPR7ahnRtemc7IpPVEsmk6JSBilpuY6JjYGNGNyU1ck+O4fHe92k9L0i5cPPB/n3HPi+/3e731/sy5Pv/fH12FZliUAAACDBPh7AAAAgB+LgAEAAMYhYAAAgHEIGAAAYBwCBgAAGIeAAQAAxiFgAACAcYL8PUBbaWlp0bFjxxQRESGHw+HvcQAAwEWwLEtfffWV3G63AgK+/zxLhw2YY8eOKS4uzt9jAACAS3D06FH17Nnze9d32ICJiIiQ9M2/gMjISD9PAwAALobH41FcXJz39/j36bABc+5lo8jISAIGAADDXOjtH7yJFwAAGIeAAQAAxiFgAACAcQgYAABgHAIGAAAYh4ABAADGIWAAAIBxCBgAAGAcAgYAABiHgAEAAMYhYAAAgHEIGAAAYBwCBgAAGIeAAQAAxiFgAACAcYL8PQDal4T8Tbbs59OFGbbsBwCA1nAGBgAAGIeAAQAAxiFgAACAcQgYAABgHAIGAAAYh4ABAADGIWAAAIBxCBgAAGCcHx0w27Zt0y233CK32y2Hw6HXXnvNZ71lWZo7d65iY2MVFhamtLQ0HTx40Geb48ePKzMzU5GRkYqKitKUKVN08uRJn20++ugj/exnP1NoaKji4uK0ePHiH390AACgQ/rRAXPq1CkNGjRIy5Yta3X94sWLVVRUpOLiYlVUVCg8PFzp6ek6c+aMd5vMzEzt3btXpaWl2rhxo7Zt26apU6d613s8Ho0ZM0a9evVSVVWVfvvb3+rxxx/XihUrLuEQAQBAR+OwLMu65Ds7HNqwYYPGjx8v6ZuzL263WzNmzNDMmTMlSQ0NDYqJidHq1as1ceJE7du3T0lJSaqsrFRycrIkqaSkROPGjdNnn30mt9ut5cuX65FHHlFtba2Cg4MlSfn5+Xrttde0f//+i5rN4/HI6XSqoaFBkZGRl3qInY5dlxKwC5ckAIDO5WJ/f9v6HpjDhw+rtrZWaWlp3mVOp1MpKSkqLy+XJJWXlysqKsobL5KUlpamgIAAVVRUeLe58cYbvfEiSenp6Tpw4IC+/PLLVh+7sbFRHo/H5wYAADomWwOmtrZWkhQTE+OzPCYmxruutrZWPXr08FkfFBSk6Ohon21a28e3H+O7CgsL5XQ6vbe4uLjLPyAAANAudZhPIc2ZM0cNDQ3e29GjR/09EgAAaCO2BozL5ZIk1dXV+Syvq6vzrnO5XKqvr/dZ//XXX+v48eM+27S2j28/xneFhIQoMjLS5wYAADomWwMmMTFRLpdLZWVl3mUej0cVFRVKTU2VJKWmpurEiROqqqrybrN582a1tLQoJSXFu822bdvU1NTk3aa0tFR9+/bV1VdfbefIAADAQD86YE6ePKnq6mpVV1dL+uaNu9XV1aqpqZHD4VBubq7mz5+v119/Xbt379bkyZPldru9n1Tq37+/xo4dq/vuu087d+7Ue++9p5ycHE2cOFFut1uS9C//8i8KDg7WlClTtHfvXr388staunSp8vLybDtwAABgrqAfe4f3339fN910k/fnc1GRlZWl1atX6+GHH9apU6c0depUnThxQqNGjVJJSYlCQ0O991m7dq1ycnI0evRoBQQEaMKECSoqKvKudzqdevvtt5Wdna2hQ4eqe/fumjt3rs93xQAAgM7rsr4Hpj3je2AuDd8DAwDwJ798DwwAAMCVQMAAAADjEDAAAMA4BAwAADAOAQMAAIxDwAAAAOMQMAAAwDgEDAAAMA4BAwAAjEPAAAAA4xAwAADAOAQMAAAwDgEDAACMQ8AAAADjEDAAAMA4BAwAADAOAQMAAIxDwAAAAOMQMAAAwDgEDAAAMA4BAwAAjEPAAAAA4xAwAADAOAQMAAAwDgEDAACMQ8AAAADjEDAAAMA4BAwAADAOAQMAAIxDwAAAAOMQMAAAwDgEDAAAME6QvwcAfkhC/iZb9vPpwgxb9gMAaB84AwMAAIxDwAAAAOMQMAAAwDgEDAAAMA4BAwAAjEPAAAAA4xAwAADAOAQMAAAwDgEDAACMQ8AAAADjEDAAAMA4BAwAADAOF3PsAOy64CEAAKbgDAwAADAOAQMAAIxDwAAAAOMQMAAAwDgEDAAAMA4BAwAAjEPAAAAA4xAwAADAOAQMAAAwDgEDAACMw6UE0CnYdbmFTxdm2LIfAMDlsf0MTHNzswoKCpSYmKiwsDBde+21evLJJ2VZlncby7I0d+5cxcbGKiwsTGlpaTp48KDPfo4fP67MzExFRkYqKipKU6ZM0cmTJ+0eFwAAGMj2gFm0aJGWL1+uF154Qfv27dOiRYu0ePFiPf/8895tFi9erKKiIhUXF6uiokLh4eFKT0/XmTNnvNtkZmZq7969Ki0t1caNG7Vt2zZNnTrV7nEBAICBHNa3T43Y4Je//KViYmL00ksveZdNmDBBYWFh+sMf/iDLsuR2uzVjxgzNnDlTktTQ0KCYmBitXr1aEydO1L59+5SUlKTKykolJydLkkpKSjRu3Dh99tlncrvdF5zD4/HI6XSqoaFBkZGRdh5iu8PVqK8cXkICgLZ1sb+/bT8DM2LECJWVlemTTz6RJH344Yd69913dfPNN0uSDh8+rNraWqWlpXnv43Q6lZKSovLycklSeXm5oqKivPEiSWlpaQoICFBFRUWrj9vY2CiPx+NzAwAAHZPtb+LNz8+Xx+NRv379FBgYqObmZj311FPKzMyUJNXW1kqSYmJifO4XExPjXVdbW6sePXr4DhoUpOjoaO8231VYWKgnnnjC7sMBAADtkO1nYP74xz9q7dq1WrdunXbt2qU1a9bo6aef1po1a+x+KB9z5sxRQ0OD93b06NE2fTwAAOA/tp+BmTVrlvLz8zVx4kRJ0sCBA3XkyBEVFhYqKytLLpdLklRXV6fY2Fjv/erq6jR48GBJksvlUn19vc9+v/76ax0/ftx7/+8KCQlRSEiI3YcDAADaIdvPwJw+fVoBAb67DQwMVEtLiyQpMTFRLpdLZWVl3vUej0cVFRVKTU2VJKWmpurEiROqqqrybrN582a1tLQoJSXF7pEBAIBhbD8Dc8stt+ipp55SfHy8BgwYoA8++EDPPvus7r33XkmSw+FQbm6u5s+fr969eysxMVEFBQVyu90aP368JKl///4aO3as7rvvPhUXF6upqUk5OTmaOHHiRX0CCQAAdGy2B8zzzz+vgoIC3X///aqvr5fb7da//du/ae7cud5tHn74YZ06dUpTp07ViRMnNGrUKJWUlCg0NNS7zdq1a5WTk6PRo0crICBAEyZMUFFRkd3jAgAAA9n+PTDtBd8Dg7bA98AAQNvy2/fAAAAAtDUCBgAAGIeAAQAAxiFgAACAcQgYAABgHAIGAAAYh4ABAADGIWAAAIBxCBgAAGAcAgYAABiHgAEAAMYhYAAAgHEIGAAAYBwCBgAAGIeAAQAAxiFgAACAcQgYAABgHAIGAAAYh4ABAADGIWAAAIBxCBgAAGAcAgYAABiHgAEAAMYhYAAAgHEIGAAAYBwCBgAAGIeAAQAAxiFgAACAcQgYAABgHAIGAAAYh4ABAADGIWAAAIBxCBgAAGAcAgYAABiHgAEAAMYhYAAAgHGC/D0AYJKE/E227OfThRm27AcAOisCBvADQggALg8vIQEAAOMQMAAAwDgEDAAAMA4BAwAAjEPAAAAA4/ApJD+y65MoAAB0NpyBAQAAxiFgAACAcQgYAABgHAIGAAAYh4ABAADGIWAAAIBxCBgAAGAcAgYAABiHgAEAAMYhYAAAgHEIGAAAYBwCBgAAGIeAAQAAxmmTgPnrX/+qu+++W926dVNYWJgGDhyo999/37vesizNnTtXsbGxCgsLU1pamg4ePOizj+PHjyszM1ORkZGKiorSlClTdPLkybYYFwAAGMb2gPnyyy81cuRIdenSRW+++aY+/vhjPfPMM7r66qu92yxevFhFRUUqLi5WRUWFwsPDlZ6erjNnzni3yczM1N69e1VaWqqNGzdq27Ztmjp1qt3jAgAAAzksy7Ls3GF+fr7ee+89/fnPf251vWVZcrvdmjFjhmbOnClJamhoUExMjFavXq2JEydq3759SkpKUmVlpZKTkyVJJSUlGjdunD777DO53e7z9tvY2KjGxkbvzx6PR3FxcWpoaFBkZKSdh2ibhPxN/h4Bhvt0YYa/RwAAW3k8Hjmdzgv+/rb9DMzrr7+u5ORk3XHHHerRo4eGDBmilStXetcfPnxYtbW1SktL8y5zOp1KSUlReXm5JKm8vFxRUVHeeJGktLQ0BQQEqKKiotXHLSwslNPp9N7i4uLsPjQAANBOBNm9w7/85S9avny58vLy9B//8R+qrKzUgw8+qODgYGVlZam2tlaSFBMT43O/mJgY77ra2lr16NHDd9CgIEVHR3u3+a45c+YoLy/P+/O5MzBAR2bXWTzO5AAwje0B09LSouTkZC1YsECSNGTIEO3Zs0fFxcXKysqy++G8QkJCFBIS0mb7BwAA7YftLyHFxsYqKSnJZ1n//v1VU1MjSXK5XJKkuro6n23q6uq861wul+rr633Wf/311zp+/Lh3GwAA0HnZHjAjR47UgQMHfJZ98skn6tWrlyQpMTFRLpdLZWVl3vUej0cVFRVKTU2VJKWmpurEiROqqqrybrN582a1tLQoJSXF7pEBAIBhbH8Jafr06RoxYoQWLFigO++8Uzt37tSKFSu0YsUKSZLD4VBubq7mz5+v3r17KzExUQUFBXK73Ro/frykb87YjB07Vvfdd5+Ki4vV1NSknJwcTZw4sdVPIAEAgM7F9oAZNmyYNmzYoDlz5mjevHlKTEzUc889p8zMTO82Dz/8sE6dOqWpU6fqxIkTGjVqlEpKShQaGurdZu3atcrJydHo0aMVEBCgCRMmqKioyO5xAQCAgWz/Hpj24mI/R+5PfA8M2gs+hQSgvfDb98AAAAC0NQIGAAAYh4ABAADGIWAAAIBxCBgAAGAcAgYAABiHgAEAAMYhYAAAgHEIGAAAYBwCBgAAGIeAAQAAxiFgAACAcQgYAABgHAIGAAAYh4ABAADGIWAAAIBxCBgAAGAcAgYAABgnyN8DAPC/hPxNtuzn04UZtuwHAC6EMzAAAMA4nIG5BHb9bRUAAFwazsAAAADjEDAAAMA4BAwAADAOAQMAAIxDwAAAAOMQMAAAwDgEDAAAMA4BAwAAjEPAAAAA4/BNvABswzWVAFwpnIEBAADGIWAAAIBxCBgAAGAcAgYAABiHgAEAAMYhYAAAgHEIGAAAYBwCBgAAGIeAAQAAxiFgAACAcQgYAABgHAIGAAAYh4ABAADGIWAAAIBxCBgAAGAcAgYAABiHgAEAAMYhYAAAgHEIGAAAYBwCBgAAGIeAAQAAxiFgAACAcQgYAABgHAIGAAAYh4ABAADGIWAAAIBx2jxgFi5cKIfDodzcXO+yM2fOKDs7W926dVPXrl01YcIE1dXV+dyvpqZGGRkZuuqqq9SjRw/NmjVLX3/9dVuPCwAADNCmAVNZWanf/e53uv76632WT58+XW+88YZeeeUVbd26VceOHdPtt9/uXd/c3KyMjAydPXtW27dv15o1a7R69WrNnTu3LccFAACGaLOAOXnypDIzM7Vy5UpdffXV3uUNDQ166aWX9Oyzz+qf/umfNHToUK1atUrbt2/Xjh07JElvv/22Pv74Y/3hD3/Q4MGDdfPNN+vJJ5/UsmXLdPbs2VYfr7GxUR6Px+cGAAA6pjYLmOzsbGVkZCgtLc1neVVVlZqamnyW9+vXT/Hx8SovL5cklZeXa+DAgYqJifFuk56eLo/Ho71797b6eIWFhXI6nd5bXFxcGxwVAABoD9okYNavX69du3apsLDwvHW1tbUKDg5WVFSUz/KYmBjV1tZ6t/l2vJxbf25da+bMmaOGhgbv7ejRozYcCQAAaI+C7N7h0aNH9dBDD6m0tFShoaF27/57hYSEKCQk5Io9HgAA8B/bz8BUVVWpvr5eP/3pTxUUFKSgoCBt3bpVRUVFCgoKUkxMjM6ePasTJ0743K+urk4ul0uS5HK5zvtU0rmfz20DAAA6L9sDZvTo0dq9e7eqq6u9t+TkZGVmZnr/uUuXLiorK/Pe58CBA6qpqVFqaqokKTU1Vbt371Z9fb13m9LSUkVGRiopKcnukQEAgGFsfwkpIiJC1113nc+y8PBwdevWzbt8ypQpysvLU3R0tCIjI/XAAw8oNTVVw4cPlySNGTNGSUlJmjRpkhYvXqza2lo9+uijys7O5mUiAABgf8BcjCVLliggIEATJkxQY2Oj0tPT9eKLL3rXBwYGauPGjZo2bZpSU1MVHh6urKwszZs3zx/jAgCAdsZhWZbl7yHagsfjkdPpVENDgyIjI23dd0L+Jlv3B8DXpwsz/D0CAD+52N/fXAsJAAAYh4ABAADGIWAAAIBxCBgAAGAcAgYAABiHgAEAAMYhYAAAgHEIGAAAYBwCBgAAGIeAAQAAxiFgAACAcQgYAABgHAIGAAAYh4ABAADGIWAAAIBxCBgAAGAcAgYAABiHgAEAAMYhYAAAgHEIGAAAYBwCBgAAGIeAAQAAxiFgAACAcQgYAABgnCB/DwAAbSUhf5Mt+/l0YYYt+wFgH87AAAAA4xAwAADAOAQMAAAwDgEDAACMQ8AAAADjEDAAAMA4BAwAADAOAQMAAIxDwAAAAOMQMAAAwDgEDAAAMA4BAwAAjEPAAAAA4xAwAADAOEH+HgAAvishf5O/RwDQznEGBgAAGIeAAQAAxiFgAACAcQgYAABgHAIGAAAYh4ABAADGIWAAAIBxCBgAAGAcAgYAABiHgAEAAMYhYAAAgHG4FhIAXIBd12b6dGGGLfsBwBkYAABgIAIGAAAYh4ABAADGIWAAAIBxbH8Tb2FhoV599VXt379fYWFhGjFihBYtWqS+fft6tzlz5oxmzJih9evXq7GxUenp6XrxxRcVExPj3aampkbTpk3Tli1b1LVrV2VlZamwsFBBQbzvGICZeDMwYB/bz8Bs3bpV2dnZ2rFjh0pLS9XU1KQxY8bo1KlT3m2mT5+uN954Q6+88oq2bt2qY8eO6fbbb/eub25uVkZGhs6ePavt27drzZo1Wr16tebOnWv3uAAAwEAOy7KstnyAL774Qj169NDWrVt14403qqGhQf/wD/+gdevW6Ve/+pUkaf/+/erfv7/Ky8s1fPhwvfnmm/rlL3+pY8eOec/KFBcXa/bs2friiy8UHBx8wcf1eDxyOp1qaGhQZGSkrcdk19+iAOBScAYGHdnF/v5u8/fANDQ0SJKio6MlSVVVVWpqalJaWpp3m379+ik+Pl7l5eWSpPLycg0cONDnJaX09HR5PB7t3bu31cdpbGyUx+PxuQEAgI6pTQOmpaVFubm5GjlypK677jpJUm1trYKDgxUVFeWzbUxMjGpra73bfDtezq0/t641hYWFcjqd3ltcXJzNRwMAANqLNg2Y7Oxs7dmzR+vXr2/Lh5EkzZkzRw0NDd7b0aNH2/wxAQCAf7TZR3pycnK0ceNGbdu2TT179vQud7lcOnv2rE6cOOFzFqaurk4ul8u7zc6dO332V1dX513XmpCQEIWEhNh8FAAAoD2y/QyMZVnKycnRhg0btHnzZiUmJvqsHzp0qLp06aKysjLvsgMHDqimpkapqamSpNTUVO3evVv19fXebUpLSxUZGamkpCS7RwYAAIax/QxMdna21q1bpz/96U+KiIjwvmfF6XQqLCxMTqdTU6ZMUV5enqKjoxUZGakHHnhAqampGj58uCRpzJgxSkpK0qRJk7R48WLV1tbq0UcfVXZ2NmdZAACA/QGzfPlySdIvfvELn+WrVq3SPffcI0lasmSJAgICNGHCBJ8vsjsnMDBQGzdu1LRp05Samqrw8HBlZWVp3rx5do8LAAAM1ObfA+MvfA8MgI6K74FBR9ZuvgcGAADAbgQMAAAwDgEDAACMQ8AAAADjEDAAAMA4BAwAADAOAQMAAIxDwAAAAOMQMAAAwDgEDAAAMA4BAwAAjEPAAAAA4xAwAADAOAQMAAAwDgEDAACMQ8AAAADjEDAAAMA4BAwAADAOAQMAAIxDwAAAAOMQMAAAwDgEDAAAME6QvwcAAPw4CfmbbNnPpwszbNkP4A+cgQEAAMYhYAAAgHEIGAAAYBwCBgAAGIeAAQAAxiFgAACAcQgYAABgHAIGAAAYh4ABAADGIWAAAIBxCBgAAGAcAgYAABiHizkCQCfFRSFhMs7AAAAA4xAwAADAOAQMAAAwDgEDAACMQ8AAAADjEDAAAMA4BAwAADAOAQMAAIxDwAAAAOMQMAAAwDgEDAAAMA4BAwAAjEPAAAAA43A1agDAZeGq1vAHzsAAAADjEDAAAMA4BAwAADAO74EBAHQovCencyBgAADtgl3hgc6Bl5AAAIBxCBgAAGCcdh0wy5YtU0JCgkJDQ5WSkqKdO3f6eyQAANAOtNuAefnll5WXl6fHHntMu3bt0qBBg5Senq76+np/jwYAAPzMYVmW5e8hWpOSkqJhw4bphRdekCS1tLQoLi5ODzzwgPLz8y94f4/HI6fTqYaGBkVGRto6G280AwBcLLs+zdTefve01ae0Lvb3d7v8FNLZs2dVVVWlOXPmeJcFBAQoLS1N5eXlrd6nsbFRjY2N3p8bGhokffMvwm4tjadt3ycAoGOKn/6Kv0doE23x+/Xb+73Q+ZV2GTB/+9vf1NzcrJiYGJ/lMTEx2r9/f6v3KSws1BNPPHHe8ri4uDaZEQCAzsz5XNvu/6uvvpLT6fze9e0yYC7FnDlzlJeX5/25paVFx48fV7du3eRwOHy29Xg8iouL09GjR21/eam96AzHKHWO4+QYO47OcJwcY8fhr+O0LEtfffWV3G73D27XLgOme/fuCgwMVF1dnc/yuro6uVyuVu8TEhKikJAQn2VRUVE/+DiRkZEd+j8+qXMco9Q5jpNj7Dg6w3FyjB2HP47zh868nNMuP4UUHBysoUOHqqyszLuspaVFZWVlSk1N9eNkAACgPWiXZ2AkKS8vT1lZWUpOTtYNN9yg5557TqdOndJvfvMbf48GAAD8rN0GzF133aUvvvhCc+fOVW1trQYPHqySkpLz3th7KUJCQvTYY4+d95JTR9IZjlHqHMfJMXYcneE4OcaOo70fZ7v9HhgAAIDv0y7fAwMAAPBDCBgAAGAcAgYAABiHgAEAAMbpdAGzbNkyJSQkKDQ0VCkpKdq5c6e/R7JVYWGhhg0bpoiICPXo0UPjx4/XgQMH/D1Wm1q4cKEcDodyc3P9PYqt/vrXv+ruu+9Wt27dFBYWpoEDB+r999/391i2am5uVkFBgRITExUWFqZrr71WTz755AWvgdKebdu2TbfccovcbrccDodee+01n/WWZWnu3LmKjY1VWFiY0tLSdPDgQf8Mexl+6Dibmpo0e/ZsDRw4UOHh4XK73Zo8ebKOHTvmv4EvwYX+LL/t3//93+VwOPTcc89dsfnscDHHuG/fPt16661yOp0KDw/XsGHDVFNTc+WH/Y5OFTAvv/yy8vLy9Nhjj2nXrl0aNGiQ0tPTVV9f7+/RbLN161ZlZ2drx44dKi0tVVNTk8aMGaNTp075e7Q2UVlZqd/97ne6/vrr/T2Krb788kuNHDlSXbp00ZtvvqmPP/5YzzzzjK6++mp/j2arRYsWafny5XrhhRe0b98+LVq0SIsXL9bzzz/v79Eu2alTpzRo0CAtW7as1fWLFy9WUVGRiouLVVFRofDwcKWnp+vMmTNXeNLL80PHefr0ae3atUsFBQXatWuXXn31VR04cEC33nqrHya9dBf6szxnw4YN2rFjxwW/+r49utAxHjp0SKNGjVK/fv30zjvv6KOPPlJBQYFCQ0Ov8KStsDqRG264wcrOzvb+3NzcbLndbquwsNCPU7Wt+vp6S5K1detWf49iu6+++srq3bu3VVpaav385z+3HnroIX+PZJvZs2dbo0aN8vcYbS4jI8O69957fZbdfvvtVmZmpp8mspcka8OGDd6fW1paLJfLZf32t7/1Ljtx4oQVEhJi/dd//ZcfJrTHd4+zNTt37rQkWUeOHLkyQ9ns+47xs88+s37yk59Ye/bssXr16mUtWbLkis9ml9aO8a677rLuvvtu/wx0AZ3mDMzZs2dVVVWltLQ077KAgAClpaWpvLzcj5O1rYaGBklSdHS0nyexX3Z2tjIyMnz+TDuK119/XcnJybrjjjvUo0cPDRkyRCtXrvT3WLYbMWKEysrK9Mknn0iSPvzwQ7377ru6+eab/TxZ2zh8+LBqa2t9/pt1Op1KSUnp0M9D0jfPRQ6H44LXqDNJS0uLJk2apFmzZmnAgAH+Hsd2LS0t2rRpk/r06aP09HT16NFDKSkpP/hS2pXUaQLmb3/7m5qbm8/7Jt+YmBjV1tb6aaq21dLSotzcXI0cOVLXXXedv8ex1fr167Vr1y4VFhb6e5Q28Ze//EXLly9X79699dZbb2natGl68MEHtWbNGn+PZqv8/HxNnDhR/fr1U5cuXTRkyBDl5uYqMzPT36O1iXPPNZ3peUiSzpw5o9mzZ+vXv/51h7r44aJFixQUFKQHH3zQ36O0ifr6ep08eVILFy7U2LFj9fbbb+u2227T7bffrq1bt/p7vPZ7KQFcvuzsbO3Zs0fvvvuuv0ex1dGjR/XQQw+ptLS0fbwO2wZaWlqUnJysBQsWSJKGDBmiPXv2qLi4WFlZWX6ezj5//OMftXbtWq1bt04DBgxQdXW1cnNz5Xa7O9RxdmZNTU268847ZVmWli9f7u9xbFNVVaWlS5dq165dcjgc/h6nTbS0tEiS/vmf/1nTp0+XJA0ePFjbt29XcXGxfv7zn/tzvM5zBqZ79+4KDAxUXV2dz/K6ujq5XC4/TdV2cnJytHHjRm3ZskU9e/b09zi2qqqqUn19vX76058qKChIQUFB2rp1q4qKihQUFKTm5mZ/j3jZYmNjlZSU5LOsf//+7eKd/3aaNWuW9yzMwIEDNWnSJE2fPr3Dnlk791zTWZ6HzsXLkSNHVFpa2qHOvvz5z39WfX294uPjvc9DR44c0YwZM5SQkODv8WzRvXt3BQUFtdvnok4TMMHBwRo6dKjKysq8y1paWlRWVqbU1FQ/TmYvy7KUk5OjDRs2aPPmzUpMTPT3SLYbPXq0du/ererqau8tOTlZmZmZqq6uVmBgoL9HvGwjR4487+Pvn3zyiXr16uWnidrG6dOnFRDg+zQUGBjo/ZtfR5OYmCiXy+XzPOTxeFRRUdGhnoek/4+XgwcP6r//+7/VrVs3f49kq0mTJumjjz7yeR5yu92aNWuW3nrrLX+PZ4vg4GANGzas3T4XdaqXkPLy8pSVlaXk5GTdcMMNeu6553Tq1Cn95je/8fdotsnOzta6dev0pz/9SREREd7X1Z1Op8LCwvw8nT0iIiLOe09PeHi4unXr1mHe6zN9+nSNGDFCCxYs0J133qmdO3dqxYoVWrFihb9Hs9Utt9yip556SvHx8RowYIA++OADPfvss7r33nv9PdolO3nypP7nf/7H+/Phw4dVXV2t6OhoxcfHKzc3V/Pnz1fv3r2VmJiogoICud1ujR8/3n9DX4IfOs7Y2Fj96le/0q5du7Rx40Y1Nzd7n4uio6MVHBzsr7F/lAv9WX43yrp06SKXy6W+ffte6VEv2YWOcdasWbrrrrt044036qabblJJSYneeOMNvfPOO/4b+hx/fwzqSnv++eet+Ph4Kzg42LrhhhusHTt2+HskW0lq9bZq1Sp/j9amOtrHqC3Lst544w3ruuuus0JCQqx+/fpZK1as8PdItvN4PNZDDz1kxcfHW6GhodY111xjPfLII1ZjY6O/R7tkW7ZsafX/waysLMuyvvkodUFBgRUTE2OFhIRYo0ePtg4cOODfoS/BDx3n4cOHv/e5aMuWLf4e/aJd6M/yu0z8GPXFHONLL71k/eM//qMVGhpqDRo0yHrttdf8N/C3OCzL4K+8BAAAnVKneQ8MAADoOAgYAABgHAIGAAAYh4ABAADGIWAAAIBxCBgAAGAcAgYAABiHgAEAAMYhYADoF7/4hXJzc3/Uffbv36/hw4crNDRUgwcPvqj7PP744z7b3nPPPVfsK/RXr16tqKgo27b97uyWZWnq1KmKjo6Ww+FQdXX1Jc8K4MI61bWQALTu1VdfVZcuXX7UfR577DGFh4frwIED6tq1axtN1n4tXbpU3/4i85KSEq1evVrvvPOOrrnmGnXv3l0Oh0MbNmww7jpHgAkIGACKjo7+0fc5dOiQMjIyruhVaZubm+VwOM67grU/OJ1On58PHTqk2NhYjRgxwk8TAZ2L/58FAPjdd19CSkhI0IIFC3TvvfcqIiJC8fHxPlfCdjgcqqqq0rx58+RwOPT4449LkmbPnq0+ffroqquu0jXXXKOCggI1NTVd8lznXsp5/fXXlZSUpJCQENXU1KixsVEzZ87UT37yE4WHhyslJeW8q+OuXr1a8fHxuuqqq3Tbbbfp73//u8/6Dz/8UDfddJMiIiIUGRmpoUOH6v333/fZ5q233lL//v3VtWtXjR07Vp9//rl33bdfQrrnnnv0wAMPqKamRg6HQwkJCUpISJAk3Xbbbd5lAOxDwABo1TPPPKPk5GR98MEHuv/++zVt2jQdOHBAkvT5559rwIABmjFjhj7//HPNnDlTkhQREaHVq1fr448/1tKlS7Vy5UotWbLksuY4ffq0Fi1apP/8z//U3r171aNHD+Xk5Ki8vFzr16/XRx99pDvuuENjx47VwYMHJUkVFRWaMmWKcnJyVF1drZtuuknz58/32W9mZqZ69uypyspKVVVVKT8/3+dltNOnT+vpp5/W73//e23btk01NTXe4/yupUuXat68eerZs6c+//xzVVZWqrKyUpK0atUq7zIA9uElJACtGjdunO6//35J35xZWbJkibZs2aK+ffvK5XIpKChIXbt2lcvl8t7n0Ucf9f5zQkKCZs6cqfXr1+vhhx++5Dmampr04osvatCgQZKkmpoarVq1SjU1NXK73ZKkmTNnqqSkRKtWrdKCBQu0dOlSjR071vu4ffr00fbt21VSUuLdb01NjWbNmqV+/fpJknr37n3e4xYXF+vaa6+VJOXk5GjevHmtzuh0OhUREaHAwECffx+SFBUVdd4yAJePgAHQquuvv977zw6HQy6XS/X19T94n5dffllFRUU6dOiQTp48qa+//lqRkZGXNUdwcLDPLLt371Zzc7P69Onjs11jY6O6desmSdq3b59uu+02n/Wpqak+AZOXl6d//dd/1e9//3ulpaXpjjvu8MaKJF111VU+P8fGxl7w+AFcObyEBKBV3/1UksPhUEtLy/duX15erszMTI0bN04bN27UBx98oEceeURnz569rDnCwsLkcDi8P588eVKBgYGqqqpSdXW197Zv3z4tXbr0ovf7+OOPa+/evcrIyNDmzZuVlJSkDRs2eNe3dvzf/tQRAP/iDAwAW2zfvl29evXSI4884l125MgR2x9nyJAham5uVn19vX72s5+1uk3//v1VUVHhs2zHjh3nbdenTx/16dNH06dP169//WutWrXqvDM3l6NLly5qbm62bX8A/h9nYADYonfv3qqpqdH69et16NAhFRUV+ZzRsEufPn2UmZmpyZMn69VXX9Xhw4e1c+dOFRYWatOmTZKkBx98UCUlJXr66ad18OBBvfDCCz4vH/3v//6vcnJy9M477+jIkSN67733VFlZqf79+9s6a0JCgsrKylRbW6svv/zS1n0DnR0BA8AWt956q6ZPn66cnBwNHjxY27dvV0FBQZs81qpVqzR58mTNmDFDffv21fjx41VZWan4+HhJ0vDhw7Vy5UotXbpUgwYN0ttvv+3zBuPAwED9/e9/1+TJk9WnTx/deeeduvnmm/XEE0/YOuczzzyj0tJSxcXFaciQIbbuG+jsHBYv6gIAAMNwBgYAABiHgAEAAMYhYAAAgHEIGAAAYBwCBgAAGIeAAQAAxiFgAACAcQgYAABgHAIGAAAYh4ABAADGIWAAAIBx/g8cdyWoCaEjoQAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGwCAYAAAC3qV8qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxpklEQVR4nO3de1xVZb7H8e8W3SAG2ysghaCWeAVNk8g0S0dUjl2nUikxTbtApWRHmSnFOolpmVaOZaV2TjraNGml5QiaMiXesD3eilFT6QLYpLJFTyiwzh/zcp/ZgRW6mc2Dn/frtV4v1vM8a63fA9X+ttaz97ZZlmUJAADAIA18XQAAAEBNEWAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIzT0NcF1JbKykp99913CgoKks1m83U5AADgV7AsSydPnlR4eLgaNDj/fZZ6G2C+++47RURE+LoMAABwAb7++mtdccUV5+2vtwEmKChI0j9/AcHBwT6uBgAA/Boul0sRERHu1/HzqbcB5txjo+DgYAIMAACG+aXlHyziBQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABinoa8LAAAAvyxqyhqvnOfwzESvnMfXuAMDAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxqlxgMnJydGwYcMUHh4um82mVatWefTbbLZqt9mzZ7vHREVFVemfOXOmx3l27dqlvn37KiAgQBEREZo1a9aFzRAAANQ7NQ4wp06dUmxsrObPn19tf2Fhoce2aNEi2Ww23XHHHR7jnn76aY9xjzzyiLvP5XJp0KBBioyMVF5enmbPnq2MjAwtXLiwpuUCAIB6qGFNDxgyZIiGDBly3v6wsDCP/ffff1833nij2rVr59EeFBRUZew5S5cu1ZkzZ7Ro0SLZ7XZ16dJFTqdTc+bM0fjx42taMgAAqGdqdQ1McXGx1qxZo7Fjx1bpmzlzplq0aKEePXpo9uzZKi8vd/fl5uaqX79+stvt7raEhATl5+fr+PHj1V6rrKxMLpfLYwMAAPVTje/A1MRbb72loKAg3X777R7tjz76qK6++mo1b95cmzdvVnp6ugoLCzVnzhxJUlFRkdq2betxTGhoqLuvWbNmVa6VmZmp6dOn19JMAABAXVKrAWbRokVKSkpSQECAR3taWpr755iYGNntdj3wwAPKzMyUv7//BV0rPT3d47wul0sREREXVjgAAKjTai3A/PWvf1V+fr5WrFjxi2Pj4uJUXl6uw4cPKzo6WmFhYSouLvYYc27/fOtm/P39Lzj8AAAAs9TaGpg333xTPXv2VGxs7C+OdTqdatCggUJCQiRJ8fHxysnJ0dmzZ91jsrKyFB0dXe3jIwAAcGmpcYApLS2V0+mU0+mUJB06dEhOp1MFBQXuMS6XS3/60590//33Vzk+NzdXc+fO1d/+9jd99dVXWrp0qSZOnKh77rnHHU5Gjhwpu92usWPHau/evVqxYoXmzZvn8YgIAABcumr8CGnHjh268cYb3fvnQkVycrKWLFkiSVq+fLksy9KIESOqHO/v76/ly5crIyNDZWVlatu2rSZOnOgRThwOh9atW6eUlBT17NlTLVu21NSpU3kLNQAAkCTZLMuyfF1EbXC5XHI4HCopKVFwcLCvywEA4KJETVnjlfMcnpnolfPUll/7+s13IQEAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgnBoHmJycHA0bNkzh4eGy2WxatWqVR//o0aNls9k8tsGDB3uMOXbsmJKSkhQcHKymTZtq7NixKi0t9Riza9cu9e3bVwEBAYqIiNCsWbNqPjsAAFAv1TjAnDp1SrGxsZo/f/55xwwePFiFhYXu7Y9//KNHf1JSkvbu3ausrCytXr1aOTk5Gj9+vLvf5XJp0KBBioyMVF5enmbPnq2MjAwtXLiwpuUCAIB6qGFNDxgyZIiGDBnys2P8/f0VFhZWbd8XX3yhtWvXavv27erVq5ck6eWXX9bQoUP1/PPPKzw8XEuXLtWZM2e0aNEi2e12denSRU6nU3PmzPEIOgAA4NJUK2tgNm7cqJCQEEVHR+uhhx7SDz/84O7Lzc1V06ZN3eFFkgYOHKgGDRpo69at7jH9+vWT3W53j0lISFB+fr6OHz9e7TXLysrkcrk8NgAAUD95PcAMHjxY//3f/63169frueee06ZNmzRkyBBVVFRIkoqKihQSEuJxTMOGDdW8eXMVFRW5x4SGhnqMObd/bsxPZWZmyuFwuLeIiAhvTw0AANQRNX6E9EuGDx/u/rlbt26KiYlR+/bttXHjRg0YMMDbl3NLT09XWlqae9/lchFiAACop2r9bdTt2rVTy5YtdeDAAUlSWFiYjh496jGmvLxcx44dc6+bCQsLU3FxsceYc/vnW1vj7++v4OBgjw0AANRPtR5gvvnmG/3www9q3bq1JCk+Pl4nTpxQXl6ee8yGDRtUWVmpuLg495icnBydPXvWPSYrK0vR0dFq1qxZbZcMAADquBoHmNLSUjmdTjmdTknSoUOH5HQ6VVBQoNLSUj3xxBPasmWLDh8+rPXr1+uWW27RlVdeqYSEBElSp06dNHjwYI0bN07btm3TZ599ptTUVA0fPlzh4eGSpJEjR8put2vs2LHau3evVqxYoXnz5nk8IgIAAJeuGgeYHTt2qEePHurRo4ckKS0tTT169NDUqVPl5+enXbt26eabb1aHDh00duxY9ezZU3/961/l7+/vPsfSpUvVsWNHDRgwQEOHDtX111/v8RkvDodD69at06FDh9SzZ089/vjjmjp1Km+hBgAAkiSbZVmWr4uoDS6XSw6HQyUlJayHAQAYL2rKGq+c5/DMRK+cp7b82tdvvgsJAAAYhwADAACMQ4ABAADG8foH2QEAgP/nrbUr8MQdGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAODUOMDk5ORo2bJjCw8Nls9m0atUqd9/Zs2c1efJkdevWTU2aNFF4eLhGjRql7777zuMcUVFRstlsHtvMmTM9xuzatUt9+/ZVQECAIiIiNGvWrAubIQAAqHdqHGBOnTql2NhYzZ8/v0rf6dOntXPnTj311FPauXOn3nvvPeXn5+vmm2+uMvbpp59WYWGhe3vkkUfcfS6XS4MGDVJkZKTy8vI0e/ZsZWRkaOHChTUtFwAA1EMNa3rAkCFDNGTIkGr7HA6HsrKyPNpeeeUV9e7dWwUFBWrTpo27PSgoSGFhYdWeZ+nSpTpz5owWLVoku92uLl26yOl0as6cORo/fnxNSwYAAPVMra+BKSkpkc1mU9OmTT3aZ86cqRYtWqhHjx6aPXu2ysvL3X25ubnq16+f7Ha7uy0hIUH5+fk6fvx4tdcpKyuTy+Xy2AAAQP1U4zswNfHjjz9q8uTJGjFihIKDg93tjz76qK6++mo1b95cmzdvVnp6ugoLCzVnzhxJUlFRkdq2betxrtDQUHdfs2bNqlwrMzNT06dPr8XZAACAuqLWAszZs2d11113ybIsLViwwKMvLS3N/XNMTIzsdrseeOABZWZmyt/f/4Kul56e7nFel8uliIiICyseAADUabUSYM6FlyNHjmjDhg0ed1+qExcXp/Lych0+fFjR0dEKCwtTcXGxx5hz++dbN+Pv73/B4QcAAJjF62tgzoWX/fv3Kzs7Wy1atPjFY5xOpxo0aKCQkBBJUnx8vHJycnT27Fn3mKysLEVHR1f7+AgAAFxaanwHprS0VAcOHHDvHzp0SE6nU82bN1fr1q3129/+Vjt37tTq1atVUVGhoqIiSVLz5s1lt9uVm5urrVu36sYbb1RQUJByc3M1ceJE3XPPPe5wMnLkSE2fPl1jx47V5MmTtWfPHs2bN08vvviil6YNAABMZrMsy6rJARs3btSNN95YpT05OVkZGRlVFt+e88knn6h///7auXOnHn74YX355ZcqKytT27Ztde+99yotLc3jEdCuXbuUkpKi7du3q2XLlnrkkUc0efLkX12ny+WSw+FQSUnJLz7CAgCgtkRNWePrEjwcnpno6xJ+1q99/a5xgDEFAQYAUBcQYGrm175+811IAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADBOjb+NGgCAS0Fd+w4jeOIODAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYJwaB5icnBwNGzZM4eHhstlsWrVqlUe/ZVmaOnWqWrdurcaNG2vgwIHav3+/x5hjx44pKSlJwcHBatq0qcaOHavS0lKPMbt27VLfvn0VEBCgiIgIzZo1q+azAwAA9VKNA8ypU6cUGxur+fPnV9s/a9YsvfTSS3r11Ve1detWNWnSRAkJCfrxxx/dY5KSkrR3715lZWVp9erVysnJ0fjx4939LpdLgwYNUmRkpPLy8jR79mxlZGRo4cKFFzBFAABQ39gsy7Iu+GCbTStXrtStt94q6Z93X8LDw/X4449r0qRJkqSSkhKFhoZqyZIlGj58uL744gt17txZ27dvV69evSRJa9eu1dChQ/XNN98oPDxcCxYs0O9//3sVFRXJbrdLkqZMmaJVq1bpyy+//FW1uVwuORwOlZSUKDg4+EKnCAC4REVNWePrEmrF4ZmJvi7hZ/3a12+vroE5dOiQioqKNHDgQHebw+FQXFyccnNzJUm5ublq2rSpO7xI0sCBA9WgQQNt3brVPaZfv37u8CJJCQkJys/P1/Hjx6u9dllZmVwul8cGAADqJ68GmKKiIklSaGioR3toaKi7r6ioSCEhIR79DRs2VPPmzT3GVHeOf73GT2VmZsrhcLi3iIiIi58QAACokxr6ugBvSU9PV1pamnvf5XIRYgAA+AlvPRrz9aMor96BCQsLkyQVFxd7tBcXF7v7wsLCdPToUY/+8vJyHTt2zGNMdef412v8lL+/v4KDgz02AABQP3k1wLRt21ZhYWFav369u83lcmnr1q2Kj4+XJMXHx+vEiRPKy8tzj9mwYYMqKysVFxfnHpOTk6OzZ8+6x2RlZSk6OlrNmjXzZskAAMBANQ4wpaWlcjqdcjqdkv65cNfpdKqgoEA2m00TJkzQf/3Xf+mDDz7Q7t27NWrUKIWHh7vfqdSpUycNHjxY48aN07Zt2/TZZ58pNTVVw4cPV3h4uCRp5MiRstvtGjt2rPbu3asVK1Zo3rx5Ho+IAADApavGa2B27NihG2+80b1/LlQkJydryZIl+s///E+dOnVK48eP14kTJ3T99ddr7dq1CggIcB+zdOlSpaamasCAAWrQoIHuuOMOvfTSS+5+h8OhdevWKSUlRT179lTLli01depUj8+KAQAAl66L+hyYuozPgQEAXIz6+jkw3lJbi3h98jkwAAAA/w4EGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcerNt1EDACDxAXSXCu7AAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOF4PMFFRUbLZbFW2lJQUSVL//v2r9D344IMe5ygoKFBiYqICAwMVEhKiJ554QuXl5d4uFQAAGKqht0+4fft2VVRUuPf37Nmj3/zmN7rzzjvdbePGjdPTTz/t3g8MDHT/XFFRocTERIWFhWnz5s0qLCzUqFGj1KhRI82YMcPb5QIAAAN5PcC0atXKY3/mzJlq3769brjhBndbYGCgwsLCqj1+3bp12rdvn7KzsxUaGqru3bvrmWee0eTJk5WRkSG73e7tkgEAgGFqdQ3MmTNn9Pbbb2vMmDGy2Wzu9qVLl6ply5bq2rWr0tPTdfr0aXdfbm6uunXrptDQUHdbQkKCXC6X9u7de95rlZWVyeVyeWwAAKB+8vodmH+1atUqnThxQqNHj3a3jRw5UpGRkQoPD9euXbs0efJk5efn67333pMkFRUVeYQXSe79oqKi814rMzNT06dP9/4kAABAnVOrAebNN9/UkCFDFB4e7m4bP368++du3bqpdevWGjBggA4ePKj27dtf8LXS09OVlpbm3ne5XIqIiLjg8wEAgLqr1gLMkSNHlJ2d7b6zcj5xcXGSpAMHDqh9+/YKCwvTtm3bPMYUFxdL0nnXzUiSv7+//P39L7JqAABgglpbA7N48WKFhIQoMTHxZ8c5nU5JUuvWrSVJ8fHx2r17t44ePeoek5WVpeDgYHXu3Lm2ygUAAAaplTswlZWVWrx4sZKTk9Ww4f9f4uDBg1q2bJmGDh2qFi1aaNeuXZo4caL69eunmJgYSdKgQYPUuXNn3XvvvZo1a5aKior05JNPKiUlhTssAFCPRU1Z4+sSYJBaCTDZ2dkqKCjQmDFjPNrtdruys7M1d+5cnTp1ShEREbrjjjv05JNPusf4+flp9erVeuihhxQfH68mTZooOTnZ43NjAADApa1WAsygQYNkWVaV9oiICG3atOkXj4+MjNRHH31UG6UBAIB6gO9CAgAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjNPQ1wUAAMwWNWWNr0vAJYg7MAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwjtcDTEZGhmw2m8fWsWNHd/+PP/6olJQUtWjRQpdddpnuuOMOFRcXe5yjoKBAiYmJCgwMVEhIiJ544gmVl5d7u1QAAGCoWvkyxy5duig7O/v/L9Lw/y8zceJErVmzRn/605/kcDiUmpqq22+/XZ999pkkqaKiQomJiQoLC9PmzZtVWFioUaNGqVGjRpoxY0ZtlAsAAAxTKwGmYcOGCgsLq9JeUlKiN998U8uWLdNNN90kSVq8eLE6deqkLVu26Nprr9W6deu0b98+ZWdnKzQ0VN27d9czzzyjyZMnKyMjQ3a7vTZKBgAABqmVNTD79+9XeHi42rVrp6SkJBUUFEiS8vLydPbsWQ0cONA9tmPHjmrTpo1yc3MlSbm5uerWrZtCQ0PdYxISEuRyubR3797zXrOsrEwul8tjAwAA9ZPXA0xcXJyWLFmitWvXasGCBTp06JD69u2rkydPqqioSHa7XU2bNvU4JjQ0VEVFRZKkoqIij/Byrv9c3/lkZmbK4XC4t4iICO9ODAAA1Blef4Q0ZMgQ988xMTGKi4tTZGSk3nnnHTVu3Njbl3NLT09XWlqae9/lchFiAACop2r9bdRNmzZVhw4ddODAAYWFhenMmTM6ceKEx5ji4mL3mpmwsLAq70o6t1/duppz/P39FRwc7LEBAID6qdYDTGlpqQ4ePKjWrVurZ8+eatSokdavX+/uz8/PV0FBgeLj4yVJ8fHx2r17t44ePeoek5WVpeDgYHXu3Lm2ywUAAAbw+iOkSZMmadiwYYqMjNR3332nadOmyc/PTyNGjJDD4dDYsWOVlpam5s2bKzg4WI888oji4+N17bXXSpIGDRqkzp07695779WsWbNUVFSkJ598UikpKfL39/d2uQBwyYqassbXJQAXzOsB5ptvvtGIESP0ww8/qFWrVrr++uu1ZcsWtWrVSpL04osvqkGDBrrjjjtUVlamhIQE/eEPf3Af7+fnp9WrV+uhhx5SfHy8mjRpouTkZD399NPeLhUAABjKZlmW5esiaoPL5ZLD4VBJSQnrYQCgGtyBwcU4PDOxVs77a1+/+S4kAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAAAADjEGAAAIBxCDAAAMA4BBgAAGAcAgwAADAOAQYAABiHAAMAAIxDgAEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4zT0dQEAgJqJmrLG1yUAPscdGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxiHAAAAA4xBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADG8XqAyczM1DXXXKOgoCCFhITo1ltvVX5+vseY/v37y2azeWwPPvigx5iCggIlJiYqMDBQISEheuKJJ1ReXu7tcgEAgIEaevuEmzZtUkpKiq655hqVl5frd7/7nQYNGqR9+/apSZMm7nHjxo3T008/7d4PDAx0/1xRUaHExESFhYVp8+bNKiws1KhRo9SoUSPNmDHD2yUDAADDeD3ArF271mN/yZIlCgkJUV5envr16+duDwwMVFhYWLXnWLdunfbt26fs7GyFhoaqe/fueuaZZzR58mRlZGTIbrd7u2wAAGCQWl8DU1JSIklq3ry5R/vSpUvVsmVLde3aVenp6Tp9+rS7Lzc3V926dVNoaKi7LSEhQS6XS3v37q32OmVlZXK5XB4bAACon7x+B+ZfVVZWasKECerTp4+6du3qbh85cqQiIyMVHh6uXbt2afLkycrPz9d7770nSSoqKvIIL5Lc+0VFRdVeKzMzU9OnT6+lmQAAgLqkVgNMSkqK9uzZo08//dSjffz48e6fu3XrptatW2vAgAE6ePCg2rdvf0HXSk9PV1pamnvf5XIpIiLiwgoHAAB1Wq09QkpNTdXq1av1ySef6IorrvjZsXFxcZKkAwcOSJLCwsJUXFzsMebc/vnWzfj7+ys4ONhjAwAA9ZPXA4xlWUpNTdXKlSu1YcMGtW3b9hePcTqdkqTWrVtLkuLj47V7924dPXrUPSYrK0vBwcHq3Lmzt0sGAACG8fojpJSUFC1btkzvv/++goKC3GtWHA6HGjdurIMHD2rZsmUaOnSoWrRooV27dmnixInq16+fYmJiJEmDBg1S586dde+992rWrFkqKirSk08+qZSUFPn7+3u7ZAAAYBiv34FZsGCBSkpK1L9/f7Vu3dq9rVixQpJkt9uVnZ2tQYMGqWPHjnr88cd1xx136MMPP3Sfw8/PT6tXr5afn5/i4+N1zz33aNSoUR6fGwMAAC5dNsuyLF8XURtcLpccDodKSkpYDwOgToiassbXJQBec3hmYq2c99e+fvNdSAAAwDgEGAAAYBwCDAAAME6tfpAdANQHrF0B6h7uwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIe3UQOoc7z1tuXa+qhzAL5HgAFQb/H5LUD9xSMkAABgHAIMAAAwDgEGAAAYhwADAACMwyJeAF7DolkA/y7cgQEAAMYhwAAAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIfPgQHA57cAMA53YAAAgHG4AwPUgLfuVByemeiV83DnBMCligAD+ADBAwAuDo+QAACAcbgDAw917RGJt3DHAwDqFwIM6jSCBwCgOgSYeqAuvsjXxZoAAPUHAeYC8OIMAIBvsYgXAAAYhwADAACMQ4ABAADGqdMBZv78+YqKilJAQIDi4uK0bds2X5cEAADqgDobYFasWKG0tDRNmzZNO3fuVGxsrBISEnT06FFflwYAAHyszgaYOXPmaNy4cbrvvvvUuXNnvfrqqwoMDNSiRYt8XRoAAPCxOvk26jNnzigvL0/p6enutgYNGmjgwIHKzc2t9piysjKVlZW590tKSiRJLpfL6/VVlp32+jkBADBJbby+/ut5Lcv62XF1MsD84x//UEVFhUJDQz3aQ0ND9eWXX1Z7TGZmpqZPn16lPSIiolZqBADgUuaYW7vnP3nypBwOx3n762SAuRDp6elKS0tz71dWVurYsWNq0aKFbDabx1iXy6WIiAh9/fXXCg4O/neXWqvq89yk+j2/+jw3qX7Prz7PTWJ+JjNxbpZl6eTJkwoPD//ZcXUywLRs2VJ+fn4qLi72aC8uLlZYWFi1x/j7+8vf39+jrWnTpj97neDgYGP+oDVVn+cm1e/51ee5SfV7fvV5bhLzM5lpc/u5Oy/n1MlFvHa7XT179tT69evdbZWVlVq/fr3i4+N9WBkAAKgL6uQdGElKS0tTcnKyevXqpd69e2vu3Lk6deqU7rvvPl+XBgAAfKzOBpi7775b33//vaZOnaqioiJ1795da9eurbKw90L4+/tr2rRpVR451Qf1eW5S/Z5ffZ6bVL/nV5/nJjE/k9XnudmsX3qfEgAAQB1TJ9fAAAAA/BwCDAAAMA4BBgAAGIcAAwAAjHPJBZj58+crKipKAQEBiouL07Zt23xdklfk5ORo2LBhCg8Pl81m06pVq3xdktdkZmbqmmuuUVBQkEJCQnTrrbcqPz/f12V5zYIFCxQTE+P+oKn4+Hh9/PHHvi6rVsycOVM2m00TJkzwdSlekZGRIZvN5rF17NjR12V51bfffqt77rlHLVq0UOPGjdWtWzft2LHD12VdtKioqCp/O5vNppSUFF+X5hUVFRV66qmn1LZtWzVu3Fjt27fXM88884vfL2SSSyrArFixQmlpaZo2bZp27typ2NhYJSQk6OjRo74u7aKdOnVKsbGxmj9/vq9L8bpNmzYpJSVFW7ZsUVZWls6ePatBgwbp1KlTvi7NK6644grNnDlTeXl52rFjh2666Sbdcsst2rt3r69L86rt27frtddeU0xMjK9L8aouXbqosLDQvX366ae+Lslrjh8/rj59+qhRo0b6+OOPtW/fPr3wwgtq1qyZr0u7aNu3b/f4u2VlZUmS7rzzTh9X5h3PPfecFixYoFdeeUVffPGFnnvuOc2aNUsvv/yyr0vzHusS0rt3byslJcW9X1FRYYWHh1uZmZk+rMr7JFkrV670dRm15ujRo5Yka9OmTb4updY0a9bMeuONN3xdhtecPHnSuuqqq6ysrCzrhhtusB577DFfl+QV06ZNs2JjY31dRq2ZPHmydf311/u6jH+Lxx57zGrfvr1VWVnp61K8IjEx0RozZoxH2+23324lJSX5qCLvu2TuwJw5c0Z5eXkaOHCgu61BgwYaOHCgcnNzfVgZaqqkpESS1Lx5cx9X4n0VFRVavny5Tp06Va++NiMlJUWJiYke//7VF/v371d4eLjatWunpKQkFRQU+Lokr/nggw/Uq1cv3XnnnQoJCVGPHj30+uuv+7osrztz5ozefvttjRkzpsqX/5rquuuu0/r16/X3v/9dkvS3v/1Nn376qYYMGeLjyrynzn4Sr7f94x//UEVFRZVP8g0NDdWXX37po6pQU5WVlZowYYL69Omjrl27+rocr9m9e7fi4+P1448/6rLLLtPKlSvVuXNnX5flFcuXL9fOnTu1fft2X5fidXFxcVqyZImio6NVWFio6dOnq2/fvtqzZ4+CgoJ8Xd5F++qrr7RgwQKlpaXpd7/7nbZv365HH31UdrtdycnJvi7Pa1atWqUTJ05o9OjRvi7Fa6ZMmSKXy6WOHTvKz89PFRUVevbZZ5WUlOTr0rzmkgkwqB9SUlK0Z8+eerXOQJKio6PldDpVUlKid999V8nJydq0aZPxIebrr7/WY489pqysLAUEBPi6HK/71/+bjYmJUVxcnCIjI/XOO+9o7NixPqzMOyorK9WrVy/NmDFDktSjRw/t2bNHr776ar0KMG+++aaGDBmi8PBwX5fiNe+8846WLl2qZcuWqUuXLnI6nZowYYLCw8Przd/ukgkwLVu2lJ+fn4qLiz3ai4uLFRYW5qOqUBOpqalavXq1cnJydMUVV/i6HK+y2+268sorJUk9e/bU9u3bNW/ePL322ms+ruzi5OXl6ejRo7r66qvdbRUVFcrJydErr7yisrIy+fn5+bBC72ratKk6dOigAwcO+LoUr2jdunWVEN2pUyf9+c9/9lFF3nfkyBFlZ2frvffe83UpXvXEE09oypQpGj58uCSpW7duOnLkiDIzM+tNgLlk1sDY7Xb17NlT69evd7dVVlZq/fr19WqtQX1kWZZSU1O1cuVKbdiwQW3btvV1SbWusrJSZWVlvi7jog0YMEC7d++W0+l0b7169VJSUpKcTme9Ci+SVFpaqoMHD6p169a+LsUr+vTpU+UjC/7+978rMjLSRxV53+LFixUSEqLExERfl+JVp0+fVoMGni/xfn5+qqys9FFF3nfJ3IGRpLS0NCUnJ6tXr17q3bu35s6dq1OnTum+++7zdWkXrbS01OP/+g4dOiSn06nmzZurTZs2Pqzs4qWkpGjZsmV6//33FRQUpKKiIkmSw+FQ48aNfVzdxUtPT9eQIUPUpk0bnTx5UsuWLdPGjRv1l7/8xdelXbSgoKAqa5WaNGmiFi1a1Is1TJMmTdKwYcMUGRmp7777TtOmTZOfn59GjBjh69K8YuLEibruuus0Y8YM3XXXXdq2bZsWLlyohQsX+ro0r6isrNTixYuVnJyshg3r18vhsGHD9Oyzz6pNmzbq0qWLPv/8c82ZM0djxozxdWne4+u3Qf27vfzyy1abNm0su91u9e7d29qyZYuvS/KKTz75xJJUZUtOTvZ1aRetunlJshYvXuzr0rxizJgxVmRkpGW3261WrVpZAwYMsNatW+frsmpNfXob9d133221bt3astvt1uWXX27dfffd1oEDB3xdlld9+OGHVteuXS1/f3+rY8eO1sKFC31dktf85S9/sSRZ+fn5vi7F61wul/XYY49Zbdq0sQICAqx27dpZv//9762ysjJfl+Y1NsuqRx/LBwAALgmXzBoYAABQfxBgAACAcQgwAADAOAQYAABgHAIMAAAwDgEGAAAYhwADAACMQ4ABAADGIcAABti4caNsNptOnDjh61I8REVFae7cuT6toX///powYUKNjvnyyy917bXXKiAgQN27d/9Vx2RkZHiMHT16tG699dbzjl+yZIlsNptsNluN67tY567btGnTf+t1gX+n+vXlD0A90L9/f3Xv3t0jGFx33XUqLCyUw+HwXWHV2L59u5o0aeLTGt577z01atSoRsdMmzZNTZo0UX5+vi677LJaqkwKDg5Wfn5+ld/RgQMHNGPGDGVnZ6u4uFgtW7ZUx44dNWbMGN19990X/b08hYWFWrFihaZNm3ZR5wHqMgIMYAC73a6wsDBfl1FFq1atfF2CmjdvXuNjDh48qMTExFr/VmWbzVbl77Zt2zYNHDhQXbp00fz589WxY0dJ0o4dOzR//nx17dpVsbGxF3S9M2fOuP9ZqWthF/A2HiEBdcjo0aO1adMmzZs3z/0Y4PDhw1UeIS1ZskRNmzbV6tWrFR0drcDAQP32t7/V6dOn9dZbbykqKkrNmjXTo48+qoqKCvf5y8rKNGnSJF1++eVq0qSJ4uLitHHjxvPWY1mWMjIy1KZNG/n7+ys8PFyPPvqou/+nj5BsNpveeOMN3XbbbQoMDNRVV12lDz74wOOce/fu1X/8x38oODhYQUFB6tu3rw4ePOjuf+ONN9SpUycFBASoY8eO+sMf/vCzv7OfPkKKiorSjBkzNGbMGAUFBalNmzYe355ss9mUl5enp59+WjabTRkZGZKkyZMnq0OHDgoMDFS7du301FNP6ezZsz977ZqyLEujR49Whw4d9Nlnn2nYsGG66qqrdNVVV2nEiBH69NNPFRMTI0m66aablJqa6nH8999/L7vdrvXr17vn+swzz2jUqFEKDg7W+PHjvVovUJcRYIA6ZN68eYqPj9e4ceNUWFiowsJCRUREVDv29OnTeumll7R8+XKtXbtWGzdu1G233aaPPvpIH330kf7nf/5Hr732mt599133MampqcrNzdXy5cu1a9cu3XnnnRo8eLD2799f7TX+/Oc/68UXX9Rrr72m/fv3a9WqVerWrdvPzmH69Om66667tGvXLg0dOlRJSUk6duyYJOnbb79Vv3795O/vrw0bNigvL09jxoxReXm5JGnp0qWaOnWqnn32WX3xxReaMWOGnnrqKb311ls1+j2+8MIL6tWrlz7//HM9/PDDeuihh5Sfny/pn49XunTposcff1yFhYWaNGmSJCkoKEhLlizRvn37NG/ePL3++ut68cUXa3TdX+J0OvXFF19o0qRJatCg+v/82mw2SdL999+vZcuWqayszN339ttv6/LLL9dNN93kbnv++ecVGxurzz//XE899ZRX6wXqNN9+GTaAn7rhhhusxx57zKPtk08+sSRZx48ftyzLshYvXmxJsg4cOOAe88ADD1iBgYHWyZMn3W0JCQnWAw88YFmWZR05csTy8/Ozvv32W49zDxgwwEpPT6+2lhdeeMHq0KGDdebMmWr7IyMjrRdffNG9L8l68skn3fulpaWWJOvjjz+2LMuy0tPTrbZt2573fO3bt7eWLVvm0fbMM89Y8fHx1Y63rKq/r8jISOuee+5x71dWVlohISHWggUL3G2xsbHWtGnTzntOy7Ks2bNnWz179nTvT5s2zYqNjXXvJycnW7fccst5j1+8eLHlcDg82pYvX25Jsnbu3OluKy4utpo0aeLe5s+fb1mWZf3v//6v1axZM2vFihXusTExMVZGRobHXG+99dZffX2gPmENDGCowMBAtW/f3r0fGhqqqKgoj0WpoaGhOnr0qCRp9+7dqqioUIcOHTzOU1ZWphYtWlR7jTvvvFNz585Vu3btNHjwYA0dOlTDhg372UWm5x6BSFKTJk0UHBzsrsHpdKpv377VLro9deqUDh48qLFjx2rcuHHu9vLy8hqv5/jXGs6tQzlXw/msWLFCL730kg4ePKjS0lKVl5crODi4Rte9EC1atJDT6ZT0z8dhZ86ckSQFBATo3nvv1aJFi3TXXXdp586d2rNnT5VHcr169ar1GoG6iAADGOqnIcBms1XbVllZKUkqLS2Vn5+f8vLy5Ofn5zHufO/EiYiIUH5+vrKzs5WVlaWHH35Ys2fP1qZNm877zp+fq6Fx48bnnU9paakk6fXXX1dcXJxH30/r/SU/V0N1cnNzlZSUpOnTpyshIUEOh0PLly/XCy+8UKPr/pKrrrpKkpSfn68ePXpI+ufcrrzySkmqEgzvv/9+de/eXd98840WL16sm266qcrCY1+/CwzwFQIMUMfY7XaPhbfe0qNHD1VUVOjo0aPq27fvrz6ucePGGjZsmIYNG6aUlBR17NhRu3fv1tVXX13jGmJiYvTWW2/p7NmzVUJGaGiowsPD9dVXXykpKanG574YmzdvVmRkpH7/+9+7244cOeL16/To0UMdO3bU888/r7vuuuu862DO6datm3r16qXXX39dy5Yt0yuvvOL1mgBTEWCAOiYqKkpbt27V4cOHddlll13Q24Sr06FDByUlJWnUqFF64YUX1KNHD33//fdav369YmJilJiYWOWYJUuWqKKiQnFxcQoMDNTbb7+txo0bX/Dbj1NTU/Xyyy9r+PDhSk9Pl8Ph0JYtW9S7d29FR0dr+vTpevTRR+VwODR48GCVlZVpx44dOn78uNLS0i72V3BeV111lQoKCrR8+XJdc801WrNmjVauXOn169hsNi1evFi/+c1v1KdPH6Wnp6tTp046e/ascnJy9P3331e523T//fcrNTVVTZo00W233eb1mgBT8S4koI6ZNGmS/Pz81LlzZ7Vq1UoFBQVeO/fixYs1atQoPf7444qOjtatt96q7du3q02bNtWOb9q0qV5//XX16dNHMTExys7O1ocffnjeNTO/pEWLFtqwYYNKS0t1ww03qGfPnnr99dfdd2Puv/9+vfHGG1q8eLG6deumG264QUuWLFHbtm0veM6/xs0336yJEycqNTVV3bt31+bNm2vtHT3XXnut8vLyFB0drZSUFHXu3FnXXXed/vjHP+rFF1/UQw895DF+xIgRatiwoUaMGKGAgIBaqQkwkc2yLMvXRQBAfbNkyRJNmDDhor/+4fDhw2rfvr22b99eo8d23ro+UFdxBwYAaklJSYkuu+wyTZ48ucbHnj17VkVFRXryySd17bXX1ii8XHbZZXrwwQdrfE3AJNyBAYBacPLkSRUXF0v656O4li1b1uj4jRs36sYbb1SHDh307rvv/uIHCP6rAwcOSPrnO5xq+/Eb4CsEGAAAYBweIQEAAOMQYAAAgHEIMAAAwDgEGAAAYBwCDAAAMA4BBgAAGIcAAwAAjEOAAQAAxvk/wK91DofytIUAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "cosmo = Cosmology()\n", + "inf = InfallDistributionGalacticus2024(0.5)\n", + "samples = [inf.z_accreted_from_zlens(0.5) for _ in range(0, 10000)]\n", + "times = [cosmo.astropy.age(0.5).value - cosmo.astropy.age(samples_i).value for samples_i in samples]\n", + "plt.hist(samples,bins=25)\n", + "plt.xlabel('infall redshift')\n", + "plt.show()\n", + "plt.hist(times,bins=25)\n", + "plt.xlabel('time since infall [Gyr]')\n", + "plt.show()" + ] + }, { "cell_type": "code", "execution_count": null, diff --git a/pyHalo/Halos/HaloModels/TNFW.py b/pyHalo/Halos/HaloModels/TNFW.py index bafd719..0057216 100644 --- a/pyHalo/Halos/HaloModels/TNFW.py +++ b/pyHalo/Halos/HaloModels/TNFW.py @@ -114,17 +114,10 @@ def bound_mass(self): Computes the mass inside the virial radius (with truncation effects included) :return: the mass inside r = c * r_s """ - def _integrand(r, rs, rhos, rt): - x = r/rs - tau = rt/rs - return 4 * tau **2 * np.pi * r ** 2 * rhos / x / (1+x)**2 / (x**2 + tau**2) - - params_physical = self.params_physical - rhos = params_physical['rhos'] - rs = params_physical['rs'] - rt = params_physical['r_trunc_kpc'] - r200 = rs * self.c - #tau = params_physical['r_trunc_kpc'] / params_physical['rs'] - #f = tnfw_mass_fraction(tau, self.c) * self._rescale_norm - bound_mass = quad(_integrand, 0, r200, args=(rs, rhos, rt))[0] - return bound_mass + if hasattr(self, '_kwargs_lenstronomy'): + tau = self._kwargs_lenstronomy[0]['r_trunc'] / self._kwargs_lenstronomy[0]['Rs'] + else: + params_physical = self.params_physical + tau = params_physical['r_trunc_kpc'] / params_physical['rs'] + f = tnfw_mass_fraction(tau, self.c) + return f * self.mass diff --git a/pyHalo/Halos/galacticus_truncation/interp_mass_loss.py b/pyHalo/Halos/galacticus_truncation/interp_mass_loss.py index d517067..c51b2d4 100644 --- a/pyHalo/Halos/galacticus_truncation/interp_mass_loss.py +++ b/pyHalo/Halos/galacticus_truncation/interp_mass_loss.py @@ -80,28 +80,33 @@ class InterpGalacticus(object): def __init__(self): from pyHalo.Halos.galacticus_truncation.johnsonSUparams import a_fit, \ b_fit - log10c_values = np.linspace(np.log10(2.0), np.log10(384), 25) - t_inf_values = np.linspace(0.0, 8.1, 25) - a_values = np.array(a_fit).reshape(25, 25) - b_values = np.array(b_fit).reshape(25, 25) - _points = (t_inf_values, log10c_values) + log10c_values = np.linspace(np.log10(2.0), np.log10(384), 10) + t_inf_values = np.linspace(0.0, 8.1, 10) + chost_values = np.linspace(3.0, 9.0, 10) + a_values = np.array(a_fit).reshape(10, 10, 10) + b_values = np.array(b_fit).reshape(10, 10, 10) + _points = (t_inf_values, log10c_values, chost_values) self._a_interp = RegularGridInterpolator(_points, a_values, bounds_error=False, fill_value=None) self._b_interp = RegularGridInterpolator(_points, b_values, bounds_error=False, fill_value=None) - def __call__(self, log10_concentration_infall, time_since_infall): + def __call__(self, log10_concentration_infall, time_since_infall, chost): """ Evaluates the prediction from galacticus for subhalo bound mass :param log10_concentration_infall: log10(c) where c is the halo concentration at infall :param time_since_infall: the time ellapsed since infall and the deflector redshift + :param chost: host halo concentration at z=0.5 :return: the log10(bound mass divided by the infall mass), plus scatter """ log10_concentration_infall = max(np.log10(2), log10_concentration_infall) log10_concentration_infall = min(np.log10(384), log10_concentration_infall) time_since_infall = max(0.0, time_since_infall) time_since_infall = min(time_since_infall, 8.1) - p = (time_since_infall, log10_concentration_infall) + #chost = max(3.0, chost) + chost = max(6.0, chost) # we're getting some unphysical behavior at chost = 3.0 + chost = min(9.0, chost) + p = (time_since_infall, log10_concentration_infall, chost) a, b = self._a_interp(p), self._b_interp(p) output = float(johnsonsu.rvs(a, b)) return output diff --git a/pyHalo/Halos/galacticus_truncation/johnsonSUparams.py b/pyHalo/Halos/galacticus_truncation/johnsonSUparams.py index 773c4fe..74c5e8f 100644 --- a/pyHalo/Halos/galacticus_truncation/johnsonSUparams.py +++ b/pyHalo/Halos/galacticus_truncation/johnsonSUparams.py @@ -1,4 +1,5 @@ import numpy as np -a_fit = np.array([1.559, 1.631, 1.703, 1.774, 1.713, 1.625, 1.665, 1.704, 1.743, 1.733, 1.712, 1.691, 1.693, 1.708, 1.723, 1.688, 1.602, 1.515, 1.455, 1.449, 1.443, 1.447, 1.504, 1.56, 1.617, 3.857, 3.756, 3.654, 3.553, 3.553, 3.551, 3.55, 3.555, 3.554, 3.55, 3.549, 3.549, 3.548, 3.548, 3.548, 3.542, 3.557, 3.573, 3.588, 3.603, 3.617, 3.629, 3.643, 3.658, 3.679, 3.416, 3.314, 3.213, 3.146, 3.153, 3.151, 3.127, 3.126, 3.128, 3.129, 3.128, 3.154, 3.183, 3.186, 3.186, 3.193, 3.208, 3.219, 3.212, 3.221, 3.236, 3.251, 3.265, 3.28, 3.289, 2.829, 2.728, 2.654, 2.605, 2.623, 2.627, 2.579, 2.578, 2.585, 2.586, 2.605, 2.659, 2.688, 2.704, 2.705, 2.712, 2.727, 2.715, 2.708, 2.712, 2.722, 2.737, 2.751, 2.761, 2.77, 2.242, 2.163, 2.113, 2.063, 2.082, 2.103, 2.037, 2.033, 2.042, 2.056, 2.11, 2.164, 2.193, 2.209, 2.225, 2.232, 2.223, 2.211, 2.204, 2.208, 2.211, 2.223, 2.233, 2.242, 2.25, 1.848, 1.801, 1.751, 1.701, 1.72, 1.745, 1.691, 1.675, 1.684, 1.734, 1.788, 1.841, 1.865, 1.881, 1.898, 1.92, 1.908, 1.897, 1.89, 1.893, 1.896, 1.891, 1.898, 1.906, 1.914, 1.824, 1.782, 1.734, 1.685, 1.703, 1.705, 1.688, 1.672, 1.694, 1.744, 1.797, 1.836, 1.859, 1.875, 1.883, 1.916, 1.958, 1.947, 1.94, 1.943, 1.936, 1.916, 1.916, 1.924, 1.932, 1.8, 1.758, 1.717, 1.668, 1.677, 1.666, 1.685, 1.67, 1.704, 1.754, 1.793, 1.83, 1.853, 1.861, 1.865, 1.898, 1.961, 1.997, 1.99, 1.985, 1.961, 1.941, 1.941, 1.941, 1.95, 1.776, 1.735, 1.693, 1.651, 1.637, 1.626, 1.655, 1.68, 1.713, 1.751, 1.787, 1.823, 1.84, 1.844, 1.848, 1.88, 1.943, 2.006, 2.035, 2.011, 1.987, 1.967, 1.967, 1.967, 1.968, 1.961, 1.92, 1.878, 1.78, 1.741, 1.73, 1.758, 1.785, 1.818, 1.904, 1.941, 1.977, 1.992, 1.996, 1.997, 2.026, 2.09, 2.153, 2.187, 2.163, 2.139, 2.132, 2.136, 2.135, 2.135, 2.176, 2.134, 2.079, 1.905, 1.865, 1.854, 1.891, 1.906, 1.936, 2.042, 2.116, 2.152, 2.167, 2.171, 2.157, 2.186, 2.26, 2.322, 2.357, 2.333, 2.317, 2.314, 2.323, 2.324, 2.324, 2.39, 2.349, 2.204, 2.031, 1.99, 1.977, 2.027, 2.042, 2.054, 2.16, 2.282, 2.328, 2.343, 2.337, 2.317, 2.346, 2.423, 2.492, 2.526, 2.505, 2.5, 2.497, 2.506, 2.513, 2.513, 2.605, 2.503, 2.33, 2.156, 2.115, 2.101, 2.163, 2.178, 2.172, 2.278, 2.4, 2.504, 2.516, 2.497, 2.477, 2.505, 2.583, 2.66, 2.696, 2.687, 2.682, 2.679, 2.689, 2.698, 2.702, 2.904, 2.744, 2.571, 2.397, 2.356, 2.314, 2.374, 2.388, 2.375, 2.473, 2.594, 2.716, 2.746, 2.727, 2.707, 2.735, 2.813, 2.89, 2.977, 2.972, 2.967, 2.964, 2.974, 2.983, 2.965, 3.568, 3.477, 3.304, 3.13, 3.089, 2.85, 2.905, 2.909, 2.897, 2.994, 3.116, 3.232, 3.274, 3.254, 3.235, 3.263, 3.34, 3.417, 3.618, 3.692, 3.687, 3.685, 3.694, 3.703, 3.542, 4.232, 4.197, 4.037, 3.864, 3.823, 3.387, 3.436, 3.43, 3.418, 3.516, 3.637, 3.742, 3.786, 3.782, 3.762, 3.79, 3.868, 3.984, 4.192, 4.413, 4.408, 4.405, 4.414, 4.395, 4.119, 4.896, 4.86, 4.77, 4.597, 4.556, 3.923, 3.967, 3.961, 3.94, 4.037, 4.153, 4.252, 4.296, 4.309, 4.29, 4.318, 4.395, 4.557, 4.765, 5.052, 5.128, 5.126, 5.135, 4.972, 4.696, 5.56, 5.524, 5.489, 5.33, 5.143, 4.46, 4.499, 4.492, 4.461, 4.559, 4.663, 4.762, 4.806, 4.821, 4.817, 4.845, 4.962, 5.131, 5.338, 5.626, 5.849, 5.846, 5.824, 5.549, 5.273, 6.223, 6.188, 6.152, 6.063, 5.68, 4.997, 5.032, 5.023, 4.982, 5.074, 5.173, 5.272, 5.316, 5.331, 5.344, 5.373, 5.536, 5.704, 5.911, 6.199, 6.486, 6.567, 6.401, 6.126, 5.85, 6.887, 6.852, 6.816, 6.781, 6.216, 5.533, 5.569, 5.554, 5.509, 5.584, 5.683, 5.783, 5.826, 5.841, 5.856, 5.94, 6.109, 6.277, 6.485, 6.772, 7.06, 7.254, 6.979, 6.703, 6.427, 11.03, 10.995, 10.959, 9.658, 6.381, 5.698, 5.783, 5.963, 5.88, 5.775, 5.84, 5.939, 6.01, 6.039, 6.054, 6.144, 6.313, 6.477, 6.675, 6.962, 7.225, 7.354, 7.418, 7.142, 6.867, 15.554, 15.519, 14.065, 12.497, 8.139, 5.822, 6.179, 6.355, 6.215, 6.109, 6.012, 6.058, 6.129, 6.185, 6.217, 6.308, 6.47, 6.625, 6.823, 7.07, 7.197, 7.326, 7.462, 7.567, 7.291, 20.123, 18.473, 16.905, 15.338, 10.979, 6.084, 6.634, 6.751, 6.611, 6.451, 6.354, 6.256, 6.259, 6.315, 6.371, 6.464, 6.619, 6.774, 6.915, 7.043, 7.17, 7.299, 7.434, 7.57, 7.706, 26.634, 22.676, 19.761, 18.193, 13.834, 11.064, 9.678, 9.794, 8.306, 7.471, 7.01, 6.912, 6.915, 6.959, 6.795, 6.734, 6.777, 6.924, 7.034, 7.055, 7.157, 7.285, 7.421, 7.555, 7.687, 33.145, 29.187, 25.229, 21.271, 18.576, 16.043, 13.632, 11.683, 9.711, 8.875, 8.288, 7.701, 7.443, 7.279, 7.116, 7.054, 7.097, 7.141, 7.184, 7.205, 7.225, 7.272, 7.404, 7.536, 7.667, ]) -b_fit = np.array([1422.879, 1473.311, 1523.743, 1574.175, 1525.652, 1457.694, 1462.549, 1486.371, 1518.919, 1518.801, 1511.551, 1504.3, 1518.096, 1542.897, 1567.698, 1542.838, 1466.798, 1390.758, 1342.678, 1351.803, 1360.927, 1367.786, 1363.109, 1358.432, 1353.755, 94.951, 97.748, 100.546, 103.343, 108.351, 112.839, 118.008, 126.049, 131.879, 136.326, 143.059, 149.792, 157.115, 164.748, 172.38, 170.321, 179.153, 187.986, 197.028, 206.499, 215.969, 225.467, 235.424, 245.38, 253.115, 21.594, 24.391, 27.189, 28.444, 30.469, 34.956, 38.328, 43.945, 48.247, 53.929, 60.662, 65.064, 66.893, 73.49, 81.122, 89.345, 98.178, 105.774, 106.98, 112.334, 121.805, 131.355, 141.312, 151.269, 154.398, 14.849, 17.647, 19.186, 19.677, 20.409, 24.287, 25.766, 31.383, 33.139, 38.82, 43.808, 45.779, 47.609, 49.363, 56.683, 64.906, 73.446, 74.601, 75.807, 77.119, 83.476, 93.027, 102.983, 107.121, 108.88, 8.105, 9.928, 10.419, 10.91, 11.642, 13.618, 13.666, 17.163, 18.03, 22.552, 24.524, 26.495, 28.324, 30.079, 32.244, 40.467, 42.273, 43.428, 44.634, 45.946, 47.257, 54.698, 59.843, 61.603, 63.362, 3.636, 4.112, 4.603, 5.095, 5.827, 6.581, 6.695, 7.271, 8.138, 9.917, 11.889, 13.86, 15.482, 17.237, 18.992, 20.628, 21.783, 22.938, 24.145, 25.456, 26.767, 28.106, 29.807, 31.567, 33.326, 3.547, 3.99, 4.471, 4.962, 5.694, 6.37, 6.844, 7.42, 8.306, 10.086, 12.057, 13.444, 15.046, 16.801, 18.24, 19.971, 21.862, 23.018, 24.224, 25.535, 26.871, 28.242, 29.604, 31.339, 33.099, 3.459, 3.902, 4.345, 4.83, 5.527, 6.158, 6.994, 7.573, 8.474, 10.254, 11.699, 13.007, 14.609, 16.086, 17.401, 19.132, 21.292, 23.097, 24.303, 25.633, 27.006, 28.377, 29.739, 31.112, 32.871, 3.371, 3.814, 4.257, 4.7, 5.315, 5.946, 6.876, 7.741, 8.643, 9.953, 11.262, 12.571, 13.933, 15.248, 16.562, 18.294, 20.454, 22.614, 24.396, 25.769, 27.142, 28.513, 29.875, 31.237, 32.644, 3.24, 3.683, 4.126, 4.471, 4.946, 5.577, 6.506, 7.274, 8.175, 9.435, 10.744, 12.053, 13.366, 14.68, 15.954, 17.352, 19.494, 21.654, 23.556, 24.929, 26.302, 27.586, 29.024, 30.386, 31.748, 3.103, 3.546, 3.966, 4.181, 4.554, 5.185, 6.097, 6.742, 7.619, 8.888, 10.214, 11.523, 12.836, 14.151, 15.102, 16.499, 18.517, 20.677, 22.579, 23.952, 25.241, 26.478, 28.016, 29.396, 30.758, 2.967, 3.41, 3.675, 3.89, 4.239, 4.793, 5.68, 6.325, 7.062, 8.331, 9.68, 10.994, 12.306, 13.413, 14.249, 15.647, 17.623, 19.7, 21.602, 22.954, 24.133, 25.37, 26.908, 28.406, 29.768, 2.83, 3.169, 3.385, 3.6, 3.948, 4.401, 5.264, 5.909, 6.506, 7.775, 9.123, 10.464, 11.724, 12.561, 13.397, 14.794, 16.77, 18.746, 20.625, 21.846, 23.024, 24.262, 25.8, 27.337, 28.779, 2.761, 2.992, 3.207, 3.423, 3.771, 4.132, 4.972, 5.626, 6.173, 7.379, 8.727, 10.076, 11.151, 11.987, 12.823, 14.221, 16.197, 18.173, 20.12, 21.298, 22.477, 23.715, 25.252, 26.79, 28.251, 2.984, 3.298, 3.513, 3.728, 4.077, 4.341, 5.21, 5.913, 6.46, 7.666, 9.014, 10.373, 11.765, 12.602, 13.438, 14.835, 16.811, 18.787, 21.409, 23.138, 24.317, 25.554, 27.092, 28.63, 29.69, 3.207, 3.587, 3.819, 4.034, 4.382, 4.55, 5.448, 6.2, 6.747, 7.953, 9.301, 10.678, 12.106, 13.216, 14.052, 15.45, 17.426, 19.574, 22.225, 24.978, 26.157, 27.394, 28.932, 30.388, 31.129, 3.429, 3.81, 4.125, 4.34, 4.688, 4.759, 5.687, 6.443, 7.034, 8.24, 9.598, 10.984, 12.412, 13.83, 14.667, 16.064, 18.04, 20.39, 23.041, 26.255, 27.996, 29.234, 30.771, 31.827, 32.568, 3.652, 4.032, 4.413, 4.646, 4.922, 4.968, 5.925, 6.681, 7.321, 8.526, 9.903, 11.289, 12.717, 14.167, 15.281, 16.679, 18.829, 21.206, 23.858, 27.071, 29.836, 31.074, 32.525, 33.266, 34.007, 3.874, 4.255, 4.635, 4.951, 5.131, 5.177, 6.149, 6.92, 7.608, 8.823, 10.209, 11.594, 13.023, 14.473, 15.896, 17.293, 19.646, 22.022, 24.674, 27.887, 31.1, 32.913, 33.964, 34.705, 35.446, 4.097, 4.477, 4.858, 5.238, 5.34, 5.386, 6.358, 7.158, 7.87, 9.129, 10.514, 11.9, 13.328, 14.778, 16.229, 18.085, 20.462, 22.838, 25.49, 28.703, 31.916, 34.662, 35.403, 36.144, 36.885, 5.522, 5.903, 6.283, 5.928, 4.709, 4.755, 5.714, 6.465, 7.047, 7.63, 8.886, 10.272, 11.426, 12.737, 14.187, 16.094, 18.471, 20.566, 22.626, 25.839, 28.82, 30.594, 32.41, 33.151, 33.892, 7.079, 7.459, 7.016, 6.507, 4.769, 4.032, 4.92, 5.657, 6.043, 6.626, 7.251, 8.432, 9.587, 10.62, 11.889, 13.796, 15.804, 17.3, 19.359, 22.193, 23.913, 25.686, 27.728, 29.674, 30.415, 8.647, 8.112, 7.603, 7.094, 5.356, 3.383, 4.166, 4.698, 5.084, 5.632, 6.258, 6.884, 7.777, 8.81, 9.843, 11.071, 12.567, 14.062, 15.611, 17.331, 19.051, 20.824, 22.867, 24.909, 26.932, 10.698, 9.445, 8.517, 8.008, 6.27, 5.394, 5.165, 5.697, 5.397, 5.439, 5.726, 6.351, 7.245, 8.254, 8.855, 9.644, 10.627, 12.086, 13.461, 14.589, 16.17, 17.943, 19.985, 21.719, 23.214, 12.749, 11.496, 10.244, 8.991, 8.158, 7.405, 6.641, 6.109, 5.563, 5.604, 5.774, 5.943, 6.46, 7.061, 7.663, 8.452, 9.435, 10.418, 11.485, 12.613, 13.74, 15.01, 16.505, 18.0, 19.495, ]) +a_fit = np.array([1.993, 1.851, 1.709, 1.587, 1.593, 1.572, 1.526, 1.645, 1.792, 1.938, 1.949, 1.927, 1.862, 1.72, 1.578, 1.582, 1.731, 1.879, 1.908, 1.912, 1.803, 1.831, 1.856, 1.837, 1.708, 1.715, 1.858, 1.889, 1.89, 1.907, 1.91, 1.905, 1.899, 1.893, 1.877, 1.872, 1.879, 1.894, 1.925, 1.955, 1.948, 1.941, 1.909, 1.839, 1.712, 1.717, 1.856, 1.917, 1.925, 1.953, 1.953, 1.92, 1.855, 1.734, 1.648, 1.615, 1.618, 1.739, 1.861, 1.917, 1.944, 1.903, 1.901, 1.86, 1.807, 1.774, 1.777, 1.793, 1.812, 1.858, 2.046, 1.965, 1.884, 1.814, 1.78, 1.792, 1.809, 1.803, 1.815, 1.85, 2.083, 1.975, 1.867, 1.769, 1.732, 1.724, 1.754, 1.79, 1.825, 1.861, 2.126, 2.019, 1.912, 1.815, 1.788, 1.791, 1.826, 1.861, 1.88, 1.891, 2.389, 2.247, 2.106, 1.964, 1.822, 1.59, 1.737, 1.883, 2.029, 2.175, 2.312, 1.989, 1.842, 1.7, 1.558, 1.512, 1.66, 1.809, 1.958, 2.107, 2.18, 2.051, 1.923, 1.795, 1.578, 1.583, 1.616, 1.759, 1.901, 2.044, 2.248, 2.242, 2.237, 2.231, 2.226, 2.028, 2.043, 2.074, 2.105, 2.135, 2.284, 2.274, 2.147, 2.02, 1.893, 1.682, 1.821, 1.959, 2.098, 2.131, 2.309, 2.072, 2.002, 1.922, 1.795, 1.752, 1.89, 1.99, 2.046, 2.102, 2.395, 2.158, 2.088, 2.011, 1.93, 1.871, 1.89, 1.922, 1.978, 2.035, 2.51, 2.401, 2.308, 2.227, 2.146, 1.924, 1.959, 1.982, 2.0, 2.019, 2.548, 2.44, 2.332, 2.224, 2.124, 1.901, 1.906, 1.942, 1.977, 2.013, 2.6, 2.493, 2.386, 2.278, 2.171, 1.97, 2.005, 2.04, 2.074, 2.058, 3.398, 3.256, 3.114, 2.907, 2.316, 1.988, 1.924, 2.04, 2.186, 2.332, 3.167, 2.601, 2.201, 2.032, 1.89, 1.852, 1.798, 1.916, 2.065, 2.214, 2.962, 2.834, 2.712, 2.554, 2.049, 1.78, 1.731, 1.833, 1.972, 2.115, 3.05, 3.044, 3.038, 2.96, 2.466, 2.177, 2.1, 2.096, 2.127, 2.158, 3.103, 3.093, 2.966, 2.788, 2.305, 2.005, 1.887, 1.989, 2.128, 2.157, 3.202, 2.663, 2.29, 2.21, 2.083, 2.083, 2.028, 2.139, 2.262, 2.318, 3.288, 2.749, 2.376, 2.299, 2.218, 2.124, 2.085, 2.109, 2.165, 2.221, 3.408, 3.295, 3.201, 3.049, 2.479, 2.244, 2.143, 2.159, 2.187, 2.205, 3.447, 3.339, 3.231, 3.071, 2.511, 2.195, 2.094, 2.109, 2.144, 2.18, 3.527, 3.42, 3.313, 3.138, 2.564, 2.295, 2.329, 2.364, 2.306, 2.197, 4.406, 4.265, 3.927, 3.335, 2.744, 2.416, 2.352, 2.288, 2.343, 2.489, 3.965, 3.399, 2.897, 2.413, 2.223, 2.184, 2.143, 2.096, 2.152, 2.301, 3.744, 3.631, 3.355, 2.856, 2.372, 2.103, 2.05, 2.0, 2.051, 2.186, 3.85, 3.845, 3.621, 3.131, 2.641, 2.352, 2.264, 2.176, 2.149, 2.18, 3.922, 3.912, 3.63, 3.147, 2.664, 2.364, 2.247, 2.13, 2.157, 2.187, 4.09, 3.55, 3.01, 2.541, 2.414, 2.4, 2.388, 2.333, 2.381, 2.52, 4.182, 3.642, 3.102, 2.848, 2.739, 2.402, 2.363, 2.359, 2.416, 2.427, 4.301, 4.188, 3.881, 3.31, 2.74, 2.521, 2.482, 2.383, 2.359, 2.392, 4.347, 4.242, 3.943, 3.393, 2.851, 2.535, 2.445, 2.349, 2.322, 2.357, 4.454, 4.347, 4.037, 3.463, 2.889, 2.619, 2.654, 2.554, 2.445, 2.335, 5.415, 4.947, 4.355, 3.764, 3.172, 2.844, 2.78, 2.716, 2.653, 2.646, 4.763, 4.197, 3.695, 3.192, 2.689, 2.517, 2.476, 2.435, 2.394, 2.389, 4.55, 4.155, 3.663, 3.179, 2.694, 2.426, 2.373, 2.321, 2.269, 2.268, 4.651, 4.287, 3.797, 3.307, 2.817, 2.528, 2.44, 2.352, 2.264, 2.211, 4.732, 4.38, 3.897, 3.414, 2.929, 2.64, 2.529, 2.412, 2.294, 2.216, 4.909, 4.369, 3.829, 3.36, 2.914, 2.688, 2.689, 2.69, 2.596, 2.549, 5.075, 4.535, 3.995, 3.571, 3.0, 2.69, 2.65, 2.641, 2.642, 2.643, 5.194, 4.749, 4.179, 3.608, 3.056, 2.782, 2.743, 2.704, 2.661, 2.571, 5.255, 4.816, 4.274, 3.732, 3.191, 2.875, 2.785, 2.695, 2.604, 2.535, 5.381, 4.935, 4.362, 3.788, 3.214, 2.912, 2.802, 2.693, 2.583, 2.474, 8.975, 8.384, 7.792, 7.2, 5.279, 4.137, 4.073, 4.009, 3.83, 3.403, 7.542, 7.04, 6.537, 6.035, 4.27, 3.543, 3.114, 3.073, 3.032, 2.975, 5.9, 5.328, 4.844, 4.36, 3.761, 3.232, 2.986, 2.934, 2.881, 2.75, 5.837, 5.26, 4.77, 4.28, 3.79, 3.292, 3.204, 3.116, 3.027, 2.633, 5.856, 5.346, 4.81, 4.312, 3.829, 3.383, 3.026, 2.909, 2.792, 2.703, 6.126, 5.586, 5.007, 4.4, 3.864, 3.507, 3.184, 3.067, 3.011, 3.012, 6.395, 5.767, 5.179, 4.635, 4.095, 3.529, 3.228, 3.224, 3.184, 3.145, 6.515, 5.922, 5.345, 4.774, 4.204, 3.699, 3.338, 3.298, 3.259, 3.22, 6.794, 6.052, 5.411, 4.868, 4.327, 3.775, 3.468, 3.378, 3.288, 3.184, 6.963, 6.389, 5.816, 5.242, 4.395, 3.769, 3.468, 3.358, 3.249, 3.139, 13.224, 12.633, 11.244, 9.025, 6.805, 5.663, 5.599, 5.176, 4.748, 4.321, 9.773, 9.27, 7.927, 6.716, 5.786, 5.058, 4.542, 4.026, 3.752, 3.695, 7.687, 6.819, 5.989, 5.505, 4.906, 4.378, 4.131, 4.079, 3.715, 3.317, 7.058, 6.468, 5.879, 5.334, 4.844, 4.346, 4.258, 3.977, 3.541, 3.106, 7.044, 6.534, 5.998, 5.463, 4.929, 4.483, 4.073, 3.663, 3.402, 3.313, 7.438, 6.817, 6.237, 5.701, 5.1, 4.696, 4.327, 3.958, 3.711, 3.594, 7.774, 7.146, 6.518, 5.902, 5.362, 4.779, 4.41, 4.025, 3.761, 3.722, 7.951, 7.358, 6.779, 6.204, 5.576, 5.039, 4.605, 4.108, 3.844, 3.805, 8.435, 7.72, 7.005, 6.242, 5.678, 5.126, 4.692, 4.278, 4.063, 3.96, 8.722, 8.148, 7.411, 6.503, 5.594, 4.968, 4.624, 4.28, 4.056, 3.947, 18.118, 15.898, 13.678, 11.458, 9.239, 7.637, 7.209, 6.782, 6.354, 5.776, 12.842, 10.288, 9.279, 8.349, 7.419, 6.517, 6.051, 5.535, 5.019, 4.502, 14.47, 8.659, 7.808, 6.939, 6.071, 5.438, 5.039, 4.641, 4.244, 3.946, 13.396, 12.768, 12.179, 11.59, 11.001, 5.286, 4.85, 4.414, 3.978, 3.649, 13.353, 12.818, 12.282, 11.746, 11.185, 5.556, 5.134, 4.724, 4.314, 3.904, 8.416, 7.801, 7.185, 6.583, 6.047, 5.653, 5.284, 4.915, 4.645, 4.235, 8.8, 8.172, 7.544, 7.062, 6.446, 5.941, 5.572, 5.187, 4.689, 4.36, 9.205, 8.49, 7.993, 7.365, 6.738, 6.218, 5.767, 5.269, 4.772, 4.442, 9.701, 8.986, 8.271, 7.556, 6.94, 6.308, 5.875, 5.441, 5.018, 4.808, 10.083, 9.432, 8.524, 7.616, 6.708, 6.082, 5.738, 5.394, 5.05, 4.81, 29.01, 26.391, 23.385, 20.379, 17.373, 15.653, 14.799, 13.944, 13.089, 12.235, 18.526, 16.94, 15.354, 13.767, 11.273, 10.178, 9.121, 8.065, 6.719, 5.592, 59.205, 52.624, 44.353, 28.853, 13.354, 5.407, 5.314, 5.221, 5.122, 4.975, 67.427, 59.3, 43.944, 28.589, 13.233, 5.372, 5.252, 5.131, 5.011, 4.878, 21.532, 21.086, 20.727, 20.369, 13.023, 5.285, 5.247, 5.182, 5.118, 4.883, 6.93, 6.509, 6.276, 6.023, 5.664, 5.302, 5.238, 5.177, 5.124, 5.071, 7.036, 6.607, 6.374, 6.135, 5.852, 5.51, 5.489, 5.468, 5.419, 5.366, 7.312, 6.916, 6.662, 6.378, 6.094, 5.738, 5.669, 5.648, 5.627, 5.606, 7.469, 6.991, 6.78, 6.568, 6.363, 6.059, 5.985, 5.911, 5.837, 5.782, 7.726, 7.546, 7.367, 7.187, 7.008, 6.39, 6.341, 6.293, 6.244, 6.196, 24.715, 22.495, 20.223, 17.217, 14.211, 12.494, 12.067, 11.397, 10.542, 9.688, 16.505, 14.918, 13.11, 10.355, 9.071, 9.102, 8.046, 6.061, 5.475, 4.958, 46.769, 31.592, 29.427, 28.559, 13.405, 5.456, 5.057, 4.832, 4.734, 4.587, 46.119, 45.492, 43.913, 28.558, 13.202, 5.339, 4.903, 4.646, 4.526, 4.393, 21.919, 21.383, 20.859, 20.5, 13.155, 5.417, 5.378, 5.314, 4.927, 4.517, 7.628, 7.013, 6.423, 6.169, 5.81, 5.432, 5.073, 4.88, 4.816, 4.762, 7.734, 7.102, 6.497, 6.265, 6.033, 5.654, 5.285, 5.113, 5.092, 5.057, 8.005, 7.418, 7.009, 6.725, 6.435, 5.932, 5.464, 5.237, 5.216, 5.195, 8.354, 7.639, 6.932, 6.72, 6.515, 6.074, 5.692, 5.545, 5.471, 5.415, 8.655, 8.476, 8.296, 7.801, 6.893, 6.267, 5.923, 5.746, 5.698, 5.649, ]) + +b_fit = np.array([565.765, 565.376, 564.987, 564.823, 566.071, 556.139, 535.028, 532.351, 532.617, 532.883, 578.669, 576.41, 574.756, 574.167, 573.62, 573.617, 573.999, 574.341, 558.817, 565.062, 517.313, 524.018, 528.168, 527.736, 527.074, 526.98, 527.456, 527.698, 523.806, 539.92, 557.156, 556.57, 555.984, 555.528, 562.263, 556.188, 533.21, 533.431, 533.803, 534.175, 586.087, 584.294, 551.596, 531.451, 530.139, 530.008, 531.057, 542.692, 561.383, 548.117, 590.642, 557.945, 537.114, 537.301, 536.431, 542.784, 556.829, 559.667, 560.625, 561.205, 574.747, 568.864, 559.003, 549.172, 541.38, 547.733, 561.777, 563.829, 563.966, 564.431, 561.185, 559.53, 557.874, 555.088, 533.238, 559.297, 559.795, 566.776, 571.04, 571.476, 564.141, 561.69, 559.406, 554.275, 531.312, 526.818, 558.397, 563.923, 564.359, 564.795, 588.007, 585.488, 582.969, 576.589, 546.021, 531.012, 531.561, 532.111, 538.244, 547.55, 4.704, 4.316, 3.927, 3.538, 3.149, 2.955, 3.221, 3.487, 3.753, 4.019, 6.459, 5.993, 5.399, 4.811, 4.222, 3.952, 4.31, 4.692, 5.074, 5.456, 8.532, 7.87, 7.208, 6.545, 5.966, 5.87, 6.057, 6.532, 7.008, 7.484, 12.021, 11.435, 10.848, 10.262, 9.676, 9.277, 9.549, 9.921, 10.293, 10.666, 16.055, 15.441, 14.129, 12.817, 11.505, 10.837, 11.886, 12.935, 13.984, 14.556, 20.693, 19.221, 17.882, 16.547, 15.235, 14.601, 15.65, 16.476, 17.055, 17.634, 26.359, 24.887, 23.547, 22.011, 20.355, 19.163, 19.301, 19.605, 20.184, 20.763, 32.651, 30.367, 28.428, 26.772, 25.117, 23.825, 24.262, 24.48, 24.617, 24.755, 38.725, 36.441, 34.156, 31.872, 29.627, 27.905, 27.909, 28.345, 28.781, 29.218, 45.207, 42.688, 40.169, 37.65, 35.131, 33.731, 34.281, 34.83, 35.379, 35.15, 4.73, 4.341, 3.952, 3.565, 3.188, 2.94, 2.82, 3.03, 3.296, 3.562, 6.293, 5.894, 5.366, 4.75, 4.161, 3.846, 3.712, 4.022, 4.404, 4.786, 8.277, 7.615, 6.968, 6.356, 5.662, 5.283, 5.01, 5.505, 6.016, 6.492, 11.655, 11.069, 10.482, 9.879, 9.136, 8.544, 8.163, 8.291, 8.663, 9.035, 15.745, 15.132, 13.82, 12.547, 11.503, 10.614, 9.881, 10.672, 11.721, 12.59, 20.686, 18.973, 17.393, 16.058, 14.746, 14.372, 13.596, 14.38, 15.342, 15.922, 26.352, 24.639, 23.058, 21.522, 19.866, 17.955, 17.345, 17.541, 18.12, 18.699, 32.748, 30.401, 28.459, 26.676, 24.145, 23.054, 22.024, 22.248, 22.553, 22.69, 38.836, 36.551, 34.267, 31.956, 28.656, 26.877, 25.62, 25.754, 26.19, 26.627, 45.64, 43.121, 40.602, 37.914, 34.226, 32.657, 33.207, 33.756, 32.888, 31.24, 4.755, 4.366, 3.982, 3.605, 3.228, 2.98, 2.86, 2.741, 2.839, 3.105, 6.055, 5.656, 5.205, 4.739, 4.1, 3.785, 3.633, 3.489, 3.653, 4.035, 8.022, 7.394, 6.751, 6.104, 5.562, 5.184, 4.968, 4.712, 4.946, 5.5, 11.276, 10.689, 10.052, 9.293, 8.533, 7.941, 7.517, 7.092, 7.032, 7.404, 15.436, 14.822, 13.627, 12.583, 11.539, 10.65, 9.917, 9.184, 9.457, 10.327, 20.655, 18.942, 17.229, 15.584, 14.272, 14.096, 13.632, 12.856, 13.115, 14.164, 26.362, 24.648, 22.935, 21.258, 19.551, 17.241, 16.63, 16.469, 17.048, 16.857, 32.779, 30.432, 28.109, 25.577, 23.046, 22.068, 21.458, 20.436, 20.235, 20.626, 38.947, 36.698, 33.944, 30.9, 28.098, 26.319, 25.564, 24.567, 24.096, 24.532, 46.073, 43.554, 40.527, 36.839, 33.152, 31.583, 32.132, 30.627, 28.978, 27.33, 4.78, 4.4, 4.023, 3.646, 3.269, 3.021, 2.901, 2.781, 2.662, 2.648, 5.817, 5.418, 4.967, 4.515, 4.063, 3.724, 3.572, 3.419, 3.267, 3.285, 7.82, 7.146, 6.546, 6.004, 5.463, 5.084, 4.868, 4.652, 4.415, 4.386, 10.896, 10.209, 9.449, 8.69, 7.93, 7.338, 6.914, 6.489, 6.065, 5.829, 15.089, 14.171, 13.127, 12.083, 11.092, 10.5, 9.825, 9.092, 8.359, 8.063, 20.345, 18.632, 16.919, 15.274, 14.202, 13.607, 13.239, 12.87, 12.145, 11.901, 26.392, 24.679, 22.966, 20.984, 18.452, 16.752, 16.141, 15.882, 15.513, 15.144, 32.809, 30.329, 27.797, 25.266, 22.561, 20.969, 20.359, 19.748, 19.289, 18.535, 39.166, 35.948, 33.145, 30.342, 27.54, 25.761, 25.006, 24.252, 23.497, 22.437, 46.506, 43.14, 39.452, 35.765, 32.078, 30.013, 28.365, 26.717, 25.069, 23.421, 5.744, 5.368, 4.991, 4.614, 3.855, 3.373, 3.253, 3.133, 2.986, 2.779, 6.356, 5.904, 5.453, 5.001, 4.291, 4.01, 3.708, 3.555, 3.403, 3.274, 7.286, 6.809, 6.268, 5.726, 5.432, 5.208, 4.929, 4.713, 4.497, 4.259, 10.101, 9.976, 9.216, 8.456, 7.697, 7.426, 7.002, 6.577, 6.153, 5.636, 14.12, 13.121, 13.17, 12.432, 11.388, 10.496, 9.759, 9.025, 8.292, 7.856, 18.973, 17.259, 15.601, 14.282, 14.331, 14.015, 13.175, 12.442, 11.895, 11.526, 25.001, 24.791, 22.995, 21.182, 19.468, 17.28, 16.507, 16.108, 15.498, 14.887, 31.418, 31.133, 29.314, 26.782, 24.251, 22.007, 20.749, 20.139, 19.528, 18.918, 38.557, 37.779, 34.966, 32.142, 29.339, 27.144, 25.716, 24.961, 24.207, 22.82, 46.049, 42.362, 38.675, 34.988, 32.828, 31.164, 29.628, 27.98, 26.332, 24.684, 6.963, 6.586, 5.98, 5.135, 4.291, 3.809, 3.689, 3.483, 3.275, 3.068, 6.559, 6.108, 5.464, 5.064, 4.788, 4.508, 4.172, 3.835, 3.593, 3.464, 7.061, 6.805, 6.52, 5.979, 5.684, 5.46, 5.181, 4.965, 4.661, 4.347, 9.099, 9.059, 9.019, 8.658, 7.898, 7.627, 7.203, 6.721, 6.192, 5.663, 12.649, 11.65, 11.7, 11.749, 11.745, 10.854, 10.115, 9.376, 8.641, 8.204, 17.331, 15.731, 14.938, 14.987, 14.877, 14.92, 14.057, 13.193, 12.416, 11.683, 23.286, 23.075, 22.865, 22.696, 20.983, 18.661, 17.798, 16.839, 15.726, 15.116, 29.643, 29.358, 27.759, 25.4, 25.19, 24.114, 22.707, 21.07, 19.956, 19.346, 37.443, 36.763, 36.082, 35.194, 31.775, 29.58, 27.998, 26.417, 25.259, 23.872, 45.352, 41.665, 38.893, 37.076, 35.258, 33.594, 32.083, 30.573, 28.992, 27.344, 8.361, 7.516, 6.672, 5.827, 4.983, 4.402, 4.195, 3.987, 3.78, 3.549, 7.041, 6.06, 5.749, 5.473, 5.197, 4.885, 4.621, 4.284, 3.948, 3.611, 8.699, 6.551, 6.309, 6.053, 5.797, 5.512, 5.198, 4.884, 4.566, 4.357, 9.896, 9.831, 9.791, 9.75, 9.71, 7.516, 6.987, 6.458, 5.929, 5.617, 13.023, 13.072, 13.122, 13.171, 13.179, 11.124, 10.483, 9.744, 9.005, 8.266, 15.498, 15.354, 15.209, 15.097, 15.146, 14.584, 13.72, 12.869, 12.581, 11.842, 20.656, 20.446, 20.431, 20.553, 20.409, 19.649, 18.785, 17.827, 16.189, 15.288, 27.029, 26.349, 26.598, 26.388, 26.178, 25.101, 23.695, 22.057, 20.42, 19.518, 34.867, 34.187, 33.506, 32.826, 32.07, 30.684, 29.102, 27.521, 25.954, 25.245, 42.743, 41.899, 40.082, 38.265, 36.447, 34.783, 33.272, 31.761, 30.251, 29.489, 11.435, 10.512, 9.514, 8.515, 7.516, 6.913, 6.639, 6.365, 6.091, 5.817, 7.95, 7.375, 6.8, 6.225, 5.271, 5.006, 4.598, 4.191, 3.623, 3.109, 24.318, 21.82, 18.652, 12.642, 6.631, 3.472, 3.536, 3.601, 3.663, 3.706, 27.462, 24.546, 18.896, 13.245, 7.595, 4.604, 4.688, 4.772, 4.855, 4.921, 11.612, 11.698, 11.82, 11.942, 9.142, 6.358, 6.744, 6.919, 7.094, 6.82, 8.14, 8.241, 8.579, 8.881, 9.003, 8.64, 8.815, 9.027, 9.309, 9.591, 10.587, 11.072, 11.41, 11.756, 12.154, 11.576, 12.019, 12.461, 12.761, 13.043, 14.046, 14.658, 15.164, 15.562, 15.96, 15.383, 15.798, 16.24, 16.682, 17.124, 18.415, 17.846, 18.506, 19.166, 19.882, 19.574, 19.985, 20.397, 20.809, 21.356, 22.311, 23.245, 24.178, 25.112, 26.045, 24.42, 25.023, 25.627, 26.231, 26.835, 10.223, 9.378, 8.523, 7.525, 6.526, 5.923, 5.715, 5.47, 5.196, 4.922, 7.683, 7.108, 6.441, 5.377, 4.945, 5.17, 4.763, 3.842, 3.538, 3.202, 20.211, 14.324, 13.558, 13.302, 7.427, 4.265, 3.95, 3.851, 3.914, 3.956, 20.534, 20.469, 20.053, 14.402, 8.752, 5.757, 5.228, 5.046, 5.13, 5.195, 13.663, 13.713, 13.767, 13.889, 11.088, 8.305, 8.691, 8.866, 8.188, 7.449, 11.521, 11.376, 11.264, 11.565, 11.687, 10.583, 10.231, 10.178, 10.353, 10.624, 13.967, 14.604, 15.22, 15.559, 15.897, 14.25, 13.387, 13.263, 13.706, 14.077, 18.647, 19.208, 19.674, 20.071, 20.458, 19.096, 17.976, 17.517, 17.959, 18.401, 25.054, 24.374, 22.964, 23.624, 24.341, 23.837, 22.325, 21.983, 22.395, 22.942, 30.369, 31.302, 32.236, 31.978, 30.161, 28.496, 26.986, 26.674, 27.278, 27.882, ]) diff --git a/pyHalo/Halos/tidal_truncation.py b/pyHalo/Halos/tidal_truncation.py index 0689deb..653ce54 100644 --- a/pyHalo/Halos/tidal_truncation.py +++ b/pyHalo/Halos/tidal_truncation.py @@ -321,12 +321,13 @@ def truncation_radius(self, halo_mass, infall_concentration, class TruncationGalacticus(object): - def __init__(self, lens_cosmo): + def __init__(self, lens_cosmo, c_host): """ - :param lens_cosmo: + :param lens_cosmo: an instance of LensCosmo + :param c_host: host halo concentration at z=0.5 """ - + self._chost = c_host self._lens_cosmo = lens_cosmo self._mass_loss_interp = InterpGalacticus() self._tau_mf_interpolation = tau_mf_interpolation() @@ -354,9 +355,8 @@ def truncation_radius_halo(self, halo): halo_mass = halo.mass infall_concentration = halo.c time_since_infall = halo.time_since_infall - log10c = np.log10(infall_concentration) - log10mbound_over_minfall = self._mass_loss_interp(log10c, time_since_infall) + log10mbound_over_minfall = self._mass_loss_interp(log10c, time_since_infall, self._chost) m_bound = halo_mass * 10 ** log10mbound_over_minfall _, rs, r200 = self._lens_cosmo.NFW_params_physical(halo_mass, infall_concentration, @@ -377,7 +377,7 @@ def truncation_radius(self, halo_mass, infall_concentration, :return: """ log10c = np.log10(infall_concentration) - log10mbound_over_minfall = self._mass_loss_interp(log10c, time_since_infall) + log10mbound_over_minfall = self._mass_loss_interp(log10c, time_since_infall, self._chost) m_bound = halo_mass * 10 ** log10mbound_over_minfall _, rs, r200 = self._lens_cosmo.NFW_params_physical(halo_mass, infall_concentration, diff --git a/pyHalo/PresetModels/cdm.py b/pyHalo/PresetModels/cdm.py index c5f6fbc..ae11661 100644 --- a/pyHalo/PresetModels/cdm.py +++ b/pyHalo/PresetModels/cdm.py @@ -17,7 +17,7 @@ def CDM(z_lens, z_source, sigma_sub=0.025, log_mlow=6., log_mhigh=10., log10_sig shmf_log_slope=-1.9, cone_opening_angle_arcsec=6., log_m_host=13.3, r_tidal=0.25, LOS_normalization=1.0, two_halo_contribution=True, delta_power_law_index=0.0, geometry_type='DOUBLE_CONE', kwargs_cosmo=None, host_scaling_factor=0.5, - redshift_scaling_factor=0.3, two_halo_Lazar_correction=True, draw_poisson=True): + redshift_scaling_factor=0.3, two_halo_Lazar_correction=True, draw_poisson=True, c_host=6.0): """ This class generates realizations of dark matter structure in Cold Dark Matter @@ -57,6 +57,7 @@ def CDM(z_lens, z_source, sigma_sub=0.025, log_mlow=6., log_mhigh=10., log10_sig :param redshift_scaling_factor: the scaling with (1+z) of the projected number density of subhalos :param two_halo_Lazar_correction: bool; if True, adds the correction to the two-halo contribution from around the main deflector presented by Lazar et al. (2021) + :param c_host: manually set host halo concentration :return: a realization of dark matter halos """ @@ -90,7 +91,8 @@ def CDM(z_lens, z_source, sigma_sub=0.025, log_mlow=6., log_mhigh=10., log10_sig model_fieldhalos, kwargs_mc_field = preset_concentration_models(concentration_model_fieldhalos, kwargs_concentration_model_fieldhalos) concentration_model_fieldhalos = model_fieldhalos(**kwargs_mc_field) - c_host = concentration_model_fieldhalos.nfw_concentration(10 ** log_m_host, z_lens) + if c_host is None: + c_host = concentration_model_fieldhalos.nfw_concentration(10 ** log_m_host, z_lens) # SET THE TRUNCATION RADIUS FOR SUBHALOS AND FIELD HALOS kwargs_truncation_model_subhalos['lens_cosmo'] = pyhalo.lens_cosmo @@ -98,7 +100,7 @@ def CDM(z_lens, z_source, sigma_sub=0.025, log_mlow=6., log_mhigh=10., log10_sig model_subhalos, kwargs_trunc_subs = truncation_models(truncation_model_subhalos) kwargs_trunc_subs.update(kwargs_truncation_model_subhalos) - if truncation_model_subhalos == 'TRUNCATION_GALACTICUS_KEELEY24': + if truncation_model_subhalos in ['TRUNCATION_GALACTICUS_KEELEY24', 'TRUNCATION_GALACTICUS']: kwargs_trunc_subs['c_host'] = c_host truncation_model_subhalos = model_subhalos(**kwargs_trunc_subs) diff --git a/pyHalo/PresetModels/wdm.py b/pyHalo/PresetModels/wdm.py index f47910d..a2c8c21 100644 --- a/pyHalo/PresetModels/wdm.py +++ b/pyHalo/PresetModels/wdm.py @@ -18,7 +18,7 @@ def WDM(z_lens, z_source, log_mc, sigma_sub=0.025, log_mlow=6., log_mhigh=10., l LOS_normalization=1.0, geometry_type='DOUBLE_CONE', kwargs_cosmo=None, mdef_subhalos='TNFW', mdef_field_halos='TNFW', kwargs_density_profile={}, host_scaling_factor=0.5, redshift_scaling_factor=0.3, two_halo_Lazar_correction=True, - draw_poisson=True): + draw_poisson=True, c_host=None): """ This class generates realizations of dark matter structure in Warm Dark Matter @@ -63,6 +63,7 @@ def WDM(z_lens, z_source, log_mc, sigma_sub=0.025, log_mlow=6., log_mhigh=10., l :param redshift_scaling_factor: the scaling with (1+z) of the projected number density of subhalos :param two_halo_Lazar_correction: bool; if True, adds the correction to the two-halo contribution from around the main deflector presented by Lazar et al. (2021) + :param c_host: manually set host halo concentration :return: a realization of dark matter halos """ # FIRST WE CREATE AN INSTANCE OF PYHALO, WHICH SETS THE COSMOLOGY @@ -104,14 +105,15 @@ def WDM(z_lens, z_source, log_mc, sigma_sub=0.025, log_mlow=6., log_mhigh=10., l concentration_model_CDM = preset_concentration_models('DIEMERJOYCE19')[0] kwargs_mc_field['concentration_cdm_class'] = concentration_model_CDM concentration_model_fieldhalos = model_fieldhalos(**kwargs_mc_field) - c_host = concentration_model_fieldhalos.nfw_concentration(10 ** log_m_host, z_lens) + if c_host is None: + c_host = concentration_model_CDM.nfw_concentration(10 ** log_m_host, z_lens) # SET THE TRUNCATION RADIUS FOR SUBHALOS AND FIELD HALOS kwargs_truncation_model_subhalos['lens_cosmo'] = pyhalo.lens_cosmo kwargs_truncation_model_fieldhalos['lens_cosmo'] = pyhalo.lens_cosmo model_subhalos, kwargs_trunc_subs = truncation_models(truncation_model_subhalos) kwargs_trunc_subs.update(kwargs_truncation_model_subhalos) - if truncation_model_subhalos == 'TRUNCATION_GALACTICUS_KEELEY24': + if truncation_model_subhalos in ['TRUNCATION_GALACTICUS_KEELEY24', 'TRUNCATION_GALACTICUS']: kwargs_trunc_subs['c_host'] = c_host kwargs_trunc_subs['lens_cosmo'] = pyhalo.lens_cosmo truncation_model_subhalos = model_subhalos(**kwargs_trunc_subs) @@ -247,7 +249,7 @@ def WDMGeneral(z_lens, z_source, log_mc, dlogT_dlogk, sigma_sub=0.025, log_mlow= shmf_log_slope=-1.9, cone_opening_angle_arcsec=6., log_m_host=13.3, r_tidal=0.25, LOS_normalization=1.0, geometry_type='DOUBLE_CONE', kwargs_cosmo=None, mdef_subhalos='TNFW', mdef_field_halos='TNFW', kwargs_density_profile={}, - host_scaling_factor=0.5, redshift_scaling_factor=0.3, two_halo_Lazar_correction=True): + host_scaling_factor=0.5, redshift_scaling_factor=0.3, two_halo_Lazar_correction=True, c_host=None): """ This preset model implements a generalized treatment of warm dark matter, or any theory that produces a cutoff in @@ -290,6 +292,7 @@ def WDMGeneral(z_lens, z_source, log_mc, dlogT_dlogk, sigma_sub=0.025, log_mlow= :param redshift_scaling_factor: the scaling with (1+z) of the projected number density of subhalos :param two_halo_Lazar_correction: bool; if True, adds the correction to the two-halo contribution from around the main deflector presented by Lazar et al. (2021) + :param c_host: manually fix the host halo concentration :return: """ # FIRST WE CREATE AN INSTANCE OF PYHALO, WHICH SETS THE COSMOLOGY @@ -320,15 +323,15 @@ def WDMGeneral(z_lens, z_source, log_mc, dlogT_dlogk, sigma_sub=0.025, log_mlow= concentration_model_fieldhalos = model_fieldhalos(**kwargs_concentration_model_fieldhalos) concentration_model_CDM = preset_concentration_models('DIEMERJOYCE19')[0](pyhalo.astropy_cosmo) - c_host = concentration_model_CDM.nfw_concentration(10 ** log_m_host, z_lens) + if c_host is None: + c_host = concentration_model_CDM.nfw_concentration(10 ** log_m_host, z_lens) # SET THE TRUNCATION RADIUS FOR SUBHALOS AND FIELD HALOS kwargs_truncation_model_subhalos['lens_cosmo'] = pyhalo.lens_cosmo kwargs_truncation_model_fieldhalos['lens_cosmo'] = pyhalo.lens_cosmo model_subhalos, kwargs_trunc_subs = truncation_models(truncation_model_subhalos) kwargs_trunc_subs.update(kwargs_truncation_model_subhalos) - if truncation_model_subhalos == 'TRUNCATION_GALACTICUS_KEELEY24': + if truncation_model_subhalos in ['TRUNCATION_GALACTICUS_KEELEY24', 'TRUNCATION_GALACTICUS']: kwargs_trunc_subs['c_host'] = c_host - truncation_model_subhalos = model_subhalos(**kwargs_trunc_subs) model_fieldhalos, kwargs_trunc_field = truncation_models(truncation_model_fieldhalos) kwargs_trunc_field.update(kwargs_truncation_model_fieldhalos) diff --git a/tests/test_halos/test_galacticus_truncation.py b/tests/test_halos/test_galacticus_truncation.py index 2991c24..75da863 100644 --- a/tests/test_halos/test_galacticus_truncation.py +++ b/tests/test_halos/test_galacticus_truncation.py @@ -12,7 +12,7 @@ def setup_method(self): self.lens_cosmo_instance = LensCosmo(0.5, 1.5) chost = 5.0 - self.tg = TruncationGalacticus(self.lens_cosmo_instance) + self.tg = TruncationGalacticus(self.lens_cosmo_instance, chost) self.tgk24 = TruncationGalacticusKeeley24(self.lens_cosmo_instance, chost) def test_interp_k24(self): @@ -113,7 +113,7 @@ def rescale_normalization(self, rescale): def test_tg(self): lens_cosmo = LensCosmo(0.5, 2.0,None) - truncation_class = TruncationGalacticus(lens_cosmo) + truncation_class = TruncationGalacticus(lens_cosmo, c_host=4.0) concentration_class = ConcentrationDiemerJoyce(lens_cosmo.cosmo, scatter=False) mass = 10 ** 8 x = 0 diff --git a/tests/test_halos/test_tnfw_halo.py b/tests/test_halos/test_tnfw_halo.py index 3f57919..0368f6a 100644 --- a/tests/test_halos/test_tnfw_halo.py +++ b/tests/test_halos/test_tnfw_halo.py @@ -92,10 +92,10 @@ def test_density_mass(self): mtheory = 4 * np.pi * rho_s * rs ** 3 * (np.log(1 + c) - c/(1+c)) npt.assert_almost_equal(mtheory, tnfw_fieldhalo.mass) rmax = c * rs + m_calculated = tnfw_fieldhalo.mass_3d(rmax) npt.assert_almost_equal(mtheory/m_calculated, 1.0) - if __name__ == '__main__': pytest.main() diff --git a/tests/test_halos/test_truncation.py b/tests/test_halos/test_truncation.py index c442f58..11d6333 100644 --- a/tests/test_halos/test_truncation.py +++ b/tests/test_halos/test_truncation.py @@ -24,7 +24,7 @@ def test_load_models(self): 'TRUNCATION_GALACTICUS', 'TRUNCATION_GALACTICUS_KEELEY24'] kwargs_model_list = [{}, {'LOS_truncation_factor': 50.}, {'RocheNorm': 1.0, 'm_power': 1./3, 'RocheNu': 2.0/3.0}, {}, - {}, {}, {'rt_arcsec': 1.0}, {}, {'c_host': 5.0}] + {}, {}, {'rt_arcsec': 1.0}, {'c_host': 5.0}, {'c_host': 9.0}] for model,kwargs in zip(model_name_list, kwargs_model_list): mod, kw = truncation_models(model) kwargs.update(kw) From e7e42fe1eb24fa951922c73c92bf27d5629e7e97 Mon Sep 17 00:00:00 2001 From: Daniel Gilman Date: Thu, 5 Sep 2024 01:21:18 -0500 Subject: [PATCH 3/3] fix preset models test --- pyHalo/PresetModels/wdm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyHalo/PresetModels/wdm.py b/pyHalo/PresetModels/wdm.py index a2c8c21..88b61a0 100644 --- a/pyHalo/PresetModels/wdm.py +++ b/pyHalo/PresetModels/wdm.py @@ -106,7 +106,7 @@ def WDM(z_lens, z_source, log_mc, sigma_sub=0.025, log_mlow=6., log_mhigh=10., l kwargs_mc_field['concentration_cdm_class'] = concentration_model_CDM concentration_model_fieldhalos = model_fieldhalos(**kwargs_mc_field) if c_host is None: - c_host = concentration_model_CDM.nfw_concentration(10 ** log_m_host, z_lens) + c_host = concentration_model_fieldhalos.nfw_concentration(10 ** log_m_host, z_lens) # SET THE TRUNCATION RADIUS FOR SUBHALOS AND FIELD HALOS kwargs_truncation_model_subhalos['lens_cosmo'] = pyhalo.lens_cosmo