diff --git a/doc/api.rst b/doc/api.rst
index 957e50a21fc2af3199323ea18f9b78abf6298492..f095c3ca8097a624a3d9115345dc7377dd247c0c 100644
--- a/doc/api.rst
+++ b/doc/api.rst
@@ -26,9 +26,7 @@ API
    :members:
    :undoc-members:
 
-.. py:data:: usadelndsoc.logger
-
-   `logging` instance used.
+.. autodata:: usadelndsoc.logger
 
 .. automodule:: usadelndsoc.plotting
    :members:
diff --git a/usadelndsoc/__init__.py b/usadelndsoc/__init__.py
index 1459a04734ab1548b6cbe40687c9cd35673e872b..d083571709dd8563a4249d5b76d57bdad7337741 100644
--- a/usadelndsoc/__init__.py
+++ b/usadelndsoc/__init__.py
@@ -5,6 +5,9 @@ import contextlib
 
 
 logger = logging.getLogger(__name__)
+"""
+`logging` instance used in usadelndsoc.
+"""
 
 
 def _init_logging():
diff --git a/usadelndsoc/core.py b/usadelndsoc/core.py
index d8c4e8ac81352c5cb2dd7b8958724c79616de941..f1899c851c5963c01d55663384e8c16f84146444 100644
--- a/usadelndsoc/core.py
+++ b/usadelndsoc/core.py
@@ -10,7 +10,7 @@ __all__ = ["Core", "MASK_NONE", "MASK_TERMINAL", "MASK_VACUUM"]
 _log = logging.getLogger(__name__)
 
 
-def _action_property(name, freeze=True):
+def _action_property(name, freeze=True, doc=None):
     def fget(self):
         return getattr(self, "_" + name)
 
@@ -25,7 +25,7 @@ def _action_property(name, freeze=True):
             setattr(self._action, name, value)
         setattr(self, "_" + name, value)
 
-    return property(fget, fset)
+    return property(fget, fset, doc=doc)
 
 
 class Core:
@@ -84,10 +84,16 @@ class Core:
 
     @property
     def x(self):
+        """
+        x-coordinates of the 2D cell grid
+        """
         return np.linspace(-self.Lx / 2, self.Lx / 2, self._shape[0])
 
     @property
     def y(self):
+        """
+        y-coordinates of the 2D cell grid
+        """
         return np.linspace(-self.Ly / 2, self.Ly / 2, self._shape[1])
 
     def eval(self, Q):
@@ -114,16 +120,33 @@ class Core:
         self._ensure_action()
         return self._action.grad_Omega(Q)
 
-    U = _action_property("U")
-    Omega = _action_property("Omega")
-    mask = _action_property("mask")
-    omega = _action_property("omega", freeze=False)
+    U = _action_property(
+        "U",
+        doc="""Wilson link matrices. Array, shape (nx, ny, 4(direction), 4, 4).
+        The direction index is one of RIGHT(0)/UP(1)/LEFT(2)/DOWN(3).
+        """,
+    )
+    Omega = _action_property(
+        "Omega", doc="Local energy matrix. Array, shape (nx, ny, 4, 4)."
+    )
+    mask = _action_property(
+        "mask", doc="""Material mask. Array of int, shape (nx, ny)"""
+    )
+    omega = _action_property(
+        "omega",
+        freeze=False,
+        doc="""Local energy matrix. Array of complex, shape (nx, nx, 4, 4)""",
+    )
     hess_computed = _action_property("hess_computed")
-    shape = property(fget=lambda self: self._shape)
-    Lx = _action_property("Lx")
-    Ly = _action_property("Ly")
-    D = _action_property("D")
-    eta = _action_property("eta")
+    shape = property(
+        fget=lambda self: self._shape, doc="""Array size of space, (nx, ny)"""
+    )
+    Lx = _action_property("Lx", doc="""Total length of system in x-direction, float.""")
+    Ly = _action_property("Ly", doc="""Total length of system in y-direction, float.""")
+    D = _action_property(
+        "D", doc="""Second gradient term prefactor (diffusion constant), float."""
+    )
+    eta = _action_property("eta", doc="""Hall term prefactor constant, float.""")
 
     @property
     def hess_computed(self):
@@ -131,14 +154,20 @@ class Core:
 
     @property
     def Delta(self):
+        """
+        Superconducting order parameter part of `Omega`. Modifiable view.
+        Array of complex, shape (nx, nx, 2, 2).
+        """
         return _DeltaArray(self.Omega)
 
     @property
     def Ux(self):
+        """Modifiable view to Wilson link matrix array, for LEFT"""
         return self.U[:, :, 0, :, :]
 
     @property
     def Uy(self):
