docs: def satellitePopulation(mask, nside_pix, n,
                          range_distance=[5., 500.],
                          range_stellar_mass=[1.e1, 1.e6],
                          range_r_physical=[1.e-3, 2.],
                          plot=False):
      """
      Create a population of n randomly placed satellites within a survey mask.
      Satellites are distributed uniformly in log(distance) (kpc), uniformly in log(stellar_mass) (M_sol), and uniformly in
      physical half-light radius log(r_physical) (kpc). The ranges can be set by the user.
  
      Returns the simulated area (deg^2) as well as the
      lon (deg), lat (deg), distance modulus, stellar mass (M_sol), and half-light radius (deg) for each satellite
  
      Parameters:
      -----------
      mask : the survey mask of available area
      nside_pix : coarse resolution npix for avoiding small gaps in survey
      n : number of satellites to simulate
  
      Returns:
      --------
      area, lon, lat, distance, stellar_mass, r_physical    
      """
      
      distance = 10**np.random.uniform(np.log10(range_distance[0]),
                                       np.log10(range_distance[1]),
                                       n)
  
      stellar_mass = 10**np.random.uniform(np.log10(range_stellar_mass[0]), 
                                           np.log10(range_stellar_mass[1]), 
                                           n)
      
      # Physical half-light radius (kpc)
      r_physical = 10**np.random.uniform(np.log10(range_r_physical[0]), 
                                         np.log10(range_r_physical[1]), 
                                         n)
  
      # Call positions last because while loop has a variable number of calls to np.random (thus not preserving seed information)
      lon, lat, simulation_area = ugali.utils.skymap.randomPositions(mask, nside_pix, n=n)
  
      #half_light_radius = np.degrees(np.arcsin(half_light_radius_physical \
      #                                         / ugali.utils.projector.distanceModulusToDistance(distance_modulus)))
  
      # One choice of theory prior
      #half_light_radius_physical = ugali.analysis.kernel.halfLightRadius(stellar_mass) # kpc
      #half_light_radius = np.degrees(np.arcsin(half_light_radius_physical \
      #                                               / ugali.utils.projector.distanceModulusToDistance(distance_modulus)))
  
      if plot:
          pylab.figure()
          #pylab.scatter(lon, lat, c=distance_modulus, s=500 * half_light_radius)
          #pylab.colorbar()
          pylab.scatter(lon, lat, edgecolors='none')
          xmin, xmax = pylab.xlim() # Reverse azimuthal axis
          pylab.xlim([xmax, xmin])
          pylab.title('Random Positions in Survey Footprint')
          pylab.xlabel('Longitude (deg)')
          pylab.ylabel('Latitude (deg)')
  
          pylab.figure()
          pylab.scatter(stellar_mass, ugali.utils.projector.distanceModulusToDistance(distance_modulus),
                        c=(60. * half_light_radius), s=500 * half_light_radius, edgecolors='none')
          pylab.xscale('log')
          pylab.yscale('log')
          pylab.xlim([0.5 * range_stellar_mass[0], 2. * range_stellar_mass[1]])
          pylab.colorbar()
          pylab.title('Half-light Radius (arcmin)')
          pylab.xlabel('Stellar Mass (arcmin)')
          pylab.ylabel('Distance (kpc)')
  
      return simulation_area, lon, lat, distance, stellar_mass, r_physical