diff --git a/Core/Export/PythonFormatting.cpp b/Core/Export/PythonFormatting.cpp
index c3bb6dea7d3894ae08f4e70386ddfbed3a21fe3a..3db9e850473808cad016b1707a81484981ff9097 100644
--- a/Core/Export/PythonFormatting.cpp
+++ b/Core/Export/PythonFormatting.cpp
@@ -28,11 +28,27 @@
 #include "Rectangle.h"
 #include "StringUtils.h"
 #include "Units.h"
+#include "FixedBinAxis.h"
 #include <iomanip>
 
 namespace PythonFormatting
 {
 
+std::string scriptPreamble()
+{
+    const std::string result = "import numpy\n"
+                               "import bornagain as ba\n"
+                               "from bornagain import deg, angstrom, nm, kvector_t\n\n\n";
+
+    return result;
+}
+
+std::string getSampleFunctionName()
+{
+    return "get_sample";
+}
+
+
 //! Returns fixed Python code snippet that defines the function "runSimulation".
 
 std::string representShape2D(const std::string& indent, const IShape2D* ishape, bool mask_value,
@@ -160,6 +176,13 @@ std::string printValue(double value, const std::string& units)
                                  + "'");
 }
 
+std::string printString(const std::string& value)
+{
+    std::ostringstream result;
+    result << "\"" << value << "\"";
+    return result.str();
+}
+
 bool isSquare(double length1, double length2, double angle)
 {
     return length1 == length2 && Numeric::areAlmostEqual(angle, M_PI_2);
@@ -278,4 +301,22 @@ std::string printParameterDistribution(const ParameterDistribution& par_distr,
     return result.str();
 }
 
+std::string printAxis(const IAxis& axis, const std::string& units)
+{
+    std::ostringstream result;
+
+    if (auto fixedAxis = dynamic_cast<const FixedBinAxis*>(&axis)) {
+        result << "ba.FixedBinAxis("
+               << printString(fixedAxis->getName()) << ", "
+               << fixedAxis->size() << ", "
+               << printValue(fixedAxis->getMin(), units) << ", "
+               << printValue(fixedAxis->getMax(), units) << ")";
+
+    } else {
+        throw std::runtime_error("PythonFormatting::printAxis() -> Error. Unsupported axis");
+    }
+
+    return result.str();
+}
+
 } // namespace PythonFormatting
diff --git a/Core/Export/PythonFormatting.h b/Core/Export/PythonFormatting.h
index 5c5fcfdfbed65cffa3e24b30789abe51d5b164b2..df9879af061370c51e3cc009f7481cd8344e65a4 100644
--- a/Core/Export/PythonFormatting.h
+++ b/Core/Export/PythonFormatting.h
@@ -25,21 +25,28 @@ class IShape2D;
 class RealParameter;
 class ParameterDistribution;
 class RealLimits;
+class IAxis;
 
 //! Utility functions for writing Python code snippets.
 
 namespace PythonFormatting
 {
 
+BA_CORE_API_ std::string scriptPreamble();
+BA_CORE_API_ std::string getSampleFunctionName();
+
 BA_CORE_API_ std::string representShape2D(const std::string& indent, const IShape2D* ishape,
                                           bool mask_value,
                                           std::function<std::string(double)> printValueFunc);
+
+BA_CORE_API_ std::string printInt(int value);
 BA_CORE_API_ std::string printBool(double value);
 BA_CORE_API_ std::string printDouble(double input);
 BA_CORE_API_ std::string printNm(double input);
 BA_CORE_API_ std::string printScientificDouble(double input);
 BA_CORE_API_ std::string printDegrees(double input);
 BA_CORE_API_ std::string printValue(double value, const std::string& units = "");
+BA_CORE_API_ std::string printString(const std::string& value);
 
 BA_CORE_API_ bool isSquare(double length1, double length2, double angle);
 BA_CORE_API_ bool isHexagonal(double length1, double length2, double angle);
@@ -58,6 +65,8 @@ BA_CORE_API_ std::string printRealLimitsArg(const RealLimits& limits,
 BA_CORE_API_ std::string printParameterDistribution(const ParameterDistribution& par_distr,
                                                     const std::string& distVarName,
                                                     const std::string& units = "");
+
+BA_CORE_API_ std::string printAxis(const IAxis& axis, const std::string& units = "");
 }
 
 #endif // PYTHONFORMATTING_H
diff --git a/Core/Export/SampleToPython.cpp b/Core/Export/SampleToPython.cpp
index b750642052371f49a0a8d5a6a4a43e5ebdc0e517..8ba52c502db7319624616e0e6ede5655a3a06be9 100644
--- a/Core/Export/SampleToPython.cpp
+++ b/Core/Export/SampleToPython.cpp
@@ -86,11 +86,11 @@ SampleToPython::~SampleToPython() = default;
 
 std::string SampleToPython::defineGetSample() const
 {
-    return "def getSample():\n" + defineMaterials() + defineLayers() + defineFormFactors()
+    return "def "+getSampleFunctionName()+"():\n" + defineMaterials() + defineLayers() + defineFormFactors()
            + defineParticles() + defineCoreShellParticles() + defineParticleCompositions()
            + defineLattices() + defineCrystals() + defineMesoCrystals()
            + defineParticleDistributions() + defineInterferenceFunctions() + defineParticleLayouts()
-           + defineRoughnesses() + addLayoutsToLayers() + defineMultiLayers() + "\n";
+           + defineRoughnesses() + addLayoutsToLayers() + defineMultiLayers() + "\n\n";
 }
 
 const std::map<MATERIAL_TYPES, std::string> factory_names{
diff --git a/Core/Export/SimulationToPython.cpp b/Core/Export/SimulationToPython.cpp
index 245d3ad574f7ec5341e04b2e03e23f1bcf03bc05..1c6bf7ce4b5a8825c254d5c2bc5f9bd23fbabbb3 100644
--- a/Core/Export/SimulationToPython.cpp
+++ b/Core/Export/SimulationToPython.cpp
@@ -25,20 +25,16 @@
 #include "ResolutionFunction2DGaussian.h"
 #include "SampleToPython.h"
 #include "SphericalDetector.h"
+#include "OffSpecSimulation.h"
 #include <iomanip>
 
 using namespace PythonFormatting;
 
 namespace
 {
-const std::string preamble = "import numpy\n"
-                             "import bornagain as ba\n"
-                             "from bornagain import deg, angstrom, nm, kvector_t\n\n";
-
 const std::string defineSimulate = "def run_simulation():\n"
-                                   "    # Run Simulation\n"
-                                   "    sample = getSample()\n"
-                                   "    simulation = getSimulation()\n"
+                                   "    sample = "+getSampleFunctionName()+"()\n"
+                                   "    simulation = get_simulation()\n"
                                    "    simulation.setSample(sample)\n"
                                    "    simulation.runSimulation()\n"
                                    "    return simulation.result()\n"
@@ -52,7 +48,7 @@ std::function<std::string(double)> printFunc(const IDetector* detector)
     if (detector->defaultAxesUnits() == AxesUnits::RADIANS)
         return PythonFormatting::printDegrees;
     throw Exceptions::RuntimeErrorException(
-        "ExportToPython::defineMasks() -> Error. Unknown detector units.");
+        "SimulationToPython::defineMasks() -> Error. Unknown detector units.");
 }
 } // namespace
 
