diff --git a/CHANGELOG.md b/CHANGELOG.md
index 03a1d366..a8e32589 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,7 +4,24 @@ All notable changes to BEAT will be documented in this file.
 
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 
-## [1.2.1] TBD 2022
+
+## [1.2.2] TBD 2022
+Contributors: Hannes Vasyura-Bathke @hvasbath
+
+### Added
+- plotting.marginals.traceplot: source_idxs argument can be slice e.g. 1:10 to take mean of patch attributes
+- heart.seis_derivative: function to numerically calculate derivatives wrt source attributes for waveforms
+- utility: functions for slice to string conversion and vice versa
+- config; NonlinearGFConfig added version attribute to specify backend version to use for GF calculation
+
+### Changed
+- plotting.marginals.traceplot: CDF plotting of multiple distributions in a single subplot marks quantiles only once
+
+### Fixed
+- plotting.marginals.traceplot: multipage output repeating variables
+
+
+## [1.2.1] 14 September 2022
 Contributors: Hannes Vasyura-Bathke @hvasbath
 
 ### Added
diff --git a/beat/apps/beat.py b/beat/apps/beat.py
index 3c4efe41..697c7abb 100644
--- a/beat/apps/beat.py
+++ b/beat/apps/beat.py
@@ -35,7 +35,7 @@
 from beat.sampler.pt import SamplingHistory
 from beat.sampler.smc import sample_factor_final_stage
 from beat.sources import MTQTSource, MTSourceWithMagnitude
-from beat.utility import list2string
+from beat.utility import list2string, string2slice
 
 logger = logging.getLogger("beat")
 
@@ -1862,7 +1862,14 @@ def setup(parser):
     if len(options.source_idxs) < 1:
         source_idxs = None
     else:
-        source_idxs = [int(idx) for idx in options.source_idxs]
+        source_idxs = []
+        for str_idx in options.source_idxs:
+            try:
+                idx = int(str_idx)
+            except ValueError:
+                idx = string2slice(str_idx)
+
+            source_idxs.append(idx)
 
     po = plotting.PlotOptions(
         plot_projection=options.plot_projection,
diff --git a/beat/config.py b/beat/config.py
index 0425bedf..a1727c07 100644
--- a/beat/config.py
+++ b/beat/config.py
@@ -380,6 +380,10 @@ class NonlinearGFConfig(GFConfig):
         "depth the velocity model is not varied based on the errors "
         "defined above!",
     )
+    version = String.T(
+        default="",
+        help="Version number of the backend codes. If not defined, default versions will be used.",
+    )
 
 
 class SeismicGFConfig(NonlinearGFConfig):
diff --git a/beat/heart.py b/beat/heart.py
index ab9fa4d8..65484bba 100644
--- a/beat/heart.py
+++ b/beat/heart.py
@@ -2160,7 +2160,12 @@ def get_fomosto_baseconfig(gfconfig, event, station, waveforms, crust_ind):
 
 
 def choose_backend(
-    fomosto_config, code, source_model, receiver_model, gf_directory="qseis2d_green"
+    fomosto_config,
+    code,
+    source_model,
+    receiver_model,
+    gf_directory="qseis2d_green",
+    version=None,
 ):
     """
     Get backend related config.
@@ -2178,7 +2183,8 @@ def choose_backend(
         raise TypeError('For near-field phases the "qseis" backend has to be used!')
 
     if code == "qseis":
-        version = "2006a"
+        default_version = "2006a"
+        code_version = version or default_version
         if "slowest" in waveforms or distances.min() < 10:
             logger.info(
                 "Receiver and source"
@@ -2204,16 +2210,17 @@ def choose_backend(
             wavelet_duration_samples=0.001,
             sw_flat_earth_transform=sw_flat_earth_transform,
             sw_algorithm=sw_algorithm,
-            qseis_version=version,
+            qseis_version=code_version,
         )
 
     elif code == "qssp":
         source_model = copy.deepcopy(receiver_model)
         receiver_model = None
-        version = "2010"
+        default_version = "2010"
+        code_version = version or default_version
 
         conf = qssp.QSSPConfig(
-            qssp_version=version,
+            qssp_version=code_version,
             slowness_max=float(num.max(slowness_taper)),
             toroidal_modes=True,
             spheroidal_modes=True,
@@ -2223,6 +2230,8 @@ def choose_backend(
     else:
         raise NotImplementedError("Backend not supported: %s" % code)
 
+    logger.info("Using modelling code: %s with version: %s", code, code_version)
+
     # fill remaining fomosto params
     fc.earthmodel_1d = source_model
     fc.earthmodel_receiver_1d = receiver_model
@@ -2322,7 +2331,12 @@ def seis_construct_gf(
             gf_directory = os.path.join(sf.store_superdir, "base_gfs_%i" % crust_ind)
 
             conf = choose_backend(
-                fomosto_config, sf.code, source_model, receiver_model, gf_directory
+                fomosto_config,
+                sf.code,
+                source_model,
+                receiver_model,
+                gf_directory,
+                version=sf.version,
             )
 
             fomosto_config.validate()
@@ -2463,7 +2477,8 @@ def geo_construct_gf(event, geodetic_config, crust_ind=0, execute=True, force=Fa
     """
     from pyrocko.fomosto import psgrn_pscmp as ppp
 
