diff --git a/Base/Axis/Frame.cpp b/Base/Axis/Frame.cpp
index 2142b974babec63e4f0f541f3098345b082fdaab..4f3ae9dd64c4e92418a8e4e8f488ed5413870b91 100644
--- a/Base/Axis/Frame.cpp
+++ b/Base/Axis/Frame.cpp
@@ -139,11 +139,12 @@ bool Frame::hasSameSizes(const Frame& o) const
     return true;
 }
 
-Frame* Frame::plottableFrame() const
+Frame* Frame::convertedFrame(const std::vector<Unit>& newUnits) const
 {
+    ASSERT(newUnits.size() == m_axes.size());
     std::vector<const Scale*> outaxes;
-    for (const Scale* s : m_axes)
-        outaxes.emplace_back(new Scale(s->plottableScale()));
+    for (size_t i = 0; i < m_axes.size(); i++)
+        outaxes.emplace_back(new Scale(m_axes[i]->convertedScale(newUnits[i])));
     return new Frame(std::move(outaxes));
 }
 
diff --git a/Base/Axis/Frame.h b/Base/Axis/Frame.h
index 1cfba5a41e5dad71410c44998c69ccfbb6cdb9f0..d3ce46aac32d38d625121013a508551838914d61 100644
--- a/Base/Axis/Frame.h
+++ b/Base/Axis/Frame.h
@@ -15,6 +15,7 @@
 #ifndef BORNAGAIN_BASE_AXIS_FRAME_H
 #define BORNAGAIN_BASE_AXIS_FRAME_H
 
+#include "Base/Const/Units.h"
 #include "Base/Types/ICloneable.h"
 #include "Base/Types/OwningVector.h"
 
@@ -81,7 +82,7 @@ public:
     std::vector<const Scale*> clonedAxes() const;
 #endif // SWIG
 
-    Frame* plottableFrame() const;
+    Frame* convertedFrame(const std::vector<Unit>& newUnits) const;
     Frame* flat() const;
 
 protected:
diff --git a/Base/Axis/MakeScale.cpp b/Base/Axis/MakeScale.cpp
index 87981964760b01bbf9dff0a59460eefd034eef42..2fdd3aee46671fd2b83157ce98effec725980794 100644
--- a/Base/Axis/MakeScale.cpp
+++ b/Base/Axis/MakeScale.cpp
@@ -16,31 +16,31 @@
 #include "Base/Axis/Scale.h"
 #include "Base/Util/Assert.h"
 
-Scale* newGenericScale(const std::string& name, const std::vector<double>& limits)
+Scale* newGenericScale(const std::string& name, Unit units, const std::vector<double>& limits)
 {
     std::vector<Bin1D> vec;
     ASSERT(!(limits.size() & 1)); // must be an even number
     for (size_t i = 0; i < limits.size(); i += 2)
         vec.emplace_back(Bin1D::FromTo(limits[i], limits[i + 1]));
-    return new Scale(name, vec);
+    return new Scale(name, units, vec);
 }
 
 
-Scale ListScan(const std::string& name, const std::vector<double>& points)
+Scale ListScan(const std::string& name, Unit units, const std::vector<double>& points)
 {
     std::vector<Bin1D> vec;
     for (double p : points)
         vec.emplace_back(Bin1D::At(p));
-    return Scale(name, vec);
+    return Scale(name, units, vec);
 }
 
-Scale* newListScan(const std::string& name, const std::vector<double>& points)
+Scale* newListScan(const std::string& name, Unit units, const std::vector<double>& points)
 {
-    return new Scale(ListScan(name, points));
+    return new Scale(ListScan(name, units, points));
 }
 
 
-Scale EquiDivision(const std::string& name, size_t N, double start, double end)
+Scale EquiDivision(const std::string& name, Unit units, size_t N, double start, double end)
 {
     ASSERT(N > 0);
     ASSERT(start <= end);
@@ -50,28 +50,28 @@ Scale EquiDivision(const std::string& name, size_t N, double start, double end)
                                        (N - i - 1) * (start / N) + (i + 1) * (end / N)));
     ASSERT(vec.size() == N);
 
-    return Scale(name, vec);
+    return Scale(name, units, vec);
 }
 
-Scale* newEquiDivision(const std::string& name, size_t N, double start, double end)
+Scale* newEquiDivision(const std::string& name, Unit units, size_t N, double start, double end)
 {
-    return new Scale(EquiDivision(name, N, start, end));
+    return new Scale(EquiDivision(name, units, N, start, end));
 }
 
-std::shared_ptr<Scale> sharedEquiDivision(const std::string& name, size_t N, double start,
-                                          double end)
+std::shared_ptr<Scale> sharedEquiDivision(const std::string& name, Unit units, size_t N,
+                                          double start, double end)
 {
-    return std::shared_ptr<Scale>(newEquiDivision(name, N, start, end));
+    return std::shared_ptr<Scale>(newEquiDivision(name, units, N, start, end));
 }
 
-std::unique_ptr<Scale> uniqueEquiDivision(const std::string& name, size_t N, double start,
-                                          double end)
+std::unique_ptr<Scale> uniqueEquiDivision(const std::string& name, Unit units, size_t N,
+                                          double start, double end)
 {
-    return std::unique_ptr<Scale>(newEquiDivision(name, N, start, end));
+    return std::unique_ptr<Scale>(newEquiDivision(name, units, N, start, end));
 }
 
 
-Scale EquiScan(const std::string& name, size_t N, double start, double end)
+Scale EquiScan(const std::string& name, Unit units, size_t N, double start, double end)
 {
     ASSERT(N);
     std::vector<Bin1D> vec;
@@ -85,20 +85,22 @@ Scale EquiScan(const std::string& name, size_t N, double start, double end)
         ASSERT(vec.size() == N);
     }
 
-    return Scale(name, vec);
+    return Scale(name, units, vec);
 }
 
-Scale* newEquiScan(const std::string& name, size_t N, double start, double end)
+Scale* newEquiScan(const std::string& name, Unit units, size_t N, double start, double end)
 {
-    return new Scale(EquiScan(name, N, start, end));
+    return new Scale(EquiScan(name, units, N, start, end));
 }
 
-std::shared_ptr<Scale> sharedEquiScan(const std::string& name, size_t N, double start, double end)
+std::shared_ptr<Scale> sharedEquiScan(const std::string& name, Unit units, size_t N, double start,
+                                      double end)
 {
-    return std::shared_ptr<Scale>(newEquiDivision(name, N, start, end));
+    return std::shared_ptr<Scale>(newEquiDivision(name, units, N, start, end));
 }
 
-std::unique_ptr<Scale> uniqueEquiScan(const std::string& name, size_t N, double start, double end)
+std::unique_ptr<Scale> uniqueEquiScan(const std::string& name, Unit units, size_t N, double start,
+                                      double end)
 {
-    return std::unique_ptr<Scale>(newEquiDivision(name, N, start, end));
+    return std::unique_ptr<Scale>(newEquiDivision(name, units, N, start, end));
 }
diff --git a/Base/Axis/MakeScale.h b/Base/Axis/MakeScale.h
index 104a1b5f69c40dbc2a4f9b5f711d0a7ee002ef93..26b98212644f4ecba024eefdec95e602e0bf6b7a 100644
--- a/Base/Axis/MakeScale.h
+++ b/Base/Axis/MakeScale.h
@@ -15,6 +15,7 @@
 #ifndef BORNAGAIN_BASE_AXIS_MAKESCALE_H
 #define BORNAGAIN_BASE_AXIS_MAKESCALE_H
 
+#include "Base/Const/Units.h"
 #include <memory>
 #include <string>
 #include <vector>
@@ -22,30 +23,32 @@
 class Scale;
 
 #ifndef SWIG
-Scale* newGenericScale(const std::string& name, const std::vector<double>& limits);
+Scale* newGenericScale(const std::string& name, Unit units, const std::vector<double>& limits);
 #endif // SWIG
 
-Scale ListScan(const std::string& name, const std::vector<double>& points);
+Scale ListScan(const std::string& name, Unit units, const std::vector<double>& points);
 #ifndef SWIG
-Scale* newListScan(const std::string& name, const std::vector<double>& points);
+Scale* newListScan(const std::string& name, Unit units, const std::vector<double>& points);
 #endif // SWIG
 
 //! Returns axis with fixed bin size.
-Scale EquiDivision(const std::string& name, size_t N, double start, double end);
+Scale EquiDivision(const std::string& name, Unit units, size_t N, double start, double end);
 #ifndef SWIG
-Scale* newEquiDivision(const std::string& name, size_t N, double start, double end);
-std::shared_ptr<Scale> sharedEquiDivision(const std::string& name, size_t N, double start,
-                                          double end);
-std::unique_ptr<Scale> uniqueEquiDivision(const std::string& name, size_t N, double start,
-                                          double end);
+Scale* newEquiDivision(const std::string& name, Unit units, size_t N, double start, double end);
+std::shared_ptr<Scale> sharedEquiDivision(const std::string& name, Unit units, size_t N,
+                                          double start, double end);
+std::unique_ptr<Scale> uniqueEquiDivision(const std::string& name, Unit units, size_t N,
+                                          double start, double end);
 #endif // SWIG
 
 //! Returns axis with equidistant points (zero bin width).
-Scale EquiScan(const std::string& name, size_t N, double start, double end);
+Scale EquiScan(const std::string& name, Unit units, size_t N, double start, double end);
 #ifndef SWIG
-Scale* newEquiScan(const std::string& name, size_t N, double start, double end);
-std::shared_ptr<Scale> sharedEquiScan(const std::string& name, size_t N, double start, double end);
-std::unique_ptr<Scale> uniqueEquiScan(const std::string& name, size_t N, double start, double end);
+Scale* newEquiScan(const std::string& name, Unit units, size_t N, double start, double end);
+std::shared_ptr<Scale> sharedEquiScan(const std::string& name, Unit units, size_t N, double start,
+                                      double end);
+std::unique_ptr<Scale> uniqueEquiScan(const std::string& name, Unit units, size_t N, double start,
+                                      double end);
 #endif // SWIG
 
 #endif // BORNAGAIN_BASE_AXIS_MAKESCALE_H
diff --git a/Base/Axis/Scale.cpp b/Base/Axis/Scale.cpp
index 3c6a8beef942c90ce8a69c3fc6af3d138f821e37..64bed9224470e3889dc08cef3adca5d6d258c66c 100644
--- a/Base/Axis/Scale.cpp
+++ b/Base/Axis/Scale.cpp
@@ -13,14 +13,32 @@
 //  ************************************************************************************************
 
 #include "Base/Axis/Scale.h"
+#include "Base/Math/Numeric.h"
 #include "Base/Util/Assert.h"
+#include <functional>
 #include <iomanip>
 #include <iostream>
 #include <numbers>
 using std::numbers::pi;
 
-Scale::Scale(std::string name, const std::vector<Bin1D>& bins)
+namespace {
+
+std::vector<Bin1D> convertBins(const std::vector<Bin1D>& bins, std::function<double(double)> f)
+{
+    std::vector<Bin1D> result;
+    for (const Bin1D b : bins) {
+        double bmi = f(b.lowerBound());
+        double bma = f(b.upperBound());
+        result.emplace_back(Bin1D::FromTo(bmi, bma));
+    }
+    return result;
+}
+
+} // namespace
+
+Scale::Scale(std::string name, Unit units, const std::vector<Bin1D>& bins)
     : m_name(name)
+    , m_units(units)
     , m_bins(bins)
 {
     ASSERT(size() > 0);
@@ -33,6 +51,11 @@ Scale::Scale(std::string name, const std::vector<Bin1D>& bins)
             ASSERT(!b.binSize()); // bin widths must be all zero or all positive
 }
 
+std::string Scale::axisLabel() const
+{
+    return m_name + " (" + Units::plainName(m_units) + ")";
+}
+
 size_t Scale::size() const
 {
     return m_bins.size();
@@ -100,8 +123,10 @@ bool Scale::isEquiDivision() const
     for (size_t i = 0; i < N; ++i) {
         const Bin1D& b = bin(i);
         // exactly replicate the computation of bin bounds in the EquiDivision factory function
-        if (b.lowerBound() != (N - i) * (min() / N) + i * (max() / N)
-            || b.upperBound() != (N - i - 1) * (min() / N) + (i + 1) * (max() / N))
+        double refLowerBound = (N - i) * (min() / N) + i * (max() / N);
+        double refUpperBound = (N - i - 1) * (min() / N) + (i + 1) * (max() / N);
+        if (!Numeric::almostEqual(b.upperBound(), refUpperBound, 1)
+            || !Numeric::almostEqual(b.lowerBound(), refLowerBound, 1))
             return false;
     }
     return true;
@@ -109,16 +134,16 @@ bool Scale::isEquiDivision() const
 
 bool Scale::isEquiScan() const
 {
+    if (!isScan())
+        return false;
+
     const size_t N = size();
-    ASSERT(N);
     if (N == 1)
-        return !bin(0).binSize();
+        return true;
     for (size_t i = 0; i < N; ++i) {
-        const Bin1D& b = bin(i);
-        if (b.binSize())
-            return false;
-        // exactly replicate the computation of bin bounds in the EquiDivision factory function
-        if (b.lowerBound() != (N - 1 - i) * (min() / (N - 1)) + i * (max() / (N - 1)))
+        // exactly replicate the computation of bin bounds in the EquiScan factory function
+        double refCenter = (N - 1 - i) * (min() / (N - 1)) + i * (max() / (N - 1));
+        if (!Numeric::almostEqual(bin(i).lowerBound(), refCenter, 1))
             return false;
     }
     return true;
@@ -138,7 +163,7 @@ Scale Scale::clipped(double lower, double upper) const
     for (const Bin1D& b : m_bins)
         if (auto bc = b.clipped_or_nil(lower, upper))
             out_bins.emplace_back(bc.value());
-    return {m_name, out_bins};
+    return {m_name, m_units, out_bins};
 }
 
 Scale Scale::clipped(std::pair<double, double> bounds) const
@@ -148,18 +173,19 @@ Scale Scale::clipped(std::pair<double, double> bounds) const
 
 bool Scale::operator==(const Scale& other) const
 {
-    return m_name == other.m_name && m_bins == other.m_bins;
+    return m_name == other.m_name && m_units == other.m_units && m_bins == other.m_bins;
 }
 
 std::ostream& operator<<(std::ostream& ostr, const Scale& ax)
 {
     size_t N = ax.size();
+    std::string name_units = ax.axisName() + " (" + Units::plainName(ax.axisUnits()) + ")";
     ASSERT(N > 0);
 
     ostr << std::setprecision(15);
 
     if (ax.isScan()) {
-        ostr << "ListScan(\"" << ax.axisName() << "\", [";
+        ostr << "ListScan(\"" << name_units << "\", [";
         for (double v : ax.binCenters())
             ostr << v << ",";
         ostr << "])";
@@ -167,43 +193,38 @@ std::ostream& operator<<(std::ostream& ostr, const Scale& ax)
     }
 
     if (ax.isEquiDivision()) {
-        ostr << "EquiDivision(\"" << ax.axisName() << "\", " << ax.size() << ", " << ax.min()
-             << ", " << ax.max() << ")";
+        ostr << "EquiDivision(\"" << name_units << "\", " << ax.size() << ", " << ax.min() << ", "
+             << ax.max() << ")";
         return ostr;
     }
 
-    ostr << "Scale(\"" << ax.axisName() << "\", [";
+    ostr << "Scale(\"" << name_units << "\", [";
     for (const Bin1D& b : ax.bins())
         ostr << b.lowerBound() << "," << b.upperBound() << ",";
     ostr << "])";
     return ostr;
 }
 
-std::string Scale::unit() const
+Scale Scale::convertedScale(Unit newUnits) const
 {
-    if (m_name.empty())
-        return "";
-    size_t i = m_name.size() - 1;
-    if (m_name[i] != ')')
-        return "";
-    for (size_t j = i - 1; j > 0; --j)
-        if (m_name[j] == '(')
-            return m_name.substr(j + 1, i - j - 1);
-    return "";
-}
+    if (newUnits == m_units) {
+        return *this;
 
-Scale Scale::plottableScale() const
-{
-    std::string u = unit();
-    if (u == "rad") {
-        std::string outname = m_name.substr(0, m_name.size() - u.size() - 2) + "(deg)";
+    } else if (newUnits == Unit::nbins) {
         std::vector<Bin1D> outvector;
-        for (const Bin1D b : m_bins) {
-            double bmi = b.lowerBound() * 180 / pi;
-            double bma = b.upperBound() * 180 / pi;
-            outvector.emplace_back(Bin1D::FromTo(bmi, bma));
-        }
-        return Scale(outname, outvector);
+        for (size_t i = 0; i < m_bins.size(); i++)
+            outvector.emplace_back(Bin1D::At(i + 1, 0.5));
+        return Scale(m_name, newUnits, outvector);
+
+    } else if (Units::inOneGroup(m_units, newUnits)) {
+        auto f = [this, newUnits](double d) { return Units::convert(d, m_units, newUnits); };
+        std::vector<Bin1D> outvector = ::convertBins(m_bins, f);
+        return Scale(m_name, newUnits, outvector);
+
+    } else { // new units have have a different dimension (="from other units group")
+        // custom conversion function should be provided
+        // TODO temporary
+        return *this;
     }
     return *this;
 }
diff --git a/Base/Axis/Scale.h b/Base/Axis/Scale.h
index 824bb2d47fa6640774905273aec34129abdab823..206863a97eb2ddff8cd87a09b6fd60a2c523c33f 100644
--- a/Base/Axis/Scale.h
+++ b/Base/Axis/Scale.h
@@ -16,6 +16,7 @@
 #define BORNAGAIN_BASE_AXIS_SCALE_H
 
 #include "Base/Axis/Bin.h"
+#include "Base/Const/Units.h"
 #include <string>
 #include <utility>
 #include <vector>
@@ -24,15 +25,18 @@
 
 class Scale {
 public:
-    Scale(std::string name, const std::vector<Bin1D>& bins);
+    Scale(std::string name, Unit units, const std::vector<Bin1D>& bins);
     Scale* clone() const { return new Scale(*this); }
 
-    //! Sets the axis label
-    void setAxisName(std::string name) { m_name = name; }
-
-    //! Returns the label of the axis
+    //! Returns the name of the axis
     std::string axisName() const { return m_name; }
 
+    //! Returns the units in which cordinate values are given
+    Unit axisUnits() const { return m_units; }
+
+    //! Returns the label of the axis to display
+    std::string axisLabel() const;
+
     //... Getter functions: range
 
     //! Returns the number of bins
@@ -77,13 +81,12 @@ public:
 
     friend std::ostream& operator<<(std::ostream& ostr, const Scale& m);
 
-    std::string unit() const;
-
-    Scale plottableScale() const;
+    Scale convertedScale(Unit newUnits) const;
 
 protected:
     // private:
     std::string m_name;
+    Unit m_units;
     std::vector<Bin1D> m_bins;
 };
 
diff --git a/Base/Const/Units.cpp b/Base/Const/Units.cpp
index 9f415042c43e356fd7f63c14db27b9891a312b80..d46bae2e271de323e9e018445178a4675a0da18a 100644
--- a/Base/Const/Units.cpp
+++ b/Base/Const/Units.cpp
@@ -63,7 +63,9 @@ UnitInfo info(Unit unit)
     return {};
 }
 
-bool inOneGroup(Unit a, Unit b)
+} // namespace
+
+bool Units::inOneGroup(Unit a, Unit b)
 {
     for (auto group : infos) {
         bool hasA = false;
@@ -80,8 +82,6 @@ bool inOneGroup(Unit a, Unit b)
     return false;
 }
 
-} // namespace
-
 std::string Units::plainName(Unit unit)
 {
     return ::info(unit).plainName;
@@ -113,8 +113,8 @@ Unit Units::unitFromPlainName(std::string units)
 
 double Units::convert(double d, Unit from, Unit to)
 {
-    if (!::inOneGroup(from, to))
+    if (!inOneGroup(from, to))
         ASSERT(false);
 
-    return d * info(from).factor / info(to).factor;
+    return d * ::info(from).factor / ::info(to).factor;
 }
diff --git a/Base/Const/Units.h b/Base/Const/Units.h
index 5260c8896347191280fa9857d124613fb27207fc..68fe3283613645996f22c4f628c80f31a585a88f 100644
--- a/Base/Const/Units.h
+++ b/Base/Const/Units.h
@@ -58,6 +58,7 @@ namespace Units {
 static constexpr Unit internalAngularUnit = Unit::radian;
 static constexpr Unit internalDepthUnit = Unit::nanometer;
 
+bool inOneGroup(Unit a, Unit b);
 std::string plainName(Unit units);
 std::string symbolicName(Unit unit);
 bool unitExists(std::string units);
diff --git a/Base/Pixel/RectangularPixel.cpp b/Base/Pixel/RectangularPixel.cpp
index a98bd9f2891d6907ecf6c2e7df5cd99a3a3d7d4f..905807198377c076c0fb6a5d2b370d4be5fb7f75 100644
--- a/Base/Pixel/RectangularPixel.cpp
+++ b/Base/Pixel/RectangularPixel.cpp
@@ -83,5 +83,5 @@ Scale* RectangularPixel::createAxis(size_t n) const
     const double alpha_f_min = (pi / 2) - R3Util::theta(k00);
     const double alpha_f_max = (pi / 2) - R3Util::theta(k01);
 
-    return newEquiDivision("alpha_f (rad)", alpha_f_min, alpha_f_max, n);
+    return newEquiDivision("alpha_f", Unit::radian, alpha_f_min, alpha_f_max, n);
 }
diff --git a/Base/Py/PyFmt.cpp b/Base/Py/PyFmt.cpp
index d0185b4e51102302d7a0b712f40e4051729bf657..c4e1b5ab5fd064ca6c1b6378cbb8705bc4d732d9 100644
--- a/Base/Py/PyFmt.cpp
+++ b/Base/Py/PyFmt.cpp
@@ -204,3 +204,41 @@ std::string Py::Fmt::indent(size_t width)
 {
     return std::string(width, ' ');
 }
+
+std::string Py::Fmt::printUnit(Unit unit)
+{
+    switch (unit) {
+    case Unit::radian:
+        return "ba.Unit_radian";
+    case Unit::degree:
+        return "ba.Unit_degree";
+    case Unit::angstrom2:
+        return "ba.Unit_angstrom2";
+    case Unit::nanometer2:
+        return "ba.Unit_nanometer2";
+    case Unit::angstrom:
+        return "ba.Unit_angstrom";
+    case Unit::nanometer:
+        return "ba.Unit_nanometer";
+    case Unit::micrometer:
+        return "ba.Unit_micrometer";
+    case Unit::millimeter:
+        return "ba.Unit_millimeter";
+    case Unit::inv_angstrom:
+        return "ba.Unit_inv_angstrom";
+    case Unit::inv_nanometer:
+        return "ba.Unit_inv_nanometer";
+    case Unit::inv_angstrom2:
+        return "ba.Unit_inv_angstrom2";
+    case Unit::inv_nanometer2:
+        return "ba.Unit_inv_nanometer2";
+    case Unit::nbins:
+        return "ba.Unit_nbins";
+    case Unit::unitless:
+        return "ba.Unit_unitless";
+    case Unit::other:
+        return "ba.Unit_other";
+    default:
+        ASSERT(false);
+    }
+}
diff --git a/Base/Py/PyFmt.h b/Base/Py/PyFmt.h
index 7072572f5550f4aa431eb7c5af4a1fb9a0a97f78..b6812f8e61da6c049bcb15503a91ef6b1075ba1f 100644
--- a/Base/Py/PyFmt.h
+++ b/Base/Py/PyFmt.h
@@ -18,6 +18,7 @@
 #ifndef BORNAGAIN_BASE_PY_PYFMT_H
 #define BORNAGAIN_BASE_PY_PYFMT_H
 
+#include "Base/Const/Units.h"
 #include <heinz/Vectors3D.h>
 #include <string>
 #include <variant>
@@ -40,6 +41,7 @@ std::string printDegrees(double input);
 std::string printValue(double value, const std::string& units = "");
 std::string printValue(std::variant<double, int> value, const std::string& units = "");
 std::string printString(const std::string& value);
+std::string printUnit(Unit unit);
 
 //! Takes pairs of value/unit and concatenates them for an argument list.
 //! Each pair's content will be processed by printValue(), so the meaning
diff --git a/Device/Data/ArrayUtil.cpp b/Device/Data/ArrayUtil.cpp
index dd496f1198aef28783baecc1d8ba30062c0e54c7..09658662bd899e14227ef9d35b67226dec8816ec 100644
--- a/Device/Data/ArrayUtil.cpp
+++ b/Device/Data/ArrayUtil.cpp
@@ -33,7 +33,7 @@ std::pair<size_t, size_t> DataUtil::Array::getShape(const std::vector<std::vecto
 std::unique_ptr<Datafield> DataUtil::Array::createPField1D(const std::vector<double>& vec)
 {
     const size_t N = vec.size();
-    std::vector<const Scale*> axes{newEquiDivision("axis0", N, 0.0, (double)N)};
+    std::vector<const Scale*> axes{newEquiDivision("axis0", Unit::other, N, 0.0, (double)N)};
     return std::make_unique<Datafield>(std::move(axes), vec);
 }
 
@@ -47,8 +47,9 @@ DataUtil::Array::createPField2D(const std::vector<std::vector<double>>& vec)
     ASSERT(nrows > 0);
     ASSERT(ncols > 0);
 
-    std::vector<const Scale*> axes{newEquiDivision("axis0", ncols, 0.0, (double)ncols),
-                                   newEquiDivision("axis1", nrows, 0.0, (double)nrows)};
+    std::vector<const Scale*> axes{
+        newEquiDivision("axis0", Unit::other, ncols, 0.0, (double)ncols),
+        newEquiDivision("axis1", Unit::other, nrows, 0.0, (double)nrows)};
 
     std::vector<double> out(nrows * ncols);
     for (size_t row = 0; row < nrows; ++row) {
diff --git a/Device/Data/DataUtil.cpp b/Device/Data/DataUtil.cpp
index c51da04e041f677c2e1c96034c42c9957d195ae8..dfc50decd99fa7df5d9adcc259680b9fab96847c 100644
--- a/Device/Data/DataUtil.cpp
+++ b/Device/Data/DataUtil.cpp
@@ -144,8 +144,8 @@ Datafield DataUtil::Data::vecvecToDatafield(const std::vector<std::vector<double
     size_t nrows = array_2d.size();
     size_t ncols = array_2d[0].size();
 
-    std::vector<const Scale*> axes{newEquiDivision("x", nrows, 0.0, double(nrows)),
-                                   newEquiDivision("y", ncols, 0.0, double(ncols))};
+    std::vector<const Scale*> axes{newEquiDivision("x", Unit::other, nrows, 0.0, double(nrows)),
+                                   newEquiDivision("y", Unit::other, ncols, 0.0, double(ncols))};
     std::vector<double> out;
     out.reserve(nrows * ncols);
     for (size_t row = 0; row < nrows; row++) {
diff --git a/Device/Data/Datafield.cpp b/Device/Data/Datafield.cpp
index 7e1673566f27f1bf46c2b37aa0d57d52101d4edb..09f12f2232c77c933228cffaad3ee12b6e4b73d0 100644
--- a/Device/Data/Datafield.cpp
+++ b/Device/Data/Datafield.cpp
@@ -329,9 +329,21 @@ Datafield* Datafield::create_yProjection(int xbinlow, int xbinup) const
     return new Datafield({yAxis().clone()}, out);
 }
 
+// Legacy function for default conversion. Used in Python modules.
 Datafield Datafield::plottableField() const
 {
-    return {title(), frame().plottableFrame(), m_values, m_errSigmas};
+    std::vector<Unit> newUnits;
+    for (size_t i = 0; i < frame().rank(); i++)
+        if (frame().axis(i).axisUnits() == Unit::radian)
+            newUnits.push_back(Unit::degree);
+        else
+            newUnits.push_back(frame().axis(i).axisUnits());
+    return plottableField(newUnits);
+}
+
+Datafield Datafield::plottableField(const std::vector<Unit> &newUnits) const
+{
+    return {title(), frame().convertedFrame(newUnits), m_values, m_errSigmas};
 }
 
 Datafield Datafield::flat() const
diff --git a/Device/Data/Datafield.h b/Device/Data/Datafield.h
index b0f0dbc921cd4b96808329540e25fca11d1eecc7..f80c5506c887927d67cfb41050adf6f14b82a1b7 100644
--- a/Device/Data/Datafield.h
+++ b/Device/Data/Datafield.h
@@ -15,6 +15,7 @@
 #ifndef BORNAGAIN_DEVICE_DATA_DATAFIELD_H
 #define BORNAGAIN_DEVICE_DATA_DATAFIELD_H
 
+#include "Base/Const/Units.h"
 #ifdef BORNAGAIN_PYTHON
 #include "PyCore/Embed/PyObjectDecl.h"
 #endif
@@ -82,6 +83,7 @@ public:
     //... modifiers
 
     Datafield plottableField() const;
+    Datafield plottableField(const std::vector<Unit>& newUnits) const;
     Datafield flat() const;
 
     //! Multiplies contents by constant factor, e.g. to shift curves in log plot
diff --git a/Device/Detector/OffspecDetector.cpp b/Device/Detector/OffspecDetector.cpp
index 45ebef085369b075b1dfd422be5785053bf55a1a..b882625b59b2fd215fe7d223048d0a40e783b4f7 100644
--- a/Device/Detector/OffspecDetector.cpp
+++ b/Device/Detector/OffspecDetector.cpp
@@ -21,8 +21,9 @@
 
 OffspecDetector::OffspecDetector(size_t n_phi, double phi_min, double phi_max, size_t n_alpha,
                                  double alpha_min, double alpha_max)
-    : m_axes{sharedEquiDivision("phi_f (rad)", n_phi, phi_min, phi_max),
-             sharedEquiDivision("alpha_f (rad)", n_alpha, alpha_min, alpha_max)}
+    : m_axes{
+        sharedEquiDivision("phi_f", Units::internalAngularUnit, n_phi, phi_min, phi_max),
+        sharedEquiDivision("alpha_f", Units::internalAngularUnit, n_alpha, alpha_min, alpha_max)}
 {
 }
 
diff --git a/Device/Detector/RectangularDetector.cpp b/Device/Detector/RectangularDetector.cpp
index f30d7655759046323561759fa66a176ac3829011..2238f0eb0aff4c393d98074a43e49f11ab247f1f 100644
--- a/Device/Detector/RectangularDetector.cpp
+++ b/Device/Detector/RectangularDetector.cpp
@@ -39,7 +39,8 @@ RectangularDetector::RectangularDetector(std::array<std::shared_ptr<Scale>, 2> a
 
 RectangularDetector::RectangularDetector(size_t nxbins, double width, size_t nybins, double height)
     : RectangularDetector(std::array<std::shared_ptr<Scale>, 2>{
-        sharedEquiDivision("u", nxbins, 0.0, width), sharedEquiDivision("v", nybins, 0.0, height)})
+        sharedEquiDivision("u", Unit::millimeter, nxbins, 0.0, width),
+        sharedEquiDivision("v", Unit::millimeter, nybins, 0.0, height)})
 {
 }
 
diff --git a/Device/Detector/SphericalDetector.cpp b/Device/Detector/SphericalDetector.cpp
index 9e42dcf179948a8fd43833c6b9973cd886dfb089..09326d0fb75a3a56d444113d7bd0a7963d5b5491 100644
--- a/Device/Detector/SphericalDetector.cpp
+++ b/Device/Detector/SphericalDetector.cpp
@@ -33,8 +33,8 @@ SphericalDetector::SphericalDetector(std::array<std::shared_ptr<Scale>, 2> axes)
 SphericalDetector::SphericalDetector(size_t n_phi, double phi_min, double phi_max, size_t n_alpha,
                                      double alpha_min, double alpha_max)
     : SphericalDetector(std::array<std::shared_ptr<Scale>, 2>{
-        sharedEquiDivision("phi_f (rad)", n_phi, phi_min, phi_max),
-        sharedEquiDivision("alpha_f (rad)", n_alpha, alpha_min, alpha_max)})
+        sharedEquiDivision("phi_f", Units::internalAngularUnit, n_phi, phi_min, phi_max),
+        sharedEquiDivision("alpha_f", Units::internalAngularUnit, n_alpha, alpha_min, alpha_max)})
 {
 }
 
diff --git a/Device/IO/ParseUtil.cpp b/Device/IO/ParseUtil.cpp
index e932c5579a072967d5b0e1ed11b6fe98f8bb7115..7688bf97e5eb493eb5025639a0f23f9608aa67d0 100644
--- a/Device/IO/ParseUtil.cpp
+++ b/Device/IO/ParseUtil.cpp
@@ -14,6 +14,7 @@
 
 #include "Device/IO/ParseUtil.h"
 #include "Base/Axis/MakeScale.h"
+#include "Base/Const/Units.h"
 #include "Base/Util/PathUtil.h"
 #include "Base/Util/StringUtil.h"
 #include "Device/Data/Datafield.h"
@@ -57,6 +58,19 @@ std::vector<double> parse_x_list(std::string text, const std::string& type)
     return result;
 }
 