@@ -62,41 +58,59 @@ std::string SimulationToPython::generateSimulationCode(const Simulation& simulat
                                                        EMainType mainType)
 {
     if (simulation.sample() == nullptr)
-        throw std::runtime_error("ExportToPython::generateSimulationCode() -> Error. "
+        throw std::runtime_error("SimulationToPython::generateSimulationCode() -> Error. "
                                  "Simulation is not initialized.");
 
     SampleToPython sampleGenerator;
 
-    return preamble + sampleGenerator.generateSampleCode(*simulation.sample())
+    return scriptPreamble() + sampleGenerator.generateSampleCode(*simulation.sample())
            + defineGetSimulation(&simulation) + defineSimulate + defineMain(mainType);
 }
 
 std::string SimulationToPython::defineGetSimulation(const Simulation* simulation) const
 {
     std::ostringstream result;
+    result << "def get_simulation():\n";
+
     if (auto gisas = dynamic_cast<const GISASSimulation*>(simulation))
         result << defineGISASSimulation(gisas);
+    else if (auto offspec = dynamic_cast<const OffSpecSimulation*>(simulation))
+        result << defineOffSpecSimulation(offspec);
     else
-        throw std::runtime_error("ExportToPython::defineGetSimulation() -> Error. Wrong simulation "
-                                 "type");
+        throw std::runtime_error("SimulationToPython::defineGetSimulation() -> Error. "
+                                 "Wrong simulation type");
 
+    result << indent() << "return simulation\n\n\n";
     return result.str();
 }
 
 std::string SimulationToPython::defineGISASSimulation(const GISASSimulation* simulation) const
 {
     std::ostringstream result;
-    result << "def getSimulation():\n";
     result << indent() << "simulation = ba.GISASSimulation()\n";
     result << defineDetector(simulation);
     result << defineDetectorResolutionFunction(simulation);
     result << defineDetectorPolarizationAnalysis(simulation);
-    result << defineBeam(simulation);
+    result << defineGISASBeam(*simulation);
+    result << defineParameterDistributions(simulation);
+    result << defineMasks(simulation);
+    result << defineSimulationOptions(simulation);
+    result << defineBackground(simulation);
+    return result.str();
+}
+
+std::string SimulationToPython::defineOffSpecSimulation(const OffSpecSimulation* simulation) const
+{
+    std::ostringstream result;
+    result << indent() << "simulation = ba.OffSpecSimulation()\n";
+    result << defineDetector(simulation);
+    result << defineDetectorResolutionFunction(simulation);
+    result << defineDetectorPolarizationAnalysis(simulation);
+    result << defineOffSpecBeam(*simulation);
     result << defineParameterDistributions(simulation);
     result << defineMasks(simulation);
     result << defineSimulationOptions(simulation);
     result << defineBackground(simulation);
-    result << indent() << "return simulation\n\n\n";
     return result.str();
 }
 
@@ -105,7 +119,7 @@ std::string SimulationToPython::defineDetector(const Simulation* simulation) con
     const IDetector* iDetector = simulation->getInstrument().getDetector();
 
     if (iDetector->dimension() != 2)
-        throw Exceptions::RuntimeErrorException("ExportToPython::defineDetector: "
+        throw Exceptions::RuntimeErrorException("SimulationToPython::defineDetector: "
                                                 "detector must be two-dimensional for GISAS");
     std::ostringstream result;
     result << std::setprecision(12);
@@ -156,11 +170,12 @@ std::string SimulationToPython::defineDetector(const Simulation* simulation) con
                    << printDouble(detector->getDirectBeamV0()) << ")\n";
         } else
             throw Exceptions::RuntimeErrorException(
-                "ExportToPython::defineDetector: unknown alignment");
+                "SimulationToPython::defineDetector() -> Error. Unknown alignment.");
 
         result << indent() << "simulation.setDetector(detector)\n";
     } else
