diff --git a/Core/Instrument/IDetector2D.cpp b/Core/Instrument/IDetector2D.cpp
index 716efb17c1b2d37a21cf897b41d82641c1ba5105..960c9b41a6f73f28f5bb8911698cbd7acdbc0f3d 100644
--- a/Core/Instrument/IDetector2D.cpp
+++ b/Core/Instrument/IDetector2D.cpp
@@ -20,6 +20,9 @@
 #include "Logger.h"
 #include "SimulationElement.h"
 #include "SimulationArea.h"
+#include "BornAgainNamespace.h"
+#include "Units.h"
+#include "Exceptions.h"
 
 IDetector2D::IDetector2D()
     : m_axes()
@@ -110,13 +113,14 @@ std::string IDetector2D::addParametersToExternalPool(
     return new_path;
 }
 
-OutputData<double>* IDetector2D::createDetectorMap(const Beam&, EAxesUnits) const
+OutputData<double>* IDetector2D::createDetectorMap(const Beam& beam, EAxesUnits units) const
 {
-    return nullptr;
+    std::unique_ptr<OutputData<double>> result(new OutputData<double>);
+    result->addAxis(*constructAxis(BornAgain::X_AXIS_INDEX, beam, units));
+    result->addAxis(*constructAxis(BornAgain::Y_AXIS_INDEX, beam, units));
+    return result.release();
 }
 
-//!
-
 void IDetector2D::initOutputData(OutputData<double> &data) const {
   data.clear();
   for (size_t i = 0; i < getDimension(); ++i)
@@ -266,6 +270,53 @@ size_t IDetector2D::getAxisBinIndex(size_t index, size_t selected_axis) const
                                           "Error! No axis with given number");
 }
 
+std::unique_ptr<IAxis> IDetector2D::constructAxis(size_t axis_index, const Beam &beam,
+                                                  IDetector2D::EAxesUnits units) const
+{
+    double amin(0), amax(0);
+    calculateAxisRange(axis_index, beam, units, amin, amax);
+    return std::unique_ptr<IAxis>(
+        new FixedBinAxis(getAxisName(axis_index), getAxis(axis_index).getSize(), amin, amax));
+}
+
+void IDetector2D::calculateAxisRange(size_t axis_index, const Beam &beam,
+        IDetector2D::EAxesUnits units, double &amin, double &amax) const
+{
+    amin = 0.0; amax=0.0;
+
+    if(units == DEFAULT) {
+        amin = getAxis(axis_index).getMin();
+        amax = getAxis(axis_index).getMax();
+
+    }else if(units == NBINS) {
+        amin = 0.0;
+        amax = static_cast<double>(getAxis(axis_index).getSize());
+
+    } else if(units == QYQZ && axis_index == BornAgain::X_AXIS_INDEX) {
+        const IAxis &aX = getAxis(BornAgain::X_AXIS_INDEX);
+        SimulationElement el_left_bottom
+            = getSimulationElement(getGlobalIndex(0, 0), beam);
+        SimulationElement el_right_bottom
+            = getSimulationElement(getGlobalIndex(aX.getSize()-1, 0), beam);
+        amin = el_left_bottom.getQ(0.0, 0.0).y();
+        amax = el_right_bottom.getQ(1.0, 0.0).y();
+
+    } else if(units == QYQZ && axis_index == BornAgain::Y_AXIS_INDEX) {
+        const IAxis &aX = getAxis(BornAgain::X_AXIS_INDEX);
+        const IAxis &aY = getAxis(BornAgain::Y_AXIS_INDEX);
+        SimulationElement el_center_bottom
+            = getSimulationElement(getGlobalIndex(aX.getSize()/2, 0), beam);
+        SimulationElement el_center_top
+            = getSimulationElement(getGlobalIndex(aX.getSize()/2, aY.getSize()-1), beam);
+        amin = -el_center_bottom.getQ(0.5, 0.0).z();
+        amax = -el_center_top.getQ(0.5, 1.0).z();
+
+    } else {
+        throw Exceptions::RuntimeErrorException("IDetector2D::calculateAxisRange() -> Error. "
+                                                "Unknown units " +std::to_string(units));
+    }
+}
+
 size_t IDetector2D::getGlobalIndex(size_t x, size_t y) const
 {
     if (getDimension()!=2) return getTotalSize();
diff --git a/Core/Instrument/IDetector2D.h b/Core/Instrument/IDetector2D.h
index 3d723889cd3e730f914517e4096ffb172b1f13ab..a09cf26f177149dd238dad12dbc54644cf7e997a 100644
--- a/Core/Instrument/IDetector2D.h
+++ b/Core/Instrument/IDetector2D.h
@@ -128,7 +128,7 @@ public:
         const std::string& path, ParameterPool* external_pool, int copy_number = -1) const;
 
     //! Returns detector map in given axes units