+std::string unitsFromLabel(std::string name)
+{
+    if (name.empty())
+        return "";
+    size_t i = name.size() - 1;
+    if (name[i] != ')')
+        return "";
+    for (size_t j = i - 1; j > 0; --j)
+        if (name[j] == '(')
+            return name.substr(j + 1, i - j - 1);
+    return "";
+}
+
 } // namespace
 
 //! Creates axis of certain type from input stream
@@ -69,8 +83,14 @@ Scale* Util::Parse::parseScale(std::istream& input_stream)
     if (matched.size() != 4)
         throw std::runtime_error("Cannot read axis from input");
     const std::string type = matched[1];
-    const std::string name = matched[2];
+    const std::string label = matched[2];
     const std::string body = matched[3];
+    const std::string units = ::unitsFromLabel(label);
+    const std::string name = label.substr(0, label.size() - units.size() - 3);
+
+    // files from pre-21 have no units information --> set units to default and update them later
+    Unit axisUnits = Units::unitExists(units) ? Units::unitFromPlainName(units) : Unit::other;
+
     if (type == "EquiDivision" || type == "FixedBinAxis" /* for compatibility with pre-21 */) {
         std::vector<std::string> arr = Base::String::split(body, ",");
         int nbins;
@@ -78,13 +98,13 @@ Scale* Util::Parse::parseScale(std::istream& input_stream)
         if (!(Base::String::to_int(arr[0], &nbins) && Base::String::to_double(arr[1], &xmi)
               && Base::String::to_double(arr[2], &xma)))
             throw std::runtime_error("Reading EquiDivision: cannot read parameters");
-        return newEquiDivision(name, nbins, xmi, xma);
+        return newEquiDivision(name, axisUnits, nbins, xmi, xma);
 
     } else if (type == "ListScan" || type == "DiscreteAxis" /* for compatibility with pre-21 */) {
-        return newListScan(name, parse_x_list(body, type));
+        return newListScan(name, axisUnits, parse_x_list(body, type));
 
     } else if (type == "Scale") {
-        return newGenericScale(name, parse_x_list(body, type));
+        return newGenericScale(name, axisUnits, parse_x_list(body, type));
 
     } else
         throw std::runtime_error("Unknown axis type '" + type + "'");
diff --git a/Device/IO/ReadReflectometry.cpp b/Device/IO/ReadReflectometry.cpp
index 2ccb4c4a8a45d8cefdc52f2ea9857e4814a2424c..dc7483b6494223cee7cd6e444508e41cd2b7cd22 100644
--- a/Device/IO/ReadReflectometry.cpp
+++ b/Device/IO/ReadReflectometry.cpp
@@ -57,7 +57,8 @@ Datafield* Util::RW::readReflectometryTable(std::istream& s, const ImportSetting
             sRVec.push_back(rowVec[p.col_sR - 1]);
     }
 
-    return new Datafield(std::vector<const Scale*>{newListScan("qVector", QVec)}, RVec, sRVec);
+    return new Datafield(
+        std::vector<const Scale*>{newListScan("qVector", Unit::inv_nanometer, QVec)}, RVec, sRVec);
 }
 
 Datafield* Util::RW::readMotofit(std::istream& s)
@@ -117,5 +118,6 @@ Datafield* Util::RW::readMotofit(std::istream& s)
             throw std::runtime_error("Error in data line");
     }
 
-    return new Datafield(std::vector<const Scale*>{newListScan("qVector", Q)}, R, wR);
+    return new Datafield(std::vector<const Scale*>{newListScan("qVector", Unit::inv_nanometer, Q)},
+                         R, wR);
 }
diff --git a/Device/IO/ReadRefsans.cpp b/Device/IO/ReadRefsans.cpp
index b65d40bfb0917c64103c07eccc3695c23b8265b2..eca0dfb2029d89c1bf41f7545250704c1c9a5895 100644
--- a/Device/IO/ReadRefsans.cpp
+++ b/Device/IO/ReadRefsans.cpp
@@ -90,7 +90,7 @@ Datafield* Util::RW::readRefsans(std::istream& input_stream)
 
     // create datafield
     std::vector<const Scale*> axes;
-    axes.push_back(newListScan("Qy (1/nm)", qy));
-    axes.push_back(newListScan("Qz (1/nm)", qz));
+    axes.push_back(newListScan("Qy", Unit::inv_nanometer, qy));
+    axes.push_back(newListScan("Qz", Unit::inv_nanometer, qz));
     return new Datafield(std::move(axes), values);
 }
diff --git a/Device/IO/ReadWriteNicos.cpp b/Device/IO/ReadWriteNicos.cpp
index cbfff5d83bd78334bd8b7341850bc4990f9ee78b..1adf37dc40e372231999610de7e1ac8a6136a5af 100644
--- a/Device/IO/ReadWriteNicos.cpp
+++ b/Device/IO/ReadWriteNicos.cpp
@@ -101,8 +101,8 @@ Datafield* Util::RW::readNicos(std::istream& input_stream)
     if (height == 0)
         throw std::runtime_error("Could not find 'DataSizeY' value");
 
-    std::vector<const Scale*> axes{newEquiDivision("x", width, 0.0, width),
-                                   newEquiDivision("y", height, 0.0, height)};
+    std::vector<const Scale*> axes{newEquiDivision("x", Unit::nbins, width, 0.0, width),
+                                   newEquiDivision("y", Unit::nbins, height, 0.0, height)};
     auto result = std::make_unique<Datafield>(std::move(axes));
 
     // -- read data
diff --git a/Device/IO/ReadWriteTiff.cpp b/Device/IO/ReadWriteTiff.cpp
index a23a07847a37727c19bdf1ef3ba1217a4160790f..0578e323215edff7f0bae4aebbb67bf66eeede27 100644
--- a/Device/IO/ReadWriteTiff.cpp
+++ b/Device/IO/ReadWriteTiff.cpp
@@ -112,8 +112,8 @@ Datafield* Util::RW::readTiff(std::istream& input_stream)
         throw std::runtime_error("Cannot read TIFF file: failed allocating buffer");
 
     auto data = std::make_unique<Datafield>(
-        std::vector<const Scale*>{newEquiDivision("x", width, 0.0, double(width)),
-                                  newEquiDivision("y", height, 0.0, double(height))});
+        std::vector<const Scale*>{newEquiDivision("x", Unit::nbins, width, 0.0, double(width)),
+                                  newEquiDivision("y", Unit::nbins, height, 0.0, double(height))});
 
     std::vector<int8_t> line_buf;
     line_buf.resize(buf_size, 0);
diff --git a/GUI/Model/Axis/AmplitudeAxisItem.cpp b/GUI/Model/Axis/AmplitudeAxisItem.cpp
index fb1d9b566ae12af47f03a1559bbac752913112e9..8408fe0c5cf480667c38172cce04a6cfed90286f 100644
--- a/GUI/Model/Axis/AmplitudeAxisItem.cpp
+++ b/GUI/Model/Axis/AmplitudeAxisItem.cpp
@@ -30,8 +30,8 @@ const QString LogRangeOrders("LogRangeOrders");
 } // namespace Tag
 } // namespace
 
-AmplitudeAxisItem::AmplitudeAxisItem(QObject* parent)
-    : BasicAxisItem(parent)
+AmplitudeAxisItem::AmplitudeAxisItem(Unit units, QObject* parent)
+    : BasicAxisItem(units, parent)
     , m_lockMinMax(false)
     , m_logScale(true)
     , m_logRangeOrders(6.0)