-        throw Exceptions::RuntimeErrorException("ExportToPython::defineDetector: unknown detector");
+        throw Exceptions::RuntimeErrorException("SimulationToPython::defineDetector() -> Error. "
+                                                "Unknown detector");
     if (iDetector->regionOfInterest()) {
         result << indent() << "simulation.setRegionOfInterest("
                << printFunc(iDetector)(iDetector->regionOfInterest()->getXlow()) << ", "
@@ -187,11 +202,11 @@ std::string SimulationToPython::defineDetectorResolutionFunction(const Simulatio
                 result << printFunc(detector)(resfunc->getSigmaY()) << "))\n";
             } else
                 throw Exceptions::RuntimeErrorException(
-                    "ExportToPython::defineDetectorResolutionFunction() -> Error. "
+                    "SimulationToPython::defineDetectorResolutionFunction() -> Error. "
                     "Unknown detector resolution function");
         } else
             throw Exceptions::RuntimeErrorException(
-                "ExportToPython::defineDetectorResolutionFunction() -> Error. "
+                "SimulationToPython::defineDetectorResolutionFunction() -> Error. "
                 "Not a ConvolutionDetectorResolution function");
     }
     return result.str();
@@ -219,14 +234,39 @@ SimulationToPython::defineDetectorPolarizationAnalysis(const Simulation* simulat
     return result.str();
 }
 
-std::string SimulationToPython::defineBeam(const Simulation* simulation) const
+std::string SimulationToPython::defineGISASBeam(const GISASSimulation& simulation) const
 {
     std::ostringstream result;
-    result << std::setprecision(12);
-    // result << indent() << "# Defining Beam Parameters\n";
-    const Beam& beam = simulation->getInstrument().getBeam();
+    const Beam& beam = simulation.getInstrument().getBeam();
+
     result << indent() << "simulation.setBeamParameters(" << printNm(beam.getWavelength()) << ", "
            << printDegrees(beam.getAlpha()) << ", " << printDegrees(beam.getPhi()) << ")\n";
+
+    result << defineBeamPolarization(beam);
+    result << defineBeamIntensity(beam);
+
+    return result.str();
+}
+
+std::string SimulationToPython::defineOffSpecBeam(const OffSpecSimulation& simulation) const
+{
+    std::ostringstream result;
+    const Beam& beam = simulation.getInstrument().getBeam();
+
+    result << indent() << "alpha_i_axis = "
+           << PythonFormatting::printAxis(*simulation.beamAxis(), BornAgain::UnitsRad) << "\n";
+
+    result << indent() << "simulation.setBeamParameters(" << printNm(beam.getWavelength()) << ", "
+           << "alpha_i_axis, " << printDegrees(beam.getPhi()) << ")\n";
+
+    result << defineBeamPolarization(beam);
+    result << defineBeamIntensity(beam);
+    return result.str();
+}
+
+std::string SimulationToPython::defineBeamPolarization(const Beam& beam) const
+{
+    std::ostringstream result;
     auto bloch_vector = beam.getBlochVector();
     if (bloch_vector.mag() > 0.0) {
         std::string beam_polarization = "beam_polarization";
@@ -235,6 +275,12 @@ std::string SimulationToPython::defineBeam(const Simulation* simulation) const
                << ")\n";
         result << indent() << "simulation.setBeamPolarization(" << beam_polarization << ")\n";
     }
+    return result.str();
+}
+
+std::string SimulationToPython::defineBeamIntensity(const Beam& beam) const
+{
+    std::ostringstream result;
     double beam_intensity = beam.getIntensity();
     if (beam_intensity > 0.0)
         result << indent() << "simulation.setBeamIntensity("
@@ -242,6 +288,7 @@ std::string SimulationToPython::defineBeam(const Simulation* simulation) const
     return result.str();
 }
 
+
 std::string SimulationToPython::defineParameterDistributions(const Simulation* simulation) const
 {
     std::ostringstream result;
@@ -339,7 +386,7 @@ std::string SimulationToPython::defineMain(SimulationToPython::EMainType mainTyp
                  "        exit(\"File name is required\")\n"
                  "    ba.IntensityDataIOFactory.writeSimulationResult(result, sys.argv[1])\n";
     } else {
-        throw std::runtime_error("ExportToPython::defineMain() -> Error. Unknown main type.");
+        throw std::runtime_error("SimulationToPython::defineMain() -> Error. Unknown main type.");
     }
     return result;
 }
