From d4a9502128fd94a07a7dfa832310de5ccedb7917 Mon Sep 17 00:00:00 2001 From: Elaine Ran Date: Mon, 14 Aug 2023 18:19:49 -0400 Subject: [PATCH] added units to function definitions --- make_sz_cluster.py | 492 +++++++++++++++++++++++++++++++-------------- 1 file changed, 344 insertions(+), 148 deletions(-) diff --git a/make_sz_cluster.py b/make_sz_cluster.py index b48cb3d..f7a1139 100644 --- a/make_sz_cluster.py +++ b/make_sz_cluster.py @@ -8,7 +8,7 @@ from pixell import enmap, powspec, enplot import camb -from astropy.constants import G, sigma_T, m_e, c, h, k_B +from astropy import constants as c from astropy import units as u class GenerateCluster(): @@ -17,111 +17,180 @@ def __init__(self): self.data =[] - def P200_Battaglia2012(self,cosmo,redshift_z,M200,R200): + def P200_Battaglia2012(self,cosmo,redshift_z,M200_SM,R200_mpc): ''' - Calculates the P200 pressure profile of a cluster, as defined in Battaglia 2012 + Calculates the P200 pressure profile of a cluster, as defined in + Battaglia 2012 Parameters: - R200, the radius of the cluster at 200 times the critical density of the universe - M200, the mass contained within R200 - redshift_z, the redshift of the cluster - cosmo, background cosmology for density calculation - - Returns P200, the thermal pressure profile of the shell defined by R200 + ----------- + cosmo: FlatLambaCDM instance + background cosmology for density calculation + redshift_z: float + the redshift of the cluster (unitless) + M200_SM: float + the mass contained within R200, in units of solar masses + R200_mpc: float + the radius of the cluster at 200 times the critical density of the + universe in units of Mpc + + Returns: + -------- + P200_kevcm3: Quantity instance + the thermal pressure profile of the shell defined by R200 in units + of keV/cm**3 ''' - P200 = G * M200*u.Msun * 200. * cosmo.critical_density(redshift_z) * (cosmo.Ob0/cosmo.Om0) / (2. * R200*u.Mpc) #From Battaglia 2012 - P200=P200.to(u.keV/u.cm**3.) #Unit conversion to keV/cm^3 + GM200 = c.G * M200_SM*u.Msun * 200. * cosmo.critical_density(redshift_z) + fbR200 = (cosmo.Ob0/cosmo.Om0) / (2. * R200_mpc*u.Mpc) #From Battaglia2012 + P200 = GM200 * fbR200 + P200_kevcm3=P200.to(u.keV/u.cm**3.) #Unit conversion to keV/cm^3 - return(P200) + return(P200_kevcm3) - def param_Battaglia2012(self,A0,alpha_m,alpha_z,M200,redshift_z): + def param_Battaglia2012(self,A0,alpha_m,alpha_z,M200_SM,redshift_z): ''' - Calculates independent params as using the formula from Battaglia 2012, Equation 11 - in order for use in the pressure profile defined in Equation 10 + Calculates independent params as using the formula from Battaglia 2012, + Equation 11 in order for use in the pressure profile defined in + Equation 10 Parameters: - A0, normalization factor - alpha_m, power law index for the mass-dependent part of the function - alpha_z, power law index for the redshift dependent part - M200, the mass of the cluster at 200 times the critical density of the universe - redshift_z, the redshift of the cluster - - Returns the parameter A given the formula from Eq 11 + ----------- + A0: float + normalization factor + alpha_m: float + power law index for the mass-dependent part of the function + alpha_z: float + power law index for the redshift dependent part + M200_SM: float + the mass of the cluster at 200 times the critical density of the + universe, in units of solar masses + redshift_z: float + the redshift of the cluster (unitless) + + Returns: + -------- + A: float + the parameter A given the formula from Eq 11 ''' - A = A0 * (M200/1e14)**alpha_m * (1.+redshift_z)**alpha_z + A = A0 * (M200_SM/1e14)**alpha_m * (1.+redshift_z)**alpha_z return(A) - def Pth_Battaglia2012(self,radius,R200,gamma,alpha,beta,xc,P0): + def Pth_Battaglia2012(self,radius_mpc,R200_mpc,gamma,alpha,beta,xc,P0): ''' - Calculates Pth using the battaglia fit profile, Battaglia 2012, Equation 10 - Pth is the thermal pressure profile normalized over P200 + Calculates Pth using the battaglia fit profile, Battaglia 2012, + Equation 10 Pth is the thermal pressure profile normalized over P200 Parameters: - P0 is the normalization factor/amplitude, - xc fits for the core-scale - beta is a power law index - gamma is a fixed paremeter defined by Battaglia 2012 - alpha is a fixed parameter defined by Battaglia 2012 - R200, the radius of the cluster at 200 times the critical density of the universe - radius, the raidus for the pressure to be calculated at - - Returns Pth, the thermal pressure profile normalized over P200 + ----------- + radius_mpc: float + the radius for the pressure to be calculated at, in units of Mpc + R200_mpc: float + the radius of the cluster at 200 times the critical density of the + universe, in units of Mpc + gamma: float + fixed paremeter defined by Battaglia 2012 + alpha: float + fixed parameter defined by Battaglia 2012 + beta: float + power law index + xc: float + fits for the core-scale + P0: float + the normalization factor/amplitude, + + Returns: + -------- + Pth: float + the thermal pressure profile normalized over P200, units of + keV/cm**3 ''' - x=radius/R200 + x=radius_mpc/R200_mpc Pth = P0 * (x/xc)**gamma * (1+(x/xc)**alpha)**(-beta) return(Pth) - def epp_to_y(self, profile, radii, P200, R200, **kwargs): + def epp_to_y(self, profile, radii_mpc, P200_kevcm3, R200_mpc, **kwargs): ''' Converts from an electron pressure profile to a compton-y profile, integrates over line of sight from -1 to 1 Mpc relative to center. - Parameters: - profile, Method to get thermal pressure profile in Kev/cm^3, accepts radius, R200 and **kwargs - radii, the array of radii corresponding to the profile in Mpc - P200, as defined in battaglia2012, needed for normalization of Battaglia profile - R200, the radius of the cluster at 200 times the critical density of the universe - Return: Compton-y profile corresponding to the radii + Parameters: + ----------- + profile: + Method to get thermal pressure profile in Kev/cm^3, accepts radius, + R200 and **kwargs + radii_mpc: array + the array of radii corresponding to the profile in Mpc + P200_kevcm3: float + as defined in battaglia2012, needed for normalization of Battaglia + profile, units of keV/cm**3 + R200_mpc: float + the radius of the cluster at 200 times the critical density of the + universe in units of Mpc + + Return: + ------- + y_pro: array + Compton-y profile corresponding to the radii ''' - radii = radii * u.Mpc - pressure_integrated = np.empty(radii.size) + radii_mpc = radii_mpc * u.Mpc + pressure_integ = np.empty(radii_mpc.size) i = 0 - l_mpc = np.linspace(0,R200,10000) # Get line of sight axis + l_mpc = np.linspace(0,R200_mpc,10000) # Get line of sight axis l = l_mpc * (1 * u.Mpc).to(u.m).value # l in meters keVcm_to_Jm = (1 * u.keV/(u.cm**3.)).to(u.J/(u.m**3.)) - thermal_to_electron_pressure = 1/1.932 # from Battaglia 2012, assumes fully ionized medium - for radius in radii: - #Multiply profile by P200 specifically for Battaglia 2012 profile, since it returns Pth/P200 instead of Pth - th_pressure = profile(np.sqrt(l_mpc**2 + radius.value**2), R200=R200, **kwargs) * P200.value #pressure as a function of l, multiplication by P20 needed ad - th_pressure = th_pressure * keVcm_to_Jm.value # Use multiplication by a precaluated factor for efficiency + thermal_to_electron_pressure = 1/1.932 # from Battaglia 2012, assumes + #fully ionized medium + for radius in radii_mpc: + #Multiply profile by P200 specifically for Battaglia 2012 profile, + # since it returns Pth/P200 instead of Pth + th_pressure = profile(np.sqrt(l_mpc**2 + radius.value**2), + R200_mpc=R200_mpc, **kwargs) + th_pressure = th_pressure * P200_kevcm3.value #pressure as a + # function of l + th_pressure = th_pressure * keVcm_to_Jm.value # Use multiplication + # by a precaluated factor for efficiency pressure = th_pressure* thermal_to_electron_pressure - integral = np.trapz(pressure, l) * 2 #integrate over pressure in J/m^3 to get J/m^2, multiply by factor of 2 to get from -R200 to R200 (assuming spherical symmetry) - pressure_integrated[i] = integral + integral = np.trapz(pressure, l) * 2 #integrate over pressure in + #J/m^3 to get J/m^2, multiply by factor of 2 to get from -R200 to + # R200 (assuming spherical symmetry) + pressure_integ[i] = integral i += 1 - y_pro = pressure_integrated * sigma_T.value/ (m_e.value * c.value**2) + y_pro = pressure_integ * c.sigma_T.value/ (c.m_e.value * c.c.value**2) return y_pro - def make_y_submap(self, profile, redshift_z, cosmo, width, pix_size, *args, **kwargs): + def make_y_submap(self, profile, redshift_z, cosmo, width, pix_size_arcmin, + *args, **kwargs): ''' Converts from an electron pressure profile to a compton-y profile, integrates over line of sight from -1 to 1 Mpc relative to center. - Parameters: - profile, Method to get thermal pressure profile in Kev/cm^3, accepts radius, and **kwargs - redshift_z, the redshift of the cluster - cosmo, background cosmology for density calculation - width, num pixels to each side of center; end shape of submap will be (2*width +1, 2*width +1) - pix_size, size of each pixel in arcmin - - Return: Compton-y submap with shape (2*width +1, 2*width +1) + Parameters: + ----------- + profile: + Method to get thermal pressure profile in Kev/cm^3, accepts radius, + R200 and **kwargs + redshift_z: float + the redshift of the cluster (unitless) + cosmo: FlatLambaCDM instance + background cosmology for density calculation + width: float + num pixels to each side of center; end shape of submap will be + (2*width +1, 2*width +1) + pix_size_arcmin: float + size of each pixel in arcmin + + Return: + ------- + y_map: array + Compton-y submap with shape (2*width +1, 2*width +1) ''' - X = np.arange(-width, width + pix_size, pix_size) + X = np.arange(-width, width + pix_size_arcmin, pix_size_arcmin) X = utils.arcmin_to_Mpc(X, redshift_z, cosmo) X[X==0] = 0.001 #Solves issues of div by 0 Y = np.transpose(X) @@ -134,103 +203,179 @@ def make_y_submap(self, profile, redshift_z, cosmo, width, pix_size, *args, **kw R.append(np.sqrt(i**2 + j**2)) R = np.array(R) - cy = self.epp_to_y(profile, R, **kwargs) #evaluate compton-y for each neccesary radius + cy = self.epp_to_y(profile, R, **kwargs) #evaluate compton-y for each + # neccesary radius for i in range(X.size): for j in range(Y.size): - y_map[i][j] = cy[np.where(R == np.sqrt(X[i]**2 + Y[j]**2))[0]][0] # assign the correct compton-y to the radius + y_map[i][j] = cy[np.where( + R == np.sqrt(X[i]**2 + Y[j]**2))[0]][0] + # assign the correct compton-y to the radius return y_map - def f_sz(self, freq, T_CMB, *args, **kwargs): + def f_sz(self, freq_ghz, T_CMB_K, *args, **kwargs): ''' - Input: Observation frequency f in GHz, Temperature of cmb T_CMB - Return: Radiation frequency + Parameters: + ---------- + freq_ghz: float + Observation frequency f, in units of GHz + T_CMB_K: Quantity instance + Temperature of CMB in K + + Return: + ------ + fsz: float + radiation frequency ''' - f=freq*u.GHz #Takes input in units of GHz + f=freq_ghz*u.GHz #Takes input in units of GHz f=f.to(1/u.s) #Unit conversion - x = h * f / k_B / T_CMB + x = c.h * f / c.k_B / T_CMB_K fsz = x * (np.exp(x) + 1) / (np.exp(x) - 1) - 4 return fsz - def generate_y_submap(self, redshift_z, M200, R200, cosmo, width, pix_size, profile="Battaglia2012"): + def generate_y_submap(self, redshift_z, M200_SM, R200_mpc, cosmo, + width, pix_size_arcmin, profile="Battaglia2012"): ''' Converts from an electron pressure profile to a compton-y profile, integrates over line of sight from -1 to 1 Mpc relative to center. - Parameters: - - profile, Name of profile, currently only supports "Battaglia2012" - redshift_z, the redshift of the cluster - cosmo, background cosmology for density calculation - width, num pixels to each side of center; end shape of submap will be (2*width +1, 2*width +1) - pix_size, size of each pixel in arcmin - R200, the radius of the cluster at 200 times the critical density of the universe - M200, the mass contained within R200 - Return: Compton-y submap with dimension (2*width +1 , 2*width +1) + Parameters: + ---------- + redshift_z: float + the redshift of the cluster (unitless) + M200_SM: + the mass contained within R200 in solar masses + R200_mpc: float + the radius of the cluster at 200 times the critical density of the + universe in Mpc + cosmo: FlatLambaCDM instance + background cosmology for density calculation + width: float + num pixels to each side of center; end shape of submap will be + (2*width +1, 2*width +1) + pix_size_arcmin: float + size of each pixel in arcmin + profile: str + Name of profile, currently only supports "Battaglia2012" + + Return: + ------ + y_map: array + Compton-y submap with dimension (2*width +1 , 2*width +1) ''' if profile != "Battaglia2012": return None - P200 = self.P200_Battaglia2012(cosmo,redshift_z,M200,R200) #P200 from Battaglia et al. 2012 - P0=self.param_Battaglia2012(18.1,0.154,-0.758,M200,redshift_z) #Parameter computation from Table 1 Battaglia et al. 2012 - xc=self.param_Battaglia2012(0.497,-0.00865,0.731,M200,redshift_z) - beta=self.param_Battaglia2012(4.35,0.0393,0.415,M200,redshift_z) - y_map = self.make_y_submap(self.Pth_Battaglia2012, redshift_z, cosmo, width, pix_size, R200=R200, gamma=-0.3,alpha=1.0,beta=beta,xc=xc,P0=P0, P200=P200) + P200 = self.P200_Battaglia2012(cosmo,redshift_z,M200_solarmasses, + R200_Mpc) #P200 from Battaglia 2012 + P0=self.param_Battaglia2012(18.1,0.154,-0.758,M200_solarmasses, + redshift_z) #Parameter computation from + #Table 1 Battaglia et al. 2012 + xc=self.param_Battaglia2012(0.497,-0.00865,0.731,M200_solarmasses, + redshift_z) + beta=self.param_Battaglia2012(4.35,0.0393,0.415,M200_solarmasses, + redshift_z) + y_map = self.make_y_submap(self.Pth_Battaglia2012, redshift_z, cosmo, + width, pix_size_arcmin, R200_mpc=R200_mpc, + gamma=-0.3,alpha=1.0,beta=beta,xc=xc,P0=P0, + P200_kevcm3=P200) return y_map ####Functions needed in this file: # 3) Convolve submap with beam - def convolve_map_with_gaussian_beam(self, pix_size, beam_size_fwhp, map): + def convolve_map_with_gaussian_beam(self, pix_size_arcmin, + beam_size_fwhp_arcmin, map): ''' - Input: pixel size, beam size in arcmin, image - Return: convolved map + Parameters: + ---------- + pix_size_arcmin: float + size of each pixel in arcmin + beam_size_fwhp_arcmin: float + beam size in arcmin + map: array + image to apply beam convolution to + + Return: + ------- + convolved_map: array Note - pixel size and beam_size need to be in the same units ''' - gaussian = utils.gaussian_kernal(pix_size, beam_size_fwhp) + gaussian = utils.gaussian_kernal(pix_size_arcmin, beam_size_fwhp_arcmin) convolved_map = scipy.signal.fftconvolve(map, gaussian, mode = 'same') return(convolved_map) - def add_cmb_map_and_convolve(self, dT_map, ps, pix_size, beam_size_fwhp_arcmin): + def add_cmb_map_and_convolve(self, dT_map_uK, ps, pix_size_armin, + beam_size_fwhp_arcmin): ''' Parameters: - - dT_map, the map to add to the CMB, in units of -uK - ps, power spectrum with shape (3, 3, lmax); clTT spectrum at ps[0][0] - pix_size, size of each pixel in arcmin - beam_size_fwhp_arcmin, beam size in arcmin - - Return: dT submap with same shape as dT_map, in units of -uK + ---------- + dT_map_uK: array + the map to add to the CMB, units of -uK + ps: array + power spectrum with shape (3, 3, lmax); clTT spectrum at ps[0][0] + pix_size_arcmin: float + size of each pixel in arcmin + beam_size_fwhp_arcmin: float + beam size in arcmin + + Return: + ------ + dT submap: array + dT submap with same shape as dT_map, in units of -uK ''' - padding_value = int(np.ceil(beam_size_fwhp_arcmin/pix_size)) - expanded_shape = (dT_map.shape[0] + 2*padding_value, dT_map.shape[1]+2*padding_value) + padding_value = int(np.ceil(beam_size_fwhp_arcmin/pix_size_arcmin)) + expanded_shape = (dT_map_uK.shape[0] + 2*padding_value, + dT_map_uK.shape[1]+2*padding_value) #print(expanded_shape) - cmb_map = self.make_cmb_map(shape=expanded_shape, pix_size=pix_size, ps=ps) - if type(dT_map) is u.Quantity: + cmb_map = self.make_cmb_map(shape=expanded_shape, + pix_size_armin=pix_size_arcmin, ps=ps) + if type(dT_map_uK) is u.Quantity: cmb_map = cmb_map *u.uK - dT_map_expanded = np.pad(dT_map, (padding_value,padding_value), constant_values=0) + dT_map_expanded = np.pad(dT_map_uK, (padding_value,padding_value), + constant_values=0) signal_map = dT_map_expanded - cmb_map - conv_map = self.convolve_map_with_gaussian_beam(pix_size, beam_size_fwhp_arcmin, signal_map) + conv_map = self.convolve_map_with_gaussian_beam( + pix_size_arcmin, beam_size_fwhp_arcmin, signal_map) - return conv_map[padding_value:-padding_value, padding_value:-padding_value] + return conv_map[padding_value:-padding_value, + padding_value:-padding_value] # 5) Generate noise map - def generate_cluster(self, radius, profile, f, noise_level, beam_size, z, nums, p = None): #SOME OF THE ABOVE WILL BE TAKEN FROM THIS OUTDATED FUNCTION + def generate_cluster(self, radius, profile, f_ghz, noise_level, + beam_size_arcmin, redshift_z, nums, p = None): """ combine all elements to generate a cluster object + + Parameters: + ----------- + radius: float + profile: + f: float + Observation frequency f, in units of GHz + noise_level: float + beam_size_arcmin: float + beam size in arcmin + redshift_z: float + redshift (unitless) + nums: list or array + + Returns: + ------- + a cluster object """ - y_con = convolve_map_with_gaussian_beam(0.5, beam_size , y_img) - fsz = f_sz(f, t_cmb) + y_con = convolve_map_with_gaussian_beam(0.5, beam_size_arcmin , y_img) + fsz = f_sz(f_ghz, t_cmb) cmb_img = y_con * fsz * t_cmb * 1e6 noise = np.random.normal(0, 1, (37, 37)) * noise_level CMB_noise = cmb_img + noise @@ -239,35 +384,36 @@ def generate_cluster(self, radius, profile, f, noise_level, beam_size, z, nums, if 1 in nums: pa = p + '1' + '.png' - plot_y(radius, profile, z, pa) + plot_y(radius, profile, redshift_z, pa) if 2 in nums: pa = p + '2' + '.png' - plot_img(y_img, z, path = pa) + plot_img(y_img, redshift_z, path = pa) if 3 in nums: pa = p + '3' + '.png' - gaussian = gaussian_kernal(0.5, beam_size) - plot_img(gaussian, z, opt = 2, path = pa) + gaussian = gaussian_kernal(0.5, beam_size_arcmin) + plot_img(gaussian, redshift_z, opt = 2, path = pa) if 4 in nums: pa = p + '4' + '.png' - plot_img(y_con, z, path = pa) + plot_img(y_con, redshift_z, path = pa) if 5 in nums: pa = p + '5' + '.png' - plot_img(cmb_img, z, opt = 1, path = pa) + plot_img(cmb_img, redshift_z, opt = 1, path = pa) if 6 in nums: pa = p + '6' + '.png' - plot_img(noise, z, opt = 1, path = pa) + plot_img(noise, redshift_z, opt = 1, path = pa) if 7 in nums: pa = p + '7' + '.png' - plot_img(CMB_noise, z, opt = 1, path = pa) + plot_img(CMB_noise, redshift_z, opt = 1, path = pa) if 8 in nums: pa = p + '8' + '.png' - plot_img(y_noise, z, path = pa) + plot_img(y_noise, redshift_z, path = pa) if 9 in nums: pa = p + '9' + '.png' - plot_img(y_noise, z, opt = 3, path = pa) #vizualization starts working from z = 0.115 + plot_img(y_noise, redshift_z, opt = 3, path = pa) #vizualization + #starts working from z = 0.115 #return tSZ_signal(z, y_con), tSZ_signal(z, y_noise) return y_img, y_con, cmb_img, noise, cmb_noise, y_noise, SZsignal, aperture @@ -275,16 +421,24 @@ def generate_cluster(self, radius, profile, f, noise_level, beam_size, z, nums, def get_cls(self, ns, cosmo, lmax=2000): ''' Makes a cmb temperature map based on the given power spectrum - Parameters: - ns, scalar spectral index of the power-spectrum - cosmo, background cosmology - Return: power spectrum as can be used in szluster.make_cmb_map + Parameters: + ---------- + ns: float + scalar spectral index of the power-spectrum + cosmo: FlatLambaCDM instance + background cosmology + + Return: + ------ + ps array + power spectrum as can be used in szluster.make_cmb_map ''' - data = camb.set_params(ns=ns, H0=cosmo.H0.value, ombh2=cosmo.Ob0, omch2=cosmo.Om0, lmax=lmax, - WantTransfer=True) + data = camb.set_params(ns=ns, H0=cosmo.H0.value, ombh2=cosmo.Ob0, + omch2=cosmo.Om0, lmax=lmax, WantTransfer=True) results = camb.get_results(data) - cls = np.swapaxes(results.get_total_cls(CMB_unit='muK', raw_cl=True),0,1) + cls = np.swapaxes(results.get_total_cls(CMB_unit='muK', raw_cl=True), + 0,1) ps = np.zeros((3,3,cls.shape[1])) # Needs to be reshaped to match input for pixell.enmap ps[0][0]= cls[0] # clTT spectrum @@ -294,18 +448,26 @@ def get_cls(self, ns, cosmo, lmax=2000): ps[2][2] = cls[2] #clBB return ps - def make_cmb_map(self, shape, pix_size, ps): + def make_cmb_map(self, shape, pix_size_arcmin, ps): ''' - Makes a flat-sky CMB temperature map based on the given power spectrum - Parameters: - shape, shape of submap in arcmin - pix_size, size of each pixel in arcmin - ps, power spectrum with shape (3, 3, lmax); clTT spectrum at ps[0][0] + Makes a cmb temperature map based on the given power spectrum - Return: cmb T map + Parameters: + ---------- + shape: tuple + shape of submap in arcmin + pix_size_arcmin: float + size of each pixel in arcmin + ps: array + power spectrum with shape (3, 3, lmax); clTT spectrum at ps[0][0] + + Return: + ------- + cmb T map: array ''' #ps[0][0] is cltt spectrum - shape,wcs = enmap.geometry(shape=shape,pos=(0,0),res=np.deg2rad(pix_size/60.)) + shape,wcs = enmap.geometry(shape=shape,pos=(0,0), + res=np.deg2rad(pix_size_arcmin/60.)) shape = (3,) + shape omap = enmap.rand_map(shape,wcs,cov=ps) #omap gives TQU maps, so for temperature, we need omap[0] @@ -314,37 +476,71 @@ def make_cmb_map(self, shape, pix_size, ps): #other helpful functions - def get_tSZ_signal(self, Map, radmax, fmax=np.sqrt(2)): + def get_tSZ_signal_aperture_photometry(self, Map, radmax_arcmin, + fmax_arcmin=np.sqrt(2)): """ Parameters: - Map - radmax: the radius of the inner radius - - Returns: The average value within an annulus of inner radius R, outer radius sqrt(2)*R, and the tSZ signal + ---------- + Map: array + radmax_arcmin: float + the radius of the inner aperture, in arcmin + fmax_arcmin: float + fmax+radmax is the radius of the outer aperture, in arcmin + + Returns: + ------- + disk_mean: float + The average value within an annulus of inner radius R + ring_mean: float + The average value within an annulus of outer radius sqrt(2)*R + tSZ signal: float + thermal SZ effect signal """ center = np.array(Map.shape) / 2 x, y = np.indices(Map.shape) r = np.sqrt((x - center[0])**2 + (y - center[1])**2) - radius_out=radmax * fmax #define outer radius - disk_mean = Map[r < radmax].mean() #average in inner radius - ring_mean = Map[(r >= radmax) & (r < radius_out)].mean() #average in outer radius + radius_out=radmax_arcmin * fmax_arcmin #define outer radius + disk_mean = Map[r < radmax_arcmin].mean() #average in inner radius + ring_mean = Map[(r >= radmax_arcmin) & (r < radius_out)].mean() #average in outer radius tSZ = disk_mean - ring_mean return disk_mean, ring_mean, tSZ - def m200_to_r200(self,cosmo,sigma8,ns,M200,z): + def get_r200_and_c200(self,cosmo,sigma8,ns,M200_SM,redshift_z): + ''' + Parameters: + ---------- + cosmo: FlatLambaCDM instance + background cosmology + sigma8: float + ns: float + M200_SM: float + the mass contained within R200, in units of solar masses + redshift_z: float + redshift of the cluster (unitless) + + Returns: + ------- + M200_SM: float + the mass contained within R200, in units of solar masses + R200_mpc: float + the radius of the cluster at 200 times the critical density of the universe in Mpc + c200: float + concentration parameter + ''' + params = {'flat': True, 'H0': cosmo.H0.value, 'Om0': cosmo.Om0, 'Ob0': cosmo.Ob0, 'sigma8':sigma8, 'ns': ns} cosmology.addCosmology('myCosmo', **params) cosmo_colossus= cosmology.setCosmology('myCosmo') - M200, R200, c200 = mass_adv.changeMassDefinitionCModel(M200/cosmo.h, z, '200c', '200c', c_model = 'ishiyama21') - M200 *= cosmo.h #From M_solar/h to M_solar - R200 = R200*cosmo.h/1000 #From kpc/h to Mpc - R200 = R200/cosmo.scale_factor(z) #From Mpc proper to Mpc comoving - return M200, R200, c200 + M200_SM, R200_mpc, c200 = mass_adv.changeMassDefinitionCModel(M200_SM/cosmo.h, z, '200c', '200c', c_model = 'ishiyama21') + M200_SM *= cosmo.h #From M_solar/h to M_solar + R200_mpc = R200_mpc*cosmo.h/1000 #From kpc/h to Mpc + R200_mpc = R200_mpc/cosmo.scale_factor(z) #From Mpc proper to Mpc comoving + return M200_SM, R200_mpc, c200 #FUNCTIONS BELOW HERE HAVE NOT BEEN TESTED OR USED RECENTLY; MIGHT BE USEFUL FOR THE ABOVE TO-DO LIST