Compute the stellar luminosity (Lsun; average per star). PDF
        comes from IMF.  The range of integration only covers the
        input isochrone data (no extrapolation used), but this seems
        like a sub-percent effect if the isochrone goes to 0.15 Msun
        for the old and metal-poor stellar populations of interest.

        Note that the stellar luminosity is very sensitive to the
        post-AGB population.

        Parameters:
        -----------
        steps : Number of steps to sample the isochrone.

        Returns:
        --------
        lum   : The stellar luminosity [Lsun]

def stellar_luminosity(self, steps=10000):
        """
        Compute the stellar luminosity (Lsun; average per star). PDF
        comes from IMF.  The range of integration only covers the
        input isochrone data (no extrapolation used), but this seems
        like a sub-percent effect if the isochrone goes to 0.15 Msun
        for the old and metal-poor stellar populations of interest.

        Note that the stellar luminosity is very sensitive to the
        post-AGB population.

        Parameters:
        -----------
        steps : Number of steps to sample the isochrone.

        Returns:
        --------
        lum   : The stellar luminosity [Lsun]
        """
        mass_min = np.min(self.mass_init)
        mass_max = self.mass_init_upper_bound
        
        d_log_mass = (np.log10(mass_max) - np.log10(mass_min)) / float(steps)
        log_mass = np.linspace(np.log10(mass_min), np.log10(mass_max), steps)
        mass = 10.**log_mass
        
        luminosity_interpolation = scipy.interpolate.interp1d(self.mass_init, self.luminosity,fill_value=0,bounds_error=False)
        luminosity = luminosity_interpolation(mass)

        return np.sum(luminosity * d_log_mass * self.imf.pdf(mass, log_mode=True))