-    virtual OutputData<double>* createDetectorMap(const Beam&, EAxesUnits) const;
+    virtual OutputData<double>* createDetectorMap(const Beam& beam, EAxesUnits units) const;
 
     //! Inits axes of OutputData to match the detector and sets values to zero.
     virtual void initOutputData(OutputData<double> &data) const;
@@ -164,6 +164,14 @@ protected:
     //! Generates an axis with correct name and default binning for given index
     virtual IAxis* createAxis(size_t index, size_t n_bins, double min, double max) const=0;
 
+    //! Constructs axis with min,max corresponding to selected units
+    std::unique_ptr<IAxis> constructAxis(size_t axis_index, const Beam& beam,
+                                         EAxesUnits units) const;
+
+    //! Calculates axis range from original detector axes in given units (mm, rad, etc)
+    virtual void calculateAxisRange(size_t axis_index, const Beam& beam, EAxesUnits units,
+                                    double &amin, double &amax) const;
+
     //! Returns the name for the axis with given index
     virtual std::string getAxisName(size_t index) const=0;
 
diff --git a/Core/Instrument/RectangularDetector.cpp b/Core/Instrument/RectangularDetector.cpp
index e263734473fbff9a5c230236e51d4e9f5267a43a..fe044b83db2a8a6e45a7685b4124527b2a54d66c 100644
--- a/Core/Instrument/RectangularDetector.cpp
+++ b/Core/Instrument/RectangularDetector.cpp
@@ -202,69 +202,6 @@ RectangularDetector::EDetectorArrangement RectangularDetector::getDetectorArrang
     return m_detector_arrangement;
 }
 