diff --git a/Core/Export/SimulationToPython.h b/Core/Export/SimulationToPython.h
index a22713185e12610c3f8933d6d49f01722aa9acb6..a28751a307dcd4e161a4c80a0910f8b8bad78571 100644
--- a/Core/Export/SimulationToPython.h
+++ b/Core/Export/SimulationToPython.h
@@ -21,6 +21,8 @@
 
 class Simulation;
 class GISASSimulation;
+class OffSpecSimulation;
+class Beam;
 
 //! Write a Python script that allows to run the current simulation.
 
@@ -38,10 +40,17 @@ private:
     std::string definePreamble() const;
     std::string defineGetSimulation(const Simulation* simulation) const;
     std::string defineGISASSimulation(const GISASSimulation* simulation) const;
+    std::string defineOffSpecSimulation(const OffSpecSimulation* simulation) const;
     std::string defineDetector(const Simulation* simulation) const;
     std::string defineDetectorResolutionFunction(const Simulation* simulation) const;
     std::string defineDetectorPolarizationAnalysis(const Simulation* simulation) const;
-    std::string defineBeam(const Simulation* simulation) const;
+
+    std::string defineGISASBeam(const GISASSimulation& simulation) const;
+    std::string defineOffSpecBeam(const OffSpecSimulation& simulation) const;
+
+    std::string defineBeamPolarization(const Beam& beam) const;
+    std::string defineBeamIntensity(const Beam& beam) const;
+
     std::string defineParameterDistributions(const Simulation* simulation) const;
     std::string defineMasks(const Simulation* simulation) const;
     std::string defineSimulationOptions(const Simulation* simulation) const;
diff --git a/Core/Simulation/OffSpecSimulation.cpp b/Core/Simulation/OffSpecSimulation.cpp
index f92a247382e0bb91c6f16c3f1ee0430c0e60b38c..c5a6e697a73d972f8e37bc6ab03c6c2b34e0d36b 100644
--- a/Core/Simulation/OffSpecSimulation.cpp
+++ b/Core/Simulation/OffSpecSimulation.cpp
@@ -73,6 +73,11 @@ void OffSpecSimulation::setBeamParameters(double wavelength, const IAxis& alpha_
     updateIntensityMap();
 }
 
