diff --git a/examples/anisotropy_3d_polar.py b/examples/anisotropy_3d_polar.py
index 409bc74..d4696ea 100644
--- a/examples/anisotropy_3d_polar.py
+++ b/examples/anisotropy_3d_polar.py
@@ -7,7 +7,6 @@
 
 
 def main():
-
     theta = np.linspace(0, np.pi, 17)
     phi = np.linspace(0, 2 * np.pi, 32)
 
@@ -25,7 +24,7 @@ def main():
         time=time,
         theta=theta,
         phi=phi,
-        B=B0,
+        B0=B0,
         D=0,
         J=0,
         kinetics=[rp.kinetics.Exponential(k)],
diff --git a/radicalpy/classical.py b/radicalpy/classical.py
index 8a6ab67..e7cae64 100644
--- a/radicalpy/classical.py
+++ b/radicalpy/classical.py
@@ -3,6 +3,7 @@
 
 from typing import Tuple
 
+import graphviz
 import numpy as np
 import scipy as sp
 
diff --git a/radicalpy/estimations.py b/radicalpy/estimations.py
index 6fc978d..5e4353b 100644
--- a/radicalpy/estimations.py
+++ b/radicalpy/estimations.py
@@ -151,30 +151,42 @@ def autocorrelation_fit(
     tau_begin: float,
     tau_end: float,
     num_exp: int = 100,
+    normalise: bool = False,
 ) -> dict:
     """Fit multiexponential to autocorrelation plot and calculate the
     effective rotational correlation time.
 
     Args:
+
         ts (np.ndarray): Time interval (x-axis of the `trajectory`)
             (s).
+
         trajectory (np.ndarray): The raw data which will be fit and
             used to calculate `tau_c`.
+
         tau_begin (float): Initial lag time (s).
+
         tau_end (float): Final lag time (s).
+
         num_exp (int): Number of exponential terms in the
             multiexponential fit (default=100).
 
+        normalise (bool): When set to true, the autocorrelation is
+            normalised (default=False).
+
     Returns:
-            dict:
-                - `fit` is the multiexponential fit to the autocorrelation.
-                - `tau_c` is the effective rotational correlation time.
+        dict:
+
+        - `fit` is the multiexponential fit to the autocorrelation.
+        - `tau_c` is the effective rotational correlation time.
 
+    Thank you, Gesa Grüning!
     """
     acf = autocorrelation(trajectory)
     zero_point_crossing = np.where(np.diff(np.sign(acf)))[0][0]
     acf = acf[0:zero_point_crossing]
-    acf /= acf[0]
+    if normalise:
+        acf /= acf[0]
     taus = np.geomspace(tau_begin, tau_end, num=num_exp)
 
     def multiexponential(x, *params):
diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py
index 4ca2140..2b305e9 100644
--- a/radicalpy/simulation.py
+++ b/radicalpy/simulation.py
@@ -287,8 +287,6 @@ def projection_operator(self, state: State):
                 Projection operator corresponding to the `State`
                 `state`.
 