-OutputData<double> *RectangularDetector::createDetectorMap(
-    const Beam& beam, IDetector2D::EAxesUnits units_type) const
-{
-    if (getDimension() != 2)
-        return 0;
-
-    std::unique_ptr<OutputData<double>> result(new OutputData<double>);
-    const IAxis& aX = getAxis(BornAgain::X_AXIS_INDEX);
-    const IAxis& aY = getAxis(BornAgain::Y_AXIS_INDEX);
-
-    result->addAxis(aX);
-    result->addAxis(aY);
-
-    if (units_type == DEFAULT)
-        return result.release();
-
-    std::vector<int> indices_left_bottom = {0, 0};
-    std::vector<int> indices_right_bottom = {(int)aX.getSize() - 1, 0};
-    std::vector<int> indices_center_bottom = {(int)aX.getSize() / 2, 0};
-    std::vector<int> indices_center_top = {(int)aX.getSize() / 2, (int)aY.getSize() - 1};
-    SimulationElement el_left_bottom
-        = getSimulationElement(result->toGlobalIndex(indices_left_bottom), beam);
-    SimulationElement el_right_bottom
-        = getSimulationElement(result->toGlobalIndex(indices_right_bottom), beam);
-    SimulationElement el_center_bottom
-        = getSimulationElement(result->toGlobalIndex(indices_center_bottom), beam);
-    SimulationElement el_center_top
-        = getSimulationElement(result->toGlobalIndex(indices_center_top), beam);
-
-    result->clear();
-
-    double xmin(aX.getMin()), xmax(aX.getMax()), ymin(aY.getMin()), ymax(aY.getMax());
-
-    if (units_type == NBINS) {
-        xmin = 0.0;
-        ymin = 0.0;
-        xmax = double(aX.getSize());
-        ymax = double(aY.getSize());
-    }
-
-    else if (units_type == RADIANS || units_type == DEGREES) {
-        double scale(1.0);
-        if (units_type == DEGREES)
-            scale = 1. / Units::degree;
-        xmin = scale * el_left_bottom.getPhi(0.0, 0.0);
-        xmax = scale * el_right_bottom.getPhi(1.0, 0.0);
-        ymin = scale * el_center_bottom.getAlpha(0.5, 0.0);
-        ymax = scale * el_center_top.getAlpha(0.5, 1.0);
-    }
-
-    else if (units_type == QYQZ) {
-        xmin = el_left_bottom.getQ(0.0, 0.0).y();
-        xmax = el_right_bottom.getQ(1.0, 0.0).y();
-        ymin = -el_center_bottom.getQ(0.5, 0.0).z();
-        ymax = -el_center_top.getQ(0.5, 1.0).z();
-    }
-
-    result->addAxis(FixedBinAxis(BornAgain::U_AXIS_NAME, aX.getSize(), xmin, xmax));
-    result->addAxis(FixedBinAxis(BornAgain::V_AXIS_NAME, aY.getSize(), ymin, ymax));
-
-    return result.release();
-}
-
 std::vector<IDetector2D::EAxesUnits> RectangularDetector::getValidAxesUnits() const
 {
     std::vector<IDetector2D::EAxesUnits> result = IDetector2D::getValidAxesUnits();
@@ -297,6 +234,43 @@ IAxis *RectangularDetector::createAxis(size_t index, size_t n_bins, double min,
     return new FixedBinAxis(getAxisName(index), n_bins, min, max);
 }
 
+void RectangularDetector::calculateAxisRange(size_t axis_index, const Beam &beam,
+    IDetector2D::EAxesUnits units, double &amin, double &amax) const
+{
+    amin = 0.0; amax=0.0;
+    if(units == MM) {
+        amin = getAxis(axis_index).getMin();
+        amax = getAxis(axis_index).getMax();
+    }else if(units == RADIANS || units == DEGREES) {
+        double scale(1.0);
+        if (units == DEGREES)
+            scale = 1. / Units::degree;
+
+        if(axis_index == BornAgain::X_AXIS_INDEX) {
+            const IAxis &aX = getAxis(BornAgain::X_AXIS_INDEX);
+            SimulationElement el_left_bottom
+                = getSimulationElement(getGlobalIndex(0, 0), beam);
+            SimulationElement el_right_bottom
+                = getSimulationElement(getGlobalIndex(aX.getSize()-1, 0), beam);
+            amin = scale * el_left_bottom.getPhi(0.0, 0.0);
+            amax = scale * el_right_bottom.getPhi(1.0, 0.0);
+        } else if(axis_index == BornAgain::Y_AXIS_INDEX) {
+            const IAxis &aX = getAxis(BornAgain::X_AXIS_INDEX);
+            const IAxis &aY = getAxis(BornAgain::Y_AXIS_INDEX);
+            SimulationElement el_center_bottom
+                = getSimulationElement(getGlobalIndex(aX.getSize()/2, 0), beam);
+            SimulationElement el_center_top
+                = getSimulationElement(getGlobalIndex(aX.getSize()/2, aY.getSize()-1), beam);
+            amin = scale * el_center_bottom.getAlpha(0.5, 0.0);
+            amax = scale * el_center_top.getAlpha(0.5, 1.0);
+        }
+
+    } else {
+        IDetector2D::calculateAxisRange(axis_index, beam, units, amin, amax);
+    }
+
+}
+
 std::string RectangularDetector::getAxisName(size_t index) const
 {
     switch (index) {
diff --git a/Core/Instrument/RectangularDetector.h b/Core/Instrument/RectangularDetector.h
index eb6f5f6834a3665be869c8f8a4dc91f0bdd1697c..dd5f5ea9042585aa149a8f18eb5a334b23659da2 100644
--- a/Core/Instrument/RectangularDetector.h
+++ b/Core/Instrument/RectangularDetector.h
@@ -76,9 +76,6 @@ public:
     double getDirectBeamV0() const;
     EDetectorArrangement getDetectorArrangment() const;
 
-    //! Returns detector map in given axes units
-    OutputData<double>* createDetectorMap(const Beam& beam, EAxesUnits units_type) const override;
-
     //! returns vector of valid axes units
     std::vector<EAxesUnits> getValidAxesUnits() const override;
 
@@ -97,6 +94,10 @@ protected:
     //! Generates an axis with correct name and default binning for given index
     IAxis* createAxis(size_t index, size_t n_bins, double min, double max) const override;
 
+    //! Calculates axis range from original detector axes in given units (mm, rad, etc)
+    virtual void calculateAxisRange(size_t axis_index, const Beam& beam, EAxesUnits units,
+                                    double &amin, double &amax) const override;
+
     //! Returns the name for the axis with given index
     std::string getAxisName(size_t index) const override;
 
diff --git a/Core/Instrument/SphericalDetector.cpp b/Core/Instrument/SphericalDetector.cpp
index 2fc96693c114f48bfe1998f15d7343b5d9362de0..e7fe65bf03f6e4bfdbba1f9a2664a5539137f0f9 100644
--- a/Core/Instrument/SphericalDetector.cpp
+++ b/Core/Instrument/SphericalDetector.cpp
@@ -69,70 +69,6 @@ std::string SphericalDetector::addParametersToExternalPool(
     return new_path;
 }
 
-OutputData<double>* SphericalDetector::createDetectorMap(const Beam& beam,
-                                                         IDetector2D::EAxesUnits units_type) const
-{
-    if (getDimension() != 2)
-        return 0;
-
-    std::unique_ptr<OutputData<double>> result(new OutputData<double>);
-
-    const IAxis& aX = getAxis(BornAgain::X_AXIS_INDEX);
-    const IAxis& aY = getAxis(BornAgain::Y_AXIS_INDEX);
-
-    result->addAxis(aX);
-    result->addAxis(aY);
-
-    if (units_type == DEFAULT)
-        return result.release();
-
-    std::vector<int> indices_left_bottom = {0, 0};
-    std::vector<int> indices_right_bottom = {(int)aX.getSize() - 1, 0};
-    std::vector<int> indices_center_bottom = {(int)aX.getSize() / 2, 0};
-    std::vector<int> indices_center_top = {(int)aX.getSize() / 2, (int)aY.getSize() - 1};
-    SimulationElement el_left_bottom
-        = getSimulationElement(result->toGlobalIndex(indices_left_bottom), beam);
-    SimulationElement el_right_bottom
-        = getSimulationElement(result->toGlobalIndex(indices_right_bottom), beam);
-    SimulationElement el_center_bottom
-        = getSimulationElement(result->toGlobalIndex(indices_center_bottom), beam);
-    SimulationElement el_center_top
-        = getSimulationElement(result->toGlobalIndex(indices_center_top), beam);
-
-    result->clear();
-
-    double xmin(aX.getMin()), xmax(aX.getMax()), ymin(aY.getMin()), ymax(aY.getMax());
-
-    if (units_type == MM)
-        return 0;
-
-    else if (units_type == NBINS) {
-        xmin = 0.0;
-        ymin = 0.0;
-        xmax = double(aX.getSize());
-        ymax = double(aY.getSize());
-    }
-
-    else if (units_type == DEGREES) {
-        xmin = aX.getMin() / Units::degree;
-        xmax = aX.getMax() / Units::degree;
-        ymin = aY.getMin() / Units::degree;
-        ymax = aY.getMax() / Units::degree;
-    }
-
-    else if (units_type == QYQZ) {
-        xmin = el_left_bottom.getQ(0.0, 0.0).y();
-        xmax = el_right_bottom.getQ(1.0, 0.0).y();
-        ymin = -el_center_bottom.getQ(0.5, 0.0).z();
-        ymax = -el_center_top.getQ(0.5, 1.0).z();
-    }
-
-    result->addAxis(FixedBinAxis(BornAgain::PHI_AXIS_NAME, aX.getSize(), xmin, xmax));
-    result->addAxis(FixedBinAxis(BornAgain::ALPHA_AXIS_NAME, aY.getSize(), ymin, ymax));
-
-    return result.release();
-}
-
 std::vector<IDetector2D::EAxesUnits> SphericalDetector::getValidAxesUnits() const
 {
     std::vector<IDetector2D::EAxesUnits> result = IDetector2D::getValidAxesUnits();
@@ -177,6 +113,21 @@ IAxis* SphericalDetector::createAxis(size_t index, size_t n_bins, double min, do
     return new FixedBinAxis(getAxisName(index), n_bins, min, max);
 }
 
+void SphericalDetector::calculateAxisRange(size_t axis_index, const Beam &beam,
+        IDetector2D::EAxesUnits units, double &amin, double &amax) const
+{
+    amin = 0.0; amax=0.0;
+    if(units == DEGREES) {
+        amin = getAxis(axis_index).getMin()/Units::degree;
+        amax = getAxis(axis_index).getMax()/Units::degree;
+    }else if(units == RADIANS) {
+        amin = getAxis(axis_index).getMin();
+        amax = getAxis(axis_index).getMax();
+    } else {
+        IDetector2D::calculateAxisRange(axis_index, beam, units, amin, amax);
+    }
+}
+
 std::string SphericalDetector::getAxisName(size_t index) const
 {
     switch (index) {
diff --git a/Core/Instrument/SphericalDetector.h b/Core/Instrument/SphericalDetector.h
index fe38df6a952315e5f5a98d1d41c0273b4b758df0..302e2f7427da904460f9d3116ac341e99b422446 100644
--- a/Core/Instrument/SphericalDetector.h
+++ b/Core/Instrument/SphericalDetector.h
@@ -49,9 +49,6 @@ public:
     std::string addParametersToExternalPool(
         const std::string& path, ParameterPool* external_pool, int copy_number = -1) const override;
 
-    //! Returns detector map in given axes units
-    OutputData<double> *createDetectorMap(const Beam& beam, EAxesUnits units_type) const override;
-
     //! returns vector of valid axes units
     std::vector<EAxesUnits> getValidAxesUnits() const override;
 
@@ -70,6 +67,10 @@ protected:
     //! Generates an axis with correct name and default binning for given index
     IAxis* createAxis(size_t index, size_t n_bins, double min, double max) const override;
 
+    //! Calculates axis range from original detector axes in given units (mm, rad, etc)
+    virtual void calculateAxisRange(size_t axis_index, const Beam& beam, EAxesUnits units,
+                                    double &amin, double &amax) const override;
+
     //! Returns the name for the axis with given index
     std::string getAxisName(size_t index) const override;
 
diff --git a/Core/Simulation/GISASSimulation.cpp b/Core/Simulation/GISASSimulation.cpp
index 32a384323f8cab4e180a64d53a8dc8ca12a8c4fc..10bbd93639325f5c0f63d63ddbcdfdc976598173 100644
--- a/Core/Simulation/GISASSimulation.cpp
+++ b/Core/Simulation/GISASSimulation.cpp
@@ -58,7 +58,7 @@ int GISASSimulation::getNumberOfSimulationElements() const
         throw Exceptions::RuntimeErrorException("GISASSimulation::getNumberOfSimulationElements: "
                                     "detector is not two-dimensional");
     const IAxis& x_axis = m_instrument.getDetectorAxis(BornAgain::X_AXIS_INDEX);
-    const IAxis& y_axis = m_instrument.getDetectorAxis(BornAgain::X_AXIS_INDEX);
+    const IAxis& y_axis = m_instrument.getDetectorAxis(BornAgain::Y_AXIS_INDEX);
     int nmasked = getInstrument().getDetector()->getNumberOfMaskedChannels();
     return x_axis.getSize()*y_axis.getSize() - nmasked;
 }