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