code: def stellar_luminosity(self, steps=10000):
          
          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))
docs: 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]