diff --git a/hugo/content/installation/py/linux.md b/hugo/content/installation/py/linux.md
index e3c18e599d6a23bba78519283f9a3de2cb526a24..568bfc094db42d71436f0ade8374680c96d561c3 100644
--- a/hugo/content/installation/py/linux.md
+++ b/hugo/content/installation/py/linux.md
@@ -44,7 +44,7 @@ Then:
 # install pyenv
 curl https://pyenv.run | bash
 
-#install Python
+#install Python (e.g. {{% recommended-python %}})
 pyenv install {{% recommended-python %}}
 pyenv global {{% recommended-python %}}
 which python # shows path in virtual environment
diff --git a/rawEx/fit/specular/PolarizedSpinAsymmetry.py b/rawEx/fit/specular/PolarizedSpinAsymmetry.py
deleted file mode 120000
index 22a16a99f0b7beba4b7ea9f1b92c54ca83437e4f..0000000000000000000000000000000000000000
--- a/rawEx/fit/specular/PolarizedSpinAsymmetry.py
+++ /dev/null
@@ -1 +0,0 @@
-../../specular/PolarizedSpinAsymmetry.py
\ No newline at end of file
diff --git a/rawEx/fit/specular/PolarizedSpinAsymmetry.py b/rawEx/fit/specular/PolarizedSpinAsymmetry.py
new file mode 100755
index 0000000000000000000000000000000000000000..56f86fc5ba316d6036c4311c56146f930463d0e1
--- /dev/null
+++ b/rawEx/fit/specular/PolarizedSpinAsymmetry.py
@@ -0,0 +1,224 @@
+#!/usr/bin/env python3
+"""
+This simulation example demonstrates how to replicate the
+fitting example "Magnetically Dead Layers in Spinel Films"
+given at the Nist website:
+https://www.nist.gov/ncnr/magnetically-dead-layers-spinel-films
+
+For simplicity, here we only reproduce the first part of that
+demonstration without the magnetically dead layer.
+"""
+
+import os
+import numpy
+import matplotlib.pyplot as plt
+import bornagain as ba
+from bornagain import angstrom, ba_plot as bp, deg, R3
+
+# q-range on which the simulation and fitting are to be performed
+qmin = 0.05997
+qmax = 1.96
+
+# number of points on which the computed result is plotted
+scan_size = 1500
+
+# The SLD of the substrate is kept constant
+sldMao = (5.377e-06, 0)
+
+# constant to convert between B and magnetic SLD
+RhoMconst = 2.910429812376859e-12
+
+####################################################################
+#                  Create Sample and Simulation                    #
+####################################################################
+
+def get_sample(P):
+    """
+    construct the sample with the given parameters
+    """
+    BMagnitude = P["rhoM_Mafo"]*1e-6/RhoMconst
+    angle = 0
+    B = R3(
+        BMagnitude*numpy.sin(angle*deg),
+        BMagnitude*numpy.cos(angle*deg), 0)
+
+    vacuum = ba.MaterialBySLD("Vacuum", 0, 0)
+    material_layer = ba.MaterialBySLD("(Mg,Al,Fe)3O4", P["rho_Mafo"]*1e-6, 0, B)
+    material_substrate = ba.MaterialBySLD("MgAl2O4", *sldMao)
+
+    ambient_layer = ba.Layer(vacuum)
+    layer = ba.Layer(material_layer, P["t_Mafo"]*angstrom)
+    substrate_layer = ba.Layer(material_substrate)
+
+    r_Mafo = ba.LayerRoughness(P["r_Mafo"]*angstrom)
+    r_substrate = ba.LayerRoughness(P["r_Mao"]*angstrom)
+
+    sample = ba.MultiLayer()
+    sample.addLayer(ambient_layer)
+    sample.addLayerWithTopRoughness(layer, r_Mafo)
+    sample.addLayerWithTopRoughness(substrate_layer, r_substrate)
+
+    return sample
+
+
+def get_simulation(sample, q_axis, parameters, polarizer_vec,
+                   analyzer_vec):
+    """
+    A simulation object.
+    Polarization, analyzer and resolution are set
+    from given parameters
+    """
+    q_axis = q_axis + parameters["q_offset"]
+    distr = ba.DistributionGaussian(0., 1., 25, 4.)
+
+    scan = ba.QzScan(q_axis)
+    scan.setAbsoluteQResolution(distr, parameters["q_res"])
+                                       # TODO CHECK not parameters["q_res"]*q_axis ??
+
+    scan.setPolarization(polarizer_vec)
+    scan.setAnalyzer(analyzer_vec)
+
+    return ba.SpecularSimulation(scan, sample)
+
+
+def run_simulation(q_axis, fitP, *, polarizer_vec, analyzer_vec):
+    """
+    Run a simulation on the given q-axis, where the sample is
+    constructed with the given parameters.
+    Vectors for polarization and analyzer need to be provided
+    """
+    parameters = dict(fitP, **fixedP)
+
+    sample = get_sample(parameters)
+    simulation = get_simulation(sample, q_axis, parameters, polarizer_vec,
+                                analyzer_vec)
+
+    return simulation.simulate()
+
+
+def qr(result):
+    """
+    Returns two arrays that hold the q-values as well as the
+    reflectivity from a given simulation result
+    """
+    q = numpy.array(result.axis(0).binCenters())
+    r = numpy.array(result.npArray())
+
+    return q, r
+
+####################################################################
+#                         Plot Handling                            #
+####################################################################
+
+def plotData(qs, rs, exps, labels, colors):
+    """
+    Plot the simulated result together with the experimental data
+    """
+    fig = plt.figure()
+    ax = fig.add_subplot(111)
+
+    for q, r, exp, l, c in zip(qs, rs, exps, labels, colors):
+        ax.errorbar(exp.xAxis().binCenters(),
+                    exp.flatVector(),
+                    # xerr=TODO i742,
+                    yerr=exp.errorSigmas(),
+                    fmt='.',
+                    markersize=0.75,
+                    linewidth=0.5,
+                    color=c[1])
+        ax.plot(q, r, label=l, color=c[0])
+
+    ax.set_yscale('log')
+    plt.legend()
+
+    plt.xlabel(r"$q$ (nm$^{-1}$)")
+    plt.ylabel("$R$")
+
+    plt.tight_layout()
+    <%= sm ? "plt.close()" : "" %>
+
+def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
+    """
+    Plot the simulated spin asymmetry as well its
+    experimental counterpart with errorbars
+    """
+    Yp = data_pp.npArray()
+    Ym = data_mm.npArray()
+    Ep = data_pp.npErrors()
+    Em = data_mm.npErrors()
+    # compute the errorbars of the spin asymmetry
+    delta = numpy.sqrt(4 * (Yp**2 * Em**2 + Ym**2 * Ep**2 ) / ( Yp + Ym )**4 )
+
+    fig = plt.figure()
+    ax = fig.add_subplot(111)
+
+    ax.errorbar(data_pp.xAxis().binCenters(),
+                (Yp - Ym) / (Yp + Ym),
+                # xerr=TODO i742,
+                yerr=delta,
+                fmt='.',
+                markersize=0.75,
+                linewidth=0.5)
+
+    ax.plot(q, (r_pp - r_mm)/(r_pp + r_mm))
+
+    plt.gca().set_ylim((-0.3, 0.5))
+
+    plt.xlabel(r"$q$ (nm$^{-1}$)")
+    plt.ylabel("Spin asymmetry")
+
+    plt.tight_layout()
+    <%= sm ? "plt.close()" : "" %>
+
+####################################################################
+#                          Data Handling                           #
+####################################################################
+
+def load_data(fname):
+    flags = ba.ImportSettings1D("q (1/nm)", "", "", 1, 2, 3, 4)
+    return ba.readData1D(fname, ba.csv1D, flags)
+
+####################################################################
+#                          Main Function                           #
+####################################################################
+
+if __name__ == '__main__':
+    datadir = os.getenv('BA_DATA_DIR', '')
+    fname_stem = os.path.join(datadir, "specular/MAFO_Saturated_")
+
+    expdata_pp = load_data(fname_stem + "pp.tab")
+    expdata_mm = load_data(fname_stem + "mm.tab")
+
+    fixedP = {
+        # parameters from our own fit run
+        'q_res': 0.010542945012551425,
+        'q_offset': 7.971243487467318e-05,
+        'rho_Mafo': 6.370140108715461,
+        'rhoM_Mafo': 0.27399566816062926,
+        't_Mafo': 137.46913056084736,
+        'r_Mao': 8.60487712674644,
+        'r_Mafo': 3.7844265311293483
+    }
+
+    def run_Simulation_pp(qzs, P):
+        return run_simulation(qzs,
+                              P,
+                              polarizer_vec=R3(0, 1, 0),
+                              analyzer_vec=R3(0, 1, 0))
+
+    def run_Simulation_mm(qzs, P):
+        return run_simulation(qzs,
+                              P,
+                              polarizer_vec=R3(0, -1, 0),
+                              analyzer_vec=R3(0, -1, 0))
+
+    qzs = numpy.linspace(qmin, qmax, scan_size)
+    q_pp, r_pp = qr(run_Simulation_pp(qzs, fixedP))
+    q_mm, r_mm = qr(run_Simulation_mm(qzs, fixedP))
+
+    plotData([q_pp, q_mm], [r_pp, r_mm], [expdata_pp, expdata_mm],
+             ["$++$", "$--$"], [['orange','red'], ['green','blue']])
+
+    plotSpinAsymmetry(expdata_pp, expdata_mm, qzs, r_pp, r_mm)
+
+    bp.show_or_export()