diff --git a/.gitignore b/.gitignore
index c2d15588c5c91608146516ca8fc1b9b357776d69..fbc8936473fe27ba37cc4ee2d4709401cbd35776 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,3 +10,8 @@ docs/source/example_gallery/**/*.ipynb
 build/
 dist/
 tests/test_output/
+
+
+pylint/
+htmlcov/
+.coverage
\ No newline at end of file
diff --git a/README.rst b/README.rst
index 086f06ff165acd60b6906617b55c8f6c24fcaaca..d71064633c58c137806d2a4428c4c3a136b39bf1 100644
--- a/README.rst
+++ b/README.rst
@@ -1,4 +1,4 @@
-.. image:: docs/source/_static/Logo_cmrsim.svg
+.. image:: https://people.ee.ethz.ch/~jweine/cmrsim/latest/_images/Logo_cmrsim_moving_light.svg
    :align: right
    :scale: 150 %
 
@@ -22,14 +22,15 @@ What makes CMRsim special?
 - Reproducibility by using versioned and documented open-source software.
 - Easily scalable to detailed digital phantoms
 
+
 Installation
 ==============
 
-
 Pip
 ^^^^
 CMRsim can be installed using the python package manager :code:`pip`. The package can be found
-in the `registry`_  of the repository hosted on the ETH-Gitlab instance
+in the `registry`_  of the repository hosted on the ETH-Gitlab instance as well as on the python
+package index (PyPI).
 
 .. _registry: https://gitlab.ethz.ch/jweine/cmrsim/-/packages
 
@@ -38,30 +39,17 @@ Stable release versions are denoted with the ususal :code:`major.minor` notation
 development versions are constantly updated on solved issued and can be installed using the
 :code:`major.minor.devX` notation
 
-Installation in command line:
+Release versions are published on PyPI, which allows installing cmrsim with:
 
 .. code-block::
 
-    pip install cmrsim --extra-index-url https://__token__:<your_personal_token>@gitlab.ethz.ch/api/v4/projects/23798/packages/pypi/simple
-
+    pip install cmrsim
 
-Requirements
-^^^^^^^^^^^^^^^^
-CMRsim does rely on functionality of 3 non-standard python stack packages:
+Development versions are only available using the extra url:
 
-1. `TensorFlow`_: Pretty much all computationally heavy calculations.
-2. `Pyvista`_: For processing and visualizing graphs (e.g. 3D Velocity fields).
-3. `CMRseq`_: For defining MR-sequences used for Fourier-Encoding calculation and Bloch-simulation.
-
-.. _TensorFlow: https://www.tensorflow.org/
-.. _Pyvista: https://docs.pyvista.org/
-.. _CMRseq: https://people.ee.ethz.ch/~jweine/cmrseq/index.html
-
-Further (already indirectly included) requirements are:
--h5py~=2.10.0
--scipy>=1.4.1
--tqdm
+.. code-block::
 
+    pip install cmrsim --extra-index-url https://gitlab.ethz.ch/api/v4/projects/23798/packages/pypi/simple
 
 Docker
 ^^^^^^^^
diff --git a/ci_cd_config/.run-test-stage-ci.yml b/ci_cd_config/.run-test-stage-ci.yml
index 3712ad07518560c083f3e459a8280411c4116566..b6eb17b041818a86a6be1f981dbfd8a843dd0db8 100644
--- a/ci_cd_config/.run-test-stage-ci.yml
+++ b/ci_cd_config/.run-test-stage-ci.yml
@@ -5,11 +5,17 @@
 run_unit_test:
   image: registry.ethz.ch/jweine/cmrsim:unittest
   stage: test
+  variables:
+    BATCH_TAG: ""
+  rules:
+    - if: $CI_COMMIT_BRANCH
+      variables:
+        BATCH_TAG: ":dev"
   script:
     - coverage run --source cmrsim -m unittest discover tests
     - coverage report
     - coverage xml
-    - genbadge coverage -i coverage.xml -o coverage.svg
+    - genbadge coverage -i coverage.xml -o coverage.svg --name "coverage${BATCH_TAG}"
   artifacts:
     paths:
       - coverage.svg
