code: def find_surface_sites_by_height(self, slab, height=0.9, xy_tol=0.05):
          
  
          # Get projection of coordinates along the miller index
          m_projs = np.array([np.dot(site.coords, self.mvec)
                              for site in slab.sites])
  
          # Mask based on window threshold along the miller index.
          mask = (m_projs - np.amax(m_projs)) >= -height
          surf_sites = [slab.sites[n] for n in np.where(mask)[0]]
          if xy_tol:
              # sort surface sites by height
              surf_sites = [s for (h, s) in zip(m_projs[mask], surf_sites)]
              surf_sites.reverse()
              unique_sites, unique_perp_fracs = [], []
              for site in surf_sites:
                  this_perp = site.coords - np.dot(site.coords, self.mvec)
                  this_perp_frac = slab.lattice.get_fractional_coords(this_perp)
                  if not in_coord_list_pbc(unique_perp_fracs, this_perp_frac):
                      unique_sites.append(site)
                      unique_perp_fracs.append(this_perp_frac)
              surf_sites = unique_sites
  
          return surf_sites
docs: This method finds surface sites by determining which sites are within
          a threshold value in height from the topmost site in a list of sites
  
          Args:
              site_list (list): list of sites from which to select surface sites
              height (float): threshold in angstroms of distance from topmost
                  site in slab along the slab c-vector to include in surface
                  site determination
              xy_tol (float): if supplied, will remove any sites which are
                  within a certain distance in the miller plane.
  
          Returns:
              list of sites selected to be within a threshold of the highest