Index: sansmodels/src/sans/models/c_extensions/schulz.c
===================================================================
--- sansmodels/src/sans/models/c_extensions/schulz.c	(revision 92ffdb6d3d3c32f163272799857a959c6f467974)
+++ sansmodels/src/sans/models/c_extensions/schulz.c	(revision 2c63e80ef62befde3fa23ebbed7e0649985d92b5)
@@ -10,5 +10,4 @@
 #include "gamma_win.h"
 #endif
-
 
 /**
Index: sansmodels/src/sans/models/c_models/gamma_win.c
===================================================================
--- sansmodels/src/sans/models/c_models/gamma_win.c	(revision 2c63e80ef62befde3fa23ebbed7e0649985d92b5)
+++ sansmodels/src/sans/models/c_models/gamma_win.c	(revision 2c63e80ef62befde3fa23ebbed7e0649985d92b5)
@@ -0,0 +1,151 @@
+// Visit http://www.johndcook.com/stand_alone_code.html for the source of this code and more like it.
+// Modified for the DANSE project.
+
+#include <math>
+
+#include "gamma_win.h"
+
+// Note that the functions Gamma and LogGamma are mutually dependent.
+
+double gamma(double x)
+{
+  // Split the function domain into three intervals:
+  // (0, 0.001), [0.001, 12), and (12, infinity)
+
+  ///////////////////////////////////////////////////////////////////////////
+  // First interval: (0, 0.001)
+  //
+  // For small x, 1/Gamma(x) has power series x + gamma x^2  - ...
+  // So in this range, 1/Gamma(x) = x + gamma x^2 with error on the order of x^3.
+  // The relative error over this interval is less than 6e-7.
+
+  const double gamma = 0.577215664901532860606512090; // Euler's gamma constant
+  double y;
+  int arg_was_less_than_one = 0;
+  int n = 0;
+  // numerator coefficients for approximation over the interval (1,2)
+  static const double p[] =
+  {
+      -1.71618513886549492533811E+0,
+       2.47656508055759199108314E+1,
+      -3.79804256470945635097577E+2,
+       6.29331155312818442661052E+2,
+       8.66966202790413211295064E+2,
+      -3.14512729688483675254357E+4,
+      -3.61444134186911729807069E+4,
+       6.64561438202405440627855E+4
+  };
+  // denominator coefficients for approximation over the interval (1,2)
+  static const double q[] =
+  {
+      -3.08402300119738975254353E+1,
+       3.15350626979604161529144E+2,
+      -1.01515636749021914166146E+3,
+      -3.10777167157231109440444E+3,
+       2.25381184209801510330112E+4,
+       4.75584627752788110767815E+3,
+      -1.34659959864969306392456E+5,
+      -1.15132259675553483497211E+5
+  };
+  double num = 0.0;
+  double den = 1.0;
+  int i;
+  double z;
+  double result;
+  double temp;
+
+  if (x < 0.001)
+      return 1.0/(x*(1.0 + gamma*x));
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Second interval: [0.001, 12)
+
+  if (x < 12.0)
+    {
+      // The algorithm directly approximates gamma over (1,2) and uses
+      // reduction identities to reduce other arguments to this interval.
+
+      y = x;
+
+      if (y < 1.0) arg_was_less_than_one = 1;
+
+      // Add or subtract integers as necessary to bring y into (1,2)
+      // Will correct for this below
+      if (arg_was_less_than_one==1) {
+          y += 1.0;
+      } else {
+          n = static_cast<int> (floor(y)) - 1;  // will use n later
+          y -= n;
+      }
+
+      z = y - 1;
+      for (i = 0; i < 8; i++) {
+          num = (num + p[i])*z;
+          den = den*z + q[i];
+      }
+      result = num/den + 1.0;
+
+      // Apply correction if argument was not initially in (1,2)
+      if (arg_was_less_than_one==1) {
+          // Use identity gamma(z) = gamma(z+1)/z
+          // The variable "result" now holds gamma of the original y + 1
+          // Thus we use y-1 to get back the orginal y.
+          result /= (y-1.0);
+      } else {
+          // Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z)
+          for (i = 0; i < n; i++)
+              result *= y++;
+      }
+
+      return result;
+    }
+
+    ///////////////////////////////////////////////////////////////////////////
+    // Third interval: [12, infinity)
+
+    if (x > 171.624)
+    {
+      // Correct answer too large to display. Force +infinity.
+      temp = DBL_MAX;
+      return temp*2.0;
+    }
+
+    return exp(lgamma(x));
+}
+
+double lgamma (double x)
+{
+  static const double c[8] =
+  {
+   1.0/12.0,
+  -1.0/360.0,
+  1.0/1260.0,
+  -1.0/1680.0,
+  1.0/1188.0,
+  -691.0/360360.0,
+  1.0/156.0,
+  -3617.0/122400.0
+  };
+  double z = 1.0/(x*x);
+  double sum = c[7];
+  static const double halfLogTwoPi = 0.91893853320467274178032973640562;
+  double series;
+  double logGamma;
+
+  if (x < 12.0) {
+     return log(fabs(gamma(x)));
+  }
+
+  // Abramowitz and Stegun 6.1.41
+  // Asymptotic series should be good to at least 11 or 12 figures
+  // For error analysis, see Whittiker and Watson
+  // A Course in Modern Analysis (1927), page 252
+
+  for (int i=6; i >= 0; i--) {
+      sum *= z;
+      sum += c[i];
+  }
+  series = sum/x;
+  logGamma = (x - 0.5)*log(x) - x + halfLogTwoPi + series;
+  return logGamma;
+}
Index: sansmodels/src/sans/models/c_models/gamma_win.cpp
===================================================================
--- sansmodels/src/sans/models/c_models/gamma_win.cpp	(revision 9a137640642d3969a386b2aea613c72d5b738adb)
+++ 	(revision )
@@ -1,177 +1,0 @@
-// Visit http://www.johndcook.com/stand_alone_code.html for the source of this code and more like it.
-
-#include <cmath>
-#include <sstream>
-#include <iostream>
-#include <stdexcept>
-
-#include "gamma_win.h"
-
-// Note that the functions Gamma and LogGamma are mutually dependent.
-
-double gamma
-(
-    double x    // We require x > 0
-)
-{
-  if (x <= 0.0)
-  {
-    std::stringstream os;
-        os << "Invalid input argument " << x <<  ". Argument must be positive.";
-        throw std::invalid_argument( os.str() );
-  }
-
-    // Split the function domain into three intervals:
-    // (0, 0.001), [0.001, 12), and (12, infinity)
-
-    ///////////////////////////////////////////////////////////////////////////
-    // First interval: (0, 0.001)
-  //
-  // For small x, 1/Gamma(x) has power series x + gamma x^2  - ...
-  // So in this range, 1/Gamma(x) = x + gamma x^2 with error on the order of x^3.
-  // The relative error over this interval is less than 6e-7.
-
-  const double gamma = 0.577215664901532860606512090; // Euler's gamma constant
-
-    if (x < 0.001)
-        return 1.0/(x*(1.0 + gamma*x));
-
-    ///////////////////////////////////////////////////////////////////////////
-    // Second interval: [0.001, 12)
-
-  if (x < 12.0)
-    {
-        // The algorithm directly approximates gamma over (1,2) and uses
-        // reduction identities to reduce other arguments to this interval.
-
-    double y = x;
-        int n = 0;
-        bool arg_was_less_than_one = (y < 1.0);
-
-        // Add or subtract integers as necessary to bring y into (1,2)
-        // Will correct for this below
-        if (arg_was_less_than_one)
-        {
-            y += 1.0;
-        }
-        else
-        {
-            n = static_cast<int> (floor(y)) - 1;  // will use n later
-            y -= n;
-        }
-
-        // numerator coefficients for approximation over the interval (1,2)
-        static const double p[] =
-        {
-            -1.71618513886549492533811E+0,
-             2.47656508055759199108314E+1,
-            -3.79804256470945635097577E+2,
-             6.29331155312818442661052E+2,
-             8.66966202790413211295064E+2,
-            -3.14512729688483675254357E+4,
-            -3.61444134186911729807069E+4,
-             6.64561438202405440627855E+4
-        };
-
-        // denominator coefficients for approximation over the interval (1,2)
-        static const double q[] =
-        {
-            -3.08402300119738975254353E+1,
-             3.15350626979604161529144E+2,
-            -1.01515636749021914166146E+3,
-            -3.10777167157231109440444E+3,
-             2.25381184209801510330112E+4,
-             4.75584627752788110767815E+3,
-            -1.34659959864969306392456E+5,
-            -1.15132259675553483497211E+5
-        };
-
-        double num = 0.0;
-        double den = 1.0;
-        int i;
-
-        double z = y - 1;
-        for (i = 0; i < 8; i++)
-        {
-            num = (num + p[i])*z;
-            den = den*z + q[i];
-        }
-        double result = num/den + 1.0;
-
-        // Apply correction if argument was not initially in (1,2)
-        if (arg_was_less_than_one)
-        {
-            // Use identity gamma(z) = gamma(z+1)/z
-            // The variable "result" now holds gamma of the original y + 1
-            // Thus we use y-1 to get back the orginal y.
-            result /= (y-1.0);
-        }
-        else
-        {
-            // Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z)
-            for (i = 0; i < n; i++)
-                result *= y++;
-        }
-
-    return result;
-    }
-
-    ///////////////////////////////////////////////////////////////////////////
-    // Third interval: [12, infinity)
-
-    if (x > 171.624)
-    {
-    // Correct answer too large to display. Force +infinity.
-    double temp = DBL_MAX;
-    return temp*2.0;
-    }
-
-    return exp(lgamma(x));
-}
-
-double lgamma
-(
-    double x    // x must be positive
-)
-{
-  if (x <= 0.0)
-  {
-    std::stringstream os;
-        os << "Invalid input argument " << x <<  ". Argument must be positive.";
-        throw std::invalid_argument( os.str() );
-  }
-
-    if (x < 12.0)
-    {
-        return log(fabs(gamma(x)));
-    }
-
-  // Abramowitz and Stegun 6.1.41
-    // Asymptotic series should be good to at least 11 or 12 figures
-    // For error analysis, see Whittiker and Watson
-    // A Course in Modern Analysis (1927), page 252
-
-    static const double c[8] =
-    {
-     1.0/12.0,
-    -1.0/360.0,
-    1.0/1260.0,
-    -1.0/1680.0,
-    1.0/1188.0,
-    -691.0/360360.0,
-    1.0/156.0,
-    -3617.0/122400.0
-    };
-    double z = 1.0/(x*x);
-    double sum = c[7];
-    for (int i=6; i >= 0; i--)
-    {
-        sum *= z;
-        sum += c[i];
-    }
-    double series = sum/x;
-
-    static const double halfLogTwoPi = 0.91893853320467274178032973640562;
-    double logGamma = (x - 0.5)*log(x) - x + halfLogTwoPi + series;
-  return logGamma;
-}
Index: sansmodels/src/sans/models/c_models/gamma_win.h
===================================================================
--- sansmodels/src/sans/models/c_models/gamma_win.h	(revision 83bf68e1b7cb622a72c76c9ab5405abc3116da24)
+++ sansmodels/src/sans/models/c_models/gamma_win.h	(revision 2c63e80ef62befde3fa23ebbed7e0649985d92b5)
@@ -5,7 +5,15 @@
 
 // Note that the functions Gamma and LogGamma are mutually dependent.
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 double lgamma(double);
 double gamma(double);
 
+#ifdef __cplusplus
+}
 #endif
 
+#endif
+
