code: def query(self, coords):
          
          # gal = coords.transform_to('galactic')
          gal = coords
          l = gal.l.deg
          b = gal.b.deg
  
          # Detect scalar input
          scalar_input = not hasattr(l, '__len__')
          if scalar_input:
              l = np.array([l])
              b = np.array([b])
  
          # Fill return array with NaNs
          ebv = np.empty(l.shape, dtype='f8')
          ebv[:] = np.nan
  
          # Fill northern cap
          idx = (b >= 65.) & (b <= 90.)
          ebv[idx] = self._lb2ebv_northcap(l[idx], b[idx])
  
          # Fill southern cap
          idx = (b <= -65.) & (b >= -90.)
          ebv[idx] = self._lb2ebv_southcap(l[idx], b[idx])
  
          # Fill northern midplane
          idx = (b < 65.) & (b >= 10.)
          ebv[idx] = self._lb2ebv_midnorth(l[idx], b[idx])
  
          # Fill southern midplane
          idx = (b > -65.) & (b <= -10.)
          ebv[idx] = self._lb2ebv_midsouth(l[idx], b[idx])
  
          if scalar_input:
              ebv = ebv[0]
  
          return ebv
docs: Returns E(B-V) at the specified location(s) on the sky.
          
          Args:
              coords (`astropy.coordinates.SkyCoord`): The coordinates to query.
  
          Returns:
              A float array of reddening, in units of E(B-V), at the given
              coordinates. The shape of the output is the same as the shape of the
              coordinates stored by `coords`.