docs: def headalongline(self, x, y, layers=None):
          """Head along line or curve
          
          Parameters
          ----------
          x : array
              x values of line
          y : array
              y values of line
          layers : integer, list or array, optional
              layers for which grid is returned
          
          Returns
          -------
          h : array size `nlayers, nx`
  
          """
          
          xg, yg = np.atleast_1d(x), np.atleast_1d(y)
          if layers is None:
              Nlayers = self.aq.find_aquifer_data(xg[0], yg[0]).naq
          else:
              Nlayers = len(np.atleast_1d(layers))
          nx = len(xg)
          if len(yg) == 1:
              yg = yg * np.ones(nx)
          h = np.zeros((Nlayers, nx))
          for i in range(nx):
              h[:, i] = self.head(xg[i], yg[i], layers)
          return h