diff --git a/GUI/Model/Axis/AmplitudeAxisItem.h b/GUI/Model/Axis/AmplitudeAxisItem.h
index a0400698239d30f302defda5e96501aa4816d5c6..a23fdc037c390ab99377fc303b6bc1f184b4d186 100644
--- a/GUI/Model/Axis/AmplitudeAxisItem.h
+++ b/GUI/Model/Axis/AmplitudeAxisItem.h
@@ -21,7 +21,7 @@ class AmplitudeAxisItem : public BasicAxisItem {
     Q_OBJECT
 
 public:
-    explicit AmplitudeAxisItem(QObject* parent = nullptr);
+    explicit AmplitudeAxisItem(Unit units, QObject* parent = nullptr);
 
     void writeTo(QXmlStreamWriter* w) const override;
     void readFrom(QXmlStreamReader* r) override;
diff --git a/GUI/Model/Axis/BasicAxisItem.cpp b/GUI/Model/Axis/BasicAxisItem.cpp
index b8ceb3ae70751541bd6b611e3eb0d2fa7d6641df..11b318da71a46e89cc93589b89209e21da00a6e6 100644
--- a/GUI/Model/Axis/BasicAxisItem.cpp
+++ b/GUI/Model/Axis/BasicAxisItem.cpp
@@ -15,13 +15,13 @@
 #include "GUI/Model/Axis/BasicAxisItem.h"
 #include "Base/Axis/MakeScale.h"
 #include "Base/Axis/Scale.h"
-#include "Base/Const/Units.h"
 #include "GUI/Support/XML/UtilXML.h"
 
 namespace {
 namespace Tag {
 
 const QString IsVisible("IsVisible");
+const QString Tag_Units("Units");
 const QString Nbins("Nbins");
 const QString MinDeg("MinDeg");
 const QString MaxDeg("MaxDeg");
@@ -32,9 +32,10 @@ const QString LogScale("LogScale");
 } // namespace Tag
 } // namespace
 
-BasicAxisItem::BasicAxisItem(QObject* parent)
+BasicAxisItem::BasicAxisItem(Unit units, QObject* parent)
     : QObject(parent)
     , m_visible(true)
+    , m_units(units)
     , m_nbins(100)
     , m_min(0)
     , m_max(-1)
@@ -78,7 +79,7 @@ void BasicAxisItem::setMax(double value)
 
 Scale BasicAxisItem::makeScale(std::string name) const
 {
-    return EquiScan(name, binCount(), Units::deg2rad(min()), Units::deg2rad(max()));
+    return EquiScan(name, m_units, m_nbins, m_min, m_max);
 }
 
 bool BasicAxisItem::isVisible() const
@@ -94,13 +95,19 @@ void BasicAxisItem::setVisible(bool b)
 
 void BasicAxisItem::writeTo(QXmlStreamWriter* w) const
 {
-    XML::writeAttribute(w, XML::Attrib::version, uint(1));
+    XML::writeAttribute(w, XML::Attrib::version, uint(2));
 
     // visibility
     w->writeStartElement(Tag::IsVisible);
     XML::writeAttribute(w, XML::Attrib::value, m_visible);
     w->writeEndElement();
 
+    // units
+    w->writeStartElement(Tag::Tag_Units);
+    QString units = QString::fromStdString(Units::plainName(m_units));
+    XML::writeAttribute(w, XML::Attrib::value, units);
+    w->writeEndElement();
+
     // nbins
     w->writeStartElement(Tag::Nbins);
     XML::writeAttribute(w, XML::Attrib::value, m_nbins);
@@ -120,7 +127,6 @@ void BasicAxisItem::writeTo(QXmlStreamWriter* w) const
 void BasicAxisItem::readFrom(QXmlStreamReader* r)
 {
     const uint version = XML::readUIntAttribute(r, XML::Attrib::version);
-    Q_UNUSED(version)
 
     while (r->readNextStartElement()) {
         QString tag = r->name().toString();
@@ -130,6 +136,13 @@ void BasicAxisItem::readFrom(QXmlStreamReader* r)
             XML::readAttribute(r, XML::Attrib::value, &m_visible);
             XML::gotoEndElementOfTag(r, tag);
 
+            // units
+        } else if (version == 2 && tag == Tag::Tag_Units) {
+            QString units;
+            XML::readAttribute(r, XML::Attrib::value, &units);
+            m_units = Units::unitFromPlainName(units.toStdString());
+            XML::gotoEndElementOfTag(r, tag);
+
             // nbins
         } else if (tag == Tag::Nbins) {
             XML::readAttribute(r, XML::Attrib::value, &m_nbins);
diff --git a/GUI/Model/Axis/BasicAxisItem.h b/GUI/Model/Axis/BasicAxisItem.h
index 66de70e31b3b8480f8f46da544a906698e06cf74..03ef877b9d0235c38d6527f3d2775eb26d42ae10 100644
--- a/GUI/Model/Axis/BasicAxisItem.h
+++ b/GUI/Model/Axis/BasicAxisItem.h
@@ -15,6 +15,7 @@
 #ifndef BORNAGAIN_GUI_MODEL_AXIS_BASICAXISITEM_H
 #define BORNAGAIN_GUI_MODEL_AXIS_BASICAXISITEM_H
 
+#include "Base/Const/Units.h"
 #include <QObject>
 
 class Scale;
@@ -25,7 +26,7 @@ class BasicAxisItem : public QObject {
     Q_OBJECT
 
 public:
-    explicit BasicAxisItem(QObject* parent = nullptr);
+    explicit BasicAxisItem(Unit units, QObject* parent = nullptr);
     ~BasicAxisItem() override;
 
     virtual void writeTo(QXmlStreamWriter* w) const;
@@ -34,6 +35,9 @@ public:
 public:
     Scale makeScale(std::string name) const;
 
+    Unit units() const { return m_units; }
+    void setUnits(Unit u) { m_units = u; }
+
     int binCount() const;
     void setBinCount(size_t value);
 
@@ -53,6 +57,7 @@ signals:
 
 private:
     bool m_visible;
+    Unit m_units;
     int m_nbins;
     double m_min; // in display units (if angle, then deg)
     double m_max;
diff --git a/GUI/Model/Axis/PointwiseAxisItem.cpp b/GUI/Model/Axis/PointwiseAxisItem.cpp
index 3f3631c180dc930292e64cb84b15893cb3532c5e..dc50a06c598ef11cf117ec9579d2ca8d5cd6ede8 100644
--- a/GUI/Model/Axis/PointwiseAxisItem.cpp
+++ b/GUI/Model/Axis/PointwiseAxisItem.cpp
@@ -32,7 +32,7 @@ const QString BaseData("BaseData");
 } // namespace
 
 PointwiseAxisItem::PointwiseAxisItem(QObject* parent)
-    : BasicAxisItem(parent)
+    : BasicAxisItem(Unit::nbins, parent)
     , m_nativeAxisUnits("nbins")
 {
 }
diff --git a/GUI/Model/Beam/GrazingScanItem.cpp b/GUI/Model/Beam/GrazingScanItem.cpp
index be38e536f5713140be19b17ff6650ed55e7b3387..7123e6d4b40c14f0da0497e761aa450617c5d954 100644
--- a/GUI/Model/Beam/GrazingScanItem.cpp
+++ b/GUI/Model/Beam/GrazingScanItem.cpp
@@ -68,7 +68,7 @@ GrazingScanItem::GrazingScanItem()
     m_distribution.initWithInitializer("Distribution", "",
                                        DistributionItemCatalog::symmetricTypes(), initDistribution);
 
-    m_uniformAlphaAxis.reset(new BasicAxisItem());
+    m_uniformAlphaAxis.reset(new BasicAxisItem(Unit::degree));
     setAxisPresentationDefaults(m_uniformAlphaAxis.get());
     m_currentAxisIsUniformAxis = true;
 }
@@ -123,7 +123,7 @@ void GrazingScanItem::readFrom(QXmlStreamReader* r)
 
             // uniform axis
         } else if (tag == Tag::UniformAxis) {
-            m_uniformAlphaAxis = std::make_unique<BasicAxisItem>();
+            m_uniformAlphaAxis = std::make_unique<BasicAxisItem>(Unit::degree);
             setAxisPresentationDefaults(m_uniformAlphaAxis.get());
             m_uniformAlphaAxis->readFrom(r);
             XML::gotoEndElementOfTag(r, tag);
diff --git a/GUI/Model/Data/DataItem.cpp b/GUI/Model/Data/DataItem.cpp
index 8b5cf51af85f94423f9924bac8dea5fa78ff9c98..67fd79cb8c1d3e150fcb08c46d2ef50434350669 100644
--- a/GUI/Model/Data/DataItem.cpp
+++ b/GUI/Model/Data/DataItem.cpp
@@ -35,14 +35,11 @@ const QString YAxis("YAxis");
 } // namespace Tag
 } // namespace
 
-const QString x_axis_default_name = "X [nbins]";
-const QString y_axis_default_name = "Y [nbins]";
-
 DataItem::DataItem(const QString& modelType)
     : TYPE(modelType)
     , m_fileName("undefined")
-    , m_xAxis(std::make_unique<BasicAxisItem>())
-    , m_yAxis(std::make_unique<AmplitudeAxisItem>())
+    , m_xAxis(std::make_unique<BasicAxisItem>(Unit::nbins))
+    , m_yAxis(std::make_unique<AmplitudeAxisItem>(Unit::nbins))
     , m_last_modified(QDateTime::currentDateTime())
     , m_last_saved(m_last_modified.addSecs(-2023))
 {
@@ -53,7 +50,10 @@ DataItem::~DataItem() = default;
 void DataItem::setDatafield(const Datafield& data)
 {
     std::unique_lock<std::mutex> lock(m_update_data_mutex);
-    m_datafield = std::make_unique<Datafield>(data.plottableField());
+    std::vector<Unit> units = {xAxisItem()->units()};
+    if (data.rank() == 2)
+        units.push_back(yAxisItem()->units());
+    m_datafield = std::make_unique<Datafield>(data.plottableField(units));
     setLastModified(QDateTime::currentDateTime());
     emit datafieldChanged();
 }
@@ -223,6 +223,17 @@ void DataItem::copyXYRangesFromItem(DataItem* sourceItem)
     copyYRangeFromItem(sourceItem);
 }
 
+void DataItem::copyUnitsFromItem(DataItem* sourceItem)
+{
+    m_xAxis->setUnits(sourceItem->xAxisItem()->units());
+    m_yAxis->setUnits(sourceItem->yAxisItem()->units());
+}
+
+QString DataItem::currentAxesUnits()
+{
+    return QString::fromStdString(Units::plainName(xAxisItem()->units()));
+}
+
 const BasicAxisItem* DataItem::xAxisItem() const
 {
     return m_xAxis.get();
@@ -245,12 +256,12 @@ AmplitudeAxisItem* DataItem::yAxisItem()
 
 QString DataItem::XaxisTitle() const
 {
-    return m_datafield ? QString::fromStdString(m_datafield->xAxis().axisName()) : "";
+    return m_datafield ? QString::fromStdString(m_datafield->xAxis().axisLabel()) : "";
 }
 
 QString DataItem::YaxisTitle() const
 {
-    return m_datafield ? QString::fromStdString(m_datafield->yAxis().axisName()) : "";
+    return m_datafield ? QString::fromStdString(m_datafield->yAxis().axisLabel()) : "";
 }
 
 void DataItem::setAxesRangeToData()
diff --git a/GUI/Model/Data/DataItem.h b/GUI/Model/Data/DataItem.h
index b4c2a231961e227612ef7bef778e88acf6332496..25b2f94aea0570e81d1813e41949033076c74c3b 100644
--- a/GUI/Model/Data/DataItem.h
+++ b/GUI/Model/Data/DataItem.h
@@ -15,6 +15,7 @@
 #ifndef BORNAGAIN_GUI_MODEL_DATA_DATAITEM_H
 #define BORNAGAIN_GUI_MODEL_DATA_DATAITEM_H
 
+#include "Base/Const/Units.h"
 #include "GUI/Model/Data/ComboProperty.h"
 #include <QDateTime>
 #include <mutex>
@@ -91,6 +92,9 @@ public:
     void copyYRangeFromItem(DataItem* sourceItem);
     void copyXYRangesFromItem(DataItem* sourceItem);
 
+    void copyUnitsFromItem(DataItem* sourceItem);
+    QString currentAxesUnits();
+
     const BasicAxisItem* xAxisItem() const;
     BasicAxisItem* xAxisItem();
     const AmplitudeAxisItem* yAxisItem() const;
diff --git a/GUI/Model/Data/IntensityDataItem.cpp b/GUI/Model/Data/IntensityDataItem.cpp
index 05072fcada5b1d95c40fc30419cac77976a9bb2f..1ec41e0d95ecadb2d239b1baee34c260f2740350 100644
--- a/GUI/Model/Data/IntensityDataItem.cpp
+++ b/GUI/Model/Data/IntensityDataItem.cpp
@@ -65,7 +65,7 @@ IntensityDataItem::IntensityDataItem()
     , m_isInterpolated(true)
     , m_gradient(
           ComboProperty::fromList(gradient_map.keys() + custom_gradient_map.keys(), startGradient))
-    , m_zAxis(std::make_unique<AmplitudeAxisItem>())
+    , m_zAxis(std::make_unique<AmplitudeAxisItem>(Unit::unitless))
 {
 }
 
diff --git a/GUI/Model/Descriptor/AxisProperty.cpp b/GUI/Model/Descriptor/AxisProperty.cpp
index 91da7727632cf2f9481167a696e138275c93c5b4..1ce0e54d758876b14ce2f3e81ad882bdfa6e8368 100644
--- a/GUI/Model/Descriptor/AxisProperty.cpp
+++ b/GUI/Model/Descriptor/AxisProperty.cpp
@@ -13,12 +13,14 @@
 //  ************************************************************************************************
 
 #include "GUI/Model/Descriptor/AxisProperty.h"
+#include "Base/Axis/MakeScale.h"
 #include "Base/Axis/Scale.h"
 #include "GUI/Support/XML/UtilXML.h"
 
 namespace {
 namespace Tag {
 
+const QString Tag_Units("Units");
 const QString Nbins("Nbins");
 const QString Min("Min");
 const QString Max("Max");
@@ -29,7 +31,11 @@ const QString ExpandGroupbox("ExpandGroupbox");
 
 using std::variant;
 
-AxisProperty::AxisProperty() = default;
+
+AxisProperty::AxisProperty(Unit units)
+    : m_units(units)
+{
+}
 
 void AxisProperty::initMin(const QString& label, const QString& tooltip, double value,
                            const variant<QString, Unit>& unit, const RealLimits& limits,
@@ -45,14 +51,20 @@ void AxisProperty::initMax(const QString& label, const QString& tooltip, double
     m_max.init(label, tooltip, value, unit, decimals, limits, "max");
 }
 
-Scale AxisProperty::createAxis(const std::string& name) const
+Scale AxisProperty::makeScale(const std::string& name) const
 {
-    return EquiDivision(name, m_nbins, m_min, m_max);
+    return EquiDivision(name, m_units, m_nbins, m_min, m_max);
 }
 
 void AxisProperty::writeTo(QXmlStreamWriter* w) const
 {
-    XML::writeAttribute(w, XML::Attrib::version, uint(1));
+    XML::writeAttribute(w, XML::Attrib::version, uint(2));
+
+    // units
+    w->writeStartElement(Tag::Tag_Units);
+    QString units = QString::fromStdString(Units::plainName(m_units));
+    XML::writeAttribute(w, XML::Attrib::value, units);
+    w->writeEndElement();
 
     // nbins
     w->writeStartElement(Tag::Nbins);
@@ -83,8 +95,15 @@ void AxisProperty::readFrom(QXmlStreamReader* r)
     while (r->readNextStartElement()) {
         QString tag = r->name().toString();
 
-        // nbins
-        if (tag == Tag::Nbins) {
+        // units
+        if (version == 2 && tag == Tag::Tag_Units) {
+            QString units;
+            XML::readAttribute(r, XML::Attrib::value, &units);
+            m_units = Units::unitFromPlainName(units.toStdString());
+            XML::gotoEndElementOfTag(r, tag);
+
+            // nbins
+        } else if (tag == Tag::Nbins) {
             XML::readAttribute(r, XML::Attrib::value, &m_nbins);
             XML::gotoEndElementOfTag(r, tag);
 
diff --git a/GUI/Model/Descriptor/AxisProperty.h b/GUI/Model/Descriptor/AxisProperty.h
index 249db1d7470f890979355f1830fd08630296c4c2..2d3d3e8e54f6a360428562a4186ce49b716a3aef 100644
--- a/GUI/Model/Descriptor/AxisProperty.h
+++ b/GUI/Model/Descriptor/AxisProperty.h
@@ -15,10 +15,11 @@
 #ifndef BORNAGAIN_GUI_MODEL_DESCRIPTOR_AXISPROPERTY_H
 #define BORNAGAIN_GUI_MODEL_DESCRIPTOR_AXISPROPERTY_H
 
-#include "Base/Axis/MakeScale.h"
 #include "GUI/Model/Descriptor/DoubleProperty.h"
 #include <memory>
 
+class Scale;
+
 //! Holds values which can be used to describe a EquiDivision.
 //!
 //! Use this as a member in your class to
@@ -34,7 +35,7 @@
 
 class AxisProperty {
 public:
-    explicit AxisProperty();
+    explicit AxisProperty(Unit units);
 
     DoubleProperty& min() { return m_min; }
     const DoubleProperty& min() const { return m_min; }
@@ -51,7 +52,10 @@ public:
                  const std::variant<QString, Unit>& unit,
                  const RealLimits& limits = RealLimits::nonnegative(), uint decimals = 3);
 
-    Scale createAxis(const std::string& name) const;
+    Scale makeScale(const std::string& name) const;
+
+    Unit units() const { return m_units; }
+    void setUnits(Unit v) { m_units = v; }
 
     uint nbins() const { return m_nbins; }
     void setNbins(uint v) { m_nbins = v; }
@@ -63,6 +67,7 @@ public:
     void readFrom(QXmlStreamReader* r);
 
 private:
+    Unit m_units;
     uint m_nbins = 100;
     DoubleProperty m_min;
     DoubleProperty m_max;
diff --git a/GUI/Model/Detector/OffspecDetectorItem.cpp b/GUI/Model/Detector/OffspecDetectorItem.cpp
index f5c127330189a6dcd9f91f9fb02025e23c4a3f25..da105df74a4ed25ecc94e073c52866c364d0c031 100644
--- a/GUI/Model/Detector/OffspecDetectorItem.cpp
+++ b/GUI/Model/Detector/OffspecDetectorItem.cpp
@@ -26,18 +26,22 @@ const QString BaseData("BaseData");
 
 } // namespace Tag
 
+constexpr Unit offspecDetUnit = Unit::degree;
+
 } // namespace
 
 OffspecDetectorItem::OffspecDetectorItem()
+    : m_phiAxis(offspecDetUnit)
+    , m_alphaAxis(offspecDetUnit)
 {
-    m_phiAxis.initMin("Min", "Lower edge of first phi-bin", -1.0, Unit::degree,
+    m_phiAxis.initMin("Min", "Lower edge of first phi-bin", -1.0, m_phiAxis.units(),
                       RealLimits::limited(-90, 90));
-    m_phiAxis.initMax("Max", "Upper edge of last phi-bin", 1.0, Unit::degree,
+    m_phiAxis.initMax("Max", "Upper edge of last phi-bin", 1.0, m_phiAxis.units(),
                       RealLimits::limited(-90, 90));
 
-    m_alphaAxis.initMin("Min", "Lower edge of first alpha-bin", 0.0, Unit::degree,
+    m_alphaAxis.initMin("Min", "Lower edge of first alpha-bin", 0.0, m_alphaAxis.units(),
                         RealLimits::limited(-90, 90));
-    m_alphaAxis.initMax("Max", "Upper edge of last alpha-bin", 2.0, Unit::degree,
+    m_alphaAxis.initMax("Max", "Upper edge of last alpha-bin", 2.0, m_alphaAxis.units(),
                         RealLimits::limited(-90, 90));
 }
 
@@ -82,12 +86,16 @@ void OffspecDetectorItem::readFrom(QXmlStreamReader* r)
 std::unique_ptr<OffspecDetector> OffspecDetectorItem::createOffspecDetector() const
 {
     const int n_x = m_phiAxis.nbins();
-    const double x_min = Units::deg2rad(m_phiAxis.min());
-    const double x_max = Units::deg2rad(m_phiAxis.max());
+    const double x_min =
+        Units::convert(m_phiAxis.min(), m_phiAxis.units(), Units::internalAngularUnit);
+    const double x_max =
+        Units::convert(m_phiAxis.max(), m_phiAxis.units(), Units::internalAngularUnit);
 
     const int n_y = m_alphaAxis.nbins();
-    const double y_min = Units::deg2rad(m_alphaAxis.min());
-    const double y_max = Units::deg2rad(m_alphaAxis.max());
+    const double y_min =
+        Units::convert(m_alphaAxis.min(), m_alphaAxis.units(), Units::internalAngularUnit);
+    const double y_max =
+        Units::convert(m_alphaAxis.max(), m_alphaAxis.units(), Units::internalAngularUnit);
 
     return std::make_unique<OffspecDetector>(n_x, x_min, x_max, n_y, y_min, y_max);
 }
diff --git a/GUI/Model/Detector/RectangularDetectorItem.cpp b/GUI/Model/Detector/RectangularDetectorItem.cpp
index 4bd41457c64e0e60c125aae47229c803188ee039..0380612be3195a0478d06fc497389206a84dc759 100644
--- a/GUI/Model/Detector/RectangularDetectorItem.cpp
+++ b/GUI/Model/Detector/RectangularDetectorItem.cpp
@@ -44,9 +44,11 @@ const QMap<RectangularDetector::EDetectorArrangement, QString> alignment_names_m
     {RectangularDetector::PERPENDICULAR_TO_DIRECT_BEAM, "Perpendicular to direct beam"},
     {RectangularDetector::PERPENDICULAR_TO_REFLECTED_BEAM, "Perpendicular to reflected beam"}};
 
+constexpr Unit rectDetUnit = Unit::millimeter;
+
 void initResolutionFunction(ResolutionFunctionItem* newFunc, const ResolutionFunctionItem*)
 {
-    newFunc->setUnit("mm");
+    newFunc->setUnit(rectDetUnit);
 }
 
 } // namespace
@@ -58,9 +60,9 @@ RectangularDetectorItem::RectangularDetectorItem()
 
     m_xSize = 100;
     m_ySize = 100;
-    m_width.init("Width (x-axis)", "Width of the detector", default_detector_width, "mm",
+    m_width.init("Width (x-axis)", "Width of the detector", default_detector_width, rectDetUnit,
                  3 /* decimals */, RealLimits::positive(), "width");
-    m_height.init("Height (y-axis)", "Height of the detector", default_detector_height, "mm",
+    m_height.init("Height (y-axis)", "Height of the detector", default_detector_height, rectDetUnit,
                   3 /* decimals */, RealLimits::positive(), "height");
 
     m_normalVector.init(
@@ -73,11 +75,11 @@ RectangularDetectorItem::RectangularDetectorItem()
                            Unit::unitless, "directionVector");
     m_directionVector.setY(-1.0);
 
-    m_u0.init("u0", "", default_detector_width / 2., "mm", 3 /* decimals */,
+    m_u0.init("u0", "", default_detector_width / 2., rectDetUnit, 3 /* decimals */,
               RealLimits::limitless(), "u0");
-    m_v0.init("v0", "", 0.0, "mm", 3 /* decimals */, RealLimits::limitless(), "v0");
+    m_v0.init("v0", "", 0.0, rectDetUnit, 3 /* decimals */, RealLimits::limitless(), "v0");
     m_distance.init("Distance", "Distance from the sample origin to the detector plane",
-                    default_detector_distance, "mm", "distance");
+                    default_detector_distance, rectDetUnit, "distance");
 
     m_detectorAlignment = ComboProperty::fromList(
         alignment_names_map.values(),
@@ -167,7 +169,7 @@ void RectangularDetectorItem::readFrom(QXmlStreamReader* r)
         // parameters from base class
         if (tag == Tag::BaseData) {
             DetectorItem::readFrom(r);
-            m_resolutionFunction->setUnit("mm");
+            m_resolutionFunction->setUnit(rectDetUnit);
             XML::gotoEndElementOfTag(r, tag);
 
             // x size
diff --git a/GUI/Model/Detector/SphericalDetectorItem.cpp b/GUI/Model/Detector/SphericalDetectorItem.cpp
index 10137fb9a19677253fcf5a385956c2c065c92dbd..ee24ade809ac3e304fb0b7f71accea095f5a7512 100644
--- a/GUI/Model/Detector/SphericalDetectorItem.cpp
+++ b/GUI/Model/Detector/SphericalDetectorItem.cpp
@@ -26,25 +26,29 @@ const QString BaseData("BaseData");
 
 } // namespace Tag
 
+constexpr Unit sphDetUnit = Unit::degree;
+
 void initResolutionFunction(ResolutionFunctionItem* newFunc, const ResolutionFunctionItem*)
 {
-    newFunc->setUnit(Unit::degree);
+    newFunc->setUnit(sphDetUnit);
 }
 } // namespace
 
 SphericalDetectorItem::SphericalDetectorItem()
+    : m_phiAxis(sphDetUnit)
+    , m_alphaAxis(sphDetUnit)
 {
     m_resolutionFunction.initWithInitializer("Resolution function", "Detector resolution function",
                                              initResolutionFunction);
 
-    m_phiAxis.initMin("Min", "Lower edge of first phi-bin", -1.0, Unit::degree,
+    m_phiAxis.initMin("Min", "Lower edge of first phi-bin", -1.0, m_phiAxis.units(),
                       RealLimits::limited(-90, 90));
-    m_phiAxis.initMax("Max", "Upper edge of last phi-bin", 1.0, Unit::degree,
+    m_phiAxis.initMax("Max", "Upper edge of last phi-bin", 1.0, m_phiAxis.units(),
                       RealLimits::limited(-90, 90));
 
-    m_alphaAxis.initMin("Min", "Lower edge of first alpha-bin", 0.0, Unit::degree,
+    m_alphaAxis.initMin("Min", "Lower edge of first alpha-bin", 0.0, m_alphaAxis.units(),
                         RealLimits::limited(-90, 90));
-    m_alphaAxis.initMax("Max", "Upper edge of last alpha-bin", 2.0, Unit::degree,
+    m_alphaAxis.initMax("Max", "Upper edge of last alpha-bin", 2.0, m_alphaAxis.units(),
                         RealLimits::limited(-90, 90));
 }
 
@@ -79,7 +83,7 @@ void SphericalDetectorItem::readFrom(QXmlStreamReader* r)
         // parameters from base class
         if (tag == Tag::BaseData) {
             DetectorItem::readFrom(r);
-            m_resolutionFunction->setUnit(Unit::degree);
+            m_resolutionFunction->setUnit(sphDetUnit);
             XML::gotoEndElementOfTag(r, tag);
 
             // phi axis
@@ -100,12 +104,16 @@ void SphericalDetectorItem::readFrom(QXmlStreamReader* r)
 std::unique_ptr<IDetector> SphericalDetectorItem::createDomainDetector() const
 {
     const int n_x = m_phiAxis.nbins();
-    const double x_min = Units::deg2rad(m_phiAxis.min());
-    const double x_max = Units::deg2rad(m_phiAxis.max());
+    const double x_min =
+        Units::convert(m_phiAxis.min(), m_phiAxis.units(), Units::internalAngularUnit);
+    const double x_max =
+        Units::convert(m_phiAxis.max(), m_phiAxis.units(), Units::internalAngularUnit);
 
     const int n_y = m_alphaAxis.nbins();
-    const double y_min = Units::deg2rad(m_alphaAxis.min());
-    const double y_max = Units::deg2rad(m_alphaAxis.max());
+    const double y_min =
+        Units::convert(m_alphaAxis.min(), m_alphaAxis.units(), Units::internalAngularUnit);
+    const double y_max =
+        Units::convert(m_alphaAxis.max(), m_alphaAxis.units(), Units::internalAngularUnit);
 
     return std::make_unique<SphericalDetector>(n_x, x_min, x_max, n_y, y_min, y_max);
 }
diff --git a/GUI/Model/Device/InstrumentItems.cpp b/GUI/Model/Device/InstrumentItems.cpp
index 2fd3c7d5f48221653196884497c4a2171236ce06..0a55927b3f0f7569d93b93ff03840581a149c5ad 100644
--- a/GUI/Model/Device/InstrumentItems.cpp
+++ b/GUI/Model/Device/InstrumentItems.cpp
@@ -118,6 +118,16 @@ bool InstrumentItem::alignedWith(const RealItem* item) const
     return shape() == item->shape();
 }
 
+std::vector<Unit> InstrumentItem::availableUnits() const
+{
+    return {Unit::nbins, Unit::radian, Unit::degree, Unit::inv_angstrom, Unit::inv_nanometer};
+}
+
+Unit InstrumentItem::defaultUnits() const
+{
+    return Unit::degree;
+}
+
 void InstrumentItem::writeTo(QXmlStreamWriter* w) const
 {
     XML::writeAttribute(w, XML::Attrib::version, uint(2));
@@ -417,7 +427,7 @@ Frame SpecularInstrumentItem::makeFrame() const
         return Frame(std::move(f));
     }
 
-    Scale axis = axis_item->makeScale("alpha_i (rad)");
+    Scale axis = axis_item->makeScale("alpha_i").convertedScale(Units::internalAngularUnit);
 
     std::vector<const Scale*> f{axis.clone()};
     return Frame(std::move(f));
@@ -482,15 +492,16 @@ void SpecularInstrumentItem::readFrom(QXmlStreamReader* r)
 
 DepthprobeInstrumentItem::DepthprobeInstrumentItem()
     : ScanningFunctionality(1e8)
+    , m_zAxis(Unit::nanometer)
 {
     auto* axisItem = scanItem()->inclinationAxisItem();
     axisItem->setMin(0.0);
     axisItem->setMax(1.0);
     axisItem->setBinCount(500);
 
-    m_zAxis.initMin("Min", "Starting value below sample horizon", -100.0, Unit::nanometer,
+    m_zAxis.initMin("Min", "Starting value below sample horizon", -100.0, m_zAxis.units(),
                     RealLimits::limitless());
-    m_zAxis.initMax("Max", "Ending value above sample horizon", 100.0, Unit::nanometer,
+    m_zAxis.initMax("Max", "Ending value above sample horizon", 100.0, m_zAxis.units(),
                     RealLimits::limitless());
 }
 
@@ -507,12 +518,22 @@ void DepthprobeInstrumentItem::updateToRealData(const RealItem*)
 Frame DepthprobeInstrumentItem::makeFrame() const
 {
     BasicAxisItem* const axis_item = scanItem()->inclinationAxisItem();
-    Scale xAxis = axis_item->makeScale("alpha_i (rad)");
-    Scale zAxis = m_zAxis.createAxis("z (nm)");
+    Scale xAxis = axis_item->makeScale("alpha_i").convertedScale(Units::internalAngularUnit);
+    Scale zAxis = m_zAxis.makeScale("z").convertedScale(Units::internalDepthUnit);
     std::vector<const Scale*> axes{xAxis.clone(), zAxis.clone()};
     return Frame(std::move(axes));
 }
 
+std::vector<Unit> DepthprobeInstrumentItem::availableDepthUnits() const
+{
+    return {Unit::angstrom, Unit::nanometer, Unit::micrometer};
+}
+
+Unit DepthprobeInstrumentItem::defaultDepthUnits() const
+{
+    return Unit::nanometer;
+}
+
 ISimulation* DepthprobeInstrumentItem::createSimulation(const MultiLayer& sample) const
 {
     const Frame frame = makeFrame();
@@ -601,7 +622,7 @@ void OffspecInstrumentItem::updateToRealData(const RealItem* dataItem)
 Frame OffspecInstrumentItem::makeFrame() const
 {
     BasicAxisItem* const axis_item = scanItem()->inclinationAxisItem();
-    Scale xAxis = axis_item->makeScale("alpha_i (rad)");
+    Scale xAxis = axis_item->makeScale("alpha_i").convertedScale(Units::internalAngularUnit);
     return Frame(std::vector<const Scale*>{
         xAxis.clone(), detectorItem()->createOffspecDetector()->axis(1).clone()});
 }
@@ -704,6 +725,21 @@ Frame GISASInstrumentItem::makeFrame() const
     return normalDetector()->scatteringCoords();
 }
 
+std::vector<Unit> GISASInstrumentItem::availableUnits() const
+{
+    auto result = InstrumentItem::availableUnits();
+    if (dynamic_cast<const RectangularDetectorItem*>(m_detector.currentItem()))
+        result.push_back(Unit::millimeter);
+    return result;
+}
+
+Unit GISASInstrumentItem::defaultUnits() const
+{
+    if (dynamic_cast<const RectangularDetectorItem*>(m_detector.currentItem()))
+        return Unit::millimeter;
+    return InstrumentItem::defaultUnits();
+}
+
 ISimulation* GISASInstrumentItem::createSimulation(const MultiLayer& sample) const
 {
     const auto beam = beamItem()->createBeam();
diff --git a/GUI/Model/Device/InstrumentItems.h b/GUI/Model/Device/InstrumentItems.h
index 2ad1fa776ee0ae6661b7968d44c845369c5787de..1eb002d782a258ca20af713b76dbaadf60cf9b39 100644
--- a/GUI/Model/Device/InstrumentItems.h
+++ b/GUI/Model/Device/InstrumentItems.h
@@ -59,6 +59,11 @@ public:
 
     virtual Frame makeFrame() const = 0;
 
+    virtual std::vector<Unit> availableUnits() const;
+
+    //! These units are used to show the simulation result
+    virtual Unit defaultUnits() const;
+
     virtual ISimulation* createSimulation(const MultiLayer& sample) const = 0;
 
     virtual void writeTo(QXmlStreamWriter* w) const;
@@ -173,6 +178,8 @@ public:
     std::vector<int> shape() const override;
     void updateToRealData(const RealItem* item) override;
     Frame makeFrame() const override;
+    std::vector<Unit> availableDepthUnits() const;
+    Unit defaultDepthUnits() const;
     ISimulation* createSimulation(const MultiLayer& sample) const override;
     void writeTo(QXmlStreamWriter* w) const override;
     void readFrom(QXmlStreamReader* r) override;
@@ -211,6 +218,8 @@ public:
     std::vector<int> shape() const override;
     void updateToRealData(const RealItem* item) override;
     Frame makeFrame() const override;
+    std::vector<Unit> availableUnits() const override;
+    Unit defaultUnits() const override;
     ISimulation* createSimulation(const MultiLayer& sample) const override;
     void writeTo(QXmlStreamWriter* w) const override;
     void readFrom(QXmlStreamReader* r) override;
diff --git a/GUI/Model/Job/JobItem.cpp b/GUI/Model/Job/JobItem.cpp
index 86319083efa68d063adb92b2e03c0edc232cf92c..53f79a2c5d4b2d880c4141ffcd591c80dccb4147 100644
--- a/GUI/Model/Job/JobItem.cpp
+++ b/GUI/Model/Job/JobItem.cpp
@@ -17,6 +17,7 @@
 #include "Device/Data/Datafield.h"
 #include "Device/Detector/IDetector.h"
 #include "Device/Detector/SimulationAreaIterator.h" // roiIndex
+#include "GUI/Model/Axis/AmplitudeAxisItem.h"
 #include "GUI/Model/Data/IntensityDataItem.h"
 #include "GUI/Model/Data/SpecularDataItem.h"
 #include "GUI/Model/Device/BackgroundItems.h"
@@ -293,6 +294,23 @@ void JobItem::createSimulatedDataItem()
 {
     ASSERT(!simulatedDataItem());
     m_simulatedDataItem.reset(createNewDataItem());
+
+    // Set default axes units for simulated data.
+    // Can be overriden by units from RealItem
+    Unit default_units = instrumentItem()->defaultUnits();
+    if (instrumentItem()->is<SpecularInstrumentItem>()) {
+        m_simulatedDataItem->xAxisItem()->setUnits(default_units);
+    } else if (instrumentItem()->is<GISASInstrumentItem>()) {
+        m_simulatedDataItem->xAxisItem()->setUnits(default_units);
+        m_simulatedDataItem->yAxisItem()->setUnits(default_units);
+    } else if (instrumentItem()->is<OffspecInstrumentItem>()) {
+        m_simulatedDataItem->xAxisItem()->setUnits(default_units);
+        m_simulatedDataItem->yAxisItem()->setUnits(default_units);
+    } else if (const auto dProbe = dynamic_cast<DepthprobeInstrumentItem*>(instrumentItem())) {
+        m_simulatedDataItem->xAxisItem()->setUnits(default_units);
+        m_simulatedDataItem->yAxisItem()->setUnits(dProbe->defaultDepthUnits());
+    } else
+        ASSERT(false);
 }
 
 IntensityDataItem* JobItem::intensityDataItem()
@@ -312,7 +330,7 @@ DataItem* JobItem::createDiffDataItem()
 
     // use the same axes units as for real data
     ASSERT(m_realItem);
-    // TODO
+    m_diffDataItem->copyUnitsFromItem(realItem()->dataItem());
 
     if (isSpecularJob())
         dynamic_cast<SpecularDataItem*>(diffDataItem())->setDiffPlotStyle();
@@ -339,7 +357,7 @@ void JobItem::copyRealItemIntoJob(const RealItem* srcRealItem)
 
     // override axes units of simulated data
     ASSERT(m_simulatedDataItem);
-    // m_simulatedDataItem->setCurrentCoord(m_realItem->dataItem()->currentCoord());
+    m_simulatedDataItem->copyUnitsFromItem(realItem()->dataItem());
 
     if (isSpecularJob())
         m_realItem->specularDataItem()->setRealPlotStyle();
diff --git a/Sim/Export/PyFmt2.cpp b/Sim/Export/PyFmt2.cpp
index c9a6c9566240aedd434af4cd44c333ed3a23c5aa..a6ca8025c2f99d99202840f5e3b335d331f522fd 100644
--- a/Sim/Export/PyFmt2.cpp
+++ b/Sim/Export/PyFmt2.cpp
@@ -93,11 +93,13 @@ std::string Py::Fmt2::printAxis(const Scale* a, const std::string& unit)
 {
     std::ostringstream result;
     if (a->isEquiDivision())
-        result << "ba.EquiDivision(" << Py::Fmt::printString(a->axisName()) << ", " << a->size()
-               << ", " << Py::Fmt::printValue(a->min(), unit) << ", "
-               << Py::Fmt::printValue(a->max(), unit) << ")";
+        result << "ba.EquiDivision(" << Py::Fmt::printString(a->axisName())
+               << Py::Fmt::printUnit(a->axisUnits()) << ", " << a->size() << ", "
+               << Py::Fmt::printValue(a->min(), unit) << ", " << Py::Fmt::printValue(a->max(), unit)
+               << ")";
     else if (a->isEquiScan())
-        result << "ba.EquiScan(" << Py::Fmt::printString(a->axisName()) << ", " << a->size() << ", "
+        result << "ba.EquiScan(" << Py::Fmt::printString(a->axisName()) << ", "
+               << Py::Fmt::printUnit(a->axisUnits()) << ", " << a->size() << ", "
                << Py::Fmt::printValue(a->min(), unit) << ", " << Py::Fmt::printValue(a->max(), unit)
                << ")";
     else if (a->isScan()) {
diff --git a/Sim/Scan/AlphaScan.cpp b/Sim/Scan/AlphaScan.cpp
index 293740ba0bd7d4370ec9679ae68b5fe4ffb80a7a..0273e7ffd46ec31b35353c08dea7d5f45d5b0cae 100644
--- a/Sim/Scan/AlphaScan.cpp
+++ b/Sim/Scan/AlphaScan.cpp
@@ -50,7 +50,7 @@ AlphaScan::AlphaScan(const Scale& alpha_axis)
 }
 
 AlphaScan::AlphaScan(int nbins, double alpha_i_min, double alpha_i_max)
-    : AlphaScan(EquiScan("alpha_i (rad)", nbins, alpha_i_min, alpha_i_max))
+    : AlphaScan(EquiScan("alpha_i", Unit::radian, nbins, alpha_i_min, alpha_i_max))
 {
 }
 
diff --git a/Sim/Scan/LambdaScan.cpp b/Sim/Scan/LambdaScan.cpp
index 145cde63ba092ff8c423ac31061b2556c3f24d78..e0cf36ece613837f2b886bd5a6de52fcea90890c 100644
--- a/Sim/Scan/LambdaScan.cpp
+++ b/Sim/Scan/LambdaScan.cpp
@@ -42,7 +42,7 @@ LambdaScan::LambdaScan(double alpha_i, Scale* lambdaScale)
 }
 
 LambdaScan::LambdaScan(double alpha_i, std::vector<double> lambdaList)
-    : LambdaScan(alpha_i, newListScan("qs", std::move(lambdaList)))
+    : LambdaScan(alpha_i, newListScan("lambda", Unit::nanometer, std::move(lambdaList)))
 {
 }
 
@@ -52,7 +52,7 @@ LambdaScan::LambdaScan(double alpha_i, const Scale& lambdaScale)
 }
 
 LambdaScan::LambdaScan(double alpha_i, int nbins, double lambda_min, double lambda_max)
-    : LambdaScan(alpha_i, newEquiScan("lambda", nbins, lambda_min, lambda_max))
+    : LambdaScan(alpha_i, newEquiScan("lambda", Unit::nanometer, nbins, lambda_min, lambda_max))
 {
 }
 
diff --git a/Sim/Scan/QzScan.cpp b/Sim/Scan/QzScan.cpp
index 2d2d4efb78c7a4ca79ab55f252a805429d6c6ec2..d596b52e7dea7845b777219100858b58728bfcb8 100644
--- a/Sim/Scan/QzScan.cpp
+++ b/Sim/Scan/QzScan.cpp
@@ -37,7 +37,7 @@ QzScan::QzScan(Scale* qs_nm)
 }
 
 QzScan::QzScan(std::vector<double> qs_nm)
-    : QzScan(newListScan("qs", std::move(qs_nm)))
+    : QzScan(newListScan("qs", Unit::inv_nanometer, std::move(qs_nm)))
 {
 }
 
@@ -47,7 +47,7 @@ QzScan::QzScan(const Scale& qs_nm)
 }
 
 QzScan::QzScan(int nbins, double qz_min, double qz_max)
-    : QzScan(newEquiScan("qs", nbins, qz_min, qz_max))
+    : QzScan(newEquiScan("qs", Unit::inv_nanometer, nbins, qz_min, qz_max))
 {
 }
 
diff --git a/Tests/SimFactory/MakeSimulations.cpp b/Tests/SimFactory/MakeSimulations.cpp
index 3a6e153f9dbf44a7b603e3511757922085beb476..5cac23c50e0bbb4ddfa85d375371b38db2446540 100644
--- a/Tests/SimFactory/MakeSimulations.cpp
+++ b/Tests/SimFactory/MakeSimulations.cpp
@@ -326,14 +326,14 @@ IBeamScan* test::makeSimulation::BasicSpecularScan(bool vsQ)
     const double max_angle = 5 * deg;
 
     if (vsQ) {
-        Scale angle_axis = EquiScan("qz (1/nm)", n, max_angle / n, max_angle);
+        Scale angle_axis = EquiScan("qz", Unit::radian, n, max_angle / n, max_angle);
         auto angles = angle_axis.binCenters();
         std::vector<double> qs(angle_axis.size(), 0.0);
         for (size_t i = 0, size = qs.size(); i < size; ++i)
             qs[i] = 4.0 * pi * std::sin(angles[i]) / wavelength;
         return new QzScan(qs);
     }
-    auto* result = new AlphaScan(EquiScan("alpha_i (rad)", n, max_angle / n, max_angle));
+    auto* result = new AlphaScan(EquiScan("alpha_i", Unit::radian, n, max_angle / n, max_angle));
     result->setWavelength(wavelength);
     return result;
 }
@@ -367,7 +367,7 @@ test::makeSimulation::SpecularWithGaussianBeam(const MultiLayer& sample)
     const int n = 2000;
     const double max_angle = 5 * deg;
     auto gaussian_ff = std::make_unique<FootprintGauss>(1.0);
-    AlphaScan scan(EquiScan("alpha_i (rad)", n, max_angle / n, max_angle));
+    AlphaScan scan(EquiScan("alpha_i", Unit::radian, n, max_angle / n, max_angle));
     scan.setWavelength(wavelength);
     scan.setFootprint(gaussian_ff.get());
 
@@ -381,7 +381,7 @@ test::makeSimulation::SpecularWithSquareBeam(const MultiLayer& sample)
     const int n = 2000;
     const double max_angle = 5 * deg;
     auto square_ff = std::make_unique<FootprintSquare>(1.0);
-    AlphaScan scan(EquiScan("alpha_i (rad)", n, max_angle / n, max_angle));
+    AlphaScan scan(EquiScan("alpha_i", Unit::radian, n, max_angle / n, max_angle));
     scan.setWavelength(wavelength);
     scan.setFootprint(square_ff.get());
 
@@ -397,7 +397,7 @@ test::makeSimulation::SpecularDivergentBeam(const MultiLayer& sample)
     const double max_angle = 5 * deg;
     const double wavelength_stddev = 0.1 * angstrom;
     const double alpha_stddev = 0.1 * deg;
-    AlphaScan scan(EquiScan("alpha_i (rad)", n, max_angle / n, max_angle));
+    AlphaScan scan(EquiScan("alpha_i", Unit::radian, n, max_angle / n, max_angle));
 
     DistributionGaussian wavelength_distr(wavelength, wavelength_stddev, n_integration_points, 2.);
     DistributionGaussian alpha_distr(0, alpha_stddev, n_integration_points, 2.);
@@ -411,7 +411,7 @@ test::makeSimulation::SpecularDivergentBeam(const MultiLayer& sample)
 std::unique_ptr<SpecularSimulation>
 test::makeSimulation::TOFRWithRelativeResolution(const MultiLayer& sample)
 {
-    Scale qs = EquiScan("alpha_i (rad)", 500, 0.0, 1.0);
+    Scale qs = EquiScan("alpha_i", Unit::radian, 500, 0.0, 1.0);
     QzScan scan(qs);
     scan.setRelativeQResolution(DistributionGaussian(0., 1., 20, 2.0), 0.03);
 
@@ -423,7 +423,7 @@ test::makeSimulation::TOFRWithRelativeResolution(const MultiLayer& sample)
 std::unique_ptr<SpecularSimulation>
 test::makeSimulation::TOFRWithPointwiseResolution(const MultiLayer& sample)
 {
-    Scale qs = EquiScan("alpha_i (rad)", 500, 0.0, 1.0);
+    Scale qs = EquiScan("alpha_i", Unit::radian, 500, 0.0, 1.0);
     QzScan scan(qs);
 
     std::vector<double> resolutions;
@@ -448,7 +448,7 @@ test::makeSimulation::BasicDepthprobe(const MultiLayer& sample)
     AlphaScan scan(20, 0.025 * deg, 0.975 * deg);
     scan.setWavelength(1.);
 
-    Scale zaxis = EquiDivision("z (nm)", 20, -100., +100.);
+    Scale zaxis = EquiDivision("z", Unit::nanometer, 20, -100., +100.);
 
     return std::make_unique<DepthprobeSimulation>(scan, sample, zaxis);
 }
diff --git a/Tests/Unit/Device/ArrayUtilsTest.cpp b/Tests/Unit/Device/ArrayUtilsTest.cpp
index 11b78e91118762629bd7883ff119f86126f1c440..f3eb6f2a83670e1afce724d8a208b3106c3384a3 100644
--- a/Tests/Unit/Device/ArrayUtilsTest.cpp
+++ b/Tests/Unit/Device/ArrayUtilsTest.cpp
@@ -20,7 +20,7 @@ TEST(ArrayUtilsTest, DatafieldFromVector1D)
 TEST(ArrayUtilsTest, DatafieldToVector1D)
 {
     const std::vector<double> expected = {10.0, 20.0, 30.0, 40.0};
-    Datafield data({newEquiDivision("axis0", 4, 10.0, 20.0)});
+    Datafield data({newEquiDivision("axis0", Unit::other, 4, 10.0, 20.0)});
     data.setVector(expected);
 
     auto vec = DataUtil::Array::createVector1D(data);
@@ -52,8 +52,8 @@ TEST(ArrayUtilsTest, DatafieldFromVector2D)
 
 TEST(ArrayUtilsTest, DatafieldToVector2D)
 {
-    Datafield data(
-        {newEquiDivision("axis0", 4, 10.0, 20.0), newEquiDivision("axis1", 3, 30.0, 40.0)});
+    Datafield data({newEquiDivision("axis0", Unit::other, 4, 10.0, 20.0),
+                    newEquiDivision("axis1", Unit::other, 3, 30.0, 40.0)});
     const std::vector<double> values = {8.0,  4.0, 0.0, 9.0,  5.0, 1.0,
                                         10.0, 6.0, 2.0, 11.0, 7.0, 3.0};
 
diff --git a/Tests/Unit/Device/FixedBinAxisTest.cpp b/Tests/Unit/Device/FixedBinAxisTest.cpp
index 663e841869313a17287c9826b0356da5b97284ca..24edf4a7979294b5e2709f9bdc3f49d2f5833352 100644
--- a/Tests/Unit/Device/FixedBinAxisTest.cpp
+++ b/Tests/Unit/Device/FixedBinAxisTest.cpp
@@ -6,7 +6,7 @@
 
 TEST(EquiDivisionTest, IndexedAccessor)
 {
-    Scale a1 = EquiDivision("length", 100, 0.0, 10.0);
+    Scale a1 = EquiDivision("length", Unit::other, 100, 0.0, 10.0);
     EXPECT_EQ(100u, a1.size());
     EXPECT_EQ(0.0, a1.min());
     EXPECT_EQ(10.0, a1.max());
@@ -15,7 +15,7 @@ TEST(EquiDivisionTest, IndexedAccessor)
     EXPECT_DOUBLE_EQ(6.55, a1.binCenter(65));
     EXPECT_DOUBLE_EQ(9.95, a1.binCenter(99));
 
-    Scale a2 = EquiDivision("name", 3, -1.5, 1.5);
+    Scale a2 = EquiDivision("name", Unit::other, 3, -1.5, 1.5);
     EXPECT_DOUBLE_EQ(-1.0, a2.binCenter(0));
     EXPECT_DOUBLE_EQ(0.0, a2.binCenter(1));
     EXPECT_DOUBLE_EQ(1.0, a2.binCenter(2));
@@ -24,7 +24,7 @@ TEST(EquiDivisionTest, IndexedAccessor)
 
 TEST(EquiDivisionTest, VectorOfUnitLength)
 {
-    Scale axis = EquiDivision("name", 1, 1.0, 2.0);
+    Scale axis = EquiDivision("name", Unit::other, 1, 1.0, 2.0);
     EXPECT_EQ(1u, axis.size());
     EXPECT_EQ(double(1.0), axis.min());
     EXPECT_EQ(double(2.0), axis.max());
@@ -33,7 +33,7 @@ TEST(EquiDivisionTest, VectorOfUnitLength)
 
 TEST(EquiDivisionTest, FindClosestIndex)
 {
-    Scale v1 = EquiDivision("name", 2, 0.0, 1.0);
+    Scale v1 = EquiDivision("name", Unit::other, 2, 0.0, 1.0);
     EXPECT_EQ(size_t(2), v1.size());
     EXPECT_EQ(size_t(0), v1.closestIndex(0.0));
     EXPECT_EQ(size_t(0), v1.closestIndex(0.25));
@@ -42,7 +42,7 @@ TEST(EquiDivisionTest, FindClosestIndex)
     //    EXPECT_THROW( v1.closestIndex(1.0), std::runtime_error);
     EXPECT_EQ(size_t(1), v1.closestIndex(1.0));
 
-    Scale v2 = EquiDivision("name", 3, -1.5, 1.5);
+    Scale v2 = EquiDivision("name", Unit::other, 3, -1.5, 1.5);
     EXPECT_EQ(size_t(0), v2.closestIndex(-1.5));
     EXPECT_EQ(size_t(0), v2.closestIndex(-1.0));
     EXPECT_EQ(size_t(1), v2.closestIndex(-0.5));
@@ -55,7 +55,7 @@ TEST(EquiDivisionTest, FindClosestIndex)
 
 TEST(EquiDivisionTest, CheckBin)
 {
-    Scale axis = EquiDivision("name", 20, 0, 10);
+    Scale axis = EquiDivision("name", Unit::other, 20, 0, 10);
 
     Bin1D bin0 = axis.bin(0);
     EXPECT_DOUBLE_EQ(0.25, bin0.center());
@@ -81,7 +81,7 @@ TEST(EquiDivisionTest, CheckBin)
 
     EXPECT_THROW(axis.bin(20), std::out_of_range);
 
-    Scale axis2 = EquiDivision("name", 3, -1, 2.0);
+    Scale axis2 = EquiDivision("name", Unit::other, 3, -1, 2.0);
     EXPECT_DOUBLE_EQ(-0.5, axis2.bin(0).center());
     EXPECT_DOUBLE_EQ(0.5, axis2.bin(1).center());
     EXPECT_DOUBLE_EQ(1.5, axis2.bin(2).center());
@@ -89,22 +89,24 @@ TEST(EquiDivisionTest, CheckBin)
 
 TEST(EquiDivisionTest, CheckEquality)
 {
-    Scale b1 = EquiDivision("axis", 99, -1.01, 3.3);
-    Scale b2 = EquiDivision("axis", 99, -1.01, 3.3);
+    Scale b1 = EquiDivision("axis", Unit::other, 99, -1.01, 3.3);
+    Scale b2 = EquiDivision("axis", Unit::other, 99, -1.01, 3.3);
     EXPECT_TRUE(b1 == b2);
-    Scale b3 = EquiDivision("axissss", 99, -1.01, 3.3);
-    Scale b4 = EquiDivision("axis", 99, -1.0, 3.3);
-    Scale b5 = EquiDivision("axis", 99, -1.01, 3.29);
-    Scale b6 = EquiDivision("axiss", 98, -1.01, 3.3);
+    Scale b3 = EquiDivision("axissss", Unit::other, 99, -1.01, 3.3);
+    Scale b4 = EquiDivision("axis", Unit::other, 99, -1.0, 3.3);
+    Scale b5 = EquiDivision("axis", Unit::other, 99, -1.01, 3.29);
+    Scale b6 = EquiDivision("axiss", Unit::other, 98, -1.01, 3.3);
+    Scale b7 = EquiDivision("axis", Unit::unitless, 99, -1.01, 3.3);
     EXPECT_FALSE(b1 == b3);
     EXPECT_FALSE(b1 == b4);
     EXPECT_FALSE(b1 == b5);
     EXPECT_FALSE(b1 == b6);
+    EXPECT_FALSE(b1 == b7);
 }
 
 TEST(EquiDivisionTest, CheckClone)
 {
-    Scale a1 = EquiDivision("axis", 99, -1.2, 5.4);
+    Scale a1 = EquiDivision("axis", Unit::other, 99, -1.2, 5.4);
     Scale* clone = a1.clone();
     EXPECT_TRUE(a1 == *clone);
     delete clone;
@@ -112,7 +114,7 @@ TEST(EquiDivisionTest, CheckClone)
 
 TEST(EquiDivisionTest, IOStream)
 {
-    Scale axis = EquiDivision("name", 99, -1.2, 5.4);
+    Scale axis = EquiDivision("name", Unit::other, 99, -1.2, 5.4);
 
     std::ostringstream oss;
     oss << axis;
@@ -124,7 +126,7 @@ TEST(EquiDivisionTest, IOStream)
 
 TEST(EquiDivisionTest, BinCenters)
 {
-    Scale axis = EquiDivision("name", 3, -1.5, 1.5);
+    Scale axis = EquiDivision("name", Unit::other, 3, -1.5, 1.5);
     std::vector<double> centers = axis.binCenters();
     EXPECT_EQ(size_t(3), centers.size());
     EXPECT_DOUBLE_EQ(-1.0, centers[0]);
diff --git a/Tests/Unit/Device/IOReaderWriterTest.cpp b/Tests/Unit/Device/IOReaderWriterTest.cpp
index df9899078cfc12b8ec67d0ed007cddd0de4db663..34e0bb885ef4212e7e95a450a60b177e542c7254 100644
--- a/Tests/Unit/Device/IOReaderWriterTest.cpp
+++ b/Tests/Unit/Device/IOReaderWriterTest.cpp
@@ -14,7 +14,8 @@ protected:
             m_model_data[i] = static_cast<double>(i);
     }
 
-    Datafield m_model_data{{newEquiDivision("x", 5, 1.0, 5.0), newEquiDivision("y", 10, 6.0, 7.0)}};
+    Datafield m_model_data{{newEquiDivision("x", Unit::other, 5, 1.0, 5.0),
+                            newEquiDivision("y", Unit::other, 10, 6.0, 7.0)}};
 };
 
 
diff --git a/Tests/Unit/Device/IntensityDataFunctionsTest.cpp b/Tests/Unit/Device/IntensityDataFunctionsTest.cpp
index b2cecd6cdf6b6c13285b1b1221c036b41d1a7b9d..014811d0681c3734ad16b89602f4b1a082e16594 100644
--- a/Tests/Unit/Device/IntensityDataFunctionsTest.cpp
+++ b/Tests/Unit/Device/IntensityDataFunctionsTest.cpp
@@ -8,8 +8,8 @@
 
 TEST(DataUtilTest, createRearrangedDataSet)
 {
-    Datafield input_data{
-        {newEquiDivision("axis0", 2, 1.0, 2.0), newEquiDivision("axis1", 3, 3.0, 4.0)}};
+    Datafield input_data{{newEquiDivision("axis0", Unit::other, 2, 1.0, 2.0),
+                          newEquiDivision("axis1", Unit::other, 3, 3.0, 4.0)}};
     input_data.setVector(std::vector<double>{1.0, 2.0, 3.0, 4.0, 5.0, 6.0});
 
     std::unique_ptr<Datafield> output_data;
@@ -55,7 +55,7 @@ TEST(DataUtilTest, createRearrangedDataSet)
 
 TEST(DataUtilTest, coordinateToFromBinf)
 {
-    Scale axis = EquiDivision("axis", 8, -5.0, 3.0);
+    Scale axis = EquiDivision("axis", Unit::other, 8, -5.0, 3.0);
     EXPECT_EQ(0.5, FrameUtil::coordinateToBinf(-4.5, axis));
     EXPECT_EQ(-4.5, FrameUtil::coordinateFromBinf(0.5, axis));
 
@@ -77,8 +77,10 @@ TEST(DataUtilTest, coordinateToFromBinf)
 
 TEST(DataUtilTest, datafieldToFromBinf)
 {
-    Frame data1(newEquiDivision("axis0", 8, -5.0, 3.0), newEquiDivision("axis1", 3, 2.0, 5.0));
-    Frame data2(newEquiDivision("axis0", 8, -10.0, 70.0), newEquiDivision("axis1", 3, -10.0, 20.0));
+    Frame data1(newEquiDivision("axis0", Unit::other, 8, -5.0, 3.0),
+                newEquiDivision("axis1", Unit::other, 3, 2.0, 5.0));
+    Frame data2(newEquiDivision("axis0", Unit::other, 8, -10.0, 70.0),
+                newEquiDivision("axis1", Unit::other, 3, -10.0, 20.0));
 
     double x(-4.5), y(2.5);
     FrameUtil::coordinatesToBinf(x, y, data1);
@@ -96,8 +98,8 @@ TEST(DataUtilTest, datafieldToFromBinf)
 
 TEST(DataUtilTest, create2DArrayfromDatafield)
 {
-    Datafield out_data{
-        {newEquiDivision("axis0", 2, 1.0, 2.0), newEquiDivision("axis1", 3, 3.0, 4.0)}};
+    Datafield out_data{{newEquiDivision("axis0", Unit::other, 2, 1.0, 2.0),
+                        newEquiDivision("axis1", Unit::other, 3, 3.0, 4.0)}};
     EXPECT_EQ(6u, out_data.size());
 
     EXPECT_EQ(2u, out_data.axis(0).size()); // no. of rows
diff --git a/Tests/Unit/Device/PointwiseAxisTest.cpp b/Tests/Unit/Device/PointwiseAxisTest.cpp
index d8a2bca2b1b4c7c58f989824f814f03a368ce853..3c8ec2945ec26bd57516ede797f87a586d46cafe 100644
--- a/Tests/Unit/Device/PointwiseAxisTest.cpp
+++ b/Tests/Unit/Device/PointwiseAxisTest.cpp
@@ -6,19 +6,22 @@
 
 TEST(PointwiseAxisTest, Construction)
 {
-    EXPECT_THROW(ListScan("length", std::vector<double>{1.0, 0.0}), std::runtime_error);
-    EXPECT_THROW(ListScan("length", std::vector<double>{0.0, 1.0, 0.5}), std::runtime_error);
-    EXPECT_THROW(ListScan("length", std::vector<double>{0.0, 1.0, 1.0}), std::runtime_error);
-    Scale a1 = ListScan("length", std::vector<double>{0.0, 1.0});
+    EXPECT_THROW(ListScan("length", Unit::other, std::vector<double>{1.0, 0.0}),
+                 std::runtime_error);
+    EXPECT_THROW(ListScan("length", Unit::other, std::vector<double>{0.0, 1.0, 0.5}),
+                 std::runtime_error);
+    EXPECT_THROW(ListScan("length", Unit::other, std::vector<double>{0.0, 1.0, 1.0}),
+                 std::runtime_error);
+    Scale a1 = ListScan("length", Unit::other, std::vector<double>{0.0, 1.0});
     std::vector<double> coordinates{0.0, 1.0};
-    Scale a2 = ListScan("length", coordinates);
+    Scale a2 = ListScan("length", Unit::other, coordinates);
     EXPECT_TRUE(a1 == a2);
 }
 
 TEST(PointwiseAxisTest, BasicProperties)
 {
     std::vector<double> coordinates{0.0, 1.0, 4.0, 8.0};
-    Scale a1 = ListScan("length", coordinates);
+    Scale a1 = ListScan("length", Unit::other, coordinates);
     EXPECT_EQ(4u, a1.size());
     EXPECT_EQ(0.0, a1.min());
     EXPECT_EQ(8.0, a1.max());
@@ -37,7 +40,7 @@ TEST(PointwiseAxisTest, BasicProperties)
 
 TEST(PointwiseAxisTest, FindClosestIndex)
 {
-    Scale v1 = ListScan("name", std::vector<double>{0.0, 1.0, 4.0, 8.0});
+    Scale v1 = ListScan("name", Unit::other, std::vector<double>{0.0, 1.0, 4.0, 8.0});
     EXPECT_EQ(4u, v1.size());
     EXPECT_EQ(v1.closestIndex(-1.0), 0u);
     EXPECT_EQ(v1.closestIndex(0.0), 0u);
@@ -49,7 +52,7 @@ TEST(PointwiseAxisTest, FindClosestIndex)
     EXPECT_EQ(3u, v1.closestIndex(8.0));
     EXPECT_EQ(3u, v1.closestIndex(11.0));
 
-    Scale v2 = ListScan("name", std::vector<double>{-2.0, -1.0});
+    Scale v2 = ListScan("name", Unit::other, std::vector<double>{-2.0, -1.0});
     EXPECT_EQ(2u, v2.size());
     EXPECT_EQ(v2.closestIndex(-3.0), 0u);
     EXPECT_EQ(v2.closestIndex(-2.0), 0u);
@@ -61,12 +64,12 @@ TEST(PointwiseAxisTest, FindClosestIndex)
 
 TEST(PointwiseAxisTest, CheckEquality)
 {
-    Scale b1 = ListScan("axis", std::vector<double>{1.0, 2.0, 5.0});
-    Scale b2 = ListScan("axis", std::vector<double>{1.0, 2.0, 5.0});
+    Scale b1 = ListScan("axis", Unit::other, std::vector<double>{1.0, 2.0, 5.0});
+    Scale b2 = ListScan("axis", Unit::other, std::vector<double>{1.0, 2.0, 5.0});
     EXPECT_TRUE(b1 == b2);
-    Scale b3 = ListScan("axissss", std::vector<double>{1.0, 2.0, 5.0});
-    Scale b4 = ListScan("axis", std::vector<double>{1.0, 2.0, 6.0});
-    Scale b6 = ListScan("axiss", std::vector<double>{1.5, 2.0, 5.0});
+    Scale b3 = ListScan("axissss", Unit::other, std::vector<double>{1.0, 2.0, 5.0});
+    Scale b4 = ListScan("axis", Unit::other, std::vector<double>{1.0, 2.0, 6.0});
+    Scale b6 = ListScan("axiss", Unit::other, std::vector<double>{1.5, 2.0, 5.0});
     EXPECT_FALSE(b1 == b3);
     EXPECT_FALSE(b1 == b4);
     EXPECT_FALSE(b1 == b6);
@@ -74,14 +77,14 @@ TEST(PointwiseAxisTest, CheckEquality)
 
 TEST(PointwiseAxisTest, CheckClone)
 {
-    Scale a1 = ListScan("axis", std::vector<double>{1.0, 2.0, 5.0});
+    Scale a1 = ListScan("axis", Unit::other, std::vector<double>{1.0, 2.0, 5.0});
     Scale* clone = a1.clone();
     EXPECT_TRUE(a1 == *clone);
 }
 
 TEST(PointwiseAxisTest, IOStream)
 {
-    Scale axis = ListScan("name", std::vector<double>{1.0, 2.0, 5.0});
+    Scale axis = ListScan("name", Unit::other, std::vector<double>{1.0, 2.0, 5.0});
 
     std::ostringstream oss;
     oss << axis;
@@ -93,7 +96,7 @@ TEST(PointwiseAxisTest, IOStream)
 
 TEST(PointwiseAxisTest, ClippedAxis)
 {
-    Scale axis = ListScan("name", std::vector<double>{1.0, 2.0, 2.5, 2.7, 5.0});
+    Scale axis = ListScan("name", Unit::other, std::vector<double>{1.0, 2.0, 2.5, 2.7, 5.0});
 
     Scale clip1 = axis.clipped(0.99, 5.1);
     EXPECT_TRUE(clip1 == axis);
@@ -113,13 +116,13 @@ TEST(PointwiseAxisTest, ClippedAxis)
 
 TEST(PointwiseAxisTest, EquiDivisionComparisonWithMask)
 {
-    Scale axis = EquiDivision("reference", 10, 0.0, 10.0);
+    Scale axis = EquiDivision("reference", Unit::other, 10, 0.0, 10.0);
 
     const std::vector<size_t> mask{0u, 2u, 3u, 4u, 7u, 8u, 9u};
     std::vector<double> coordinates;
     for (auto index : mask)
         coordinates.push_back(axis.binCenter(index));
-    Scale pointwise_axis = ListScan("pointwise", coordinates);
+    Scale pointwise_axis = ListScan("pointwise", Unit::other, coordinates);
 
     // comparing on-axis values
     for (size_t i = 0; i < mask.size(); ++i)
diff --git a/Tests/Unit/Device/PowerfieldTest.cpp b/Tests/Unit/Device/PowerfieldTest.cpp
index 94dc1ef78359fd3183b8611cbe875f2f36eeddca..ba9210572dfd5dbe2e252dbaac1cc85abde0737b 100644
--- a/Tests/Unit/Device/PowerfieldTest.cpp
+++ b/Tests/Unit/Device/PowerfieldTest.cpp
@@ -8,7 +8,8 @@
 class DatafieldTest : public ::testing::Test {
 protected:
     DatafieldTest()
-        : db_data{{newEquiDivision("angle", 20, 0.0, 0.1), newEquiDivision("length", 10, 0.0, 0.5)}}
+        : db_data{{newEquiDivision("angle", Unit::other, 20, 0.0, 0.1),
+                   newEquiDivision("length", Unit::other, 10, 0.0, 0.5)}}
     {
         for (size_t i = 0; i < 200; ++i)
             db_data[i] = (double)i;
diff --git a/Tests/Unit/Device/RegionOfInterestTest.cpp b/Tests/Unit/Device/RegionOfInterestTest.cpp
index acc5ef9aad1fa0cf028955975a89346cad9806b7..1f6ae8d3d039caba3fa351a2540d6dc98cd770d9 100644
--- a/Tests/Unit/Device/RegionOfInterestTest.cpp
+++ b/Tests/Unit/Device/RegionOfInterestTest.cpp
@@ -8,7 +8,8 @@
 TEST(RegionOfInterestTest, constructor)
 {
     SphericalDetector detector(std::array<std::shared_ptr<Scale>, 2>{
-        sharedEquiDivision("axis0", 8, -3.0, 5.0), sharedEquiDivision("axis1", 4, 0.0, 4.0)});
+        sharedEquiDivision("axis0", Unit::other, 8, -3.0, 5.0),
+        sharedEquiDivision("axis1", Unit::other, 4, 0.0, 4.0)});
 
     // creating region of interest
     double xlow(-1.9), ylow(1.1), xup(2.9), yup(2.85);
@@ -34,7 +35,8 @@ TEST(RegionOfInterestTest, constructor)
 TEST(RegionOfInterestTest, largeArea)
 {
     SphericalDetector detector(std::array<std::shared_ptr<Scale>, 2>{
-        sharedEquiDivision("axis0", 8, -3.0, 5.0), sharedEquiDivision("axis1", 4, 0.0, 4.0)});
+        sharedEquiDivision("axis0", Unit::other, 8, -3.0, 5.0),
+        sharedEquiDivision("axis1", Unit::other, 4, 0.0, 4.0)});
 
     // creating region of interest
     double xlow(-3.9), ylow(-1.1), xup(6.9), yup(5.85);
@@ -56,7 +58,8 @@ TEST(RegionOfInterestTest, largeArea)
 TEST(RegionOfInterestTest, clone)
 {
     SphericalDetector detector(std::array<std::shared_ptr<Scale>, 2>{
-        sharedEquiDivision("axis0", 8, -3.0, 5.0), sharedEquiDivision("axis1", 4, 0.0, 4.0)});
+        sharedEquiDivision("axis0", Unit::other, 8, -3.0, 5.0),
+        sharedEquiDivision("axis1", Unit::other, 4, 0.0, 4.0)});
 
     // creating region of interest
     double xlow(-1.9), ylow(1.1), xup(2.9), yup(2.85);
diff --git a/Tests/Unit/Device/SphericalDetectorTest.cpp b/Tests/Unit/Device/SphericalDetectorTest.cpp
index 9a005539ffb263536e80223eb4968a3d4c8c0b77..25657e29bd5b7aeffee688a35ba206072ff00fab 100644
--- a/Tests/Unit/Device/SphericalDetectorTest.cpp
+++ b/Tests/Unit/Device/SphericalDetectorTest.cpp
@@ -15,8 +15,8 @@
 // Construction of the detector with axes.
 TEST(SphericalDetectorTest, constructionWithAxes)
 {
-    auto axis0 = sharedEquiDivision("axis0", 10, 0.0, 10.0);
-    auto axis1 = sharedEquiDivision("axis1", 20, 0.0, 20.0);
+    auto axis0 = sharedEquiDivision("axis0", Unit::other, 10, 0.0, 10.0);
+    auto axis1 = sharedEquiDivision("axis1", Unit::other, 20, 0.0, 20.0);
     SphericalDetector detector(std::array<std::shared_ptr<Scale>, 2>{axis0, axis1});
 
     // checking dimension and axes
@@ -58,7 +58,8 @@ TEST(SphericalDetectorTest, createDetectorMap)
 TEST(SphericalDetectorTest, regionOfInterest)
 {
     SphericalDetector detector(std::array<std::shared_ptr<Scale>, 2>{
-        sharedEquiDivision("axis0", 8, -3.0, 5.0), sharedEquiDivision("axis1", 4, 0.0, 4.0)});
+        sharedEquiDivision("axis0", Unit::other, 8, -3.0, 5.0),
+        sharedEquiDivision("axis1", Unit::other, 4, 0.0, 4.0)});
 
     // creating region of interest
     double xlow(-2.0), ylow(1.0), xup(4.0), yup(3.0);
diff --git a/Tests/Unit/GUI/Utils.cpp b/Tests/Unit/GUI/Utils.cpp
index bda0ff1a1385c189b9cab0ab0f80654309fc45af..507b85c8592041dabf7781f5a2f57f550ef1b6c3 100644
--- a/Tests/Unit/GUI/Utils.cpp
+++ b/Tests/Unit/GUI/Utils.cpp
@@ -42,9 +42,9 @@ void UTest::GUI::create_dir(const QString& dir_name)
 Datafield UTest::GUI::createData(double value, DIM n_dim)
 {
     std::vector<const Scale*> axes;
-    axes.emplace_back(newEquiDivision("x", nxsize, -1.0, 1.0));
+    axes.emplace_back(newEquiDivision("x", Unit::other, nxsize, -1.0, 1.0));
     if (n_dim == DIM::D2)
-        axes.emplace_back(newEquiDivision("y", nysize, 0.0, 2.0));
+        axes.emplace_back(newEquiDivision("y", Unit::other, nysize, 0.0, 2.0));
 
     Datafield result(std::move(axes));
     result.setAllTo(value);
diff --git a/Tests/Unit/Sim/FittingTestHelper.h b/Tests/Unit/Sim/FittingTestHelper.h
index cb8a51d37d6eb00d713aa503c216675d214ae3cf..3f52dc5542812017045a442da9bbc4a1a6c53100 100644
--- a/Tests/Unit/Sim/FittingTestHelper.h
+++ b/Tests/Unit/Sim/FittingTestHelper.h
@@ -51,8 +51,8 @@ public:
     std::unique_ptr<Datafield> createTestData(double value)
     {
         std::unique_ptr<Datafield> result(
-            new Datafield({newEquiDivision("phi_f (rad)", m_nx, m_xmin, m_xmax),
-                           newEquiDivision("alpha_f (rad)", m_ny, m_ymin, m_ymax)}));
+            new Datafield({newEquiDivision("phi_f", Unit::radian, m_nx, m_xmin, m_xmax),
+                           newEquiDivision("alpha_f", Unit::radian, m_ny, m_ymin, m_ymax)}));
         result->setAllTo(value);
         return result;
     }
diff --git a/Tests/Unit/Sim/SpecularScanTest.cpp b/Tests/Unit/Sim/SpecularScanTest.cpp
index b3b93b0f79d17b28d38eb0a154c103658e6cb2af..9ab39a532e553b00b626a1f301d69646677454e8 100644
--- a/Tests/Unit/Sim/SpecularScanTest.cpp
+++ b/Tests/Unit/Sim/SpecularScanTest.cpp
@@ -15,11 +15,12 @@ TEST(SpecularScanTest, AngularScanInit)
         EXPECT_EQ(scan.nScan(), axis.size());
     };
 
-    const Scale pointwise_axis = ListScan("inc_angles", std::vector<double>{0.1, 0.2, 0.3});
+    const Scale pointwise_axis =
+        ListScan("inc_angles", Unit::other, std::vector<double>{0.1, 0.2, 0.3});
     AlphaScan scan(pointwise_axis);
     check(scan, pointwise_axis);
 
-    const Scale fixed_axis = EquiDivision("inc_angles", 3, 0.1, 0.3);
+    const Scale fixed_axis = EquiDivision("inc_angles", Unit::other, 3, 0.1, 0.3);
     AlphaScan scan3(fixed_axis);
     check(scan3, fixed_axis);
 }
diff --git a/Wrap/Python/ba_plot.py b/Wrap/Python/ba_plot.py
index b612a0f7b3719fcb9a389f56f6af4300aeeec4b0..e8ff571c68cc31d38a2ee70c6a8e507e4938d9f7 100644
--- a/Wrap/Python/ba_plot.py
+++ b/Wrap/Python/ba_plot.py
@@ -169,7 +169,7 @@ def get_axes_labels(result):
     """
     labels = []
     for i in range(result.rank()):
-        labels.append(translate_axis_label(result.axis(i).axisName()))
+        labels.append(translate_axis_label(result.axis(i).axisLabel()))
 
     return labels
 
@@ -390,9 +390,9 @@ def plot_histogram(field, **kwargs):
     """
 
     if 'xlabel' not in kwargs:
-        kwargs['xlabel'] = translate_axis_label(field.xAxis().axisName())
+        kwargs['xlabel'] = translate_axis_label(field.xAxis().axisLabel())
     if 'ylabel' not in kwargs:
-        kwargs['ylabel'] = translate_axis_label(field.yAxis().axisName())
+        kwargs['ylabel'] = translate_axis_label(field.yAxis().axisLabel())
 
     axes_limits = [
         field.xAxis().min(),
diff --git a/auto/Examples/fit/specular/Specular1Par.py b/auto/Examples/fit/specular/Specular1Par.py
index c491a89e7aa733128b6fe154dcfaa7314379b79f..a35b3babaac7a16e9ae9866679fda45f4fb6a02e 100755
--- a/auto/Examples/fit/specular/Specular1Par.py
+++ b/auto/Examples/fit/specular/Specular1Par.py
@@ -47,7 +47,8 @@ def get_sample(P):
 
 
 def get_simulation(P):
-    scan = ba.AlphaScan(ba.ListScan("", exp_x))
+    alpha_units = ba.Unit_radian
+    scan = ba.AlphaScan(ba.ListScan("", alpha_units, exp_x))
     scan.setWavelength(1.54*angstrom)
     sample = get_sample(P)
 
diff --git a/auto/Examples/varia/Depthprobe1.py b/auto/Examples/varia/Depthprobe1.py
index 2d9c1a5f828f20e5f8ad92c5784ff0080a7ec415..a502f6f94c0eef42c27f93172ebcd77ff87f9f0b 100755
--- a/auto/Examples/varia/Depthprobe1.py
+++ b/auto/Examples/varia/Depthprobe1.py
@@ -27,7 +27,8 @@ def get_simulation(sample, flags):
     scan = ba.AlphaScan(n, 0*deg, 1*deg)
     scan.setWavelength(0.3*nm)
 
-    z_axis = ba.EquiDivision("z (nm)", n, -130*nm, 30*nm)
+    z_units = ba.Unit_nanometer
+    z_axis = ba.EquiDivision("z", z_units, n, -130*nm, 30*nm)
     simulation = ba.DepthprobeSimulation(scan, sample, z_axis, flags)
 
     return simulation
diff --git a/auto/Examples/varia/Resonator.py b/auto/Examples/varia/Resonator.py
index a4a83daa7939cef46f96904e9cad229115210af1..cfafd376604f620c9bf33df68d851eefa38a7db5 100755
--- a/auto/Examples/varia/Resonator.py
+++ b/auto/Examples/varia/Resonator.py
@@ -27,6 +27,7 @@ d_ang = 0.01*ba.deg  # spread width for incident angle
 #  depth position span
 z_min = -100*nm
 z_max = 100*nm
+z_units = ba.Unit_nanometer
 
 
 def get_sample():
@@ -73,7 +74,7 @@ def get_simulation(sample):
     footprint = ba.FootprintSquare(0.01)
     scan.setFootprint(footprint)
 
-    z_axis = ba.EquiDivision("z (nm)", nz, z_min, z_max)
+    z_axis = ba.EquiDivision("z", z_units, nz, z_min, z_max)
     simulation = ba.DepthprobeSimulation(scan, sample, z_axis)
 
     alpha_distr = ba.DistributionGaussian(0, d_ang, 25, 3.)
diff --git a/auto/Examples/varia/Transmission.py b/auto/Examples/varia/Transmission.py
index 2d17a31fcaca6ee0535bb6043870db45a9cd23d3..350090a9109331e877333f9fc16b99a3855d7feb 100755
--- a/auto/Examples/varia/Transmission.py
+++ b/auto/Examples/varia/Transmission.py
@@ -26,7 +26,8 @@ def simulate():
     scan = ba.AlphaScan(1, alpha, alpha)
     scan.setWavelength(0.3*nm)
 
-    z_axis = ba.EquiDivision("z (nm)", 1, -depth, -depth)
+    z_units = ba.Unit_nanometer
+    z_axis = ba.EquiDivision("z", z_units, 1, -depth, -depth)
     simulation = ba.DepthprobeSimulation(scan, get_sample(), z_axis, 0)
 
     result = simulation.simulate().flat()
diff --git a/auto/Examples/varia/TransmissionVsAlpha.py b/auto/Examples/varia/TransmissionVsAlpha.py
index 6015baba29ce6c746d36602ddf968e6820c6137b..3db173329e501b0c519b707570313124041a8a87 100755
--- a/auto/Examples/varia/TransmissionVsAlpha.py
+++ b/auto/Examples/varia/TransmissionVsAlpha.py
@@ -26,7 +26,8 @@ def simulate(depth):
     scan = ba.AlphaScan(n, alpha_max / n, alpha_max)
     scan.setWavelength(0.3*nm)
 
-    z_axis = ba.ListScan("z (nm)", [-depth])
+    z_units = ba.Unit_nanometer
+    z_axis = ba.ListScan("z", z_units, [-depth])
     simulation = ba.DepthprobeSimulation(scan, sample, z_axis, 0)
 
     result = simulation.simulate().flat()
diff --git a/auto/Examples/varia/TransmissionVsDepth.py b/auto/Examples/varia/TransmissionVsDepth.py
index 46378c88a3a60e7ecf8bbc1a089cdf8ace9d0066..14d79f541e037bb7538ff80019cbf1a1ec1850f5 100755
--- a/auto/Examples/varia/TransmissionVsDepth.py
+++ b/auto/Examples/varia/TransmissionVsDepth.py
@@ -24,7 +24,8 @@ def simulate(sample, alpha_i_deg):
     scan = ba.AlphaScan(1, alpha, alpha)
     scan.setWavelength(0.3*nm)
 
-    z_axis = ba.EquiDivision("z (nm)", 500, -130*nm, 30*nm)
+    z_units = ba.Unit_nanometer
+    z_axis = ba.EquiDivision("z", z_units, 500, -130*nm, 30*nm)
     simulation = ba.DepthprobeSimulation(scan, sample, z_axis, 0)
 
     result = simulation.simulate().flat()
diff --git a/auto/MiniExamples/fit/specular/Specular1Par.py b/auto/MiniExamples/fit/specular/Specular1Par.py
index c491a89e7aa733128b6fe154dcfaa7314379b79f..a35b3babaac7a16e9ae9866679fda45f4fb6a02e 100755
--- a/auto/MiniExamples/fit/specular/Specular1Par.py
+++ b/auto/MiniExamples/fit/specular/Specular1Par.py
@@ -47,7 +47,8 @@ def get_sample(P):
 
 
 def get_simulation(P):
-    scan = ba.AlphaScan(ba.ListScan("", exp_x))
+    alpha_units = ba.Unit_radian
+    scan = ba.AlphaScan(ba.ListScan("", alpha_units, exp_x))
     scan.setWavelength(1.54*angstrom)
     sample = get_sample(P)
 
diff --git a/auto/MiniExamples/varia/Depthprobe1.py b/auto/MiniExamples/varia/Depthprobe1.py
index af40d95162e62c8a76a32a081454764ace41d3be..47207358743d8ecd2ab690cf7c44f61397a45657 100755
--- a/auto/MiniExamples/varia/Depthprobe1.py
+++ b/auto/MiniExamples/varia/Depthprobe1.py
@@ -27,7 +27,8 @@ def get_simulation(sample, flags):
     scan = ba.AlphaScan(n, 0*deg, 1*deg)
     scan.setWavelength(0.3*nm)
 
-    z_axis = ba.EquiDivision("z (nm)", n, -130*nm, 30*nm)
+    z_units = ba.Unit_nanometer
+    z_axis = ba.EquiDivision("z", z_units, n, -130*nm, 30*nm)
     simulation = ba.DepthprobeSimulation(scan, sample, z_axis, flags)
 
     return simulation
diff --git a/auto/MiniExamples/varia/Resonator.py b/auto/MiniExamples/varia/Resonator.py
index 5b533c17f982bf2bffd2d22b40a71f460768e877..0acb35fd09614d06cdc089556f5c28955ae99ab8 100755
--- a/auto/MiniExamples/varia/Resonator.py
+++ b/auto/MiniExamples/varia/Resonator.py
@@ -27,6 +27,7 @@ d_ang = 0.01*ba.deg  # spread width for incident angle
 #  depth position span
 z_min = -100*nm
 z_max = 100*nm
+z_units = ba.Unit_nanometer
 
 
 def get_sample():
@@ -73,7 +74,7 @@ def get_simulation(sample):
     footprint = ba.FootprintSquare(0.01)
     scan.setFootprint(footprint)
 
-    z_axis = ba.EquiDivision("z (nm)", nz, z_min, z_max)
+    z_axis = ba.EquiDivision("z", z_units, nz, z_min, z_max)
     simulation = ba.DepthprobeSimulation(scan, sample, z_axis)
 
     alpha_distr = ba.DistributionGaussian(0, d_ang, 25, 3.)
diff --git a/auto/MiniExamples/varia/Transmission.py b/auto/MiniExamples/varia/Transmission.py
index f268b102b09f6ad972d0f9932a1dd840cb86f0db..e1c62e33723d2046563a7c7e1054a94ecfd92fbe 100755
--- a/auto/MiniExamples/varia/Transmission.py
+++ b/auto/MiniExamples/varia/Transmission.py
@@ -26,7 +26,8 @@ def simulate():
     scan = ba.AlphaScan(1, alpha, alpha)
     scan.setWavelength(0.3*nm)
 
-    z_axis = ba.EquiDivision("z (nm)", 1, -depth, -depth)
+    z_units = ba.Unit_nanometer
+    z_axis = ba.EquiDivision("z", z_units, 1, -depth, -depth)
     simulation = ba.DepthprobeSimulation(scan, get_sample(), z_axis, 0)
 
     result = simulation.simulate().flat()
diff --git a/auto/MiniExamples/varia/TransmissionVsAlpha.py b/auto/MiniExamples/varia/TransmissionVsAlpha.py
index 6e25a0381f57511b22665ff10c0f8464ddc6741f..9cbbd4fb200f13c1a40473d5fd85ebe5967e38ca 100755
--- a/auto/MiniExamples/varia/TransmissionVsAlpha.py
+++ b/auto/MiniExamples/varia/TransmissionVsAlpha.py
@@ -26,7 +26,8 @@ def simulate(depth):
     scan = ba.AlphaScan(n, alpha_max / n, alpha_max)
     scan.setWavelength(0.3*nm)
 
-    z_axis = ba.ListScan("z (nm)", [-depth])
+    z_units = ba.Unit_nanometer
+    z_axis = ba.ListScan("z", z_units, [-depth])
     simulation = ba.DepthprobeSimulation(scan, sample, z_axis, 0)
 
     result = simulation.simulate().flat()
diff --git a/auto/MiniExamples/varia/TransmissionVsDepth.py b/auto/MiniExamples/varia/TransmissionVsDepth.py
index 5d41f4bd2d6c27f86c1215dbd13f784eaa9bfc45..3b8fadfc848c7585f7366547a268802a41eb802f 100755
--- a/auto/MiniExamples/varia/TransmissionVsDepth.py
+++ b/auto/MiniExamples/varia/TransmissionVsDepth.py
@@ -24,7 +24,8 @@ def simulate(sample, alpha_i_deg):
     scan = ba.AlphaScan(1, alpha, alpha)
     scan.setWavelength(0.3*nm)
 
-    z_axis = ba.EquiDivision("z (nm)", 20, -130*nm, 30*nm)
+    z_units = ba.Unit_nanometer
+    z_axis = ba.EquiDivision("z", z_units, 20, -130*nm, 30*nm)
     simulation = ba.DepthprobeSimulation(scan, sample, z_axis, 0)
 
     result = simulation.simulate().flat()
diff --git a/auto/Wrap/libBornAgainBase.py b/auto/Wrap/libBornAgainBase.py
index 8ebee9537a1d65baca8a847db056f84b95687a18..e8038f0b2119f7ae37159c89e577660d66d0caed 100644
--- a/auto/Wrap/libBornAgainBase.py
+++ b/auto/Wrap/libBornAgainBase.py
@@ -1787,6 +1787,10 @@ Unit_unitless = _libBornAgainBase.Unit_unitless
 Unit_other = _libBornAgainBase.Unit_other
 
 
+def inOneGroup(a, b):
+    r"""inOneGroup(Unit a, Unit b) -> bool"""
+    return _libBornAgainBase.inOneGroup(a, b)
+
 def plainName(units):
     r"""plainName(Unit units) -> std::string"""
     return _libBornAgainBase.plainName(units)
@@ -1912,22 +1916,26 @@ class Scale(object):
     thisown = property(lambda x: x.this.own(), lambda x, v: x.this.own(v), doc="The membership flag")
     __repr__ = _swig_repr
 
-    def __init__(self, name, bins):
-        r"""__init__(Scale self, std::string name, std::vector< Bin1D,std::allocator< Bin1D > > const & bins) -> Scale"""
-        _libBornAgainBase.Scale_swiginit(self, _libBornAgainBase.new_Scale(name, bins))
+    def __init__(self, name, units, bins):
+        r"""__init__(Scale self, std::string name, Unit units, std::vector< Bin1D,std::allocator< Bin1D > > const & bins) -> Scale"""
+        _libBornAgainBase.Scale_swiginit(self, _libBornAgainBase.new_Scale(name, units, bins))
 
     def clone(self):
         r"""clone(Scale self) -> Scale"""
         return _libBornAgainBase.Scale_clone(self)
 
-    def setAxisName(self, name):
-        r"""setAxisName(Scale self, std::string name)"""
-        return _libBornAgainBase.Scale_setAxisName(self, name)
-
     def axisName(self):
         r"""axisName(Scale self) -> std::string"""
         return _libBornAgainBase.Scale_axisName(self)
 
+    def axisUnits(self):
+        r"""axisUnits(Scale self) -> Unit"""
+        return _libBornAgainBase.Scale_axisUnits(self)
+
+    def axisLabel(self):
+        r"""axisLabel(Scale self) -> std::string"""
+        return _libBornAgainBase.Scale_axisLabel(self)
+
     def size(self):
         r"""size(Scale self) -> size_t"""
         return _libBornAgainBase.Scale_size(self)
@@ -1999,29 +2007,25 @@ class Scale(object):
         r"""__eq__(Scale self, Scale right) -> bool"""
         return _libBornAgainBase.Scale___eq__(self, right)
 
-    def unit(self):
-        r"""unit(Scale self) -> std::string"""
-        return _libBornAgainBase.Scale_unit(self)
-
-    def plottableScale(self):
-        r"""plottableScale(Scale self) -> Scale"""
-        return _libBornAgainBase.Scale_plottableScale(self)
+    def convertedScale(self, newUnits):
+        r"""convertedScale(Scale self, Unit newUnits) -> Scale"""
+        return _libBornAgainBase.Scale_convertedScale(self, newUnits)
     __swig_destroy__ = _libBornAgainBase.delete_Scale
 
 # Register Scale in _libBornAgainBase:
 _libBornAgainBase.Scale_swigregister(Scale)
 
-def ListScan(name, points):
-    r"""ListScan(std::string const & name, vdouble1d_t points) -> Scale"""
-    return _libBornAgainBase.ListScan(name, points)
+def ListScan(name, units, points):
+    r"""ListScan(std::string const & name, Unit units, vdouble1d_t points) -> Scale"""
+    return _libBornAgainBase.ListScan(name, units, points)
 
-def EquiDivision(name, N, start, end):
-    r"""EquiDivision(std::string const & name, size_t N, double start, double end) -> Scale"""
-    return _libBornAgainBase.EquiDivision(name, N, start, end)
+def EquiDivision(name, units, N, start, end):
+    r"""EquiDivision(std::string const & name, Unit units, size_t N, double start, double end) -> Scale"""
+    return _libBornAgainBase.EquiDivision(name, units, N, start, end)
 
-def EquiScan(name, N, start, end):
-    r"""EquiScan(std::string const & name, size_t N, double start, double end) -> Scale"""
-    return _libBornAgainBase.EquiScan(name, N, start, end)
+def EquiScan(name, units, N, start, end):
+    r"""EquiScan(std::string const & name, Unit units, size_t N, double start, double end) -> Scale"""
+    return _libBornAgainBase.EquiScan(name, units, N, start, end)
 class Frame(ICloneable):
     r"""Proxy of C++ Frame class."""
 
@@ -2089,9 +2093,9 @@ class Frame(ICloneable):
         r"""__eq__(Frame self, Frame arg2) -> bool"""
         return _libBornAgainBase.Frame___eq__(self, arg2)
 
-    def plottableFrame(self):
-        r"""plottableFrame(Frame self) -> Frame"""
-        return _libBornAgainBase.Frame_plottableFrame(self)
+    def convertedFrame(self, newUnits):
+        r"""convertedFrame(Frame self, std::vector< Unit,std::allocator< Unit > > const & newUnits) -> Frame"""
+        return _libBornAgainBase.Frame_convertedFrame(self, newUnits)
 
     def flat(self):
         r"""flat(Frame self) -> Frame"""
diff --git a/auto/Wrap/libBornAgainBase_wrap.cpp b/auto/Wrap/libBornAgainBase_wrap.cpp
index b0964be3eb451ee928d9920ce6834c7b2a123c7a..994a800e2c09f73b2b8df29686bb6ea520cf363e 100644
--- a/auto/Wrap/libBornAgainBase_wrap.cpp
+++ b/auto/Wrap/libBornAgainBase_wrap.cpp
@@ -24679,6 +24679,36 @@ SWIGINTERN PyObject *Swig_var_internalDepthUnit_get(void) {
 }
 
 
+SWIGINTERN PyObject *_wrap_inOneGroup(PyObject *self, PyObject *args) {
+  PyObject *resultobj = 0;
+  Unit arg1 ;
+  Unit arg2 ;
+  int val1 ;
+  int ecode1 = 0 ;
+  int val2 ;
+  int ecode2 = 0 ;
+  PyObject *swig_obj[2] ;
+  bool result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "inOneGroup", 2, 2, swig_obj)) SWIG_fail;
+  ecode1 = SWIG_AsVal_int(swig_obj[0], &val1);
+  if (!SWIG_IsOK(ecode1)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "inOneGroup" "', argument " "1"" of type '" "Unit""'");
+  } 
+  arg1 = static_cast< Unit >(val1);
+  ecode2 = SWIG_AsVal_int(swig_obj[1], &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "inOneGroup" "', argument " "2"" of type '" "Unit""'");
+  } 
+  arg2 = static_cast< Unit >(val2);
+  result = (bool)Units::inOneGroup(arg1,arg2);
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_plainName(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Unit arg1 ;
@@ -25581,13 +25611,16 @@ SWIGINTERN PyObject *Bin1D_swigregister(PyObject *SWIGUNUSEDPARM(self), PyObject
 SWIGINTERN PyObject *_wrap_new_Scale(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   std::string arg1 ;
-  std::vector< Bin1D,std::allocator< Bin1D > > *arg2 = 0 ;
-  void *argp2 = 0 ;
-  int res2 = 0 ;
-  PyObject *swig_obj[2] ;
+  Unit arg2 ;
+  std::vector< Bin1D,std::allocator< Bin1D > > *arg3 = 0 ;
+  int val2 ;
+  int ecode2 = 0 ;
+  void *argp3 = 0 ;
+  int res3 = 0 ;
+  PyObject *swig_obj[3] ;
   Scale *result = 0 ;
   
-  if (!SWIG_Python_UnpackTuple(args, "new_Scale", 2, 2, swig_obj)) SWIG_fail;
+  if (!SWIG_Python_UnpackTuple(args, "new_Scale", 3, 3, swig_obj)) SWIG_fail;
   {
     std::string *ptr = (std::string *)0;
     int res = SWIG_AsPtr_std_string(swig_obj[0], &ptr);
@@ -25597,15 +25630,20 @@ SWIGINTERN PyObject *_wrap_new_Scale(PyObject *self, PyObject *args) {
     arg1 = *ptr;
     if (SWIG_IsNewObj(res)) delete ptr;
   }
-  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_std__vectorT_Bin1D_std__allocatorT_Bin1D_t_t,  0  | 0);
-  if (!SWIG_IsOK(res2)) {
-    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "new_Scale" "', argument " "2"" of type '" "std::vector< Bin1D,std::allocator< Bin1D > > const &""'"); 
+  ecode2 = SWIG_AsVal_int(swig_obj[1], &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "new_Scale" "', argument " "2"" of type '" "Unit""'");
+  } 
+  arg2 = static_cast< Unit >(val2);
+  res3 = SWIG_ConvertPtr(swig_obj[2], &argp3, SWIGTYPE_p_std__vectorT_Bin1D_std__allocatorT_Bin1D_t_t,  0  | 0);
+  if (!SWIG_IsOK(res3)) {
+    SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "new_Scale" "', argument " "3"" of type '" "std::vector< Bin1D,std::allocator< Bin1D > > const &""'"); 
   }
-  if (!argp2) {
-    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "new_Scale" "', argument " "2"" of type '" "std::vector< Bin1D,std::allocator< Bin1D > > const &""'"); 
+  if (!argp3) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "new_Scale" "', argument " "3"" of type '" "std::vector< Bin1D,std::allocator< Bin1D > > const &""'"); 
   }
-  arg2 = reinterpret_cast< std::vector< Bin1D,std::allocator< Bin1D > > * >(argp2);
-  result = (Scale *)new Scale(arg1,(std::vector< Bin1D,std::allocator< Bin1D > > const &)*arg2);
+  arg3 = reinterpret_cast< std::vector< Bin1D,std::allocator< Bin1D > > * >(argp3);
+  result = (Scale *)new Scale(arg1,arg2,(std::vector< Bin1D,std::allocator< Bin1D > > const &)*arg3);
   resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Scale, SWIG_POINTER_NEW |  0 );
   return resultobj;
 fail:
@@ -25636,38 +25674,53 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_Scale_setAxisName(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_Scale_axisName(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Scale *arg1 = (Scale *) 0 ;
-  std::string arg2 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  PyObject *swig_obj[2] ;
+  PyObject *swig_obj[1] ;
+  std::string result;
   
-  if (!SWIG_Python_UnpackTuple(args, "Scale_setAxisName", 2, 2, swig_obj)) SWIG_fail;
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Scale, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Scale_setAxisName" "', argument " "1"" of type '" "Scale *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Scale_axisName" "', argument " "1"" of type '" "Scale const *""'"); 
   }
   arg1 = reinterpret_cast< Scale * >(argp1);
-  {
-    std::string *ptr = (std::string *)0;
-    int res = SWIG_AsPtr_std_string(swig_obj[1], &ptr);
-    if (!SWIG_IsOK(res) || !ptr) {
-      SWIG_exception_fail(SWIG_ArgError((ptr ? res : SWIG_TypeError)), "in method '" "Scale_setAxisName" "', argument " "2"" of type '" "std::string""'"); 
-    }
-    arg2 = *ptr;
-    if (SWIG_IsNewObj(res)) delete ptr;
+  result = ((Scale const *)arg1)->axisName();
+  resultobj = SWIG_From_std_string(static_cast< std::string >(result));
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_Scale_axisUnits(PyObject *self, PyObject *args) {
+  PyObject *resultobj = 0;
+  Scale *arg1 = (Scale *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  Unit result;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Scale, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Scale_axisUnits" "', argument " "1"" of type '" "Scale const *""'"); 
   }
-  (arg1)->setAxisName(arg2);
-  resultobj = SWIG_Py_Void();
+  arg1 = reinterpret_cast< Scale * >(argp1);
+  result = (Unit)((Scale const *)arg1)->axisUnits();
+  resultobj = SWIG_From_int(static_cast< int >(result));
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_Scale_axisName(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_Scale_axisLabel(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Scale *arg1 = (Scale *) 0 ;
   void *argp1 = 0 ;
@@ -25679,10 +25732,10 @@ SWIGINTERN PyObject *_wrap_Scale_axisName(PyObject *self, PyObject *args) {
   swig_obj[0] = args;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Scale, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Scale_axisName" "', argument " "1"" of type '" "Scale const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Scale_axisLabel" "', argument " "1"" of type '" "Scale const *""'"); 
   }
   arg1 = reinterpret_cast< Scale * >(argp1);
-  result = ((Scale const *)arg1)->axisName();
+  result = ((Scale const *)arg1)->axisLabel();
   resultobj = SWIG_From_std_string(static_cast< std::string >(result));
   return resultobj;
 fail:
@@ -26218,45 +26271,29 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_Scale_unit(PyObject *self, PyObject *args) {
-  PyObject *resultobj = 0;
-  Scale *arg1 = (Scale *) 0 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  PyObject *swig_obj[1] ;
-  std::string result;
-  
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Scale, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Scale_unit" "', argument " "1"" of type '" "Scale const *""'"); 
-  }
-  arg1 = reinterpret_cast< Scale * >(argp1);
-  result = ((Scale const *)arg1)->unit();
-  resultobj = SWIG_From_std_string(static_cast< std::string >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_Scale_plottableScale(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_Scale_convertedScale(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Scale *arg1 = (Scale *) 0 ;
+  Unit arg2 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  PyObject *swig_obj[1] ;
+  int val2 ;
+  int ecode2 = 0 ;
+  PyObject *swig_obj[2] ;
   SwigValueWrapper< Scale > result;
   
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
+  if (!SWIG_Python_UnpackTuple(args, "Scale_convertedScale", 2, 2, swig_obj)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Scale, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Scale_plottableScale" "', argument " "1"" of type '" "Scale const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Scale_convertedScale" "', argument " "1"" of type '" "Scale const *""'"); 
   }
   arg1 = reinterpret_cast< Scale * >(argp1);
-  result = ((Scale const *)arg1)->plottableScale();
+  ecode2 = SWIG_AsVal_int(swig_obj[1], &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "Scale_convertedScale" "', argument " "2"" of type '" "Unit""'");
+  } 
+  arg2 = static_cast< Unit >(val2);
+  result = ((Scale const *)arg1)->convertedScale(arg2);
   resultobj = SWIG_NewPointerObj((new Scale(result)), SWIGTYPE_p_Scale, SWIG_POINTER_OWN |  0 );
   return resultobj;
 fail:
@@ -26300,13 +26337,16 @@ SWIGINTERN PyObject *Scale_swiginit(PyObject *SWIGUNUSEDPARM(self), PyObject *ar
 SWIGINTERN PyObject *_wrap_ListScan(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   std::string *arg1 = 0 ;
-  std::vector< double,std::allocator< double > > *arg2 = 0 ;
+  Unit arg2 ;
+  std::vector< double,std::allocator< double > > *arg3 = 0 ;
   int res1 = SWIG_OLDOBJ ;
-  int res2 = SWIG_OLDOBJ ;
-  PyObject *swig_obj[2] ;
+  int val2 ;
+  int ecode2 = 0 ;
+  int res3 = SWIG_OLDOBJ ;
+  PyObject *swig_obj[3] ;
   SwigValueWrapper< Scale > result;
   
-  if (!SWIG_Python_UnpackTuple(args, "ListScan", 2, 2, swig_obj)) SWIG_fail;
+  if (!SWIG_Python_UnpackTuple(args, "ListScan", 3, 3, swig_obj)) SWIG_fail;
   {
     std::string *ptr = (std::string *)0;
     res1 = SWIG_AsPtr_std_string(swig_obj[0], &ptr);
@@ -26318,25 +26358,30 @@ SWIGINTERN PyObject *_wrap_ListScan(PyObject *self, PyObject *args) {
     }
     arg1 = ptr;
   }
+  ecode2 = SWIG_AsVal_int(swig_obj[1], &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "ListScan" "', argument " "2"" of type '" "Unit""'");
+  } 
+  arg2 = static_cast< Unit >(val2);
   {
     std::vector< double,std::allocator< double > > *ptr = (std::vector< double,std::allocator< double > > *)0;
-    res2 = swig::asptr(swig_obj[1], &ptr);
-    if (!SWIG_IsOK(res2)) {
-      SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "ListScan" "', argument " "2"" of type '" "std::vector< double,std::allocator< double > > const &""'"); 
+    res3 = swig::asptr(swig_obj[2], &ptr);
+    if (!SWIG_IsOK(res3)) {
+      SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "ListScan" "', argument " "3"" of type '" "std::vector< double,std::allocator< double > > const &""'"); 
     }
     if (!ptr) {
-      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "ListScan" "', argument " "2"" of type '" "std::vector< double,std::allocator< double > > const &""'"); 
+      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "ListScan" "', argument " "3"" of type '" "std::vector< double,std::allocator< double > > const &""'"); 
     }
-    arg2 = ptr;
+    arg3 = ptr;
   }
-  result = ListScan((std::string const &)*arg1,(std::vector< double,std::allocator< double > > const &)*arg2);
+  result = ListScan((std::string const &)*arg1,arg2,(std::vector< double,std::allocator< double > > const &)*arg3);
   resultobj = SWIG_NewPointerObj((new Scale(result)), SWIGTYPE_p_Scale, SWIG_POINTER_OWN |  0 );
   if (SWIG_IsNewObj(res1)) delete arg1;
-  if (SWIG_IsNewObj(res2)) delete arg2;
+  if (SWIG_IsNewObj(res3)) delete arg3;
   return resultobj;
 fail:
   if (SWIG_IsNewObj(res1)) delete arg1;
-  if (SWIG_IsNewObj(res2)) delete arg2;
+  if (SWIG_IsNewObj(res3)) delete arg3;
   return NULL;
 }
 
@@ -26344,20 +26389,23 @@ fail:
 SWIGINTERN PyObject *_wrap_EquiDivision(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   std::string *arg1 = 0 ;
-  size_t arg2 ;
-  double arg3 ;
+  Unit arg2 ;
+  size_t arg3 ;
   double arg4 ;
+  double arg5 ;
   int res1 = SWIG_OLDOBJ ;
-  size_t val2 ;
+  int val2 ;
   int ecode2 = 0 ;
-  double val3 ;
+  size_t val3 ;
   int ecode3 = 0 ;
   double val4 ;
   int ecode4 = 0 ;
-  PyObject *swig_obj[4] ;
+  double val5 ;
+  int ecode5 = 0 ;
+  PyObject *swig_obj[5] ;
   SwigValueWrapper< Scale > result;
   
-  if (!SWIG_Python_UnpackTuple(args, "EquiDivision", 4, 4, swig_obj)) SWIG_fail;
+  if (!SWIG_Python_UnpackTuple(args, "EquiDivision", 5, 5, swig_obj)) SWIG_fail;
   {
     std::string *ptr = (std::string *)0;
     res1 = SWIG_AsPtr_std_string(swig_obj[0], &ptr);
@@ -26369,22 +26417,27 @@ SWIGINTERN PyObject *_wrap_EquiDivision(PyObject *self, PyObject *args) {
     }
     arg1 = ptr;
   }
-  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
+  ecode2 = SWIG_AsVal_int(swig_obj[1], &val2);
   if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "EquiDivision" "', argument " "2"" of type '" "size_t""'");
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "EquiDivision" "', argument " "2"" of type '" "Unit""'");
   } 
-  arg2 = static_cast< size_t >(val2);
-  ecode3 = SWIG_AsVal_double(swig_obj[2], &val3);
+  arg2 = static_cast< Unit >(val2);
+  ecode3 = SWIG_AsVal_size_t(swig_obj[2], &val3);
   if (!SWIG_IsOK(ecode3)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "EquiDivision" "', argument " "3"" of type '" "double""'");
+    SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "EquiDivision" "', argument " "3"" of type '" "size_t""'");
   } 
-  arg3 = static_cast< double >(val3);
+  arg3 = static_cast< size_t >(val3);
   ecode4 = SWIG_AsVal_double(swig_obj[3], &val4);
   if (!SWIG_IsOK(ecode4)) {
     SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "EquiDivision" "', argument " "4"" of type '" "double""'");
   } 
   arg4 = static_cast< double >(val4);
-  result = EquiDivision((std::string const &)*arg1,SWIG_STD_MOVE(arg2),arg3,arg4);
+  ecode5 = SWIG_AsVal_double(swig_obj[4], &val5);
+  if (!SWIG_IsOK(ecode5)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode5), "in method '" "EquiDivision" "', argument " "5"" of type '" "double""'");
+  } 
+  arg5 = static_cast< double >(val5);
+  result = EquiDivision((std::string const &)*arg1,arg2,SWIG_STD_MOVE(arg3),arg4,arg5);
   resultobj = SWIG_NewPointerObj((new Scale(result)), SWIGTYPE_p_Scale, SWIG_POINTER_OWN |  0 );
   if (SWIG_IsNewObj(res1)) delete arg1;
   return resultobj;
@@ -26397,20 +26450,23 @@ fail:
 SWIGINTERN PyObject *_wrap_EquiScan(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   std::string *arg1 = 0 ;
-  size_t arg2 ;
-  double arg3 ;
+  Unit arg2 ;
+  size_t arg3 ;
   double arg4 ;
+  double arg5 ;
   int res1 = SWIG_OLDOBJ ;
-  size_t val2 ;
+  int val2 ;
   int ecode2 = 0 ;
-  double val3 ;
+  size_t val3 ;
   int ecode3 = 0 ;
   double val4 ;
   int ecode4 = 0 ;
-  PyObject *swig_obj[4] ;
+  double val5 ;
+  int ecode5 = 0 ;
+  PyObject *swig_obj[5] ;
   SwigValueWrapper< Scale > result;
   
-  if (!SWIG_Python_UnpackTuple(args, "EquiScan", 4, 4, swig_obj)) SWIG_fail;
+  if (!SWIG_Python_UnpackTuple(args, "EquiScan", 5, 5, swig_obj)) SWIG_fail;
   {
     std::string *ptr = (std::string *)0;
     res1 = SWIG_AsPtr_std_string(swig_obj[0], &ptr);
@@ -26422,22 +26478,27 @@ SWIGINTERN PyObject *_wrap_EquiScan(PyObject *self, PyObject *args) {
     }
     arg1 = ptr;
   }
-  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
+  ecode2 = SWIG_AsVal_int(swig_obj[1], &val2);
   if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "EquiScan" "', argument " "2"" of type '" "size_t""'");
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "EquiScan" "', argument " "2"" of type '" "Unit""'");
   } 
-  arg2 = static_cast< size_t >(val2);
-  ecode3 = SWIG_AsVal_double(swig_obj[2], &val3);
+  arg2 = static_cast< Unit >(val2);
+  ecode3 = SWIG_AsVal_size_t(swig_obj[2], &val3);
   if (!SWIG_IsOK(ecode3)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "EquiScan" "', argument " "3"" of type '" "double""'");
+    SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "EquiScan" "', argument " "3"" of type '" "size_t""'");
   } 
-  arg3 = static_cast< double >(val3);
+  arg3 = static_cast< size_t >(val3);
   ecode4 = SWIG_AsVal_double(swig_obj[3], &val4);
   if (!SWIG_IsOK(ecode4)) {
     SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "EquiScan" "', argument " "4"" of type '" "double""'");
   } 
   arg4 = static_cast< double >(val4);
-  result = EquiScan((std::string const &)*arg1,SWIG_STD_MOVE(arg2),arg3,arg4);
+  ecode5 = SWIG_AsVal_double(swig_obj[4], &val5);
+  if (!SWIG_IsOK(ecode5)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode5), "in method '" "EquiScan" "', argument " "5"" of type '" "double""'");
+  } 
+  arg5 = static_cast< double >(val5);
+  result = EquiScan((std::string const &)*arg1,arg2,SWIG_STD_MOVE(arg3),arg4,arg5);
   resultobj = SWIG_NewPointerObj((new Scale(result)), SWIGTYPE_p_Scale, SWIG_POINTER_OWN |  0 );
   if (SWIG_IsNewObj(res1)) delete arg1;
   return resultobj;
@@ -27009,22 +27070,32 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_Frame_plottableFrame(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_Frame_convertedFrame(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Frame *arg1 = (Frame *) 0 ;
+  std::vector< Unit,std::allocator< Unit > > *arg2 = 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  PyObject *swig_obj[1] ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
   Frame *result = 0 ;
   
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
+  if (!SWIG_Python_UnpackTuple(args, "Frame_convertedFrame", 2, 2, swig_obj)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Frame, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Frame_plottableFrame" "', argument " "1"" of type '" "Frame const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Frame_convertedFrame" "', argument " "1"" of type '" "Frame const *""'"); 
   }
   arg1 = reinterpret_cast< Frame * >(argp1);
-  result = (Frame *)((Frame const *)arg1)->plottableFrame();
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_std__vectorT_Unit_std__allocatorT_Unit_t_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "Frame_convertedFrame" "', argument " "2"" of type '" "std::vector< Unit,std::allocator< Unit > > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "Frame_convertedFrame" "', argument " "2"" of type '" "std::vector< Unit,std::allocator< Unit > > const &""'"); 
+  }
+  arg2 = reinterpret_cast< std::vector< Unit,std::allocator< Unit > > * >(argp2);
+  result = (Frame *)((Frame const *)arg1)->convertedFrame((std::vector< Unit,std::allocator< Unit > > const &)*arg2);
   resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Frame, 0 |  0 );
   return resultobj;
 fail:
@@ -29052,6 +29123,7 @@ static PyMethodDef SwigMethods[] = {
 	 { "delete_ThreadInfo", _wrap_delete_ThreadInfo, METH_O, "delete_ThreadInfo(ThreadInfo self)"},
 	 { "ThreadInfo_swigregister", ThreadInfo_swigregister, METH_O, NULL},
 	 { "ThreadInfo_swiginit", ThreadInfo_swiginit, METH_VARARGS, NULL},
+	 { "inOneGroup", _wrap_inOneGroup, METH_VARARGS, "inOneGroup(Unit a, Unit b) -> bool"},
 	 { "plainName", _wrap_plainName, METH_O, "plainName(Unit units) -> std::string"},
 	 { "symbolicName", _wrap_symbolicName, METH_O, "symbolicName(Unit unit) -> std::string"},
 	 { "unitExists", _wrap_unitExists, METH_O, "unitExists(std::string units) -> bool"},
@@ -29078,10 +29150,11 @@ static PyMethodDef SwigMethods[] = {
 	 { "Bin1D_clipped_or_nil", _wrap_Bin1D_clipped_or_nil, METH_VARARGS, "Bin1D_clipped_or_nil(Bin1D self, double lower, double upper) -> std::optional< Bin1D >"},
 	 { "delete_Bin1D", _wrap_delete_Bin1D, METH_O, "delete_Bin1D(Bin1D self)"},
 	 { "Bin1D_swigregister", Bin1D_swigregister, METH_O, NULL},
-	 { "new_Scale", _wrap_new_Scale, METH_VARARGS, "new_Scale(std::string name, std::vector< Bin1D,std::allocator< Bin1D > > const & bins) -> Scale"},
+	 { "new_Scale", _wrap_new_Scale, METH_VARARGS, "new_Scale(std::string name, Unit units, std::vector< Bin1D,std::allocator< Bin1D > > const & bins) -> Scale"},
 	 { "Scale_clone", _wrap_Scale_clone, METH_O, "Scale_clone(Scale self) -> Scale"},
-	 { "Scale_setAxisName", _wrap_Scale_setAxisName, METH_VARARGS, "Scale_setAxisName(Scale self, std::string name)"},
 	 { "Scale_axisName", _wrap_Scale_axisName, METH_O, "Scale_axisName(Scale self) -> std::string"},
+	 { "Scale_axisUnits", _wrap_Scale_axisUnits, METH_O, "Scale_axisUnits(Scale self) -> Unit"},
+	 { "Scale_axisLabel", _wrap_Scale_axisLabel, METH_O, "Scale_axisLabel(Scale self) -> std::string"},
 	 { "Scale_size", _wrap_Scale_size, METH_O, "Scale_size(Scale self) -> size_t"},
 	 { "Scale_min", _wrap_Scale_min, METH_O, "Scale_min(Scale self) -> double"},
 	 { "Scale_max", _wrap_Scale_max, METH_O, "Scale_max(Scale self) -> double"},
@@ -29102,14 +29175,13 @@ static PyMethodDef SwigMethods[] = {
 		"Scale_clipped(Scale self, pvacuum_double_t bounds) -> Scale\n"
 		""},
 	 { "Scale___eq__", _wrap_Scale___eq__, METH_VARARGS, "Scale___eq__(Scale self, Scale right) -> bool"},
-	 { "Scale_unit", _wrap_Scale_unit, METH_O, "Scale_unit(Scale self) -> std::string"},
-	 { "Scale_plottableScale", _wrap_Scale_plottableScale, METH_O, "Scale_plottableScale(Scale self) -> Scale"},
+	 { "Scale_convertedScale", _wrap_Scale_convertedScale, METH_VARARGS, "Scale_convertedScale(Scale self, Unit newUnits) -> Scale"},
 	 { "delete_Scale", _wrap_delete_Scale, METH_O, "delete_Scale(Scale self)"},
 	 { "Scale_swigregister", Scale_swigregister, METH_O, NULL},
 	 { "Scale_swiginit", Scale_swiginit, METH_VARARGS, NULL},
-	 { "ListScan", _wrap_ListScan, METH_VARARGS, "ListScan(std::string const & name, vdouble1d_t points) -> Scale"},
-	 { "EquiDivision", _wrap_EquiDivision, METH_VARARGS, "EquiDivision(std::string const & name, size_t N, double start, double end) -> Scale"},
-	 { "EquiScan", _wrap_EquiScan, METH_VARARGS, "EquiScan(std::string const & name, size_t N, double start, double end) -> Scale"},
+	 { "ListScan", _wrap_ListScan, METH_VARARGS, "ListScan(std::string const & name, Unit units, vdouble1d_t points) -> Scale"},
+	 { "EquiDivision", _wrap_EquiDivision, METH_VARARGS, "EquiDivision(std::string const & name, Unit units, size_t N, double start, double end) -> Scale"},
+	 { "EquiScan", _wrap_EquiScan, METH_VARARGS, "EquiScan(std::string const & name, Unit units, size_t N, double start, double end) -> Scale"},
 	 { "new_Frame", _wrap_new_Frame, METH_VARARGS, "\n"
 		"Frame(std::vector< Scale const *,std::allocator< Scale const * > > && axes)\n"
 		"Frame(Scale const *&& ax0)\n"
@@ -29129,7 +29201,7 @@ static PyMethodDef SwigMethods[] = {
 	 { "Frame_toGlobalIndex", _wrap_Frame_toGlobalIndex, METH_VARARGS, "Frame_toGlobalIndex(Frame self, std::vector< unsigned int,std::allocator< unsigned int > > const & axes_indices) -> size_t"},
 	 { "Frame_hasSameSizes", _wrap_Frame_hasSameSizes, METH_VARARGS, "Frame_hasSameSizes(Frame self, Frame arg2) -> bool"},
 	 { "Frame___eq__", _wrap_Frame___eq__, METH_VARARGS, "Frame___eq__(Frame self, Frame arg2) -> bool"},
-	 { "Frame_plottableFrame", _wrap_Frame_plottableFrame, METH_O, "Frame_plottableFrame(Frame self) -> Frame"},
+	 { "Frame_convertedFrame", _wrap_Frame_convertedFrame, METH_VARARGS, "Frame_convertedFrame(Frame self, std::vector< Unit,std::allocator< Unit > > const & newUnits) -> Frame"},
 	 { "Frame_flat", _wrap_Frame_flat, METH_O, "Frame_flat(Frame self) -> Frame"},
 	 { "Frame_swigregister", Frame_swigregister, METH_O, NULL},
 	 { "Frame_swiginit", Frame_swiginit, METH_VARARGS, NULL},
@@ -29239,7 +29311,7 @@ static swig_type_info _swigt__p_std__optionalT_Bin1D_t = {"_p_std__optionalT_Bin
 static swig_type_info _swigt__p_std__pairT_double_double_t = {"_p_std__pairT_double_double_t", "std::pair< double,double > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_Bin1D_std__allocatorT_Bin1D_t_t = {"_p_std__vectorT_Bin1D_std__allocatorT_Bin1D_t_t", "std::vector< Bin1D,std::allocator< Bin1D > > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_Scale_const_p_std__allocatorT_Scale_const_p_t_t = {"_p_std__vectorT_Scale_const_p_std__allocatorT_Scale_const_p_t_t", "std::vector< Scale const *,std::allocator< Scale const * > > *", 0, 0, (void*)0, 0};
-static swig_type_info _swigt__p_std__vectorT_Unit_std__allocatorT_Unit_t_t = {"_p_std__vectorT_Unit_std__allocatorT_Unit_t_t", "Units::Vec *|std::vector< enum Unit,std::allocator< enum Unit > > *", 0, 0, (void*)0, 0};
+static swig_type_info _swigt__p_std__vectorT_Unit_std__allocatorT_Unit_t_t = {"_p_std__vectorT_Unit_std__allocatorT_Unit_t_t", "std::vector< Unit,std::allocator< Unit > > *|std::vector< enum Unit,std::allocator< enum Unit > > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_double_std__allocatorT_double_t_t = {"_p_std__vectorT_double_std__allocatorT_double_t_t", "std::vector< double,std::allocator< double > > *|std::vector< double > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_int_std__allocatorT_int_t_t = {"_p_std__vectorT_int_std__allocatorT_int_t_t", "std::vector< int,std::allocator< int > > *|std::vector< int > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_std__complexT_double_t_std__allocatorT_std__complexT_double_t_t_t = {"_p_std__vectorT_std__complexT_double_t_std__allocatorT_std__complexT_double_t_t_t", "std::vector< std::complex< double >,std::allocator< std::complex< double > > > *|std::vector< std::complex< double > > *", 0, 0, (void*)0, 0};
diff --git a/auto/Wrap/libBornAgainDevice.py b/auto/Wrap/libBornAgainDevice.py
index e4b7dc675390a63fa27657443abbc5e7478f3c93..502b5b32798ce791e24ff3f02f46acaffa1dc5b4 100644
--- a/auto/Wrap/libBornAgainDevice.py
+++ b/auto/Wrap/libBornAgainDevice.py
@@ -2137,9 +2137,12 @@ class Datafield(object):
         r"""minVal(Datafield self) -> double"""
         return _libBornAgainDevice.Datafield_minVal(self)
 
-    def plottableField(self):
-        r"""plottableField(Datafield self) -> Datafield"""
-        return _libBornAgainDevice.Datafield_plottableField(self)
+    def plottableField(self, *args):
+        r"""
+        plottableField(Datafield self) -> Datafield
+        plottableField(Datafield self, std::vector< Unit,std::allocator< Unit > > const & newUnits) -> Datafield
+        """
+        return _libBornAgainDevice.Datafield_plottableField(self, *args)
 
     def flat(self):
         r"""flat(Datafield self) -> Datafield"""
diff --git a/auto/Wrap/libBornAgainDevice_wrap.cpp b/auto/Wrap/libBornAgainDevice_wrap.cpp
index ef3d413f166627be84f2fbb4de3823b51557efac..032ee49088d0ecbc86abcd421f74e84d6a5d165b 100644
--- a/auto/Wrap/libBornAgainDevice_wrap.cpp
+++ b/auto/Wrap/libBornAgainDevice_wrap.cpp
@@ -3459,24 +3459,25 @@ namespace Swig {
 #define SWIGTYPE_p_std__vectorT_INode_const_p_std__allocatorT_INode_const_p_t_t swig_types[70]
 #define SWIGTYPE_p_std__vectorT_ParaMeta_std__allocatorT_ParaMeta_t_t swig_types[71]
 #define SWIGTYPE_p_std__vectorT_Scale_const_p_std__allocatorT_Scale_const_p_t_t swig_types[72]
-#define SWIGTYPE_p_std__vectorT_Vec3T_double_t_std__allocatorT_Vec3T_double_t_t_t swig_types[73]
-#define SWIGTYPE_p_std__vectorT_double_std__allocatorT_double_t_t swig_types[74]
-#define SWIGTYPE_p_std__vectorT_int_std__allocatorT_int_t_t swig_types[75]
-#define SWIGTYPE_p_std__vectorT_size_t_std__allocatorT_size_t_t_t swig_types[76]
-#define SWIGTYPE_p_std__vectorT_std__complexT_double_t_std__allocatorT_std__complexT_double_t_t_t swig_types[77]
-#define SWIGTYPE_p_std__vectorT_std__pairT_double_double_t_std__allocatorT_std__pairT_double_double_t_t_t swig_types[78]
-#define SWIGTYPE_p_std__vectorT_std__string_std__allocatorT_std__string_t_t swig_types[79]
-#define SWIGTYPE_p_std__vectorT_std__vectorT_double_std__allocatorT_double_t_t_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t_t swig_types[80]
-#define SWIGTYPE_p_std__vectorT_std__vectorT_int_std__allocatorT_int_t_t_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t_t swig_types[81]
-#define SWIGTYPE_p_std__vectorT_unsigned_long_std__allocatorT_unsigned_long_t_t swig_types[82]
-#define SWIGTYPE_p_swig__SwigPyIterator swig_types[83]
-#define SWIGTYPE_p_unsigned_char swig_types[84]
-#define SWIGTYPE_p_unsigned_int swig_types[85]
-#define SWIGTYPE_p_unsigned_long_long swig_types[86]
-#define SWIGTYPE_p_unsigned_short swig_types[87]
-#define SWIGTYPE_p_value_type swig_types[88]
-static swig_type_info *swig_types[90];
-static swig_module_info swig_module = {swig_types, 89, 0, 0, 0, 0};
+#define SWIGTYPE_p_std__vectorT_Unit_std__allocatorT_Unit_t_t swig_types[73]
+#define SWIGTYPE_p_std__vectorT_Vec3T_double_t_std__allocatorT_Vec3T_double_t_t_t swig_types[74]
+#define SWIGTYPE_p_std__vectorT_double_std__allocatorT_double_t_t swig_types[75]
+#define SWIGTYPE_p_std__vectorT_int_std__allocatorT_int_t_t swig_types[76]
+#define SWIGTYPE_p_std__vectorT_size_t_std__allocatorT_size_t_t_t swig_types[77]
+#define SWIGTYPE_p_std__vectorT_std__complexT_double_t_std__allocatorT_std__complexT_double_t_t_t swig_types[78]
+#define SWIGTYPE_p_std__vectorT_std__pairT_double_double_t_std__allocatorT_std__pairT_double_double_t_t_t swig_types[79]
+#define SWIGTYPE_p_std__vectorT_std__string_std__allocatorT_std__string_t_t swig_types[80]
+#define SWIGTYPE_p_std__vectorT_std__vectorT_double_std__allocatorT_double_t_t_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t_t swig_types[81]
+#define SWIGTYPE_p_std__vectorT_std__vectorT_int_std__allocatorT_int_t_t_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t_t swig_types[82]
+#define SWIGTYPE_p_std__vectorT_unsigned_long_std__allocatorT_unsigned_long_t_t swig_types[83]
+#define SWIGTYPE_p_swig__SwigPyIterator swig_types[84]
+#define SWIGTYPE_p_unsigned_char swig_types[85]
+#define SWIGTYPE_p_unsigned_int swig_types[86]
+#define SWIGTYPE_p_unsigned_long_long swig_types[87]
+#define SWIGTYPE_p_unsigned_short swig_types[88]
+#define SWIGTYPE_p_value_type swig_types[89]
+static swig_type_info *swig_types[91];
+static swig_module_info swig_module = {swig_types, 90, 0, 0, 0, 0};
 #define SWIG_TypeQuery(name) SWIG_TypeQueryModule(&swig_module, &swig_module, name)
 #define SWIG_MangledTypeQuery(name) SWIG_MangledTypeQueryModule(&swig_module, &swig_module, name)
 
@@ -28837,16 +28838,14 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_Datafield_plottableField(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_Datafield_plottableField__SWIG_0(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   Datafield *arg1 = (Datafield *) 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  PyObject *swig_obj[1] ;
   SwigValueWrapper< Datafield > result;
   
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
+  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Datafield, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
     SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Datafield_plottableField" "', argument " "1"" of type '" "Datafield const *""'"); 
@@ -28860,6 +28859,78 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_Datafield_plottableField__SWIG_1(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
+  PyObject *resultobj = 0;
+  Datafield *arg1 = (Datafield *) 0 ;
+  std::vector< Unit,std::allocator< Unit > > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  SwigValueWrapper< Datafield > result;
+  
+  if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Datafield, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Datafield_plottableField" "', argument " "1"" of type '" "Datafield const *""'"); 
+  }
+  arg1 = reinterpret_cast< Datafield * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_std__vectorT_Unit_std__allocatorT_Unit_t_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "Datafield_plottableField" "', argument " "2"" of type '" "std::vector< Unit,std::allocator< Unit > > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "Datafield_plottableField" "', argument " "2"" of type '" "std::vector< Unit,std::allocator< Unit > > const &""'"); 
+  }
+  arg2 = reinterpret_cast< std::vector< Unit,std::allocator< Unit > > * >(argp2);
+  result = ((Datafield const *)arg1)->plottableField((std::vector< Unit,std::allocator< Unit > > const &)*arg2);
+  resultobj = SWIG_NewPointerObj((new Datafield(result)), SWIGTYPE_p_Datafield, SWIG_POINTER_OWN |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_Datafield_plottableField(PyObject *self, PyObject *args) {
+  Py_ssize_t argc;
+  PyObject *argv[3] = {
+    0
+  };
+  
+  if (!(argc = SWIG_Python_UnpackTuple(args, "Datafield_plottableField", 0, 2, argv))) SWIG_fail;
+  --argc;
+  if (argc == 1) {
+    int _v = 0;
+    void *vptr = 0;
+    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_Datafield, 0);
+    _v = SWIG_CheckState(res);
+    if (_v) {
+      return _wrap_Datafield_plottableField__SWIG_0(self, argc, argv);
+    }
+  }
+  if (argc == 2) {
+    int _v = 0;
+    void *vptr = 0;
+    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_Datafield, 0);
+    _v = SWIG_CheckState(res);
+    if (_v) {
+      int res = SWIG_ConvertPtr(argv[1], 0, SWIGTYPE_p_std__vectorT_Unit_std__allocatorT_Unit_t_t, SWIG_POINTER_NO_NULL | 0);
+      _v = SWIG_CheckState(res);
+      if (_v) {
+        return _wrap_Datafield_plottableField__SWIG_1(self, argc, argv);
+      }
+    }
+  }
+  
+fail:
+  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'Datafield_plottableField'.\n"
+    "  Possible C/C++ prototypes are:\n"
+    "    Datafield::plottableField() const\n"
+    "    Datafield::plottableField(std::vector< Unit,std::allocator< Unit > > const &) const\n");
+  return 0;
+}
+
+
 SWIGINTERN PyObject *_wrap_Datafield_flat(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Datafield *arg1 = (Datafield *) 0 ;
@@ -37826,7 +37897,10 @@ static PyMethodDef SwigMethods[] = {
 	 { "Datafield_flatVector", _wrap_Datafield_flatVector, METH_O, "Datafield_flatVector(Datafield self) -> vdouble1d_t"},
 	 { "Datafield_maxVal", _wrap_Datafield_maxVal, METH_O, "Datafield_maxVal(Datafield self) -> double"},
 	 { "Datafield_minVal", _wrap_Datafield_minVal, METH_O, "Datafield_minVal(Datafield self) -> double"},
-	 { "Datafield_plottableField", _wrap_Datafield_plottableField, METH_O, "Datafield_plottableField(Datafield self) -> Datafield"},
+	 { "Datafield_plottableField", _wrap_Datafield_plottableField, METH_VARARGS, "\n"
+		"Datafield_plottableField(Datafield self) -> Datafield\n"
+		"Datafield_plottableField(Datafield self, std::vector< Unit,std::allocator< Unit > > const & newUnits) -> Datafield\n"
+		""},
 	 { "Datafield_flat", _wrap_Datafield_flat, METH_O, "Datafield_flat(Datafield self) -> Datafield"},
 	 { "Datafield_scale", _wrap_Datafield_scale, METH_VARARGS, "Datafield_scale(Datafield self, double factor)"},
 	 { "Datafield_crop", _wrap_Datafield_crop, METH_VARARGS, "\n"
@@ -38294,6 +38368,7 @@ static swig_type_info _swigt__p_std__pairT_double_double_t = {"_p_std__pairT_dou
 static swig_type_info _swigt__p_std__vectorT_INode_const_p_std__allocatorT_INode_const_p_t_t = {"_p_std__vectorT_INode_const_p_std__allocatorT_INode_const_p_t_t", "std::vector< INode const *,std::allocator< INode const * > > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_ParaMeta_std__allocatorT_ParaMeta_t_t = {"_p_std__vectorT_ParaMeta_std__allocatorT_ParaMeta_t_t", "std::vector< ParaMeta,std::allocator< ParaMeta > > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_Scale_const_p_std__allocatorT_Scale_const_p_t_t = {"_p_std__vectorT_Scale_const_p_std__allocatorT_Scale_const_p_t_t", "std::vector< Scale const *,std::allocator< Scale const * > > *", 0, 0, (void*)0, 0};
+static swig_type_info _swigt__p_std__vectorT_Unit_std__allocatorT_Unit_t_t = {"_p_std__vectorT_Unit_std__allocatorT_Unit_t_t", "std::vector< Unit,std::allocator< Unit > > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_Vec3T_double_t_std__allocatorT_Vec3T_double_t_t_t = {"_p_std__vectorT_Vec3T_double_t_std__allocatorT_Vec3T_double_t_t_t", "std::vector< Vec3< double >,std::allocator< Vec3< double > > > *|std::vector< Vec3< double > > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_double_std__allocatorT_double_t_t = {"_p_std__vectorT_double_std__allocatorT_double_t_t", "std::vector< double,std::allocator< double > > *|std::vector< double > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_int_std__allocatorT_int_t_t = {"_p_std__vectorT_int_std__allocatorT_int_t_t", "std::vector< int,std::allocator< int > > *|std::vector< int > *", 0, 0, (void*)0, 0};
@@ -38385,6 +38460,7 @@ static swig_type_info *swig_type_initial[] = {
   &_swigt__p_std__vectorT_INode_const_p_std__allocatorT_INode_const_p_t_t,
   &_swigt__p_std__vectorT_ParaMeta_std__allocatorT_ParaMeta_t_t,
   &_swigt__p_std__vectorT_Scale_const_p_std__allocatorT_Scale_const_p_t_t,
+  &_swigt__p_std__vectorT_Unit_std__allocatorT_Unit_t_t,
   &_swigt__p_std__vectorT_Vec3T_double_t_std__allocatorT_Vec3T_double_t_t_t,
   &_swigt__p_std__vectorT_double_std__allocatorT_double_t_t,
   &_swigt__p_std__vectorT_int_std__allocatorT_int_t_t,
@@ -38476,6 +38552,7 @@ static swig_cast_info _swigc__p_std__pairT_double_double_t[] = {  {&_swigt__p_st
 static swig_cast_info _swigc__p_std__vectorT_INode_const_p_std__allocatorT_INode_const_p_t_t[] = {  {&_swigt__p_std__vectorT_INode_const_p_std__allocatorT_INode_const_p_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__vectorT_ParaMeta_std__allocatorT_ParaMeta_t_t[] = {  {&_swigt__p_std__vectorT_ParaMeta_std__allocatorT_ParaMeta_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__vectorT_Scale_const_p_std__allocatorT_Scale_const_p_t_t[] = {  {&_swigt__p_std__vectorT_Scale_const_p_std__allocatorT_Scale_const_p_t_t, 0, 0, 0},{0, 0, 0, 0}};
+static swig_cast_info _swigc__p_std__vectorT_Unit_std__allocatorT_Unit_t_t[] = {  {&_swigt__p_std__vectorT_Unit_std__allocatorT_Unit_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__vectorT_Vec3T_double_t_std__allocatorT_Vec3T_double_t_t_t[] = {  {&_swigt__p_std__vectorT_Vec3T_double_t_std__allocatorT_Vec3T_double_t_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__vectorT_double_std__allocatorT_double_t_t[] = {  {&_swigt__p_std__vectorT_double_std__allocatorT_double_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__vectorT_int_std__allocatorT_int_t_t[] = {  {&_swigt__p_std__vectorT_int_std__allocatorT_int_t_t, 0, 0, 0},{0, 0, 0, 0}};
@@ -38567,6 +38644,7 @@ static swig_cast_info *swig_cast_initial[] = {
   _swigc__p_std__vectorT_INode_const_p_std__allocatorT_INode_const_p_t_t,
   _swigc__p_std__vectorT_ParaMeta_std__allocatorT_ParaMeta_t_t,
   _swigc__p_std__vectorT_Scale_const_p_std__allocatorT_Scale_const_p_t_t,
+  _swigc__p_std__vectorT_Unit_std__allocatorT_Unit_t_t,
   _swigc__p_std__vectorT_Vec3T_double_t_std__allocatorT_Vec3T_double_t_t_t,
   _swigc__p_std__vectorT_double_std__allocatorT_double_t_t,
   _swigc__p_std__vectorT_int_std__allocatorT_int_t_t,
diff --git a/rawEx/fit/specular/Specular1Par.py b/rawEx/fit/specular/Specular1Par.py
index c491a89e7aa733128b6fe154dcfaa7314379b79f..a35b3babaac7a16e9ae9866679fda45f4fb6a02e 100755
--- a/rawEx/fit/specular/Specular1Par.py
+++ b/rawEx/fit/specular/Specular1Par.py
@@ -47,7 +47,8 @@ def get_sample(P):
 
 
 def get_simulation(P):
-    scan = ba.AlphaScan(ba.ListScan("", exp_x))
+    alpha_units = ba.Unit_radian
+    scan = ba.AlphaScan(ba.ListScan("", alpha_units, exp_x))
     scan.setWavelength(1.54*angstrom)
     sample = get_sample(P)
 
diff --git a/rawEx/varia/Depthprobe1.py b/rawEx/varia/Depthprobe1.py
index a372812fd9f69594a01ba921cf956652b57a96eb..323162b1942e114c733c841c0663fda0220dfaa8 100755
--- a/rawEx/varia/Depthprobe1.py
+++ b/rawEx/varia/Depthprobe1.py
@@ -27,7 +27,8 @@ def get_simulation(sample, flags):
     scan = ba.AlphaScan(n, 0*deg, 1*deg)
     scan.setWavelength(0.3*nm)
 
-    z_axis = ba.EquiDivision("z (nm)", n, -130*nm, 30*nm)
+    z_units = ba.Unit_nanometer
+    z_axis = ba.EquiDivision("z", z_units, n, -130*nm, 30*nm)
     simulation = ba.DepthprobeSimulation(scan, sample, z_axis, flags)
 
     return simulation
diff --git a/rawEx/varia/Resonator.py b/rawEx/varia/Resonator.py
index 9856b4b4191094d5084d70a36ebdd74405055115..e52ab8b2734d0b70351a51ec3463364d40ceec2e 100755
--- a/rawEx/varia/Resonator.py
+++ b/rawEx/varia/Resonator.py
@@ -27,6 +27,7 @@ d_ang = 0.01*ba.deg  # spread width for incident angle
 #  depth position span
 z_min = -100*nm
 z_max = 100*nm
+z_units = ba.Unit_nanometer
 
 
 def get_sample():
@@ -73,7 +74,7 @@ def get_simulation(sample):
     footprint = ba.FootprintSquare(0.01)
     scan.setFootprint(footprint)
 
-    z_axis = ba.EquiDivision("z (nm)", nz, z_min, z_max)
+    z_axis = ba.EquiDivision("z", z_units, nz, z_min, z_max)
     simulation = ba.DepthprobeSimulation(scan, sample, z_axis)
 
     alpha_distr = ba.DistributionGaussian(0, d_ang, 25, 3.)
diff --git a/rawEx/varia/Transmission.py b/rawEx/varia/Transmission.py
index 55f3af0cc2773f7ee4493688a0e0fd11d3ad50d1..d4ed4343b2813e0cc9b840c0202a7d34476504f4 100755
--- a/rawEx/varia/Transmission.py
+++ b/rawEx/varia/Transmission.py
@@ -26,7 +26,8 @@ def simulate():
     scan = ba.AlphaScan(1, alpha, alpha)
     scan.setWavelength(0.3*nm)
 
-    z_axis = ba.EquiDivision("z (nm)", 1, -depth, -depth)
+    z_units = ba.Unit_nanometer
+    z_axis = ba.EquiDivision("z", z_units, 1, -depth, -depth)
     simulation = ba.DepthprobeSimulation(scan, get_sample(), z_axis, 0)
 
     result = simulation.simulate().flat()
diff --git a/rawEx/varia/TransmissionVsAlpha.py b/rawEx/varia/TransmissionVsAlpha.py
index 28b3ebb58e3de29fa978ca3910c09782738dbdcc..55bed708e6a2d953ad0584069ebf48db49486c8e 100755
--- a/rawEx/varia/TransmissionVsAlpha.py
+++ b/rawEx/varia/TransmissionVsAlpha.py
@@ -26,7 +26,8 @@ def simulate(depth):
     scan = ba.AlphaScan(n, alpha_max / n, alpha_max)
     scan.setWavelength(0.3*nm)
 
-    z_axis = ba.ListScan("z (nm)", [-depth])
+    z_units = ba.Unit_nanometer
+    z_axis = ba.ListScan("z", z_units, [-depth])
     simulation = ba.DepthprobeSimulation(scan, sample, z_axis, 0)
 
     result = simulation.simulate().flat()
diff --git a/rawEx/varia/TransmissionVsDepth.py b/rawEx/varia/TransmissionVsDepth.py
index 796b1bdce0aaf6155bb8beedc553f90bc59eae51..7d583d0e1243d5646a1e95741c7ebb0f3714a403 100755
--- a/rawEx/varia/TransmissionVsDepth.py
+++ b/rawEx/varia/TransmissionVsDepth.py
@@ -24,7 +24,8 @@ def simulate(sample, alpha_i_deg):
     scan = ba.AlphaScan(1, alpha, alpha)
     scan.setWavelength(0.3*nm)
 
-    z_axis = ba.EquiDivision("z (nm)", <%= sm ? 20 : 500 %>, -130*nm, 30*nm)
+    z_units = ba.Unit_nanometer
+    z_axis = ba.EquiDivision("z", z_units, <%= sm ? 20 : 500 %>, -130*nm, 30*nm)
     simulation = ba.DepthprobeSimulation(scan, sample, z_axis, 0)
 
     result = simulation.simulate().flat()