+const IAxis* OffSpecSimulation::beamAxis() const
+{
+    return mP_alpha_i_axis.get();
+}
+
 OffSpecSimulation::OffSpecSimulation(const OffSpecSimulation& other)
     : Simulation2D(other)
 {
diff --git a/Core/Simulation/OffSpecSimulation.h b/Core/Simulation/OffSpecSimulation.h
index 61e17033b9305676c4911a8a5668d2869cf05ec0..86a7c9b7ab0c1418ea293736f7a561ccf165fc0e 100644
--- a/Core/Simulation/OffSpecSimulation.h
+++ b/Core/Simulation/OffSpecSimulation.h
@@ -49,6 +49,9 @@ public:
     //! Sets beam parameters from here (forwarded to Instrument)
     void setBeamParameters(double wavelength, const IAxis& alpha_axis, double phi_i);
 
+    //! Returns axis of the beam.
+    const IAxis* beamAxis() const;
+
 private:
     OffSpecSimulation(const OffSpecSimulation& other);
 
diff --git a/Core/StandardSamples/ResonatorBuilder.cpp b/Core/StandardSamples/ResonatorBuilder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..31526bc1d6f53d320a9adc8485c73e34047ccdb6
--- /dev/null
+++ b/Core/StandardSamples/ResonatorBuilder.cpp
@@ -0,0 +1,60 @@
+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      Core/StandardSamples/ResonatorBuilder.cpp
+//! @brief     Implements ResonatorBuilder class.
+//!
+//! @homepage  http://www.bornagainproject.org
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************** //
+
+#include "ResonatorBuilder.h"
+#include "Layer.h"
+#include "LayerRoughness.h"
+#include "MaterialFactoryFuncs.h"
+#include "MultiLayer.h"
+#include "Units.h"
+#include <memory>
+
+MultiLayer* ResonatorBuilder::buildSample() const
+{
+    std::unique_ptr<MultiLayer> result(new MultiLayer);
+
+    Material m_Si = HomogeneousMaterial("Si", 8.25218379931e-06, 0.0);
+    Material m_Ti = HomogeneousMaterial("Ti", -7.6593316363e-06, 3.81961616312e-09);
+    Material m_TiO2 = HomogeneousMaterial("TiO2", 1.04803530026e-05, 2.03233519385e-09);
+    Material m_Pt = HomogeneousMaterial("Pt", 2.52936993309e-05, 7.54553992473e-09);
+    Material m_D2O = HomogeneousMaterial("D2O", 2.52897204573e-05, 4.5224432814e-13);
+
+    Layer l_TiO2(m_TiO2, 3.0 * Units::nm);
+    Layer l_Ti_top(m_Ti, 10.0 * Units::nm);
+    Layer l_Ti(m_Ti, 13.0 * Units::nm);
+    Layer l_Si(m_Si);
+    Layer l_Pt(m_Pt, 32.0 * Units::nm);
+    Layer l_D2O(m_D2O);
+
+    LayerRoughness roughness;
+    roughness.setSigma(2.0 * Units::nm);
+    roughness.setHurstParameter(0.8);
+    roughness.setLatteralCorrLength(10.0 * Units::micrometer);
+
+    result->addLayer(l_Si);
+
+    const int nlayers = 3;
+    for (size_t i = 0; i < nlayers; ++i) {
+        result->addLayerWithTopRoughness(l_Ti, roughness);
+        result->addLayerWithTopRoughness(l_Pt, roughness);
+    }
+
+    result->addLayerWithTopRoughness(l_Ti_top, roughness);
+    result->addLayerWithTopRoughness(l_TiO2, roughness);
+    result->addLayerWithTopRoughness(l_D2O, roughness);
+
+    result->setCrossCorrLength(400 * Units::nm);
+
+    return result.release();
+}
diff --git a/Core/StandardSamples/ResonatorBuilder.h b/Core/StandardSamples/ResonatorBuilder.h
new file mode 100644
index 0000000000000000000000000000000000000000..9330f7289038c5e7e5b00c23d4f03ab53e58c70c
--- /dev/null
+++ b/Core/StandardSamples/ResonatorBuilder.h
@@ -0,0 +1,30 @@
+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      Core/StandardSamples/ResonatorBuilder.h
+//! @brief     Defines ResonatorBuilder class.
+//!
+//! @homepage  http://www.bornagainproject.org
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************** //
+
+#ifndef RESONATORBUILDER_H
+#define RESONATORBUILDER_H
+
+#include "IMultiLayerBuilder.h"
+
+//! Builds sample: multilayer with Ti/Pt layers sequence.
+//! @ingroup standard_samples
+
+class BA_CORE_API_ ResonatorBuilder : public IMultiLayerBuilder
+{
+public:
+    ResonatorBuilder() = default;
+    MultiLayer* buildSample() const;
+};
+
+#endif  // RESONATORBUILDER_H
diff --git a/Core/StandardSamples/SampleBuilderFactory.cpp b/Core/StandardSamples/SampleBuilderFactory.cpp
index a4f028ba4578e6ef6f905505e5923c8a8e903f4c..f65230e57d57a02cfa064ab8930b3b477471bcd2 100644
--- a/Core/StandardSamples/SampleBuilderFactory.cpp
+++ b/Core/StandardSamples/SampleBuilderFactory.cpp
@@ -40,6 +40,7 @@
 #include "TwoDimLatticeBuilder.h"
 #include "TwoLayerRoughnessBuilder.h"
 #include "SlicedParticleBuilder.h"
+#include "ResonatorBuilder.h"
 
 SampleBuilderFactory::SampleBuilderFactory()
 {
@@ -300,6 +301,12 @@ SampleBuilderFactory::SampleBuilderFactory()
         "MaterialBySLDBuilder",
         create_new<MaterialBySLDBuilder>,
         "Alternating homogeneous layers of Ti and Ni on silicone substrate (wavelength-independent).");
+
+    registerItem(
+        "ResonatorBuilder",
+        create_new<ResonatorBuilder>,
+        "Multilayer with bi-layer sequence for OffSpec testing.");
+
 }
 
 //! Retrieves a SampleBuilder from the registry, does the build, and returns the result.
diff --git a/Core/StandardSamples/SimulationFactory.cpp b/Core/StandardSamples/SimulationFactory.cpp
index 26f5ec460bc0df305e07e6dce6191d1c56b7859b..c93ef9e43a50f92d2d60c9f7f41752ca0eac80bf 100644
--- a/Core/StandardSamples/SimulationFactory.cpp
+++ b/Core/StandardSamples/SimulationFactory.cpp
@@ -17,6 +17,7 @@
 #include "SpecularSimulation.h"
 #include "RealParameter.h"
 #include "StandardSimulations.h"
+#include "OffSpecSimulation.h"
 
 SimulationFactory::SimulationFactory()
 {
@@ -147,4 +148,7 @@ SimulationFactory::SimulationFactory()
     registerItem("SpecularDivergentBeam", StandardSimulations::SpecularDivergentBeam,
                  "Simulation implies beam divergence both in wavelength and "
                  "inclination angle.");
+
+    registerItem("OffSpecMini", StandardSimulations::MiniOffSpec,
+                 "Mini OffSpecular simulation for resonator experiment.");
 }
diff --git a/Core/StandardSamples/StandardSimulations.cpp b/Core/StandardSamples/StandardSimulations.cpp
index 7c053629c05996466f145010b5f4af0e3c89b525..7eef2cd19e913b4db9f06f0d84f3158a43f66f2e 100644
--- a/Core/StandardSamples/StandardSimulations.cpp
+++ b/Core/StandardSamples/StandardSimulations.cpp
@@ -30,6 +30,7 @@
 #include "ResolutionFunction2DGaussian.h"
 #include "SampleBuilderFactory.h"
 #include "SpecularSimulation.h"
