docs: def equation(self):
          """Mix-in class that returns matrix rows for leaky wall condition.
          Qnormal = resfac * (headin - headout)
          Returns matrix part (nunknowns,neq)
          Returns rhs part nunknowns
          """
          mat = np.empty((self.nunknowns, self.model.neq))
          rhs = np.zeros(self.nunknowns)  # Needs to be initialized to zero
          for icp in range(self.ncp):
              istart = icp * self.nlayers
              ieq = 0
              for e in self.model.elementlist:
                  if e.nunknowns > 0:
                      qx, qy = e.disvecinflayers(self.xc[icp], self.yc[icp], self.layers)
                      mat[istart:istart + self.nlayers, ieq:ieq + e.nunknowns] = \
                          qx * self.cosnorm[icp] + qy * self.sinnorm[icp] - self.resfac[:, np.newaxis] * \
                          (e.potinflayers(self.xcin[icp], self.ycin[icp], self.layers, aq=self.aq) / self.aq.Tcol[
                              self.layers] - \
                           e.potinflayers(self.xcout[icp], self.ycout[icp], self.layers, aq=self.aq) / self.aq.Tcol[
                               self.layers])
                      ieq += e.nunknowns
                  else:
                      qx, qy = e.disveclayers(self.xc[icp], self.yc[icp], self.layers)
                      rhs[istart:istart + self.nlayers] -= qx * self.cosnorm[icp] + qy * self.sinnorm[icp] + self.resfac * \
                                                            (e.potentiallayers(self.xcin[icp], self.ycin[icp], self.layers,
                                                                               aq=self.aq) / self.aq.T[self.layers] -
                                                             e.potentiallayers(self.xcout[icp], self.ycout[icp],
                                                                               self.layers, aq=self.aq) / self.aq.T[
                                                                self.layers])
          return mat, rhs