deleted unused functions and fixed tests

make_sz_cluster.py

@@ -445,73 +445,6 @@ def simulate_submap(M200_dist, z_dist, id_dist=None,
 
     return clusters
 
-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_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
-    y_noise = CMB_noise / fsz / t_cmb / 1e6
-    y_img = make_proj_image_new(radius,profile,extrapolate=True)
-
-    if 1 in nums:
-        pa = p + '1' + '.png'
-        plot_y(radius, profile, redshift_z, pa)
-
-    if 2 in nums:
-        pa = p + '2' + '.png'
-        plot_img(y_img, redshift_z, path = pa)
-
-    if 3 in nums:
-        pa = p + '3' + '.png'
-        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, redshift_z, path = pa)
-    if 5 in nums:
-        pa = p + '5' + '.png'
-        plot_img(cmb_img, redshift_z, opt = 1, path = pa)
-    if 6 in nums:
-        pa = p + '6' + '.png'
-        plot_img(noise, redshift_z, opt = 1, path = pa)
-    if 7 in nums:
-        pa = p + '7' + '.png'
-        plot_img(CMB_noise, redshift_z, opt = 1, path = pa)
-    if 8 in nums:
-        pa = p + '8' + '.png'
-        plot_img(y_noise, redshift_z, path = pa)
-    if 9 in nums:
-        pa = p + '9' + '.png'
-        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
-
 def get_r200_and_c200(cosmo,sigma8,ns,M200_SM,redshift_z):
     '''
     Parameters:
@@ -547,63 +480,4 @@ def get_r200_and_c200(cosmo,sigma8,ns,M200_SM,redshift_z):
     R200_mpc = R200_mpc*cosmo.h/1000 #From kpc/h to Mpc
     #From Mpc proper to Mpc comoving
     R200_mpc = R200_mpc/cosmo.scale_factor(redshift_z) 
-    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
-
-
-def make_proj_image_new(radius, profile,range=18,pixel_scale=0.5,extrapolate=False):
-    '''
-    Input: Profile as function of Radius, range (default to 18) & pixel_scale (default to 0.5) in Mpc
-    Return: 2D profile
-    '''
-    image_size = range/pixel_scale+1
-
-    if extrapolate:
-        profile_spline = interp1d(radius, profile, kind = 3, bounds_error=False, fill_value="extrapolate")
-    else:
-        profile_spline = interp1d(radius, profile, bounds_error=False)
-
-    x,y=np.meshgrid(np.arange(image_size),np.arange(image_size))
-    r = np.sqrt((x-image_size//2)**2+(y-image_size//2)**2)*pixel_scale
-    image = profile_spline(r)
-
-    return image, x, y, r
-
-def battaglia_profile(r, Mvir, z, cosmo): #THIS IS OLD; WILL LIKELY DELETE SOON
-    '''
-    Using Battaglia et al (2012). Eq. 10.
-    Input: Virial Mass in solar mass and Radius in Mpc
-    Return: Pressure profile in keV/cm^3 at radius r
-    '''
-
-    a = cosmo.scale_factor(z)
-    rho_critical = cosmo.critical_density(0.5).to(u.M_sun/u.Mpc**3) #In M_sun/Mpc^3
-
-
-
-    M200, R200, c200 = mass_adv.changeMassDefinitionCModel(Mvir/cosmo.h, z, 'vir', '200c', c_model = 'ishiyama21')
-    #M200, R200, c200 = mass_adv.changeMassDefinitionCModel(Mvir/cosmo_h, z, 'vir', '200c', c_model = 'ishiyama21')
-    #M200, R200, c200 = mass_adv.changeMassDefinitionCModel(M500/cosmo_h, z, '500c', '200c', c_model = 'ishiyama21')
-    #cvir = concentration.concentration(Mvir, 'vir', z, model = 'ishiyama21')      #Ishiyama et al. (2021)
-    #Option to customize concentration, currently default, using Bullock et al. (2001)
-
-    M200 *= cosmo.h
-    R200 = R200 / 1000 * cosmo.h * (1.+z)
-    x = r / R200
-
-    G2 = G * 1e-6 / float(Mpc_to_m) * m_sun
-    gamma = -0.3
-
-    P200 =  G2 * M200 * 200. * rho_critical * (omega_b / omega_m) / 2. / (R200 / (1. + z))    # Msun km^2 / s^2 / Mpc^3
-
-    P0 = 18.1 * ((M200 / 1e14)**0.154 * (1. + z)**-0.758)
-    xc = 0.497 * ((M200 / 1e14)**-0.00865 * (1. + z)**0.731)
-    beta = 4.35 * ((M200 / 1e14)**0.0393 * (1. + z)**0.415)
-    pth = P200 * P0 * (x / xc)**gamma * (1. + (x/xc))**(-1. * beta)      # Msun km^2 / s^2 / Mpc^3
-
-    pth *= (m_sun * 1e6 * j_to_kev  / ((Mpc_to_m*100)**3))       # keV/cm^3
-    p_e = pth * 0.518       # Vikram et al (2016)
-
-    return p_e, M200, R200, c200
+    return M200_SM, R200_mpc, c200

simsz/make_dm_halo.py

@@ -6,17 +6,4 @@ def flatdist_halo(zmin,zmax,m500min,m500max,size):
     zdist=np.random.uniform(low=zmin, high=zmax, size=size)
     mdist=np.random.uniform(low=m500min, high=m500max, size=size)
 
-    return zdist, mdist
-
-def vir_to_200_colossus(cosmo,sigma8,ns,Mvir,z):
-    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(Mvir/cosmo.h, z, 'vir', '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)
+    return zdist, mdist

simsz/read_yaml.py

@@ -26,19 +26,3 @@ def parse_yaml(self):
                 return yaml.safe_load(stream)
             except yaml.YAMLError as exc:
                 print(exc)
-
-    def load_vars(self, ref):
-        survey=[key for key in ref['SURVEYS'].keys()][0]
-        survey_freq=[key for key in ref['SURVEYS'][survey].keys()][0]
-        beam_size=ref['SURVEYS'][survey][survey_freq]['beam_size']
-        noise_level=ref['SURVEYS'][survey][survey_freq]['noise_level']
-        image_size = ref['IMAGES']['image_size']
-        pixel_scale = ref['IMAGES']['pixel_scale']
-        for key in ref['COSMOLOGY'].keys():
-            cosmo_dict=ref['COSMOLOGY'][key] #Read in cosmological parameters
-
-        sigma8=cosmo_dict['sigma8']
-        ns=cosmo_dict['ns']
-
-        cosmo=FlatLambdaCDM(cosmo_dict['H0'], cosmo_dict['Omega_m0'], Tcmb0=cosmo_dict['t_cmb'], Ob0=cosmo_dict['Omega_b0'])
-        return(survey,survey_freq,beam_size,noise_level,image_size,pixel_scale,cosmo,sigma8,ns)

tests/test_szcluster.py

@@ -0,0 +1,82 @@
+import pytest
+from make_sz_cluster import P200_Battaglia2012, param_Battaglia2012, Pth_Battaglia2012
+import simsz.utils as utils
+import numpy as np
+import astropy.units as u
+from astropy.cosmology import FlatLambdaCDM
+
+'''
+Tests for the make_sz_cluster.py file, only testing functions used in simulation of cluster using 
+Battaglia 2012 profile. 
+
+Constants used in testing are similiar to those used in Battaglia 2012 paper
+In order to run tests, use the command `PYTHONPATH=.. pytest` while in the 'tests' directory.
+'''
+
+def get_mock_cosmology():
+    '''
+    Generates a mock cosmology for testing purposes
+    Returns a tuple of an Astropy FlatLambdaCDM cosmology,
+    a sigma8 value, and a ns (the scalar spectral index)
+    '''
+    cosmo = FlatLambdaCDM(70, 0.25, Tcmb0=2.725, Ob0=0.043)
+    sigma8 = 0.8
+    ns = 0.96
+    return (cosmo,sigma8,ns)
+
+class TestSZCluster:
+    def test_P200_Battaglia2012(self):
+        '''
+        Added testing for the method P200_Battaglia2012,
+        which calculates the P200 param as defined in Battaglia 2012
+        M200 and R200 are the mass and radius of the cluster at 200 times the critical density of the universe
+        P200 is the thermal pressure profile of the shell defined by R200
+        '''
+        redshift_z = 0
+        (cosmo,sigma8,ns) = get_mock_cosmology()
+        M200 = 1.3e13
+        R200 = 0.386
+        P200_expected = 0.00014312182 * (u.keV/u.cm**3.)
+        P200_calculated = P200_Battaglia2012(cosmo, redshift_z, M200, R200)
+        assert u.isclose(P200_calculated,P200_expected), f"Expected {P200_expected}, but got {P200_calculated}"
+
+    def test_param_Battaglia2012(self):
+        '''
+        Test for the method param_Battaglia2012,
+        which calculates independent params as defined in Battaglia 2012, Equation 11 
+        These parameters are used to make the profile as described in Eq 10
+        P0 is the normalization factor/amplitude,
+        xc fits for the core-scale
+        beta is a power law index
+        '''
+        redshift_z = 0
+        M200 = 1.3 * 10e13 # in solar masses
+        P0_expected = 18.84628919814473
+        xc_expected = 0.49587336181740654
+        beta_expected = 4.395084514715711
+        assert u.isclose(param_Battaglia2012(18.1,0.154,-0.758,M200,redshift_z), P0_expected), "Incorrect param calculation: P0"
+        assert u.isclose(param_Battaglia2012(0.497,-0.00865,0.731,M200,redshift_z), xc_expected), "Incorrect param calculation: xc"
+        assert u.isclose(param_Battaglia2012(4.35,0.0393,0.415,M200,redshift_z), beta_expected), "Incorrect param calculation: beta"
+
+    def test_Pth_Battaglia2012(self):
+        '''
+        Test for the method Pth_Battaglia2012,
+        which calculates Pth using the battaglia fit profile, Battaglia 2012, Equation 10
+        Pth is the thermal pressure profile normalized over P200
+        P0 is the normalization factor/amplitude,
+        xc fits for the core-scale
+        beta is a power law index
+        M200 and R200 are the mass and radius of the cluster at 200 times the critical density of the universe
+        P200 is the thermal pressure profile of the shell defined by R200
+        '''
+        radii=np.linspace(0.01,10,10000) #Generate a space of radii in arcmin
+        (cosmo, sigma8, ns) = get_mock_cosmology()
+        radii=utils.arcmin_to_Mpc(radii,0.5,cosmo)
+        P0 = 18.84628919814473
+        xc = 0.49587336181740654
+        beta = 4.395084514715711
+        R200 = 0.386
+        x = radii/R200 #As defined in Battaglia 2012
+        Pth_expected = P0 * (x/xc)**(-0.3) * (1 + (x/xc))**(-beta)
+        result = Pth_Battaglia2012(radii,R200,-0.3,1.0,beta,xc,P0)
+        assert np.allclose(result, Pth_expected), f"Expected {Pth_expected}, but got {result}"