Correction to crowding metric

[gpastorelli-astro]

The implementation of the crowding error calculation in the function _comp_crowd_error requires differential luminosity functions (LFs) as input.

However, the current LFs files are stored as cumulative LFs and they are directly used as input to that function.

This applies to the case of LF files calculated from the TRILEGAL simulations (rubin_sim_data/maps/TriMaps/TRIstarDensity*.npz), and also to the Galfast LFs (rubin_sim_data/maps/StarMaps/starDensity*.npz).

As a quick fix, differential LFs can be computed from the cumulative ones, once the latter are loaded in the function _comp_crowd_error. This should work with current LF files.
See below for a suggested fix.

Please, let me know if further details are needed.
@rhiannonlynne @knutago

Suggested fix:

def _comp_crowd_error(mag_vector, lum_func, seeing, single_mag=None):
    """
    Compute the photometric crowding error given the luminosity
    function and best seeing.

    Equation from Olsen, Blum, & Rigaut 2003, AJ, 126, 452

    Parameters
    ----------
    mag_vector : `np.array` (N,)
        Stellar magnitudes.
    lum_func : `np.array` (N,)
        Stellar luminosity function.
    seeing : `float`
        The best seeing conditions.
        Assuming forced-photometry can use the best seeing conditions
        to help with confusion errors.
    single_mag : `float` or None
        If single_mag is None, the crowding error is calculated
        for each mag in mag_vector. If single_mag is a float,
        the crowding error is interpolated to that single value.

    Returns
    -------
    mag_uncert : `np.array` (N,)
        Magnitude uncertainties.
    """
    lum_area_arcsec = 3600.0**2

    # ==> Suggested fix starts
    lum_func = np.diff(lum_func) # NEW =>  cumulative LF to differential LF
    mag_vector = mag_vector[1:] # NEW => first point of differential LF cannot be recovered
    #<== Suggested fix ends

    lum_vector = 10 ** (-0.4 * mag_vector)
    coeff = np.sqrt(np.pi / lum_area_arcsec) * seeing / 2.0
    my_int = (np.add.accumulate((lum_vector**2 * lum_func)[::-1]))[::-1]
    temp = np.sqrt(my_int) / lum_vector
    if single_mag is not None:
        interp = interp1d(mag_vector, temp)
        temp = interp(single_mag)
    crowd_error = coeff * temp
    return crowd_error