diff --git a/Core/Algorithms/inc/FTDistributions.h b/Core/Algorithms/inc/FTDistributions.h
index f2193a3984e9dbd83061e7619e216bbea7a67f39..4f98bb8b7c8438614d04a5cd3299046fdc0efdc0 100644
--- a/Core/Algorithms/inc/FTDistributions.h
+++ b/Core/Algorithms/inc/FTDistributions.h
@@ -18,8 +18,7 @@
 #define FTDISTRIBUTIONS_H_
 
 #include "IParameterized.h"
-#include <cmath>
-
+#include "Units.h"
 #include "Numeric.h"
 
 
@@ -168,7 +167,7 @@ public:
         : m_coherence_length_x(coherence_length_x)
         , m_coherence_length_y(coherence_length_y)
         , m_gamma(0.0)
-        , m_delta(M_PI/2.0) {}
+        , m_delta(Units::PI/2.0) {}
     virtual ~IFTDistribution2D() {}
 
     virtual IFTDistribution2D *clone() const=0;
diff --git a/Core/Algorithms/src/DWBASimulation.cpp b/Core/Algorithms/src/DWBASimulation.cpp
index 7888ddbaa76e75c4574b2e249f4e1d44a004cc3a..16c9528eed4a071bce4329db6bbcb90bb57234c2 100644
--- a/Core/Algorithms/src/DWBASimulation.cpp
+++ b/Core/Algorithms/src/DWBASimulation.cpp
@@ -108,7 +108,7 @@ bool DWBASimulation::checkPolarizationPresent() const
 double DWBASimulation::getWaveLength() const
 {
     kvector_t real_ki(m_ki.x().real(), m_ki.y().real(), m_ki.z().real());
-    return 2.0*M_PI/real_ki.mag();
+    return Units::PI2/real_ki.mag();
 }
 
 const OutputData<double>& DWBASimulation::getPolarizationData() const