+#include "OffSpecSimulation.h"
 #include "Units.h"
 #include <memory>
 
@@ -441,3 +442,30 @@ SpecularSimulation* StandardSimulations::SpecularDivergentBeam()
 
     return result.release();
 }
+
+// OffSpec simulation used in ResonatorOffSpecSetup.py
+OffSpecSimulation* StandardSimulations::MiniOffSpec()
+{
+    std::unique_ptr<OffSpecSimulation> result(new OffSpecSimulation());
+
+    const int n_alpha(19);
+    const double alpha_min(0.0*Units::deg);
+    const double alpha_max(4.0*Units::deg);
+    const int n_phi(9);
+    const double phi_min(-0.1*Units::deg);
+    const double phi_max(0.1*Units::deg);
+
+    result->setDetectorParameters(n_phi, phi_min, phi_max, n_alpha, alpha_min, alpha_max);
+
+    const int n_scan_points(n_alpha);
+    const double alpha_i_min(alpha_min);
+    const double alpha_i_max(alpha_max);
+
+    FixedBinAxis alpha_i_axis("alpha_i", n_scan_points, alpha_i_min, alpha_i_max);
+    result->setBeamParameters(5.0*Units::angstrom, alpha_i_axis, 0.0);
+
+    result->setBeamIntensity(1e9);
+    result->getOptions().setIncludeSpecular(true);
+
+    return result.release();
+}
diff --git a/Core/StandardSamples/StandardSimulations.h b/Core/StandardSamples/StandardSimulations.h
index 996bc2a970872cc68411b000b004aeec66efdc37..ba8f98c7fc4448b38f5f424a8780a22ee5e0351b 100644
--- a/Core/StandardSamples/StandardSimulations.h
+++ b/Core/StandardSamples/StandardSimulations.h
@@ -17,6 +17,7 @@
 
 class GISASSimulation;
 class SpecularSimulation;
+class OffSpecSimulation;
 
 //! Standard pre-defined simulations.
 
@@ -57,6 +58,10 @@ SpecularSimulation* SpecularWithGaussianBeam();
 SpecularSimulation* SpecularWithSquareBeam();
 SpecularSimulation* SpecularDivergentBeam();
 
+// OffSpec simulations
+OffSpecSimulation* MiniOffSpec();
+
+
 } // namespace StandardSimulations
 
 #endif // STANDARDSIMULATIONS_H
diff --git a/Examples/python/simulation/ex05_BeamAndDetector/OffSpecularSimulation.py b/Examples/python/simulation/ex05_BeamAndDetector/OffSpecularSimulation.py
index 8467e81a8ef3c89d26d043fd95cf79cf5633185a..7d45d07e4fdfa31b3f408cc99328c71a1025d86b 100644
--- a/Examples/python/simulation/ex05_BeamAndDetector/OffSpecularSimulation.py
+++ b/Examples/python/simulation/ex05_BeamAndDetector/OffSpecularSimulation.py
@@ -6,7 +6,6 @@ from bornagain import deg, angstrom, nm
 
 phi_f_min, phi_f_max = -1.0, 1.0
 alpha_f_min, alpha_f_max = 0.0, 10.0
-
 alpha_i_min, alpha_i_max = 0.0, 10.0  # incoming beam
 
 
@@ -74,6 +73,4 @@ def run_simulation():
 
 if __name__ == '__main__':
     result = run_simulation()
