From 0e49962eadda6acfd2661624f3b86c3cc608603a Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (o)" <j.wuttke@fz-juelich.de>
Date: Thu, 6 Apr 2023 10:20:11 +0200
Subject: [PATCH] adapt to libheinz post 1.0: unit -> unit_or_throw

---
 Base/Pixel/RectangularPixel.cpp         |  2 +-
 Device/Coord/CoordSystem2D.cpp          |  8 ++--
 Device/Detector/RectangularDetector.cpp |  8 ++--
 Device/Pol/PolFilter.cpp                |  4 +-
 Sample/Correlations/IPeakShape.cpp      | 20 ++++----
 Sample/Lattice/Lattice3D.cpp            |  2 +-
 Tests/Unit/Numeric/BisectFF.cpp         |  4 +-
 Tests/Unit/Numeric/MultiQTest.cpp       |  2 +-
 auto/Wrap/libBornAgainBase.py           | 20 +++++---
 auto/Wrap/libBornAgainBase_wrap.cpp     | 64 +++++++++++++++++++++----
 auto/Wrap/libBornAgainDevice.py         | 20 +++++---
 auto/Wrap/libBornAgainDevice_wrap.cpp   | 64 +++++++++++++++++++++----
 auto/Wrap/libBornAgainResample.py       | 20 +++++---
 auto/Wrap/libBornAgainResample_wrap.cpp | 64 +++++++++++++++++++++----
 auto/Wrap/libBornAgainSample.py         | 20 +++++---
 auto/Wrap/libBornAgainSample_wrap.cpp   | 64 +++++++++++++++++++++----
 auto/Wrap/libBornAgainSim.py            | 20 +++++---
 auto/Wrap/libBornAgainSim_wrap.cpp      | 64 +++++++++++++++++++++----
 18 files changed, 375 insertions(+), 95 deletions(-)

diff --git a/Base/Pixel/RectangularPixel.cpp b/Base/Pixel/RectangularPixel.cpp
index 73d147f90d4..986f4aa3753 100644
--- a/Base/Pixel/RectangularPixel.cpp
+++ b/Base/Pixel/RectangularPixel.cpp
@@ -65,7 +65,7 @@ double RectangularPixel::solidAngle() const
 
 R3 RectangularPixel::normalizeLength(const R3 direction, double length) const
 {
-    return direction.unit() * length;
+    return direction.unit_or_throw() * length;
 }
 
 double RectangularPixel::calculateSolidAngle() const