diff --git a/Core/Algorithms/src/Distributions.cpp b/Core/Algorithms/src/Distributions.cpp
index 9613cc079b9f23d5e0142c0dd35ce3c1ce0906d5..4d9f14b8e219a3920fdc13459085485968375d55 100644
--- a/Core/Algorithms/src/Distributions.cpp
+++ b/Core/Algorithms/src/Distributions.cpp
@@ -173,7 +173,7 @@ DistributionLorentz::DistributionLorentz(double mean, double hwhm)
 double DistributionLorentz::probabilityDensity(double x) const
 {
     if (m_hwhm == 0.0) return x==m_mean ? 1.0 : 0.0;
-    return m_hwhm/(m_hwhm*m_hwhm + (x-m_mean)*(x-m_mean))/M_PI;
+    return m_hwhm/(m_hwhm*m_hwhm + (x-m_mean)*(x-m_mean))/Units::PI;
 }
 
 std::vector<double> DistributionLorentz::generateValueList(size_t nbr_samples,
@@ -226,7 +226,7 @@ double DistributionGaussian::probabilityDensity(double x) const
     if (m_std_dev == 0.0) return x==m_mean ? 1.0 : 0.0;
     double exponential = std::exp(-(x-m_mean)*(x-m_mean)
             /(2.0*m_std_dev*m_std_dev));
-    return exponential/m_std_dev/std::sqrt(2.0*M_PI);
+    return exponential/m_std_dev/std::sqrt(Units::PI2);
 }
 
 std::vector<double> DistributionGaussian::generateValueList(size_t nbr_samples,
@@ -278,7 +278,7 @@ double DistributionLogNormal::probabilityDensity(double x) const
 {
     if (m_scale_param==0.0) return x==m_median ? 1.0 : 0.0;
     double t = std::log(x/m_median)/m_scale_param;
-    return std::exp(-t*t/2.0)/(x*m_scale_param*std::sqrt(2.0*M_PI));
+    return std::exp(-t*t/2.0)/(x*m_scale_param*std::sqrt(Units::PI2));
 }
 
 double DistributionLogNormal::getMean() const
@@ -342,16 +342,16 @@ DistributionCosine::DistributionCosine(double mean, double sigma)
 double DistributionCosine::probabilityDensity(double x) const
 {
     if (m_sigma == 0.0) return x==m_mean ? 1.0 : 0.0;
-    if (std::abs(x-m_mean)>M_PI*m_sigma) return 0.0;
-    return (1.0 + std::cos((x-m_mean)/m_sigma))/(m_sigma*2.0*M_PI);
+    if (std::abs(x-m_mean)>Units::PI*m_sigma) return 0.0;
+    return (1.0 + std::cos((x-m_mean)/m_sigma))/(m_sigma*Units::PI2);
 }
 
 std::vector<double> DistributionCosine::generateValueList(size_t nbr_samples,
         double sigma_factor, const AttLimits &limits) const
 {
     if (sigma_factor <= 0.0 || sigma_factor > 2.0) sigma_factor = 2.0;
-    double xmin = m_mean - sigma_factor*m_sigma*M_PI/2.0;
-    double xmax = m_mean + sigma_factor*m_sigma*M_PI/2.0;
+    double xmin = m_mean - sigma_factor*m_sigma*Units::PID2;
+    double xmax = m_mean + sigma_factor*m_sigma*Units::PID2;
     adjustMinMaxForLimits(xmin, xmax, limits);
     return generateValues(nbr_samples, xmin, xmax);
 }
diff --git a/Core/Algorithms/src/FTDistributions.cpp b/Core/Algorithms/src/FTDistributions.cpp
index 70ca7477999a43176e42846d6786b2ef3b645005..340002fb46e4d704d85b06b98b1dcfc647375859 100644
--- a/Core/Algorithms/src/FTDistributions.cpp
+++ b/Core/Algorithms/src/FTDistributions.cpp
@@ -118,11 +118,11 @@ FTDistribution1DCosine *FTDistribution1DCosine::clone() const
 double FTDistribution1DCosine::evaluate(double q) const
 {
     double qw = std::abs(q*m_omega);
-    if (std::abs(qw/M_PI-1.0) < Numeric::double_epsilon) {
+    if (std::abs(qw/Units::PI-1.0) < Numeric::double_epsilon) {
         return 0.5;
     }
     else {
-        return MathFunctions::Sinc(qw)/(1.0-qw*qw/M_PI/M_PI);
+        return MathFunctions::Sinc(qw)/(1.0-qw*qw/Units::PI/Units::PI);
     }
 }
 
@@ -172,7 +172,7 @@ double IFTDistribution2D::evaluateLattice(double qx, double qy) const
 void IFTDistribution2D::transformToStarBasis(double qX, double qY,
         double alpha, double a, double b, double& qa, double& qb) const
 {
-    double prefactor = 1.0/(2*M_PI); // divide by sin(m_delta)
+    double prefactor = 1.0/Units::PI2; // divide by sin(m_delta)
                                      // for unnormalized X*,Y* basis
     qa = a*prefactor*( std::sin(m_gamma+m_delta)*qX - std::sin(m_gamma)*qY );
     qb = b*prefactor*( -std::sin(alpha-m_gamma-m_delta)*qX + std::sin(alpha-m_gamma)*qY );
diff --git a/Core/Algorithms/src/IInterferenceFunctionStrategy.cpp b/Core/Algorithms/src/IInterferenceFunctionStrategy.cpp
index d1251cf9d0964daa63dd7a00adfc73102c1d26a0..1e59d02e6f125e264f411bb5f58890a5be18f972 100644
--- a/Core/Algorithms/src/IInterferenceFunctionStrategy.cpp
+++ b/Core/Algorithms/src/IInterferenceFunctionStrategy.cpp
@@ -181,7 +181,7 @@ IInterferenceFunctionStrategy::getIntegrationParams(const cvector_t &k_i,
 
     IntegrationParamsAlpha result;
     result.k_i = k_i;
-    result.wavelength = 2.0*M_PI/real_ki.mag();
+    result.wavelength = Units::PI2/real_ki.mag();
     result.alpha_bin = alpha_f_bin;
     result.phi_bin = phi_f_bin;
     return result;
diff --git a/Core/Algorithms/src/LayerStrategyBuilder.cpp b/Core/Algorithms/src/LayerStrategyBuilder.cpp
index 4253f4c26ec1569d6c1bdbd282867bd1ab97cdf2..dfd30fdfcb69f7eb62aa52bfde1a222d97c0aad9 100644
--- a/Core/Algorithms/src/LayerStrategyBuilder.cpp
+++ b/Core/Algorithms/src/LayerStrategyBuilder.cpp
@@ -109,7 +109,7 @@ void LayerStrategyBuilder::collectFormFactorInfos()
     double wavelength = getWavelength();
     double total_abundance = mp_layer->getTotalAbundance();
     if (total_abundance<=0.0) total_abundance = 1.0;
-    complex_t wavevector_scattering_factor = M_PI/wavelength/wavelength;
+    complex_t wavevector_scattering_factor = Units::PI/wavelength/wavelength;
     size_t number_of_particles = p_layout->getNumberOfParticles();
     for (size_t particle_index =
              0; particle_index<number_of_particles; ++particle_index) {
@@ -139,7 +139,7 @@ double LayerStrategyBuilder::getWavelength()
 {
     cvector_t ki = mp_simulation->getInstrument().getBeam().getCentralK();
     kvector_t ki_real(ki.x().real(), ki.y().real(), ki.z().real());
-    return 2*M_PI/ki_real.mag();
+    return Units::PI2/ki_real.mag();
 }
 
 FormFactorInfo *LayerStrategyBuilder::createFormFactorInfo(
diff --git a/Core/FormFactors/inc/IFormFactor.h b/Core/FormFactors/inc/IFormFactor.h
index e4677c829010d9d1039d9b9aba068baad3956e92..77cc743c00d20e54e021a32cb12b4ab91d28deac 100644
--- a/Core/FormFactors/inc/IFormFactor.h
+++ b/Core/FormFactors/inc/IFormFactor.h
@@ -16,6 +16,7 @@
 #ifndef IFORMFACTOR_H
 #define IFORMFACTOR_H
 
+#include "Units.h"
 #include "IMaterial.h"
 #include "ISample.h"
 #include "Bin.h"
@@ -114,7 +115,7 @@ inline double IFormFactor::getHeight() const
 
 inline double IFormFactor::getRadius() const
 {
-    return std::sqrt(getVolume()/getHeight()/M_PI);
+    return std::sqrt(getVolume()/getHeight()/Units::PI);
 }
 
 #endif // IFORMFACTOR_H
diff --git a/Core/FormFactors/inc/IFormFactorBornSeparable.h b/Core/FormFactors/inc/IFormFactorBornSeparable.h
index 193eb86281cbc426bdae7a78fc581de607cce3ec..e778057aaeaaee1eb3d6a29c45d57ed2b6d7f4a3 100644
--- a/Core/FormFactors/inc/IFormFactorBornSeparable.h
+++ b/Core/FormFactors/inc/IFormFactorBornSeparable.h
@@ -50,13 +50,13 @@ inline complex_t IFormFactorBornSeparable::evaluate(const cvector_t& k_i,
     complex_t radial, zpart;
     Bin1DCVector q_bin(k_i - k_f_bin.m_q_lower, k_i - k_f_bin.m_q_upper);
     double delta_qr = std::abs( q_bin.getDelta().magxy() );
-    if ( delta_qr > M_PI/2.0/getRadius() ) {
+    if ( delta_qr > Units::PID2/getRadius() ) {
         radial = bigRadialPart(q_bin);
     } else {
     radial = evaluate_for_q_radial(q_bin.getMidPoint());
     }
     double delta_qz = std::abs( q_bin.getDelta().z() );
-    if ( delta_qz > M_PI/2.0/getHeight() ) {
+    if ( delta_qz > Units::PID2/getHeight() ) {
         zpart = bigZPart(q_bin);
     } else {
         zpart = evaluate_for_q_z(q_bin.getMidPoint());
diff --git a/Core/FormFactors/src/FormFactorCone.cpp b/Core/FormFactors/src/FormFactorCone.cpp
index b5888a63a0288d302d19917666474472694bcac8..f352f027320b5d963a9f0c91fa52ca21da891f3e 100644
--- a/Core/FormFactors/src/FormFactorCone.cpp
+++ b/Core/FormFactors/src/FormFactorCone.cpp
@@ -89,14 +89,14 @@ complex_t FormFactorCone::evaluate_for_q(const cvector_t& q) const
         double tga = std::tan(m_alpha);
         double HdivRtga = H/tga/R;
 
-        return  M_PI/3.0*tga*R*R*R*
+        return  Units::PI/3.0*tga*R*R*R*
                 (1.0 - (1.0 - HdivRtga)*(1.0 - HdivRtga)*(1.0 - HdivRtga));
 
     } else {
 
         complex_t integral = m_integrator->integrate(0., m_height);
 
-        return 2.0*M_PI*integral;
+        return Units::PI2*integral;
     }
 }
 
diff --git a/Core/FormFactors/src/FormFactorCrystal.cpp b/Core/FormFactors/src/FormFactorCrystal.cpp
index 9f4403e23f72f3c4a34e51bef4d2309f2c3daaa7..91ed23a672ea7fb2974e48c3d3b7d8e5f8e1c221 100644
--- a/Core/FormFactors/src/FormFactorCrystal.cpp
+++ b/Core/FormFactors/src/FormFactorCrystal.cpp
@@ -149,6 +149,6 @@ void FormFactorCrystal::calculateLargestReciprocalDistance()
     kvector_t a2 = m_lattice.getBasisVectorB();
     kvector_t a3 = m_lattice.getBasisVectorC();
 
-    m_max_rec_length = std::max(M_PI / a1.mag(), M_PI / a2.mag());
-    m_max_rec_length = std::max(m_max_rec_length, M_PI / a3.mag());
+    m_max_rec_length = std::max(Units::PI / a1.mag(), Units::PI / a2.mag());
+    m_max_rec_length = std::max(m_max_rec_length, Units::PI / a3.mag());
 }
diff --git a/Core/FormFactors/src/FormFactorCylinder.cpp b/Core/FormFactors/src/FormFactorCylinder.cpp
index 80c0f8ee3c31a36656d45e4e8c8bdda8d519f5a6..fb5e2c1b2ea704d248c1b5b47acc0c644fc80bd5 100644
--- a/Core/FormFactors/src/FormFactorCylinder.cpp
+++ b/Core/FormFactors/src/FormFactorCylinder.cpp
@@ -59,7 +59,7 @@ complex_t FormFactorCylinder::evaluate_for_q(const cvector_t& q) const
     complex_t J1_qrR_div_qrR = std::abs(qrR) > Numeric::double_epsilon ?
         MathFunctions::Bessel_C1(qrR) :
         0.5;
-    complex_t radial_part = 2*M_PI*R*R*J1_qrR_div_qrR;
+    complex_t radial_part = Units::PI2*R*R*J1_qrR_div_qrR;
 
     complex_t result = radial_part*z_part;
 
diff --git a/Core/FormFactors/src/FormFactorDecoratorMaterial.cpp b/Core/FormFactors/src/FormFactorDecoratorMaterial.cpp
index b0fea9075fd8ea457b964efcbc139a2a23c34790..8345771da511c9bb73aae10cb037e7352e05079d 100644
--- a/Core/FormFactors/src/FormFactorDecoratorMaterial.cpp
+++ b/Core/FormFactors/src/FormFactorDecoratorMaterial.cpp
@@ -74,7 +74,7 @@ Eigen::Matrix2cd FormFactorDecoratorMaterial::evaluatePol(const cvector_t& k_i,
     time_reverse_conj(0,1) = 1.0;
     time_reverse_conj(1,0) = -1.0;
     // the interaction and time reversal taken together:
-    double k_mag2 = 4.0 * M_PI * m_wavevector_scattering_factor.real();
+    double k_mag2 = 4.0 * Units::PI * m_wavevector_scattering_factor.real();
     Eigen::Matrix2cd V_eff = m_wavevector_scattering_factor * time_reverse_conj
             * (mP_material->getScatteringMatrix(k_mag2) -
                mP_ambient_material->getScatteringMatrix(k_mag2));
diff --git a/Core/FormFactors/src/FormFactorEllipsoidalCylinder.cpp b/Core/FormFactors/src/FormFactorEllipsoidalCylinder.cpp
index 437f3dee5fd6592ede3573c57f78d7be1d138852..60903930f42090b4d3c2f2cd02064c315517740e 100644
--- a/Core/FormFactors/src/FormFactorEllipsoidalCylinder.cpp
+++ b/Core/FormFactors/src/FormFactorEllipsoidalCylinder.cpp
@@ -44,7 +44,7 @@ complex_t FormFactorEllipsoidalCylinder::evaluate_for_q(const cvector_t& q) cons
     complex_t gamma  = std::sqrt((qxRa)*(qxRa) + (qyRb)*(qyRb));
     complex_t J1_gamma_div_gamma = MathFunctions::Bessel_C1(gamma);
 
-    return 2.*M_PI *m_radius_a*m_radius_b*m_height * Fz*J1_gamma_div_gamma;
+    return Units::PI2 *m_radius_a*m_radius_b*m_height * Fz*J1_gamma_div_gamma;
 
 }
 
diff --git a/Core/FormFactors/src/FormFactorFullSphere.cpp b/Core/FormFactors/src/FormFactorFullSphere.cpp
index 61b47f4d8076662731f89cc471b17a9a35342b91..a3b77b3292e7c27fae191853856f4c56c66eaf18 100644
--- a/Core/FormFactors/src/FormFactorFullSphere.cpp
+++ b/Core/FormFactors/src/FormFactorFullSphere.cpp
@@ -52,7 +52,7 @@ complex_t FormFactorFullSphere::evaluate_for_q(const cvector_t& q) const
     complex_t z_part = std::exp(iqzR);
 
     complex_t qR = std::sqrt( std::norm(q.x()) + std::norm(q.y()) + std::norm(qz) )*R;
-    double volume = 4*M_PI*R*R*R/3;
+    double volume = 4*Units::PI*R*R*R/3;
     complex_t radial;
     if (std::abs(qR) < Numeric::double_epsilon) {
         radial = volume;
diff --git a/Core/FormFactors/src/FormFactorFullSpheroid.cpp b/Core/FormFactors/src/FormFactorFullSpheroid.cpp
index e2dfa6bbbe0c2f99a8a74cf23874daab8866bf27..0a805937915a567e858fe202a730f722f3f2ab56 100644
--- a/Core/FormFactors/src/FormFactorFullSpheroid.cpp
+++ b/Core/FormFactors/src/FormFactorFullSpheroid.cpp
@@ -76,14 +76,14 @@ complex_t FormFactorFullSpheroid::evaluate_for_q(const cvector_t& q) const
 
     if (std::abs(m_q.mag()) <= Numeric::double_epsilon) {
 
-        return 2.*M_PI*R*R*H/3.;
+        return Units::PI2*R*R*H/3.;
 
     } else {
 
         complex_t qzH_half  = m_q.z()*H/2.0;
         complex_t z_part    =  std::exp(complex_t(0.0, 1.0)*qzH_half);
 
-        return 4.0* M_PI * z_part *m_integrator->integrate(0.0, H/2.0);
+        return 4.0* Units::PI * z_part *m_integrator->integrate(0.0, H/2.0);
 
     }
 }
diff --git a/Core/FormFactors/src/FormFactorGauss.cpp b/Core/FormFactors/src/FormFactorGauss.cpp
index e01ab0a6da2297dbaadfa833f334635d94a14cc5..960932dfe6ffa686c19940404df801fa7d79ff30 100644
--- a/Core/FormFactors/src/FormFactorGauss.cpp
+++ b/Core/FormFactors/src/FormFactorGauss.cpp
@@ -25,7 +25,7 @@ FormFactorGauss::FormFactorGauss(double volume)
     m_width = m_height;
     check_initialization();
     init_parameters();
-    m_max_ql = std::sqrt(-4.0*M_PI*std::log(Numeric::double_epsilon)
+    m_max_ql = std::sqrt(-4.0*Units::PI*std::log(Numeric::double_epsilon)
                / 3.0);
 }
 
@@ -36,7 +36,7 @@ FormFactorGauss::FormFactorGauss(double width, double height)
     m_height = height;
     check_initialization();
     init_parameters();
-    m_max_ql = std::sqrt(-4.0*M_PI*std::log(Numeric::double_epsilon)
+    m_max_ql = std::sqrt(-4.0*Units::PI*std::log(Numeric::double_epsilon)
                / 3.0);
 }
 
@@ -71,9 +71,9 @@ complex_t FormFactorGauss::evaluate_for_q(const cvector_t& q) const
 
 
     complex_t z_part = std::exp(complex_t(0.,1.)*qzHdiv2) * m_height
-            * std::exp(-qzh*qzh/4.0/M_PI);
+            * std::exp(-qzh*qzh/4.0/Units::PI);
     double radial_part = m_width * m_width
-            * std::exp(-(qxr*qxr+qyr*qyr)/4.0/M_PI);
+            * std::exp(-(qxr*qxr+qyr*qyr)/4.0/Units::PI);
     complex_t result = radial_part * z_part;
     return result;
 }
diff --git a/Core/FormFactors/src/FormFactorHemiEllipsoid.cpp b/Core/FormFactors/src/FormFactorHemiEllipsoid.cpp
index b22becd1f99ffdacfeda676ff0e6e923f2470728..cd297044423c760073d775072be04f2c9febfbfe 100644
--- a/Core/FormFactors/src/FormFactorHemiEllipsoid.cpp
+++ b/Core/FormFactors/src/FormFactorHemiEllipsoid.cpp
@@ -91,11 +91,11 @@ complex_t FormFactorHemiEllipsoid::evaluate_for_q(const cvector_t& q) const
 
      if (std::abs(m_q.mag()) <= Numeric::double_epsilon) {
 
-         return 2.0*M_PI*R*W*H/3.;
+         return Units::PI2*R*W*H/3.;
 
      } else {
 
-         return 2.0*M_PI*m_integrator->integrate(0.,H );
+         return Units::PI2*m_integrator->integrate(0.,H );
      }
 }
 
diff --git a/Core/FormFactors/src/FormFactorInfLongBox.cpp b/Core/FormFactors/src/FormFactorInfLongBox.cpp
index 5bd7fa9a45bae649f14c977c8fc587792eb47b81..34f4258e9ec5525d476089418c0e61f4a362c133 100644
--- a/Core/FormFactors/src/FormFactorInfLongBox.cpp
+++ b/Core/FormFactors/src/FormFactorInfLongBox.cpp
@@ -43,7 +43,7 @@ complex_t FormFactorInfLongBox::evaluate(const cvector_t& k_i,
     complex_t qyWdiv2 = m_width*q.y()/2.0;
     complex_t qzHdiv2 = m_height*q.z()/2.0;
 
-    return 2*M_PI*m_height*m_width*
+    return Units::PI2*m_height*m_width*
         std::exp(complex_t(0.,1.)*qzHdiv2) *
         MathFunctions::Sinc(qyWdiv2) *
         MathFunctions::Sinc(qzHdiv2);
@@ -69,7 +69,6 @@ complex_t FormFactorInfLongBox::evaluate_for_q(const cvector_t& q) const
 }
 
 double FormFactorInfLongBox::getVolume() const {
-    // return 2*M_PI*m_height*m_width;
     // volume of the infinite object is infinite
     throw NotImplementedException(
                 "FormFactorInfLongBox::getVolume() -> Error: not implemented exception. Volume of the infinite object is infinite.");
diff --git a/Core/FormFactors/src/FormFactorInfLongRipple1.cpp b/Core/FormFactors/src/FormFactorInfLongRipple1.cpp
index 0ffe28467eeedfea8e0c6af2aa3bc121a11e3ed5..1a71ff95a4cb3e5ecd3385e5cb0d70c70ef1174b 100644
--- a/Core/FormFactors/src/FormFactorInfLongRipple1.cpp
+++ b/Core/FormFactors/src/FormFactorInfLongRipple1.cpp
@@ -70,7 +70,7 @@ complex_t FormFactorInfLongRipple1::Integrand(double Z, void* params) const
     (void)params;  // to avoid unused-variable warning
     complex_t iqZ = complex_t(0.0, 1.0)*m_q.z()*Z;
     complex_t aa = std::acos(2.0*Z/m_height - 1.0);
-    return std::exp(iqZ)*aa*MathFunctions::Sinc(aa*m_q.y()*m_width/(2*M_PI));
+    return std::exp(iqZ)*aa*MathFunctions::Sinc(aa*m_q.y()*m_width/(Units::PI2));
 }
 
 //! Complex formfactor.
@@ -104,7 +104,6 @@ complex_t FormFactorInfLongRipple1::evaluate(const cvector_t& k_i,
 }
 
 double FormFactorInfLongRipple1::getVolume() const {
-    // return 2*M_PI*m_height*m_width;
     // volume of the infinite object is infinite
     throw NotImplementedException(
                 "FormFactorInfLongRipple1::getVolume() -> Error: not implemented exception. Volume of the infinite object is infinite.");
diff --git a/Core/FormFactors/src/FormFactorInfLongRipple2.cpp b/Core/FormFactors/src/FormFactorInfLongRipple2.cpp
index 4768b32f3644b35d5a111d88e4696d79c88f6297..e1a27f726e631599318590bf9215955b8c1fa44c 100644
--- a/Core/FormFactors/src/FormFactorInfLongRipple2.cpp
+++ b/Core/FormFactors/src/FormFactorInfLongRipple2.cpp
@@ -75,7 +75,6 @@ FormFactorInfLongRipple2 *FormFactorInfLongRipple2::clone() const
 }
 
 double FormFactorInfLongRipple2::getVolume() const {
-    // return 2*M_PI*m_height*m_width;
     // volume of the infinite object is infinite
     throw NotImplementedException(
                 "FormFactorInfLongRipple2::getVolume() -> Error: not implemented exception. Volume of the infinite object is infinite.");
@@ -106,7 +105,7 @@ complex_t FormFactorInfLongRipple2::evaluate(const cvector_t& k_i,
 
     m_q = k_i - k_f_bin.getMidPoint();
 
-    complex_t factor = 2.0*M_PI*m_width;
+    complex_t factor = Units::PI2*m_width;
     complex_t result = 0;
     complex_t iqzH = complex_t(0.0, 1.0)*m_q.z()*m_height;
     complex_t iqyW = complex_t(0.0, 1.0)*m_q.y()*m_width;
diff --git a/Core/FormFactors/src/FormFactorLorentz.cpp b/Core/FormFactors/src/FormFactorLorentz.cpp
index 663f0a9c8971f3f80c3c5bf68893d3258b57fd82..7b9098d8faf14dd9b914d6e64db678dfaaec66a3 100644
--- a/Core/FormFactors/src/FormFactorLorentz.cpp
+++ b/Core/FormFactors/src/FormFactorLorentz.cpp
@@ -58,7 +58,7 @@ FormFactorLorentz* FormFactorLorentz::clone() const
 
 complex_t FormFactorLorentz::evaluate_for_q(const cvector_t& q) const
 {
-    static const double sigma2 = 4.0*std::pow(M_PI, 2.0/3.0);
+    static const double sigma2 = 4.0*std::pow(Units::PI, 2.0/3.0);
     double R = m_width;
     double H = m_height;
 
diff --git a/Core/FormFactors/src/FormFactorRipple1.cpp b/Core/FormFactors/src/FormFactorRipple1.cpp
index a0165371b920f7f003bebaab2685f10722707ce5..a1bea2bf9a8af395f05d6024d6ea75b1dc21b8aa 100644
--- a/Core/FormFactors/src/FormFactorRipple1.cpp
+++ b/Core/FormFactors/src/FormFactorRipple1.cpp
@@ -73,7 +73,7 @@ complex_t FormFactorRipple1::Integrand(double Z, void* params) const
     (void)params;  // to avoid unused-variable warning
     complex_t iqZ = complex_t(0.0, 1.0)*m_q.z()*Z;
     complex_t aa = std::acos(2.0*Z/m_height - 1.0);
-    return std::exp(iqZ)*aa*MathFunctions::Sinc(aa*m_q.y()*m_width/(2*M_PI));
+    return std::exp(iqZ)*aa*MathFunctions::Sinc(aa*m_q.y()*m_width/(Units::PI2));
 }
 
 //! Complex formfactor.
@@ -81,17 +81,17 @@ complex_t FormFactorRipple1::Integrand(double Z, void* params) const
 complex_t FormFactorRipple1::evaluate_for_q(const cvector_t& q) const
 {
     m_q = q;
-    complex_t factor = m_length*MathFunctions::Sinc(m_q.x()*m_length*0.5)*m_width/M_PI;
-    complex_t aaa = m_q.y()*m_width/(2*M_PI);
+    complex_t factor = m_length*MathFunctions::Sinc(m_q.x()*m_length*0.5)*m_width/Units::PI;
+    complex_t aaa = m_q.y()*m_width/(Units::PI2);
     complex_t aaa2 = aaa*aaa;
 
     // analytical expressions for some particular cases
     if (0.0==m_q.y() && 0.0==m_q.z())
-        return factor*M_PI_2*m_height;
+        return factor*Units::PID2*m_height;
     else if (0.0==m_q.z() && 1.0 == aaa2)
-        return factor*M_PI_4*m_height;
+        return factor*Units::PID4*m_height;
     else if (0.0==m_q.z())
-        return factor*M_PI_2*m_height*MathFunctions::Sinc(m_q.y()*m_width*0.5)/(1.0-aaa2);
+        return factor*Units::PID2*m_height*MathFunctions::Sinc(m_q.y()*m_width*0.5)/(1.0-aaa2);
 
     // numerical integration otherwise
     complex_t integral = m_integrator->integrate(0, m_height);
diff --git a/Core/FormFactors/src/FormFactorSphereUniformRadius.cpp b/Core/FormFactors/src/FormFactorSphereUniformRadius.cpp
index d132d58a7a9a1d8bf0488b1efbada46c72c16363..d27527f3133949f6098317469d540d0264f8fc0f 100644
--- a/Core/FormFactors/src/FormFactorSphereUniformRadius.cpp
+++ b/Core/FormFactors/src/FormFactorSphereUniformRadius.cpp
@@ -55,11 +55,11 @@ complex_t FormFactorSphereUniformRadius::evaluate_for_q(
     double q2 = std::norm(q.x()) + std::norm(q.y()) + std::norm(q.z());
     double q_r = std::sqrt(q2);
     if (q_r*R < Numeric::double_epsilon) {
-        return (4.0*M_PI*R*R*R + M_PI*R*W*W)/3.0;
+        return (4.0*Units::PI*R*R*R + Units::PI*R*W*W)/3.0;
     }
     double qR = q_r*R;
     double qW = q_r*W;
-    double nominator = 4*M_PI*( 4*std::sin(qR)*std::sin(qW/2.0)
+    double nominator = 4*Units::PI*( 4*std::sin(qR)*std::sin(qW/2.0)
                              - qW*std::cos(qW/2.0)*std::sin(qR)
                              - 2.0*qR*std::cos(qR)*std::sin(qW/2.0));
     return nominator/(q2*q2*W);
diff --git a/Core/FormFactors/src/FormFactorTruncatedSphere.cpp b/Core/FormFactors/src/FormFactorTruncatedSphere.cpp
index 868b198dc1ce5ee8655277d248e1d4d7324d429a..c5eb57a94d88620d32696d59f9573b13ab897711 100644
--- a/Core/FormFactors/src/FormFactorTruncatedSphere.cpp
+++ b/Core/FormFactors/src/FormFactorTruncatedSphere.cpp
@@ -78,13 +78,13 @@ complex_t FormFactorTruncatedSphere::evaluate_for_q(const cvector_t& q) const
 {   m_q = q;
     if ( std::abs(m_q.mag()) < Numeric::double_epsilon) {
         double HdivR = m_height/m_radius;
-        return M_PI/3.*m_radius*m_radius*m_radius
+        return Units::PI/3.*m_radius*m_radius*m_radius
                 *(3.*HdivR -1. - (HdivR - 1.)*(HdivR - 1.)*(HdivR - 1.));
     }
     else {
     complex_t iqzR = complex_t(0.0, 1.0)*m_q.z()*(m_height-m_radius);
     complex_t integral = m_integrator->integrate(m_radius-m_height, m_radius);
-    return 2*M_PI*integral*std::exp(iqzR);
+    return Units::PI2*integral*std::exp(iqzR);
     }
 }
 
diff --git a/Core/FormFactors/src/FormFactorTruncatedSpheroid.cpp b/Core/FormFactors/src/FormFactorTruncatedSpheroid.cpp
index 681142f06dc787a033be12c9b364d469835e6fb9..889e835674dee7ce6fae2c82a64211c6f442bc0a 100644
--- a/Core/FormFactors/src/FormFactorTruncatedSpheroid.cpp
+++ b/Core/FormFactors/src/FormFactorTruncatedSpheroid.cpp
@@ -88,13 +88,13 @@ complex_t FormFactorTruncatedSpheroid::evaluate_for_q(const cvector_t& q) const
 
     if (std::abs(m_q.mag()) <= Numeric::double_epsilon) {
 
-        return M_PI*R*H*H/fp*(1.-H/(3.*fp*R));
+        return Units::PI*R*H*H/fp*(1.-H/(3.*fp*R));
 
     } else {
 
         complex_t z_part    =  std::exp(complex_t(0.0, 1.0)*m_q.z()*(H-fp*R));
 
-        return 2.0* M_PI * z_part *m_integrator->integrate(fp*R-H,fp*R );
+        return Units::PI2 * z_part *m_integrator->integrate(fp*R-H,fp*R );
     }
 }
 
diff --git a/Core/FormFactors/src/IFormFactorBorn.cpp b/Core/FormFactors/src/IFormFactorBorn.cpp
index 915282d0102227dd24001c7daea76c7e1fcff358..be7c153b5b36780bf00a6fa46cb41ba4432686df 100644
--- a/Core/FormFactors/src/IFormFactorBorn.cpp
+++ b/Core/FormFactors/src/IFormFactorBorn.cpp
@@ -86,9 +86,9 @@ bool IFormFactorBorn::useLargeBinApproximation(const Bin1DCVector& q_bin) const
     double delta_qz = std::abs( q_bin.getDelta().z() );
     if (delta_qr == 0 || delta_qz == 0)
         return false;
-    if ( delta_qr > M_PI/2./getRadius() )
+    if ( delta_qr > Units::PID2/getRadius() )
         return true;
-    if ( delta_qz > M_PI/2./getHeight() )
+    if ( delta_qz > Units::PID2/getHeight() )
         return true;
     return false;
 }
diff --git a/Core/Geometry/inc/BasicVector3D.h b/Core/Geometry/inc/BasicVector3D.h
index f9f5b4a08d8d1c014e7c6f5ba0d15859c8b80f5c..24351b78dd44be42a1552c84145ae28e16e56b41 100644
--- a/Core/Geometry/inc/BasicVector3D.h
+++ b/Core/Geometry/inc/BasicVector3D.h
@@ -22,6 +22,7 @@
 #define GEOMETRY_BASICVECTOR3D_H
 
 #include "Numeric.h"
+#include "Units.h"
 #include "Exceptions.h"
 #include "WinDllMacros.h"
 #include <complex>
@@ -223,7 +224,7 @@ public:
     inline void setLambdaAlphaPhi(
         const T& _lambda, const T& _alpha, const T& _phi)
         {
-            T k = 2*M_PI/_lambda;
+            T k = Units::PI2/_lambda;
             v_[0] = k*std::cos(_alpha) * std::cos(_phi);
             v_[1] = k*std::cos(_alpha) * std::sin(_phi);
             v_[2] = k*std::sin(_alpha);
diff --git a/Core/StandardSamples/IsGISAXS06Builder.cpp b/Core/StandardSamples/IsGISAXS06Builder.cpp
index 48952688790539d97189be4000d9273c36bf7f43..86a3b7881bfb86ebb7e46f6751893ea376b07ed5 100644
--- a/Core/StandardSamples/IsGISAXS06Builder.cpp
+++ b/Core/StandardSamples/IsGISAXS06Builder.cpp
@@ -44,8 +44,8 @@ ISample *IsGISAXS06Lattice1Builder::buildSample() const
 
     InterferenceFunction2DLattice *p_interference_function =
         InterferenceFunction2DLattice::createSquare(10.0*Units::nanometer);
-    FTDistribution2DCauchy pdf(300.0*Units::nanometer/2.0/M_PI,
-                               100.0*Units::nanometer/2.0/M_PI);
+    FTDistribution2DCauchy pdf(300.0*Units::nanometer/2.0/Units::PI,
+                               100.0*Units::nanometer/2.0/Units::PI);
     p_interference_function->setProbabilityDistribution(pdf);
 
     // particles
@@ -80,9 +80,9 @@ ISample *IsGISAXS06Lattice2Builder::buildSample() const
     Layer substrate_layer(substrate_material);
 
     InterferenceFunction2DLattice interference_function(10.0*Units::nanometer,
-            10.0*Units::nanometer, M_PI/2.0);
-    FTDistribution2DCauchy pdf(300.0*Units::nanometer/2.0/M_PI,
-                               100.0*Units::nanometer/2.0/M_PI);
+            10.0*Units::nanometer, Units::PI/2.0);
+    FTDistribution2DCauchy pdf(300.0*Units::nanometer/2.0/Units::PI,
+                               100.0*Units::nanometer/2.0/Units::PI);
     interference_function.setProbabilityDistribution(pdf);
 
     FormFactorCylinder ff_cyl(5.0*Units::nanometer, 5.0*Units::nanometer);
@@ -125,8 +125,8 @@ ISample *IsGISAXS06Lattice3Builder::buildSample() const
     InterferenceFunction2DLattice *p_interference_function =
         InterferenceFunction2DLattice::createSquare(10.0*Units::nanometer,
                                                     30.0*Units::degree);
-    FTDistribution2DCauchy pdf(300.0*Units::nanometer/2.0/M_PI,
-                               100.0*Units::nanometer/2.0/M_PI);
+    FTDistribution2DCauchy pdf(300.0*Units::nanometer/2.0/Units::PI,
+                               100.0*Units::nanometer/2.0/Units::PI);
     pdf.setGamma(30.0*Units::degree);
     p_interference_function->setProbabilityDistribution(pdf);
 
@@ -177,8 +177,8 @@ ISample *IsGISAXS06Lattice4Builder::buildSample() const
     InterferenceFunction2DLattice *p_interference_function =
         InterferenceFunction2DLattice::createSquare(10.0*Units::nanometer,
                                                     m_xi);
-    FTDistribution2DCauchy pdf(300.0*Units::nanometer/2.0/M_PI,
-                               100.0*Units::nanometer/2.0/M_PI);
+    FTDistribution2DCauchy pdf(300.0*Units::nanometer/2.0/Units::PI,
+                               100.0*Units::nanometer/2.0/Units::PI);
     p_interference_function->setProbabilityDistribution(pdf);
 
     ParticleLayout particle_layout;
diff --git a/Core/StandardSamples/IsGISAXS08Builder.cpp b/Core/StandardSamples/IsGISAXS08Builder.cpp
index 33d6bf94108eaf44766db3a9bfe1980db2640f43..ec052cb1c81fd47be7636a3be020dbfa2993bc17 100644
--- a/Core/StandardSamples/IsGISAXS08Builder.cpp
+++ b/Core/StandardSamples/IsGISAXS08Builder.cpp
@@ -38,7 +38,7 @@ ISample *IsGISAXS08ABuilder::buildSample() const
 
     InterferenceFunction2DParaCrystal *p_interference_function =
             new InterferenceFunction2DParaCrystal(10.0*Units::nanometer,
-                    10.0*Units::nanometer, M_PI/2.0, 0.0, 0.0);
+                    10.0*Units::nanometer, Units::PI/2.0, 0.0, 0.0);
     p_interference_function->setDomainSizes(20.0*Units::micrometer,
             20.0*Units::micrometer);
     FTDistribution2DCauchy pdf1(0.5*Units::nanometer, 2.0*Units::nanometer);
@@ -81,7 +81,7 @@ ISample *IsGISAXS08BBuilder::buildSample() const
 
     InterferenceFunction2DParaCrystal *p_interference_function =
             new InterferenceFunction2DParaCrystal(10.0*Units::nanometer,
-                    10.0*Units::nanometer, M_PI/2.0, 0.0, 0.0);
+                    10.0*Units::nanometer, Units::PI/2.0, 0.0, 0.0);
     p_interference_function->setDomainSizes(20.0*Units::micrometer,
             20.0*Units::micrometer);
     FTDistribution2DCauchy pdf1(0.5*Units::nanometer, 0.5*Units::nanometer);
diff --git a/Core/StandardSamples/MesoCrystal01Builder.cpp b/Core/StandardSamples/MesoCrystal01Builder.cpp
index 6a037e0248680ba9dfcf8562c7aac6865864e112..e50d5e67799a17fb12ed2a67ae5e2fe6f22a910a 100644
--- a/Core/StandardSamples/MesoCrystal01Builder.cpp
+++ b/Core/StandardSamples/MesoCrystal01Builder.cpp
@@ -68,7 +68,7 @@ void MesoCrystal01Builder::init_parameters()
 ISample* MesoCrystal01Builder::buildSample() const
 {
     // create mesocrystal
-    double surface_density = m_surface_filling_ratio/M_PI/m_meso_radius/m_meso_radius;
+    double surface_density = m_surface_filling_ratio/Units::PI/m_meso_radius/m_meso_radius;
 //    complex_t n_particle(1.0-1.55e-5, 1.37e-6); // data from Artur
     complex_t n_particle(1.0-2.84e-5, 4.7e-7); // data from http://henke.lbl.gov/optical_constants/getdb2.html
     complex_t avg_n_squared_meso = 0.7886*n_particle*n_particle + 0.2114;
@@ -101,7 +101,7 @@ ISample* MesoCrystal01Builder::buildSample() const
 //    double alpha_step = 5.0*Units::degree/n_alpha_rotation_steps;
 //    double alpha_start = - (n_alpha_rotation_steps/2.0)*alpha_step;
 
-    double phi_step = 2*M_PI/3.0/n_max_phi_rotation_steps;
+    double phi_step = 2*Units::PI/3.0/n_max_phi_rotation_steps;
     double phi_start = 0.0;
     for (size_t i=0; i<n_max_phi_rotation_steps; ++i) {
         for (size_t j=0; j<n_alpha_rotation_steps; ++j) {
diff --git a/Core/Tools/inc/MathFunctions.h b/Core/Tools/inc/MathFunctions.h
index 9912ef73597d41e5d1b32cfa5889edf3a34e2700..17ec3fef9d6c042197246c9d30726167c183a90e 100644
--- a/Core/Tools/inc/MathFunctions.h
+++ b/Core/Tools/inc/MathFunctions.h
@@ -212,7 +212,7 @@ inline double MathFunctions::Si(double value)  // int_0^x du Sin(u)/u
 
 inline double MathFunctions::Sinc(double value)  // Sin(x)/x
 {
-    return gsl_sf_sinc(value/M_PI);
+    return gsl_sf_sinc(value/Units::PI);
 }
 
 inline complex_t MathFunctions::Sinc(const complex_t &value)  // Sin(x)/x
@@ -236,7 +236,7 @@ inline double MathFunctions::FastSin(const double& x) {
     const double P = 0.225f;
     const double A = 16 * std::sqrt(P);
     const double B = (1 - P) / std::sqrt(P);
-    double y = x / (2 * M_PI);
+    double y = x / (2 * Units::PI);
     y = y - std::floor(y + 0.5);  // y in range -0.5..0.5
     y = A * y * (0.5 - std::abs(y));
     return y * (B + std::abs(y));
@@ -247,7 +247,7 @@ inline double MathFunctions::FastCos(const double& x) {
     const double P = 0.225f;
     const double A = 16 * std::sqrt(P);
     const double B = (1 - P) / std::sqrt(P);
-    double y = x / (2 * M_PI)+0.25; // pi/2. shift
+    double y = x / (2 * Units::PI)+0.25; // pi/2. shift
     y = y - std::floor(y + 0.5);  // y in range -0.5..0.5
     y = A * y * (0.5 - std::abs(y));
     return y * (B + std::abs(y));
diff --git a/Core/Tools/inc/Numeric.h b/Core/Tools/inc/Numeric.h
index bd27e8efdbcd1054d1df1bdb3f740b1f2debec8b..e0af341e31a32a9b69736ec66cc5fe897c5a296f 100644
--- a/Core/Tools/inc/Numeric.h
+++ b/Core/Tools/inc/Numeric.h
@@ -21,10 +21,6 @@
 
 //! Floating-point epsilon, tolerances, almost-equal.
 
-//#ifndef M_PI
-//#define M_PI		3.14159265358979323846
-//#endif
-
 namespace Numeric {
 
 static const double required_precision = 1.e-4;
diff --git a/Core/Tools/inc/Units.h b/Core/Tools/inc/Units.h
index be32d6d0eb8b9537da939db75f539d1c68d0de94..51b54b6fb9a48d54f9218c4e89a70446089eb089 100644
--- a/Core/Tools/inc/Units.h
+++ b/Core/Tools/inc/Units.h
@@ -16,7 +16,7 @@
 #ifndef UNITS_H
 #define UNITS_H
 
-#include "Types.h"
+//#include "Types.h"
 
 namespace Units {  // BornAgain namespace
 
@@ -24,6 +24,7 @@ namespace Units {  // BornAgain namespace
 static const double PI          = 3.14159265358979323846264338327950288;
 static const double PI2         = 6.28318530717958647692528676655900577;
 static const double PID2        = 1.57079632679489661923132169163975144;
+static const double PID4        = PI/4.0;
 static const double PI_SQR      = 9.86960440108935861883449099987615114;
 
 // Length
diff --git a/Core/Tools/src/MathFunctions.cpp b/Core/Tools/src/MathFunctions.cpp
index f40c9122f564c649e7a27d106c26efcc633ecec5..5aae437ce1554c7afbbdac20bece622ba752c413 100644
--- a/Core/Tools/src/MathFunctions.cpp
+++ b/Core/Tools/src/MathFunctions.cpp
@@ -36,7 +36,7 @@ double MathFunctions::IntegratedGaussian(double value, double average, double st
 
 double MathFunctions::StandardNormal(double value)
 {
-    return std::exp(-value * value / 2.0) / std::sqrt(2 * M_PI);
+    return std::exp(-value * value / 2.0) / std::sqrt(Units::PI2);
 }
 
 double MathFunctions::GenerateStandardNormalRandom() // using GSL
@@ -185,7 +185,7 @@ complex_t MathFunctions::crbond_bessel_J0(const complex_t &z)
             kz = 10; //   "      "     "  12
         else
             kz = 12; //   "      "     "  14
-        complex_t ct1 = z1 - M_PI_4;
+        complex_t ct1 = z1 - Units::PID4;
         complex_t cp0 = cone;
         complex_t cq0 = -0.125 / z1;
         complex_t ptmp = std::pow(z1, -2.0);
@@ -258,7 +258,7 @@ complex_t MathFunctions::crbond_bessel_J1(const complex_t &z)
             cq1 += b1[k] * ptmp / z1;
             ptmp /= (z1 * z1);
         }
-        complex_t ct2 = z1 - 0.75 * M_PI;
+        complex_t ct2 = z1 - 0.75 * Units::PI;
         cj1 = std::sqrt(M_2_PI / z1) * (cp1 * std::cos(ct2) - cq1 * std::sin(ct2));
     }
     if (std::real(z) < 0.0)
diff --git a/Core/Tools/src/PyGenTools.cpp b/Core/Tools/src/PyGenTools.cpp
index 1a1b09ae201859a4b6772c89bce20c6188d5a705..a808681bce2c0e8045fe561bf207e04c5cd23b80 100644
--- a/Core/Tools/src/PyGenTools.cpp
+++ b/Core/Tools/src/PyGenTools.cpp
@@ -44,7 +44,7 @@ std::string PyGenTools::genPyScript(Simulation *simulation)
 
 bool PyGenTools::isSquare(double length1, double length2, double angle)
 {
-    if(length1 == length2 && Numeric::areAlmostEqual(angle, M_PI/2.0)) {
+    if(length1 == length2 && Numeric::areAlmostEqual(angle, Units::PI/2.0)) {
         return true;
     }
     return false;
@@ -53,7 +53,7 @@ bool PyGenTools::isSquare(double length1, double length2, double angle)
 
 bool PyGenTools::isHexagonal(double length1, double length2, double angle)
 {
-    if(length1 == length2 && Numeric::areAlmostEqual(angle, 2*M_PI/3.0)) {
+    if(length1 == length2 && Numeric::areAlmostEqual(angle, 2*Units::PI/3.0)) {
         return true;
     }
     return false;
diff --git a/GUI/coregui/Models/TransformFromDomain.cpp b/GUI/coregui/Models/TransformFromDomain.cpp
index ead2f1495ebbea340b6b6cf0a01de10fdbea9866..849dede8e42120e93557d0309228409c1c14afba 100644
--- a/GUI/coregui/Models/TransformFromDomain.cpp
+++ b/GUI/coregui/Models/TransformFromDomain.cpp
@@ -522,7 +522,7 @@ bool TransformFromDomain::isValidRoughness(const LayerRoughness *roughness)
 bool TransformFromDomain::isSquareLattice(double length1, double length2,
                                           double angle)
 {
-    if(length1 == length2 && Numeric::areAlmostEqual(angle, M_PI/2.0)) {
+    if(length1 == length2 && Numeric::areAlmostEqual(angle, Units::PI/2.0)) {
         return true;
     }
     return false;
@@ -532,7 +532,7 @@ bool TransformFromDomain::isSquareLattice(double length1, double length2,
 bool TransformFromDomain::isHexagonalLattice(double length1, double length2,
                                              double angle)
 {
-    if(length1 == length2 && Numeric::areAlmostEqual(angle, 2*M_PI/3.0)) {
+    if(length1 == length2 && Numeric::areAlmostEqual(angle, 2*Units::PI/3.0)) {
         return true;
     }
     return false;
diff --git a/GUI/coregui/Models/TransformToDomain.cpp b/GUI/coregui/Models/TransformToDomain.cpp
index 9bd3878754eee16049a33cb4f8f1445aa10e5d9c..6db3935a2ba448040c78ac76a8c047b1f8022dba 100644
--- a/GUI/coregui/Models/TransformToDomain.cpp
+++ b/GUI/coregui/Models/TransformToDomain.cpp
@@ -193,13 +193,13 @@ IInterferenceFunction *TransformToDomain::createInterferenceFunction(
             length_1 = latticeItem->getRegisteredProperty(
                         SquareLatticeTypeItem::P_LATTICE_LENGTH).toDouble();
             length_2 = length_1;
-            alpha_lattice = M_PI/2.0;
+            alpha_lattice = Units::PI/2.0;
         }
         else if(latticeItem->modelType() == Constants::HexagonalLatticeType) {
             length_1 = latticeItem->getRegisteredProperty(
                         HexagonalLatticeTypeItem::P_LATTICE_LENGTH).toDouble();
             length_2 = length_1;
-            alpha_lattice = 2*M_PI/3.0;
+            alpha_lattice = 2*Units::PI/3.0;
         }
         else {
             throw GUIHelpers::Error("TransformToDomain::createInterferenceFunction() -> Error");
@@ -250,13 +250,13 @@ IInterferenceFunction *TransformToDomain::createInterferenceFunction(
             length_1 = latticeItem->getRegisteredProperty(
                         SquareLatticeTypeItem::P_LATTICE_LENGTH).toDouble();
             length_2 = length_1;
-            angle = M_PI/2.0;
+            angle = Units::PI/2.0;
         }
         else if(latticeItem->modelType() == Constants::HexagonalLatticeType) {
             length_1 = latticeItem->getRegisteredProperty(
                         HexagonalLatticeTypeItem::P_LATTICE_LENGTH).toDouble();
             length_2 = length_1;
-            angle = 2*M_PI/3.0;
+            angle = 2*Units::PI/3.0;
         }
         else {
             throw GUIHelpers::Error("TransformToDomain::createInterferenceFunction() -> Error");
diff --git a/Tests/UnitTests/TestCore/BeamTest.h b/Tests/UnitTests/TestCore/BeamTest.h
index fc111aef4ed282b0ab153da71547fc5b0f60a6c8..58ddb1dd3dffaf0bc88f0aaec7127477992527a0 100644
--- a/Tests/UnitTests/TestCore/BeamTest.h
+++ b/Tests/UnitTests/TestCore/BeamTest.h
@@ -4,7 +4,6 @@
 #include "Beam.h"
 #include "Units.h"
 #include "ParameterPool.h"
-
 #include "gtest/gtest.h"
 
 
@@ -30,7 +29,7 @@ BeamTest::~BeamTest()
 
 TEST_F(BeamTest, BeamInitialState)
 {
-    EXPECT_DOUBLE_EQ(2.0*M_PI, emptyBeam.getCentralK()[0].real());
+    EXPECT_DOUBLE_EQ(2.0*Units::PI, emptyBeam.getCentralK()[0].real());
     EXPECT_DOUBLE_EQ(0.0, emptyBeam.getCentralK()[0].imag());
     EXPECT_EQ(complex_t(0.0,0.0), emptyBeam.getCentralK()[1]);
     EXPECT_EQ(complex_t(0.0,0.0), emptyBeam.getCentralK()[2]);
diff --git a/Tests/UnitTests/TestCore/DistributionsTest.h b/Tests/UnitTests/TestCore/DistributionsTest.h
index 2ec95b8adc083fbd0dd4fee41801c96a6f52fd54..ba71961926b9da68c9ba68d18bd504d5fb6014ef 100644
--- a/Tests/UnitTests/TestCore/DistributionsTest.h
+++ b/Tests/UnitTests/TestCore/DistributionsTest.h
@@ -1,6 +1,7 @@
 #ifndef DISTRIBUTIONSTEST_H
 #define DISTRIBUTIONSTEST_H
 
+#include "Units.h"
 #include "Distributions.h"
 #include "ParameterSample.h"
 #include <cmath>
@@ -102,7 +103,7 @@ TEST_F(DistributionsTest, DistributionLorentzDefaultConstructor)
     EXPECT_EQ(0.0, distr->getMean());
     EXPECT_EQ(1.0, distr->getHWHM());
     EXPECT_EQ("DistributionLorentz", distr->getName());
-    EXPECT_EQ(1/(2*M_PI), distr->probabilityDensity(1.0));
+    EXPECT_EQ(1/(2*Units::PI), distr->probabilityDensity(1.0));
 
     std::vector<double> list1 = distr->generateValueList(1, 0.0);
     EXPECT_EQ(distr->getMean(), list1[0]);
@@ -126,7 +127,7 @@ TEST_F(DistributionsTest, DistributionLorentzConstructor)
     EXPECT_EQ(1.0, distr2.getMean());
     EXPECT_EQ(1.0, distr2.getHWHM());
     EXPECT_EQ("DistributionLorentz", distr2.getName());
-    EXPECT_EQ(1.0/M_PI, distr2.probabilityDensity(1.0));
+    EXPECT_EQ(1.0/Units::PI, distr2.probabilityDensity(1.0));
 
     std::vector<double> list2 = distr2.generateValueList(1, 0.0);
     EXPECT_EQ(distr2.getMean(), list2[0]);
@@ -198,7 +199,7 @@ TEST_F(DistributionsTest, DistributionGaussianDefaultConstructor)
     DistributionGaussian* id1D = new DistributionGaussian();
     EXPECT_EQ(0.0, id1D->getMean());
     EXPECT_EQ(1.0, id1D->getStdDev());
-    EXPECT_EQ(std::exp(-1.0/2.0)/std::sqrt(2.0*M_PI), id1D->probabilityDensity(1.0));
+    EXPECT_EQ(std::exp(-1.0/2.0)/std::sqrt(2.0*Units::PI), id1D->probabilityDensity(1.0));
     EXPECT_EQ("DistributionGaussian", id1D->getName());
 
     std::vector<double> list1 = id1D->generateValueList(1, 0.0);
@@ -223,7 +224,7 @@ TEST_F(DistributionsTest, DistributionGaussianConstructor)
     DistributionGaussian distr2(1.0, 1.0);
     EXPECT_EQ(1.0, distr2.getMean());
     EXPECT_EQ(1.0, distr2.getStdDev());
-    EXPECT_EQ(1/std::sqrt(2.0*M_PI), distr2.probabilityDensity(1.0));
+    EXPECT_EQ(1/std::sqrt(2.0*Units::PI), distr2.probabilityDensity(1.0));
     EXPECT_EQ("DistributionGaussian", distr2.getName());
 
     std::vector<double> list2 = distr2.generateValueList(1, 0.0);
@@ -247,7 +248,7 @@ TEST_F(DistributionsTest, DistributionGaussianClone)
     DistributionGaussian* id1DClone = id1D->clone();
     EXPECT_EQ(1.0, id1DClone->getMean());
     EXPECT_EQ(1.0, id1DClone->getStdDev());
-    EXPECT_EQ(1/std::sqrt(2.0*M_PI), id1DClone->probabilityDensity(1.0));
+    EXPECT_EQ(1/std::sqrt(2.0*Units::PI), id1DClone->probabilityDensity(1.0));
     EXPECT_EQ("DistributionGaussian", id1DClone->getName());
 
     std::vector<double> list1 = id1DClone->generateValueList(1, 0.0);
@@ -277,7 +278,7 @@ TEST_F(DistributionsTest, DistributionLogNormalConstructorWithOneParameter)
     EXPECT_EQ(1.0, distr2.getMedian());
     EXPECT_EQ(1.0, distr2.getScalePar());
     EXPECT_EQ(std::exp(0.5), distr2.getMean());
-    EXPECT_EQ(1.0/std::sqrt(2.0*M_PI), distr2.probabilityDensity(1.0));
+    EXPECT_EQ(1.0/std::sqrt(2.0*Units::PI), distr2.probabilityDensity(1.0));
     EXPECT_EQ("DistributionLogNormal", distr2.getName());
 
     std::vector<double> list2 = distr2.generateValueList(1, 0.0);
@@ -294,7 +295,7 @@ TEST_F(DistributionsTest, DistributionLogNormalConstructorWithTwoParameter)
     EXPECT_EQ(1.0, id1D->getMedian());
     EXPECT_EQ(1.0, id1D->getScalePar());
     EXPECT_EQ(std::exp(0.5), id1D->getMean());
-    EXPECT_EQ(1.0/std::sqrt(2.0*M_PI), id1D->probabilityDensity(1.0));
+    EXPECT_EQ(1.0/std::sqrt(2.0*Units::PI), id1D->probabilityDensity(1.0));
     EXPECT_EQ("DistributionLogNormal", id1D->getName());
 
     std::vector<double> list1 = id1D->generateValueList(1, 0.0);
@@ -321,7 +322,7 @@ TEST_F(DistributionsTest, DistributionLogNormalClone)
     EXPECT_EQ(1.0, id1D->getMedian());
     EXPECT_EQ(1.0, id1D->getScalePar());
     EXPECT_EQ(std::exp(0.5), id1D->getMean());
-    EXPECT_EQ(1/std::sqrt(2.0*M_PI), id1D->probabilityDensity(1.0));
+    EXPECT_EQ(1/std::sqrt(2.0*Units::PI), id1D->probabilityDensity(1.0));
     EXPECT_EQ("DistributionLogNormal", id1D->getName());
 
     std::vector<double> list1 = id1D->generateValueList(1.0, 0.0);
@@ -342,7 +343,7 @@ TEST_F(DistributionsTest, DistributionCosineDefaultConstructor)
     DistributionCosine* id1D = new DistributionCosine();
     EXPECT_EQ(0.0, id1D->getMean());
     EXPECT_EQ(1.0, id1D->getSigma());
-    EXPECT_DOUBLE_EQ((1.0+std::cos(1.0))/(2.0*M_PI), id1D->probabilityDensity(1.0));
+    EXPECT_DOUBLE_EQ((1.0+std::cos(1.0))/(2.0*Units::PI), id1D->probabilityDensity(1.0));
     EXPECT_EQ(0, id1D->probabilityDensity(100.0));
     EXPECT_EQ("DistributionCosine", id1D->getName());
 
@@ -350,8 +351,8 @@ TEST_F(DistributionsTest, DistributionCosineDefaultConstructor)
     EXPECT_EQ(id1D->getMean(), list1[0]);
 
     std::vector<double> list2 = id1D->generateValueList(2, 0.0);
-    EXPECT_EQ(-M_PI, list2[0]);
-    EXPECT_EQ(M_PI, list2[1]);
+    EXPECT_EQ(-Units::PI, list2[0]);
+    EXPECT_EQ(Units::PI, list2[1]);
 
     delete id1D;
 }
@@ -368,7 +369,7 @@ TEST_F(DistributionsTest, DistributionCosineConstructor)
     DistributionCosine distr2(1.0,1.0);
     EXPECT_EQ(1.0, distr2.getMean());
     EXPECT_EQ(1.0, distr2.getSigma());
-    EXPECT_EQ(2.0/(2.0*M_PI), distr2.probabilityDensity(1.0));
+    EXPECT_EQ(2.0/(2.0*Units::PI), distr2.probabilityDensity(1.0));
     EXPECT_EQ(0, distr2.probabilityDensity(100.0));
     EXPECT_EQ("DistributionCosine", distr2.getName());
 
@@ -376,8 +377,8 @@ TEST_F(DistributionsTest, DistributionCosineConstructor)
     EXPECT_EQ(distr2.getMean(), list2[0]);
 
     std::vector<double> list3 = distr2.generateValueList(2, 0.0);
-    EXPECT_EQ(1-M_PI, list3[0]);
-    EXPECT_EQ(1+M_PI, list3[1]);
+    EXPECT_EQ(1-Units::PI, list3[0]);
+    EXPECT_EQ(1+Units::PI, list3[1]);
 }
 
 TEST_F(DistributionsTest, DistributionCosineParameters)
@@ -393,7 +394,7 @@ TEST_F(DistributionsTest, DistributionCosineClone)
     DistributionCosine* id1DClone = id1D->clone();
     EXPECT_EQ(1.0, id1DClone->getMean());
     EXPECT_EQ(1.0, id1DClone->getSigma());
-    EXPECT_EQ(2.0/(2.0*M_PI), id1DClone->probabilityDensity(1.0));
+    EXPECT_EQ(2.0/(2.0*Units::PI), id1DClone->probabilityDensity(1.0));
     EXPECT_EQ(0, id1D->probabilityDensity(100.0));
     EXPECT_EQ("DistributionCosine", id1DClone->getName());
 
@@ -401,8 +402,8 @@ TEST_F(DistributionsTest, DistributionCosineClone)
     EXPECT_EQ(id1DClone->getMean(), list1[0]);
 
     std::vector<double> list2 = id1DClone->generateValueList(2, 0.0);
-    EXPECT_EQ(1-M_PI, list2[0]);
-    EXPECT_EQ(1+M_PI, list2[1]);
+    EXPECT_EQ(1-Units::PI, list2[0]);
+    EXPECT_EQ(1+Units::PI, list2[1]);
 
     delete id1D;
     delete id1DClone;
diff --git a/Tests/UnitTests/TestCore/FTDistributionsTest.h b/Tests/UnitTests/TestCore/FTDistributionsTest.h
index 46518b565290a2d8ad4537992bf890523888a0c9..16ec4f586014287852e5e4572395cddab38496b0 100644
--- a/Tests/UnitTests/TestCore/FTDistributionsTest.h
+++ b/Tests/UnitTests/TestCore/FTDistributionsTest.h
@@ -179,7 +179,7 @@ TEST_F(FTDistributionsTest, FTDistribution2DCauchyConstructor)
     IFTDistribution2D * iftd2D = new FTDistribution2DCauchy(1.0,-2.0);
     EXPECT_EQ(1.0, iftd2D->getCoherenceLengthX());
     EXPECT_EQ(-2.0, iftd2D->getCoherenceLengthY());
-    EXPECT_EQ(M_PI/2.0, iftd2D->getDelta());
+    EXPECT_EQ(Units::PI/2.0, iftd2D->getDelta());
     EXPECT_EQ(0.0, iftd2D->getGamma());
     EXPECT_EQ("2DDistributionCauchy", iftd2D->getName());
     EXPECT_NEAR(0.343206, iftd2D->evaluate(0.2, 0.5),0.000001);
@@ -202,7 +202,7 @@ TEST_F(FTDistributionsTest, FTDistribution2DCauchyClone)
 
     EXPECT_EQ(-5.0, iftd2DClone->getCoherenceLengthX());
     EXPECT_EQ(2.3, iftd2DClone->getCoherenceLengthY());
-    EXPECT_EQ(M_PI/2.0, iftd2DClone->getDelta());
+    EXPECT_EQ(Units::PI/2.0, iftd2DClone->getDelta());
     EXPECT_EQ(0.0, iftd2DClone->getGamma());
     EXPECT_EQ("2DDistributionCauchy", iftd2DClone->getName());
     EXPECT_NEAR(0.165121078, iftd2DClone->evaluate(0.2, 0.5),0.000001);
@@ -218,7 +218,7 @@ TEST_F(FTDistributionsTest, FTDistribution2DGaussConstructor)
     IFTDistribution2D * iftd2D = new FTDistribution2DGauss(1.0,-2.0);
     EXPECT_EQ(1.0, iftd2D->getCoherenceLengthX());
     EXPECT_EQ(-2.0, iftd2D->getCoherenceLengthY());
-    EXPECT_EQ(M_PI/2.0, iftd2D->getDelta());
+    EXPECT_EQ(Units::PI/2.0, iftd2D->getDelta());
     EXPECT_EQ(0.0, iftd2D->getGamma());
     EXPECT_EQ("2DDistributionGauss", iftd2D->getName());
     EXPECT_NEAR(0.5945205, iftd2D->evaluate(0.2, 0.5),0.000001);
@@ -236,7 +236,7 @@ TEST_F(FTDistributionsTest, FTDistribution2DGaussClone)
 
     EXPECT_EQ(-5.0, iftd2DClone->getCoherenceLengthX());
     EXPECT_EQ(2.3, iftd2DClone->getCoherenceLengthY());
-    EXPECT_EQ(M_PI/2.0, iftd2DClone->getDelta());
+    EXPECT_EQ(Units::PI/2.0, iftd2DClone->getDelta());
     EXPECT_EQ(0.0, iftd2DClone->getGamma());
     EXPECT_EQ("2DDistributionGauss", iftd2DClone->getName());
     EXPECT_NEAR(0.3130945, iftd2DClone->evaluate(0.2, 0.5),0.000001);
@@ -250,7 +250,7 @@ TEST_F(FTDistributionsTest, FTDistribution2DGateConstructor)
     IFTDistribution2D * iftd2D = new FTDistribution2DGate(1.0,-2.0);
     EXPECT_EQ(1.0, iftd2D->getCoherenceLengthX());
     EXPECT_EQ(-2.0, iftd2D->getCoherenceLengthY());
-    EXPECT_EQ(M_PI/2.0, iftd2D->getDelta());
+    EXPECT_EQ(Units::PI/2.0, iftd2D->getDelta());
     EXPECT_EQ(0.0, iftd2D->getGamma());
     EXPECT_EQ("2DDistributionGate", iftd2D->getName());
     EXPECT_NEAR(0.875513, iftd2D->evaluate(0.2, 0.5),0.000001);
@@ -268,7 +268,7 @@ TEST_F(FTDistributionsTest, FTDistribution2DGateClone)
 
     EXPECT_EQ(-5.0, iftd2DClone->getCoherenceLengthX());
     EXPECT_EQ(2.3, iftd2DClone->getCoherenceLengthY());
-    EXPECT_EQ(M_PI/2.0, iftd2DClone->getDelta());
+    EXPECT_EQ(Units::PI/2.0, iftd2DClone->getDelta());
     EXPECT_EQ(0.0, iftd2DClone->getGamma());
     EXPECT_EQ("2DDistributionGate", iftd2DClone->getName());
     EXPECT_NEAR(0.736461, iftd2DClone->evaluate(0.2, 0.5),0.000001);
@@ -282,7 +282,7 @@ TEST_F(FTDistributionsTest, FTDistribution2DConeConstructor)
     IFTDistribution2D * iftd2D = new FTDistribution2DCone(1.0,-2.0);
     EXPECT_EQ(1.0, iftd2D->getCoherenceLengthX());
     EXPECT_EQ(-2.0, iftd2D->getCoherenceLengthY());
-    EXPECT_EQ(M_PI/2.0, iftd2D->getDelta());
+    EXPECT_EQ(Units::PI/2.0, iftd2D->getDelta());
     EXPECT_EQ(0.0, iftd2D->getGamma());
     EXPECT_EQ("2DDistributionCone", iftd2D->getName());
     EXPECT_NEAR(0.924374, iftd2D->evaluate(0.2, 0.5),0.000001);
@@ -300,7 +300,7 @@ TEST_F(FTDistributionsTest, FTDistribution2DConeClone)
 
     EXPECT_EQ(-5.0, iftd2DClone->getCoherenceLengthX());
     EXPECT_EQ(2.3, iftd2DClone->getCoherenceLengthY());
-    EXPECT_EQ(M_PI/2.0, iftd2DClone->getDelta());
+    EXPECT_EQ(Units::PI/2.0, iftd2DClone->getDelta());
     EXPECT_EQ(0.0, iftd2DClone->getGamma());
     EXPECT_EQ("2DDistributionCone", iftd2DClone->getName());
     EXPECT_NEAR(0.837410, iftd2DClone->evaluate(0.2, 0.5),0.000001);
@@ -314,7 +314,7 @@ TEST_F(FTDistributionsTest, FTDistribution2DVoigtConstructor)
     IFTDistribution2D * iftd2D = new FTDistribution2DVoigt(1.0,-2.0,3.5);
     EXPECT_EQ(1.0, iftd2D->getCoherenceLengthX());
     EXPECT_EQ(-2.0, iftd2D->getCoherenceLengthY());
-    EXPECT_EQ(M_PI/2.0, iftd2D->getDelta());
+    EXPECT_EQ(Units::PI/2.0, iftd2D->getDelta());
     EXPECT_EQ(0.0, iftd2D->getGamma());
     EXPECT_EQ("2DDistributionVoigt", iftd2D->getName());
     EXPECT_NEAR(1.2228072, iftd2D->evaluate(0.2, 0.5),0.000001);
@@ -332,7 +332,7 @@ TEST_F(FTDistributionsTest, FTDistribution2DVoigtClone)
 
     EXPECT_EQ(-5.0, iftd2DClone->getCoherenceLengthX());
     EXPECT_EQ(2.3, iftd2DClone->getCoherenceLengthY());
-    EXPECT_EQ(M_PI/2.0, iftd2DClone->getDelta());
+    EXPECT_EQ(Units::PI/2.0, iftd2DClone->getDelta());
     EXPECT_EQ(0.0, iftd2DClone->getGamma());
     EXPECT_EQ("2DDistributionVoigt", iftd2DClone->getName());
     EXPECT_NEAR(-0.6635305, iftd2DClone->evaluate(0.2, 0.5),0.000001);
diff --git a/Tests/UnitTests/TestCore/FormFactorTest.h b/Tests/UnitTests/TestCore/FormFactorTest.h
index b2d115dc6c27ac649f1a27fd4e68980068b40b63..1d8480f762fc45f00c0bfe3964f45970c486854e 100644
--- a/Tests/UnitTests/TestCore/FormFactorTest.h
+++ b/Tests/UnitTests/TestCore/FormFactorTest.h
@@ -16,30 +16,8 @@
 #ifndef FORMFACTORTEST_H
 #define FORMFACTORTEST_H
 
-#include "FormFactorAnisoPyramid.h"
-#include "FormFactorBox.h"
-#include "FormFactorCone.h"
-#include "FormFactorCone6.h"
-#include "FormFactorCrystal.h"
-#include "FormFactorCylinder.h"
-#include "FormFactorCuboctahedron.h"
-#include "FormFactorEllipsoidalCylinder.h"
-#include "FormFactorFullSphere.h"
-#include "FormFactorFullSpheroid.h"
-#include "FormFactorHemiEllipsoid.h"
-#include "FormFactorPrism3.h"
-#include "FormFactorPrism6.h"
-#include "FormFactorPyramid.h"
-#include "FormFactorTruncatedSphere.h"
-#include "FormFactorTruncatedSpheroid.h"
-#include "FormFactorTetrahedron.h"
-#include "FormFactorRipple2.h"
-#include "FormFactorRipple1.h"
-#include "FormFactorInfLongBox.h"
-#include "FormFactorInfLongRipple2.h"
-#include "FormFactorInfLongRipple1.h"
-#include "FormFactorTruncatedCube.h"
-
+#include "Units.h"
+#include "FormFactors.h"
 #include "gtest/gtest.h"
 
 class FormFactorTest : public ::testing::Test
@@ -89,7 +67,7 @@ TEST_F(FormFactorTest, HemiEllipsoid)
     double radius_b = 7.;
     double height = 5.;
 
-    double volume = 2.*M_PI*radius_a*radius_b*height/3.;
+    double volume = 2.*Units::PI*radius_a*radius_b*height/3.;
 
     FormFactorHemiEllipsoid hemiellipsoid(radius_a, radius_b, height);
 
@@ -143,7 +121,7 @@ TEST_F(FormFactorTest, Cone)
     double alpha = 0.8;
     double tga = std::tan(alpha);
     double HdivRtga = height/tga/radius;
-    double volume = M_PI/3.*tga*radius*radius*radius*
+    double volume = Units::PI/3.*tga*radius*radius*radius*
             (1. - (1.- HdivRtga)*(1.- HdivRtga)*(1.- HdivRtga));
 
     FormFactorCone cone(radius, height, alpha);
@@ -235,7 +213,7 @@ TEST_F(FormFactorTest, Cylinder)
 {
     double radius = 3.;
     double height = 5.;
-    double volume = M_PI*radius*radius*height;
+    double volume = Units::PI*radius*radius*height;
 
     FormFactorCylinder cylinder(radius,height);
     EXPECT_EQ("FormFactorCylinder",cylinder.getName());
@@ -258,7 +236,7 @@ TEST_F(FormFactorTest, EllipsoidalCylinder)
     double radius_a = 3.;
     double radius_b = 5.;
     double height = 4;
-    double volume = M_PI*radius_a*radius_b*height;
+    double volume = Units::PI*radius_a*radius_b*height;
 
     FormFactorEllipsoidalCylinder ellipscyl(radius_a, radius_b, height);
 
@@ -283,7 +261,7 @@ TEST_F(FormFactorTest, EllipsoidalCylinder)
 TEST_F(FormFactorTest, FullSphere)
 {
     double radius = 5.;
-    double volume = 4./3.*M_PI*radius*radius*radius;
+    double volume = 4./3.*Units::PI*radius*radius*radius;
 
     FormFactorFullSphere fullsphere(radius);
 
@@ -306,7 +284,7 @@ TEST_F(FormFactorTest, FullSpheroid)
 {
     double radius = 3.;
     double height = 5.;
-    double volume = 2./3.*M_PI*radius*radius*height;
+    double volume = 2./3.*Units::PI*radius*radius*height;
 
     FormFactorFullSpheroid fullspheroid(radius,height);
 
@@ -407,7 +385,7 @@ TEST_F(FormFactorTest, TruncatedSphere)
     double radius = 5.;
     double height = 3.;
     double HdivR = height/radius;
-    double volume = M_PI/3.*radius*radius*radius
+    double volume = Units::PI/3.*radius*radius*radius
             *(3.*HdivR -1. - (HdivR - 1.)*(HdivR - 1.)*(HdivR - 1.));
 
     FormFactorTruncatedSphere trsphere(radius, height);
@@ -430,7 +408,7 @@ TEST_F(FormFactorTest, TruncatedSpheroid)
     double radius = 3.;
     double flattening = 1.5;
     double total_height =2.*flattening *radius;
-    double volume = M_PI*radius*height*height/flattening
+    double volume = Units::PI*radius*height*height/flattening
             *(1.-height/(3.*flattening*radius));
 
     FormFactorTruncatedSpheroid trspheroid(radius, height,flattening);
@@ -532,73 +510,73 @@ TEST_F(FormFactorTest, Ripple1)
 TEST_F(FormFactorTest, InfLongBox)
 {
     double height = 15.;
-    double width = 100.0/M_PI;
+    double width = 100.0/Units::PI;
 
     FormFactorInfLongBox ilbox(width, height);
 
     EXPECT_EQ("FormFactorInfLongBox",ilbox.getName());
-    EXPECT_DOUBLE_EQ(100./M_PI, ilbox.getWidth());
+    EXPECT_DOUBLE_EQ(100./Units::PI, ilbox.getWidth());
     EXPECT_EQ(15., ilbox.getHeight());
-    EXPECT_DOUBLE_EQ(50./M_PI, ilbox.getRadius());
+    EXPECT_DOUBLE_EQ(50./Units::PI, ilbox.getRadius());
     EXPECT_THROW(ilbox.getVolume(), NotImplementedException);
     EXPECT_EQ(2, ilbox.getNumberOfStochasticParameters());
 
     FormFactorInfLongBox *ilboxclone = ilbox.clone();
     EXPECT_EQ("FormFactorInfLongBox",ilboxclone->getName());
-    EXPECT_DOUBLE_EQ(100./M_PI, ilboxclone->getWidth());
+    EXPECT_DOUBLE_EQ(100./Units::PI, ilboxclone->getWidth());
     EXPECT_EQ(15., ilboxclone->getHeight());
     EXPECT_THROW(ilboxclone->getVolume(),NotImplementedException);
     EXPECT_EQ(2, ilboxclone->getNumberOfStochasticParameters());
-    EXPECT_DOUBLE_EQ(50./M_PI, ilboxclone->getRadius());
+    EXPECT_DOUBLE_EQ(50./Units::PI, ilboxclone->getRadius());
 }
 
 // Test form factor of a infinite long ripple1
 TEST_F(FormFactorTest, InfLongRipple1)
 {
     double height = 15.;
-    double width = 100.0/M_PI;
+    double width = 100.0/Units::PI;
 
     FormFactorInfLongRipple1 ilripple1(width, height);
 
     EXPECT_EQ("FormFactorInfLongRipple1",ilripple1.getName());
-    EXPECT_DOUBLE_EQ(100./M_PI, ilripple1.getWidth());
+    EXPECT_DOUBLE_EQ(100./Units::PI, ilripple1.getWidth());
     EXPECT_EQ(15., ilripple1.getHeight());
-    EXPECT_DOUBLE_EQ(50./M_PI, ilripple1.getRadius());
+    EXPECT_DOUBLE_EQ(50./Units::PI, ilripple1.getRadius());
     EXPECT_THROW(ilripple1.getVolume(), NotImplementedException);
     EXPECT_EQ(2, ilripple1.getNumberOfStochasticParameters());
 
     FormFactorInfLongRipple1 *ilripple1clone = ilripple1.clone();
     EXPECT_EQ("FormFactorInfLongRipple1",ilripple1clone->getName());
-    EXPECT_DOUBLE_EQ(100./M_PI, ilripple1clone->getWidth());
+    EXPECT_DOUBLE_EQ(100./Units::PI, ilripple1clone->getWidth());
     EXPECT_EQ(15., ilripple1clone->getHeight());
     EXPECT_THROW(ilripple1clone->getVolume(), NotImplementedException);
     EXPECT_EQ(2, ilripple1clone->getNumberOfStochasticParameters());
-    EXPECT_DOUBLE_EQ(50./M_PI, ilripple1clone->getRadius());
+    EXPECT_DOUBLE_EQ(50./Units::PI, ilripple1clone->getRadius());
 }
 
 // Test form factor of a infinite long ripple2
 TEST_F(FormFactorTest, InfLongRipple2)
 {
     double height = 15.;
-    double width = 100.0/M_PI;
+    double width = 100.0/Units::PI;
     double d = 0.3; // asymetry
 
     FormFactorInfLongRipple2 ilripple2(width, height, d);
 
     EXPECT_EQ("FormFactorInfLongRipple2",ilripple2.getName());
-    EXPECT_DOUBLE_EQ(100./M_PI, ilripple2.getWidth());
+    EXPECT_DOUBLE_EQ(100./Units::PI, ilripple2.getWidth());
     EXPECT_EQ(15., ilripple2.getHeight());
-    EXPECT_DOUBLE_EQ(50./M_PI, ilripple2.getRadius());
+    EXPECT_DOUBLE_EQ(50./Units::PI, ilripple2.getRadius());
     EXPECT_THROW(ilripple2.getVolume(), NotImplementedException);
     EXPECT_EQ(3, ilripple2.getNumberOfStochasticParameters());
 
     FormFactorInfLongRipple2 *ilripple2clone = ilripple2.clone();
     EXPECT_EQ("FormFactorInfLongRipple2",ilripple2clone->getName());
-    EXPECT_DOUBLE_EQ(100./M_PI, ilripple2clone->getWidth());
+    EXPECT_DOUBLE_EQ(100./Units::PI, ilripple2clone->getWidth());
     EXPECT_EQ(15., ilripple2clone->getHeight());
     EXPECT_THROW(ilripple2clone->getVolume(), NotImplementedException);
     EXPECT_EQ(3, ilripple2clone->getNumberOfStochasticParameters());
-    EXPECT_DOUBLE_EQ(50./M_PI, ilripple2clone->getRadius());
+    EXPECT_DOUBLE_EQ(50./Units::PI, ilripple2clone->getRadius());
 }
 
 // Test form factor of a truncated cube
diff --git a/Tests/UnitTests/TestCore/InstrumentTest.h b/Tests/UnitTests/TestCore/InstrumentTest.h
index 387f5eacdaaaea7b798c4a833124e998ade0fba2..86719f8818b63b2cc34470f3315822f60a296756 100644
--- a/Tests/UnitTests/TestCore/InstrumentTest.h
+++ b/Tests/UnitTests/TestCore/InstrumentTest.h
@@ -1,6 +1,7 @@
 #ifndef INSTRUMENTTEST_H_
 #define INSTRUMENTTEST_H_
 
+#include "Units.h"
 #include "Instrument.h"
 #include "BornAgainNamespace.h"
 
@@ -34,7 +35,7 @@ TEST_F(InstrumentTest, InstrumentInitialState)
 TEST_F(InstrumentTest, BeamManipulation)
 {
     double lambda(1), alpha(-1), phi(1);
-    double k = 2.*M_PI/lambda;
+    double k = 2.*Units::PI/lambda;
     double x = k*std::cos(alpha) * std::cos(phi);
     double y = k*std::cos(alpha) * std::sin(phi);
     double z = k*std::sin(alpha);
diff --git a/Tests/UnitTests/TestCore/KVectorTest.h b/Tests/UnitTests/TestCore/KVectorTest.h
index f0f810d182fcad5e840f0dca0d8ffc4b883bf0ba..f2aa5a1a672e34ac30e18c961c7c91693adfbe3d 100644
--- a/Tests/UnitTests/TestCore/KVectorTest.h
+++ b/Tests/UnitTests/TestCore/KVectorTest.h
@@ -1,9 +1,9 @@
 #ifndef KVECTORTEST_H
 #define KVECTORTEST_H
 
+#include "Units.h"
 #include "Types.h"
 
-
 class KVectorTest : public ::testing::Test
 {
  protected:
@@ -109,14 +109,14 @@ TEST_F(KVectorTest, BasicTransformation)
 
     // rotation via transformation
     a = kvector_t(std::sqrt(3.)/2., 2., 0.5);
-    Geometry::Transform3D m2 = Geometry::Transform3D::createRotateY(M_PI/6.);
+    Geometry::Transform3D m2 = Geometry::Transform3D::createRotateY(Units::PI/6.);
     v = m2.transformed(a);
     ASSERT_NEAR(      v.x(), 1.0, epsilon );
     EXPECT_DOUBLE_EQ( v.y(), 2.0 );
     ASSERT_NEAR(      v.z(), 0.0, epsilon );
 
     a = kvector_t(0.5, std::sqrt(3.)/2., 2.);
-    Geometry::Transform3D m3 = Geometry::Transform3D::createRotateZ(M_PI/6.);
+    Geometry::Transform3D m3 = Geometry::Transform3D::createRotateZ(Units::PI/6.);
     v = m3.transformed(a);
     ASSERT_NEAR(      v.x(), 0.0, epsilon );
     ASSERT_NEAR(      v.y(), 1.0, epsilon );
diff --git a/Tests/UnitTests/TestCore/MatrixSpecularInfoMapTest.h b/Tests/UnitTests/TestCore/MatrixSpecularInfoMapTest.h
index c42b17793cab85bc5d4044eddd3ff1c65a37c758..397eefb52c1ef356b1f246f93db0c2f584eb017e 100644
--- a/Tests/UnitTests/TestCore/MatrixSpecularInfoMapTest.h
+++ b/Tests/UnitTests/TestCore/MatrixSpecularInfoMapTest.h
@@ -1,6 +1,7 @@
 #ifndef MATRIXSPECULARINFOMAPTEST_H
 #define MATRIXSPECULARINFOMAPTEST_H
 
+#include "Units.h"
 #include "MatrixRTCoefficients.h"
 #include "MatrixSpecularInfoMap.h"
 #include "gtest/gtest.h"
@@ -34,7 +35,7 @@ MatrixSpecularInfoMapTest::MatrixSpecularInfoMapTest()
 
 TEST_F(MatrixSpecularInfoMapTest, getCoefficients)
 {
-    MatrixSpecularInfoMap map(mp_multilayer, 0, 2.0*M_PI);
+    MatrixSpecularInfoMap map(mp_multilayer, 0, 2.0*Units::PI);
     boost::scoped_ptr<const MatrixRTCoefficients> rt_coeffs(
                 map.getCoefficients(1.0, 1.0));
     complex_t R0 = complex_t(0.1750375, -0.0222467);
diff --git a/Tests/UnitTests/TestCore/ScalarSpecularInfoMapTest.h b/Tests/UnitTests/TestCore/ScalarSpecularInfoMapTest.h
index d0365d4cc9c1b4de6f66ec3a64244f89f397d2a6..268550371240665fcc24df16188f2676ac4e1179 100644
--- a/Tests/UnitTests/TestCore/ScalarSpecularInfoMapTest.h
+++ b/Tests/UnitTests/TestCore/ScalarSpecularInfoMapTest.h
@@ -1,6 +1,7 @@
 #ifndef SCALARSPECULARINFOMAPTEST_H
 #define SCALARSPECULARINFOMAPTEST_H
 
+#include "Units.h"
 #include "ScalarRTCoefficients.h"
 #include "ScalarSpecularInfoMap.h"
 #include "gtest/gtest.h"
@@ -35,7 +36,7 @@ ScalarSpecularInfoMapTest::ScalarSpecularInfoMapTest()
 
 TEST_F(ScalarSpecularInfoMapTest, getCoefficients)
 {
-    ScalarSpecularInfoMap map(mp_multilayer, 0, 2.0*M_PI);
+    ScalarSpecularInfoMap map(mp_multilayer, 0, 2.0*Units::PI);
     boost::scoped_ptr<const ScalarRTCoefficients> rt_coeffs(
                 map.getCoefficients(1.0, 1.0));
     complex_t R0 = complex_t(0.1750375, -0.0222467);