-    ba.plot_colormap(result, xlabel=r'$\alpha_i ^{\circ}$')
-    from matplotlib import pyplot as plt
-    plt.show()
+    ba.plot_simulation_result(result, zmin=1.0)
diff --git a/Examples/python/simulation/ex05_BeamAndDetector/ResonatorOffSpecSetup.py b/Examples/python/simulation/ex05_BeamAndDetector/ResonatorOffSpecSetup.py
new file mode 100644
index 0000000000000000000000000000000000000000..cc6597095e60804fd4caca9d952b650ade30ab76
--- /dev/null
+++ b/Examples/python/simulation/ex05_BeamAndDetector/ResonatorOffSpecSetup.py
@@ -0,0 +1,95 @@
+"""
+Resonator in off-specular experiment.
+"""
+import bornagain as ba
+from bornagain import deg, nm, micrometer, angstrom
+
+
+def get_sample(nlayers = 3):
+    """
+    Construct resonator with given number of Ti/Pt double layers nlayers.
+    """
+
+    # define materials
+    m_Si = ba.HomogeneousMaterial("Si", 8.25218379931e-06, 0.0)
+    m_Ti = ba.HomogeneousMaterial("Ti", -7.6593316363e-06, 3.81961616312e-09)
+    m_TiO2 = ba.HomogeneousMaterial("TiO2", 1.04803530026e-05, 2.03233519385e-09)
+    m_Pt = ba.HomogeneousMaterial("Pt", 2.52936993309e-05, 7.54553992473e-09)
+    m_D2O = ba.HomogeneousMaterial("D2O", 2.52897204573e-05, 4.5224432814e-13)
+
+    # create layers
+    l_TiO2 = ba.Layer(m_TiO2, 3.0*nm)
+    l_Ti_top = ba.Layer(m_Ti, 10.0*nm)
+    l_Ti = ba.Layer(m_Ti, 13.0*nm)
+    l_Si = ba.Layer(m_Si)    # consider it as an ambient layer
+    l_Pt = ba.Layer(m_Pt, 32.0*nm)
+    l_D2O = ba.Layer(m_D2O)  # thickness is not given, seems to be like a substrate
+
+    # describe layer roughness
+    roughness = ba.LayerRoughness()
+    roughness.setSigma(2.0*nm)
+    roughness.setHurstParameter(0.8)
+    roughness.setLatteralCorrLength(10.0*micrometer)
+
+    # assemble multilayer
+    sample = ba.MultiLayer()
+    sample.addLayer(l_Si)  # Assume huge Si block to be infinite
+
+    for i in range(nlayers):
+        sample.addLayerWithTopRoughness(l_Ti, roughness)
+        sample.addLayerWithTopRoughness(l_Pt, roughness)
+
+    sample.addLayerWithTopRoughness(l_Ti_top, roughness)
+    sample.addLayerWithTopRoughness(l_TiO2, roughness)
+    sample.addLayerWithTopRoughness(l_D2O, roughness)
+
+    sample.setCrossCorrLength(400*nm)
+
+    return sample
+
+
+def get_offspec_simulation():
+    """
+    characterizing the input beam and output detector
+    """
+
+    # create OffSpecular simulation
+    simulation = ba.OffSpecSimulation()
+    simulation.setTerminalProgressMonitor()
+
+    # define detector parameters
+    n_alpha, alpha_min, alpha_max = 300, 0.0*deg, 4.0*deg
+    n_phi, phi_min, phi_max = 10, -0.1*deg, 0.1*deg
+    simulation.setDetectorParameters(
+        n_phi, phi_min, phi_max, n_alpha, alpha_min, alpha_max)
+
+    # define the beam with alpha_i varied between alpha_i_min and alpha_i_max
+    n_scan_points, alpha_i_min, alpha_i_max = n_alpha, alpha_min, alpha_max
+    alpha_i_axis = ba.FixedBinAxis(
+        "alpha_i", n_scan_points, alpha_i_min, alpha_i_max)
+    simulation.setBeamParameters(5.0*angstrom, alpha_i_axis, 0.0)
+
+    simulation.setBeamIntensity(1e9)
+    simulation.getOptions().setIncludeSpecular(True)
+
+    # define detector resolution function with smearing depending on bin size
+    d_alpha = (alpha_max - alpha_min)/n_alpha
+    d_phi = (phi_max-phi_min)/n_phi
+    sigma_factor = 1.0
+    simulation.setDetectorResolutionFunction(
+        ba.ResolutionFunction2DGaussian(sigma_factor*d_alpha, sigma_factor*d_phi))
+
+    return simulation
+
+
+def run_simulation():
+    sample = get_sample(nlayers=3)
+    simulation = get_offspec_simulation()
+    simulation.setSample(sample)
+    simulation.runSimulation()
+    return simulation.result()
+
+
+if __name__ == '__main__':
+    result = run_simulation()
+    ba.plot_simulation_result(result, zmin=1e-03)
diff --git a/Tests/Functional/Core/CoreStandardTest/CMakeLists.txt b/Tests/Functional/Core/CoreStandardTest/CMakeLists.txt
index 9f490363257c2c49b613c075ec3f34e4b786144b..181768c0dcebde60838e82dade7bedeccac32066 100644
--- a/Tests/Functional/Core/CoreStandardTest/CMakeLists.txt
+++ b/Tests/Functional/Core/CoreStandardTest/CMakeLists.txt
@@ -66,6 +66,7 @@ set(test_cases
     GaussianBeamFootprint
     SquareBeamFootprint
     SpecularDivergentBeam
+    OffSpecularResonator
     )
 
 add_executable(CoreStandardTest main.cpp CoreStandardTest.h CoreStandardTest.cpp)
diff --git a/Tests/Functional/Python/PyEmbedded/TestCases.cpp b/Tests/Functional/Python/PyEmbedded/TestCases.cpp
index c695d31a878031c66708f883142fbcdd4064ff22..b4764bbe8c291ff32adf4ce53cd8733319c880b9 100644
--- a/Tests/Functional/Python/PyEmbedded/TestCases.cpp
+++ b/Tests/Functional/Python/PyEmbedded/TestCases.cpp
@@ -23,6 +23,7 @@
 #include "BornAgainNamespace.h"
 #include "SampleBuilderFactory.h"
 #include "ExportToPython.h"
+#include "PythonFormatting.h"
 #include <iostream>
 #include <sstream>
 
@@ -363,12 +364,10 @@ bool ExportToPythonAndBack::runTest()
     auto code = ExportToPython::generateSampleCode(*sample);
 
     std::stringstream snippet;
-    snippet << "import bornagain as ba                                        \n";
-    snippet << "from bornagain import deg, angstrom, nm                     \n\n";
-    snippet << code;
+    snippet << PythonFormatting::scriptPreamble() << code;
 