+        """Modifiable view to Wilson link matrix array, for UP"""
         return self.U[:, :, 1, :, :]
 
 
diff --git a/usadelndsoc/solver.py b/usadelndsoc/solver.py
index 86883303a68bb02bb0082267d86054cb2442b034..61d9511dd771a5905b1f73c4726a6f3a8281595d 100644
--- a/usadelndsoc/solver.py
+++ b/usadelndsoc/solver.py
@@ -41,10 +41,10 @@ __all__ = [
     "Solver",
     "Result",
     "cpr",
+    "tr",
     "MASK_NONE",
     "MASK_TERMINAL",
     "MASK_VACUUM",
-    "tr",
     "LEFT",
     "UP",
     "RIGHT",
@@ -65,26 +65,56 @@ S_z = np.array([[1, 0], [0, -1]])
 S_0 = np.array([[1, 0], [0, 1]])
 
 RIGHT = 0
+"""
+Direction index used for right.
+"""
+
 UP = 1
+"""
+Direction index used for up.
+"""
+
 LEFT = 2
+"""
+Direction index used for left.
+"""
+
 DOWN = 3
+"""
+Direction index used for down.
+"""
+
+MASK_NONE = MASK_NONE
+"""
+Material mask indicating free material cell.
+"""
+
+MASK_TERMINAL = MASK_TERMINAL
+"""
+Material mask indicating rigid terminal boundary condition cell.
+"""
 
+MASK_VACUUM = MASK_VACUUM
+"""
+Material mask indicating vacuum cell.
+"""
 
-def _core_property(name):
+
+def _core_property(name, doc=None):
     def fget(self):
         return getattr(self._core, name)
 
     def fset(self, value):
         setattr(self._core, name, value)
 
-    return property(fget, fset)
+    return property(fget, fset, doc=doc or getattr(Core, name).__doc__)
 
 
-def _array_property(name):
+def _array_property(name, doc=None):
     def fget(self):
         return getattr(self, name)
 
-    return property(fget)
+    return property(fget, doc=doc)
 
 
 def _norm(z):
@@ -116,6 +146,9 @@ class Solver:
         self._t_solve_min = np.inf
 
     def set_shape(self, nx, ny=1):
+        """
+        Set system spatial size in cells.
+        """
         self._core.set_shape(nx, ny)
         self._Phi = np.zeros((nx, ny, 2, 2, 2), dtype=np.complex_)
         self._J_tot = np.zeros((nx, ny, 4, 4, 4), dtype=np.complex_)
@@ -135,10 +168,22 @@ class Solver:
     Ux = _core_property("Ux")
     Uy = _core_property("Uy")
 
-    Phi = _array_property("_Phi")
-    J_tot = _array_property("_J_tot")
+    Phi = _array_property(
+        "_Phi",
+        doc="""Current saddle-point solution in Riccati parametrization. Array of complex, shape (nx, ny, 2, 2, 2).""",
+    )
+    J_tot = _array_property(
+        "_J_tot",
+        doc="""Matsubara-summed matrix current between cells. Array of complex, shape (nx, ny, 4(direction), 4, 4)""",
+    )
 
     def reset(self):
+        """Reset solver.
+
+        Allows modification of material mask and other parameters,
+        that become frozen when solver is run.
+
+        """
         Solver.__global_id += 1
         self._id = Solver.__global_id
         self._core.reset()
@@ -615,6 +660,46 @@ class Solver:
         return Result(self, omega=omega0, Phi=Phi)
 
     def self_consistency(self, T, T_c0, continuation=False, **kw):
+        """
+        Solve self-consistency equations for Delta.
+
+        Assumes simple singlet weak coupling.
+
+        Supports additional constraint variables and equations.
+
+        Parameters
+        ----------
+        T : float
+            Temperature
+        T_c0 : float or array of float, shape (nx, ny)
+            Critical temperature.
+        continuation : bool, optional
+            Whether to continue from previous solution for Delta.
+        maxiter : int
+        tol : float
+        xtol : float
+        workers : int
+        constraint_fun : function(Delta, x) -> residual
+            Function returning additional components of the self-consistency
+            residual, where `x` are additional constraint variables.
+            This can be used e.g. for pseudoarclength continuation.
+        constraint_x : array
+            Initial values of the additional constraint variables.
+        plot : bool
+        **solver_kw
+
+        Returns
+        -------
+        Delta
+            Order parameter matrix
+        J_tot
+            Matsubara-summed matrix current
+        cx
+            Final constraint variable values
+        succes : bool
+            Whether convergence was obtained.
+
+        """
         nx, ny = self.shape
         T_c0 = np.broadcast_to(T_c0, (nx, ny))
         self.Delta[self.mask == MASK_VACUUM] = 0