-    version = "2008a"
+    default_version = "2008a"
+
     gfc = geodetic_config.gf_config
 
     # extract source crustal profile and check for water layer
@@ -2477,9 +2492,9 @@ def geo_construct_gf(event, geodetic_config, crust_ind=0, execute=True, force=Fa
 
     c = ppp.PsGrnPsCmpConfig()
 
-    c.pscmp_config.version = version
+    c.pscmp_config.version = gfc.version or default_version
 
-    c.psgrn_config.version = version
+    c.psgrn_config.version = gfc.version or default_version
     c.psgrn_config.sampling_interval = gfc.sampling_interval
     c.psgrn_config.gf_depth_spacing = gfc.medium_depth_spacing
     c.psgrn_config.gf_distance_spacing = gfc.medium_distance_spacing
@@ -3797,6 +3812,131 @@ def seis_synthetics(
         raise TypeError("Outmode %s not supported!" % outmode)
 
 
+spatial_derivative_parameters = {"dn": "north_shift", "de": "east_shift", "dd": "depth"}
+
+
+def seis_derivative(
+    engine,
+    sources,
+    targets,
+    arrival_taper,
+    arrival_times,
+    wavename,
+    filterer,
+    h,
+    parameter,
+    stencil_order=3,
+    outmode="tapered_data",
+):
+    """
+    Calculate numerical derivative with respect to source or spatial parameter
+    Parameters
+    ----------
+    engine : :class:`pyrocko.gf.seismosizer.LocalEngine`
+    sources : list
+        containing :class:`pyrocko.gf.seismosizer.Source` Objects
+        reference source is the first in the list!!!
+    targets : list
+        containing :class:`pyrocko.gf.seismosizer.Target` Objects
+    arrival_taper : :class:`ArrivalTaper`
+    arrival_times : list or:class:`numpy.ndarray`
+        containing the start times [s] since 1st.January 1970 to start
+        tapering
+    wavename : string
+        of the tabulated phase that determines the phase arrival
+    filterer : :class:`Filterer`
+    h : float
+        distance for derivative calculation
+    parameter : str
+        parameter with respect to which the derivative
+        is being calculated e.g. 'strike', 'dip', 'depth'
+    stencil_order : int
+        order N of numerical stencil differentiation, available; 3 or 5
+    Returns
+    -------
+    :class:`num.array` ntargets x nsamples with the first derivative
+    """
+
+    ntargets = len(targets)
+
+    available_params = list(sources[0].keys()) + list(
+        spatial_derivative_parameters.keys()
+    )
+    if parameter not in available_params:
+        raise AttributeError(
+            "Parameter for which the derivative was requested is neither"
+            " represented by the source nor the target. Supported parameters:"
+            " %s" % utility.list2string(available_params)
+        )
+
+    calc_sources = copy.deepcopy(sources)
+    store = engine.get_store(targets[0].store_id)
+    nsamples = int(num.ceil(store.config.sample_rate * arrival_taper.duration()))
+
+    stencil = utility.StencilOperator(h=h, order=stencil_order)
+
+    # loop over stencil steps
+    n_stencil_steps = len(stencil)
+    tmp = num.zeros((n_stencil_steps, ntargets, nsamples), dtype="float64")
+    for i, hstep in enumerate(stencil.hsteps):
+
+        if parameter in spatial_derivative_parameters:
+            target_param_name = spatial_derivative_parameters[parameter]
+            diff_targets = []
+            diff_sources = sources
+            for target in targets:
+                target_diff = copy.deepcopy(target)
+                target_param = getattr(target, target_param_name)
+                setattr(target_diff, target_param_name, target_param + hstep)
+                diff_targets.append(target_diff)
+
+            arrival_times = num.repeat(arrival_times, n_stencil_steps)
+        else:
+            diff_targets = targets
+            diff_sources = []
+            for source in calc_sources:
+                source_diff = source.clone()
+                source_param = source[parameter]
+                setattr(source_diff, parameter, source_param + hstep)
+                diff_sources.append(source_diff)
+
+        tmp[i, :, :], tmins = seis_synthetics(
+            engine=engine,
+            sources=diff_sources,
+            targets=diff_targets,
+            arrival_taper=arrival_taper,
+            wavename=wavename,
+            filterer=filterer,
+            arrival_times=arrival_times,
+            pre_stack_cut=True,
+            outmode="array",
+            chop_bounds=["b", "c"],
+        )
+
+    diff_array = (tmp * stencil.coefficients).sum(axis=0) / stencil.denominator
+
+    if outmode == "array":
+        return diff_array
+    elif outmode == "tapered_data":
+        out_traces = []
+        for i, target in enumerate(targets):
+            store = engine.get_store(target.store_id)
+            network, station, location, channel = target.codes
+            strain_trace = trace.Trace(
+                tmin=tmins[i],
+                ydata=diff_array[i, :],
+                network=network,
+                station=station,
+                location="{}{}".format(parameter, location),
+                channel=channel,
+                deltat=store.config.deltat,
+            )
+            out_traces.append(strain_trace)
+        return out_traces
+    else:
+        raise IOError("Outmode %s not supported!" % outmode)
+
+
 def radiation_weights_p(takeoff_angles, azimuths):
     """
     Station dependent propagation coefficients for P waves
diff --git a/beat/plotting/common.py b/beat/plotting/common.py
index 281e5f10..7aadad8f 100644
--- a/beat/plotting/common.py
+++ b/beat/plotting/common.py
@@ -310,6 +310,8 @@ def histplot_op(
     """
 
     cumulative = kwargs.pop("cumulative", False)
+    nsources = kwargs.pop("nsources", False)
+    isource = kwargs.pop("isource", 0)
 
     if color is not None and cmap is not None:
         logger.debug("Using color for histogram edgecolor ...")
@@ -396,30 +398,39 @@ def histplot_op(
 
             for quantile, value in sigma_quants.items():
                 quantile /= 100.0
-                x = [leftb, value, value]
-                y = [quantile, quantile, 0.0]
+                if nsources == 1:
+                    x = [leftb, value, value]
+                    y = [quantile, quantile, 0.0]
+                else:
+                    x = [leftb, rightb]
+                    y = [quantile, quantile]
 
-                ax.plot(x, y, "--k", linewidth=0.5)
                 fontsize = 6
-                xval = (value - num.abs(leftb)) / 2 + leftb
-
-                ax.text(
-                    xval,
-                    quantile,
-                    "{}%".format(int(quantile * 100)),
-                    fontsize=fontsize,
-                    horizontalalignment="center",
-                    verticalalignment="bottom",
-                )
-                # if quantile > 0.3:
-                ax.text(
-                    value,
-                    quantile / 2,
-                    "%.3f" % value,
-                    fontsize=fontsize,
-                    horizontalalignment="left",
-                    verticalalignment="bottom",
-                )
+
+                if isource + 1 == nsources:
+                    # plot for last hist in axis
+                    ax.plot(x, y, "--k", linewidth=0.5)
+
+                    xval = (value - leftb) / 2 + leftb
+
+                    ax.text(
+                        xval,
+                        quantile,
+                        "{}%".format(int(quantile * 100)),
+                        fontsize=fontsize,
+                        horizontalalignment="center",
+                        verticalalignment="bottom",
+                    )
+
+                if nsources == 1:
+                    ax.text(
+                        value,
+                        quantile / 2,
+                        "%.3f" % value,
+                        fontsize=fontsize,
+                        horizontalalignment="left",
+                        verticalalignment="bottom",
+                    )
 
 
 def kde2plot_op(ax, x, y, grid=200, **kwargs):
diff --git a/beat/plotting/marginals.py b/beat/plotting/marginals.py
index 5c718dd5..330d8095 100644
--- a/beat/plotting/marginals.py
+++ b/beat/plotting/marginals.py
@@ -284,6 +284,7 @@ def remove_var(varnames, varname):
 
     figs = []
     fig_axs = []
+    var_idx = 0
     for nsubplots in nsubplots_page:
 
         width, height = mpl_papersize("a4", "portrait")
@@ -298,7 +299,7 @@ def remove_var(varnames, varname):
             coli, rowi = utility.mod_i(i, nrow)
             ax = axs[rowi, coli]
 
-            if i > n - 1:
+            if var_idx > n - 1:
                 try:
                     fig.delaxes(ax)
                 except KeyError:
@@ -307,7 +308,8 @@ def remove_var(varnames, varname):
                 if nvar == 1:
                     source_idxs = [backup_source_idxs[i]]
 
-                v = varnames[i]
+                v = varnames[var_idx]
+                var_idx += 1
 
                 color = copy.deepcopy(input_color)
 
@@ -320,7 +322,9 @@ def remove_var(varnames, varname):
                     if v in dist_vars:
                         if source_idxs is None:
                             source_idx_step = int(num.floor(d.shape[1] / 6))
-                            logger.info("No patches defined using 1 every %i!")
+                            logger.info(
+                                "No patches defined using 1 every %i!", source_idx_step
+                            )
                             source_idxs = num.arange(
                                 0, d.shape[1], source_idx_step
                             ).tolist()
@@ -329,13 +333,21 @@ def remove_var(varnames, varname):
                             "Plotting patches: %s" % utility.list2string(source_idxs)
                         )
 
-                        try:
-                            selected = num.atleast_2d(d.T[source_idxs])
-                        except IndexError:
-                            raise IndexError(
-                                "One or several patches do not exist! "
-                                "Patch idxs: %s" % utility.list2string(source_idxs)
-                            )
+                        selected = []
+                        for s_idx in source_idxs:
+                            try:
+                                if isinstance(s_idx, slice):
+                                    d_sel = num.atleast_2d(d.T[s_idx].mean(0))
+                                else:
+                                    d_sel = num.atleast_2d(d.T[s_idx])
+                            except IndexError:
+                                raise IndexError(
+                                    "One or several patches do not exist! "
+                                    "Patch idxs: %s" % utility.list2string([s_idx])
+                                )
+                            selected.append(d_sel)
+
+                        selected = num.vstack(selected)
                     else:
                         selected = d.T
 
@@ -384,6 +396,9 @@ def remove_var(varnames, varname):
                         elif plot_style in ["pdf", "cdf"]:
 
                             kwargs["label"] = source_idxs
+                            # following determine quantile annotations in cdf
+                            kwargs["nsources"] = nsources
+                            kwargs["isource"] = isource
                             if plot_style == "cdf":
                                 kwargs["cumulative"] = True
                             else:
@@ -409,8 +424,8 @@ def remove_var(varnames, varname):
 
                             if v in dist_vars:
                                 try:  # variable bounds
-                                    lower = param.lower[source_idxs]
-                                    upper = param.upper[source_idxs]
+                                    lower = param.lower[tuple(source_idxs)]
+                                    upper = param.upper[tuple(source_idxs)]
                                 except IndexError:
                                     lower, upper = param.lower, param.upper
 
diff --git a/beat/sources.py b/beat/sources.py
index 363ad3bc..a128190f 100644
--- a/beat/sources.py
+++ b/beat/sources.py
@@ -58,13 +58,6 @@ def dipvector(self):
         """
         Get 3 dimensional dip-vector of the planar fault.
 
-        Parameters
-        ----------
-        dip : scalar, float
-            dip-angle [deg] of the fault
-        strike : scalar, float
-            strike-abgle [deg] of the fault
-
         Returns
         -------
         :class:`numpy.ndarray`
@@ -83,11 +76,6 @@ def strikevector(self):
         """
         Get 3 dimensional strike-vector of the planar fault.
 
-        Parameters
-        ----------
-        strike : scalar, float
-            strike-abgle [deg] of the fault
-
         Returns
         -------
         :class:`numpy.ndarray`
@@ -95,6 +83,17 @@ def strikevector(self):
 
         return num.array([num.sin(self.strike * d2r), num.cos(self.strike * d2r), 0.0])
 
+    @property
+    def normalvector(self):
+        """
+        Get 3 dimensional normal-vector of the planar fault.
+        Returns
+        -------
+        :class:`numpy.ndarray`
+        """
+
+        return num.cross(self.strikevector, self.dipvector)
+
     @property
     def center(self):
         """
diff --git a/beat/utility.py b/beat/utility.py
index dd2531c6..bf93fa8c 100644
--- a/beat/utility.py
+++ b/beat/utility.py
@@ -20,6 +20,7 @@
 from pyrocko import catalog, orthodrome, util
 from pyrocko.cake import LayeredModel, m2d, read_nd_model_str
 from pyrocko.gf.seismosizer import RectangularSource
+from pyrocko.guts import Float, Int, Object
 from theano import config as tconfig
 
 logger = logging.getLogger("utility")
@@ -1093,6 +1094,23 @@ def running_window_rms(data, window_size, mode="valid"):
     return num.sqrt(num.convolve(data2, window, mode))
 
 
+def slice2string(slice_obj):
+    """
+    Wrapper for better formatted string method for slices.
+
+    Returns
+    -------
+    str
+    """
+    if isinstance(slice_obj, slice):
+        if slice_obj.step:
+            return "{}:{}:{}".format(slice_obj.start, slice_obj.stop, slice_obj.step)
+        else:
+            return "{}:{}".format(slice_obj.start, slice_obj.stop)
+    else:
+        return slice_obj
+
+
 def list2string(l, fill=", "):
     """
     Convert list of string to single string.
@@ -1102,7 +1120,20 @@ def list2string(l, fill=", "):
     l: list
         of strings
     """
-    return fill.join("%s" % listentry for listentry in l)
+    return fill.join("%s" % slice2string(listentry) for listentry in l)
+
+
+def string2slice(slice_string):
+    """
+    Convert string of slice form to python slice object.
+
+    Parameters
+    ----------
+    slice_string: str
+        of form "0:2" i.e. two integer numbers separated by colon
+    """
+
+    return slice(*[int(idx) for idx in slice_string.split(":")])
 
 
 def unique_list(l):
@@ -1556,3 +1587,42 @@ def find_elbow(data, theta=None, rotate_left=False):
     # rotate data vector
     rotated_data = data.dot(rotation_matrix)
     return rotated_data[:, 1].argmin(), rotated_data
+
+
+class StencilOperator(Object):
+
+    h = Float.T(default=0.1, help="step size left and right of the reference value")
+    order = Int.T(default=3, help="number of points of central differences")
+
+    def __init__(self, **kwargs):
+
+        stencil_order = kwargs["order"]
+        if stencil_order not in [3, 5]:
+            raise ValueError(
+                "Only stencil orders 3 and 5 implemented."
+                " Requested: %i" % stencil_order
+            )
+
+        self._coeffs = {3: num.array([1.0, -1.0]), 5: num.array([1.0, 8.0, -8.0, -1.0])}
+
+        self._denominator = {3: 2.0, 5: 12.0}
+
+        self._hsteps = {3: num.array([-1, 1]), 5: num.array([-2, -1, 1, 2])}
+
+        Object.__init__(self, **kwargs)
+
+    @property
+    def coefficients(self):
+        coeffs = self._coeffs[self.order]
+        return coeffs.reshape((coeffs.size, 1, 1))
+
+    def __len__(self):
+        return self.coefficients.size
+
+    @property
+    def denominator(self):
+        return self._denominator[self.order] * self.h
+
+    @property
+    def hsteps(self):
+        return self._hsteps[self.order] * self.h
diff --git a/docs/anaconda_installation.rst b/docs/anaconda_installation.rst
index 7e1cb64e..7470ad09 100644
--- a/docs/anaconda_installation.rst
+++ b/docs/anaconda_installation.rst
@@ -10,9 +10,9 @@ A general advice when dealing with anaconda is that the "sudo" command must NOT
 instead of the respective anaconda environment.
 Below are a series of commands that might be able to get you up and running using anaconda3 (thanks to Rebecca Salvage).
 
-Create and activate a new conda environment e.g. called "beat" using python3.6 (minimum required is 3.5)::
+Create and activate a new conda environment e.g. called "beat" using python3.8 (minimum required is 3.7)::
 
-  conda create -n beat python=3.6
+  conda create -n beat python=3.8
   conda activate beat
   cd ~/src  # or wherever you keep the packages
 
diff --git a/test/test_utility.py b/test/test_utility.py
index 9095ccce..4a7ba0d8 100644
--- a/test/test_utility.py
+++ b/test/test_utility.py
@@ -70,7 +70,7 @@ def test_list_ordering(self):
     def test_window_rms(self):
 
         data = num.random.randn(5000)
-        ws = data.size / 5
+        ws = int(data.size / 5)
         t0 = time()
         data_stds = utility.running_window_rms(data, window_size=ws)
         t1 = time()
@@ -87,6 +87,15 @@ def test_window_rms(self):
         data_stds = utility.running_window_rms(data, window_size=ws, mode="same")
         print(data_stds.shape)
 
+    def test_stencil(self):
+        for order in [3, 5]:
+            so = utility.StencilOperator(order=order, h=0.001)
+            print(so)
+            print(len(so))
+            print("hsteps", so.hsteps)
+            print("coeffs", so.coefficients)
+            print("denom", so.denominator)
+
 
 if __name__ == "__main__":
     util.setup_logging("test_utility", "warning")