diff --git a/pyHalo/Halos/HaloModels/gaussian.py b/pyHalo/Halos/HaloModels/gaussianhalo.py
similarity index 86%
rename from pyHalo/Halos/HaloModels/gaussian.py
rename to pyHalo/Halos/HaloModels/gaussianhalo.py
index e717f92..49e4124 100644
--- a/pyHalo/Halos/HaloModels/gaussian.py
+++ b/pyHalo/Halos/HaloModels/gaussianhalo.py
@@ -2,7 +2,7 @@
 from lenstronomy.LensModel.Profiles.gaussian import Gaussian
 import numpy as np
 
-class Gaussian(Halo):
+class GaussianHalo(Halo):
     """
     The base class for a Gaussian fluctuation
 
@@ -19,8 +19,8 @@ def __init__(self, mass, x, y, r3d, z,
         self._concentration_class = concentration_class
         self._truncation_class = truncation_class
         mdef = 'GAUSSIAN_KAPPA'
-        super(Gaussian, self).__init__(mass, x, y, r3d, mdef, z, sub_flag,
-                                            lens_cosmo_instance, args, unique_tag, fixed_position=True)
+        super(GaussianHalo, self).__init__(mass, x, y, r3d, mdef, z, sub_flag,
+                                           lens_cosmo_instance, args, unique_tag, fixed_position=True)
 
     @property
     def profile_args(self):
diff --git a/pyHalo/Halos/HaloModels/powerlaw.py b/pyHalo/Halos/HaloModels/powerlaw.py
index 037c6f1..7fc84bd 100644
--- a/pyHalo/Halos/HaloModels/powerlaw.py
+++ b/pyHalo/Halos/HaloModels/powerlaw.py
@@ -116,3 +116,75 @@ def z_eval(self):
         Returns the redshift at which to evaluate the concentration-mass relation
         """
         return self.z
