docs: def separation(self, mag_1, mag_2):
          """ 
          Calculate the separation between a specific point and the
          isochrone in magnitude-magnitude space. Uses an interpolation
  
          ADW: Could speed this up...
  
          Parameters:
          -----------
          mag_1 : The magnitude of the test points in the first band
          mag_2 : The magnitude of the test points in the second band
  
          Returns:
          --------
          sep : Minimum separation between test points and isochrone interpolation
          """
  
          iso_mag_1 = self.mag_1 + self.distance_modulus
          iso_mag_2 = self.mag_2 + self.distance_modulus
          
          def interp_iso(iso_mag_1,iso_mag_2,mag_1,mag_2):
              interp_1 = scipy.interpolate.interp1d(iso_mag_1,iso_mag_2,bounds_error=False)
              interp_2 = scipy.interpolate.interp1d(iso_mag_2,iso_mag_1,bounds_error=False)
  
              dy = interp_1(mag_1) - mag_2
              dx = interp_2(mag_2) - mag_1
  
              dmag_1 = np.fabs(dx*dy) / (dx**2 + dy**2) * dy
              dmag_2 = np.fabs(dx*dy) / (dx**2 + dy**2) * dx
  
              return dmag_1, dmag_2
  
          # Separate the various stellar evolution stages
          if np.issubdtype(self.stage.dtype,np.number):
              sel = (self.stage < self.hb_stage)
          else:
              sel = (self.stage != self.hb_stage)
  
          # First do the MS/RGB
          rgb_mag_1 = iso_mag_1[sel]
          rgb_mag_2 = iso_mag_2[sel]
          dmag_1,dmag_2 = interp_iso(rgb_mag_1,rgb_mag_2,mag_1,mag_2)
  
          # Then do the HB (if it exists)
          if not np.all(sel):
              hb_mag_1 = iso_mag_1[~sel]
              hb_mag_2 = iso_mag_2[~sel]
  
              hb_dmag_1,hb_dmag_2 = interp_iso(hb_mag_1,hb_mag_2,mag_1,mag_2)
  
              dmag_1 = np.nanmin([dmag_1,hb_dmag_1],axis=0)
              dmag_2 = np.nanmin([dmag_2,hb_dmag_2],axis=0)
  
          #return dmag_1,dmag_2
          return np.sqrt(dmag_1**2 + dmag_2**2)