docs: def sample(self, mode='data', mass_steps=1000, mass_min=0.1, full_data_range=False):
          """Sample the isochrone in steps of mass interpolating between
          the originally defined isochrone points.
  
          Parameters:
          -----------
          mode : 
          mass_steps : 
          mass_min : Minimum mass [Msun]
          full_data_range :
          
          Returns:
          --------
          mass_init : Initial mass of each point
          mass_pdf : PDF of number of stars in each point
          mass_act : Actual (current mass) of each stellar point
          mag_1 : Array of absolute magnitudes in first band (no distance modulus applied)
          mag_2 : Array of absolute magnitudes in second band (no distance modulus applied)
          """
  
          if full_data_range:
              # ADW: Might be depricated 02/10/2015
              # Generate points over full isochrone data range
              select = slice(None)
          else:
              # Not generating points for the post-AGB stars,
              # but still count those stars towards the normalization
              select = slice(self.index)
  
          mass_steps = int(mass_steps)
  
          mass_init = self.mass_init[select]
          mass_act = self.mass_act[select]
          mag_1 = self.mag_1[select]
          mag_2 = self.mag_2[select]
          
          # ADW: Assume that the isochrones are pre-sorted by mass_init
          # This avoids some numerical instability from points that have the same
          # mass_init value (discontinuities in the isochrone).
          # ADW: Might consider using np.interp for speed
          mass_act_interpolation = scipy.interpolate.interp1d(mass_init, mass_act,assume_sorted=True)
          mag_1_interpolation = scipy.interpolate.interp1d(mass_init, mag_1,assume_sorted=True)
          mag_2_interpolation = scipy.interpolate.interp1d(mass_init, mag_2,assume_sorted=True)
  
          # ADW: Any other modes possible?
          if mode=='data':
              # Mass interpolation with uniform coverage between data points from isochrone file 
              mass_interpolation = scipy.interpolate.interp1d(np.arange(len(mass_init)), mass_init)
              mass_array = mass_interpolation(np.linspace(0, len(mass_init)-1, mass_steps+1))
              d_mass = mass_array[1:] - mass_array[:-1]
              mass_init_array = np.sqrt(mass_array[1:] * mass_array[:-1])
              mass_pdf_array = d_mass * self.imf.pdf(mass_init_array, log_mode=False)
              mass_act_array = mass_act_interpolation(mass_init_array)
              mag_1_array = mag_1_interpolation(mass_init_array)
              mag_2_array = mag_2_interpolation(mass_init_array)
  
          # Horizontal branch dispersion
          if self.hb_spread and (self.stage==self.hb_stage).any():
              logger.debug("Performing dispersion of horizontal branch...")
              mass_init_min = self.mass_init[self.stage==self.hb_stage].min()
              mass_init_max = self.mass_init[self.stage==self.hb_stage].max()
              cut = (mass_init_array>mass_init_min)&(mass_init_array<mass_init_max)
              if isinstance(self.hb_spread,collections.Iterable):
                  # Explicit dispersion spacing
                  dispersion_array = self.hb_spread
                  n = len(dispersion_array)
              else:
                  # Default dispersion spacing
                  dispersion = self.hb_spread
                  spacing = 0.025
                  n = int(round(2.0*self.hb_spread/spacing))
                  if n % 2 != 1: n += 1
                  dispersion_array = np.linspace(-dispersion, dispersion, n)
  
              # Reset original values
              mass_pdf_array[cut] = mass_pdf_array[cut] / float(n)
  
              # Isochrone values for points on the HB
              mass_init_hb = mass_init_array[cut]
              mass_pdf_hb = mass_pdf_array[cut]
              mass_act_hb = mass_act_array[cut]
              mag_1_hb = mag_1_array[cut]
              mag_2_hb = mag_2_array[cut]
  
              # Add dispersed values
              for dispersion in dispersion_array:
                  if dispersion == 0.: continue
                  msg = 'Dispersion=%-.4g, HB Points=%i, Iso Points=%i'%(dispersion,cut.sum(),len(mass_init_array))
                  logger.debug(msg)
  
                  mass_init_array = np.append(mass_init_array, mass_init_hb) 
                  mass_pdf_array = np.append(mass_pdf_array, mass_pdf_hb)
                  mass_act_array = np.append(mass_act_array, mass_act_hb) 
                  mag_1_array = np.append(mag_1_array, mag_1_hb + dispersion)
                  mag_2_array = np.append(mag_2_array, mag_2_hb + dispersion)
  
          # Note that the mass_pdf_array is not generally normalized to unity
          # since the isochrone data range typically covers a different range
          # of initial masses
          #mass_pdf_array /= np.sum(mass_pdf_array) # ORIGINAL
          # Normalize to the number of stars in the satellite with mass > mass_min
          mass_pdf_array /= self.imf.integrate(mass_min, self.mass_init_upper_bound)
          out = np.vstack([mass_init_array,mass_pdf_array,mass_act_array,mag_1_array,mag_2_array])
          return out