docs: 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))