def mass_within_region(gals, x_bound, y_bound):
    """
    Calculate the total mass and number of galaxies within a specified region.

    Parameters
    ----------

    gals: list of ``Galaxy`` class instances.
        Galaxies that we're calculating the mass for.

    x_bound, y_bound: [float, float]
        The minimum and maximum bounds that define the region we're averaging
        inside.

    Returns
    -------

    mass_in_region: float
        The total galaxy mass within the specified region.

    num_gals_in_region: int
        The number of galaxies within the specified region.
    """

    # Initialize our counters.
    mass_in_region = 0.0
    num_gals_in_region = 0
    region_bounds = [x_bound, y_bound]

    for gal in gals:

        gal_pos = [gal.x, gal.y]
        in_region = True

        # We're going to go through each dimension of the galaxy and ask if the position is
        # inside the region.
        for region_bound, dim_pos in zip(region_bounds, gal_pos):

            # Galaxy is outside the region. Flag it and move to the next galaxy.
            if dim_pos < region_bound[0] or dim_pos > region_bound[1]:
                in_region = False
                break

        # Galaxy was in the region, add it.
        if in_region:
            mass_in_region += gal.mass
            num_gals_in_region += 1

    return mass_in_region, num_gals_in_region