diff --git a/cmrsim/analytic/contrast/base.py b/cmrsim/analytic/contrast/base.py
index f16c8c0ca6837c2bdce1ac76b75ae0636296c7cb..4e6b59fe1fca9e920f3ba67ece6583fca0363bbf 100644
--- a/cmrsim/analytic/contrast/base.py
+++ b/cmrsim/analytic/contrast/base.py
@@ -250,7 +250,7 @@ class LookUpTableModule(BaseSignalModel):
 
         :param r_vectors: (-1, 3) in meter
         :param batch_index: (#batch) indices of the first dimension (T) used for spatial lookup
-        :return: (#batch_idx, -1, self._n_channels) trilinear interpolated map values for specifiec
+        :return: (#batch_idx, -1, self._n_channels) trilinear interpolated map values for specified
                     locations
         """
         with tf.name_scope("linear_interpolation_lookup") as scope:
diff --git a/cmrsim/analytic/contrast/coil_sensitivities.py b/cmrsim/analytic/contrast/coil_sensitivities.py
index 8a0ec0d393f7397b69cdb415a41b70953f3ffe73..162bf1f1c483963a8a19371a51781e239af6050d 100644
--- a/cmrsim/analytic/contrast/coil_sensitivities.py
+++ b/cmrsim/analytic/contrast/coil_sensitivities.py
@@ -116,6 +116,21 @@ class CoilSensitivity(LookUpTableModule):
             result = tf.reshape(product, new_shape)
             return result
 
+    def lookup_complex_coil_factors(self, r_vectors: tf.Tensor):
+        """ Looks the sensitivity values for all coils at specified locations.
+        
+        :param r_vectors: [..., 3] dtype: tf.float32
+        :return: coil_factors [..., n_coils] dtype: tf.complex64
+        """
+        old_shape = tf.shape(r_vectors)
+        new_shape = tf.concat([old_shape[:-1], [self.expansion_factor, ]], axis=0)
+        r_flat = tf.reshape(r_vectors, [-1, 3])
+        sense_factors = self.linear_interpolation(r_vectors=r_flat)
+        sense_factors = tf.reshape(sense_factors, [-1, self.expansion_factor, 2])
+        sense_factors = tf.complex(sense_factors[..., 0], sense_factors[..., 1])
+        sense_factors = tf.reshape(sense_factors, new_shape)
+        return sense_factors
+
     @property
     # pylint: disable=invalid-name
     def coil_maps(self):
@@ -137,6 +152,46 @@ class CoilSensitivity(LookUpTableModule):
         coils = tf.reshape(coils, [*coil_sensitivities.shape[1:-1], -1])
         self.look_up_map3d.assign(coils)
 
+    @classmethod
+    def from_3d_birdcage(cls, map_shape: Union['np.ndarray', tf.Tensor, Tuple],
+                         map_dimensions: np.ndarray, relative_coil_radius: float = 1.5,
+                         device: str = 'CPU:0', **kwargs):
+        """Lookup for ideal bird-cage coil sensitivities.
+
+        .. note::
+
+            expand_repetitions is always set to True for coil-sensitivities, which means that
+            the output hast a #repetions * #coils dimension on the repetitions axis
+
+        :param map_shape: (X, Y, Z, #coils) array shape of the resulting coil-maps
+        :param map_dimensions: (3, 2) -> [(xlow, xhigh), (ylow, yhigh), (zlow, zhigh)] Extend of
+        :param relative_coil_radius: distance of the birdcage ring from the map center
+            the map in scanner coordinate system. Used to sort iso-chromat positions into bins.
+            This is assumed to be the *length* in meter in each spatial direction
+        :param device: (str) Device on which all operations except for the look-up are
+                        placed (default: GPU:0)
+        :param device_lookup: (str) Device where the lookup table is places (defaults to CPU:0)
+
+        """
+        nx, ny, nz, nc = map_shape
+        c, z, y, x = np.mgrid[:nc, :nz, :ny, :nx]
+        coilx = relative_coil_radius * np.cos(c * (2 * np.pi / nc))
+        coily = relative_coil_radius * np.sin(c * (2 * np.pi / nc))
+        coilz = np.floor(c / nc) - 0.25  # * nc
+        coil_phs = -(c + np.floor(c / nc)) * (2 * np.pi / nc)
+
+        x_co = (x - nx / 2.0) / (nx / 2.0) - coilx
+        y_co = (y - ny / 2.0) / (ny / 2.0) - coily
+        z_co = (z - nz / 2.0) / (nz / 2.0) - coilz
+        rr = (x_co ** 2 + y_co ** 2 + z_co ** 2) ** 0.5
+        phi = np.arctan2(x_co, -y_co) + coil_phs
+        out = (1 / rr) * np.exp(1j * phi)
+
+        rss = sum(abs(out) ** 2, 0) ** 0.5
+        out /= rss
+        coil_maps3d = tf.constant(out.transpose(3, 2, 1, 0), dtype=tf.complex64)
+        return cls(coil_maps3d, map_dimensions, device, **kwargs)
+
 
 class BirdcageCoilSensitivities(CoilSensitivity):
     """ Lookup for ideal bird-cage coil sensitivities.
diff --git a/cmrsim/analytic/contrast/phase_tracking.py b/cmrsim/analytic/contrast/phase_tracking.py
index adcd8938edfeb75ad780dcf80360691154e13dd2..81b704a5aa078f50b63d6d84dbb8b026e6bfe08b 100644
--- a/cmrsim/analytic/contrast/phase_tracking.py
+++ b/cmrsim/analytic/contrast/phase_tracking.py
@@ -78,7 +78,8 @@ class PhaseTracking(BaseSignalModel):
             # Case 2: repetition-axis of signal_tensor must match self.expansion_factor
             else:
                 tf.Assert(input_shape[1] == self.expansion_factor,
-                          ["Shape missmatch for input argument signal_tensor for case no-expand"])
+                          ["Shape missmatch for input argument signal_tensor for case no-expand",
+                           input_shape, self.expansion_factor])
                 result = signal_tensor * complex_phase
             return result
 
diff --git a/cmrsim/analytic/contrast/sequences.py b/cmrsim/analytic/contrast/sequences.py
index c10e6fa2fccf1e45f4bea7f13a7b04f5c9ccc7ef..3720180e0f2fa9e7be468fe0f980d5887ab9759f 100644
--- a/cmrsim/analytic/contrast/sequences.py
+++ b/cmrsim/analytic/contrast/sequences.py
@@ -124,7 +124,7 @@ class BSSFP(BaseSignalModel):
     :param expand_repetitions: See BaseSignalModel for explanation
     """
     # pylint: disable=invalid-name
-    required_quantities = ('T1', 'T2', 'T2star', 'theta')
+    required_quantities = ('T1', 'T2', 'T2star', 'off_res')
     #: Echo time in ms
     TE: tf.Variable = None
     #: Repetition time in ms
@@ -161,7 +161,7 @@ class BSSFP(BaseSignalModel):
 
     # pylint: disable=too-many-arguments
     def __call__(self, signal_tensor: tf.Tensor, T1: tf.Tensor, T2: tf.Tensor,
-                 T2star: tf.Tensor, theta: tf.Tensor, **kwargs):
+                 T2star: tf.Tensor, off_res: tf.Tensor, **kwargs):
         # pylint: disable=invalid-name, too-many-locals,
         """ Evaluates the signal Equation.
 
@@ -169,10 +169,13 @@ class BSSFP(BaseSignalModel):
         :param T1: in ms
         :param T2: in ms
         :param T2star: in ms
-        :param theta: phase in radians acquired over one TR due to off-resonance
+        :param off_res: b0 inhomogeneity in rad/ms (angular frequency difference)
         :return:  (#voxel, #repetitions, #k-space-samples)
         """
 
+        #  phase in radians acquired over one TR due to off-resonance
+        theta = off_res * tf.reshape(self.TR, [1, -1, 1])
+
         T2p = 1. / (1. / T2star - 1./T2)
         E1 = tf.exp(-tf.reshape(self.TR, [1, -1, 1]) / T1)
         E2 = tf.exp(-tf.reshape(self.TR, [1, -1, 1]) / T2)
@@ -183,7 +186,7 @@ class BSSFP(BaseSignalModel):
         b = tf.math.cos(tf.reshape(self.FA, [1, -1, 1])) - E1
 
         # Equation 5
-        Lambda =  (a - E2**2 * b - tf.sqrt((a**2 - E2**2 * b**2) * (1 - E2**2))) / (a - b)
+        Lambda = (a - E2**2 * b - tf.sqrt((a**2 - E2**2 * b**2) * (1 - E2**2))) / (a - b)
 
         # Equation 3
         A = ((a + b) * E2 / 2.) / (a - Lambda * (a - b) / 2.)
diff --git a/cmrsim/analytic/contrast/t2_star.py b/cmrsim/analytic/contrast/t2_star.py
index 9122e87e011d373e2afe676d721061eeb97aef93..72d5507852c5deda3cc350d7d1a338b90c32f083 100644
--- a/cmrsim/analytic/contrast/t2_star.py
+++ b/cmrsim/analytic/contrast/t2_star.py
@@ -94,7 +94,7 @@ class StaticT2starDecay(BaseSignalModel):
 
             if self.expand_repetitions or self.expansion_factor == 1:
                 result = signal_tensor[:, :, tf.newaxis, :] * decay_factor[:, tf.newaxis, :, :]
-                new_shape = tf.concat([input_shape[0], input_shape[1] * self.expansion_factor,
+                new_shape = tf.stack([input_shape[0], input_shape[1] * self.expansion_factor,
                                        n_samples], axis=0)
                 result = tf.reshape(result, new_shape)
             else:
diff --git a/cmrsim/analytic/simulation.py b/cmrsim/analytic/simulation.py
index a3de5877f3acc6d0695a9a324b47ba1c21e83cb0..d046de0710bcd3aca05fe369cf8f6a42e45cda43 100644
--- a/cmrsim/analytic/simulation.py
+++ b/cmrsim/analytic/simulation.py
@@ -162,6 +162,8 @@ class AnalyticSimulation(tf.Module):
         s_of_k_segments = tf.TensorArray(dtype=tf.complex64, size=n_segments, dynamic_size=False)
         self.progress_bar.reset_segments()
 
+        m_transverse = self.forward_model(batch_dict["M0"], segment_index=0, **batch_dict)
+        
         # Loop over k-space segments, defined by k-space trajectory in self.encoding_module
         m_transverse = self.forward_model(batch_dict["M0"], segment_index=0, **batch_dict)
         for segment_index in tf.range(n_segments):
diff --git a/cmrsim/bloch/__init__.py b/cmrsim/bloch/__init__.py
index 7b4c1ad5639b6a04b0318fe4e2db4176f5bcd540..949228d9fa806dda8cffee6aca717c8062532615 100644
--- a/cmrsim/bloch/__init__.py
+++ b/cmrsim/bloch/__init__.py
@@ -3,10 +3,11 @@
 from cmrsim.bloch._base import *
 from cmrsim.bloch._generic import *
 from cmrsim.bloch._ideal import *
+from cmrsim.bloch._multi_coil import *
 import cmrsim.bloch.submodules
 
 import importlib
-_submodules = ["_base", "_generic", "_ideal", "_ideal"]
+_submodules = ["_base", "_generic", "_ideal", "_multi_coil"]
 mod_handles = [importlib.import_module(f"cmrsim.bloch.{m}") for m in _submodules]
 __all__ = [item for m in mod_handles for item in getattr(m, '__all__')] + ["submodules"]
 
diff --git a/cmrsim/bloch/_generic.py b/cmrsim/bloch/_generic.py
index e7ff204d9b25c8ef8cf1ce894279456ae3b95481..6e7494560b7b67bd5c331a51fdb5d3707c8b45cb 100644
--- a/cmrsim/bloch/_generic.py
+++ b/cmrsim/bloch/_generic.py
@@ -328,118 +328,158 @@ class GeneralBlochOperator(BaseBlochOperator):
         :param n_repetitions: actual n_reps for par, 1 for sequential
         :param n_samples:
         :param trajectory_module:
-        :param initial_position: (n_repetitions/1, #batch, 3)
+        :param initial_position: (#batch, 3)
         :param magnetization: (n_repetitions/1, #batch, 3)
         :param T1: (n_repetitions/1, #batch)
         :param T2: (n_repetitions/1, #batch)
         :param M0: (n_repetitions/1, #batch)
         :param kwargs:
-        :return: magnetization, r_next, signal at adc_on==True
+        :return:  - magnetization
+                  - r_next
+                  - signal at adc_on==True with shape (n_repetitions, n_samples)
+
         """
         with tf.device(self.device):
-            batch_signal = tf.TensorArray(dtype=tf.complex64, size=n_samples + tf.constant(1),
-                                          element_shape=(n_repetitions,))
-            result_writing_idx = 0
-            r_next = initial_position
-            repetitions = tf.range(repetition_index, repetition_index + n_repetitions, 1)
-
-            # If gradient waveforms are specified, gather all required repetitions (either all or 1)
-            # and calculate the left value of for trapezoidal phase integration
-            if self.gradient_waveforms is not None:
-                current_gradient_waveforms = tf.gather(self.gradient_waveforms, repetitions,
-                                                       axis=0)[:, tf.newaxis]
-                gdotr_l = tf.reduce_sum(current_gradient_waveforms[..., 0, :] * r_next, axis=-1)
-            else:
-                gdotr_l = tf.constant(0, dtype=tf.float32)
-            gdotr_r = tf.zeros_like(gdotr_l)
-
-            # If rf waveforms are specified, gather all required repetitions (either all or 1)
-            if self.rf_waveforms is not None:
-                current_rf_waveforms = tf.gather(self.rf_waveforms, repetitions, axis=0)
-
-            # If adc events are specified, gather all required repetitions (either all or 1)
-            # or otherwise create a default zero return value for the batch signal as this it is
-            # required autographing this function by tensorflow
-            if self.adc_events is None:
-                for i in tf.range(n_samples):
-                    batch_signal = batch_signal.write(i, tf.zeros(n_repetitions,
-                                                                  dtype=tf.complex64))
-            else:
-                current_adc_events = tf.gather(self.adc_events, repetitions, axis=0)
-                current_adc_phase = tf.gather(self.adc_phase, repetitions, axis=0)
+            rf_waveforms, grad_waveforms, adc_on, adc_phase, adc_widx, batch_signal = \
+                self._init_tensors(n_samples, n_repetitions, repetition_index)
+
+            r_prev = tf.reshape(initial_position, [1, -1, 3])
+            r_next = r_prev
 
             # Iterate over all time intervals to obtain the magnetization evolution
-            for t_idx in tf.range(0, tf.shape(self.dt_grid)[0]):
+            for t_idx in tf.range(1, tf.shape(self.dt_grid)[0]+1):
 
                 # The shape of the magnetization tensor must always be (#reps, #batch, 3), with
                 # needs to be specified for tf.autograph as:
                 tf.autograph.experimental.set_loop_options(
                     shape_invariants=[(magnetization, tf.TensorShape([None, None, 3])),
-                                      (r_next, tf.TensorShape([None, None, 3])),
-                                      (gdotr_r, tf.TensorShape([None, None])),
-                                      (gdotr_l, tf.TensorShape([None, None]))
-                                      ])
+                                      (r_next, tf.TensorShape([None, None, 3]))])
 
                 # Increment particle positions and if available field definitions at given locations
                 # r_next.shape = (#particles, #nsteps=1)
-                r_next, additional_fields = trajectory_module.increment_particles(
-                                                                    r_next[0], self.dt_grid[t_idx])
+                r_next, additional_fields = trajectory_module.increment_particles(r_prev[0], self.dt_grid[t_idx-1])
                 r_next = r_next[tf.newaxis]   # add dummy #reps-dim for subsequent broadcasting
 
                 # Evaluate effective flip angle for given interval and apply the hard-pulse rotation
                 if self.rf_waveforms is not None:
-                    rf = current_rf_waveforms[:, t_idx:t_idx + 2]
-                    alpha = (tf.complex(self.dt_grid[t_idx] * self.gamma, 0.)
+                    rf = rf_waveforms[:, t_idx-1:t_idx+1]
+                    alpha = (tf.complex(self.dt_grid[t_idx-1] * self.gamma, 0.)
                              * (rf[:, 1] + rf[:, 0]) / 2)
                     magnetization = self.hard_pulse(alpha[:, tf.newaxis], magnetization)
 
                 # Evaluate effective phase increment for given interval including the submodule
                 # effects and apply the corresponding rotation of the magnetization vectors
+                delta_phi = tf.zeros_like(M0, dtype=tf.float32)
                 if self.gradient_waveforms is not None:
-                    gdotr_r = tf.reduce_sum(current_gradient_waveforms[..., t_idx + 1, :] * r_next,
-                                            axis=-1)
-                    delta_phi = (gdotr_l + gdotr_r) / 2. * self.gamma * self.dt_grid[t_idx]
-                    gdotr_l = gdotr_r
-                    for submod in self._submodules:
-                        phi_s = submod(
-                            gradient_wave_form=current_gradient_waveforms[:, 0, t_idx:t_idx + 1, :],
-                            trajectories=r_next[:, :, tf.newaxis], dt=self.dt_grid[t_idx],
-                            **additional_fields, **kwargs)
-                        if tf.shape(phi_s)[0] < tf.shape(delta_phi)[0]:
-                            phi_s = tf.tile(phi_s, [tf.shape(delta_phi)[0], 1, 1])
-                        # phi_s shape = (#reps, #particles, #nsteps=1)
-                        tf.ensure_shape(phi_s, (None, None, 1))
-                        delta_phi = tf.reshape(delta_phi + phi_s[:, :, 0], [n_repetitions, -1])
-
-                    # delta phi for current interval with (nreps, nbatch) dimension
-                    tf.ensure_shape(delta_phi, (None, None))
-                    magnetization = self.phase_increment(delta_phi[:, tf.newaxis], magnetization)
+                    grad = grad_waveforms[..., t_idx-1:t_idx+1, :]
+                    gdotr_l = tf.reduce_sum(grad[..., 0, :] * r_prev, axis=-1)
+                    gdotr_r = tf.reduce_sum(grad[..., 1, :] * r_next, axis=-1)
+                    delta_phi += (gdotr_l + gdotr_r) / 2. * self.gamma * self.dt_grid[t_idx-1]
+
+                for submod in self._submodules:
+                    phi_s = submod(
+                        gradient_wave_form=grad_waveforms[:, 0, t_idx:t_idx+1, :],
+                        trajectories=r_next[:, :, tf.newaxis], dt=self.dt_grid[t_idx-1],
+                        **additional_fields, **kwargs)
+                    if tf.shape(phi_s)[0] < tf.shape(delta_phi)[0]:
+                        phi_s = tf.tile(phi_s, [tf.shape(delta_phi)[0], 1, 1])
+                    # phi_s shape = (#reps, #particles, #nsteps=1)
+                    tf.ensure_shape(phi_s, (None, None, 1))
+                    delta_phi = tf.reshape(delta_phi + phi_s[:, :, 0], [n_repetitions, -1])
+
+                # delta phi for current interval with (nreps, nbatch) dimension
+                tf.ensure_shape(delta_phi, (None, None))
+                magnetization = self.phase_increment(delta_phi[:, tf.newaxis], magnetization)
 
                 # Relax magnetization for the time interval
-                magnetization = self.relax(T1, T2, self.dt_grid[t_idx], magnetization)
+                magnetization = self.relax(T1, T2, self.dt_grid[t_idx-1], magnetization)
 
                 # If adc-events are specified and the current interval is set to 'on' calculate
                 # the weighed sum over all particles and demodulate according to the current ADC
                 # phase. The result is then written into the signal return array.
-                # Note: To make this work with tf.function, the conditional needs an 'else' that
-                # also writes into the batch_signal return array at the last position. Therefore,
-                # the last entry is cut-off in the end, which implicitly assumes the last interval
-                # is never an active ADC-sample
                 if self.adc_events is not None:
-                    n_active_adc_channels = tf.reduce_sum(current_adc_events[:, t_idx + 1], axis=0)
-                    if n_active_adc_channels >= 1.:
-                        signal = tf.reduce_sum(magnetization[..., 0] * tf.complex(M0, 0.), axis=1)
-                        signal = signal * tf.complex(current_adc_events[:, t_idx + 1], 0.)
-                        demodulated_signal = signal * tf.exp(
-                            tf.complex(0., -current_adc_phase[:, t_idx + 1]))
-                        batch_signal = batch_signal.write(result_writing_idx, demodulated_signal)
-                        result_writing_idx += 1
-                    else:
-                        batch_signal = batch_signal.write(result_writing_idx,
-                                                          tf.zeros(n_repetitions,
-                                                                   dtype=tf.complex64))
-            batch_signal = batch_signal.stack()
-            return magnetization, r_next, tf.transpose(batch_signal[:-1], [1, 0])
+                    batch_signal = self._sample(adc_on[:, t_idx], adc_phase[:, t_idx],
+                                                adc_widx[:, t_idx], magnetization, M0,
+                                                r_next, batch_signal)
+                    
+                r_prev = r_next
+                
+            transpose_indices = tf.range(tf.shape(tf.shape(batch_signal))[0])
+            transpose_indices = tf.tensor_scatter_nd_update(transpose_indices, [[0], [1]], [1, 0])
+            return magnetization, r_next, tf.transpose(batch_signal, transpose_indices)
+
+    def _init_tensors(self, n_samples: int, n_repetitions: int, repetition_index: int):
+        """ Initializes / slices the waveforms according to the specified current repetitions and
+        number of samples.
+
+        :param n_samples:
+        :param n_repetitions:
+        :param repetition_index:
+        :return: - rf_waveform (#reps, #steps) or None
+                 - grad_waveform (#reps, #steps, 3) or None
+                 - adc_events (#reps, #steps) or None
+                 - adc_phase (#reps, #steps) or None
+                 - adc_writing_idx (#reps, #steps) or None
+                 - batch_signal (#reps, #samples) -> used for signal accumulation
+        """
+        batch_signal = tf.zeros(shape=(n_samples, n_repetitions), dtype=tf.complex64)
+        repetitions = tf.range(repetition_index, repetition_index + n_repetitions, 1)
+
+        # If rf waveforms are specified, gather all required repetitions (either all or 1)
+        if self.rf_waveforms is not None:
+            current_rf_waveforms = tf.gather(self.rf_waveforms, repetitions, axis=0)
+        else:
+            current_rf_waveforms = None
+
+        # If gradient waveforms are specified, gather all required repetitions (either all or 1)
+        # and calculate the left value of for trapezoidal phase integration
+        if self.gradient_waveforms is not None:
+            current_gradient_waveforms = tf.gather(self.gradient_waveforms, repetitions,
+                                                   axis=0)[:, tf.newaxis]
+        else:
+            current_gradient_waveforms = None
+
+        # If adc events are specified, gather all required repetitions (either all or 1)
+        # or otherwise create a default zero return value for the batch signal as this it is
+        # required autographing this function by tensorflow
+        if self.adc_events is not None:
+            current_adc_events = tf.gather(self.adc_events, repetitions, axis=0)
+            current_adc_phase = tf.gather(self.adc_phase, repetitions, axis=0)
+            adc_writing_indices = tf.cast(tf.math.cumsum(current_adc_events, axis=1) - 1, tf.int32)
+        else:
+            current_adc_events, current_adc_phase, adc_writing_indices = None, None, None
+
+        return (current_rf_waveforms, current_gradient_waveforms, current_adc_events,
+                current_adc_phase, adc_writing_indices, batch_signal)
+
+    @staticmethod
+    def _sample(adc_events: tf.Tensor, adc_phase: tf.Tensor, adc_writing_idx: tf.Tensor,
+                magnetization: tf.Tensor, M0: tf.Tensor, r: tf.Tensor, batch_signal: tf.Tensor):
+        """Sums up the :math:`m_{xy}^{+}` magnetization weighted with the associated proton-density
+        of the particles and adds the signal to the batch-signal accumulator ad index specified
+        in adc_writing_index.
+
+        :param adc_events: (#reps, )
+        :param adc_phase: (#reps, )
+        :param adc_writing_idx: (#reps, )
+        :param magnetization: (#batch, )
+        :param M0: (n_repetitions/1, #batch)
+        :param r: (n_repetitions/1, #batch, 3)
+        :param batch_signal: (#samples, #reps)
+        :return:
+        """
+        n_active_adc_channels = tf.reduce_sum(adc_events, axis=0)
+        if n_active_adc_channels >= 1.:
+            signal = tf.reduce_sum(magnetization[..., 0] * tf.complex(M0, 0.), axis=1)
+            signal = signal * tf.complex(adc_events, 0.)
+            demodulated_signal = signal * tf.exp(tf.complex(0., -adc_phase))
+            n_reps = tf.shape(adc_writing_idx)[0]
+            indices = tf.stack([adc_writing_idx, tf.range(n_reps, dtype=tf.int32)], axis=1)
+            in_bounds = tf.where(adc_writing_idx > -1)
+            indices = tf.gather(indices, in_bounds)[0]
+            demodulated_signal = tf.gather(demodulated_signal, in_bounds)[0]
+            batch_signal = tf.tensor_scatter_nd_add(batch_signal, indices, demodulated_signal)
+        return batch_signal
 
     def reset(self):
         """ Resets the time-signal accumulators, i.e. sets all entries to 0"""
diff --git a/cmrsim/bloch/_multi_coil.py b/cmrsim/bloch/_multi_coil.py
new file mode 100644
index 0000000000000000000000000000000000000000..8157cbe288c758e890fb2cccc74c6489f6b36398
--- /dev/null
+++ b/cmrsim/bloch/_multi_coil.py
@@ -0,0 +1,121 @@
+"""Contains implementation of Bloch-solving operator including multiple receive coils"""
+
+__all__ = ["ParallelReceiveBlochOperator"]
+
+from typing import List, Tuple, Union, TYPE_CHECKING
+
+import tensorflow as tf
+import numpy as np
+
+from cmrsim.bloch._generic import GeneralBlochOperator
+from cmrsim.bloch.submodules import BaseSubmodule
+from cmrsim.trajectory import StaticParticlesTrajectory
+
+if TYPE_CHECKING:
+    import cmrsim.analytic.contrast.coil_sensitivities
+
+
+class ParallelReceiveBlochOperator(GeneralBlochOperator):
+    """Extends the GeneralBlochOperator with parallel receive functionality. The only additionally
+    required input is an instance of a CoilSensitivity contrast module from the
+    cmrsim.analytic.contrast submodule, which is used internally to apply a spatially dependent
+    weight to the transverse magnetization before summing, in the process of acquiring signals.
+
+    The resulting signal is also stored in the 'time_signal_acc' attribute, where for the
+    parallel receive case each entry has an additional axis corresponding to the number of coils.
+
+    :param name: Verbose name of the module (non-functionally relevant)
+    :param gamma: gyromagnetic ratio in rad/ms/m
+    :param coil_module:
+    :param time_grid: (n_steps) - tf.float32 - time definition for numerical integration steps
+    :param gradient_waveforms: (#repetitions, #steps, 3) - tf.float32 - Definition of gradients
+                        in mT/m
+    :param rf_waveforms: (#repetitions, #steps) - tf.complex64 -
+    :param adc_events: (#repetition, #steps, 2) - tf.float32 - last axis contains the adc-on
+                            definitions and the adc-phase
+    :param submodules: List[tf.Modules] implementing the __call__ function returning the
+                            resulting phase increment contribution
+    :param device: str
+    """
+    #: Instance of the cmrsim.analytic.contrast.CoilSensitivity class that is used to compute
+    #: the positionally dependent sensitivity factors prior to sampling
+    coil_lookup_module: 'cmrsim.analytic.contrast.coil_sensitivities.CoilSensitivity'
+
+    #:
+    n_coils: int
+
+    def __init__(self, name: str,
+                 gamma: float,
+                 coil_module: Union['cmrsim.analytic.contrast.coil_sensitivities.CoilSensitivity'],
+                 time_grid: tf.Tensor,
+                 adc_events: tf.Tensor,
+                 gradient_waveforms: tf.Tensor = None,
+                 rf_waveforms: tf.Tensor = None,
+                 submodules: Tuple[BaseSubmodule, ...] = (),
+                 device: str = None):
+        """
+
+
+        """
+        super().__init__(name=name, gamma=gamma, time_grid=time_grid,
+                         gradient_waveforms=gradient_waveforms, rf_waveforms=rf_waveforms,
+                         adc_events=adc_events, submodules=submodules, device=device)
+        self.coil_lookup_module = coil_module
+        self.n_coils = self.coil_lookup_module.expansion_factor
+
+        # Expand the dimension of the signal accumulator to incorporate coil-maps
+        with tf.device("CPU"):
+            self.time_signal_acc = [
+                tf.Variable(initial_value=tf.zeros([v.read_value().shape[0], self.n_coils], dtype=tf.complex64),
+                            shape=(None, self.n_coils), dtype=tf.complex64, trainable=False,
+                            name=f"k_segment_{idx}") for idx, v in enumerate(self.time_signal_acc)]
+
+    def _init_tensors(self, n_samples: int, n_repetitions: int, repetition_index: int):
+        """ Initializes / slices the waveforms according to the specified current repetitions and
+        number of samples.
+
+        :param n_samples:
+        :param n_repetitions:
+        :param repetition_index:
+        :return: - rf_waveform (#reps, #steps) or None
+                 - grad_waveform (#reps, #steps, 3) or None
+                 - adc_events (#reps, #steps) or None
+                 - adc_phase (#reps, #steps) or None
+                 - adc_writing_idx (#reps, #steps) or None
+                 - batch_signal (#reps, #samples) -> used for signal accumulation
+        """
+        (rf_waveforms, gradient_waveforms, adc_events, adc_phase, adc_writing_indices,
+         batch_signal) = super()._init_tensors(n_samples, n_repetitions, repetition_index)
+        batch_signal = tf.zeros(shape=(n_samples, n_repetitions, self.n_coils), dtype=tf.complex64)
+        return (rf_waveforms, gradient_waveforms, adc_events, adc_phase,
+                adc_writing_indices, batch_signal)
+
+    def _sample(self, adc_events: tf.Tensor, adc_phase: tf.Tensor, adc_writing_idx: tf.Tensor,
+                magnetization: tf.Tensor, M0: tf.Tensor, r: tf.Tensor, batch_signal: tf.Tensor):
+        """Sums up the :math:`m_{xy}^{+}` magnetization weighted with the associated proton-density
+        of the particles and adds the signal to the batch-signal accumulator ad index specified
+        in adc_writing_index. Before summation, a coil-sensitivity weighting is applied
+
+        :param adc_events: (#reps, )
+        :param adc_phase: (#reps, )
+        :param adc_writing_idx: (#reps, )
+        :param magnetization: (#batch, )
+        :param M0: (n_repetitions/1, #batch)
+        :param r: (n_repetitions/1, #batch, 3) particle positions
+        :param batch_signal: (#samples, #reps, #coils)
+        :return:
+        """
+        n_active_adc_channels = tf.reduce_sum(adc_events, axis=0)
+        if n_active_adc_channels >= 1.:
+            sense_factors = self.coil_lookup_module.lookup_complex_coil_factors(r)
+            weighted_magnetization = magnetization[..., 0, tf.newaxis] * sense_factors
+            signal = tf.reduce_sum(weighted_magnetization * tf.complex(M0[..., tf.newaxis], 0.), axis=1)
+            signal = signal * tf.complex(adc_events[..., tf.newaxis], 0.)
+            demodulated_signal = signal * tf.exp(tf.complex(0., -adc_phase[..., tf.newaxis]))
+            n_reps = tf.shape(adc_writing_idx)[0]
+            indices = tf.stack([adc_writing_idx, tf.range(n_reps, dtype=tf.int32)], axis=1)
+            in_bounds = tf.where(adc_writing_idx > -1)
+            indices = tf.gather(indices, in_bounds)[0]
+            demodulated_signal = tf.gather(demodulated_signal, in_bounds)[0]
+            batch_signal = tf.tensor_scatter_nd_add(batch_signal, indices, demodulated_signal)
+        return batch_signal
diff --git a/cmrsim/bloch/submodules.py b/cmrsim/bloch/submodules.py
index 9df756633425ba0145e2516f602376bb3637b47c..bc656184b66416f9436c4bddd4b14906e3d53fce 100644
--- a/cmrsim/bloch/submodules.py
+++ b/cmrsim/bloch/submodules.py
@@ -122,14 +122,14 @@ class ConcomitantFields(BaseSubmodule):
         """
         with tf.device(self.device):
             tf.Assert(tf.shape(tf.shape(gradient_wave_form))[0] == 3,
-                      ["Shape missmatch in ConcomittantFields for specified gradient",
+                      ["Shape missmatch in ConcomitantFields for specified gradient",
                        tf.shape(gradient_wave_form)])
             gradient_wave_form = tf.expand_dims(gradient_wave_form, axis=1)
 
             if tf.shape(tf.shape(trajectories))[0] == 3:
                 trajectories = tf.expand_dims(trajectories, axis=0)
             tf.Assert(tf.shape(tf.shape(gradient_wave_form))[0] == 4,
-                      ["Shape missmatch in ConcomittantFields for trajectories",
+                      ["Shape missmatch in ConcomitantFields for trajectories",
                       tf.shape(trajectories)])
 
             max_shape = tf.reduce_max(tf.stack([tf.shape(gradient_wave_form),
diff --git a/cmrsim/datasets/__init__.py b/cmrsim/datasets/__init__.py
index 549dc1007dc69cfb84b89a43fcab33c0d1247ccd..b6b1c70a49a585466a6d7150a2a3c4e8e1c69c51 100644
--- a/cmrsim/datasets/__init__.py
+++ b/cmrsim/datasets/__init__.py
@@ -10,4 +10,5 @@ from cmrsim.datasets._analytic import AnalyticDataset
 from cmrsim.datasets._flow import RefillingFlowDataset
 from cmrsim.datasets._bloch import BlochDataset
 from cmrsim.datasets._cardiac_mesh import MeshDataset, CardiacMeshDataset
+from cmrsim.datasets._regular_grid import RegularGridDataset
 
diff --git a/cmrsim/datasets/_flow.py b/cmrsim/datasets/_flow.py
index faf8e9e46c75a4cdf315ff1adde321f0828ac8db..9526c19b53c77460db43dfaaa3a62709206fb1d0 100644
--- a/cmrsim/datasets/_flow.py
+++ b/cmrsim/datasets/_flow.py
@@ -12,6 +12,7 @@ import pyvista
 
 from scipy.ndimage import gaussian_filter
 
+
 # pylint: disable=too-many-instance-attributes
 class RefillingFlowDataset(tf.Module):
     """ Wraps a vtk-mesh of containing a velocity field. Functionality tackles the problem of in-/
diff --git a/cmrsim/datasets/_regular_grid.py b/cmrsim/datasets/_regular_grid.py
new file mode 100644
index 0000000000000000000000000000000000000000..5b9c185ce10bd3bdc346475ee532efd14d000b32
--- /dev/null
+++ b/cmrsim/datasets/_regular_grid.py
@@ -0,0 +1,299 @@
+__all__ = ["RegularGridDataset", ]
+
+from typing import Dict, Tuple, Union, Sequence
+
+import tensorflow as tf
+import numpy as np
+import pyvista
+from pint import Quantity
+
+
+class RegularGridDataset:
+    """Dataset class containing functionality, that requires a regular-grid representation of the
+    data. Multiple class-methods define inputs from which an instance of this class can be
+    created.
+
+    :param mesh: Instance of a pyvista.UniformGrid which is copied on instantiation
+    """
+    #: Instance of a pyvista.UniformGrid used as basis for all contained computations
+    mesh: pyvista.UniformGrid
+
+    def __init__(self, mesh: pyvista.UniformGrid):
+        self.mesh = mesh.copy()
+
+    def get_phantom_def(self, filter_by: str = None, keys: Sequence[str] = None,
+                        perturb_positions_std: float = None) -> Dict[str, np.ndarray]:
+        """Extracts points and their corresponding physical properties from the contained mesh
+        and returns them as dictionary of flattened array
+
+        :param filter_by: used to determine which points to extract from the mesh by using the
+                            indices = np.where(self.mesh[filter_by])
+        :param keys: sequence of strings defining the field arrays returned as dictionary keys
+        :param perturb_positions_std: if specified the positions of the extracted mesh-points are
+                    perturbed by sampling a normal distribution with zero mean and given standard
+                    deviation to reduce descretization artifacts in simulation
+        :return:
+        """
+
+        if filter_by is not None:
+            indices = np.where(self.mesh[filter_by])
+        else:
+            indices = np.where(np.ones_like(self.mesh.points[:, 0]))
+
+        if keys is None:
+            keys = self.mesh.array_names
+
+        r_vectors = np.array(self.mesh.points[indices], dtype=np.float32).reshape(-1, 3)
+        if perturb_positions_std is not None:
+            r_vectors += np.random.normal(0, perturb_positions_std, size=r_vectors.shape)
+
+        phantom = dict(r_vectors=r_vectors)
+        phantom.update({k: np.squeeze(self.mesh[k][indices].reshape(r_vectors.shape[0], -1))
+                        for k in keys})
+        return phantom
+
+    def add_field(self, key: str, data: np.ndarray):
+        """Adds field to self.mesh from a numpy array (takes care of correct ordering from
+        default numpy C-order to F-order used in pyvista).
+
+        :param key: name of the new field
+        :param data: (x, y, z, N) dimensional array containing arbitrary data
+        """
+        n_channels = np.prod(data.shape) // np.prod(self.mesh.dimensions)
+        flat_data = data.reshape([-1, n_channels], order="F")
+        self.mesh[key] = flat_data
+
+    def resample_spacing(self, spacing: Quantity):
+        """Resamples the contained self.mesh to a new spacing, which might be necessary to reduce
+        discretization artifacts on simulation.
+
+        :param spacing: New target spacing (3, )
+        """
+        r_max = np.max(self.mesh.points, axis=0)
+        r_min = np.min(self.mesh.points, axis=0)
+        nbins = (Quantity(r_max - r_min, "m") / spacing.to("m")).m.astype(int)
+        new_mesh = pyvista.UniformGrid(dimensions=nbins, spacing=spacing.m_as("m"), origin=r_min)
+        self.mesh = new_mesh.sample(self.mesh)
+
+    def select_slice(self, slice_normal: np.array, spacing: Quantity,
+                     slice_position: Quantity, field_of_view: Quantity,
+                     readout_direction: np.array = None,
+                     in_mps: bool = False) -> pyvista.UniformGrid:
+        """Creates a uniform grid (slice) according to the specified parameters, and uses it to
+        probe the self.mesh (see example). Slice (coordinates) can either be returned in mps
+        or in global coordinates.
+
+        .. Dropdown:: Example Visualization
+
+            .. image:: ../../_static/api/dataset_regular_grid_select_slice.png
+
+        :param slice_normal: (3, ) normal vector of the extracted slice
+        :param spacing: (3, ) dimension of voxels in all directions
+        :param slice_position: (3, ) vector pointing to the slice-center
+        :param field_of_view: (3, ) total extent of the extracted slice in 3 dimensions
+        :param readout_direction: (3, ) if not provided, is assumed to be x (part orthogonal to
+                specified slice normal). Used to rotate the extracted slice about the slice-normal
+        :param in_mps: if true the slice is returned in MPS coordinates
+        :return:
+        """
+        from scipy.spatial.transform import Rotation
+        slice_normal /= np.linalg.norm(slice_normal)
+
+        if readout_direction is None:
+            if np.allclose(slice_normal, np.array([1., 0., 0.]), atol=1e-3):
+                readout_direction = np.array([0., 1., 0.])
+            else:
+                readout_direction = np.array([1., 0., 0.])
+        else:
+
+            if np.isclose(
+                    np.dot(slice_normal, readout_direction / np.linalg.norm(readout_direction)),
+                    1., atol=1e-3):
+                raise ValueError("Slice normal and readout-direction are parallel")
+
+        readout_direction = readout_direction - np.dot(slice_normal,
+                                                       readout_direction) * slice_normal
+        readout_direction /= np.linalg.norm(readout_direction)
+        pe_direction = np.cross(readout_direction, slice_normal)
+        pe_direction /= np.linalg.norm(pe_direction)
+        original_basis = np.eye(3, 3)
+        new_basis = np.stack([readout_direction, pe_direction, slice_normal])
+
+        rot, _ = Rotation.align_vectors(original_basis, new_basis)
+
+        total_transform = np.eye(4, 4)
+        total_transform[:3, :3] = rot.as_matrix()
+        total_transform[:3, 3] = slice_position.m_as("m")
+        nbins = (field_of_view / spacing).m_as("dimensionless").astype(int)
+        new_slice = pyvista.UniformGrid(dimensions=nbins, spacing=spacing.m_as("m"),
+                                        origin=-field_of_view.m_as("m") / 2)
+        new_slice = new_slice.transform(total_transform, inplace=False)
+        new_slice = self.mesh.probe(new_slice)
+
+        if in_mps:
+            inv_transform = np.eye(4, 4)
+            inv_transform[:3, :3] = rot.as_matrix().T
+            inv_transform[:3, 3] = -slice_position.m_as("m")
+            new_slice = new_slice.transform(inv_transform, inplace=False,
+                                            transform_all_input_vectors=True)
+        return new_slice
+
+    def compute_offresonance(self, b0: Quantity, susceptibility_key: str) -> Quantity:
+        """Approximates the magnetic field variations in z-direction due to interfaces with varying
+        susceptibilities by applying a convolution (Fourier Product) with a dipole kernel to the
+        susceptibility field contained in self.mesh. The name of the mesh array is specified by
+        argument.
+
+        Assumes that the main magnetic field is pointing in z-direction of the contained mesh.
+        The susceptibility os assumed to be specified as parts per million
+        (e.g. :math:`\chi_{air} = 0.36ppm`).
+
+        The result is store in self.mesh["offres"] and is also returned as Quantity.
+
+        .. Dropdown:: Literature Reference
+
+            This function is based on the work of Frank-Zijlstra, published on Mathworks:
+
+            Frank Zijlstra (2022). MRI simulation using FORECAST: Fourier-based Off-REsonanCe
+            Artifact simulation in the STeady-state
+            (https://www.mathworks.com/matlabcentral/fileexchange/56680-mri-simulation-using-
+            forecast-fourier-based-off-resonance-artifact-simulation-in-the-steady-state),
+
+            MATLAB Central File Exchange. Retrieved December 8, 2022.
+
+        :param b0: Field strength of main magnetic field pointing in z-direction in Tesla
+        :param susceptibility_key: name of the array contained in self.mesh, used as susceptibility
+        :return: (X, Y, Z) array wrapped as Quantity with dimension ("T")
+        """
+        chi_3d = self.mesh[susceptibility_key].reshape(self.mesh.dimensions, order="F")
+        # Pad the chi-map to a multiple of four:
+        residual_per_dim = np.mod(self.mesh.dimensions, 4)
+        chi_3d = np.pad(chi_3d,
+                        np.stack([np.zeros_like(residual_per_dim), residual_per_dim], axis=1),
+                        mode="constant", constant_values=0.)
+        voxel_size = Quantity(self.mesh.spacing, "m")
+        field_of_view = voxel_size * Quantity(chi_3d.shape, "dimensionless")
+
+        # The fourier based convolution (from original implementation):
+        # 1) The forward Fourier Transform (FT): Instead of doing an extra convolution to calculate
+        # the effect of aliasing (as proposed in the article), here we use the 'imaginary channel'
+        # to do this. This is done by downscaling the artificial environment, and adding this as
+        # an imaginary component to the original (unpadded) susceptibility distribution.
+        # Together this is referred to as the "DUAL" distribution:
+        # down_sampled_chi = tf.nn.avg_pool3d(chi_3d[np.newaxis, ..., np.newaxis], ksize=(2, 2, 2),
+        #                                     strides=(2, 2, 2), padding="VALID")
+        # imaginary_channel = self._zc(down_sampled_chi[0, ..., 0])
+        # ft_chi_3d_dual = tf.signal.fft3d(chi_3d + 1j * imaginary_channel)
+        ft_chi_3d_dual = tf.signal.fft3d(chi_3d)
+        # 2) multiplication with the dipole function
+        dipole_kernel = self._compute_fourier_dipole_kernel(field_of_view,
+                                                       np.array(ft_chi_3d_dual.shape))
+        ft_field_3d = ft_chi_3d_dual * dipole_kernel
+        # 3) inverse FT to spatial domain
+        field_3d_dual = tf.signal.ifft3d(ft_field_3d)
+        # 4) Subtract the up-scaled cropped center of the aliasing
+        # field_3d = tf.math.real(field_3d_dual) - UC(tf.math.imag(field_3d_dual))
+        field_3d = tf.math.real(field_3d_dual)
+        i, j, k = [_ if _ != 0 else None for _ in -residual_per_dim]
+        field_3d = field_3d.numpy()[:i, :j, :k]
+        self.mesh["offres"] = b0.to("mT") * field_3d.reshape(-1, order="f")
+        return field_3d * b0.to("mT")
+
+    @staticmethod
+    def _compute_fourier_dipole_kernel(fov: Quantity, n_samples: np.ndarray):
+        """ Computes the dipole kernel for a given FOV and 3D array
+        :param fov: (3, ) length - quantity
+        :param n_samples: (3, ) integer number of regular points per dimension
+        :return:
+        """
+        k_squared = [np.arange(-n/2, n/2, 1) / fov[i].m_as("m") for i, n in enumerate(n_samples)]
+        k_squared = [np.fft.ifftshift(k)**2 for k in k_squared]
+        [kx2, ky2, kz2] = np.meshgrid(*k_squared, indexing="ij")
+        kernel = 1/3 - kz2 / (kx2 + ky2 + kz2)
+        kernel[0, 0, 0] = 0
+        return kernel
+
+    @classmethod
+    def from_unstructured_grid(cls, input_mesh: pyvista.DataSet, pixel_spacing: Quantity,
+                               padding: Quantity = None) -> 'RegularGridDataset':
+        """Resamples a given unstructured mesh onto a regular grid and instantiates a
+        RegularGridDataset from it.
+
+        :param input_mesh: pyvista.Dataset
+        :param pixel_spacing: Quantity (3,) specifying the 3D spacing of the uniform grid
+        :param padding: Space added symmetrically around the bounding box of the input_mesh,
+                        defaults to (0, 0, 0) "meter"
+        """
+        if padding is None:
+            padding = Quantity([0., 0., 0.], "m")
+        r_max = np.max(input_mesh.points, axis=0) + padding.m_as("m") / 2
+        r_min = np.min(input_mesh.points, axis=0) - padding.m_as("m") / 2
+        nbins = (Quantity(r_max - r_min, "m") / pixel_spacing.to("m")).m.astype(int)
+        mesh = pyvista.UniformGrid(dimensions=nbins, spacing=pixel_spacing.m_as("m"), origin=r_min)
+        mesh = input_mesh.probe(mesh)
+        return cls(mesh)
+
+    @classmethod
+    def from_numpy(cls, data: Union[np.ndarray, Dict[str, np.ndarray]],
+                   origin_offset: Quantity = Quantity([0, 0, 0], "m"), extent: Quantity = None,
+                   pixel_spacing: Quantity = None) -> 'RegularGridDataset':
+        """Constructs a uniform regular 3D grid of the same shape as the data provided as argument.
+
+        Assumes that the origin is at the center of the provided array.
+
+        :raises: ValueError:
+                    - if neither extent nor pixel_spacing is specified
+                    - data-argument is a dictionary containing arrays with non-identical shapes
+
+        :param data: (X, Y, Z, ...)
+        :param origin_offset: (3, ) If not provided, the grid is positioned symmetrically around the
+                                point (0., 0., 0.). Otherwise this is vector is added onto the
+                                center-point.
+        :param extent: (3, ) Length-Quantity defining the extend of the provided data in 3 dimen
+        :param pixel_spacing: (3, )
+        """
+
+        if (all(v is None for v in [pixel_spacing, extent]) or
+                all(v is not None for v in [pixel_spacing, extent])):
+            raise ValueError("Either 'extent'/'pixel-spacing' keyword argument must be specified")
+
+        if isinstance(data, dict):
+            shapes = [v.shape for v in data.values()]
+            if not all(shapes[0][0:3] == x[0:3] for x in shapes):
+                raise ValueError("Data arrays specified in argument 'data' are not identical:\n"
+                                 + "\n".join([str(s) for s in shapes]))
+        else:
+            shapes = [data.shape, ]
+
+        if extent is not None:
+            pixel_spacing = extent.to("m") / np.array(shapes[0])
+
+        shape = np.array(shapes[0][0:3])
+        origin = - pixel_spacing * shape / 2 + origin_offset
+
+        mesh = pyvista.UniformGrid(dimensions=shape, spacing=pixel_spacing.m_as("m"),
+                                   origin=origin.m_as("m"))
+
+        for name, arr in data.items():
+            mesh[name] = arr.reshape([shape.prod(), *arr.shape[3:]], order="F")
+
+        return cls(mesh)
+
+    @classmethod
+    def from_shepp_logan(cls, map_size: Tuple[int, int, int], extent: Quantity) \
+            -> 'RegularGridDataset':
+        """Creates a 3D Shepp-Logan phantom and stores it in a RegularGridDataset.
+
+        :param map_size: Target number of voxels per dimension
+        :param extent: Target spatial extent of the resulting shepp-logan phantom.
+        """
+        try:
+            import phantominator
+        except ImportError:
+            raise ImportError("The python-package phantominator is necessary to create a shepp-"
+                              "logan phantom. To install it, type:\n pip install phantominator")
+
+        M0, T1, T2 = phantominator.shepp_logan(map_size, MR=True)
+        M0, T1, T2 = [t.transpose(2, 1, 0) for t in [M0, T1, T2]]
+        data = dict(rho=M0, T1=T1 * 1000, T2=T2*1000)
+        return cls.from_numpy(data=data, extent=extent)
diff --git a/cmrsim/trajectory/_flow.py b/cmrsim/trajectory/_flow.py
index c5449b9804cd6e821e3f2924e79509112cb0db43..a1d73649791f768db250a3441b3ae52acb8fa613 100644
--- a/cmrsim/trajectory/_flow.py
+++ b/cmrsim/trajectory/_flow.py
@@ -49,11 +49,11 @@ class FlowTrajectory(BaseTrajectoryModule, LookUpTableModule):
     :param map_dimensions: (3, 2) -> [(xlow, xhigh), (ylow, yhigh), (zlow, zhigh)]
                                      Physical extend of the gridded velocity fields.
     :param additional_dimension_mapping: List of name/len pairs to unstack the looked up
-                                          values. If not None this mapping must explain all
-                                          'n' additional dimension in the parameter
-                                          velocity-field e.g. ('off_res', 1) meaning that
-                                          the 1 index (starting from 3) in velocity field
-                                          describes the off-resonance at all time points.
+                                        values. If not None this mapping must explain all
+                                        'n' additional dimension in the parameter
+                                        velocity-field e.g. [('off_res', 1)] meaning that
+                                        off-resonance for all time points has a size of 1,
+                                        starting from index 3 in velocity field.
     :param kwargs: - device: str defaults to CPU:0
 
     """
@@ -70,9 +70,9 @@ class FlowTrajectory(BaseTrajectoryModule, LookUpTableModule):
         :param additional_dimension_mapping: List of name/len pairs to unstack the looked up
                                               values. If not None this mapping must explain all
                                               'n' additional dimension in the parameter
-                                              velocity-field e.g. ('off_res', 1) meaning that
-                                              the 1 index (starting from 3) in velocity field
-                                              describes the off-resonance at all time points.
+                                              velocity-field e.g. [('off_res', 1)] meaning that
+                                              off-resonance for all time points has a size of 1,
+                                              starting from index 3 in velocity field.
         :param kwargs: - device: str defaults to CPU:0
         """
 
@@ -88,6 +88,8 @@ class FlowTrajectory(BaseTrajectoryModule, LookUpTableModule):
             field_dims.insert(0, 0)
             self._field_names = {k: (field_dims[i], field_dims[i + 1])
                                  for i, k in enumerate(field_labels)}
+        else:
+            self._field_names = {}
 
         super().__init__(velocity_field, map_dimensions, name="flow_trajectory_evaluation",
                          **kwargs)
@@ -175,19 +177,75 @@ class FlowTrajectory(BaseTrajectoryModule, LookUpTableModule):
 
 # pylint: disable=abstract-method
 class TurbulentTrajectory(FlowTrajectory):
-    """ Implements the path integration for particles in a temporally static velocity field
+    """ Numerically integrates the path of particles in a mean velocity field and incorporates
+    random turbulent velocity
+    fluctuations based on Langevin updates from sampled RST values.
+    Looks up and returns additional fields for the previous particle positions.
 
+        .. dropdown:: Example usage
 
+            .. code-block:: python
+                :caption: Instantiation
 
-    """
-    def __init__(self, velocity_field: np.ndarray, map_dimensions: np.ndarray, **kwargs):
-        """
+                velocity_field_3d = ...  # shape - (X, Y, Z, 10)
+                map_dimensions = [(x_lo, x_hi), (y_lo, y_hi), (z_lo, z_hi)]
+                field_list = []  # ... no addition a field mapping
+                trajectory_module = TurbulentTrajectory(velocity_field_3d, map_dimensions, field_list)
+
+            .. code-block:: python
+                :caption: Call function
+
+                r_init = ...  # initial_positions of shape (N, 3)
+                timing = np.arange(0., 1.01, 0.01)
+                r_at_timing, fields = trajectory_module(r_init, timing, dt_max=0.01,
+                                                       return_velocities=True)
+
+            .. code-block:: python
+                :caption: increment particle positions
 
-        :param flow_field: (X, Y, Z, 10) 3D map storing gridded velocities in m/s at positions
-                                    X, Y, Z, Chol11, Chol12, Chol13, Chol22, Chol23, Chol33, Tau.
+                r_init = ...  # initial_positions of shape (N, 3)
+                delta_t = tf.constant(0.01, tf.float32)
+                r_new, fields = self.trajectory_module.increment_particles(r_init, dt=delta_t)
+
+
+        :param velocity_field: (X, Y, Z, 10+n) 3D map storing gridded fields at positions X, Y, Z.
+                                            Fields 0-2 : Mean velocities (X, Y, Z).  Unit should be m/ms
+                                            Fields 3-8 : Lower triangular Cholesky decomposition of the turbulent
+                                                         Reynolds stress tensor, L. Ordering is
+                                                         (L11, L21, L22, L31, L32, L33)
+                                            Field 9 : Langevin timescale tau, in Hz.
+                                            Additional dimensions past 9 are treated as additional
+                                            data fields and returned
         :param map_dimensions: (3, 2) -> [(xlow, xhigh), (ylow, yhigh), (zlow, zhigh)]
                                          Physical extend of the gridded velocity fields.
+        :param additional_dimension_mapping: List of name/len pairs to unstack the looked up
+                                              values. If not None this mapping must explain all
+                                              'n' additional dimension in the parameter
+                                              velocity-field e.g. [('off_res', 1)] meaning that
+                                              off-resonance for all time points has a size of 1,
+                                              starting from index 10 in velocity field.
+        :param kwargs: - device: str defaults to CPU:0
+        """
 
+    def __init__(self, velocity_field: np.ndarray, map_dimensions: np.ndarray,
+                 additional_dimension_mapping: List[Tuple[str, int]] = None, **kwargs):
+        """
+        :param velocity_field: (X, Y, Z, 10+n) 3D map storing gridded fields at positions X, Y, Z.
+                                            Fields 0-2 : Mean velocities (X, Y, Z).  Unit should be m/ms
+                                            Fields 3-8 : Lower triangular Cholesky decomposition of the turbulent
+                                                         Reynolds stress tensor, L. Ordering is
+                                                         (L11, L21, L22, L31, L32, L33)
+                                            Field 9 : Langevin timescale tau, in Hz.
+                                            Additional dimensions past 9 are treated as additional
+                                            data fields and returned
+        :param map_dimensions: (3, 2) -> [(xlow, xhigh), (ylow, yhigh), (zlow, zhigh)]
+                                         Physical extend of the gridded velocity fields.
+        :param additional_dimension_mapping: List of name/len pairs to unstack the looked up
+                                              values. If not None this mapping must explain all
+                                              'n' additional dimension in the parameter
+                                              velocity-field e.g. [('off_res', 1)] meaning that
+                                              off-resonance for all time points has a size of 1,
+                                              starting from index 10 in velocity field.
         :param kwargs: - device: str defaults to CPU:0
         """
         self.prev_v_turb = tf.Variable(tf.zeros([1, 3]), shape=(None, 3), dtype=tf.float32)
@@ -206,11 +264,29 @@ class TurbulentTrajectory(FlowTrajectory):
     @tf.function(jit_compile=False, reduce_retracing=True)
     def increment_particles(self, particle_positions: tf.Tensor, dt: tf.Tensor,
                             return_velocities: bool = False):
-        data = self.linear_interpolation(r_vectors=particle_positions)
-        impuls = tf.random.normal([tf.shape(particle_positions)[0], 3])
+        """ For given particle positions and specified time-interval returns the values of all
+        additional fields at the initial particle position and moves the particles according
+        to:
+
+        .. math::
+
+            r_{t+\\delta t} = r_{t} + \\delta t (U(r_{t}) + u_{t})
+
+            u_{t} = u_{t-\\delta t} e^{-\\frac{\\delta t}{\\tau}} + \\zeta \\sqrt{1-e^{-2 \\frac{\\delta t}{\\tau}}}
 
-        # Langevin equation
+        where U is the mean velocity in m/ms, u is fluctuating velocity in m/ms, :math:`\\tau` is the Langevin timescale
+        in ms, and :math:`\\zeta` is a random Gaussian sample from the Reynolds stress tensor.
+
+        :param particle_positions: (-1, 3) batch of 3D particle coordinates
+        :param dt: scalar value in ms
+        :param return_velocities: if True the additional field lookup contains the field "velocity"
+        :return: new particle positions (-1, 3), dict(\*\*aditional_fields[-1, c])
+        """
+        data = self.linear_interpolation(r_vectors=particle_positions)
         v_mean = data[:, 0:3]
+
+        # Random sample from RST (Cholesky decomposition)
+        impuls = tf.random.normal([tf.shape(particle_positions)[0], 3])
         vtu = impuls[:, 0] * data[:, 3]
         vtv = tf.reduce_sum(impuls[:, :2] * data[:, 4:6], axis=1)
         vtw = tf.reduce_sum(impuls * data[:, 6:9], axis=1)
@@ -227,12 +303,12 @@ class TurbulentTrajectory(FlowTrajectory):
         # New turbulence velocity
         v_turb = damping_update[:, tf.newaxis] * tf.stack([vtu, vtv, vtw], axis=1)
 
-        refmap = tf.reshape(tf.floor(data[:, 10]), [tf.shape(v_turb)[0], 1])
+        #refmap = tf.reshape(tf.floor(data[:, 10]), [tf.shape(v_turb)[0], 1])
 
         # Old + new turbulence velocity
         v_turb = self.prev_v_turb * tf.reshape(damping, [tf.shape(v_turb)[0], 1]) + v_turb
 
-        v_turb = refmap * v_turb + (refmap - 1) * self.prev_v_turb
+        #v_turb = refmap * v_turb + (refmap - 1) * self.prev_v_turb
 
         v_total = v_mean + v_turb
         additional_fields = self._unstack_add_dims(data[:, 10:])
@@ -243,8 +319,8 @@ class TurbulentTrajectory(FlowTrajectory):
         return particle_positions + v_total * dt, additional_fields
 
     def warmup_step(self, particle_positions: tf.Tensor):
-        """ Updates turbulence velocity of all particles based on sampled RST values at their
-        positions.
+        """ Updates turbulence velocity of all particles based on sampled RST values at their positions.
+        Used to initialize turbulence velocity to prevent transient states due to long Langevin timescales
 
         :param particle_positions: (-1, 3) batch of 3D particle coordinates
         """
diff --git a/cmrsim/utils/__init__.py b/cmrsim/utils/__init__.py
index c517baa1a665162e66f4ba1507cadaa570b2aff6..ba1478cabec16d7dade84f31d4291b0ab380cd57 100644
--- a/cmrsim/utils/__init__.py
+++ b/cmrsim/utils/__init__.py
@@ -6,5 +6,4 @@ __all__ = ["coordinates", "snr", "particle_properties", "display"]
 import cmrsim.utils.coordinates
 import cmrsim.utils.snr
 import cmrsim.utils.particle_properties
-import cmrsim.utils.offresonance
 import cmrsim.utils.display
diff --git a/cmrsim/utils/offresonance.py b/cmrsim/utils/offresonance.py
deleted file mode 100644
index 2266aa9898a22268da827818f82e2683784bd28c..0000000000000000000000000000000000000000
--- a/cmrsim/utils/offresonance.py
+++ /dev/null
@@ -1,44 +0,0 @@
-""" This modules contains utility to create offresonance related parameters"""
-__all__ = ["generate_inhomogeneity_map"]
-
-import numpy as np
-from pint import Quantity
-
-
-def generate_inhomogeneity_map(labels, B0: Quantity, sus_maps: dict, padsize):
-    """ Calculates inhomogeneity maps on discrete 3D maps of suseptibility values
-
-    :param labels: 3D volume oof labels
-    :param B0: Field strength in
-    :param sus_maps: dictionary that maps a suszeptibility for all values in labels
-    :param padsize:
-    :return:
-    """
-    labels_pad = np.pad(labels, padsize, constant_values=0)
-
-    phantom = np.ones_like(labels_pad).astype(np.float64) * -7e-7
-    for key, val in sus_maps.items():
-        phantom[np.where(labels_pad == key)] = val
-    # phantom = Quantity(phantom, "")
-
-    mu0 = Quantity(4. * np.pi * 1e-7, "H/m")
-
-    # Equation 5
-    center = np.ceil((np.array(phantom.shape) + 1) / 2)
-    x, z, y = np.meshgrid(*[np.arange(phantom.shape[i]) - center[i] for i in range(3)], indexing="ij")
-    corterm = 3 * z ** 2 / (x ** 2 + y ** 2 + z ** 2) - 1
-
-    # Equation 4
-    Mz = B0.to("T") * phantom / mu0
-    Mk = np.fft.fftshift(np.fft.fftn(Mz.m))
-    Bk = (-mu0 * Mk / 3. * corterm).m
-    Bk[int(center[0]), int(center[1]), int(center[2])] = -B0.m_as("T") * 2.9e-8 / 3
-
-    result = np.fft.ifftn(np.fft.ifftshift(Bk))
-    inhomogeneity = np.real(result[padsize:-padsize, padsize:-padsize, padsize:-padsize])
-    background_mean = np.mean(inhomogeneity[np.where(labels == 0)])
-    inhomogeneity -= background_mean
-
-    pmax = np.percentile(np.abs(inhomogeneity[np.where(labels != 0)]), 95)
-    inhomogeneity /= pmax
-    return inhomogeneity
diff --git a/docker/Dockerfile_cpu b/docker/Dockerfile_cpu
index 3093a3fcd8d8f540de1dd9a8e8b8cf029619282f..41a6579092dfffff4aff7d8200ba4106b320129a 100644
--- a/docker/Dockerfile_cpu
+++ b/docker/Dockerfile_cpu
@@ -2,7 +2,7 @@
 # docker build -t registry.ethz.ch/jweine/cmrsim/cpu -f Dockerfile_cpu .
 
 # Use the tensorflow image as base
-FROM tensorflow/tensorflow:2.9.1
+FROM tensorflow/tensorflow:2.11.0
 
 # maintainter from IBT
 LABEL maintainers="Jonathan Weine (weine@biomed.ee.ethz.ch)"
@@ -16,13 +16,10 @@ RUN pip3 install h5py scipy tqdm deprecated
 RUN pip3 install pyvista ffmpeg ipyvtklink panel
 
 # Install cmrseq for sequence definitions. Latest version whl needs to be provided in context
-#COPY $cmrseq_file $cmrseq_file
-#RUN pip install $cmrseq_file
-# When cmrseq will be public, this is not necessary anymore
+RUN pip install cmrseq
 
 # Install cmrsim
-COPY cmrsim-latest-py3-none-any.whl cmrsim-latest-py3-none-any.whl
-RUN pip install cmrsim-latest-py3-none-any.whl
+RUN pip install cmrsim
 
 ##
 # Environment
diff --git a/docker/Dockerfile_gpu b/docker/Dockerfile_gpu
index 4ce0bc75c714c19278955025f19111541728ea2f..d14f4e2cf15f35e687b8e5c3fa809c26d0344715 100644
--- a/docker/Dockerfile_gpu
+++ b/docker/Dockerfile_gpu
@@ -1,5 +1,7 @@
+## docker build -t registry.ethz.ch/jweine/cmrsim/gpu:latest -f ./Dockerfile_gpu .
+
 # take the GPU python3 tensorflow image as base
-FROM tensorflow/tensorflow:2.9.1-gpu
+FROM tensorflow/tensorflow:2.11.0-gpu
 
 # maintainter from IBT
 LABEL maintainers="Jonathan Weine (weine@biomed.ee.ethz.ch)"
@@ -12,14 +14,11 @@ RUN apt-get install -y libx11-6 libgl1-mesa-dev libxrender-dev
 RUN pip3 install h5py scipy tqdm deprecated
 RUN pip3 install pyvista ffmpeg ipyvtklink panel
 
-# Install cmrseq for sequence definitions. Latest version whl needs to be provided in context
-#COPY $cmrseq_file $cmrseq_file
-#RUN pip install $cmrseq_file
-# When cmrseq will be public, this is not necessary anymore
+# Install cmrseq for sequence definitions.
+RUN pip install cmrseq
 
 # Install cmrsim
-COPY cmrsim-latest-py3-none-any.whl cmrsim-latest-py3-none-any.whl
-RUN pip install cmrsim-latest-py3-none-any.whl
+RUN pip install cmrsim
 
 ##
 # Environment
diff --git a/docker/Dockerfile_jupyter b/docker/Dockerfile_jupyter
index 7c214d6fdef5ffad9ec3a0658222548d04685cc7..599357abcf7c9a494fb93698d2875f944ec0d1b5 100644
--- a/docker/Dockerfile_jupyter
+++ b/docker/Dockerfile_jupyter
@@ -1,3 +1,4 @@
+## docker build -t registry.ethz.ch/jweine/cmrsim/jupyter:latest -f ./Dockerfile_jupyter .
 FROM registry.ethz.ch/jweine/cmrsim/gpu:latest
 
 # maintainter from IBT
@@ -7,7 +8,7 @@ LABEL Description="This Docker image is meant to start MR simulations with the C
 # Set up jupyter lab with interactive plots and variable inspector
 RUN apt-get -y update
 RUN pip3 install jupyterlab
-RUN pip3 install ipympl
+RUN pip3 install ipympl==0.8.8 matplotlib==3.5
 RUN apt-get install -y curl
 RUN curl -sL https://deb.nodesource.com/setup_14.x | bash -E
 RUN apt-get install -y nodejs
@@ -16,8 +17,8 @@ RUN jupyter labextension install jupyter-matplotlib
 RUN pip install lckr-jupyterlab-variableinspector
 
 # Install pyvista and dependencies
-RUN apt install -y libgl1-mesa-glx libxrender1
-RUN pip install vtk pyvista
+RUN apt install -y libgl1-mesa-glx libxrender1 xvfb
+RUN pip install vtk pyvista phantominator
 
 RUN apt-get install -y git
 RUN apt-get install apt-transport-https ca-certificates -y
@@ -26,5 +27,4 @@ RUN update-ca-certificates
 ## Environment ##
 ENV LD_LIBRARY_PATH /usr/local/cuda/extras/CUPTI/lib64:$LD_LIBRARY_PATH
 RUN echo $LD_LIBRARY_PATH
-ENV PATH /opt/conda/bin:$PATH
 
diff --git a/docker/Dockerfile_unittest b/docker/Dockerfile_unittest
index 7082a90c75a07ac46829764d54a33a8d6b4dd6af..9e752278e3484c66287c74cffa01fbac793435fe 100644
--- a/docker/Dockerfile_unittest
+++ b/docker/Dockerfile_unittest
@@ -11,8 +11,5 @@ RUN pip install coverage anybadge phantominator genbadge
 RUN pip install pylint pylint-exit
 RUN apt install -y libgl1-mesa-glx xvfb
 
-COPY cmrseq-0.15.dev4-py3-none-any.whl cmrseq-0.15.dev4-py3-none-any.whl
-RUN pip install cmrseq-0.15.dev4-py3-none-any.whl
-
 # To push image:
 # docker push registry.ethz.ch/jweine/cmrsim:unittest
\ No newline at end of file
diff --git a/docker/cmrseq-0.15.dev4-py3-none-any.whl b/docker/cmrseq-0.15.dev4-py3-none-any.whl
deleted file mode 100644
index ec4cc3610f303d909bb39d9614bd4743af47ad2b..0000000000000000000000000000000000000000
Binary files a/docker/cmrseq-0.15.dev4-py3-none-any.whl and /dev/null differ
diff --git a/docs/source/_static/api/dataset_regular_grid_select_slice.png b/docs/source/_static/api/dataset_regular_grid_select_slice.png
new file mode 100644
index 0000000000000000000000000000000000000000..6325cc53a6b56d583bc41466348c97aea5e293b9
--- /dev/null
+++ b/docs/source/_static/api/dataset_regular_grid_select_slice.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f4d9883a1c80a8c2e212ee817d7184f77be2d51c5ea130b5296d7c4fa89e5120
+size 85981
diff --git a/docs/source/examples.rst b/docs/source/examples.rst
index 4b3c27c6440a3efab5aee2b8dbaa87b107ab2ee4..58aede87acb1ce5c4040c2df0e3fdec7210c907d 100644
--- a/docs/source/examples.rst
+++ b/docs/source/examples.rst
@@ -15,6 +15,11 @@ Datasets
 
         .. nbgallery:: example_gallery/datasets/cardiac_mesh_dataset
 
+    .. card:: Volumes on regular Grids
+
+        .. nbgallery:: example_gallery/datasets/regular_grid_dataset
+
+
 
 Trajectory Modules
 ====================
@@ -68,7 +73,7 @@ Static Phantoms
 
     .. card:: Off-resonance with BSSFP
 
-        .. nbgallery:: example_gallery/bloch_simulation/static_BSSFP_submodules.ipynb
+        .. nbgallery:: example_gallery/bloch_simulation/static_BSSFP_offresonance
 
 Flowing particles
 ^^^^^^^^^^^^^^^^^^^
@@ -92,7 +97,7 @@ Flowing particles
         .. nbgallery:: example_gallery/bloch_simulation/flow_bssfp_submodules_aorta
 
 
-Mesh Motions
+Mesh Motion
 ^^^^^^^^^^^^^
 
 .. card-carousel:: 2
@@ -120,6 +125,10 @@ Analytic Simulation
 
         .. nbgallery:: example_gallery/analytic_simulation/3_multiple_coils
 
+    .. card:: BSSFP on spherical off-resonance Phantom
+
+        .. nbgallery:: example_gallery/analytic_simulation/BSSFP_offresonance
+
     .. card:: BSSFP on texturized XCAT phantom
 
         .. nbgallery:: example_gallery/analytic_simulation/simulate_mrxcat_cmrsim
\ No newline at end of file
diff --git a/docs/source/format_datasets.py b/docs/source/format_datasets.py
index 86b05ba2fef53cdac1d28299ed4a6d92be8158f6..8e9d3b5c64281a8883af8d06a9d79414347dcc33 100644
--- a/docs/source/format_datasets.py
+++ b/docs/source/format_datasets.py
@@ -2,6 +2,7 @@ import cmrsim
 import glob
 import os
 
+
 def format_dataset_index():
     all_modules = glob.glob(f"api/datasets/cmrsim.datasets.*.rst")
     all_modules = [s.replace("api/", "") for s in all_modules]
diff --git a/notebooks/analytic_simulation/BSSFP_offresonance.ipynb b/notebooks/analytic_simulation/BSSFP_offresonance.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..cd7b41c50bdaef93b493912fa805de6ca23c3226
--- /dev/null
+++ b/notebooks/analytic_simulation/BSSFP_offresonance.ipynb
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:bb017045678b9a3a2299fa7512712d6d045d67a3edd5158adc6ca6d7cbb57964
+size 539217
diff --git a/notebooks/bloch_simulation/static_BSSFP_offresonance.ipynb b/notebooks/bloch_simulation/static_BSSFP_offresonance.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..94eb5d2c75ae82a0bc06064a4480c15c70ccf88b
--- /dev/null
+++ b/notebooks/bloch_simulation/static_BSSFP_offresonance.ipynb
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:2c0a3d4bc78511d90712facf3dd3de544e92fe84c721638e5fa6ded3f943b954
+size 919308
diff --git a/notebooks/bloch_simulation/static_BSSFP_submodules.ipynb b/notebooks/bloch_simulation/static_BSSFP_submodules.ipynb
deleted file mode 100644
index 6dd7c66edc7a8ccd5be35517eecfc2f50872402d..0000000000000000000000000000000000000000
--- a/notebooks/bloch_simulation/static_BSSFP_submodules.ipynb
+++ /dev/null
@@ -1,3 +0,0 @@
-version https://git-lfs.github.com/spec/v1
-oid sha256:78ee646487c91de39dabe4fc4f65feb6271627c392654a05f21c099b546ca101
-size 15550702
diff --git a/notebooks/bloch_simulation/static_multi_coil_bssfp.ipynb b/notebooks/bloch_simulation/static_multi_coil_bssfp.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..fb272988219409e330da1b4f8c8d2051ff42c818
--- /dev/null
+++ b/notebooks/bloch_simulation/static_multi_coil_bssfp.ipynb
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:e2b13574fc91146c330f7e71c18a632e31e39f054e6960dc0446c60757205696
+size 2106036
diff --git a/notebooks/datasets/regular_grid_dataset.ipynb b/notebooks/datasets/regular_grid_dataset.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..1e4b62e1c44bfb96c7e3db4a5e15efad74f720ce
--- /dev/null
+++ b/notebooks/datasets/regular_grid_dataset.ipynb
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:389468d74ee83518a1045a2b2bc11189c8c333839623de4bc2d54c31125da0b3
+size 324270
diff --git a/notebooks/local_functions.py b/notebooks/local_functions.py
index b71f1db5277a42671d7c76b1aec33c8b22dd4028..a4d93a0351ff94c91eb6b8c2074bfb3b0ee73505 100644
--- a/notebooks/local_functions.py
+++ b/notebooks/local_functions.py
@@ -1,22 +1,105 @@
+__all__ = ["transformed_box", "create_spherical_phantom", "add_custom_axes", 
+           "add_refilling_box", "animate_trajectories", "get_custom_theme"]
 from typing import Sequence
 
 import pyvista 
 import numpy as np
+from pint import Quantity
 from tqdm.notebook import tqdm
 
 
-def add_custom_axes(plotter):
+def transformed_box(reference_basis: np.ndarray, slice_normal: np.ndarray,
+                    readout_direction: np.ndarray, slice_position: Quantity,
+                    extend: Quantity = Quantity([30, 30, 1], "cm")):
+    """
+    
+    :param reference_basis: (3x3) matrix where row-index reference the x, y, z
+                            basis vectors
+    
+    """
+    from scipy.spatial.transform import Rotation
+    slice_normal /= np.linalg.norm(slice_normal)
+
+    if readout_direction is None:
+        if np.allclose(slice_normal, np.array([1., 0., 0.]), atol=1e-3):
+            readout_direction = np.array([0., 1., 0.])
+        else:
+            readout_direction = np.array([1., 0., 0.])
+    else:
+
+        if np.isclose(
+                np.dot(slice_normal, readout_direction / np.linalg.norm(readout_direction)),
+                1., atol=1e-3):
+            raise ValueError("Slice normal and readout-direction are parallel")
+
+    readout_direction = readout_direction - np.dot(slice_normal,
+                                                   readout_direction) * slice_normal
+    readout_direction /= np.linalg.norm(readout_direction)
+    pe_direction = np.cross(readout_direction, slice_normal)
+    pe_direction /= np.linalg.norm(pe_direction)
+    new_basis = np.stack([readout_direction, pe_direction, slice_normal])
+    rot, _ = Rotation.align_vectors(reference_basis, new_basis)
+
+    total_transform = np.eye(4, 4)
+    total_transform[:3, :3] = rot.as_matrix()
+    total_transform[:3, 3] = slice_position.m_as("m")
+    bounds = np.stack([-extend.m_as("m")/2, extend.m_as("m")/2], axis=1).reshape(-1)    
+    
+    box = pyvista.Box(bounds=bounds)
+    box.transform(total_transform, inplace=True)
+    return box
+
+
+def create_spherical_phantom(dimensions=(120, 120, 120), spacing=(0.002, 0.002, 0.002),
+                             big_radius=0.1) -> pyvista.UnstructuredGrid:
+    """ Creates a spherical phantom with 6 spherical holes on the coordinate axes and
+     a rectangular hole in the center.
+
+    :param dimensions: (int, int, int)
+    :param spacing: (float, float, float) in meter
+    :param big_radius: float in meter
+    :return: pyvista.UnstructuredGrid
+    """
+    raster = pyvista.UniformGrid(dimensions=dimensions, spacing=spacing,
+                                 origin=-np.array(dimensions) * np.array(spacing) / 2)
+    sphere_big = pyvista.Sphere(radius=big_radius).extract_surface()
+    box = pyvista.Box(bounds=(-big_radius * 0.3, big_radius * 0.3, -big_radius/5,
+                              big_radius/5, -big_radius/10, big_radius/10))
+
+    o = big_radius * 0.7
+    spheres_small = [pyvista.Sphere(radius=big_radius*0.05, center=c).extract_surface() for c in
+                     ((-o, 0., 0.), (o, 0., 0.), (0., -o, 0.), (0., o, 0.), (0., 0., -o),
+                      (0., 0., o))]
+    del o
+
+    filled_sphere = raster.clip_surface(sphere_big, progress_bar=True)
+    filled_sphere = filled_sphere.clip_surface(box, progress_bar=True, invert=False)
+    for sms in spheres_small:
+        filled_sphere = filled_sphere.clip_surface(sms, progress_bar=True, invert=False)
+    return filled_sphere
+
+
+def add_custom_axes(plotter: pyvista.Plotter, scale: float = 0.1):
     """ Add scaled arrows to the plotter at origin"""
-    origin_arrows = [pyvista.Arrow(start=(0., 0., 0.), direction=direction, tip_length=0.1, tip_radius=0.025, 
-                                   tip_resolution=10, shaft_radius=0.005, shaft_resolution=10, scale=0.1) 
+    origin_arrows = [pyvista.Arrow(start=(0., 0., 0.), direction=direction, tip_length=0.1,
+                                   tip_radius=0.025, tip_resolution=10, shaft_radius=0.005,
+                                   shaft_resolution=10, scale=scale)
                      for direction in [(1., 0., 0.), (0., 1., 0.), (0., 0., 1.)]]
     plotter.add_mesh(origin_arrows[0], color="blue", label="x")
     plotter.add_mesh(origin_arrows[1], color="green", label="y")
     plotter.add_mesh(origin_arrows[2], color="red", label="z")
     
     
-def add_refilling_box(plotter, bounding_box, thickness, position):
-    """ Adds an orange box oriented in (x,y,z) to the plotter, to illustrate a reseeding volume"""
+def add_refilling_box(plotter: pyvista.Plotter, bounding_box: Quantity,
+                      thickness: Quantity, position: Quantity):
+    """ Adds an orange box oriented in (x,y,z) to the plotter, to illustrate a reseeding volume
+
+    :param plotter: -
+    :param bounding_box: (2, )
+    :param thickness: ( )
+    :param position: (3, )
+    :return:
+    """
     slice_mesh = pyvista.Box(bounds=(-bounding_box[0].m_as("m")/2, bounding_box[0].m_as("m")/2,
                                      -bounding_box[1].m_as("m")/2, bounding_box[1].m_as("m") /2,
                                      -thickness.m_as("m")/2, thickness.m_as("m")/2))
@@ -24,8 +107,8 @@ def add_refilling_box(plotter, bounding_box, thickness, position):
     plotter.add_mesh(slice_mesh, color="orange", opacity=0.5, show_edges=True)
     
     
-def animate_trajectories(plotter, trajectories: Sequence, timing: np.ndarray,
-                         filename: str, scalars: Sequence = None, 
+def animate_trajectories(plotter: pyvista.Plotter, trajectories: Sequence,
+                         timing: np.ndarray, filename: str, scalars: Sequence = None,
                          vectors: Sequence = None,
                          multi_vecs: bool = False,
                          vector_pos: Sequence = None,
@@ -33,66 +116,83 @@ def animate_trajectories(plotter, trajectories: Sequence, timing: np.ndarray,
                          text_kwargs: dict = None,
                          vector_kwargs: dict = None,
                          visualize_trajectories: bool = False,
-                         trajectory_kwargs: dict = None
-                        ):
-
-        if mesh_kwargs is None:
-            mesh_kwargs = dict()
-        
-        if text_kwargs is None:
-            text_kwargs = dict(color='white', shadow=False, font_size=26,
-                               position='upper_right')
-        if vector_kwargs is None:
-            vector_kwargs = dict(mag=1e-3)
-        
-        if vector_pos is None:
-            vector_pos = trajectories
-            
-        if trajectory_kwargs is None:
-            trajectory_kwargs = dict(color="red", show_scalar_bar=False)
-        
-        plotter.open_gif(f"{filename}")        
+                         trajectory_kwargs: dict = None):
+    """Animates Sequence of trajectories and saves the animation as gif
+
+    :param plotter:
+    :param trajectories: sequence containing positions of arbitrary points. Each entry corresponds
+            to a single snapshot
+    :param timing: sequence of time-values corresponding to the snapshots in trajectories
+    :param filename: full path of the saved file (must include .gif ending)
+    :param scalars: scalar values corresponding to each entry of trajectories (shapes must match)
+    :param vectors: vectors corresponding to each entry of vector pos (shapes must match)
+    :param multi_vecs: bool - if true it is assumed that each entry of 'vectors' argument contains
+                    multiple vectors per time-step / per point
+    :param vector_pos: positions of vectors specified in argument 'vectors'
+    :param mesh_kwargs: dictionary containing arguments to plotter.add_mesh(.., **mesh_kwargs)
+    :param text_kwargs: dictionary containing arguments to plotter.add_text(..., **text_kwargs)
+    :param vector_kwargs: dictionary containing arguments to plotter.add_arrows(..., **vector_kwargs)
+    :param visualize_trajectories: if true, particle trajectories are traced with a line
+    :param trajectory_kwargs: dictionary containing argument for lines tracing the trajectories
+    """
+
+    if mesh_kwargs is None:
+        mesh_kwargs = dict()
+
+    if text_kwargs is None:
+        text_kwargs = dict(color='white', shadow=False, font_size=26,
+                           position='upper_right')
+    if vector_kwargs is None:
+        vector_kwargs = dict(mag=1e-3)
+
+    if vector_pos is None:
+        vector_pos = trajectories
+
+    if trajectory_kwargs is None:
+        trajectory_kwargs = dict(color="red", show_scalar_bar=False)
+
+    plotter.open_gif(f"{filename}")
+    if timing is not None:
+        desc = f"Rendering Mesh-Animation from {timing[0]:2.3f}ms to {timing[-1]:2.3f} ms"
+
+    pbar = tqdm(range(len(trajectories)), desc=desc)
+    for step in pbar:
+        plot_mesh = pyvista.PolyData(trajectories[step])
+
         if timing is not None:
-            desc = f"Rendering Mesh-Animation from {timing[0]:2.3f}ms to {timing[-1]:2.3f} ms"
-            
-        pbar = tqdm(range(len(trajectories)), desc=desc)
-        for step in pbar:
-            plot_mesh = pyvista.PolyData(trajectories[step])
-            
-            if timing is not None:
-                text_actor = plotter.add_text(f"{timing[step]:>2.2f} ms", **text_kwargs)
-            if scalars is not None:
-                point_actor = plotter.add_mesh(plot_mesh, scalars=scalars[step], **mesh_kwargs)
+            text_actor = plotter.add_text(f"{timing[step]:>2.2f} ms", **text_kwargs)
+        if scalars is not None:
+            point_actor = plotter.add_mesh(plot_mesh, scalars=scalars[step], **mesh_kwargs)
+        else:
+            point_actor = plotter.add_mesh(plot_mesh, **mesh_kwargs)
+
+        if vectors is not None:
+            if multi_vecs:
+                arrow_actors = [plotter.add_arrows(cent=vector_pos[step], direction=vectors[i][step],
+                                                   **vector_kwargs[i]) for i in range(len(vectors))]
             else:
-                point_actor = plotter.add_mesh(plot_mesh, **mesh_kwargs)
-
-            if vectors is not None:
-                if multi_vecs:
-                    arrow_actors = [plotter.add_arrows(cent=vector_pos[step], direction=vectors[i][step], 
-                                                       **vector_kwargs[i]) for i in range(len(vectors))]
-                else:
-                    arrow_actor = plotter.add_arrows(cent=vector_pos[step],
-                                                     direction=vectors[step],
-                                                     **vector_kwargs)
-                    
-            if visualize_trajectories:
-                if step >= 1:
-                    plot_mesh["disp"] = trajectories[step]-trajectories[step-1]
-                    glyphs = plot_mesh.glyph(orient="disp", geom=pyvista.Line())
-                    plotter.add_mesh(glyphs, **trajectory_kwargs)
-    
-            plotter.render()
-            plotter.write_frame()
-            if timing is not None:
-                plotter.remove_actor(text_actor)
-            if vectors is not None:
-                if multi_vecs:
-                    for act in arrow_actors:
-                        plotter.remove_actor(act)
-                else:
-                    plotter.remove_actor(arrow_actor)
-
-            plotter.remove_actor(point_actor)
+                arrow_actor = plotter.add_arrows(cent=vector_pos[step],
+                                                 direction=vectors[step],
+                                                 **vector_kwargs)
+
+        if visualize_trajectories:
+            if step >= 1:
+                plot_mesh["disp"] = trajectories[step]-trajectories[step-1]
+                glyphs = plot_mesh.glyph(orient="disp", geom=pyvista.Line())
+                plotter.add_mesh(glyphs, **trajectory_kwargs)
+
+        plotter.render()
+        plotter.write_frame()
+        if timing is not None:
+            plotter.remove_actor(text_actor)
+        if vectors is not None:
+            if multi_vecs:
+                for act in arrow_actors:
+                    plotter.remove_actor(act)
+            else:
+                plotter.remove_actor(arrow_actor)
+
+        plotter.remove_actor(point_actor)
     
     
 def get_custom_theme():
diff --git a/requirements.txt b/requirements.txt
index b6a6a5ecfd7e66d8b1929e41e6aef3862ed1fe26..c2efbd39f879cb8c33d6de4b3d77fd979f10ad46 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -8,5 +8,5 @@ tqdm
 scipy
 setuptools
 deprecated
-cmrseq>=0.14
+cmrseq>=0.16
 
diff --git a/setup.py b/setup.py
index 716b991a20173fe7388c36d519147b0e28cb9ed4..79ca2eccb4e10b5cf51b3c52c220971242726fe7 100644
--- a/setup.py
+++ b/setup.py
@@ -16,6 +16,12 @@ author_list = "Jonathan Weine, Charles McGrath"
 author_email_list = "weine@biomed.ee.ethz.ch, mcgrath@biomed.ee.ethz.ch"
 url = "https://gitlab.ethz.ch/jweine/cmrsim"
 
+project_urls = {
+    'Documentation': 'https://people.ee.ethz.ch/~jweine/cmrsim/latest/index.html',
+    'Source': 'https://gitlab.ethz.ch/jweine/cmrsim',
+    'Institute': 'https://cmr.ethz.ch/'
+}
+
 # Get version tag
 if '--version' in sys.argv:
     tag_index = sys.argv.index('--version') + 1
@@ -33,13 +39,14 @@ with open(f'{project_name}/__init__.py', 'r+') as module_header_file:
 setuptools.setup(
     name=project_name,
     url=url,
+    project_urls=project_urls,
     version=current_version,
     author=author_list,
     author_email_=author_email_list,
     long_description=long_description,
     packages=setuptools.find_packages(exclude=["*.tests", "*.tests.*", "tests.*", "tests"]),
     include_package_data=True,
-    install_requires=["tensorflow>=2.9", "pyvista>=0.34", "numpy>=0.19", "pint>=0.18",
-                      "Deprecated>=1.2.13"],
+    install_requires=["tensorflow>=2.11", "pyvista>=0.37", "numpy>=1.19", "pint>=0.18",
+                      "Deprecated>=1.2.13", "cmrseq>=0.17"],
     python_requires=">=3.6"
 )
diff --git a/tests/test_bloch_generic.py b/tests/test_bloch_generic.py
index ecfc3285751b750b4126c318f424ba1b8a9ff6da..e20cf0ee30730caaa7508b4078e2974f153f1f81 100644
--- a/tests/test_bloch_generic.py
+++ b/tests/test_bloch_generic.py
@@ -92,47 +92,46 @@ class TestGenericBlochModule(unittest.TestCase):
                                          max_grad=Quantity(80., "mT/m"),
                                          grad_raster_time=Quantity(0.005, "ms"),
                                          rf_raster_time=Quantity(0.005, "ms"))
-        seq = cmrseq.seqdefs.diffusion.bipolar(system_specs=system_specs,
-                                               dt=Quantity(0.1, "ms"),
-                                               Dt=Quantity(10., "ms"),
-                                               direction=np.array([1., 0., 0.]),
-                                               amplitude=Quantity(1., "mT/m"),
-                                               rise_time=Quantity(0.03,  "ms"))
+        seq = cmrseq.seqdefs.velocity.bipolar(system_specs=system_specs,
+                                              venc=Quantity(200, "cm/s"),
+                                              duration=Quantity(2.5, "ms"),
+                                              direction=np.array([-1., 0., 0.]))
         t, gridded_seq = seq.gradients_to_grid()
         wf = gridded_seq.T
         dt = system_specs.grad_raster_time.m_as("ms")
         module = cmrsim.bloch.GeneralBlochOperator("test_grad",
-                                                   system_specs.gamma_rad.m_as("rad/ms/mT"),
+                                                   gamma=system_specs.gamma_rad.m_as("rad/ms/mT"),
                                                    time_grid=t,
                                                    gradient_waveforms=wf[np.newaxis])
 
         class TrajMod(tf.Module):
-            def increment_particles(self, r, dt):
-                v = Quantity([25, 0., 0.], "cm/s")
-                v = tf.constant(v.m_as("m/ms"), dtype=tf.float32, shape=(1, 3))
-                return r + v * dt, dict()
+            def increment_particles(self, particle_positions: tf.Tensor, dt: tf.Tensor):
+                v = tf.constant(Quantity([25, 0., 0.], "cm/s").m_as("m/ms"),
+                                dtype=tf.float32, shape=(1, 3))
+                return particle_positions + v * dt, dict()
+    
+
 
         with self.subTest("Mean Phase increase due to motion in direction of bipolar gradient"):
             # Set up phantom
             n_particles = 100
-            r_init = tf.random.normal((n_particles, 3), 0., 0.05)
+            r_init = tf.random.normal((n_particles, 3), 0., 0.005)
             m = tf.tile(tf.constant([[1., 0., 0.], ], dtype=tf.complex64), [n_particles, 1])
             properties = OrderedDict(
                 magnetization=m,
                 T1=tf.ones((n_particles,), dtype=tf.float32) * Quantity(1000, "ms").m,
-                T2=tf.ones((n_particles,), dtype=tf.float32) * Quantity(300, "ms").m,
+                T2=tf.ones((n_particles,), dtype=tf.float32) * Quantity(500, "ms").m,
                 M0=tf.ones((n_particles,), dtype=tf.float32) * 10)
 
             # call module
-            m_new, _ = module(trajectory_module=TrajMod(),
+            m_new, r = module(trajectory_module=TrajMod(),
                               initial_position=r_init,
                               repetition_index=0,
                               run_parallel=True, **properties)
-
+            
             mean_phase = tf.reduce_mean(tf.math.angle(m_new[0, :, 0]))
             relaxation_factor = 1. - np.exp(-(wf.shape[0]-1)*dt/1000)  # must match T1 definition
             self.assertGreater(tf.abs(mean_phase), tf.reduce_mean(tf.math.angle(m[:, 0])))
-            self.assertTrue(np.allclose(tf.math.angle(m_new[0, :, 0]).numpy(), mean_phase.numpy(),
-                                        atol=1e-4))
+            self.assertTrue(np.allclose(tf.math.angle(m_new[0, :, 0]).numpy(), mean_phase.numpy(), atol=1e-3))
             self.assertTrue(np.allclose(tf.abs(m_new[0, :, 2]), relaxation_factor, atol=1e-4))
 
diff --git a/tests/test_bloch_multi_coil.py b/tests/test_bloch_multi_coil.py
new file mode 100644
index 0000000000000000000000000000000000000000..3f284f9a28b48208d5e6ab31303b73c87e5dd92a
--- /dev/null
+++ b/tests/test_bloch_multi_coil.py
@@ -0,0 +1,95 @@
+import unittest
+from collections import OrderedDict
+from copy import deepcopy
+
+import numpy as np
+import tensorflow as tf
+tf.config.run_functions_eagerly(True)
+
+from pint import Quantity
+import cmrseq
+
+import cmrsim
+
+
+class TestParallelReceiveBlochSimulation(unittest.TestCase):
+    def setUp(self) -> None:
+        self.n_reps = 10
+        self.n_steps = 100
+        self.time = tf.range(0., self.n_steps, dtype=tf.float32) * 0.01
+        self.grads = tf.ones([self.n_reps, self.n_steps, 3], dtype=tf.float32)
+        self.rf = tf.zeros((self.n_reps, self.n_steps), dtype=tf.complex64)
+        adc_on = tf.roll(tf.eye(self.n_reps, self.n_steps,
+                                dtype=tf.float32).numpy(), 1, axis=1)
+        adc_phase = tf.zeros_like(adc_on)
+        self.adc = tf.stack([adc_on, adc_phase], axis=-1)
+
+        n_particles = 1000
+        self.simulation_input = dict(
+            magnetization=tf.zeros([n_particles, 3], dtype=tf.complex64),
+            T1=tf.ones([n_particles, ], dtype=tf.float32),
+            T2=tf.ones([n_particles, ], dtype=tf.float32),
+            M0=tf.ones([n_particles, ], dtype=tf.float32),
+            initial_position=tf.random.normal([n_particles, 3], 0., 0.005),
+            off_res=tf.random.normal([n_particles, 1], 0, 1, dtype=tf.float32),
+            omega_t2s=tf.random.normal([n_particles, ], 0., 0.1, dtype=tf.float32)
+        )
+        self.n_particles = n_particles
+
+        self.coil_module = cmrsim.analytic.contrast.CoilSensitivity(
+                            coil_sensitivities=np.ones([1, 20, 20, 20, 3], dtype=np.complex64),
+                            map_dimensions=np.array([[-0.1, 0.1], [-0.1, 0.1], [-0.1, 0.1]]))
+
+    def test_input_shapes(self):
+
+        def _test_sequential_call(module: cmrsim.bloch.ParallelReceiveBlochOperator):
+            m, r = module(repetition_index=tf.constant(0),
+                          run_parallel=False, **self.simulation_input)
+            self.assertTrue(all([i == j for i, j in zip(m.shape, [self.n_particles, 3])]))
+            self.assertTrue(m.dtype is tf.complex64)
+            self.assertTrue(r.dtype is tf.float32)
+
+        def _test_parallel_call(module: cmrsim.bloch.ParallelReceiveBlochOperator):
+            m, r = module(run_parallel=True, **self.simulation_input)
+            self.assertTrue(all([i == j for i, j in
+                                 zip(m.shape, [self.n_reps, self.n_particles, 3])]))
+            self.assertTrue(m.dtype is tf.complex64)
+            self.assertTrue(r.dtype is tf.float32)
+
+            second_iteration_input = deepcopy(self.simulation_input)
+            [second_iteration_input.pop(k) for k in ("magnetization", "initial_position")]
+            m, r = module(run_parallel=True, magnetization=m, initial_position=r,
+                          **second_iteration_input)
+            self.assertTrue(all([i == j for i, j in
+                                 zip(m.shape, [self.n_reps, self.n_particles, 3])]))
+            self.assertTrue(m.dtype is tf.complex64)
+            self.assertTrue(r.dtype is tf.float32)
+
+        # test __call__ without submodules
+        module_ = cmrsim.bloch.ParallelReceiveBlochOperator("test", 1., self.coil_module,
+                                                            time_grid=self.time,
+                                                            gradient_waveforms=self.grads,
+                                                            adc_events=self.adc,
+                                                            rf_waveforms=self.rf)
+        with self.subTest("No submodules, sequential simulation"):
+            _test_sequential_call(module_)
+
+        with self.subTest("No submodules, parallel simulation"):
+            _test_parallel_call(module_)
+
+        # test __call__ with submodules
+        submodules = (cmrsim.bloch.submodules.OffResonance(1),
+                      cmrsim.bloch.submodules.T2starDistributionModel(),
+                      cmrsim.bloch.submodules.ConcomitantFields(1, 1),
+                      )
+        module_ = cmrsim.bloch.ParallelReceiveBlochOperator("test", 1., self.coil_module,
+                                                            time_grid=self.time,
+                                                            gradient_waveforms=self.grads,
+                                                            adc_events=self.adc,
+                                                            rf_waveforms=self.rf,
+                                                            submodules=submodules)
+        with self.subTest("With all submodules sequential call"):
+            _test_sequential_call(module_)
+
+        with self.subTest("With all submodule parallel call"):
+            _test_parallel_call(module_)
diff --git a/tests/test_datasets.py b/tests/test_datasets.py
index e252a8ff3d72b04c6914f5117788d7e7da7e8d0a..1b1ad2036ce3ba622baf2734396359ed0371131f 100644
--- a/tests/test_datasets.py
+++ b/tests/test_datasets.py
@@ -203,7 +203,6 @@ class TestMeshDatasets(unittest.TestCase):
         slice_mesh.add_data_per_time(np.random.normal(0., 1., size=(3, npoints, 3)),
                                      name="random_vector")
 
-
         from sys import platform
         if platform == "linux" or platform == "linux2":
             pyvista.start_xvfb()
@@ -220,3 +219,64 @@ class TestMeshDatasets(unittest.TestCase):
         mesh_dataset.evaluate_cylinder_coordinates(mesh_dataset.mesh.points, 4, 40*30+1)
         mesh_dataset.compute_local_basis(time_stamp=Quantity(0., "ms"))
         mesh_dataset.compute_cylinder_coords(time_stamp=Quantity(0., "ms"))
+
+
+class TestRegularGridDataset(unittest.TestCase):
+    def test_shepp_logan(self):
+        dataset = cmrsim.datasets.RegularGridDataset.from_shepp_logan(
+            map_size=(100, 100, 5), extent=Quantity([10, 10, 5], "cm"))
+
+    def test_from_unstructured_grid(self):
+        sphere = pyvista.Sphere(radius=0.1, center=(0., 0., 0.))
+        ds = cmrsim.datasets.RegularGridDataset.from_unstructured_grid(
+            sphere, pixel_spacing=Quantity([4., 4., 4.], "mm"),
+            padding=Quantity([1., 1., 1.], "cm"))
+
+    def test_general_functionality(self):
+        sphere = pyvista.Sphere(radius=0.1, center=(0., 0., 0.))
+        ds = cmrsim.datasets.RegularGridDataset.from_unstructured_grid(
+            sphere, pixel_spacing=Quantity([4., 4., 4.], "mm"),
+            padding=Quantity([1., 1., 1.], "cm"))
+
+        ds.mesh["M0"] = np.where(ds.mesh["vtkValidPointMask"],
+                                  np.ones_like(ds.mesh["vtkValidPointMask"]),
+                                  np.zeros_like(ds.mesh["vtkValidPointMask"]))
+
+        with self.subTest("compute_offres"):
+            # susceptibilty in ppm
+            ds.mesh["chi"] = np.where(ds.mesh["vtkValidPointMask"],
+                                      np.ones_like(ds.mesh["vtkValidPointMask"]) * (-9.05),
+                                      np.ones_like(ds.mesh["vtkValidPointMask"]) * 0.36) * 1e-6
+            ds.compute_offresonance(b0=Quantity(3, "T"), susceptibility_key="chi")
+
+        with self.subTest("select_slice"):
+            slice_normal = np.array([0., 1., 1.])
+            readout_direction = np.array([1., 1., 0.])
+            slice_position = Quantity([0., 0., 1], "cm")
+            field_of_view = Quantity([30, 30, 1], "cm")
+            spacing = Quantity([1, 1, 1], "mm")
+
+            global_slice = ds.select_slice(slice_normal, spacing, slice_position,
+                                           field_of_view, readout_direction, in_mps=False)
+            mps_slice = ds.select_slice(slice_normal, spacing, slice_position, field_of_view,
+                                        readout_direction, in_mps=True)
+            global_slice = ds.select_slice(slice_normal, spacing, slice_position,
+                                           field_of_view, in_mps=False)
+            self.assertRaises(ValueError, lambda: ds.select_slice(slice_normal, spacing,
+                                        slice_position, field_of_view, slice_normal, in_mps=False))
+
+        with self.subTest("get_phantom_dict"):
+            ds.mesh["is_non_trivial"] = ds.mesh["M0"] > 0.
+            phantom_dict = ds.get_phantom_def(filter_by="is_non_trivial", keys=["M0", "offres"],
+                                              perturb_positions_std=Quantity(0.1, "mm").m_as("m"))
+
+        with self.subTest("add_field"):
+            dat = np.ones(ds.mesh.dimensions)
+            ds.add_field(key="dummy", data=dat)
+            self.assertTrue("dummy" in ds.mesh.array_names)
+
+        with self.subTest("resample_spacing"):
+            ds.resample_spacing(spacing=Quantity([3, 3, 3], "mm"))
+
+
+
diff --git a/tests/test_signal_so_phasetracking.py b/tests/test_signal_so_phasetracking.py
index 6864d00cdd96dbe8efb8771d49a21fb328d17979..51968eb055c01524b0ab39380c6dcc2a2ee8722a 100644
--- a/tests/test_signal_so_phasetracking.py
+++ b/tests/test_signal_so_phasetracking.py
@@ -26,7 +26,9 @@ class TestDiffusionWaveFormTracking(unittest.TestCase):
         self.system_specs = cmrseq.SystemSpec(max_grad=Quantity(80., "mT/m"),
                                               max_slew=Quantity(80., "mT/m/ms"),
                                               grad_raster_time=Quantity(0.01, "ms"))
-        seq = cmrseq.seqdefs.diffusion.m012(self.system_specs, Quantity(100., "s/mm**2"),
+        seq = cmrseq.seqdefs.diffusion.m012(self.system_specs, lambda_=Quantity(5.6, "ms"),
+                                            zeta=Quantity(0.9, "ms"),
+                                            bvalue=Quantity(100., "s/mm**2"),
                                             direction=np.array([0., 1., 0.]))
         t, wf = seq.gradients_to_grid()
         wf = np.concatenate([t[np.newaxis], wf], axis=0).T
@@ -53,7 +55,6 @@ class TestDiffusionWaveFormTracking(unittest.TestCase):
         signal_tensor = np.ones((rad_cords.shape[0], 1, 1), dtype=np.complex64)
         signal_tensor = self.mod(signal_tensor, trajectory)
         phase = np.angle(signal_tensor)
-        print(phase)
         self.assertTrue(np.allclose(phase, 0., atol=np.pi/1000),
                         msg="Sum of all phases is greater than numerical precision "
                             "tolerance for static phantom")
diff --git a/tests/test_signal_so_sequences.py b/tests/test_signal_so_sequences.py
index a7e6b3f212b4bbc2ffbcd99d47a7517b5e030c8b..1b59a3b4eb2fdf0d89dd8fcc49c78e56c0e1d7ea 100644
--- a/tests/test_signal_so_sequences.py
+++ b/tests/test_signal_so_sequences.py
@@ -86,18 +86,20 @@ class TestBSSFP(unittest.TestCase):
         T2star = tf.ones(shape=(10, 1, 1), dtype=tf.float32)
         offres = tf.zeros(shape=(10, 1, 1), dtype=tf.float32)
 
-        a = bssfp_module(m0, T1=T1, T2=T2, T2star=T2star, theta=offres)
+        a = bssfp_module(m0, T1=T1, T2=T2, T2star=T2star, off_res=offres)
 
     def test_banding(self):
-        theta_grid = tf.range(-2 * np.pi, 2 * np.pi, np.pi / 50, dtype=tf.float32)
+        TR = 6  # ms
+        theta_grid = tf.range(-2 * np.pi, 2 * np.pi, np.pi / 50, dtype=tf.float32) / TR
+
         sequence_module = BSSFP(flip_angle=np.array([35, ]), echo_time=[3, ],
-                                repetition_time=[6, ], expand_repetitions=True)
-        inp1 = dict(signal_tensor=tf.ones([1, 1, 1], dtype=tf.complex64), T1=tf.ones([1, 1, 1], dtype=tf.float32) * 600,
+                                repetition_time=[TR, ], expand_repetitions=True)
+        inp1 = dict(signal_tensor=tf.ones([1, 1, 1], dtype=tf.complex64), T1=tf.ones([1, 1, 1] ,dtype=tf.float32) * 600,
                     T2=tf.ones([1, 1, 1], dtype=tf.float32) * 100, T2star=tf.ones([1, 1, 1], dtype=tf.float32) * 1000,
-                    theta=tf.reshape(theta_grid, [1, -1, 1]))
+                    off_res=tf.reshape(theta_grid, [1, -1, 1]))
         inp2 = dict(signal_tensor=tf.ones([1, 1, 1], dtype=tf.complex64), T1=tf.ones([1, 1, 1], dtype=tf.float32) * 900,
                     T2=tf.ones([1, 1, 1], dtype=tf.float32) * 50, T2star=tf.ones([1, 1, 1], dtype=tf.float32) * 1000,
-                    theta=tf.reshape(theta_grid, [1, -1, 1]))
+                    off_res=tf.reshape(theta_grid, [1, -1, 1]))
 
         res1 = sequence_module(**inp1)
         res2 = sequence_module(**inp2)