+
+class GlobularCluster(Halo):
+
+    def __init__(self, mass, x, y, z, lens_cosmo_instance, args, unique_tag):
+        """
+        A mass profile model for a globular cluster
+        :param mass: the total mass of the GC
+        :param x: x coordinate [arcsec]
+        :param y: y coordinate [arcsec]
+        :param z: redshift
+        :param lens_cosmo_instance: an instance of LensCosmo class
+        :param args: keyword arguments for the profile, must contain 'gamma', 'r_core_fraction', 'gc_size_lightyear'
+        which are the logarithmic profile slope outside the core, the core size in units of the GC size, and the gc size
+        in light years. The size and mass are related by
+        gc_size = gc_size_lightyear (m / 10^5)^1/3
+        The mass is the total mass inside a sphere with radius gc_size
+        :param unique_tag:
+        """
+        self._prof = SPLCORE()
+        self._lens_cosmo = lens_cosmo_instance
+        mdef = 'SPL_CORE'
+        super(GlobularCluster, self).__init__(mass, x, y, None, mdef, z, True,
+                                              lens_cosmo_instance, args, unique_tag)
+
+    @property
+    def lenstronomy_params(self):
+        """
+        See documentation in base class (Halos/halo_base.py)
+        """
+        if not hasattr(self, '_lenstronomy_args'):
+
+            profile_args = self.profile_args
+            rho0 = profile_args['rho0']
+            r_core_arcsec = profile_args['r_core_arcsec']
+            gamma = profile_args['gamma']
+            sigma0 = rho0 * r_core_arcsec
+            self._lenstronomy_args = [{'sigma0': sigma0, 'gamma': gamma, 'center_x': self.x, 'center_y': self.y,
+                                       'r_core': r_core_arcsec}]
+        return self._lenstronomy_args, None
+
+    @property
+    def profile_args(self):
+        """
+        See documentation in base class (Halos/halo_base.py)
+        """
+        if not hasattr(self, '_profile_args'):
+            kpc_per_arcsec = self._lens_cosmo.cosmo.kpc_proper_per_asec(self.z)
+            gamma = self._args['gamma']
+            r_core_fraction = self._args['r_core_fraction']
+            gc_size_lightyear_0 = self._args['gc_size_lightyear']
+            sigma_crit_mpc = self._lens_cosmo.get_sigma_crit_lensing(self.z, self._lens_cosmo.z_source)
+            sigma_crit_arcsec = sigma_crit_mpc * (0.001 * kpc_per_arcsec) ** 2
+            kpc_per_lightyear = 0.3 * 1e-3
+            gc_size = gc_size_lightyear_0 * (self.mass / 10 ** 5) ** (1 / 3) * kpc_per_lightyear
+            r_core_arcsec = r_core_fraction * gc_size
+            rho0 = self.mass / self._prof.mass_3d(gc_size, sigma_crit_arcsec, r_core_arcsec, gamma)
+            self._profile_args = {'rho0': rho0, 'gc_size': gc_size, 'gamma': gamma, 'r_core_arcsec': r_core_arcsec}
+        return self._profile_args
+
+    @property
+    def lenstronomy_ID(self):
+        """
+        See documentation in base class (Halos/halo_base.py)
+        """
+        return ['SPL_CORE']
+
+    @property
+    def z_eval(self):
+        """
+        Returns the redshift at which to evalate the concentration-mass relation
+        """
+        return self.z
diff --git a/pyHalo/Halos/lens_cosmo.py b/pyHalo/Halos/lens_cosmo.py
index 426ac3e..1667a85 100644
--- a/pyHalo/Halos/lens_cosmo.py
+++ b/pyHalo/Halos/lens_cosmo.py
@@ -40,14 +40,12 @@ def __init__(self, z_lens=None, z_source=None, cosmology=None,
         if cosmology is None:
             from pyHalo.Cosmology.cosmology import Cosmology
             cosmology = Cosmology()
-
         self.cosmo = cosmology
         self._arcsec = 2 * numpy.pi / 360 / 3600
         self.h = self.cosmo.h
         # critical density of the universe in M_sun h^2 Mpc^-3
         rhoc = un.Quantity(self.cosmo.astropy.critical_density(0), unit=un.Msun / un.Mpc ** 3).value
         self.rhoc = rhoc / self.cosmo.h ** 2
-
         if z_lens is not None and z_source is not None:
             # critical density for lensing in units M_sun * Mpc ^ -2
             self.sigma_crit_lensing = self.get_sigma_crit_lensing(z_lens, z_source)
diff --git a/pyHalo/PresetModels/wdm.py b/pyHalo/PresetModels/wdm.py
index 7c9d96c..2023edf 100644
--- a/pyHalo/PresetModels/wdm.py
+++ b/pyHalo/PresetModels/wdm.py
@@ -49,8 +49,8 @@ def WDM(z_lens, z_source, log_mc, sigma_sub=0.025, log_mlow=6., log_mhigh=10., l
     :param truncation_model_fieldhalos: the truncation model applied to field halos, see truncation_models for a
     complete list
     :param kwargs_truncation_model_fieldhalos: keyword arguments for the truncation model applied to field halos
-    :param infall_redshift_model: a string that specifies that infall redshift sampling distribution, see the LensCosmo
-    class for details
+    :param infall_redshift_model: a string that specifies that infall redshift sampling distribution, options are
+    HYBRID_INFALL (accounts for subhalos and sub-subhalos) and DIRECT_INFALL (only subhalos)
     :param kwargs_infall_model: keyword arguments for the infall redshift model
     :param shmf_log_slope: the logarithmic slope of the subhalo mass function pivoting around 10^8 M_sun
     :param cone_opening_angle_arcsec: the opening angle of the rendering volume in arcsec
diff --git a/pyHalo/Rendering/MassFunctions/gaussian.py b/pyHalo/Rendering/MassFunctions/gaussian.py
new file mode 100644
index 0000000..01e2359
--- /dev/null
+++ b/pyHalo/Rendering/MassFunctions/gaussian.py
@@ -0,0 +1,41 @@
+from copy import deepcopy
+import numpy as np
+
+
+class Gaussian(object):
+    """
+    This class generates masses from a delta function normalized with respect to a
+    background density, a mass, and a volume
+
+    """
+    name = 'GAUSSIAN'
+    def __init__(self, n, mean, sigma, draw_poisson=True, *args, **kwargs):
+        """
+        :param n: normalization, also equal to the number of objects
+        :param mean: mass of objects to render
+        :param sigma: rendering volume
+        :param draw poisson: whether or not to draw from a poisson distribution
+        """
+        self.mean = mean
+        self.sigma = sigma
+        self.n_mean = n
+        self.first_moment = self.mean
+        self.draw_poisson = draw_poisson
+
+    def draw(self):
+
+        """
+        :return: an array of masses
+        """
+
+        if self.draw_poisson:
+            n = int(np.random.poisson(self.n_mean))
+        else:
+            n = int(np.round(self.n_mean))
+
+        if n > 0:
+            log10_mass = np.random.normal(self.mean, self.sigma, n)
+            m = 10 ** log10_mass
+            return np.array(m)
+        else:
+            return np.array([])
diff --git a/pyHalo/Rendering/MassFunctions/stucker.py b/pyHalo/Rendering/MassFunctions/stucker.py
index 0cc5a1f..d577587 100644
--- a/pyHalo/Rendering/MassFunctions/stucker.py
+++ b/pyHalo/Rendering/MassFunctions/stucker.py
@@ -64,3 +64,4 @@ def _supp_mscale(a, b, c, frac=0.5):
         # This is just here, to tell the fitting function when it goes out of bounds
         return np.nan
     return a / (frac**(1./c) - 1.)**(1./b)
+
diff --git a/pyHalo/concentration_models.py b/pyHalo/concentration_models.py
index 5c41f97..a9af0e5 100644
--- a/pyHalo/concentration_models.py
+++ b/pyHalo/concentration_models.py
@@ -15,7 +15,6 @@ def preset_concentration_models(model_name, kwargs_model=None):
         kwargs_model_return = {}
     else:
         kwargs_model_return = deepcopy(kwargs_model)
-
     if model_name == 'DIEMERJOYCE19':
         return ConcentrationDiemerJoyce, kwargs_model_return
     elif model_name == 'LUDLOW2016':
diff --git a/pyHalo/realization_extensions.py b/pyHalo/realization_extensions.py
index 4ac1e44..9bf2647 100644
--- a/pyHalo/realization_extensions.py
+++ b/pyHalo/realization_extensions.py
@@ -1,10 +1,11 @@
 import numpy as np
-from pyHalo.Halos.HaloModels.powerlaw import PowerLawSubhalo, PowerLawFieldHalo
+from pyHalo.Halos.HaloModels.powerlaw import PowerLawSubhalo, PowerLawFieldHalo, GlobularCluster
 from pyHalo.Halos.HaloModels.generalized_nfw import GeneralNFWSubhalo, GeneralNFWFieldHalo
 from pyHalo.single_realization import Realization
-from pyHalo.Halos.HaloModels.gaussian import Gaussian
+from pyHalo.Halos.HaloModels.gaussianhalo import GaussianHalo
 from pyHalo.Rendering.correlated_structure import CorrelatedStructure
 from pyHalo.Rendering.MassFunctions.delta_function import DeltaFunction
+from pyHalo.Rendering.MassFunctions.gaussian import Gaussian
 from pyHalo.Cosmology.geometry import Geometry
 from pyHalo.utilities import generate_lens_plane_redshifts, mask_annular
 from pyHalo.Rendering.SpatialDistributions.uniform import Uniform
@@ -32,6 +33,66 @@ def __init__(self, realization):
 
         self._realization = realization
 
+    def add_globular_clusters(self, log10_mgc_mean, log10_mgc_sigma, rendering_radius_arcsec, gc_profile_args,
+                              galaxy_Re=6.0, host_halo_mass=10**13.3, f=3.4e-5, center_x=0, center_y=0):
+        """
+        Add globular clusters at z = z_lens to a realization sampling from a log-normal mass distribution
+        :param log10_mgc_mean: the median log10(mass) of the GC's
+        :param log10_mgc_sigma: the standard deviation of the Gaussian mass function for log10(m)
+        :param rendering_radius_arcsec [arcsec]: sets the area around (center_x, center_y) where the GC's appear; GC's are
+        distributed uniformly in a circule centered at (center_x, center_y) with this radius
+        :param gc_profile_args: the keyword arguments for the GC mass profile, must include "gamma", "gc_size_lightyear",
+        "r_core_fraction", or the logarithmic power law slope, the size of the gc in light years, and the size of the core
+        relative to the size of the GC
+        :param galaxy_Re: galaxy half-light radius
+        :param host_halo_mass: total mass of the host DM halo
+        :param f: the fraction M_gc / M_host from https://arxiv.org/pdf/1602.01105
+        :param coordinate_center_x: center of rendering area in arcsec
+        :param coordinate_center_y: center of rendering area in arcsec
+        :return: an instance of Realization that includes the GC's
+        """
+
+        n = self.number_globular_clusters(log10_mgc_mean, log10_mgc_sigma, rendering_radius_arcsec, galaxy_Re,
+                                          host_halo_mass, f)
+        mfunc = Gaussian(n, log10_mgc_mean, log10_mgc_sigma)
+        m = mfunc.draw()
+        z = self._realization.lens_cosmo.z_lens
+        uniform_spatial_distribution = Uniform(rendering_radius_arcsec, self._realization.geometry)
+        x, y = uniform_spatial_distribution.draw(len(m), z, center_x, center_y)
+        lens_cosmo = self._realization.lens_cosmo
+        GCS = []
+        for (m_gc, x_center_gc, y_center_gc) in zip(m, x, y):
+            unique_tag = np.random.rand()
+            profile = GlobularCluster(m_gc, x_center_gc, y_center_gc, lens_cosmo.z_lens, lens_cosmo,
+                                      gc_profile_args, unique_tag)
+            GCS.append(profile)
+        gcs_realization = Realization.from_halos(GCS, self._realization.lens_cosmo,
+                                                 self._realization.kwargs_halo_model,
+                                                 self._realization.apply_mass_sheet_correction,
+                                                 self._realization.rendering_classes,
+                                                 self._realization._rendering_center_x,
+                                                 self._realization._rendering_center_y,
+                                                 self._realization.geometry)
+        new_realization = self._realization.join(gcs_realization)
+        return new_realization
+
+    def number_globular_clusters(self, log10_mgc_mean, log10_mgc_sigma, rendering_radius_arcsec,
+                              galaxy_Re=10.0, host_halo_mass=10**13.3, f=3.4e-5):
+        """
+
+        :return:
+        """
+        m_gc = f * host_halo_mass
+        area = np.pi * (6 * galaxy_Re) ** 2 # physical kpc^2
+        surface_mass_density = m_gc / area
+        z = self._realization.lens_cosmo.z_lens
+        kpc_per_arcsec = self._realization.lens_cosmo.cosmo.kpc_proper_per_asec(z)
+        # assuming log-normal mass function
+        integral = np.exp(np.log(10 ** log10_mgc_mean) + np.log(10 ** log10_mgc_sigma) ** 2 / 2)
+        mass_in_gc = np.pi * surface_mass_density * (rendering_radius_arcsec * kpc_per_arcsec) ** 2
+        n = int(mass_in_gc / integral)
+        return n
+
     def toSIDM(self, cross_section):
         """
 
@@ -496,14 +557,14 @@ def _get_fluctuation_halos(realization, fluctuation_amplitude, fluctuation_size,
 
     args_fluc=[{'amp': amps[i], 'sigma': sigs[i], 'center_x': xs[i], 'center_y': ys[i]} for i in range(len(amps))]
     masses = np.absolute(amps)
-    fluctuations = [Gaussian(masses[i], xs[i], ys[i], None, realization.lens_cosmo.z_lens,
-                             True, realization.lens_cosmo, args_fluc[i],
-                             truncation_class=None, concentration_class=None,
-                             unique_tag=np.random.rand()) for i in range(len(amps))]
+    fluctuations = [GaussianHalo(masses[i], xs[i], ys[i], None, realization.lens_cosmo.z_lens,
+                                 True, realization.lens_cosmo, args_fluc[i],
+                                 truncation_class=None, concentration_class=None,
+                                 unique_tag=np.random.rand()) for i in range(len(amps))]
 
     return fluctuations
 
-def corr_kappa_with_mask(kappa_map, map_size, r, mu, apply_mask=True, r_min=0.5, r_max=1.5, normalization=False): 
+def corr_kappa_with_mask(kappa_map, map_size, r, mu, apply_mask=True, r_min=0.5, r_max=1.5, normalization=False):
 
     """
     This function computes the two-point correlation function from a convergence map.
@@ -515,7 +576,7 @@ def corr_kappa_with_mask(kappa_map, map_size, r, mu, apply_mask=True, r_min=0.5,
             of the angles between 0 and 180 degrees.
     :param apply_mask: if True, apply the mask on the convergence map.
     :param r_min: inner radius of mask in units of grid coordinates
-    :param r_max: outer radius of mask in units of grid coordinates. If r_max = None, the size of the convergence map's outside 
+    :param r_max: outer radius of mask in units of grid coordinates. If r_max = None, the size of the convergence map's outside
             boundary becomes the mask's outer boundary.
     :param normalization: if True, apply normalization to the correlation function.
     :return: the two-point correlation function on the (mu, r) coordinate grid.
@@ -553,7 +614,7 @@ def corr_kappa_with_mask(kappa_map, map_size, r, mu, apply_mask=True, r_min=0.5,
         mask_interp = RectBivariateSpline(X_, Y_, mask, kx=1, ky=1, s=0)
     else:
         mask = np.ones(XX_.shape)
-        
+
     kappa_interp = RectBivariateSpline(X_, Y_, kappa_map, kx=5, ky=5, s=0)
 
     corr = np.zeros((r.shape[0], mu.shape[0]))
diff --git a/pyHalo/single_realization.py b/pyHalo/single_realization.py
index 3851586..7a96691 100644
--- a/pyHalo/single_realization.py
+++ b/pyHalo/single_realization.py
@@ -9,7 +9,7 @@
 from pyHalo.Halos.HaloModels.PTMass import PTMass
 from pyHalo.Halos.HaloModels.ULDM import ULDMFieldHalo, ULDMSubhalo
 from pyHalo.Halos.HaloModels.NFW_core_trunc import TNFWCFieldHaloSIDM, TNFWCSubhaloSIDM
-from pyHalo.Halos.HaloModels.gaussian import Gaussian
+from pyHalo.Halos.HaloModels.gaussianhalo import GaussianHalo
 import numpy as np
 from copy import deepcopy
 
@@ -677,7 +677,7 @@ def _load_halo_model(mass, x, y, r3d, mdef, z, is_subhalo,
             else:
                 model = ULDMFieldHalo
         elif mdef == 'GAUSSIAN_KAPPA':
-            model = Gaussian
+            model = GaussianHalo
         elif mdef == 'GNFW':
             if is_subhalo:
                 model = GeneralNFWSubhalo
diff --git a/pyHalo/utilities.py b/pyHalo/utilities.py
index cc7c2c9..8909eb7 100644
--- a/pyHalo/utilities.py
+++ b/pyHalo/utilities.py
@@ -203,7 +203,7 @@ def compute_comoving_ray_path(x_coordinate, y_coordinate, lens_model, kwargs_len
         alpha_x_start, alpha_y_start = x_coordinate, y_coordinate
         for zi in all_redshifts_sorted:
 
-            x_start, y_start, alpha_x_start, alpha_y_start = lens_model.lens_model.ray_shooting_partial(x_start, y_start,
+            x_start, y_start, alpha_x_start, alpha_y_start = lens_model.lens_model.ray_shooting_partial_comoving(x_start, y_start,
                                                                                 alpha_x_start, alpha_y_start,
                                                                                 z_start, zi, kwargs_lens)
             d = float(comoving_distance_calc(zi))
diff --git a/tests/test_halos/test_gaussian.py b/tests/test_halos/test_gaussian.py
index a81a91e..14f838e 100644
--- a/tests/test_halos/test_gaussian.py
+++ b/tests/test_halos/test_gaussian.py
@@ -1,6 +1,6 @@
 import numpy.testing as npt
 import numpy as np
-from pyHalo.Halos.HaloModels.gaussian import Gaussian
+from pyHalo.Halos.HaloModels.gaussianhalo import GaussianHalo
 from pyHalo.Halos.lens_cosmo import LensCosmo
 from pyHalo.Cosmology.cosmology import Cosmology
 from astropy.cosmology import FlatLambdaCDM
@@ -22,9 +22,9 @@ def setup_method(self):
         lens_cosmo = LensCosmo(z, 2., cosmo)
         profile_args = {'amp':1,'sigma':1,'center_x':1.0,'center_y':1.0}
         sub_flag = False
-        self.halo = Gaussian(mass, x, y, r3d, z,
-                               sub_flag, lens_cosmo,
-                               profile_args, None, None, unique_tag=np.random.rand())
+        self.halo = GaussianHalo(mass, x, y, r3d, z,
+                                 sub_flag, lens_cosmo,
+                                 profile_args, None, None, unique_tag=np.random.rand())
 
     def test_lenstronomy_params(self):
 
diff --git a/tests/test_halos/test_globular_clusters.py b/tests/test_halos/test_globular_clusters.py
new file mode 100644
index 0000000..28693c0
--- /dev/null
+++ b/tests/test_halos/test_globular_clusters.py
@@ -0,0 +1,52 @@
+import numpy.testing as npt
+from pyHalo.Halos.lens_cosmo import LensCosmo
+from lenstronomy.LensModel.Profiles.splcore import SPLCORE
+from pyHalo.Halos.HaloModels.powerlaw import GlobularCluster
+import pytest
+import numpy as np
+
+
+class TestSPLCORE(object):
+
+    def setup_method(self):
+
+        self.zhalo = 0.5
+        self.zsource = 2.0
+        self.lens_cosmo = LensCosmo(self.zhalo, self.zsource, None)
+        self.splcore = SPLCORE()
+
+    def test_lenstronomy_ID(self):
+
+        mass = 10 ** 5
+        args = {'gamma': 2.5,
+                'r_core_fraction': 0.05,
+                'gc_size_lightyear': 100}
+        profile = GlobularCluster(mass, 0.0, 0.0, self.zhalo, self.lens_cosmo,
+                                  args, 1)
+        lenstronomy_ID = profile.lenstronomy_ID
+        npt.assert_string_equal(lenstronomy_ID[0], 'SPL_CORE')
+
+    def test_lenstronomy_args(self):
+
+        logM = 5.0
+        mass = 10 ** logM
+        args = {'gamma': 2.5,
+                'r_core_fraction': 0.05,
+                'gc_size_lightyear': 100}
+        profile = GlobularCluster(mass, 0.0, 0.0, self.zhalo, self.lens_cosmo,
+                                  args, 1)
+        lenstronomy_args, _ = profile.lenstronomy_params
+        profile_lenstronomy = SPLCORE()
+        profile_args = profile.profile_args
+        gc_size = profile_args['gc_size']
+        rho0 = profile_args['rho0']
+        r_core = profile_args['r_core_arcsec']
+        gamma = profile_args['gamma']
+        mass = profile_lenstronomy.mass_3d(gc_size, rho0, r_core, gamma)
+        sigma_crit_mpc = self.lens_cosmo.get_sigma_crit_lensing(profile.z, self.lens_cosmo.z_source)
+        kpc_per_arcsec = self.lens_cosmo.cosmo.kpc_proper_per_asec(profile.z)
+        sigma_crit_arcsec = sigma_crit_mpc * (0.001 * kpc_per_arcsec) ** 2
+        npt.assert_almost_equal(logM, np.log10(mass * sigma_crit_arcsec))
+
+if __name__ == '__main__':
+    pytest.main()
diff --git a/tests/test_realization_extensions.py b/tests/test_realization_extensions.py
index 1ce471b..b0f4d5e 100644
--- a/tests/test_realization_extensions.py
+++ b/tests/test_realization_extensions.py
@@ -8,11 +8,64 @@
 from pyHalo.Halos.concentration import ConcentrationDiemerJoyce
 from pyHalo.Halos.tidal_truncation import TruncationRoche, TruncationRN
 from pyHalo.Halos.lens_cosmo import LensCosmo
+from pyHalo.PresetModels.cdm import CDM
 from scipy.special import eval_chebyt
 
 
 class TestRealizationExtensions(object):
 
+    def test_globular_clusters(self):
+
+        cdm = CDM(0.5, 2.0, sigma_sub=0.01, LOS_normalization=0.1)
+        n_halos_cdm = len(cdm.halos)
+        ext = RealizationExtensions(cdm)
+        log10_mgc_mean = 4.5
+        log10_mgc_sigma = 0.2
+        rendering_radius_arcsec = 10.0
+        galaxy_Re = 10  # kpc
+        center_x = 1.0
+        center_y = 0.0
+        cdm_with_GCs = ext.add_globular_clusters(log10_mgc_mean,
+                                                   log10_mgc_sigma,
+                                                   rendering_radius_arcsec,
+                                                   galaxy_Re,
+                                                   center_x=center_x,
+                                                   center_y=center_y)
+        n_halos_cdm_plus_gcs = len(cdm_with_GCs.halos)
+        npt.assert_equal(n_halos_cdm_plus_gcs>n_halos_cdm, True)
+
+        cdm = CDM(0.5, 2.0, sigma_sub=0.0, LOS_normalization=0.0)
+        ext = RealizationExtensions(cdm)
+        log10_mgc_mean = 4.5
+        log10_mgc_sigma = 0.2
+        rendering_radius_arcsec = 10.0
+        galaxy_Re = 10 # kpc
+        center_x = 1.0
+        center_y = 0.0
+
+        cdm_with_GCs_1 = ext.add_globular_clusters(log10_mgc_mean,
+                                                 log10_mgc_sigma,
+                                                 rendering_radius_arcsec,
+                                                 galaxy_Re,
+                                                 f=3.4e-4,
+                                                 center_x=center_x,
+                                                   center_y=center_y)
+        scale_re = 3
+        cdm_with_GCs_2 = ext.add_globular_clusters(log10_mgc_mean,
+                                                   log10_mgc_sigma,
+                                                   rendering_radius_arcsec,
+                                                   scale_re*galaxy_Re,
+                                                   )
+        number_1 = len(cdm_with_GCs_1.halos)
+        number_2 = len(cdm_with_GCs_2.halos)
+        npt.assert_almost_equal(number_1/number_2 / scale_re**2, 1.0, 1)
+        masses = []
+        for i in range(0, len(cdm_with_GCs_1.halos)):
+            masses.append(cdm_with_GCs_1.halos[i].mass)
+        log10mean_mass = np.mean(np.log10(masses))
+        npt.assert_almost_equal(log10mean_mass/log10_mgc_mean, 1, 1)
+        npt.assert_almost_equal(np.std(np.log10(masses))/log10_mgc_sigma, 1, 1)
+
     def test_toSIDM(self):
 
         halo_mass = 10 ** 8
@@ -446,7 +499,7 @@ def test_add_pbh(self):
 #         corr = np.ones((r.shape[0], mu.shape[0]))
 #         r, xi_0 = xi_l(0, corr, r, mu)
 #         npt.assert_almost_equal(xi_0_real, xi_0)
-# 
+#
 #     def test_xi_l_to_Pk_l(self):
 #         l = 0
 #         x = np.logspace(-3, 3, num=60, endpoint=False)
@@ -467,6 +520,8 @@ def test_add_pbh(self):
 #         npt.assert_array_almost_equal(As, As_fit)
 #         npt.assert_array_almost_equal(n, n_fit)
 
+t = TestRealizationExtensions()
+t.test_globular_clusters()
 
-if __name__ == '__main__':
-      pytest.main()
+# if __name__ == '__main__':
+#       pytest.main()
diff --git a/tests/test_rendering/test_mass_functions/test_gaussian.py b/tests/test_rendering/test_mass_functions/test_gaussian.py
new file mode 100644
index 0000000..ca876a2
--- /dev/null
+++ b/tests/test_rendering/test_mass_functions/test_gaussian.py
@@ -0,0 +1,36 @@
+import pytest
+import numpy as np
+import numpy.testing as npt
+from pyHalo.Rendering.MassFunctions.gaussian import Gaussian
+
+class TestGaussian(object):
+
+    def setup_method(self):
+        n0 = 5000
+        mean = 5.0
+        sigma = 0.5
+        mfunc = Gaussian(n0, mean, sigma, draw_poisson=False)
+        mfunc_piosson = Gaussian(n0, mean, sigma, draw_poisson=True)
+        self.mean = mean
+        self.sigma = sigma
+        self.n0 = n0
+        self.mass_function_poisson = mfunc_piosson
+        self.mfunc = mfunc
+
+    def test_mass_function(self):
+
+        m = self.mass_function_poisson.draw()
+        log10_m = np.log10(m)
+        mu, sigma = np.mean(log10_m), np.std(log10_m)
+        npt.assert_almost_equal(mu, self.mean, 2)
+        npt.assert_almost_equal(sigma, self.sigma, 2)
+        npt.assert_almost_equal(mu, self.mass_function_poisson.first_moment,2)
+
+    def test_draw_poisson(self):
+
+        m = self.mfunc.draw()
+        m_poisson = self.mass_function_poisson.draw()
+        npt.assert_equal(True, len(m)!=len(m_poisson))
+
+if __name__ == '__main__':
+   pytest.main()