diff --git a/Core/Core.pro b/Core/Core.pro
index 05f7ce0e86c21d63deb359b555a65dd50ae8110f..38b301de89fc349c777e77ce8a3ae4b11a48f09b 100644
--- a/Core/Core.pro
+++ b/Core/Core.pro
@@ -242,6 +242,7 @@ HEADERS += \
     Tools/inc/LLData.h \
     Tools/inc/Macros.h \
     Tools/inc/MathFunctions.h \
+    Tools/inc/MemberFunctionIntegrator.h \
     Tools/inc/NamedVector.h \
     Tools/inc/Numeric.h \
     Tools/inc/OutputData.h \
diff --git a/Core/FormFactors/inc/FormFactorBigCylinder.h b/Core/FormFactors/inc/FormFactorBigCylinder.h
index f1ac05a3aad99ce1a8514e87694184aac8a487b0..0713f1d5301405e10aa8b35356a2a38b7ccb385a 100644
--- a/Core/FormFactors/inc/FormFactorBigCylinder.h
+++ b/Core/FormFactors/inc/FormFactorBigCylinder.h
@@ -46,7 +46,7 @@ private:
     virtual void init_parameters();
 
     //! approximation to radial function for integration
-    double iTilde(double qR) const;
+    double iTilde(double qR, void *params) const;
 
     //! print class
     void print(std::ostream &ostr) const;
diff --git a/Core/FormFactors/src/FormFactorBigCylinder.cpp b/Core/FormFactors/src/FormFactorBigCylinder.cpp
index 9bfecb8e807cc64f88bf379af35d28f52431a0af..9037ab922e0dfeff373aeab321b806083f08def9 100644
--- a/Core/FormFactors/src/FormFactorBigCylinder.cpp
+++ b/Core/FormFactors/src/FormFactorBigCylinder.cpp
@@ -2,15 +2,9 @@
 #include "StochasticDiracDelta.h"
 
 #include "MathFunctions.h"
+#include "MemberFunctionIntegrator.h"
 #include "Numeric.h"
 
-double wrapFunctionForGSL2(double qR, void* p_unary_function)
-{
-    std::binder1st<std::const_mem_fun1_t<double, FormFactorBigCylinder, double> > *p_f =
-            (std::binder1st<std::const_mem_fun1_t<double, FormFactorBigCylinder, double> > *)p_unary_function;
-    return (*p_f)(qR);
-}
-
 FormFactorBigCylinder::FormFactorBigCylinder(double height, double radius)
 : m_bin_size(1.0)
 {
@@ -63,12 +57,9 @@ complex_t FormFactorBigCylinder::evaluate_for_q(const cvector_t &q) const
     double qRmin = qrR - effective_bin_size/2.0;
     double qRmax = qrR + effective_bin_size/2.0;
 
-    std::binder1st<std::const_mem_fun1_t<double, FormFactorBigCylinder, double> >
-        f_base = std::bind1st(std::mem_fun(&FormFactorBigCylinder::iTilde), this);
-    gsl_function f;
-    f.function = &wrapFunctionForGSL2;
-    f.params = &f_base;
-    double average_intensity = MathFunctions::Integrate1D(&f, qRmin, qRmax)/effective_bin_size;
+    MemberFunctionIntegrator<FormFactorBigCylinder>::mem_function p_mf = &FormFactorBigCylinder::iTilde;
+    MemberFunctionIntegrator<FormFactorBigCylinder> integrator(p_mf, this);
+    double average_intensity = integrator.integrate(qRmin, qRmax, (void*)0);
 
     double J1_qrR_div_qrR = std::sqrt(average_intensity);
     double radial_part = 2.0*M_PI*R*R*J1_qrR_div_qrR;
@@ -76,8 +67,9 @@ complex_t FormFactorBigCylinder::evaluate_for_q(const cvector_t &q) const
     return radial_part*z_part;
 }
 
