diff --git a/addons/gm2calc/MSSMNoFV_onshell.cpp b/addons/gm2calc/MSSMNoFV_onshell.cpp
index 9876a8b58..93cd7a993 100644
--- a/addons/gm2calc/MSSMNoFV_onshell.cpp
+++ b/addons/gm2calc/MSSMNoFV_onshell.cpp
@@ -25,6 +25,7 @@
 #include <cmath>
 #include <complex>
 #include <iostream>
+#include <sstream>
 
 #define WARNING(message) std::cerr << "Warning: " << message << '\n';
 
@@ -45,7 +46,6 @@ namespace {
    }
 }
 
-namespace flexiblesusy {
 namespace gm2calc {
 
 MSSMNoFV_onshell::MSSMNoFV_onshell()
@@ -82,6 +82,15 @@ void MSSMNoFV_onshell::set_alpha_thompson(double alpha)
    EL0 = calculate_e(alpha);
 }
 
+void MSSMNoFV_onshell::set_TB(double tanb)
+{
+   const double vev = std::sqrt(sqr(vu) + sqr(vd));
+   const double sinb = tanb / std::sqrt(1 + tanb*tanb);
+   const double cosb = 1.   / std::sqrt(1 + tanb*tanb);
+   set_vd(vev * cosb);
+   set_vu(vev * sinb);
+}
+
 /**
  * Returns the electromagnetig gauge coupling in the Thompson limit.
  */
@@ -110,7 +119,11 @@ void MSSMNoFV_onshell::convert_to_onshell(double precision) {
    convert_yukawa_couplings(); // first guess of resummed yukawas
    convert_mf2(precision, 1000);
    convert_yukawa_couplings();
+
+   // final mass spectrum
+   get_problems().clear();
    calculate_DRbar_masses();
+   check_problems();
 }
 
 /**
@@ -127,7 +140,11 @@ void MSSMNoFV_onshell::calculate_masses() {
    convert_vev();
    convert_yukawa_couplings_treelevel();
    convert_yukawa_couplings(); // tan(beta) resummation in Yukawas
-   calculate_DRbar_masses();   // final mass spectrum
+
+   // final mass spectrum
+   get_problems().clear();
+   calculate_DRbar_masses();
+   check_problems();
 
    // copy SUSY masses to physical struct
    copy_susy_masses_to_pole();
@@ -143,6 +160,21 @@ void MSSMNoFV_onshell::check_input()
       throw EInvalidInput("Z mass is zero");
    if (is_zero(get_MM()))
       throw EInvalidInput("Muon mass is zero");
+   if (is_zero(get_MassB()))
+      throw EInvalidInput("Bino mass M1 is zero");
+   if (is_zero(get_MassWB()))
+      throw EInvalidInput("Bino mass M2 is zero");
+   if (is_zero(get_MassG()))
+      throw EInvalidInput("Gluino mass M3 is zero");
+}
+
+void MSSMNoFV_onshell::check_problems()
+{
+   if (get_problems().have_problem()) {
+      std::ostringstream sstr;
+      sstr << get_problems();
+      throw EPhysicalProblem(sstr.str());
+   }
 }
 
 void MSSMNoFV_onshell::copy_susy_masses_to_pole()
@@ -461,4 +493,3 @@ std::ostream& operator<<(std::ostream& os, const MSSMNoFV_onshell& model)
 }
 
 } // gm2calc
-} // namespace flexiblesusy
diff --git a/addons/gm2calc/MSSMNoFV_onshell.hpp b/addons/gm2calc/MSSMNoFV_onshell.hpp
index 9789d2423..98b796e5e 100644
--- a/addons/gm2calc/MSSMNoFV_onshell.hpp
+++ b/addons/gm2calc/MSSMNoFV_onshell.hpp
@@ -23,7 +23,6 @@
 #include <limits>
 #include <Eigen/Core>
 
-namespace flexiblesusy {
 namespace gm2calc {
 
 class MSSMNoFV_onshell : public MSSMNoFV_onshell_mass_eigenstates {
@@ -44,6 +43,7 @@ class MSSMNoFV_onshell : public MSSMNoFV_onshell_mass_eigenstates {
    void set_Au(unsigned i, unsigned k, double a) { Au(i,k) = a; }
    void set_Ad(unsigned i, unsigned k, double a) { Ad(i,k) = a; }
    void set_MA0(double m) { get_physical().MAh(1) = m; }
+   void set_TB(double);
 
    double get_MUDIM() const {return get_scale();}
    double get_EL0() const {return EL0;}
@@ -105,6 +105,7 @@ class MSSMNoFV_onshell : public MSSMNoFV_onshell_mass_eigenstates {
    unsigned find_bino_like_neutralino();
 
    void check_input();
+   void check_problems();
    void convert_gauge_couplings();
    void convert_BMu();
    void convert_mf2(double, unsigned);
@@ -115,6 +116,5 @@ class MSSMNoFV_onshell : public MSSMNoFV_onshell_mass_eigenstates {
 };
 
 } // namespace gm2calc
-} // namespace flexiblesusy
 
 #endif
diff --git a/addons/gm2calc/MSSMNoFV_onshell_mass_eigenstates.cpp b/addons/gm2calc/MSSMNoFV_onshell_mass_eigenstates.cpp
index c032e3e15..4bee0185f 100644
--- a/addons/gm2calc/MSSMNoFV_onshell_mass_eigenstates.cpp
+++ b/addons/gm2calc/MSSMNoFV_onshell_mass_eigenstates.cpp
@@ -32,7 +32,6 @@
 #include "eigen_utils.hpp"
 #include "linalg2.hpp"
 #include "numerics2.hpp"
-#include "error.hpp"
 #include "ffunctions.hpp"
 
 #include <cmath>
@@ -65,9 +64,9 @@ void Hermitianize(Eigen::MatrixBase<Derived>& m)
 }
 } // anonymous namespace
 
-namespace flexiblesusy {
+namespace gm2calc {
 
-using namespace gm2calc;
+using namespace flexiblesusy;
 
 #define CLASSNAME MSSMNoFV_onshell_mass_eigenstates
 #define PHYSICAL(parameter) physical.parameter
@@ -78,6 +77,7 @@ CLASSNAME::MSSMNoFV_onshell_mass_eigenstates()
    , force_output(false)
    , precision(1.0e-3)
    , physical()
+   , problems()
    , MVG(0), MGlu(0), MVP(0), MVZ(0), MFd(0), MFs(0), MFb(0), MFu(0), MFc(0),
       MFt(0), MFve(0), MFvm(0), MFvt(0), MFe(0), MFm(0), MFtau(0), MSveL(0), MSvmL
       (0), MSvtL(0), MSd(Eigen::Array<double,2,1>::Zero()), MSu(Eigen::Array<
@@ -133,6 +133,16 @@ void CLASSNAME::set_physical(const MSSMNoFV_onshell_physical& physical_)
    physical = physical_;
 }
 
+const MSSMNoFV_onshell_problems& CLASSNAME::get_problems() const
+{
+   return problems;
+}
+
+MSSMNoFV_onshell_problems& CLASSNAME::get_problems()
+{
+   return problems;
+}
+
 int CLASSNAME::solve_ewsb_tree_level()
 {
    return solve_ewsb_tree_level_via_soft_higgs_masses();
@@ -232,6 +242,7 @@ void CLASSNAME::print(std::ostream& ostr) const
    ostr << "UP = " << UP << '\n';
 
    physical.print(ostr);
+   problems.print(ostr);
 }
 
 /**
@@ -423,6 +434,7 @@ void CLASSNAME::clear()
    MSSMNoFV_onshell_soft_parameters::clear();
    clear_DRbar_parameters();
    physical.clear();
+   problems.clear();
 }
 
 std::string CLASSNAME::name() const
@@ -670,7 +682,7 @@ void CLASSNAME::calculate_MSveL()
    MSveL = std::abs(mass_matrix_SveL);
 
    // if (MSveL < 0.)
-   //    problems.flag_tachyon(MSSMNoFV_onshell_info::SveL);
+   //    problems.flag_tachyon("SveL");
 
    MSveL = std::sqrt(std::abs(MSveL));
 }
@@ -688,8 +700,8 @@ void CLASSNAME::calculate_MSvmL()
    const auto mass_matrix_SvmL = get_mass_matrix_SvmL();
    MSvmL = std::abs(mass_matrix_SvmL);
 
-   // if (MSvmL < 0.)
-   //    problems.flag_tachyon(MSSMNoFV_onshell_info::SvmL);
+   if (MSvmL < 0.)
+      problems.flag_tachyon("SvmL");
 
    MSvmL = std::sqrt(std::abs(MSvmL));
 }
@@ -708,7 +720,7 @@ void CLASSNAME::calculate_MSvtL()
    MSvtL = std::abs(mass_matrix_SvtL);
 
    // if (MSvtL < 0.)
-   //    problems.flag_tachyon(MSSMNoFV_onshell_info::SvtL);
+   //    problems.flag_tachyon("SvtL");
 
    MSvtL = std::sqrt(std::abs(MSvtL));
 }
@@ -736,7 +748,7 @@ void CLASSNAME::calculate_MSd()
    fs_diagonalize_hermitian(mass_matrix_Sd, MSd, ZD);
 
    // if (MSd.minCoeff() < 0.)
-   //    problems.flag_tachyon(MSSMNoFV_onshell_info::Sd);
+   //    problems.flag_tachyon("Sd");
 
    MSd = sqrt(MSd.cwiseAbs());
 }
@@ -764,7 +776,7 @@ void CLASSNAME::calculate_MSu()
    fs_diagonalize_hermitian(mass_matrix_Su, MSu, ZU);
 
    // if (MSu.minCoeff() < 0.)
-   //    problems.flag_tachyon(MSSMNoFV_onshell_info::Su);
+   //    problems.flag_tachyon("Su");
 
    MSu = sqrt(MSu.cwiseAbs());
 }
@@ -792,7 +804,7 @@ void CLASSNAME::calculate_MSe()
    fs_diagonalize_hermitian(mass_matrix_Se, MSe, ZE);
 
    // if (MSe.minCoeff() < 0.)
-   //    problems.flag_tachyon(MSSMNoFV_onshell_info::Se);
+   //    problems.flag_tachyon("Se");
 
    MSe = sqrt(MSe.cwiseAbs());
 }
@@ -819,8 +831,8 @@ void CLASSNAME::calculate_MSm()
    const auto mass_matrix_Sm(get_mass_matrix_Sm());
    fs_diagonalize_hermitian(mass_matrix_Sm, MSm, ZM);
 
-   // if (MSm.minCoeff() < 0.)
-   //    problems.flag_tachyon(MSSMNoFV_onshell_info::Sm);
+   if (MSm.minCoeff() < 0.)
+      problems.flag_tachyon("Sm");
 
    MSm = sqrt(MSm.cwiseAbs());
 }
@@ -847,8 +859,8 @@ void CLASSNAME::calculate_MStau()
    const auto mass_matrix_Stau(get_mass_matrix_Stau());
    fs_diagonalize_hermitian(mass_matrix_Stau, MStau, ZTau);
 
-   // if (MStau.minCoeff() < 0.)
-   //    problems.flag_tachyon(MSSMNoFV_onshell_info::Stau);
+   if (MStau.minCoeff() < 0.)
+      problems.flag_tachyon("Stau");
 
    MStau = sqrt(MStau.cwiseAbs());
 }
@@ -876,7 +888,7 @@ void CLASSNAME::calculate_MSs()
    fs_diagonalize_hermitian(mass_matrix_Ss, MSs, ZS);
 
    // if (MSs.minCoeff() < 0.)
-   //    problems.flag_tachyon(MSSMNoFV_onshell_info::Ss);
+   //    problems.flag_tachyon("Ss");
 
    MSs = sqrt(MSs.cwiseAbs());
 }
@@ -904,7 +916,7 @@ void CLASSNAME::calculate_MSc()
    fs_diagonalize_hermitian(mass_matrix_Sc, MSc, ZC);
 
    // if (MSc.minCoeff() < 0.)
-   //    problems.flag_tachyon(MSSMNoFV_onshell_info::Sc);
+   //    problems.flag_tachyon("Sc");
 
    MSc = sqrt(MSc.cwiseAbs());
 }
@@ -931,8 +943,8 @@ void CLASSNAME::calculate_MSb()
    const auto mass_matrix_Sb(get_mass_matrix_Sb());
    fs_diagonalize_hermitian(mass_matrix_Sb, MSb, ZB);
 
-   // if (MSb.minCoeff() < 0.)
-   //    problems.flag_tachyon(MSSMNoFV_onshell_info::Sb);
+   if (MSb.minCoeff() < 0.)
+      problems.flag_tachyon("Sb");
 
    MSb = sqrt(MSb.cwiseAbs());
 }
@@ -959,8 +971,8 @@ void CLASSNAME::calculate_MSt()
    const auto mass_matrix_St(get_mass_matrix_St());
    fs_diagonalize_hermitian(mass_matrix_St, MSt, ZT);
 
-   // if (MSt.minCoeff() < 0.)
-   //    problems.flag_tachyon(MSSMNoFV_onshell_info::St);
+   if (MSt.minCoeff() < 0.)
+      problems.flag_tachyon("St");
 
    MSt = sqrt(MSt.cwiseAbs());
 }
@@ -986,8 +998,8 @@ void CLASSNAME::calculate_Mhh()
    const auto mass_matrix_hh(get_mass_matrix_hh());
    fs_diagonalize_hermitian(mass_matrix_hh, Mhh, ZH);
 
-   // if (Mhh.minCoeff() < 0.)
-   //    problems.flag_tachyon(MSSMNoFV_onshell_info::hh);
+   if (Mhh.minCoeff() < 0.)
+      problems.flag_tachyon("hh");
 
    Mhh = sqrt(Mhh.cwiseAbs());
 }
@@ -1018,8 +1030,8 @@ void CLASSNAME::calculate_MAh()
    const auto mass_matrix_Ah(get_mass_matrix_Ah());
    fs_diagonalize_hermitian(mass_matrix_Ah, MAh, ZA);
 
-   // if (MAh.minCoeff() < 0.)
-   //    problems.flag_tachyon(MSSMNoFV_onshell_info::Ah);
+   if (MAh.minCoeff() < 0.)
+      problems.flag_tachyon("Ah");
 
    MAh = sqrt(MAh.cwiseAbs());
 }
@@ -1044,8 +1056,8 @@ void CLASSNAME::calculate_MHpm()
    const auto mass_matrix_Hpm(get_mass_matrix_Hpm());
    fs_diagonalize_hermitian(mass_matrix_Hpm, MHpm, ZP);
 
-   // if (MHpm.minCoeff() < 0.)
-   //    problems.flag_tachyon(MSSMNoFV_onshell_info::Hpm);
+   if (MHpm.minCoeff() < 0.)
+      problems.flag_tachyon("Hpm");
 
    MHpm = sqrt(MHpm.cwiseAbs());
 }
@@ -1106,8 +1118,8 @@ void CLASSNAME::calculate_MVWm()
    const auto mass_matrix_VWm = get_mass_matrix_VWm();
    MVWm = std::abs(mass_matrix_VWm);
 
-   // if (MVWm < 0.)
-   //    problems.flag_tachyon(MSSMNoFV_onshell_info::VWm);
+   if (MVWm < 0.)
+      problems.flag_tachyon("VWm");
 
    MVWm = std::sqrt(std::abs(MVWm));
 }
@@ -1148,4 +1160,4 @@ std::ostream& operator<<(std::ostream& ostr, const MSSMNoFV_onshell_mass_eigenst
    return ostr;
 }
 
-} // namespace flexiblesusy
+} // namespace gm2calc
diff --git a/addons/gm2calc/MSSMNoFV_onshell_mass_eigenstates.hpp b/addons/gm2calc/MSSMNoFV_onshell_mass_eigenstates.hpp
index 44a88db53..27654177b 100644
--- a/addons/gm2calc/MSSMNoFV_onshell_mass_eigenstates.hpp
+++ b/addons/gm2calc/MSSMNoFV_onshell_mass_eigenstates.hpp
@@ -32,13 +32,13 @@
 
 #include "MSSMNoFV_onshell_soft_parameters.hpp"
 #include "MSSMNoFV_onshell_physical.hpp"
+#include "MSSMNoFV_onshell_problems.hpp"
 
 #include <iosfwd>
 #include <string>
 
 #include <Eigen/Core>
 
-namespace flexiblesusy {
 namespace gm2calc {
 
 /**
@@ -64,13 +64,14 @@ class MSSMNoFV_onshell_mass_eigenstates : public MSSMNoFV_onshell_soft_parameter
    void set_physical(const MSSMNoFV_onshell_physical&);
    const MSSMNoFV_onshell_physical& get_physical() const;
    MSSMNoFV_onshell_physical& get_physical();
+   const MSSMNoFV_onshell_problems& get_problems() const;
+   MSSMNoFV_onshell_problems& get_problems();
    int solve_ewsb_tree_level();
    int solve_ewsb();
 
    std::string name() const;
    void print(std::ostream&) const;
 
-
    double get_MVG() const { return MVG; }
    double get_MGlu() const { return MGlu; }
    double get_MVP() const { return MVP; }
@@ -238,6 +239,7 @@ class MSSMNoFV_onshell_mass_eigenstates : public MSSMNoFV_onshell_soft_parameter
    bool force_output;             ///< switch to force output of pole masses
    double precision;              ///< RG running precision
    MSSMNoFV_onshell_physical physical; ///< contains the pole masses and mixings
+   MSSMNoFV_onshell_problems problems; ///< problems
 
    int solve_ewsb_tree_level_via_soft_higgs_masses();
 
@@ -301,6 +303,5 @@ class MSSMNoFV_onshell_mass_eigenstates : public MSSMNoFV_onshell_soft_parameter
 std::ostream& operator<<(std::ostream&, const MSSMNoFV_onshell_mass_eigenstates&);
 
 } // namespace gm2calc
-} // namespace flexiblesusy
 
 #endif
diff --git a/addons/gm2calc/MSSMNoFV_onshell_physical.cpp b/addons/gm2calc/MSSMNoFV_onshell_physical.cpp
index b6aadf22e..83fba708f 100644
--- a/addons/gm2calc/MSSMNoFV_onshell_physical.cpp
+++ b/addons/gm2calc/MSSMNoFV_onshell_physical.cpp
@@ -21,9 +21,10 @@
 
 #include <iostream>
 
-namespace flexiblesusy {
 namespace gm2calc {
 
+using namespace flexiblesusy;
+
 MSSMNoFV_onshell_physical::MSSMNoFV_onshell_physical()
    :
     MVG(0), MGlu(0), MVP(0), MVZ(0), MFd(0), MFs(0), MFb(0), MFu(0), MFc(0),
@@ -228,4 +229,3 @@ std::ostream& operator<<(std::ostream& ostr, const MSSMNoFV_onshell_physical& ph
 }
 
 } // namespace gm2calc
-} // namespace flexiblesusy
diff --git a/addons/gm2calc/MSSMNoFV_onshell_physical.hpp b/addons/gm2calc/MSSMNoFV_onshell_physical.hpp
index 7fb81467c..084e0e15d 100644
--- a/addons/gm2calc/MSSMNoFV_onshell_physical.hpp
+++ b/addons/gm2calc/MSSMNoFV_onshell_physical.hpp
@@ -23,7 +23,6 @@
 
 #include <iosfwd>
 
-namespace flexiblesusy {
 namespace gm2calc {
 
 struct MSSMNoFV_onshell_physical {
@@ -89,6 +88,5 @@ struct MSSMNoFV_onshell_physical {
 std::ostream& operator<<(std::ostream&, const MSSMNoFV_onshell_physical&);
 
 } // namespace gm2calc
-} // namespace flexiblesusy
 
 #endif
diff --git a/addons/gm2calc/MSSMNoFV_onshell_problems.cpp b/addons/gm2calc/MSSMNoFV_onshell_problems.cpp
new file mode 100644
index 000000000..e6f1577ea
--- /dev/null
+++ b/addons/gm2calc/MSSMNoFV_onshell_problems.cpp
@@ -0,0 +1,61 @@
+// ====================================================================
+// This file is part of GM2Calc.
+//
+// GM2Calc is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published
+// by the Free Software Foundation, either version 3 of the License,
+// or (at your option) any later version.
+//
+// GM2Calc is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with GM2Calc.  If not, see
+// <http://www.gnu.org/licenses/>.
+// ====================================================================
+
+#include "MSSMNoFV_onshell_problems.hpp"
+#include "numerics2.hpp"
+
+#include <iostream>
+
+namespace gm2calc {
+
+MSSMNoFV_onshell_problems::MSSMNoFV_onshell_problems()
+   : have_tachyon(false)
+   , tachyonic_particle("")
+{
+}
+
+void MSSMNoFV_onshell_problems::clear()
+{
+   have_tachyon = false;
+   tachyonic_particle.clear();
+}
+
+void MSSMNoFV_onshell_problems::flag_tachyon(const std::string& particle_name)
+{
+   have_tachyon = true;
+   tachyonic_particle = particle_name;
+}
+
+bool MSSMNoFV_onshell_problems::have_problem() const
+{
+   return have_tachyon;
+}
+
+void MSSMNoFV_onshell_problems::print(std::ostream& ostr) const
+{
+   if (have_tachyon)
+      ostr << "Problem: " << tachyonic_particle << " tachyon";
+}
+
+std::ostream& operator<<(std::ostream& ostr, const MSSMNoFV_onshell_problems& problems)
+{
+   problems.print(ostr);
+   return ostr;
+}
+
+} // namespace gm2calc
diff --git a/addons/gm2calc/MSSMNoFV_onshell_problems.hpp b/addons/gm2calc/MSSMNoFV_onshell_problems.hpp
new file mode 100644
index 000000000..2bbb815c3
--- /dev/null
+++ b/addons/gm2calc/MSSMNoFV_onshell_problems.hpp
@@ -0,0 +1,44 @@
+// ====================================================================
+// This file is part of GM2Calc.
+//
+// GM2Calc is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published
+// by the Free Software Foundation, either version 3 of the License,
+// or (at your option) any later version.
+//
+// GM2Calc is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with GM2Calc.  If not, see
+// <http://www.gnu.org/licenses/>.
+// ====================================================================
+
+#ifndef MSSMNoFV_onshell_PROBLEMS_H
+#define MSSMNoFV_onshell_PROBLEMS_H
+
+#include <iosfwd>
+#include <string>
+
+namespace gm2calc {
+
+class MSSMNoFV_onshell_problems {
+public:
+   MSSMNoFV_onshell_problems();
+   void clear();
+   void flag_tachyon(const std::string&);
+   bool have_problem() const;
+   void print(std::ostream&) const;
+
+private:
+   bool have_tachyon;
+   std::string tachyonic_particle;
+};
+
+std::ostream& operator<<(std::ostream&, const MSSMNoFV_onshell_problems&);
+
+} // namespace gm2calc
+
+#endif
diff --git a/addons/gm2calc/MSSMNoFV_onshell_soft_parameters.cpp b/addons/gm2calc/MSSMNoFV_onshell_soft_parameters.cpp
index 17e6b3cc1..452b68ce0 100644
--- a/addons/gm2calc/MSSMNoFV_onshell_soft_parameters.cpp
+++ b/addons/gm2calc/MSSMNoFV_onshell_soft_parameters.cpp
@@ -20,7 +20,6 @@
 
 #include <iostream>
 
-namespace flexiblesusy {
 namespace gm2calc {
 
 MSSMNoFV_onshell_soft_parameters::MSSMNoFV_onshell_soft_parameters()
@@ -100,4 +99,3 @@ std::ostream& operator<<(std::ostream& ostr, const MSSMNoFV_onshell_soft_paramet
 }
 
 } // namespace flexiblesusy
-} // namespace gm2calc
diff --git a/addons/gm2calc/MSSMNoFV_onshell_soft_parameters.hpp b/addons/gm2calc/MSSMNoFV_onshell_soft_parameters.hpp
index 6a0ce4293..4e0b66ef9 100644
--- a/addons/gm2calc/MSSMNoFV_onshell_soft_parameters.hpp
+++ b/addons/gm2calc/MSSMNoFV_onshell_soft_parameters.hpp
@@ -23,7 +23,6 @@
 
 #include <iosfwd>
 
-namespace flexiblesusy {
 namespace gm2calc {
 
 class MSSMNoFV_onshell_soft_parameters : public MSSMNoFV_onshell_susy_parameters {
@@ -111,6 +110,5 @@ class MSSMNoFV_onshell_soft_parameters : public MSSMNoFV_onshell_susy_parameters
 std::ostream& operator<<(std::ostream&, const MSSMNoFV_onshell_soft_parameters&);
 
 } // namespace gm2calc
-} // namespace flexiblesusy
 
 #endif
diff --git a/addons/gm2calc/MSSMNoFV_onshell_susy_parameters.cpp b/addons/gm2calc/MSSMNoFV_onshell_susy_parameters.cpp
index fea3b0f85..370fe8a64 100644
--- a/addons/gm2calc/MSSMNoFV_onshell_susy_parameters.cpp
+++ b/addons/gm2calc/MSSMNoFV_onshell_susy_parameters.cpp
@@ -20,7 +20,6 @@
 
 #include <iostream>
 
-namespace flexiblesusy {
 namespace gm2calc {
 
 MSSMNoFV_onshell_susy_parameters::MSSMNoFV_onshell_susy_parameters()
@@ -79,4 +78,3 @@ std::ostream& operator<<(std::ostream& ostr, const MSSMNoFV_onshell_susy_paramet
 }
 
 } // namespace gm2calc
-} // namespace flexiblesusy
diff --git a/addons/gm2calc/MSSMNoFV_onshell_susy_parameters.hpp b/addons/gm2calc/MSSMNoFV_onshell_susy_parameters.hpp
index 288e796f7..5009d80b2 100644
--- a/addons/gm2calc/MSSMNoFV_onshell_susy_parameters.hpp
+++ b/addons/gm2calc/MSSMNoFV_onshell_susy_parameters.hpp
@@ -22,7 +22,6 @@
 #include <iosfwd>
 #include <Eigen/Core>
 
-namespace flexiblesusy {
 namespace gm2calc {
 
 class MSSMNoFV_onshell_susy_parameters {
@@ -85,6 +84,5 @@ class MSSMNoFV_onshell_susy_parameters {
 std::ostream& operator<<(std::ostream&, const MSSMNoFV_onshell_susy_parameters&);
 
 } // namespace gm2calc
-} // namespace flexiblesusy
 
 #endif
diff --git a/addons/gm2calc/ffunctions.cpp b/addons/gm2calc/ffunctions.cpp
index d3df7c8a2..c829c5bd3 100644
--- a/addons/gm2calc/ffunctions.cpp
+++ b/addons/gm2calc/ffunctions.cpp
@@ -26,9 +26,11 @@
 
 #define ERROR(message) std::cerr << "Error: " << message << '\n';
 
-namespace flexiblesusy {
+namespace gm2calc {
+
+using namespace flexiblesusy;
 
-namespace {
+const double precision = 1e-10;
 
 double dilog(double x) {
    // The DiLogarithm function
@@ -110,10 +112,6 @@ std::complex<double> dilog(const std::complex<double>& x) {
    return std::complex<double>(ansreal, ansimag);
 }
 
-} // anonymous namespace
-
-namespace gm2calc {
-
 double sqr(double x) { return x*x; }
 double cube(double x) { return x*x*x; }
 double quad(double x) { return sqr(x)*sqr(x); }
@@ -124,23 +122,38 @@ double signed_abs_sqrt(double x) {
 }
 
 double F1C(double x) {
-   if(is_equal(x, 1.))
-      ERROR("F1C: x must not be 1 !");
+   if (is_equal(x, 1., 0.01)) {
+      return 1. - 0.6*(-1 + x) + 0.4*sqr(-1 + x)
+         - 0.2857142857142857*cube(-1 + x)
+         + 0.21428571428571427*quad(-1 + x)
+         - 0.16666666666666666*std::pow(-1 + x,5)
+         + 0.13333333333333333*std::pow(-1 + x,6);
+   }
 
    return 2. / quad(1. - x) * (2. + 3. * x - 6. * sqr(x)
                            + cube(x) + 6. * x * log(x));
 }
 
 double F2C(double x) {
-   if(is_equal(x, 1.))
-      ERROR("F2C: x must not be 1 !");
+   if (is_equal(x, 1., 0.01)) {
+      return 1. - 0.75*(-1 + x) + 0.6*sqr(-1 + x) - 0.5*cube(-1 + x)
+         + 0.42857142857142855*quad(-1 + x)
+         - 0.375*std::pow(-1 + x,5)
+         + 0.3333333333333333*std::pow(-1 + x,6);
+   }
 
    return 3. / (2. * cube(1. - x)) * (- 3. + 4. * x - sqr(x) - 2. * log(x));
 }
 
 double F3C(double x) {
-   if(is_equal(x, 1.))
-      ERROR("F3C: x must not be 1 !");
+   if (is_equal(x, 1., 0.07)) {
+      return 1. + 0.9012765957446807*(-1 + x)
+         - 1.2235460992907807*sqr(-1 + x)
+         + 1.2279808944854547*cube(-1 + x)
+         - 1.1530221450282228*quad(-1 + x)
+         + 1.0620712114633095*std::pow(-1 + x,5)
+         - 0.9740314565542522*std::pow(-1 + x,6);
+   }
 
    return ( 4. / (141. * quad(1. - x)) * ((1. - x) * (151. * sqr(x) - 335. * x + 592.)
              + 6. * (21. * cube(x) - 108. * sqr(x) - 93. * x + 50.) * log(x)
@@ -149,8 +162,14 @@ double F3C(double x) {
 }
 
 double F4C(double x) {
-   if(is_equal(x, 1.))
-      ERROR("F4C: x must not be 1 !");
+   if (is_equal(x, 1., 0.07)) {
+      return 1. - 0.3688524590163934*(-1 + x)
+         + 0.15426229508196715*sqr(-1 + x)
+         - 0.05573770491803279*cube(-1 + x)
+         + 0.0037738374038140845*quad(-1 + x)
+         + 0.025907494145199085*std::pow(-1 + x,5)
+         - 0.04369818965837698*std::pow(-1 + x,6);
+   }
 
    return ( - 9. / (122. * cube(1. - x)) * (8. * (sqr(x) - 3. * x + 2.)
              +  (11. * sqr(x) - 40. * x + 5.) * log(x)
@@ -159,23 +178,38 @@ double F4C(double x) {
 }
 
 double F1N(double x) {
-   if(is_equal(x, 1.))
-      ERROR("F1N: x must not be 1 !");
+   if (is_equal(x, 1., 0.01)) {
+      return 1. - 0.4*(-1 + x) + 0.2*sqr(-1 + x)
+         - 0.1142857142857143*cube(-1 + x)
+         + 0.07142857142857142*quad(-1 + x)
+         - 0.047619047619047616*std::pow(-1 + x,5)
+         + 0.03333333333333333*std::pow(-1 + x,6);
+   }
 
    return 2. / quad(1. - x) * (1. - 6. * x + 3. * sqr(x)
                           + 2. * cube(x) - 6. * sqr(x) * log(x));
 }
 
 double F2N(double x) {
-   if(is_equal(x, 1.))
-      ERROR("F2N: x must not be 1 !");
+   if(is_equal(x, 1., 0.01)){
+      return 1. - 0.5*(-1 + x) + 0.3*sqr(-1 + x)
+         - 0.2*cube(-1 + x) + 0.14285714285714285*quad(-1 + x)
+         - 0.10714285714285714*std::pow(-1 + x,5)
+         + 0.08333333333333333*std::pow(-1 + x,6);
+   }
 
    return 3. / cube(1. - x) * (1. - sqr(x) + 2. * x * log(x));
 }
 
 double F3N(double x) {
-   if(is_equal(x, 1.))
-      ERROR("F3N: x must not be 1 !");
+   if (is_equal(x, 1., 0.07)) {
+      return 1. + 0.08685714285714365*(-1 + x)
+         - 0.1641904761904751*sqr(-1 + x)
+         + 0.13662973760932912*cube(-1 + x)
+         - 0.1038192419825078*quad(-1 + x)
+         + 0.07823129251700683*std::pow(-1 + x,5)
+         - 0.05960544217687109*std::pow(-1 + x,6);
+   }
 
    return 4. / (105. * quad(1. - x)) * ((1. - x) * (- 97. * sqr(x) - 529. * x + 2.)
             + 6. * sqr(x) * (13. * x + 81.) * log(x)
@@ -183,27 +217,75 @@ double F3N(double x) {
 }
 
 double F4N(double x) {
-   if(is_equal(x, 1.))
-      ERROR("F4N: x must not be 1 !");
+   if (is_equal(x, 1., 0.07)) {
+      return 1. + 6.245004513516506e-17*(-1 + x)
+         - 0.13875*sqr(-1 + x)
+         + 0.1475*cube(-1 + x)
+         - 0.13163265306122446*quad(-1 + x)
+         + 0.1128826530612245*std::pow(-1 + x,5)
+         - 0.09610615079365077*std::pow(-1 + x,6);
+   }
 
    return - 2.25 / cube(1. - x) * ((x + 3.) * (x * log(x) + x - 1.)
                                   + (6. * x + 2.) * dilog(1. - x));
 }
 
-double Fa(double x, double y) {
-   if(is_equal(x, y))
-      ERROR("Fa: x must not be equal y!");
+/// Fb(x,1)
+double Fb1(double x) {
+   return (cube(x) - 6*sqr(x) + 3*x + 2. + 3*x*log(sqr(x)))
+      / (6*quad(x - 1.));
+}
 
-   return - (G3(x) - G3(y)) / (x - y);
+/// Fb(x,x)
+double Fbx(double x) {
+   return 0.5 * (sqr(x) + 4*x - 5. - (2*x + 1.)*log(sqr(x)))
+      / quad(x - 1.);
 }
 
 double Fb(double x, double y) {
-   if(is_equal(x, y))
-      ERROR("Fb: x must not be equal y!");
+   if (is_equal(x, 1., precision) && is_equal(y, 1., precision))
+      return 1./12.;
+
+   if (is_equal(x, 1., precision))
+      return Fb1(y);
+
+   if (is_equal(y, 1., precision))
+      return Fb1(x);
+
+   if (is_equal(x, y, precision))
+      return Fbx(x);
 
    return - (G4(x) - G4(y)) / (x - y);
 }
 
+/// Fa(x,1)
+double Fa1(double x) {
+   return 0.25 * (cube(x) - 4*sqr(x) + 11*x - 8. - (x + 2.)*log(sqr(x)))
+      / quad(x - 1.) + 0.5 * Fb1(x);
+}
+
+/// Fa(x,x)
+double Fax(double x) {
+   return 0.25 * (cube(x) - 16*sqr(x) + 11*x + 4. + x*(2*x + 7.)*log(sqr(x)))
+      / (x * quad(x - 1.)) + 0.5 * Fbx(x);
+}
+
+double Fa(double x, double y) {
+   if (is_equal(x, 1., precision) && is_equal(y, 1., precision))
+      return 0.25;
+
+   if (is_equal(x, 1., precision))
+      return Fa1(y);
+
+   if (is_equal(y, 1., precision))
+      return Fa1(x);
+
+   if (is_equal(x, y, precision))
+      return Fax(x);
+
+   return - (G3(x) - G3(y)) / (x - y);
+}
+
 double G3(double x) {
    if(is_equal(x, 1.))
       ERROR("G3: x must not be 1 !");
@@ -225,13 +307,38 @@ double G4(double x) {
  * \f$\frac{1}{z} H_2(\frac{x}{z},\frac{y}{z}) = - I(x,y,z)\f$.
  */
 double H2(double x, double y) {
-   if (is_equal(x,1.) || is_equal(y,1.) || is_equal(x,y))
-      ERROR("H2(" << x << "," << y << ") is not well-defined!");
+   if (is_equal(x, 1., precision) && is_equal(y, 1., precision))
+      return -0.5;
+
+   if (is_equal(x, 1., precision))
+      return (-1. + y - y*log(y))/sqr(-1. + y);
+
+   if (is_equal(y, 1., precision))
+      return (-1. + x - x*log(x))/sqr(-1. + x);
+
+   if (is_equal(x, y, precision))
+      return (1 - y + log(y))/sqr(-1 + y);
 
    return x * log(x) / ((1-x)*(x-y)) + y * log(y) / ((1-y)*(y-x));
 }
 
+/// Iabc(a,a,c)
+double Iaac(double a, double c) {
+   return (sqr(a) - sqr(c) - sqr(c)*log(sqr(a/c))) / sqr(sqr(a) - sqr(c));
+}
+
 double Iabc(double a, double b, double c) {
+   if (is_equal(a,b) && is_equal(b,c))
+      return 0.5 / sqr(a);
+
+   if (is_equal(a,b))
+      return Iaac(a,c);
+
+   if (is_equal(b,c))
+      return Iaac(b,a);
+
+   if (is_equal(a,c))
+      return Iaac(a,b);
 
    return ( (sqr(a * b) * log(sqr(a / b))
            + sqr(b * c) * log(sqr(b / c))
@@ -277,4 +384,3 @@ double f_sferm(double z) {
 }
 
 } // namespace gm2calc
-} // namespace flexiblesusy
diff --git a/addons/gm2calc/ffunctions.hpp b/addons/gm2calc/ffunctions.hpp
index 497448afb..4513d78ee 100644
--- a/addons/gm2calc/ffunctions.hpp
+++ b/addons/gm2calc/ffunctions.hpp
@@ -19,7 +19,6 @@
 #ifndef GM2_FFUNCTIONS_H
 #define GM2_FFUNCTIONS_H
 
-namespace flexiblesusy {
 namespace gm2calc {
 
 /// \f$F_1^C(x)\f$, Eq. (54) arXiv:hep-ph/0609168
@@ -60,6 +59,5 @@ double signed_sqr(double);
 double sqr(double);
 
 } // namespace gm2calc
-} // namespace flexiblesusy
 
 #endif
diff --git a/addons/gm2calc/gm2_1loop.cpp b/addons/gm2calc/gm2_1loop.cpp
index 52797acc9..ece7f0944 100644
--- a/addons/gm2calc/gm2_1loop.cpp
+++ b/addons/gm2calc/gm2_1loop.cpp
@@ -22,10 +22,18 @@
 
 #include <complex>
 
-namespace flexiblesusy {
+/**
+ * \file gm2_1loop.cpp
+ *
+ * Contains functions necessary to calculate the SUSY contributions
+ * for g-2 at the 1-loop level.
+ */
+
 namespace gm2calc {
 
 /**
+ * \fn calculate_amu_1loop_non_tan_beta_resummed
+ *
  * Calculates full 1-loop SUSY contribution to (g-2), Eq. (45) of
  * arXiv:hep-ph/0609168.
  *
@@ -504,4 +512,3 @@ double delta_bottom_correction(const MSSMNoFV_onshell& model)
 }
 
 } // namespace gm2calc
-} // namespace flexiblesusy
diff --git a/addons/gm2calc/gm2_1loop.hpp b/addons/gm2calc/gm2_1loop.hpp
index ca48227fc..b45b279f9 100644
--- a/addons/gm2calc/gm2_1loop.hpp
+++ b/addons/gm2calc/gm2_1loop.hpp
@@ -22,7 +22,6 @@
 #include <cmath>
 #include <Eigen/Core>
 
-namespace flexiblesusy {
 namespace gm2calc {
 
 class MSSMNoFV_onshell;
@@ -77,6 +76,5 @@ Eigen::Matrix<double,4,2> x_im(const MSSMNoFV_onshell&);
 Eigen::Array<double,2,1> x_k(const MSSMNoFV_onshell&);
 
 } // namespace gm2calc
-} // namespace flexiblesusy
 
 #endif
diff --git a/addons/gm2calc/gm2_2loop.cpp b/addons/gm2calc/gm2_2loop.cpp
index f71fed148..f31ad4f33 100644
--- a/addons/gm2calc/gm2_2loop.cpp
+++ b/addons/gm2calc/gm2_2loop.cpp
@@ -24,9 +24,15 @@
 #include <complex>
 #include <cmath>
 
-namespace flexiblesusy {
 namespace gm2calc {
 
+/**
+ * \file gm2_2loop.cpp
+ *
+ * Contains functions necessary to calculate the SUSY contributions
+ * for g-2 at the 2-loop level.
+ */
+
 // fermion/sfermion corrections, log-approximations
 
 /**
@@ -282,8 +288,7 @@ double amuChi0Photonic(const MSSMNoFV_onshell& model) {
 
 /**
  * The following functions include resummation of 1/(1+Deltamu) within
- * muon Yukawa coupling but not in other Yukawa couplings (in
- * particular not in tau, bottom)
+ * the muon, tau and bottom Yukawa couplings.
  */
 
 /**
@@ -510,4 +515,3 @@ double amua2LCha(const MSSMNoFV_onshell& model) {
 }
 
 } // namespace gm2calc
-} // namespace flexiblesusy
diff --git a/addons/gm2calc/gm2_2loop.hpp b/addons/gm2calc/gm2_2loop.hpp
index 2877d31f6..d2571ad85 100644
--- a/addons/gm2calc/gm2_2loop.hpp
+++ b/addons/gm2calc/gm2_2loop.hpp
@@ -21,7 +21,6 @@
 
 #include <Eigen/Core>
 
-namespace flexiblesusy {
 namespace gm2calc {
 
 class MSSMNoFV_onshell;
@@ -64,6 +63,5 @@ Eigen::Matrix<std::complex<double>,2,2> lambda_sbot(const MSSMNoFV_onshell&);
 Eigen::Matrix<std::complex<double>,2,2> lambda_stau(const MSSMNoFV_onshell&);
 
 } // namespace gm2calc
-} // namespace flexiblesusy
 
 #endif
diff --git a/addons/gm2calc/gm2_error.hpp b/addons/gm2calc/gm2_error.hpp
index a52fd8faf..ee20b0795 100644
--- a/addons/gm2calc/gm2_error.hpp
+++ b/addons/gm2calc/gm2_error.hpp
@@ -19,11 +19,27 @@
 #ifndef GM2_ERROR_H
 #define GM2_ERROR_H
 
-#include "error.hpp"
-
-namespace flexiblesusy {
 namespace gm2calc {
 
+class Error {
+public:
+   virtual ~Error() {}
+   virtual std::string what() const = 0;
+};
+
+/**
+ * @class SetupError
+ * @brief Spectrum generator was not setup correctly
+ */
+class SetupError : public Error {
+public:
+   explicit SetupError(const std::string& message_) : message(message_) {}
+   virtual ~SetupError() {}
+   virtual std::string what() const { return message; }
+private:
+   std::string message;
+};
+
 class EInvalidInput : public Error {
 public:
    explicit EInvalidInput(const std::string& message_) : message(message_) {}
@@ -33,7 +49,24 @@ class EInvalidInput : public Error {
    std::string message;
 };
 
+class EPhysicalProblem : public Error {
+public:
+   explicit EPhysicalProblem(const std::string& message_) : message(message_) {}
+   virtual ~EPhysicalProblem() {}
+   virtual std::string what() const { return message; }
+private:
+   std::string message;
+};
+
+class ReadError : public Error {
+public:
+   ReadError(const std::string& message_) : message(message_) {}
+   virtual ~ReadError() {}
+   virtual std::string what() const { return message; }
+private:
+   std::string message;
+};
+
 } // namespace gm2calc
-} // namespace flexiblesusy
 
 #endif
diff --git a/addons/gm2calc/gm2_slha_io.cpp b/addons/gm2calc/gm2_slha_io.cpp
index 7142f80a7..00f733ee7 100644
--- a/addons/gm2calc/gm2_slha_io.cpp
+++ b/addons/gm2calc/gm2_slha_io.cpp
@@ -26,9 +26,10 @@
 
 #include <Eigen/Core>
 
-namespace flexiblesusy {
 namespace gm2calc {
 
+using namespace flexiblesusy;
+
 #define ERROR(message) std::cerr << "Error: " << message << '\n';
 #define WARNING(message) std::cerr << "Warning: " << message << '\n';
 
@@ -231,6 +232,35 @@ void GM2_slha_io::fill_block_entry(const std::string& block_name,
    }
 }
 
+/**
+ * Fills a block entry with a string.  If the block or the entry do
+ * not exist, the block / entry is created.
+ *
+ * @param block_name block name
+ * @param entry number of the entry
+ * @param description comment
+ */
+void GM2_slha_io::fill_block_entry(const std::string& block_name,
+                                   unsigned entry,
+                                   const std::string& description)
+{
+   std::ostringstream sstr;
+   sstr << FORMAT_SPINFO(entry, description);
+
+   SLHAea::Coll::const_iterator block =
+      data.find(data.cbegin(), data.cend(), block_name);
+
+   if (block == data.cend()) {
+      // create new block
+      std::ostringstream block;
+      block << "Block " << block_name << '\n'
+            << sstr.str();
+      set_block(block, GM2_slha_io::back);
+   } else {
+      data[block_name][SLHAea::to<double>(entry)] = sstr.str();
+   }
+}
+
 double read_scale(const GM2_slha_io& slha_io)
 {
    char const * const drbar_blocks[] =
@@ -578,4 +608,3 @@ void GM2_slha_io::process_gm2calcconfig_tuple(Config_options& config_options,
 }
 
 } // namespace gm2calc
-} // namespace flexiblesusy
diff --git a/addons/gm2calc/gm2_slha_io.hpp b/addons/gm2calc/gm2_slha_io.hpp
index 22d42aab4..3c45dbebf 100644
--- a/addons/gm2calc/gm2_slha_io.hpp
+++ b/addons/gm2calc/gm2_slha_io.hpp
@@ -25,10 +25,9 @@
 #include <Eigen/Core>
 #include <boost/format.hpp>
 #include "slhaea.h"
-#include "error.hpp"
+#include "gm2_error.hpp"
 #include "numerics2.hpp"
 
-namespace flexiblesusy {
 namespace gm2calc {
 
 class MSSMNoFV_onshell;
@@ -68,18 +67,10 @@ class GM2_slha_io {
    typedef std::function<void(int, double)> Tuple_processor;
    enum Position { front, back };
 
-   class ReadError : public Error {
-   public:
-      ReadError(const std::string& message_) : message(message_) {}
-      virtual ~ReadError() {}
-      virtual std::string what() const { return message; }
-   private:
-      std::string message;
-   };
-
    void clear();
 
    // reading functions
+   const SLHAea::Coll& get_data() const { return data; }
    bool block_exists(const std::string&) const;
    void read_from_file(const std::string&);
    void read_from_source(const std::string&);
@@ -95,8 +86,8 @@ class GM2_slha_io {
    void set_block(const std::ostringstream&, Position position = back);
    void write_to_file(const std::string&);
    void write_to_stream(std::ostream& = std::cout);
-
    void fill_block_entry(const std::string&, unsigned, double, const std::string&);
+   void fill_block_entry(const std::string&, unsigned, const std::string&);
 
    static void process_gm2calcconfig_tuple(Config_options&, int, double);
    static void process_gm2calcinput_tuple(MSSMNoFV_onshell&, int, double);
@@ -187,6 +178,5 @@ void fill_slha(const GM2_slha_io&, MSSMNoFV_onshell&);
 void fill(const GM2_slha_io&, Config_options&);
 
 } // namespace gm2calc
-} // namespace flexiblesusy
 
 #endif
diff --git a/addons/gm2calc/gm2calc.cpp b/addons/gm2calc/gm2calc.cpp
index 2ec31eec9..7466f4d56 100644
--- a/addons/gm2calc/gm2calc.cpp
+++ b/addons/gm2calc/gm2calc.cpp
@@ -18,7 +18,7 @@
 
 #include "gm2_1loop.hpp"
 #include "gm2_2loop.hpp"
-#include "error.hpp"
+#include "gm2_error.hpp"
 #include "config.h"
 
 #include "gm2_slha_io.hpp"
@@ -126,28 +126,11 @@ void setup_model(gm2calc::MSSMNoFV_onshell& model,
       model.calculate_masses();
       break;
    default:
-      throw SetupError("Unknown input option");
+      throw gm2calc::SetupError("Unknown input option");
       break;
    }
 }
 
-/**
- * Calculate most precise value of a_mu for a given set op parameters
- */
-double calculate_amu_best(gm2calc::MSSMNoFV_onshell& model)
-{
-   const double gm2_1l_tan_beta_resummed =
-      gm2calc::calculate_amu_1loop(model);
-   const double gm2_2l_tan_beta_resummed =
-      + gm2calc::amu2LFSfapprox(model)
-      + gm2calc::amuChipmPhotonic(model)
-      + gm2calc::amuChi0Photonic(model);
-
-   const double gm2_best = gm2_1l_tan_beta_resummed + gm2_2l_tan_beta_resummed;
-
-   return gm2_best;
-}
-
 /**
  * Prints detailed a_mu calculation (1-loop w/ and w/o tan(beta)
  * resumation, 2-loop, and differenc contributions).
@@ -248,10 +231,10 @@ void print_amu_detailed(const gm2calc::MSSMNoFV_onshell& model)
  *
  * @return a_mu
  */
-double calculate_amu(gm2calc::MSSMNoFV_onshell& model,
+double calculate_amu(const gm2calc::MSSMNoFV_onshell& model,
                      const gm2calc::Config_options& config_options)
 {
-   double result = 0.;
+   double result = std::numeric_limits<double>::signaling_NaN();
 
    switch (config_options.loop_order) {
    case 0:
@@ -280,35 +263,23 @@ double calculate_amu(gm2calc::MSSMNoFV_onshell& model,
    return result;
 }
 
-int main(int argc, const char* argv[])
+/**
+ * Calculates amu and prints it to std::cout.  The output format
+ * depends on the \a config_options .
+ *
+ * @param model model
+ * @param slha_io SLHA object
+ * @param config_options configuration options
+ */
+void print_amu(const gm2calc::MSSMNoFV_onshell& model,
+               gm2calc::GM2_slha_io& slha_io,
+               const gm2calc::Config_options& config_options)
 {
-   Gm2_cmd_line_options options(get_cmd_line_options(argc, argv));
-
-   if (options.input_source.empty()) {
-      ERROR("No input source given!\n"
-            "   Please provide an SLHA input via the option --slha-input-file=\n"
-            "   or a GM2Calc input via the option --gm2calc-input-file=");
-      return EXIT_FAILURE;
-   }
-
-   gm2calc::MSSMNoFV_onshell model;
-   gm2calc::GM2_slha_io slha_io;
-   gm2calc::Config_options config_options;
-
-   try {
-      slha_io.read_from_source(options.input_source);
-      fill(slha_io, config_options);
-      setup_model(model, slha_io, options);
-   } catch (const Error& error) {
-      ERROR(error.what());
-      return EXIT_FAILURE;
-   }
-
    switch (config_options.output_format) {
    case gm2calc::Config_options::Minimal:
       std::cout << std::setprecision(std::numeric_limits<double>::digits10)
                 << std::scientific
-                << calculate_amu_best(model) << '\n';
+                << calculate_amu(model, config_options) << '\n';
       break;
    case gm2calc::Config_options::Detailed:
       print_amu_detailed(model);
@@ -327,9 +298,72 @@ int main(int argc, const char* argv[])
       break;
    default:
       ERROR("Unknown output format: " << config_options.output_format);
-      return EXIT_FAILURE;
       break;
    }
+}
+
+/**
+ * Prints output if an error has occured.
+ *
+ * @param error error object
+ * @param model model
+ * @param slha_io SLHA object
+ * @param config_options configuration options
+ */
+void print_error(const gm2calc::Error& error,
+                 const gm2calc::MSSMNoFV_onshell& model,
+                 gm2calc::GM2_slha_io& slha_io,
+                 const gm2calc::Config_options& config_options)
+{
+   ERROR(error.what());
+
+   switch (config_options.output_format) {
+   case gm2calc::Config_options::Detailed:
+      // in detailed mode print amu even if an error has occured
+      print_amu_detailed(model);
+      std::cout << '\n'
+                << "================================\n"
+                << "   " << error.what() << '\n'
+                << "================================\n";
+      break;
+   case gm2calc::Config_options::NMSSMTools:
+   case gm2calc::Config_options::SPheno:
+      // print SPINFO block with error description
+      slha_io.fill_block_entry("SPINFO", 1, "GM2Calc");
+      slha_io.fill_block_entry("SPINFO", 2, GM2CALC_VERSION);
+      slha_io.fill_block_entry("SPINFO", 4, error.what());
+      slha_io.write_to_stream(std::cout);
+      break;
+   default:
+      break;
+   }
+}
+
+int main(int argc, const char* argv[])
+{
+   Gm2_cmd_line_options options(get_cmd_line_options(argc, argv));
+
+   if (options.input_source.empty()) {
+      ERROR("No input source given!\n"
+            "   Please provide an SLHA input via the option --slha-input-file=\n"
+            "   or a GM2Calc input via the option --gm2calc-input-file=");
+      return EXIT_FAILURE;
+   }
+
+   gm2calc::MSSMNoFV_onshell model;
+   gm2calc::GM2_slha_io slha_io;
+   gm2calc::Config_options config_options;
+   int exit_code = 0;
+
+   try {
+      slha_io.read_from_source(options.input_source);
+      fill(slha_io, config_options);
+      setup_model(model, slha_io, options);
+      print_amu(model, slha_io, config_options);
+   } catch (const gm2calc::Error& error) {
+      print_error(error, model, slha_io, config_options);
+      exit_code = EXIT_FAILURE;
+   }
 
-   return 0;
+   return exit_code;
 }
diff --git a/addons/gm2calc/gm2scan.cpp b/addons/gm2calc/gm2scan.cpp
new file mode 100644
index 000000000..a5718b741
--- /dev/null
+++ b/addons/gm2calc/gm2scan.cpp
@@ -0,0 +1,140 @@
+// ====================================================================
+// This file is part of GM2Calc.
+//
+// GM2Calc is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published
+// by the Free Software Foundation, either version 3 of the License,
+// or (at your option) any later version.
+//
+// GM2Calc is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with GM2Calc.  If not, see
+// <http://www.gnu.org/licenses/>.
+// ====================================================================
+
+#include "gm2_1loop.hpp"
+#include "gm2_2loop.hpp"
+#include "gm2_error.hpp"
+#include "MSSMNoFV_onshell.hpp"
+
+#include <cstdio>
+#include <iostream>
+
+double calculate_amu(gm2calc::MSSMNoFV_onshell& model)
+{
+   return gm2calc::calculate_amu_1loop(model)
+      + gm2calc::amu2LFSfapprox(model)
+      + gm2calc::amuChipmPhotonic(model)
+      + gm2calc::amuChi0Photonic(model);
+}
+
+gm2calc::MSSMNoFV_onshell setup()
+{
+   gm2calc::MSSMNoFV_onshell model;
+
+   const double g3 = 1.06459;
+   const double MW = 80.385;
+   const double MZ = 91.1876;
+   const double ME = 0.00051;
+   const double MM = 0.105658;
+   const double ML = 1.777;
+   const double MT = 173.5;
+   const double MB = 3.;
+   const double cW = MW / MZ;
+   const double EL = model.get_EL(); // default
+   const double g2 = EL / std::sqrt(1. - cW*cW);
+   const double vev = 2. * MW / g2;
+   Eigen::Matrix<double,3,3> Ae;
+   Eigen::Matrix<double,3,3> Au;
+   Eigen::Matrix<double,3,3> Ad;
+   Eigen::Matrix<double,3,3> mq2;
+   Eigen::Matrix<double,3,3> ml2;
+   Eigen::Matrix<double,3,3> md2;
+   Eigen::Matrix<double,3,3> mu2;
+   Eigen::Matrix<double,3,3> me2;
+   double Mu, M1, M2, M3;
+   double scale;
+   double MA0;
+   double TB;
+
+   Mu = 350;
+   M1 = 150;
+   M2 = 300;
+   M3 = 1;
+   TB = 10.;
+   scale = 454.7;
+   MA0 = 1500;
+
+   Au.setZero();
+   Ad.setZero();
+   Ae.setZero();
+
+   mq2 << 400*400,           0,           0,
+                    0, 400*400,           0,
+                    0,       0,      400*400;
+
+   ml2 = md2 = mu2 = me2 = mq2;
+
+   // set parameters
+   model.set_g1(std::sqrt(5. / 3.) * EL / cW);
+   model.set_g2(g2);
+   model.set_g3(g3);
+   model.set_vu(vev / sqrt(1. + 1. / (TB*TB)));
+   model.set_vd(model.get_vu() / TB);
+   model.set_Ae(Ae);
+   model.set_Ad(Ad);
+   model.set_Au(Au);
+   model.set_Mu(Mu);
+   model.set_mq2(mq2);
+   model.set_ml2(ml2);
+   model.set_md2(md2);
+   model.set_mu2(mu2);
+   model.set_me2(me2);
+   model.set_MassB(M1);
+   model.set_MassWB(M2);
+   model.set_MassG(M3);
+   model.set_scale(scale);
+   model.set_MA0(MA0);
+
+   model.get_physical().MVWm = MW;
+   model.get_physical().MVZ = MZ;
+   model.get_physical().MFe = ME;
+   model.get_physical().MFm = MM;
+   model.get_physical().MFtau = ML;
+   model.get_physical().MFt = MT;
+   model.get_physical().MFb = MB;
+
+   return model;
+}
+
+int main()
+{
+   const double tanb_start = 2.;
+   const double tanb_stop = 100.;
+
+   printf("# %14s %16s\n", "tan(beta)", "amu");
+
+   for (unsigned n = 0; n < 100; n++) {
+      double amu;
+      const double tanb = tanb_start + (tanb_stop - tanb_start) * n / 100;
+
+      gm2calc::MSSMNoFV_onshell model(setup());
+      model.set_TB(tanb);
+
+      try {
+         model.calculate_masses();
+         amu = calculate_amu(model);
+      } catch (const gm2calc::Error& e) {
+         std::cerr << "Error: " << e.what() << std::endl;
+         amu = std::numeric_limits<double>::signaling_NaN();
+      }
+
+      printf("%16.8e %16.8e\n", tanb, amu);
+   }
+
+   return 0;
+}
diff --git a/addons/gm2calc/module.mk b/addons/gm2calc/module.mk
index db08511b7..f58ad78f4 100644
--- a/addons/gm2calc/module.mk
+++ b/addons/gm2calc/module.mk
@@ -18,12 +18,14 @@ LIBgm2calc_SRC := \
 		$(DIR)/MSSMNoFV_onshell.cpp \
 		$(DIR)/MSSMNoFV_onshell_mass_eigenstates.cpp \
 		$(DIR)/MSSMNoFV_onshell_physical.cpp \
+		$(DIR)/MSSMNoFV_onshell_problems.cpp \
 		$(DIR)/MSSMNoFV_onshell_soft_parameters.cpp \
 		$(DIR)/MSSMNoFV_onshell_susy_parameters.cpp
 
 # main()
 EXEgm2calc_SRC := \
-		$(DIR)/gm2calc.cpp
+		$(DIR)/gm2calc.cpp \
+		$(DIR)/gm2scan.cpp
 
 # header files
 LIBgm2calc_HDR := \
@@ -35,6 +37,7 @@ LIBgm2calc_HDR := \
 		$(DIR)/MSSMNoFV_onshell.hpp \
 		$(DIR)/MSSMNoFV_onshell_mass_eigenstates.hpp \
 		$(DIR)/MSSMNoFV_onshell_physical.hpp \
+		$(DIR)/MSSMNoFV_onshell_problems.hpp \
 		$(DIR)/MSSMNoFV_onshell_soft_parameters.hpp \
 		$(DIR)/MSSMNoFV_onshell_susy_parameters.hpp
 
@@ -118,7 +121,7 @@ $(LIBgm2calc_DEP) $(EXEgm2calc_DEP) $(LIBgm2calc_OBJ) $(EXEgm2calc_OBJ): CPPFLAG
 $(LIBgm2calc): $(LIBgm2calc_OBJ)
 		$(MAKELIB) $@ $^
 
-$(EXEgm2calc_EXE): $(EXEgm2calc_OBJ) $(LIBgm2calc) $(LIBMSSMNoFV_onshell) $(LIBFLEXI) $(LIBLEGACY) $(filter-out -%,$(LOOPFUNCLIBS))
+$(DIR)/%.x: $(DIR)/%.o $(LIBgm2calc) $(LIBMSSMNoFV_onshell) $(LIBFLEXI) $(LIBLEGACY) $(filter-out -%,$(LOOPFUNCLIBS))
 		$(CXX) -o $@ $(call abspathx,$^) $(filter -%,$(LOOPFUNCLIBS)) $(GSLLIBS) $(BOOSTTHREADLIBS) $(THREADLIBS) $(LAPACKLIBS) $(BLASLIBS) $(FLIBS)
 
 ALLDEP += $(LIBgm2calc_DEP) $(EXEgm2calc_DEP)
diff --git a/configure b/configure
index 19ba08050..98919d6b1 100755
--- a/configure
+++ b/configure
@@ -10,7 +10,7 @@ FLEXIBLESUSY_PATCH=1
 FLEXIBLESUSY_EXTRA=""
 FLEXIBLESUSY_VERSION="${FLEXIBLESUSY_MAJOR}.${FLEXIBLESUSY_MINOR}.${FLEXIBLESUSY_PATCH}${FLEXIBLESUSY_EXTRA}"
 GIT_COMMIT="`git describe --tags 2> /dev/null || echo unknown`"
-GM2CALC_VERSION="0.2.3"
+GM2CALC_VERSION="0.2.7"
 
 # directory of this script
 BASEDIR=$(dirname $0)