diff --git a/Device/Coord/CoordSystem2D.cpp b/Device/Coord/CoordSystem2D.cpp
index a499d0a9d0d..d1ceefe2fad 100644
--- a/Device/Coord/CoordSystem2D.cpp
+++ b/Device/Coord/CoordSystem2D.cpp
@@ -126,7 +126,7 @@ double SphericalCoords::calculateValue(size_t i_axis, Coords units, double value
         if (i_axis == 1) {
             // v axis is perpendicular to ki and y.
             const R3 kf = vecOfKAlphaPhi(m_ki.mag(), value);
-            static const R3 unit_v = (m_ki.cross(R3(0, 1, 0))).unit();
+            static const R3 unit_v = (m_ki.cross(R3(0, 1, 0))).unit_or_throw();
             return (kf - m_ki).dot(unit_v);
         }
         ASSERT(false);
@@ -212,8 +212,8 @@ double ImageCoords::calculateValue(size_t i_axis, Coords units, double value) co
     const auto k10 = m_detector_pixel->getPosition(1.0, 0.0);
     const auto& max_pos = i_axis == 0 ? k10 : k01; // position of max along given axis
     const double shift = value - m_axes[i_axis]->min();
-    const R3 out_dir = k00 + shift * (max_pos - k00).unit();
-    const R3 kf = out_dir.unit() * m_ki.mag();
+    const R3 out_dir = k00 + shift * (max_pos - k00).unit_or_throw();
+    const R3 kf = out_dir.unit_or_throw() * m_ki.mag();
 
     switch (units) {
     case Coords::RADIANS:
@@ -226,7 +226,7 @@ double ImageCoords::calculateValue(size_t i_axis, Coords units, double value) co
             return (m_ki - kf).y();
         if (i_axis == 1) {
             // v axis is perpendicular to ki and y.
-            static const R3 unit_v = (m_ki.cross(R3(0, 1, 0))).unit();
+            static const R3 unit_v = (m_ki.cross(R3(0, 1, 0))).unit_or_throw();
             return (kf - m_ki).dot(unit_v);
         }
     } break;
diff --git a/Device/Detector/RectangularDetector.cpp b/Device/Detector/RectangularDetector.cpp
index 551632aea18..0178f0762ba 100644
--- a/Device/Detector/RectangularDetector.cpp
+++ b/Device/Detector/RectangularDetector.cpp
@@ -183,7 +183,7 @@ const IPixel* RectangularDetector::createPixel(size_t index) const
 size_t RectangularDetector::indexOfSpecular(const Beam& beam) const
 {
     const R3 k_spec = beam.k_reflected();
-    const R3 normal_unit = m_normal_to_detector.unit();
+    const R3 normal_unit = m_normal_to_detector.unit_or_throw();
     const double kd = k_spec.dot(normal_unit);
     if (kd <= 0.0)
         return totalSize();
@@ -213,7 +213,7 @@ void RectangularDetector::setDistanceAndOffset(double distance, double u0, doubl
 
 void RectangularDetector::initNormalVector(const R3 central_k)
 {
-    R3 central_k_unit = central_k.unit();
+    R3 central_k_unit = central_k.unit_or_throw();
 
     if (m_detector_arrangement == GENERIC) {
         // do nothing
@@ -241,8 +241,8 @@ void RectangularDetector::initUandV()
     double d2 = m_normal_to_detector.dot(m_normal_to_detector);
     R3 u_direction =
         d2 * m_direction - m_direction.dot(m_normal_to_detector) * m_normal_to_detector;
-    m_u_unit = u_direction.unit();
-    m_v_unit = m_u_unit.cross(m_normal_to_detector).unit();
+    m_u_unit = u_direction.unit_or_throw();
+    m_v_unit = m_u_unit.cross(m_normal_to_detector).unit_or_throw();
 }
 
 const CoordSystem2D* RectangularDetector::scatteringCoords(const Beam& beam) const
diff --git a/Device/Pol/PolFilter.cpp b/Device/Pol/PolFilter.cpp
index e203ba99568..4e508ca55ea 100644
--- a/Device/Pol/PolFilter.cpp
+++ b/Device/Pol/PolFilter.cpp
@@ -26,7 +26,7 @@ PolFilter::PolFilter(R3 direction, double efficiency, double total_transmission)
         m_efficiency = 0;
         m_total_transmission = total_transmission;
     } else {
-        m_direction = direction.unit();
+        m_direction = direction.unit_or_throw();
         m_efficiency = efficiency;
         m_total_transmission = total_transmission;
     }
@@ -41,7 +41,7 @@ SpinMatrix PolFilter::matrix() const
 {
     if (m_direction.mag() == 0.0 || m_efficiency == 0.0)
         return m_total_transmission * SpinMatrix::One();
-    R3 v = m_direction.unit() * m_efficiency;
+    R3 v = m_direction.unit_or_throw() * m_efficiency;
     return 2 * m_total_transmission * SpinMatrix::FromBlochVector(v);
 }
 
diff --git a/Sample/Correlations/IPeakShape.cpp b/Sample/Correlations/IPeakShape.cpp
index 03597a27449..778628ab59f 100644
--- a/Sample/Correlations/IPeakShape.cpp
+++ b/Sample/Correlations/IPeakShape.cpp
@@ -211,7 +211,7 @@ MisesFisherGaussPeakShape::MisesFisherGaussPeakShape(double max_intensity, doubl
                                                      R3 zenith, double kappa_1, double kappa_2)
     : m_max_intensity(max_intensity)
     , m_radial_size(radial_size)
-    , m_zenith(zenith.unit())
+    , m_zenith(zenith.unit_or_throw())
     , m_kappa_1(kappa_1)
     , m_kappa_2(kappa_2)
 {
@@ -238,17 +238,17 @@ double MisesFisherGaussPeakShape::peakDistribution(const R3 q, const R3 q_lattic
     // angular part
     const R3 vy = m_zenith.cross(q_lattice_point);
     const R3 zxq = m_zenith.cross(q);
-    const R3 up = q_lattice_point.unit();
+    const R3 up = q_lattice_point.unit_or_throw();
     if (vy.mag2() <= 0.0 || zxq.mag2() <= 0.0) {
-        const double x = q.unit().dot(up);
+        const double x = q.unit_or_throw().dot(up);
         const double angular_part = FisherDistribution(x, m_kappa_1);
         return m_max_intensity * radial_part * angular_part;
     }
-    const R3 uy = vy.unit();
+    const R3 uy = vy.unit_or_throw();
     const R3 ux = uy.cross(m_zenith);
     const R3 q_ortho = q - q.dot(m_zenith) * m_zenith;
-    const double phi0 = std::acos(q_ortho.unit().dot(ux));
-    const double theta = std::acos(q.unit().dot(m_zenith));
+    const double phi0 = std::acos(q_ortho.unit_or_throw().dot(ux));
+    const double theta = std::acos(q.unit_or_throw().dot(m_zenith));
     const double pre_1 = FisherPrefactor(m_kappa_1);
     const double pre_2 = MisesPrefactor(m_kappa_2);
     const double integral = RealIntegrator().integrate(
@@ -271,7 +271,7 @@ MisesGaussPeakShape::MisesGaussPeakShape(double max_intensity, double radial_siz
                                          double kappa)
     : m_max_intensity(max_intensity)
     , m_radial_size(radial_size)
-    , m_zenith(zenith.unit())
+    , m_zenith(zenith.unit_or_throw())
     , m_kappa(kappa)
 {
 }
@@ -293,11 +293,11 @@ double MisesGaussPeakShape::peakDistribution(const R3 q, const R3 q_lattice_poin
     }
     const double m_qr = q.mag();
     const R3 m_p = q_lattice_point;
-    const R3 uy = vy.unit();
+    const R3 uy = vy.unit_or_throw();
     const R3 ux = uy.cross(m_zenith);
     const R3 q_ortho = q - q.dot(m_zenith) * m_zenith;
-    const double phi0 = std::acos(q_ortho.unit().dot(ux));
-    const double theta = std::acos(q.unit().dot(m_zenith));
+    const double phi0 = std::acos(q_ortho.unit_or_throw().dot(ux));
+    const double theta = std::acos(q.unit_or_throw().dot(m_zenith));
     const double pre = MisesPrefactor(m_kappa);
     const double integral = RealIntegrator().integrate(
         [&](double phi) -> double {
diff --git a/Sample/Lattice/Lattice3D.cpp b/Sample/Lattice/Lattice3D.cpp
index 50692cabc86..6d39c77bcca 100644
--- a/Sample/Lattice/Lattice3D.cpp
+++ b/Sample/Lattice/Lattice3D.cpp
@@ -49,7 +49,7 @@ Lattice3D Lattice3D::rotated(const RotMatrix& rotMatrix) const
 R3 Lattice3D::getMillerDirection(double h, double k, double l) const
 {
     R3 direction = h * m_ra + k * m_rb + l * m_rc;
-    return direction.unit();
+    return direction.unit_or_throw();
 }
 
 double Lattice3D::unitCellVolume() const
diff --git a/Tests/Unit/Numeric/BisectFF.cpp b/Tests/Unit/Numeric/BisectFF.cpp
index 8cd7da8dad8..530fd9f4cf5 100644
--- a/Tests/Unit/Numeric/BisectFF.cpp
+++ b/Tests/Unit/Numeric/BisectFF.cpp
@@ -102,8 +102,8 @@ void run_test(IFormFactor& ff)
     ::testing::internal::ParamGenerator<std::tuple<C3, C3, double, double, double>> gen = qlist;
     int ifail = 0;
     for (auto it : gen) {
-        const C3 q_dir0 = std::get<0>(it).unit();
-        const C3 q_dir1 = std::get<1>(it).unit();
+        const C3 q_dir0 = std::get<0>(it).unit_or_throw();
+        const C3 q_dir1 = std::get<1>(it).unit_or_throw();
         const double qrealmag = std::get<2>(it);
         const double qrel1 = std::get<3>(it);
         const double qimagrel = std::get<4>(it);
diff --git a/Tests/Unit/Numeric/MultiQTest.cpp b/Tests/Unit/Numeric/MultiQTest.cpp
index 3bf7fef0ac7..3f97a7e8a89 100644
--- a/Tests/Unit/Numeric/MultiQTest.cpp
+++ b/Tests/Unit/Numeric/MultiQTest.cpp
@@ -47,7 +47,7 @@ int run_test_for_many_q(std::function<complex_t(C3)> fff0, std::function<complex
             continue;
         if (q_maindir == q_sidedir)
             continue;
-        const C3 q = qmag * (q_maindir + qsidemag * q_sidedir).unit();
+        const C3 q = qmag * (q_maindir + qsidemag * q_sidedir).unit_or_throw();
 
         const auto [f0, f1, avge, abserr, deviation] = evaluate(q);
 
diff --git a/auto/Wrap/libBornAgainBase.py b/auto/Wrap/libBornAgainBase.py
index a00485c698a..8050ea3c773 100644
--- a/auto/Wrap/libBornAgainBase.py
+++ b/auto/Wrap/libBornAgainBase.py
@@ -2257,9 +2257,13 @@ class R3(object):
         r"""magxy(R3 self) -> double"""
         return _libBornAgainBase.R3_magxy(self)
 
-    def unit(self):
-        r"""unit(R3 self) -> R3"""
-        return _libBornAgainBase.R3_unit(self)
+    def unit_or_throw(self):
+        r"""unit_or_throw(R3 self) -> R3"""
+        return _libBornAgainBase.R3_unit_or_throw(self)
+
+    def unit_or_null(self):
+        r"""unit_or_null(R3 self) -> R3"""
+        return _libBornAgainBase.R3_unit_or_null(self)
 
     def complex(self):
         r"""complex(R3 self) -> C3"""
@@ -2353,9 +2357,13 @@ class C3(object):
         r"""magxy(C3 self) -> double"""
         return _libBornAgainBase.C3_magxy(self)
 
-    def unit(self):
-        r"""unit(C3 self) -> C3"""
-        return _libBornAgainBase.C3_unit(self)
+    def unit_or_throw(self):
+        r"""unit_or_throw(C3 self) -> C3"""
+        return _libBornAgainBase.C3_unit_or_throw(self)
+
+    def unit_or_null(self):
+        r"""unit_or_null(C3 self) -> C3"""
+        return _libBornAgainBase.C3_unit_or_null(self)
 
     def complex(self):
         r"""complex(C3 self) -> C3"""
diff --git a/auto/Wrap/libBornAgainBase_wrap.cpp b/auto/Wrap/libBornAgainBase_wrap.cpp
index e2aad77e7d6..01f380567d8 100644
--- a/auto/Wrap/libBornAgainBase_wrap.cpp
+++ b/auto/Wrap/libBornAgainBase_wrap.cpp
@@ -28190,7 +28190,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_R3_unit(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_R3_unit_or_throw(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< double > *arg1 = (Vec3< double > *) 0 ;
   void *argp1 = 0 ;
@@ -28202,10 +28202,33 @@ SWIGINTERN PyObject *_wrap_R3_unit(PyObject *self, PyObject *args) {
   swig_obj[0] = args;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3_unit" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3_unit_or_throw" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
   }
   arg1 = reinterpret_cast< Vec3< double > * >(argp1);
-  result = ((Vec3< double > const *)arg1)->unit();
+  result = ((Vec3< double > const *)arg1)->unit_or_throw();
+  resultobj = SWIG_NewPointerObj((new Vec3< double >(result)), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_OWN |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_R3_unit_or_null(PyObject *self, PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< double > *arg1 = (Vec3< double > *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  Vec3< double > result;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3_unit_or_null" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< double > * >(argp1);
+  result = ((Vec3< double > const *)arg1)->unit_or_null();
   resultobj = SWIG_NewPointerObj((new Vec3< double >(result)), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_OWN |  0 );
   return resultobj;
 fail:
@@ -28858,7 +28881,30 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_C3_unit(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_C3_unit_or_throw(PyObject *self, PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  Vec3< std::complex< double > > result;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3_unit_or_throw" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
+  result = ((Vec3< std::complex< double > > const *)arg1)->unit_or_throw();
+  resultobj = SWIG_NewPointerObj((new Vec3< std::complex< double > >(result)), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_OWN |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_C3_unit_or_null(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
   void *argp1 = 0 ;
@@ -28870,10 +28916,10 @@ SWIGINTERN PyObject *_wrap_C3_unit(PyObject *self, PyObject *args) {
   swig_obj[0] = args;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3_unit" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3_unit_or_null" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
   }
   arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
-  result = ((Vec3< std::complex< double > > const *)arg1)->unit();
+  result = ((Vec3< std::complex< double > > const *)arg1)->unit_or_null();
   resultobj = SWIG_NewPointerObj((new Vec3< std::complex< double > >(result)), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_OWN |  0 );
   return resultobj;
 fail:
@@ -29830,7 +29876,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "R3_mag", _wrap_R3_mag, METH_O, "R3_mag(R3 self) -> double"},
 	 { "R3_magxy2", _wrap_R3_magxy2, METH_O, "R3_magxy2(R3 self) -> double"},
 	 { "R3_magxy", _wrap_R3_magxy, METH_O, "R3_magxy(R3 self) -> double"},
-	 { "R3_unit", _wrap_R3_unit, METH_O, "R3_unit(R3 self) -> R3"},
+	 { "R3_unit_or_throw", _wrap_R3_unit_or_throw, METH_O, "R3_unit_or_throw(R3 self) -> R3"},
+	 { "R3_unit_or_null", _wrap_R3_unit_or_null, METH_O, "R3_unit_or_null(R3 self) -> R3"},
 	 { "R3_complex", _wrap_R3_complex, METH_O, "R3_complex(R3 self) -> C3"},
 	 { "R3_real", _wrap_R3_real, METH_O, "R3_real(R3 self) -> R3"},
 	 { "R3___eq__", _wrap_R3___eq__, METH_VARARGS, "R3___eq__(R3 self, R3 other) -> bool"},
@@ -29857,7 +29904,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "C3_mag", _wrap_C3_mag, METH_O, "C3_mag(C3 self) -> double"},
 	 { "C3_magxy2", _wrap_C3_magxy2, METH_O, "C3_magxy2(C3 self) -> double"},
 	 { "C3_magxy", _wrap_C3_magxy, METH_O, "C3_magxy(C3 self) -> double"},
-	 { "C3_unit", _wrap_C3_unit, METH_O, "C3_unit(C3 self) -> C3"},
+	 { "C3_unit_or_throw", _wrap_C3_unit_or_throw, METH_O, "C3_unit_or_throw(C3 self) -> C3"},
+	 { "C3_unit_or_null", _wrap_C3_unit_or_null, METH_O, "C3_unit_or_null(C3 self) -> C3"},
 	 { "C3_complex", _wrap_C3_complex, METH_O, "C3_complex(C3 self) -> C3"},
 	 { "C3_real", _wrap_C3_real, METH_O, "C3_real(C3 self) -> R3"},
 	 { "C3___eq__", _wrap_C3___eq__, METH_VARARGS, "C3___eq__(C3 self, C3 other) -> bool"},
diff --git a/auto/Wrap/libBornAgainDevice.py b/auto/Wrap/libBornAgainDevice.py
index ce5883910b5..a64ff77216b 100644
--- a/auto/Wrap/libBornAgainDevice.py
+++ b/auto/Wrap/libBornAgainDevice.py
@@ -1731,9 +1731,13 @@ class R3(object):
         r"""magxy(R3 self) -> double"""
         return _libBornAgainDevice.R3_magxy(self)
 
-    def unit(self):
-        r"""unit(R3 self) -> R3"""
-        return _libBornAgainDevice.R3_unit(self)
+    def unit_or_throw(self):
+        r"""unit_or_throw(R3 self) -> R3"""
+        return _libBornAgainDevice.R3_unit_or_throw(self)
+
+    def unit_or_null(self):
+        r"""unit_or_null(R3 self) -> R3"""
+        return _libBornAgainDevice.R3_unit_or_null(self)
 
     def complex(self):
         r"""complex(R3 self) -> C3"""
@@ -1827,9 +1831,13 @@ class C3(object):
         r"""magxy(C3 self) -> double"""
         return _libBornAgainDevice.C3_magxy(self)
 
-    def unit(self):
-        r"""unit(C3 self) -> C3"""
-        return _libBornAgainDevice.C3_unit(self)
+    def unit_or_throw(self):
+        r"""unit_or_throw(C3 self) -> C3"""
+        return _libBornAgainDevice.C3_unit_or_throw(self)
+
+    def unit_or_null(self):
+        r"""unit_or_null(C3 self) -> C3"""
+        return _libBornAgainDevice.C3_unit_or_null(self)
 
     def complex(self):
         r"""complex(C3 self) -> C3"""
diff --git a/auto/Wrap/libBornAgainDevice_wrap.cpp b/auto/Wrap/libBornAgainDevice_wrap.cpp
index 98e10fcd5c6..7954367e655 100644
--- a/auto/Wrap/libBornAgainDevice_wrap.cpp
+++ b/auto/Wrap/libBornAgainDevice_wrap.cpp
@@ -24590,7 +24590,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_R3_unit(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_R3_unit_or_throw(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< double > *arg1 = (Vec3< double > *) 0 ;
   void *argp1 = 0 ;
@@ -24602,10 +24602,33 @@ SWIGINTERN PyObject *_wrap_R3_unit(PyObject *self, PyObject *args) {
   swig_obj[0] = args;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3_unit" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3_unit_or_throw" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
   }
   arg1 = reinterpret_cast< Vec3< double > * >(argp1);
-  result = ((Vec3< double > const *)arg1)->unit();
+  result = ((Vec3< double > const *)arg1)->unit_or_throw();
+  resultobj = SWIG_NewPointerObj((new Vec3< double >(result)), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_OWN |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_R3_unit_or_null(PyObject *self, PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< double > *arg1 = (Vec3< double > *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  Vec3< double > result;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3_unit_or_null" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< double > * >(argp1);
+  result = ((Vec3< double > const *)arg1)->unit_or_null();
   resultobj = SWIG_NewPointerObj((new Vec3< double >(result)), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_OWN |  0 );
   return resultobj;
 fail:
@@ -25258,7 +25281,30 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_C3_unit(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_C3_unit_or_throw(PyObject *self, PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  Vec3< std::complex< double > > result;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3_unit_or_throw" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
+  result = ((Vec3< std::complex< double > > const *)arg1)->unit_or_throw();
+  resultobj = SWIG_NewPointerObj((new Vec3< std::complex< double > >(result)), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_OWN |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_C3_unit_or_null(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
   void *argp1 = 0 ;
@@ -25270,10 +25316,10 @@ SWIGINTERN PyObject *_wrap_C3_unit(PyObject *self, PyObject *args) {
   swig_obj[0] = args;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3_unit" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3_unit_or_null" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
   }
   arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
-  result = ((Vec3< std::complex< double > > const *)arg1)->unit();
+  result = ((Vec3< std::complex< double > > const *)arg1)->unit_or_null();
   resultobj = SWIG_NewPointerObj((new Vec3< std::complex< double > >(result)), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_OWN |  0 );
   return resultobj;
 fail:
@@ -38164,7 +38210,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "R3_mag", _wrap_R3_mag, METH_O, "R3_mag(R3 self) -> double"},
 	 { "R3_magxy2", _wrap_R3_magxy2, METH_O, "R3_magxy2(R3 self) -> double"},
 	 { "R3_magxy", _wrap_R3_magxy, METH_O, "R3_magxy(R3 self) -> double"},
-	 { "R3_unit", _wrap_R3_unit, METH_O, "R3_unit(R3 self) -> R3"},
+	 { "R3_unit_or_throw", _wrap_R3_unit_or_throw, METH_O, "R3_unit_or_throw(R3 self) -> R3"},
+	 { "R3_unit_or_null", _wrap_R3_unit_or_null, METH_O, "R3_unit_or_null(R3 self) -> R3"},
 	 { "R3_complex", _wrap_R3_complex, METH_O, "R3_complex(R3 self) -> C3"},
 	 { "R3_real", _wrap_R3_real, METH_O, "R3_real(R3 self) -> R3"},
 	 { "R3___eq__", _wrap_R3___eq__, METH_VARARGS, "R3___eq__(R3 self, R3 other) -> bool"},
@@ -38191,7 +38238,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "C3_mag", _wrap_C3_mag, METH_O, "C3_mag(C3 self) -> double"},
 	 { "C3_magxy2", _wrap_C3_magxy2, METH_O, "C3_magxy2(C3 self) -> double"},
 	 { "C3_magxy", _wrap_C3_magxy, METH_O, "C3_magxy(C3 self) -> double"},
-	 { "C3_unit", _wrap_C3_unit, METH_O, "C3_unit(C3 self) -> C3"},
+	 { "C3_unit_or_throw", _wrap_C3_unit_or_throw, METH_O, "C3_unit_or_throw(C3 self) -> C3"},
+	 { "C3_unit_or_null", _wrap_C3_unit_or_null, METH_O, "C3_unit_or_null(C3 self) -> C3"},
 	 { "C3_complex", _wrap_C3_complex, METH_O, "C3_complex(C3 self) -> C3"},
 	 { "C3_real", _wrap_C3_real, METH_O, "C3_real(C3 self) -> R3"},
 	 { "C3___eq__", _wrap_C3___eq__, METH_VARARGS, "C3___eq__(C3 self, C3 other) -> bool"},
diff --git a/auto/Wrap/libBornAgainResample.py b/auto/Wrap/libBornAgainResample.py
index 36576948488..7c69e8c52c0 100644
--- a/auto/Wrap/libBornAgainResample.py
+++ b/auto/Wrap/libBornAgainResample.py
@@ -1763,9 +1763,13 @@ class R3(object):
         r"""magxy(R3 self) -> double"""
         return _libBornAgainResample.R3_magxy(self)
 
-    def unit(self):
-        r"""unit(R3 self) -> R3"""
-        return _libBornAgainResample.R3_unit(self)
+    def unit_or_throw(self):
+        r"""unit_or_throw(R3 self) -> R3"""
+        return _libBornAgainResample.R3_unit_or_throw(self)
+
+    def unit_or_null(self):
+        r"""unit_or_null(R3 self) -> R3"""
+        return _libBornAgainResample.R3_unit_or_null(self)
 
     def complex(self):
         r"""complex(R3 self) -> C3"""
@@ -1859,9 +1863,13 @@ class C3(object):
         r"""magxy(C3 self) -> double"""
         return _libBornAgainResample.C3_magxy(self)
 
-    def unit(self):
-        r"""unit(C3 self) -> C3"""
-        return _libBornAgainResample.C3_unit(self)
+    def unit_or_throw(self):
+        r"""unit_or_throw(C3 self) -> C3"""
+        return _libBornAgainResample.C3_unit_or_throw(self)
+
+    def unit_or_null(self):
+        r"""unit_or_null(C3 self) -> C3"""
+        return _libBornAgainResample.C3_unit_or_null(self)
 
     def complex(self):
         r"""complex(C3 self) -> C3"""
diff --git a/auto/Wrap/libBornAgainResample_wrap.cpp b/auto/Wrap/libBornAgainResample_wrap.cpp
index 4ae506abf18..1e369af586b 100644
--- a/auto/Wrap/libBornAgainResample_wrap.cpp
+++ b/auto/Wrap/libBornAgainResample_wrap.cpp
@@ -24839,7 +24839,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_R3_unit(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_R3_unit_or_throw(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< double > *arg1 = (Vec3< double > *) 0 ;
   void *argp1 = 0 ;
@@ -24851,10 +24851,33 @@ SWIGINTERN PyObject *_wrap_R3_unit(PyObject *self, PyObject *args) {
   swig_obj[0] = args;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3_unit" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3_unit_or_throw" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
   }
   arg1 = reinterpret_cast< Vec3< double > * >(argp1);
-  result = ((Vec3< double > const *)arg1)->unit();
+  result = ((Vec3< double > const *)arg1)->unit_or_throw();
+  resultobj = SWIG_NewPointerObj((new Vec3< double >(result)), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_OWN |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_R3_unit_or_null(PyObject *self, PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< double > *arg1 = (Vec3< double > *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  Vec3< double > result;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3_unit_or_null" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< double > * >(argp1);
+  result = ((Vec3< double > const *)arg1)->unit_or_null();
   resultobj = SWIG_NewPointerObj((new Vec3< double >(result)), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_OWN |  0 );
   return resultobj;
 fail:
@@ -25507,7 +25530,30 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_C3_unit(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_C3_unit_or_throw(PyObject *self, PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  Vec3< std::complex< double > > result;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3_unit_or_throw" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
+  result = ((Vec3< std::complex< double > > const *)arg1)->unit_or_throw();
+  resultobj = SWIG_NewPointerObj((new Vec3< std::complex< double > >(result)), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_OWN |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_C3_unit_or_null(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
   void *argp1 = 0 ;
@@ -25519,10 +25565,10 @@ SWIGINTERN PyObject *_wrap_C3_unit(PyObject *self, PyObject *args) {
   swig_obj[0] = args;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3_unit" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3_unit_or_null" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
   }
   arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
-  result = ((Vec3< std::complex< double > > const *)arg1)->unit();
+  result = ((Vec3< std::complex< double > > const *)arg1)->unit_or_null();
   resultobj = SWIG_NewPointerObj((new Vec3< std::complex< double > >(result)), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_OWN |  0 );
   return resultobj;
 fail:
@@ -28274,7 +28320,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "R3_mag", _wrap_R3_mag, METH_O, "R3_mag(R3 self) -> double"},
 	 { "R3_magxy2", _wrap_R3_magxy2, METH_O, "R3_magxy2(R3 self) -> double"},
 	 { "R3_magxy", _wrap_R3_magxy, METH_O, "R3_magxy(R3 self) -> double"},
-	 { "R3_unit", _wrap_R3_unit, METH_O, "R3_unit(R3 self) -> R3"},
+	 { "R3_unit_or_throw", _wrap_R3_unit_or_throw, METH_O, "R3_unit_or_throw(R3 self) -> R3"},
+	 { "R3_unit_or_null", _wrap_R3_unit_or_null, METH_O, "R3_unit_or_null(R3 self) -> R3"},
 	 { "R3_complex", _wrap_R3_complex, METH_O, "R3_complex(R3 self) -> C3"},
 	 { "R3_real", _wrap_R3_real, METH_O, "R3_real(R3 self) -> R3"},
 	 { "R3___eq__", _wrap_R3___eq__, METH_VARARGS, "R3___eq__(R3 self, R3 other) -> bool"},
@@ -28301,7 +28348,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "C3_mag", _wrap_C3_mag, METH_O, "C3_mag(C3 self) -> double"},
 	 { "C3_magxy2", _wrap_C3_magxy2, METH_O, "C3_magxy2(C3 self) -> double"},
 	 { "C3_magxy", _wrap_C3_magxy, METH_O, "C3_magxy(C3 self) -> double"},
-	 { "C3_unit", _wrap_C3_unit, METH_O, "C3_unit(C3 self) -> C3"},
+	 { "C3_unit_or_throw", _wrap_C3_unit_or_throw, METH_O, "C3_unit_or_throw(C3 self) -> C3"},
+	 { "C3_unit_or_null", _wrap_C3_unit_or_null, METH_O, "C3_unit_or_null(C3 self) -> C3"},
 	 { "C3_complex", _wrap_C3_complex, METH_O, "C3_complex(C3 self) -> C3"},
 	 { "C3_real", _wrap_C3_real, METH_O, "C3_real(C3 self) -> R3"},
 	 { "C3___eq__", _wrap_C3___eq__, METH_VARARGS, "C3___eq__(C3 self, C3 other) -> bool"},
diff --git a/auto/Wrap/libBornAgainSample.py b/auto/Wrap/libBornAgainSample.py
index f707827c767..eb89faa3fcc 100644
--- a/auto/Wrap/libBornAgainSample.py
+++ b/auto/Wrap/libBornAgainSample.py
@@ -1730,9 +1730,13 @@ class R3(object):
         r"""magxy(R3 self) -> double"""
         return _libBornAgainSample.R3_magxy(self)
 
-    def unit(self):
-        r"""unit(R3 self) -> R3"""
-        return _libBornAgainSample.R3_unit(self)
+    def unit_or_throw(self):
+        r"""unit_or_throw(R3 self) -> R3"""
+        return _libBornAgainSample.R3_unit_or_throw(self)
+
+    def unit_or_null(self):
+        r"""unit_or_null(R3 self) -> R3"""
+        return _libBornAgainSample.R3_unit_or_null(self)
 
     def complex(self):
         r"""complex(R3 self) -> C3"""
@@ -1826,9 +1830,13 @@ class C3(object):
         r"""magxy(C3 self) -> double"""
         return _libBornAgainSample.C3_magxy(self)
 
-    def unit(self):
-        r"""unit(C3 self) -> C3"""
-        return _libBornAgainSample.C3_unit(self)
+    def unit_or_throw(self):
+        r"""unit_or_throw(C3 self) -> C3"""
+        return _libBornAgainSample.C3_unit_or_throw(self)
+
+    def unit_or_null(self):
+        r"""unit_or_null(C3 self) -> C3"""
+        return _libBornAgainSample.C3_unit_or_null(self)
 
     def complex(self):
         r"""complex(C3 self) -> C3"""
diff --git a/auto/Wrap/libBornAgainSample_wrap.cpp b/auto/Wrap/libBornAgainSample_wrap.cpp
index 7f1fa602c97..ac6755c4b6c 100644
--- a/auto/Wrap/libBornAgainSample_wrap.cpp
+++ b/auto/Wrap/libBornAgainSample_wrap.cpp
@@ -25369,7 +25369,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_R3_unit(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_R3_unit_or_throw(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< double > *arg1 = (Vec3< double > *) 0 ;
   void *argp1 = 0 ;
@@ -25381,10 +25381,33 @@ SWIGINTERN PyObject *_wrap_R3_unit(PyObject *self, PyObject *args) {
   swig_obj[0] = args;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3_unit" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3_unit_or_throw" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
   }
   arg1 = reinterpret_cast< Vec3< double > * >(argp1);
-  result = ((Vec3< double > const *)arg1)->unit();
+  result = ((Vec3< double > const *)arg1)->unit_or_throw();
+  resultobj = SWIG_NewPointerObj((new Vec3< double >(result)), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_OWN |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_R3_unit_or_null(PyObject *self, PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< double > *arg1 = (Vec3< double > *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  Vec3< double > result;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3_unit_or_null" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< double > * >(argp1);
+  result = ((Vec3< double > const *)arg1)->unit_or_null();
   resultobj = SWIG_NewPointerObj((new Vec3< double >(result)), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_OWN |  0 );
   return resultobj;
 fail:
@@ -26037,7 +26060,30 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_C3_unit(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_C3_unit_or_throw(PyObject *self, PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  Vec3< std::complex< double > > result;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3_unit_or_throw" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
+  result = ((Vec3< std::complex< double > > const *)arg1)->unit_or_throw();
+  resultobj = SWIG_NewPointerObj((new Vec3< std::complex< double > >(result)), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_OWN |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_C3_unit_or_null(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
   void *argp1 = 0 ;
@@ -26049,10 +26095,10 @@ SWIGINTERN PyObject *_wrap_C3_unit(PyObject *self, PyObject *args) {
   swig_obj[0] = args;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3_unit" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3_unit_or_null" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
   }
   arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
-  result = ((Vec3< std::complex< double > > const *)arg1)->unit();
+  result = ((Vec3< std::complex< double > > const *)arg1)->unit_or_null();
   resultobj = SWIG_NewPointerObj((new Vec3< std::complex< double > >(result)), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_OWN |  0 );
   return resultobj;
 fail:
@@ -57641,7 +57687,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "R3_mag", _wrap_R3_mag, METH_O, "R3_mag(R3 self) -> double"},
 	 { "R3_magxy2", _wrap_R3_magxy2, METH_O, "R3_magxy2(R3 self) -> double"},
 	 { "R3_magxy", _wrap_R3_magxy, METH_O, "R3_magxy(R3 self) -> double"},
-	 { "R3_unit", _wrap_R3_unit, METH_O, "R3_unit(R3 self) -> R3"},
+	 { "R3_unit_or_throw", _wrap_R3_unit_or_throw, METH_O, "R3_unit_or_throw(R3 self) -> R3"},
+	 { "R3_unit_or_null", _wrap_R3_unit_or_null, METH_O, "R3_unit_or_null(R3 self) -> R3"},
 	 { "R3_complex", _wrap_R3_complex, METH_O, "R3_complex(R3 self) -> C3"},
 	 { "R3_real", _wrap_R3_real, METH_O, "R3_real(R3 self) -> R3"},
 	 { "R3___eq__", _wrap_R3___eq__, METH_VARARGS, "R3___eq__(R3 self, R3 other) -> bool"},
@@ -57668,7 +57715,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "C3_mag", _wrap_C3_mag, METH_O, "C3_mag(C3 self) -> double"},
 	 { "C3_magxy2", _wrap_C3_magxy2, METH_O, "C3_magxy2(C3 self) -> double"},
 	 { "C3_magxy", _wrap_C3_magxy, METH_O, "C3_magxy(C3 self) -> double"},
-	 { "C3_unit", _wrap_C3_unit, METH_O, "C3_unit(C3 self) -> C3"},
+	 { "C3_unit_or_throw", _wrap_C3_unit_or_throw, METH_O, "C3_unit_or_throw(C3 self) -> C3"},
+	 { "C3_unit_or_null", _wrap_C3_unit_or_null, METH_O, "C3_unit_or_null(C3 self) -> C3"},
 	 { "C3_complex", _wrap_C3_complex, METH_O, "C3_complex(C3 self) -> C3"},
 	 { "C3_real", _wrap_C3_real, METH_O, "C3_real(C3 self) -> R3"},
 	 { "C3___eq__", _wrap_C3___eq__, METH_VARARGS, "C3___eq__(C3 self, C3 other) -> bool"},
diff --git a/auto/Wrap/libBornAgainSim.py b/auto/Wrap/libBornAgainSim.py
index 54c85365fb2..e76bfc952aa 100644
--- a/auto/Wrap/libBornAgainSim.py
+++ b/auto/Wrap/libBornAgainSim.py
@@ -1731,9 +1731,13 @@ class R3(object):
         r"""magxy(R3 self) -> double"""
         return _libBornAgainSim.R3_magxy(self)
 
-    def unit(self):
-        r"""unit(R3 self) -> R3"""
-        return _libBornAgainSim.R3_unit(self)
+    def unit_or_throw(self):
+        r"""unit_or_throw(R3 self) -> R3"""
+        return _libBornAgainSim.R3_unit_or_throw(self)
+
+    def unit_or_null(self):
+        r"""unit_or_null(R3 self) -> R3"""
+        return _libBornAgainSim.R3_unit_or_null(self)
 
     def complex(self):
         r"""complex(R3 self) -> C3"""
@@ -1843,9 +1847,13 @@ class C3(object):
         r"""magxy(C3 self) -> double"""
         return _libBornAgainSim.C3_magxy(self)
 
-    def unit(self):
-        r"""unit(C3 self) -> C3"""
-        return _libBornAgainSim.C3_unit(self)
+    def unit_or_throw(self):
+        r"""unit_or_throw(C3 self) -> C3"""
+        return _libBornAgainSim.C3_unit_or_throw(self)
+
+    def unit_or_null(self):
+        r"""unit_or_null(C3 self) -> C3"""
+        return _libBornAgainSim.C3_unit_or_null(self)
 
     def complex(self):
         r"""complex(C3 self) -> C3"""
diff --git a/auto/Wrap/libBornAgainSim_wrap.cpp b/auto/Wrap/libBornAgainSim_wrap.cpp
index b00b37eca66..548c4a71fd4 100644
--- a/auto/Wrap/libBornAgainSim_wrap.cpp
+++ b/auto/Wrap/libBornAgainSim_wrap.cpp
@@ -24892,7 +24892,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_R3_unit(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_R3_unit_or_throw(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< double > *arg1 = (Vec3< double > *) 0 ;
   void *argp1 = 0 ;
@@ -24904,10 +24904,33 @@ SWIGINTERN PyObject *_wrap_R3_unit(PyObject *self, PyObject *args) {
   swig_obj[0] = args;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3_unit" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3_unit_or_throw" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
   }
   arg1 = reinterpret_cast< Vec3< double > * >(argp1);
-  result = ((Vec3< double > const *)arg1)->unit();
+  result = ((Vec3< double > const *)arg1)->unit_or_throw();
+  resultobj = SWIG_NewPointerObj((new Vec3< double >(result)), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_OWN |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_R3_unit_or_null(PyObject *self, PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< double > *arg1 = (Vec3< double > *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  Vec3< double > result;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3_unit_or_null" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< double > * >(argp1);
+  result = ((Vec3< double > const *)arg1)->unit_or_null();
   resultobj = SWIG_NewPointerObj((new Vec3< double >(result)), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_OWN |  0 );
   return resultobj;
 fail:
@@ -25684,7 +25707,30 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_C3_unit(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_C3_unit_or_throw(PyObject *self, PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  Vec3< std::complex< double > > result;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3_unit_or_throw" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
+  result = ((Vec3< std::complex< double > > const *)arg1)->unit_or_throw();
+  resultobj = SWIG_NewPointerObj((new Vec3< std::complex< double > >(result)), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_OWN |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_C3_unit_or_null(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
   void *argp1 = 0 ;
@@ -25696,10 +25742,10 @@ SWIGINTERN PyObject *_wrap_C3_unit(PyObject *self, PyObject *args) {
   swig_obj[0] = args;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3_unit" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3_unit_or_null" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
   }
   arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
-  result = ((Vec3< std::complex< double > > const *)arg1)->unit();
+  result = ((Vec3< std::complex< double > > const *)arg1)->unit_or_null();
   resultobj = SWIG_NewPointerObj((new Vec3< std::complex< double > >(result)), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_OWN |  0 );
   return resultobj;
 fail:
@@ -36043,7 +36089,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "R3_mag", _wrap_R3_mag, METH_O, "R3_mag(R3 self) -> double"},
 	 { "R3_magxy2", _wrap_R3_magxy2, METH_O, "R3_magxy2(R3 self) -> double"},
 	 { "R3_magxy", _wrap_R3_magxy, METH_O, "R3_magxy(R3 self) -> double"},
-	 { "R3_unit", _wrap_R3_unit, METH_O, "R3_unit(R3 self) -> R3"},
+	 { "R3_unit_or_throw", _wrap_R3_unit_or_throw, METH_O, "R3_unit_or_throw(R3 self) -> R3"},
+	 { "R3_unit_or_null", _wrap_R3_unit_or_null, METH_O, "R3_unit_or_null(R3 self) -> R3"},
 	 { "R3_complex", _wrap_R3_complex, METH_O, "R3_complex(R3 self) -> C3"},
 	 { "R3_real", _wrap_R3_real, METH_O, "R3_real(R3 self) -> R3"},
 	 { "R3___eq__", _wrap_R3___eq__, METH_VARARGS, "R3___eq__(R3 self, R3 other) -> bool"},
@@ -36074,7 +36121,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "C3_mag", _wrap_C3_mag, METH_O, "C3_mag(C3 self) -> double"},
 	 { "C3_magxy2", _wrap_C3_magxy2, METH_O, "C3_magxy2(C3 self) -> double"},
 	 { "C3_magxy", _wrap_C3_magxy, METH_O, "C3_magxy(C3 self) -> double"},
-	 { "C3_unit", _wrap_C3_unit, METH_O, "C3_unit(C3 self) -> C3"},
+	 { "C3_unit_or_throw", _wrap_C3_unit_or_throw, METH_O, "C3_unit_or_throw(C3 self) -> C3"},
+	 { "C3_unit_or_null", _wrap_C3_unit_or_null, METH_O, "C3_unit_or_null(C3 self) -> C3"},
 	 { "C3_complex", _wrap_C3_complex, METH_O, "C3_complex(C3 self) -> C3"},
 	 { "C3_real", _wrap_C3_real, METH_O, "C3_real(C3 self) -> R3"},
 	 { "C3___eq__", _wrap_C3___eq__, METH_VARARGS, "C3___eq__(C3 self, C3 other) -> bool"},
-- 
GitLab