docs: def catsimSatellite(config, lon_centroid, lat_centroid, distance, stellar_mass, r_physical, 
                      m_maglim_1, m_maglim_2, m_ebv,
                      plot=False, title='test'):
      """
      Simulate a single satellite. This is currently only valid for band_1 = g and band_2 = r.
      r_physical is azimuthally averaged half-light radius, kpc
      """
  
      # Probably don't want to parse every time
      completeness = getCompleteness(config)
      log_photo_error = getPhotoError(config)
  
      s = ugali.analysis.source.Source()
  
      # Following McConnachie 2012, ellipticity = 1 - (b/a) , where a is semi-major axis and b is semi-minor axis
      
      r_h = np.degrees(np.arcsin(r_physical / distance)) # Azimuthally averaged half-light radius
      #ellipticity = 0.3 # Semi-arbitrary default for testing purposes
      # See http://iopscience.iop.org/article/10.3847/1538-4357/833/2/167/pdf
      # Based loosely on https://arxiv.org/abs/0805.2945
      ellipticity = np.random.uniform(0.1, 0.8)
      position_angle = np.random.uniform(0., 180.) # Random position angle (deg)
      a_h = r_h / np.sqrt(1. - ellipticity) # semi-major axis (deg)
      
          
      # Elliptical kernels take the "extension" as the semi-major axis
      ker = ugali.analysis.kernel.EllipticalPlummer(lon=lon_centroid, lat=lat_centroid, ellipticity=ellipticity, position_angle=position_angle)
  
      flag_too_extended = False
      if a_h >= 1.0:
          print 'Too extended: a_h = %.2f'%(a_h)
          a_h = 1.0
          flag_too_extended = True
      ker.setp('extension', value=a_h, bounds=[0.0,1.0])
      s.set_kernel(ker)
      
      age = np.random.choice([10., 12.0, 13.5])
      metal_z = np.random.choice([0.0001, 0.0002])
      distance_modulus = ugali.utils.projector.distanceToDistanceModulus(distance)
      iso = isochrone_factory('Bressan2012', survey=config['survey'], age=age, z=metal_z, distance_modulus=distance_modulus)
      s.set_isochrone(iso)
      # Simulate takes stellar mass as an argument, NOT richness
      mag_1, mag_2 = s.isochrone.simulate(stellar_mass) 
  
      lon, lat = s.kernel.sample_lonlat(len(mag_2))
  
      nside = healpy.npix2nside(len(m_maglim_1)) # Assuming that the two maglim maps have same resolution
      pix = ugali.utils.healpix.angToPix(nside, lon, lat)
      maglim_1 = m_maglim_1[pix]
      maglim_2 = m_maglim_2[pix]
      if config['survey'] == 'des':
          # DES Y3 Gold fiducial
          mag_extinction_1 = 3.186 * m_ebv[pix]
          mag_extinction_2 = 2.140 * m_ebv[pix]
      elif config['survey'] == 'ps1':
          # From Table 6 in Schlafly 2011 with Rv = 3.1
          # http://iopscience.iop.org/article/10.1088/0004-637X/737/2/103/pdf
          mag_extinction_1 = 3.172 * m_ebv[pix]
          mag_extinction_2 = 2.271 * m_ebv[pix]
      
      # Photometric uncertainties are larger in the presence of interstellar dust reddening
      mag_1_error = 0.01 + 10**(log_photo_error((mag_1 + mag_extinction_1) - maglim_1))
      mag_2_error = 0.01 + 10**(log_photo_error((mag_2 + mag_extinction_2) - maglim_2))
  
      # It would be better to convert to a flux uncertainty and then transform back to a magnitude
      #mag_1_meas = mag_1 + np.random.normal(scale=mag_1_error)
      #mag_2_meas = mag_2 + np.random.normal(scale=mag_2_error)
      flux_1_meas = magToFlux(mag_1) + np.random.normal(scale=getFluxError(mag_1, mag_1_error))
      mag_1_meas = np.where(flux_1_meas > 0., fluxToMag(flux_1_meas), 99.)
      flux_2_meas = magToFlux(mag_2) + np.random.normal(scale=getFluxError(mag_2, mag_2_error))
      mag_2_meas = np.where(flux_2_meas > 0., fluxToMag(flux_2_meas), 99.)
  
      # In the HSC SXDS ultra-deep field:
      # mean maglim_r_sof_gold_2.0 = 23.46
      # median maglim_r_sof_gold_2.0 = 23.47
      # m = healpy.read_map('/Users/keithbechtol/Documents/DES/projects/mw_substructure/des/y3a1/data/maps/y3a2_gold_1.0_cmv02-001_v1_nside4096_nest_r_depth.fits.gz')
      # np.mean(m[ugali.utils.healpix.angToDisc(4096, 34.55, -4.83, 0.75)])
      # np.median(m[ugali.utils.healpix.angToDisc(4096, 34.55, -4.83, 0.75)])
  
      # Includes penalty for interstellar extinction and also include variations in depth
      if config['survey'] == 'des':
          cut_detect = (np.random.uniform(size=len(mag_2)) < completeness(mag_2 + mag_extinction_2 + (23.46 - np.clip(maglim_2, 20., 26.))))
      elif config['survey'] == 'ps1':
          cut_detect = (np.random.uniform(size=len(mag_2)) < completeness(mag_2 + mag_extinction_2))
  
      n_g22 = np.sum(cut_detect & (mag_1 < 22.))
      n_g24 = np.sum(cut_detect & (mag_1 < 24.))
      print '  n_sim = %i, n_detect = %i, n_g24 = %i, n_g22 = %i'%(len(mag_1),np.sum(cut_detect),n_g24,n_g22)
      
      richness = stellar_mass / s.isochrone.stellarMass()
      #abs_mag = s.isochrone.absolute_magnitude()
      #abs_mag_martin = s.isochrone.absolute_magnitude_martin(richness=richness, n_trials=10)[0] # 100 trials seems to be sufficient for rough estimate
      #print 'abs_mag_martin = %.2f mag'%(abs_mag_martin)
  
      # The more clever thing to do would be to sum up the actual simulated stars
      if config['survey'] == 'des':
          v = mag_1 - 0.487*(mag_1 - mag_2) - 0.0249 # See https://github.com/DarkEnergySurvey/ugali/blob/master/ugali/isochrone/model.py
      elif config['survey'] == 'ps1':
          # https://arxiv.org/pdf/1706.06147.pdf
          # V - g = C_0 + C_1 * (g - r)
          C_0 = -0.017
          C_1 = -0.508
          v = mag_1 + C_0 + C_1 * (mag_1 - mag_2)
      flux = np.sum(10**(-v/2.5))
      abs_mag = -2.5*np.log10(flux) - distance_modulus
  
      #print abs_mag, abs_mag_martin
  
      #distance = ugali.utils.projector.distanceModulusToDistance(distance_modulus)
      #r_h = extension * np.sqrt(1. - ellipticity) # Azimuthally averaged half-light radius
      r_physical = distance * np.tan(np.radians(r_h)) # Azimuthally averaged half-light radius, kpc
      #print 'distance = %.3f kpc'%(distance)
      #print 'r_physical = %.3f kpc'%(r_physical)
      surface_brightness = ugali.analysis.results.surfaceBrightness(abs_mag, r_physical, distance) # Average within azimuthally averaged half-light radius
      #print 'surface_brightness = %.3f mag arcsec^-2'%(surface_brightness)
      
      if plot:
          import pylab
          pylab.ion()
  
          n_sigma_p = np.sum(cut_detect & (mag_1 < 23.))
  
          pylab.figure(figsize=(6., 6.))
          pylab.scatter(mag_1_meas[cut_detect] - mag_2_meas[cut_detect], mag_1_meas[cut_detect], edgecolor='none', c='black', s=5)
          pylab.xlim(-0.5, 1.)
          pylab.ylim(26., 16.)
          pylab.xlabel('g - r')
          pylab.ylabel('g')
          pylab.title('Number of stars with g < 23: %i'%(n_sigma_p))
          pylab.savefig('y3_sat_sim_cmd_%s.png'%(title), dpi=150.)
          
          print 'n_Sigma_p = %i'%(n_sigma_p)
          raw_input('WAIT')
          
      #if flag_too_extended:
      #    # This is a kludge to remove these satellites. fragile!!
      #    n_g24 = 1.e6
  
      return lon[cut_detect], lat[cut_detect], mag_1_meas[cut_detect], mag_2_meas[cut_detect], mag_1_error[cut_detect], mag_2_error[cut_detect], mag_extinction_1[cut_detect], mag_extinction_2[cut_detect], n_g22, n_g24, abs_mag, surface_brightness, ellipticity, position_angle, age, metal_z, flag_too_extended