def catsimPopulation(tag, mc_source_id_start=1, n=5000, n_chunk=100, config='simulate_population.yaml'):
    """
    n = Number of satellites to simulation
    n_chunk = Number of satellites in a file chunk
    """

    assert mc_source_id_start >= 1, "Starting mc_source_id must be >= 1" 
    assert n % n_chunk == 0, "Total number of satellites must be divisible by the chunk size"
    nside_pix = 256 # NSIDE = 128 -> 27.5 arcmin, NSIDE = 256 -> 13.7 arcmin 
    
    if not os.path.exists(tag): os.makedirs(tag)

    if isinstance(config,str): config = yaml.load(open(config))
    assert config['survey'] in ['des', 'ps1']

    infile_ebv = config['ebv']
    infile_fracdet = config['fracdet']
    infile_maglim_g = config['maglim_g']
    infile_maglim_r = config['maglim_r']
    infile_density = config['stellar_density']

    range_distance = config.get('range_distance',[5., 500.])
    range_stellar_mass = config.get('range_stellar_mass',[1.e1, 1.e6])
    range_r_physical = config.get('range_r_physical',[1.e-3, 2.0])
    
    m_density = np.load(infile_density)
    nside_density = healpy.npix2nside(len(m_density))
    m_fracdet = read_map(infile_fracdet, nest=False) #.astype(np.float16)
    nside_fracdet = healpy.npix2nside(len(m_fracdet))

    m_maglim_g = read_map(infile_maglim_g, nest=False) #.astype(np.float16)
    m_maglim_r = read_map(infile_maglim_r, nest=False) #.astype(np.float16)

    m_ebv = read_map(infile_ebv, nest=False) #.astype(np.float16)
    
    #m_foreground = healpy.read_map(infile_foreground)

    mask = (m_fracdet > 0.5)

    kwargs = dict(range_distance = range_distance,
                  range_stellar_mass = range_stellar_mass,
                  range_r_physical = range_r_physical)
    print kwargs
    # r_physical is azimuthally-averaged half-light radius, kpc
    simulation_area, lon_population, lat_population, distance_population, stellar_mass_population, r_physical_population = ugali.simulation.population.satellitePopulation(mask, nside_pix, n, **kwargs)
    n_g22_population = np.tile(np.nan, n)
    n_g24_population = np.tile(np.nan, n)
    abs_mag_population = np.tile(np.nan, n)
    surface_brightness_population = np.tile(np.nan, n)
    ellipticity_population = np.tile(np.nan, n)
    position_angle_population = np.tile(np.nan, n)
    age_population = np.tile(np.nan, n)
    metal_z_population = np.tile(np.nan, n)
    mc_source_id_population = np.arange(mc_source_id_start, mc_source_id_start + n)
    #cut_difficulty_population = np.tile(False, n)
    difficulty_population = np.tile(0, n)

    lon_array = []
    lat_array = []
    mag_1_array = []
    mag_2_array = []
    mag_1_error_array = []
    mag_2_error_array = []
    mag_extinction_1_array = []
    mag_extinction_2_array = []
    mc_source_id_array = []
    for ii, mc_source_id in enumerate(mc_source_id_population):
        print '  Simulating satellite (%i/%i) ... MC_SOURCE_ID = %i'%(ii + 1, n, mc_source_id)
        print '  distance=%.2e, stellar_mass=%.2e, rhalf=%.2e'%(distance_population[ii],stellar_mass_population[ii],r_physical_population[ii])
        lon, lat, mag_1, mag_2, mag_1_error, mag_2_error, mag_extinction_1, mag_extinction_2, n_g22, n_g24, abs_mag, surface_brightness, ellipticity, position_angle, age, metal_z, flag_too_extended = catsimSatellite(config,
                                                                                                                                                                             lon_population[ii], 
                                                                                                                                                                             lat_population[ii], 
                                                                                                                                                                             distance_population[ii], 
                                                                                                                                                                             stellar_mass_population[ii], 
                                                                                                                                                                             r_physical_population[ii],
                                                                                                                                                                             m_maglim_g,
                                                                                                                                                                             m_maglim_r,
                                                                                                                                                                             m_ebv)
        print '  ', len(lon)
        
        n_g22_population[ii] = n_g22
        n_g24_population[ii] = n_g24
        abs_mag_population[ii] = abs_mag
        surface_brightness_population[ii] = surface_brightness
        ellipticity_population[ii] = ellipticity
        position_angle_population[ii] = position_angle
        age_population[ii] = age
        metal_z_population[ii] = metal_z

        #print "Difficulty masking..."  

        # These objects are too extended and are not simulated
        if (flag_too_extended):
            difficulty_population[ii] |= 0b0001

        # We assume that these objects would be easily detected and
        # remove them to reduce data volume
        if (surface_brightness_population[ii]<25.)&(n_g22_population[ii]>1e2):
            difficulty_population[ii] |= 0b0010
        if (surface_brightness_population[ii]<28.)&(n_g22_population[ii]>1e4):
            difficulty_population[ii] |= 0b0100
        if (surface_brightness_population[ii]<30.)&(n_g22_population[ii]>1e5):
            difficulty_population[ii] |= 0b1000
        
        #cut_easy = (surface_brightness_population[ii]<25.)&(n_g22_population[ii]>1.e2) \
        #           | ((surface_brightness_population[ii] < 30.) & (n_g24_population[ii] > 1.e4)) \
        #           | ((surface_brightness_population[ii] < 31.) & (n_g24_population[ii] > 1.e5))
        #cut_hard = (surface_brightness_population[ii] > 35.) | (n_g24_population[ii] < 1.)
        #cut_difficulty_population[ii] = ~cut_easy & ~cut_hard
        #if cut_easy:
        #    difficulty_population[ii] += 1 # TOO EASY
        #if cut_hard:
        #    difficulty_population[ii] += 2 # TOO HARD
        #if flag_too_extended:
        #    difficulty_population[ii] += 3 # TOO EXTENDED

        if difficulty_population[ii] == 0:
            lon_array.append(lon)
            lat_array.append(lat)
            mag_1_array.append(mag_1)
            mag_2_array.append(mag_2)
            mag_1_error_array.append(mag_1_error)
            mag_2_error_array.append(mag_2_error)
            mag_extinction_1_array.append(mag_extinction_1)
            mag_extinction_2_array.append(mag_extinction_2)
            mc_source_id_array.append(np.tile(mc_source_id, len(lon)))

    # Concatenate all the arrays

    print "Concatenating arrays..."
    lon_array = np.concatenate(lon_array)
    lat_array = np.concatenate(lat_array)
    mag_1_array = np.concatenate(mag_1_array)
    mag_2_array = np.concatenate(mag_2_array)
    mag_1_error_array = np.concatenate(mag_1_error_array)
    mag_2_error_array = np.concatenate(mag_2_error_array)
    mag_extinction_1_array = np.concatenate(mag_extinction_1_array)
    mag_extinction_2_array = np.concatenate(mag_extinction_2_array)
    mc_source_id_array = np.concatenate(mc_source_id_array)

    # Now do the masking all at once

    print "Fracdet masking..."
    pix_array = ugali.utils.healpix.angToPix(nside_fracdet, lon_array, lat_array)
    cut_fracdet = (np.random.uniform(size=len(lon_array)) < m_fracdet[pix_array])

    lon_array = lon_array[cut_fracdet]
    lat_array = lat_array[cut_fracdet]
    mag_1_array = mag_1_array[cut_fracdet]
    mag_2_array = mag_2_array[cut_fracdet]
    mag_1_error_array = mag_1_error_array[cut_fracdet]
    mag_2_error_array = mag_2_error_array[cut_fracdet]
    mag_extinction_1_array = mag_extinction_1_array[cut_fracdet]
    mag_extinction_2_array = mag_extinction_2_array[cut_fracdet]
    mc_source_id_array = mc_source_id_array[cut_fracdet]

    # Cut out the entries that are easily detectable
    """
    lon_population = lon_population[cut_difficulty_population]
    lat_population = lat_population[cut_difficulty_population]
    distance_population = distance_population[cut_difficulty_population]
    stellar_mass_population = stellar_mass_population[cut_difficulty_population]
    r_physical_population = r_physical_population[cut_difficulty_population]
    n_g24_population = n_g24_population[cut_difficulty_population]
    abs_mag_population = abs_mag_population[cut_difficulty_population]
    surface_brightness_population = surface_brightness_population[cut_difficulty_population]
    ellipticity_population = ellipticity_population[cut_difficulty_population]
    position_angle_population = position_angle_population[cut_difficulty_population]
    age_population = age_population[cut_difficulty_population]
    metal_z_population = metal_z_population[cut_difficulty_population]
    mc_source_id_population = mc_source_id_population[cut_difficulty_population]
    """
    
    # Create bonus columns
    
    print "Creating bonus columns..."
    distance_modulus_population = ugali.utils.projector.distanceToDistanceModulus(distance_population)
    hpix_32_population = ugali.utils.healpix.angToPix(32, lon_population, lat_population) # Make sure this matches the dataset

    # Local stellar density
    pixarea = healpy.nside2pixarea(nside_density, degrees=True) * 60.**2 # arcmin^2
    density_population = m_density[ugali.utils.healpix.angToPix(nside_density, lon_population, lat_population)] / pixarea # arcmin^-2

    # Average fracdet within the azimuthally averaged half-light radius
    #m_fracdet_zero = np.where(m_fracdet >= 0., m_fracdet, 0.)
    #m_fracdet_zero = m_fracdet
    r_half = np.degrees(np.arctan2(r_physical_population, distance_population)) # Azimuthally averaged half-light radius in degrees
    fracdet_half_population = meanFracdet(m_fracdet, lon_population, lat_population, r_half)
    fracdet_core_population = meanFracdet(m_fracdet, lon_population, lat_population, 0.1)
    fracdet_wide_population = meanFracdet(m_fracdet, lon_population, lat_population, 0.5)

    # Magnitude limits
    nside_maglim = healpy.npix2nside(len(m_maglim_g))
    pix_population = ugali.utils.healpix.angToPix(nside_maglim, lon_population, lat_population)
    maglim_g_population = m_maglim_g[pix_population]
    maglim_r_population = m_maglim_r[pix_population]
    
    # E(B-V)
    nside_ebv = healpy.npix2nside(len(m_ebv))
    pix_population = ugali.utils.healpix.angToPix(nside_ebv, lon_population, lat_population)
    ebv_population = m_ebv[pix_population]

    # Survey
    survey_population = np.tile(config['survey'], len(lon_population))

    # Number of surviving catalog stars
    n_catalog_population = np.histogram(mc_source_id_array, bins=np.arange(mc_source_id_population[0] - 0.5, mc_source_id_population[-1] + 0.51))[0]

    # Faked-up coadd_object_ids
    coadd_object_id_array = []
    for mc_source_id in mc_source_id_population:
        coadd_object_id_array.append((1000000 * mc_source_id) + 1 + np.arange(np.sum(mc_source_id == mc_source_id_array)))
    coadd_object_id_array = -1 * np.concatenate(coadd_object_id_array) # Assign negative numbers to distinguish from real objects

    # Catalog output file

    # for ii in range(0, len(d.formats)): print '\'%s\': [ , \'%s\'],'%(d.names[ii], d.formats[ii])
    
    # See: 
    # https://github.com/sidneymau/simple/blob/master/search_algorithm.py 
    # https://github.com/sidneymau/simple/blob/master/config.yaml
    # /home/s1/kadrlica/projects/y3a2/dsphs/v2/skim/ , e.g., /home/s1/kadrlica/projects/y3a2/dsphs/v2/skim/y3a2_ngmix_cm_11755.fits

    #default_array = np.tile(np.nan, len(mc_source_id_array)) # To recognize that those values are synthetic filler
    default_array = np.tile(-9999., len(mc_source_id_array))

    """
    # Column name, data, fits format
    # Y3A2 pre-Gold
    key_map = {'CM_MAG_ERR_G': [mag_1_error_array, 'D'],
               'CM_MAG_ERR_R': [mag_2_error_array, 'D'],
               'CM_MAG_G': [mag_1_array, 'D'],
               'CM_MAG_R': [mag_2_array, 'D'],
               'CM_T': [default_array, 'D'],
               'CM_T_ERR': [default_array, 'D'],
               'COADD_OBJECT_ID': [coadd_object_id_array, 'K'],
               'DEC': [lat_array, 'D'],
               'FLAGS': [default_array, 'K'],
               'PSF_MAG_ERR_G': [mag_1_error_array, 'D'],
               'PSF_MAG_ERR_R': [mag_2_error_array, 'D'],
               'PSF_MAG_G': [mag_1_array, 'D'],
               'PSF_MAG_R': [mag_2_array, 'D'],
               'RA': [lon_array, 'D'],
               'SEXTRACTOR_FLAGS_G': [np.tile(0, len(mc_source_id_array)), 'I'],
               'SEXTRACTOR_FLAGS_R': [np.tile(0, len(mc_source_id_array)), 'I'],
               'WAVG_MAG_PSF_G': [mag_1_array, 'E'],
               'WAVG_MAG_PSF_R': [mag_2_array, 'E'],
               'WAVG_MAGERR_PSF_G': [mag_1_error_array, 'E'],
               'WAVG_MAGERR_PSF_R': [mag_2_error_array, 'E'],
               'WAVG_SPREAD_MODEL_I': [default_array, 'E'],
               'WAVG_SPREADERR_MODEL_I': [default_array, 'E'],
               'EXT_SFD98_G': [default_array, 'E'],
               'EXT_SFD98_R': [default_array, 'E'],
               'CM_MAG_SFD_G': [mag_1_array, 'D'],
               'CM_MAG_SFD_R': [mag_2_array, 'D'],
               'FLAG_FOOTPRINT': [np.tile(1, len(mc_source_id_array)), 'J'],
               'FLAG_FOREGROUND': [np.tile(0, len(mc_source_id_array)), 'J'],
               'EXTENDED_CLASS_MASH': [np.tile(0, len(mc_source_id_array)), 'K'],
               'PSF_MAG_SFD_G': [mag_1_array, 'D'],
               'PSF_MAG_SFD_R': [mag_2_array, 'D'],
               'WAVG_MAG_PSF_SFD_G': [mag_1_array, 'E'],
               'WAVG_MAG_PSF_SFD_R': [mag_2_array, 'E']}
    """
    
    if config['survey'] == 'des':
        # Y3 Gold v2.0
        key_map = odict([
                ('COADD_OBJECT_ID', [coadd_object_id_array, 'K']),
                ('RA', [lon_array, 'D']),
                ('DEC', [lat_array, 'D']),
                ('SOF_PSF_MAG_CORRECTED_G', [mag_1_array, 'D']),
                ('SOF_PSF_MAG_CORRECTED_R', [mag_2_array, 'D']),
                ('SOF_PSF_MAG_ERR_G', [mag_1_error_array, 'D']),
                ('SOF_PSF_MAG_ERR_R', [mag_2_error_array, 'D']),
                ('A_SED_SFD98_G', [mag_extinction_1_array, 'E']),
                ('A_SED_SFD98_R', [mag_extinction_2_array, 'E']),
                ('WAVG_MAG_PSF_G', [mag_1_array+mag_extinction_1_array, 'E']),
                ('WAVG_MAG_PSF_R', [mag_2_array+mag_extinction_2_array, 'E']),
                ('WAVG_MAGERR_PSF_G', [mag_1_error_array, 'E']),
                ('WAVG_MAGERR_PSF_R', [mag_2_error_array, 'E']),
                ('WAVG_SPREAD_MODEL_I', [default_array, 'E']),
                ('WAVG_SPREADERR_MODEL_I', [default_array, 'E']),
                ('SOF_CM_T', [default_array, 'D']),
                ('SOF_CM_T_ERR', [default_array, 'D']),
                ('FLAGS_GOLD', [np.tile(0, len(mc_source_id_array)), 'J']),
                ('EXTENDED_CLASS_MASH_SOF', [np.tile(0, len(mc_source_id_array)), 'I']),
                ])
    elif config['survey'] == 'ps1':
        # PS1
        key_map = odict([
                ('OBJID', [coadd_object_id_array, 'K']),
                ('RA', [lon_array, 'D']),
                ('DEC', [lat_array, 'D']),
                #('UNIQUEPSPSOBID', [coadd_object_id_array, 'K']),
                #('OBJINFOFLAG', [default_array, 'E']),
                #('QUALITYFLAG', [np.tile(16, len(mc_source_id_array)), 'I']),
                #('NSTACKDETECTIONS', [np.tile(99, len(mc_source_id_array)), 'I']),
                #('NDETECTIONS', [np.tile(99, len(mc_source_id_array)), 'I']),
                #('NG', [default_array, 'E']),
                #('NR', [default_array, 'E']),
                #('NI', [default_array, 'E']),
                ('GFPSFMAG', [mag_1_array+mag_extinction_1_array, 'E']),
                ('RFPSFMAG', [mag_2_array+mag_extinction_2_array, 'E']),
                #('IFPSFMAG', [np.tile(0., len(mc_source_id_array)), 'E'], # Too pass star selection
                ('GFPSFMAGERR', [mag_1_error_array, 'E']),
                ('RFPSFMAGERR', [mag_2_error_array, 'E']),
                #('IFPSFMAGERR', [default_array, 'E']),
                #('GFKRONMAG', [mag_1_array, 'E']),
                #('RFKRONMAG', [mag_2_array, 'E']),
                #('IFKRONMAG', [np.tile(0., len(mc_source_id_array)), 'E'], # Too pass star selection
                #('GFKRONMAGERR', [mag_1_error_array, 'E']),
                #('RFKRONMAGERR', [mag_2_error_array, 'E']),
                #('IFKRONMAGERR', [default_array, 'E']),
                #('GFLAGS', [np.tile(0, len(mc_source_id_array)), 'I']),
                #('RFLAGS', [np.tile(0, len(mc_source_id_array)), 'I']),
                #('IFLAGS', [np.tile(0, len(mc_source_id_array)), 'I']),
                #('GINFOFLAG', [np.tile(0, len(mc_source_id_array)), 'I']),
                #('RINFOFLAG', [np.tile(0, len(mc_source_id_array)), 'I']),
                #('IINFOFLAG', [np.tile(0, len(mc_source_id_array)), 'I']),
                #('GINFOFLAG2', [np.tile(0, len(mc_source_id_array)), 'I']),
                #('RINFOFLAG2', [np.tile(0, len(mc_source_id_array)), 'I']),
                #('IINFOFLAG2', [np.tile(0, len(mc_source_id_array)), 'I']),
                #('GINFOFLAG3', [np.tile(0, len(mc_source_id_array)), 'I']),
                #('RINFOFLAG3', [np.tile(0, len(mc_source_id_array)), 'I']),
                #('IINFOFLAG3', [np.tile(0, len(mc_source_id_array)), 'I']),
                #('PRIMARYDETECTION', [default_array, 'E']),
                #('BESTDETECTION', [default_array, 'E']),
                #('EBV', [default_array, 'E']),
                #('EXTSFD_G', [mag_extinction_1_array 'E']),
                #('EXTSFD_R', [mag_extinction_2_array, 'E']),
                #('EXTSFD_I', [default_array, 'E']),
                ('GFPSFMAG_SFD', [mag_1_array, 'E']),
                ('RFPSFMAG_SFD', [mag_2_array, 'E']),
                ('EXTENDED_CLASS', [np.tile(0, len(mc_source_id_array)), 'I']),
                ])
    key_map['MC_SOURCE_ID'] = [mc_source_id_array, 'K']

    print "Writing catalog files..."
    columns = []
    for key in key_map:
        columns.append(pyfits.Column(name=key, format=key_map[key][1], array=key_map[key][0]))
    tbhdu = pyfits.BinTableHDU.from_columns(columns)
    tbhdu.header.set('AREA', simulation_area, 'Simulation area (deg^2)')

    for mc_source_id_chunk in np.split(np.arange(mc_source_id_start, mc_source_id_start + n), n / n_chunk):
        print '  writing MC_SOURCE_ID values from %i to %i'%(mc_source_id_chunk[0], mc_source_id_chunk[-1])
        cut_chunk = np.in1d(mc_source_id_array, mc_source_id_chunk)
        outfile = '%s/sim_catalog_%s_mc_source_id_%07i-%07i.fits'%(tag, tag, mc_source_id_chunk[0], mc_source_id_chunk[-1])
        header = copy.deepcopy(tbhdu.header)
        header.set('IDMIN',mc_source_id_chunk[0], 'Minimum MC_SOURCE_ID')
        header.set('IDMAX',mc_source_id_chunk[-1], 'Maximum MC_SOURCE_ID')
        pyfits.writeto(outfile, tbhdu.data[cut_chunk], header, clobber=True)

    # Population metadata output file
    
    print "Writing population metadata file..."
    tbhdu = pyfits.BinTableHDU.from_columns([
        pyfits.Column(name='RA', format='E', array=lon_population, unit='deg'),
        pyfits.Column(name='DEC', format='E', array=lat_population, unit='deg'),
        pyfits.Column(name='DISTANCE', format='E', array=distance_population, unit='kpc'),
        pyfits.Column(name='DISTANCE_MODULUS', format='E', array=distance_modulus_population, unit='kpc'),
        pyfits.Column(name='STELLAR_MASS', format='E', array=stellar_mass_population, unit='m_solar'),
        pyfits.Column(name='R_PHYSICAL', format='E', array=r_physical_population, unit='kpc'),
        pyfits.Column(name='N_G22', format='J', array=n_g22_population, unit=''),
        pyfits.Column(name='N_G24', format='J', array=n_g24_population, unit=''),
        pyfits.Column(name='N_CATALOG', format='J', array=n_catalog_population, unit=''),
        pyfits.Column(name='DIFFICULTY', format='J', array=difficulty_population, unit=''),
        pyfits.Column(name='ABS_MAG', format='E', array=abs_mag_population, unit='mag'),
        pyfits.Column(name='SURFACE_BRIGHTNESS', format='E', array=surface_brightness_population, unit='mag arcsec^-2'),
        pyfits.Column(name='ELLIPTICITY', format='E', array=ellipticity_population, unit=''),
        pyfits.Column(name='POSITION_ANGLE', format='E', array=position_angle_population, unit='deg'),
        pyfits.Column(name='AGE', format='E', array=age_population, unit='deg'),
        pyfits.Column(name='METAL_Z', format='E', array=metal_z_population, unit=''),
        pyfits.Column(name='MC_SOURCE_ID', format='K', array=mc_source_id_population, unit=''),
        pyfits.Column(name='HPIX_32', format='E', array=hpix_32_population, unit=''),
        pyfits.Column(name='DENSITY', format='E', array=density_population, unit='arcmin^-2'),
        pyfits.Column(name='FRACDET_HALF', format='E', array=fracdet_half_population, unit=''),
        pyfits.Column(name='FRACDET_CORE', format='E', array=fracdet_core_population, unit=''),
        pyfits.Column(name='FRACDET_WIDE', format='E', array=fracdet_wide_population, unit=''),
        pyfits.Column(name='MAGLIM_G', format='E', array=maglim_g_population, unit='mag'),
        pyfits.Column(name='MAGLIM_R', format='E', array=maglim_r_population, unit='mag'),
        pyfits.Column(name='EBV', format='E', array=ebv_population, unit='mag'),
        pyfits.Column(name='SURVEY', format='A12', array=survey_population, unit=''),
    ])
    tbhdu.header.set('AREA', simulation_area, 'Simulation area (deg^2)')
    tbhdu.writeto('%s/sim_population_%s_mc_source_id_%07i-%07i.fits'%(tag, tag, mc_source_id_start, mc_source_id_start + n - 1), clobber=True)

    # 5284.2452461023322

    # Mask output file

    print "Writing population mask file..."
    outfile_mask = '%s/sim_mask_%s_cel_nside_%i.fits'%(tag, tag, healpy.npix2nside(len(mask)))
    if not os.path.exists(outfile_mask):
        healpy.write_map(outfile_mask, mask.astype(int), nest=True, coord='C', overwrite=True)
        os.system('gzip -f %s'%(outfile_mask))