docs: def generate(self, trilegal_filename, ra=None, dec=None,
                   n=2e4, ichrone='mist', MAfn=None,
                   mags=None, maxrad=None, f_binary=0.4, **kwargs):
          """
          Generate population.
          """
          n = int(n)
  
          #generate/load BG primary stars from TRILEGAL simulation
          bgpop = BGStarPopulation_TRILEGAL(trilegal_filename,
                                          ra=ra, dec=dec, mags=mags,
                                          maxrad=maxrad, **kwargs)
  
          # Make sure that
          # properties of stars are within allowable range for isochrone.
          # This is a bit hacky, admitted.
          mass = bgpop.stars['m_ini'].values
          age = bgpop.stars['logAge'].values
          feh = bgpop.stars['[M/H]'].values
  
          ichrone = get_ichrone(ichrone)
  
          pct = 0.05 #pct distance from "edges" of ichrone interpolation
          mass[mass < ichrone.minmass*(1+pct)] = ichrone.minmass*(1+pct)
          mass[mass > ichrone.maxmass*(1-pct)] = ichrone.maxmass*(1-pct)
          age[age < ichrone.minage*(1+pct)] = ichrone.minage*(1+pct)
          age[age > ichrone.maxage*(1-pct)] = ichrone.maxage*(1-pct)
          feh[feh < ichrone.minfeh+0.05] = ichrone.minfeh+0.05
          feh[feh > ichrone.maxfeh-0.05] = ichrone.maxfeh-0.05
  
          distance = bgpop.stars['distance'].values
  
          #Generate binary population to draw eclipses from
          pop = MultipleStarPopulation(mA=mass, age=age, feh=feh,
                                              f_triple=0, f_binary=1,
                                              distance=distance,
                                              ichrone=ichrone)
  
          all_stars = pop.stars.dropna(subset=['mass_A'])
          all_stars.reset_index(inplace=True)
  
          #generate eclipses
          stars = pd.DataFrame()
          df_orbpop = pd.DataFrame()
          tot_prob = None; tot_dprob=None; prob_norm=None
  
          n_adapt = n
          while len(stars) < n:
              n_adapt = int(n_adapt)
              inds = np.random.randint(len(all_stars), size=n_adapt)
  
              s = all_stars.iloc[inds]
  
              #calculate limb-darkening coefficients
              u1A, u2A = ldcoeffs(s['Teff_A'], s['logg_A'])
              u1B, u2B = ldcoeffs(s['Teff_B'], s['logg_B'])
  
              inds, df, (prob,dprob) = calculate_eclipses(s['mass_A'], s['mass_B'],
                                                          s['radius_A'], s['radius_B'],
                                                          s['{}_mag_A'.format(self.band)],
                                                          s['{}_mag_B'.format(self.band)],
                                                          u11s=u1A, u21s=u2A,
                                                          u12s=u1B, u22s=u2B,
                                                          band=self.band,
                                                          period=self.period,
                                                          calc_mininc=True,
                                                          return_indices=True,
                                                          MAfn=MAfn)
              s = s.iloc[inds].copy()
              s.reset_index(inplace=True)
              for col in df.columns:
                  s[col] = df[col]
              stars = pd.concat((stars, s))
  
              #new_df_orbpop = pop.orbpop.orbpop_long.dataframe.iloc[inds].copy()
              #new_df_orbpop.reset_index(inplace=True)
  
              #df_orbpop = pd.concat((df_orbpop, new_df_orbpop))
  
              logging.info('{} BEB systems generated (target {})'.format(len(stars),n))
              #logging.debug('{} nans in stars[dpri]'.format(np.isnan(stars['dpri']).sum()))
              #logging.debug('{} nans in df[dpri]'.format(np.isnan(df['dpri']).sum()))
  
              if tot_prob is None:
                  prob_norm = (1/dprob**2)
                  tot_prob = prob
                  tot_dprob = dprob
              else:
                  prob_norm = (1/tot_dprob**2 + 1/dprob**2)
                  tot_prob = (tot_prob/tot_dprob**2 + prob/dprob**2)/prob_norm
                  tot_dprob = 1/np.sqrt(prob_norm)
  
              n_adapt = min(int(1.2*(n-len(stars)) * n_adapt//len(s)), 5e5)
              #logging.debug('n_adapt = {}'.format(n_adapt))
              n_adapt = max(n_adapt, 100)
              n_adapt = int(n_adapt)
  
          stars = stars.iloc[:n]
  
          if 'level_0' in stars:
              stars.drop('level_0', axis=1, inplace=True) #dunno where this came from
          stars = stars.reset_index()
          stars.drop('index', axis=1, inplace=True)
  
          stars['mass_1'] = stars['mass_A']
          stars['radius_1'] = stars['radius_A']
          stars['mass_2'] = stars['mass_B']
          stars['radius_2'] = stars['radius_B']
  
          MultipleStarPopulation.__init__(self, stars=stars,
                                          #orbpop=orbpop,
                                          f_triple=0, f_binary=f_binary,
                                          period_long=self.period)
  
          priorfactors = {'f_binary':f_binary}
  
          #attributes needed for BGStarPopulation
          self.density = bgpop.density
          self.trilegal_args = bgpop.trilegal_args
          self._maxrad = bgpop._maxrad
  
          #create an OrbitPopulation here?
  
          EclipsePopulation.__init__(self, stars=stars, #orbpop=orbpop,
                                     period=self.period, cadence=self.cadence,
                                     model=self.model,
                                     lhoodcachefile=self.lhoodcachefile,
                                     priorfactors=priorfactors, prob=tot_prob)
  
          #add Rsky property
          self.stars['Rsky'] = randpos_in_circle(len(self.stars),
                                                 self._maxrad, return_rad=True)