code: def simulate(self, stellar_mass, distance_modulus=None, **kwargs):
          
          if distance_modulus is None: distance_modulus = self.distance_modulus
          # Total number of stars in system
          n = int(round(stellar_mass / self.stellar_mass()))
          f_1 = scipy.interpolate.interp1d(self.mass_init, self.mag_1)
          f_2 = scipy.interpolate.interp1d(self.mass_init, self.mag_2)
          mass_init_sample = self.imf.sample(n, np.min(self.mass_init), np.max(self.mass_init), **kwargs)
          mag_1_sample, mag_2_sample = f_1(mass_init_sample), f_2(mass_init_sample) 
          return mag_1_sample + distance_modulus, mag_2_sample + distance_modulus
docs: Simulate a set of stellar magnitudes (no uncertainty) for a
          satellite of a given stellar mass and distance.
  
          Parameters:
          -----------
          stellar_mass : the total stellar mass of the system (Msun)
          distance_modulus : distance modulus of the system (if None takes from isochrone)
          kwargs : passed to iso.imf.sample
  
          Returns:
          --------
          mag_1, mag_2 : simulated magnitudes with length stellar_mass/iso.stellar_mass()