-double FormFactorBigCylinder::iTilde(double qR) const
+double FormFactorBigCylinder::iTilde(double qR, void *params) const
 {
+    (void)params;
     static double a = 1.0/4.0;
     static double b = std::sqrt(M_PI/3.0/std::sqrt(3.0));
 
diff --git a/Core/Samples/inc/InterferenceFunction2DParaCrystal.h b/Core/Samples/inc/InterferenceFunction2DParaCrystal.h
index 809629a30da81542bf20653faedd8acf71c0c708..4da321a940ba7ed4f3695e152373068dc35bf5d6 100644
--- a/Core/Samples/inc/InterferenceFunction2DParaCrystal.h
+++ b/Core/Samples/inc/InterferenceFunction2DParaCrystal.h
@@ -66,7 +66,7 @@ private:
     virtual void init_parameters();
 
     //! Calculate interference function for fixed rotation xi
-    double interferenceForXi(double xi) const;
+    double interferenceForXi(double xi, void *params) const;
 
     //! calculate interference function for fixed xi in 1d
     double interference1D(double qx, double qy, double xi, size_t index) const;
diff --git a/Core/Samples/src/InterferenceFunction2DParaCrystal.cpp b/Core/Samples/src/InterferenceFunction2DParaCrystal.cpp
index 5b2c7793534658e4cd3a7586b873d20b94c9d9d9..c6d125ef60f023c16f27c220539a8959ea10269f 100644
--- a/Core/Samples/src/InterferenceFunction2DParaCrystal.cpp
+++ b/Core/Samples/src/InterferenceFunction2DParaCrystal.cpp
@@ -1,16 +1,10 @@
 #include "InterferenceFunction2DParaCrystal.h"
 #include "MathFunctions.h"
+#include "MemberFunctionIntegrator.h"
 #include "Exceptions.h"
 
 #include <functional>
 
-double wrapFunctionForGSL(double xi, void* p_unary_function)
-{
-    std::binder1st<std::const_mem_fun1_t<double, InterferenceFunction2DParaCrystal, double> > *p_f =
-            (std::binder1st<std::const_mem_fun1_t<double, InterferenceFunction2DParaCrystal, double> > *)p_unary_function;
-    return (*p_f)(xi);
-}
-
 InterferenceFunction2DParaCrystal::InterferenceFunction2DParaCrystal(double length_1, double length_2, double alpha_lattice, double xi, double corr_length)
 : m_alpha_lattice(alpha_lattice)
 , m_xi(xi)
@@ -54,14 +48,12 @@ double InterferenceFunction2DParaCrystal::evaluate(const cvector_t &q) const
     m_qy = q.y().real();
     double result;
     if (m_integrate_xi) {
-        std::binder1st<std::const_mem_fun1_t<double, InterferenceFunction2DParaCrystal, double> > f_base = std::bind1st(std::mem_fun(&InterferenceFunction2DParaCrystal::interferenceForXi), this);
-        gsl_function f;
-        f.function = &wrapFunctionForGSL;
-        f.params = &f_base;
-        result = MathFunctions::Integrate1D(&f, 0.0, M_PI)/M_PI;
-    }
+        MemberFunctionIntegrator<InterferenceFunction2DParaCrystal>::mem_function p_member_function = &InterferenceFunction2DParaCrystal::interferenceForXi;
+        MemberFunctionIntegrator<InterferenceFunction2DParaCrystal> integrator(p_member_function, this);
+        result = integrator.integrate(0.0, M_PI, (void*)0);
+   }
     else {
-        result = interferenceForXi(m_xi);
+        result = interferenceForXi(m_xi, (void*)0);
     }
     return result;
 }
@@ -119,8 +111,9 @@ void InterferenceFunction2DParaCrystal::init_parameters()
     getParameterPool()->registerParameter("domain_size_2", &m_domain_sizes[1]);
 }
 
-double InterferenceFunction2DParaCrystal::interferenceForXi(double xi) const
+double InterferenceFunction2DParaCrystal::interferenceForXi(double xi, void *params) const
 {
+    (void)params;
     double result = interference1D(m_qx, m_qy, xi, 0)*interference1D(m_qx, m_qy, xi + m_alpha_lattice, 1);
     return result;
 }
diff --git a/Core/Tools/inc/MemberFunctionIntegrator.h b/Core/Tools/inc/MemberFunctionIntegrator.h
new file mode 100644
index 0000000000000000000000000000000000000000..e046c903c72e0b3687269cca153a1112248d309a
--- /dev/null
+++ b/Core/Tools/inc/MemberFunctionIntegrator.h
@@ -0,0 +1,73 @@
+#ifndef MEMBERFUNCTIONINTEGRATOR_H_
+#define MEMBERFUNCTIONINTEGRATOR_H_
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   MemberFunctionIntegrator.h
+//! @brief  Definition of MemberFunctionIntegrator template
+//! @author Scientific Computing Group at FRM II
+//! @date   Nov 26, 2012
+
+#include "gsl/gsl_integration.h"
+
+template <class C> class MemberFunctionIntegrator
+{
+public:
+    //! member function type
+    typedef double (C::*mem_function)(double, void*) const;
+
+    //! structure holding the object and possible extra parameters
+    struct CallBackHolder {
+        const C *m_object_pointer;
+        mem_function m_member_function;
+        void *m_data;
+    };
+
+    //! constructor taking a member function and the object whose member function to integrate
+    MemberFunctionIntegrator(mem_function p_member_function, const C *const p_object);
+
+    //! perform the actual integration over the range [lmin, lmax]
+    double integrate(double lmin, double lmax, void *params);
+
+private:
+    //! static function that can be passed to gsl integrator
+    static double StaticCallBack(double d, void *v) {
+        CallBackHolder *p_cb = static_cast<CallBackHolder *>(v);
+        mem_function mf = p_cb->m_member_function;
+        return (p_cb->m_object_pointer->*mf)(d, p_cb->m_data);
+    }
+    mem_function m_member_function; //!< the member function to integrate
+    const C *mp_object; //!< the object whose member function to integrate
+};
+
+template<class C> MemberFunctionIntegrator<C>::MemberFunctionIntegrator(
+        mem_function p_member_function, const C *const p_object)
+: m_member_function(p_member_function)
+, mp_object(p_object)
+{
+}
+
+template<class C> inline double MemberFunctionIntegrator<C>::integrate(
+        double lmin, double lmax, void* params)
+{
+    CallBackHolder cb = { mp_object, m_member_function, params };
+
+    gsl_function f;
+    f.function = StaticCallBack;
+    f.params = &cb;
+
+    gsl_integration_workspace *ws = gsl_integration_workspace_alloc(200);
+    double result, error;
+    gsl_integration_qag(&f, lmin, lmax, 1e-10, 1e-8, 50, 1, ws, &result, &error);
+    gsl_integration_workspace_free(ws);
+
+    return result;
+}
+
+#endif /* MEMBERFUNCTIONINTEGRATOR_H_ */