-        .. todo:: Write proper docs.
-
         """
         # Spin operators
         SAx, SAy, SAz = [self.spin_operator(0, ax) for ax in "xyz"]
@@ -678,8 +676,8 @@ def product_probability(self, obs: State, rhos: np.ndarray) -> np.ndarray:
     @staticmethod
     def product_yield(product_probability, time, k):
         """Calculate the product yield and the product yield sum."""
-        product_yield = sp.integrate.cumtrapz(product_probability, time, initial=0) * k
-        product_yield_sum = np.trapz(product_probability, dx=time[1]) * k
+        product_yield = k * sp.integrate.cumtrapz(product_probability, time, initial=0)
+        product_yield_sum = k * np.trapz(product_probability, dx=time[1])
         return product_yield, product_yield_sum
 
     def apply_liouville_hamiltonian_modifiers(self, H, modifiers):
@@ -797,17 +795,50 @@ def anisotropy_loop(
         self,
         init_state: State,
         time: np.ndarray,
-        B: float,
+        B0: float,
         H_base: np.ndarray,
-        theta: float | np.ndarray,
-        phi: float | np.ndarray,
-    ):
+        theta: np.ndarray,
+        phi: np.ndarray,
+    ) -> np.ndarray:
+        """Inner loop of anisotropy experiment.
+
+        Args:
+
+            init_state (State): Initial `State` of the density matrix.
+
+            time (np.ndarray): An sequence of (uniform) time points,
+                usually created using `np.arange` or `np.linspace`.
+
+            B0 (float): External magnetic field intensity (milli
+                Tesla) (see `zeeman_hamiltonian`).
+
+            H_base (np.ndarray): A "base" Hamiltonian, i.e., the
+                Zeeman Hamiltonian will be added to this base, usually
+                obtained with `total_hamiltonian` and `B0=0`.
+
+            theta (np.ndarray): rotation (polar) angle between the
+                external magnetic field and the fixed molecule. See
+                `zeeman_hamiltonian_3d`.
+
+            phi (np.ndarray): rotation (azimuth) angle between the
+                external magnetic field and the fixed molecule. See
+                `zeeman_hamiltonian_3d`.
+
+        Returns:
+            np.ndarray:
+
+            A tensor which has a series of density matrices for each
+            angle `theta` and `phi` obtained by running
+            `time_evolution` for each of them (with `time`
+            time\-steps, `B0` magnetic intensity).
+
+        """
         shape = self._get_rho_shape(H_base.shape[0])
         rhos = np.zeros((len(theta), len(phi), len(time), *shape), dtype=complex)
 
         iters = itertools.product(enumerate(theta), enumerate(phi))
         for (i, th), (j, ph) in tqdm(list(iters)):
-            H_zee = self.zeeman_hamiltonian(B, th, ph)
+            H_zee = self.zeeman_hamiltonian(B0, th, ph)
             H = H_base + self.convert(H_zee)
             rhos[i, j] = self.time_evolution(init_state, time, H)
         return rhos
@@ -819,17 +850,70 @@ def anisotropy(
         time: np.ndarray,
         theta: np.ndarray | float,
         phi: np.ndarray | float,
-        B: float,
+        B0: float,
         D: np.ndarray,
         J: float,
         kinetics: list[HilbertIncoherentProcessBase] = [],
         relaxations: list[HilbertIncoherentProcessBase] = [],
     ) -> dict:
+        """Anisotropy experiment.
+
+        Args:
+
+            init_state (State): Initial `State` of the density matrix.
+
+            obs_state (State): Observable `State` of the density matrix.
+
+            time (np.ndarray): An sequence of (uniform) time points,
+                usually created using `np.arange` or `np.linspace`.
+
+            H_base (np.ndarray): A "base" Hamiltonian, i.e., the
+                Zeeman Hamiltonian will be added to this base, usually
+                obtained with `total_hamiltonian` and `B0=0`.
+
+            theta (np.ndarray): rotation (polar) angle between the
+                external magnetic field and the fixed molecule. See
+                `zeeman_hamiltonian_3d`.
+
+            B0 (float): External magnetic field intensity (milli
+                Tesla) (see `zeeman_hamiltonian`).
+
+            phi (np.ndarray): rotation (azimuth) angle between the
+                external magnetic field and the fixed molecule. See
+                `zeeman_hamiltonian_3d`.
+
+            D (np.ndarray): Dipolar coupling constant (see
+                `dipolar_hamiltonian`).
+
+            J (float): Exchange coupling constant (see
+                `exchange_hamiltonian`).
+
+            kinetics (list): A list of kinetic (super)operators of
+                type `radicalpy.kinetics.HilbertKineticsBase` or
+                `radicalpy.kinetics.LiouvilleKineticsBase`.
+
+            relaxations (list): A list of relaxation superoperators of
+                type `radicalpy.relaxation.LiouvilleRelaxationBase`.
+
+        Returns:
+            dict:
+
+            - time: the original `time` object
+            - B0: `B0` parameter
+            - theta: `theta` parameter
+            - phi: `phi` parameter
+            - rhos: tensor of sequences of time evolution of density
+              matrices
+            - time_evolutions: product probabilities
+            - product_yields: product yields
+            - product_yield_sums: product yield sums
+
+        """
         H = self.total_hamiltonian(B0=0, D=D, J=J, hfc_anisotropy=True)
 
         self.apply_liouville_hamiltonian_modifiers(H, kinetics + relaxations)
         theta, phi = utils.anisotropy_check(theta, phi)
-        rhos = self.anisotropy_loop(init_state, time, B, H, theta=theta, phi=phi)
+        rhos = self.anisotropy_loop(init_state, time, B0, H, theta=theta, phi=phi)
         product_probabilities = self.product_probability(obs_state, rhos)
 
         self.apply_hilbert_kinetics(time, product_probabilities, kinetics)
@@ -841,7 +925,7 @@ def anisotropy(
 
         return dict(
             time=time,
-            B=B,
+            B0=B0,
             theta=theta,
             phi=phi,
             rhos=rhos,
diff --git a/requirements.txt b/requirements.txt
index 6aeb7fb..9f155ca 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -7,3 +7,4 @@ scipy
 seaborn
 sympy
 tqdm
+graphviz