-    auto multilayer = PyImport::createFromPython(snippet.str(), "getSample",
-                                                 BABuild::buildLibDir());
+    auto multilayer = PyImport::createFromPython(snippet.str(),
+            PythonFormatting::getSampleFunctionName(), BABuild::buildLibDir());
     auto new_code = ExportToPython::generateSampleCode(*multilayer);
 
     return code == new_code;
diff --git a/Tests/Functional/Python/PyStandard/CMakeLists.txt b/Tests/Functional/Python/PyStandard/CMakeLists.txt
index 36d8db852dfa16f4fd18e8800e65fda111958bcd..2c0b77e20872d1b8d7f2ae7c98cfcb574f83933f 100644
--- a/Tests/Functional/Python/PyStandard/CMakeLists.txt
+++ b/Tests/Functional/Python/PyStandard/CMakeLists.txt
@@ -56,6 +56,7 @@ set(test_cases
     RotatedPyramidsDistribution
     SpheresWithLimitsDistribution
     ConesWithLimitsDistribution
+    OffSpecularResonator
 )
 
 add_executable(PyStandardTest main.cpp PyStandardTest.h PyStandardTest.cpp)
diff --git a/Tests/Functional/Python/PyStandard/PyStandardTest.cpp b/Tests/Functional/Python/PyStandard/PyStandardTest.cpp
index 411a0ce1764da72a04258caa2a8c265e64cecad1..a79bfedd679c95aafc962c984f3450c167c04082 100644
--- a/Tests/Functional/Python/PyStandard/PyStandardTest.cpp
+++ b/Tests/Functional/Python/PyStandard/PyStandardTest.cpp
@@ -30,7 +30,7 @@ bool PyStandardTest::runTest()
         = FileSystemUtils::jointPath(BATesting::PyStandardOutputDir(), getName());
     std::string output_path = output_name + ".ref.int.gz";
     std::remove(output_path.c_str());
-    std::cout << "Removed old output " << output_path << "n";
+    std::cout << "Removed old output " << output_path << "\n";
 
     // Generate Python script
     std::string pyscript_filename
diff --git a/Tests/Functional/TestMachinery/StandardTestCatalogue.cpp b/Tests/Functional/TestMachinery/StandardTestCatalogue.cpp
index 1e94321dc98f896d6a397a6622cd95726fc0fcff..84dbbeaf195081dbe18aeb5ecda6ecc3d185abcc 100644
--- a/Tests/Functional/TestMachinery/StandardTestCatalogue.cpp
+++ b/Tests/Functional/TestMachinery/StandardTestCatalogue.cpp
@@ -408,6 +408,12 @@ StandardTestCatalogue::StandardTestCatalogue()
         "SpecularDivergentBeam",
         "HomogeneousMultilayerBuilder",
         1e-10);
+
+    add("OffSpecularResonator",
+        "Simulates resonator in OffSpec setup",
+        "OffSpecMini",
+        "ResonatorBuilder",
+        1e-10);
 }
 
 //! Adds test description to the catalogue.
diff --git a/Tests/ReferenceData/Core/OffSpecularResonator.int.gz b/Tests/ReferenceData/Core/OffSpecularResonator.int.gz
new file mode 100644
index 0000000000000000000000000000000000000000..5c34fd0725fb1d6036e0ab72d873fda7142782d9
Binary files /dev/null and b/Tests/ReferenceData/Core/OffSpecularResonator.int.gz differ
diff --git a/Tests/UnitTests/Core/ExportToPython/PythonFormattingTest.h b/Tests/UnitTests/Core/ExportToPython/PythonFormattingTest.h
index 6571aef8fc6259ce2af5ba9b7957bba94efee882..b44c527d57b697da052590823b5dd00ac4c0db50 100644
--- a/Tests/UnitTests/Core/ExportToPython/PythonFormattingTest.h
+++ b/Tests/UnitTests/Core/ExportToPython/PythonFormattingTest.h
@@ -5,6 +5,7 @@
 #include "PythonFormatting.h"
 #include "RealLimits.h"
 #include "Units.h"
+#include "FixedBinAxis.h"
 
 class PythonFormattingTest : public ::testing::Test
 {
@@ -106,3 +107,13 @@ TEST_F(PythonFormattingTest, printParameterDistribution)
               "ba.ParameterDistribution(\"/Particle/ZRotation/Angle\", "
               "distr_1, 5, 2.0, ba.RealLimits.limited(1.0*deg, 2.0*deg))");
 }
+
+TEST_F(PythonFormattingTest, printAxis)
+{
+    FixedBinAxis axis1("axis0", 10, -1.0, 2.0);
+    EXPECT_EQ(PythonFormatting::printAxis(axis1), "ba.FixedBinAxis(\"axis0\", 10, -1.0, 2.0)");
+
+    FixedBinAxis axis2("axis0", 10, -1.0*Units::deg, 2.0*Units::deg);
+    EXPECT_EQ(PythonFormatting::printAxis(axis2, BornAgain::UnitsRad),
+              "ba.FixedBinAxis(\"axis0\", 10, -1.0*deg, 2.0*deg)");
+}