Added some documentation to the numerical integration classes and applied small cosmetic changes to msr2data.

src/classes/PMsr2Data.cpp

@@ -1959,7 +1959,7 @@ int PMsr2Data::WriteOutput(const string &outfile, bool db, unsigned int withHead
         outFile.seekg(-i, ios::end);
         getline(outFile, s);
         trim(s); // remove whitespace
-        if (s.empty() || (s == "\n")) { // (s == "\n") check is for M$-systems using "\r\n" linebreaks
+        if (s.empty()) { // trim cuts off also characters like '\n', therefore this should work also with M$-DOS linebreaks
           if (i == size) {
             outFile.seekp(0);
             fHeaderWritten = false;  // if the file contained only empty lines, default to writing the header

src/external/BMWIntegrator/BMWIntegrator.cpp

@@ -5,7 +5,7 @@
   Author: Bastian M. Wojek
   e-mail: bastian.wojek@psi.ch
 
-  $Id$
+  $Id
 
 ***************************************************************************/
 
@@ -38,6 +38,12 @@
 
 std::vector<double> TDWaveGapIntegralCuhre::fPar;
 
+/**
+ * <p>Integrate the function using the Cuhre interface
+ *
+ * <p><b>return:</b>
+ * - value of the integral
+ */
 double TDWaveGapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
@@ -61,6 +67,18 @@ double TDWaveGapIntegralCuhre::IntegrateFunc()
   return integral[0];
 }
 
+/**
+ * <p>Calculate the function value for the use with Cuhre---actual implementation of the function
+ *
+ * <p><b>return:</b>
+ * - 0
+ *
+ * \param ndim number of dimensions of the integral (2 here)
+ * \param x point where the function should be evaluated
+ * \param ncomp number of components of the integrand (1 here)
+ * \param f function value
+ * \param userdata additional user parameters (required by the interface, NULL here)
+ */
 int TDWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
                       const int *ncomp, double f[], void *userdata) // x = {E, phi}, fPar = {twokBT, Delta(T), Ec, phic}
 {
@@ -71,6 +89,12 @@ int TDWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
 
 std::vector<double> TCosSqDWaveGapIntegralCuhre::fPar;
 
+/**
+ * <p>Integrate the function using the Cuhre interface
+ *
+ * <p><b>return:</b>
+ * - value of the integral
+ */
 double TCosSqDWaveGapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
@@ -94,6 +118,18 @@ double TCosSqDWaveGapIntegralCuhre::IntegrateFunc()
   return integral[0];
 }
 
+/**
+ * <p>Calculate the function value for the use with Cuhre---actual implementation of the function
+ *
+ * <p><b>return:</b>
+ * - 0
+ *
+ * \param ndim number of dimensions of the integral (2 here)
+ * \param x point where the function should be evaluated
+ * \param ncomp number of components of the integrand (1 here)
+ * \param f function value
+ * \param userdata additional user parameters (required by the interface, NULL here)
+ */
 int TCosSqDWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
                       const int *ncomp, double f[], void *userdata) // x = {E, phi}, fPar = {twokBT, DeltaD(T), Ec, phic, DeltaS(T)}
 {
@@ -104,6 +140,12 @@ int TCosSqDWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
 
 std::vector<double> TSinSqDWaveGapIntegralCuhre::fPar;
 
+/**
+ * <p>Integrate the function using the Cuhre interface
+ *
+ * <p><b>return:</b>
+ * - value of the integral
+ */
 double TSinSqDWaveGapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
@@ -127,6 +169,18 @@ double TSinSqDWaveGapIntegralCuhre::IntegrateFunc()
   return integral[0];
 }
 
+/**
+ * <p>Calculate the function value for the use with Cuhre---actual implementation of the function
+ *
+ * <p><b>return:</b>
+ * - 0
+ *
+ * \param ndim number of dimensions of the integral (2 here)
+ * \param x point where the function should be evaluated
+ * \param ncomp number of components of the integrand (1 here)
+ * \param f function value
+ * \param userdata additional user parameters (required by the interface, NULL here)
+ */
 int TSinSqDWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
                       const int *ncomp, double f[], void *userdata) // x = {E, phi}, fPar = {twokBT, DeltaD(T), Ec, phic, DeltaS(T)}
 {
@@ -135,9 +189,14 @@ int TSinSqDWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
   return 0;
 }
 
-
 std::vector<double> TAnSWaveGapIntegralCuhre::fPar;
 
+/**
+ * <p>Integrate the function using the Cuhre interface
+ *
+ * <p><b>return:</b>
+ * - value of the integral
+ */
 double TAnSWaveGapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
@@ -161,6 +220,18 @@ double TAnSWaveGapIntegralCuhre::IntegrateFunc()
   return integral[0];
 }
 
+/**
+ * <p>Calculate the function value for the use with Cuhre---actual implementation of the function
+ *
+ * <p><b>return:</b>
+ * - 0
+ *
+ * \param ndim number of dimensions of the integral (2 here)
+ * \param x point where the function should be evaluated
+ * \param ncomp number of components of the integrand (1 here)
+ * \param f function value
+ * \param userdata additional user parameters (required by the interface, NULL here)
+ */
 int TAnSWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
                       const int *ncomp, double f[], void *userdata) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
 {
@@ -171,6 +242,12 @@ int TAnSWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
 
 std::vector<double> TAnSWaveGapIntegralDivonne::fPar;
 
+/**
+ * <p>Integrate the function using the Divonne interface
+ *
+ * <p><b>return:</b>
+ * - value of the integral
+ */
 double TAnSWaveGapIntegralDivonne::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
@@ -202,6 +279,18 @@ double TAnSWaveGapIntegralDivonne::IntegrateFunc()
   return integral[0];
 }
 
+/**
+ * <p>Calculate the function value for the use with Divonne---actual implementation of the function
+ *
+ * <p><b>return:</b>
+ * - 0
+ *
+ * \param ndim number of dimensions of the integral (2 here)
+ * \param x point where the function should be evaluated
+ * \param ncomp number of components of the integrand (1 here)
+ * \param f function value
+ * \param userdata additional user parameters (required by the interface, NULL here)
+ */
 int TAnSWaveGapIntegralDivonne::Integrand(const int *ndim, const double x[],
                       const int *ncomp, double f[], void *userdata) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
 {
@@ -212,6 +301,12 @@ int TAnSWaveGapIntegralDivonne::Integrand(const int *ndim, const double x[],
 
 std::vector<double> TAnSWaveGapIntegralSuave::fPar;
 
+/**
+ * <p>Integrate the function using the Suave interface
+ *
+ * <p><b>return:</b>
+ * - value of the integral
+ */
 double TAnSWaveGapIntegralSuave::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
@@ -236,6 +331,18 @@ double TAnSWaveGapIntegralSuave::IntegrateFunc()
   return integral[0];
 }
 
+/**
+ * <p>Calculate the function value for the use with Suave---actual implementation of the function
+ *
+ * <p><b>return:</b>
+ * - 0
+ *
+ * \param ndim number of dimensions of the integral (2 here)
+ * \param x point where the function should be evaluated
+ * \param ncomp number of components of the integrand (1 here)
+ * \param f function value
+ * \param userdata additional user parameters (required by the interface, NULL here)
+ */
 int TAnSWaveGapIntegralSuave::Integrand(const int *ndim, const double x[],
                       const int *ncomp, double f[], void *userdata) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
 {
@@ -246,6 +353,12 @@ int TAnSWaveGapIntegralSuave::Integrand(const int *ndim, const double x[],
 
 std::vector<double> TNonMonDWave1GapIntegralCuhre::fPar;
 
+/**
+ * <p>Integrate the function using the Cuhre interface
+ *
+ * <p><b>return:</b>
+ * - value of the integral
+ */
 double TNonMonDWave1GapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
@@ -269,6 +382,18 @@ double TNonMonDWave1GapIntegralCuhre::IntegrateFunc()
   return integral[0];
 }
 
+/**
+ * <p>Calculate the function value for the use with Cuhre---actual implementation of the function
+ *
+ * <p><b>return:</b>
+ * - 0
+ *
+ * \param ndim number of dimensions of the integral (2 here)
+ * \param x point where the function should be evaluated
+ * \param ncomp number of components of the integrand (1 here)
+ * \param f function value
+ * \param userdata additional user parameters (required by the interface, NULL here)
+ */
 int TNonMonDWave1GapIntegralCuhre::Integrand(const int *ndim, const double x[],
                       const int *ncomp, double f[], void *userdata) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
 {
@@ -279,6 +404,12 @@ int TNonMonDWave1GapIntegralCuhre::Integrand(const int *ndim, const double x[],
 
 std::vector<double> TNonMonDWave2GapIntegralCuhre::fPar;
 
+/**
+ * <p>Integrate the function using the Cuhre interface
+ *
+ * <p><b>return:</b>
+ * - value of the integral
+ */
 double TNonMonDWave2GapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
@@ -302,6 +433,18 @@ double TNonMonDWave2GapIntegralCuhre::IntegrateFunc()
   return integral[0];
 }
 
+/**
+ * <p>Calculate the function value for the use with Cuhre---actual implementation of the function
+ *
+ * <p><b>return:</b>
+ * - 0
+ *
+ * \param ndim number of dimensions of the integral (2 here)
+ * \param x point where the function should be evaluated
+ * \param ncomp number of components of the integrand (1 here)
+ * \param f function value
+ * \param userdata additional user parameters (required by the interface, NULL here)
+ */
 int TNonMonDWave2GapIntegralCuhre::Integrand(const int *ndim, const double x[],
                       const int *ncomp, double f[], void *userdata) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
 {

src/external/BMWIntegrator/BMWIntegrator.h

@@ -5,7 +5,7 @@
   Author: Bastian M. Wojek
   e-mail: bastian.wojek@psi.ch
 
-  $Id:
+  $Id
 
 ***************************************************************************/
 
@@ -41,8 +41,11 @@
 
 using namespace std;
 
-// 1D Integrator base class - the function to be integrated have to be implemented in a derived class
+/**
-
+ * <p>Base class for 1D integrations using the GNU Scientific Library integrator.
+ *    The function which should be integrated has to be implemented in a derived class.
+ *    Note: The purpose of this is to offer an easy-to-use interface---not the most efficient integration routine.
+ */
 class TIntegrator {
   public:
     TIntegrator();
@@ -52,37 +55,67 @@ class TIntegrator {
     double IntegrateFunc(double, double);
 
   protected:
-    mutable vector<double> fPar;
+    mutable vector<double> fPar; ///< parameters of the integrand
 
   private:
     static double FuncAtXgsl(double, void *);
-    ROOT::Math::GSLIntegrator *fIntegrator;
+    ROOT::Math::GSLIntegrator *fIntegrator; ///< pointer to the GSL integrator
-    mutable double (*fFunc)(double, void *);
+    mutable double (*fFunc)(double, void *); ///< pointer to the integrand function
 };
 
+/**
+ * <p>Constructor of the base class for 1D integrations
+ *    Allocation of memory for an integration using the adaptive 31 point Gauss-Kronrod rule
+ */
 inline TIntegrator::TIntegrator() : fFunc(0) {
   fIntegrator = new ROOT::Math::GSLIntegrator(ROOT::Math::Integration::kADAPTIVE,ROOT::Math::Integration::kGAUSS31);
 }
 
+/**
+ * <p>Destructor of the base class for 1D integrations
+ *    Clean up.
+ */
 inline TIntegrator::~TIntegrator(){
+  fPar.clear();
   delete fIntegrator;
   fIntegrator=0;
   fFunc=0;
 }
 
+/**
+ * <p>Method for passing the integrand function value to the integrator.
+ *
+ * <p><b>return:</b>
+ * - function value of the integrand
+ *
+ * \param x point at which the function value is calculated
+ * \param obj pointer to the integrator
+ */
 inline double TIntegrator::FuncAtXgsl(double x, void *obj)
 {
   return ((TIntegrator*)obj)->FuncAtX(x);
 }
 
+/**
+ * <p>Calculate the integral of the function between the given boundaries
+ *
+ * <p><b>return:</b>
+ * - value of the integral
+ *
+ * \param x1 lower boundary
+ * \param x2 upper boundary
+ */
 inline double TIntegrator::IntegrateFunc(double x1, double x2)
 {
   fFunc = &TIntegrator::FuncAtXgsl;
   return fIntegrator->Integral(fFunc, (this), x1, x2);
 }
 
-// Multi dimensional GSL Monte Carlo Integrations
+/**
-
+ * <p>Base class for multidimensional Monte-Carlo integrations using the GNU Scientific Library integrator.
+ *    The function which should be integrated has to be implemented in a derived class.
+ *    Note: The purpose of this is to offer an easy-to-use interface---not the most efficient integration routine.
+ */
 class TMCIntegrator {
   public:
     TMCIntegrator();
@@ -92,37 +125,69 @@ class TMCIntegrator {
     double IntegrateFunc(size_t, double *, double *);
 
   protected:
-    mutable vector<double> fPar;
+    mutable vector<double> fPar; ///< parameters of the integrand
 
   private:
     static double FuncAtXgsl(double *, size_t, void *);
-    ROOT::Math::GSLMCIntegrator *fMCIntegrator;
+    ROOT::Math::GSLMCIntegrator *fMCIntegrator; ///< pointer to the GSL integrator
-    mutable double (*fFunc)(double *, size_t, void *);
+    mutable double (*fFunc)(double *, size_t, void *); ///< pointer to the integrand function
 };
 
+/**
+ * <p>Constructor of the base class for multidimensional Monte-Carlo integrations
+ *    Allocation of memory for an integration using the MISER algorithm of Press and Farrar
+ */
 inline TMCIntegrator::TMCIntegrator() : fFunc(0) {
   fMCIntegrator = new ROOT::Math::GSLMCIntegrator(ROOT::Math::MCIntegration::kMISER, 1.E-6, 1.E-4, 500000);
 }
 
+/**
+ * <p>Destructor of the base class for 1D integrations
+ *    Clean up.
+ */
 inline TMCIntegrator::~TMCIntegrator(){
+  fPar.clear();
   delete fMCIntegrator;
   fMCIntegrator=0;
   fFunc=0;
 }
 
+/**
+ * <p>Method for passing the integrand function value to the integrator.
+ *
+ * <p><b>return:</b>
+ * - function value of the integrand
+ *
+ * \param x point at which the function value is calculated
+ * \param dim number of dimensions
+ * \param obj pointer to the integrator
+ */
 inline double TMCIntegrator::FuncAtXgsl(double *x, size_t dim, void *obj)
 {
   return ((TMCIntegrator*)obj)->FuncAtX(x);
 }
 
+/**
+ * <p>Calculate the integral of the function between the given boundaries
+ *
+ * <p><b>return:</b>
+ * - value of the integral
+ *
+ * \param dim number of dimensions
+ * \param x1 lower boundary array
+ * \param x2 upper boundary array
+ */
 inline double TMCIntegrator::IntegrateFunc(size_t dim, double *x1, double *x2)
 {
   fFunc = &TMCIntegrator::FuncAtXgsl;
   return fMCIntegrator->Integral(fFunc, dim, x1, x2, (this));
 }
 
-// Multidimensional Integrator class for a d-wave gap integral using the Cuhre algorithm
+/**
-
+ * <p>Two-dimensional integrator class for the efficient calculation of the superfluid density within the semi-classical model
+ *    assuming a cylindrical Fermi surface and a d_{x^2-y^2} symmetry of the superconducting order parameter.
+ *    The integration uses the Cuhre algorithm of the Cuba library.
+ */
 class TDWaveGapIntegralCuhre {
   public:
     TDWaveGapIntegralCuhre() : fNDim(2) {}
@@ -132,11 +197,16 @@ class TDWaveGapIntegralCuhre {
     double IntegrateFunc();
 
   protected:
-    static vector<double> fPar;
+    static vector<double> fPar; ///< parameters of the integrand
-    unsigned int fNDim;
+    unsigned int fNDim; ///< dimension of the integral
-
 };
 
+/**
+ * <p>Two-dimensional integrator class for the efficient calculation of the superfluid density along the a-axis
+ *    within the semi-classical model assuming a cylindrical Fermi surface and a mixed d_{x^2-y^2} + s symmetry of the
+ *    superconducting order parameter (effectively: d_{x^2-y^2} with shifted nodes and a-b-anisotropy).
+ *    The integration uses the Cuhre algorithm of the Cuba library.
+ */
 class TCosSqDWaveGapIntegralCuhre {
   public:
     TCosSqDWaveGapIntegralCuhre() : fNDim(2) {}
@@ -146,11 +216,16 @@ class TCosSqDWaveGapIntegralCuhre {
     double IntegrateFunc();
 
   protected:
-    static vector<double> fPar;
+    static vector<double> fPar; ///< parameters of the integrand
-    unsigned int fNDim;
+    unsigned int fNDim; ///< dimension of the integral
-
 };
 
+/**
+ * <p>Two-dimensional integrator class for the efficient calculation of the superfluid density along the b-axis
+ *    within the semi-classical model assuming a cylindrical Fermi surface and a mixed d_{x^2-y^2} + s symmetry of the
+ *    superconducting order parameter (effectively: d_{x^2-y^2} with shifted nodes and a-b-anisotropy).
+ *    The integration uses the Cuhre algorithm of the Cuba library.
+ */
 class TSinSqDWaveGapIntegralCuhre {
   public:
     TSinSqDWaveGapIntegralCuhre() : fNDim(2) {}
@@ -160,11 +235,15 @@ class TSinSqDWaveGapIntegralCuhre {
     double IntegrateFunc();
 
   protected:
-    static vector<double> fPar;
+    static vector<double> fPar; ///< parameters of the integrand
-    unsigned int fNDim;
+    unsigned int fNDim; ///< dimension of the integral
-
 };
 
+/**
+ * <p>Two-dimensional integrator class for the efficient calculation of the superfluid density within the semi-classical model
+ *    assuming a cylindrical Fermi surface and an "anisotropic s-wave" symmetry of the superconducting order parameter.
+ *    The integration uses the Cuhre algorithm of the Cuba library.
+ */
 class TAnSWaveGapIntegralCuhre {
   public:
     TAnSWaveGapIntegralCuhre() : fNDim(2) {}
@@ -174,11 +253,15 @@ class TAnSWaveGapIntegralCuhre {
     double IntegrateFunc();
 
   protected:
-    static vector<double> fPar;
+    static vector<double> fPar; ///< parameters of the integrand
-    unsigned int fNDim;
+    unsigned int fNDim; ///< dimension of the integral
-
 };
 
+/**
+ * <p>Two-dimensional integrator class for the efficient calculation of the superfluid density within the semi-classical model
+ *    assuming a cylindrical Fermi surface and an "anisotropic s-wave" symmetry of the superconducting order parameter.
+ *    The integration uses the Divonne algorithm of the Cuba library.
+ */
 class TAnSWaveGapIntegralDivonne {
   public:
     TAnSWaveGapIntegralDivonne() : fNDim(2) {}
@@ -188,11 +271,15 @@ class TAnSWaveGapIntegralDivonne {
     double IntegrateFunc();
 
   protected:
-    static vector<double> fPar;
+    static vector<double> fPar; ///< parameters of the integrand
-    unsigned int fNDim;
+    unsigned int fNDim; ///< dimension of the integral
-
 };
 
+/**
+ * <p>Two-dimensional integrator class for the efficient calculation of the superfluid density within the semi-classical model
+ *    assuming a cylindrical Fermi surface and an "anisotropic s-wave" symmetry of the superconducting order parameter.
+ *    The integration uses the Suave algorithm of the Cuba library.
+ */
 class TAnSWaveGapIntegralSuave {
   public:
     TAnSWaveGapIntegralSuave() : fNDim(2) {}
@@ -202,11 +289,15 @@ class TAnSWaveGapIntegralSuave {
     double IntegrateFunc();
 
   protected:
-    static vector<double> fPar;
+    static vector<double> fPar; ///< parameters of the integrand
-    unsigned int fNDim;
+    unsigned int fNDim; ///< dimension of the integral
-
 };
 
+/**
+ * <p>Two-dimensional integrator class for the efficient calculation of the superfluid density within the semi-classical model
+ *    assuming a cylindrical Fermi surface and an "non-monotonic d-wave" symmetry of the superconducting order parameter.
+ *    The integration uses the Cuhre algorithm of the Cuba library.
+ */
 class TNonMonDWave1GapIntegralCuhre {
   public:
     TNonMonDWave1GapIntegralCuhre() : fNDim(2) {}
@@ -216,11 +307,15 @@ class TNonMonDWave1GapIntegralCuhre {
     double IntegrateFunc();
 
   protected:
-    static vector<double> fPar;
+    static vector<double> fPar; ///< parameters of the integrand
-    unsigned int fNDim;
+    unsigned int fNDim; ///< dimension of the integral
-
 };
 
+/**
+ * <p>Two-dimensional integrator class for the efficient calculation of the superfluid density within the semi-classical model
+ *    assuming a cylindrical Fermi surface and an "non-monotonic d-wave" symmetry of the superconducting order parameter.
+ *    The integration uses the Cuhre algorithm of the Cuba library.
+ */
 class TNonMonDWave2GapIntegralCuhre {
   public:
     TNonMonDWave2GapIntegralCuhre() : fNDim(2) {}
@@ -230,13 +325,14 @@ class TNonMonDWave2GapIntegralCuhre {
     double IntegrateFunc();
 
   protected:
-    static vector<double> fPar;
+    static vector<double> fPar; ///< parameters of the integrand
-    unsigned int fNDim;
+    unsigned int fNDim; ///< dimension of the integral
-
 };
 
-// To be integrated: x*y dx dy
+/**
-
+ * <p>Test class for the 2D MC integration
+ *    Integral: x*y dx dy
+ */
 class T2DTest : public TMCIntegrator {
   public:
     T2DTest() {}
@@ -244,13 +340,24 @@ class T2DTest : public TMCIntegrator {
     double FuncAtX(double *) const;
 };
 
+/**
+ * <p>Calculate the function value---actual implementation of the function x*y
+ *
+ * <p><b>return:</b>
+ * - function value
+ *
+ * \param x point where the function should be evaluated
+ */
 inline double T2DTest::FuncAtX(double *x) const
 {
   return x[0]*x[1];
 }
 
-// To be integrated: d wave gap integral
+/**
-
+ * <p>Class for the 2D Monte-Carlo integration for the calculation of the superfluid density within the semi-classical model
+ *    assuming a cylindrical Fermi surface and a d_{x^2-y^2} symmetry of the superconducting order parameter.
+ *    The integration uses the GSL integration routines.
+ */
 class TDWaveGapIntegral : public TMCIntegrator {
   public:
     TDWaveGapIntegral() {}
@@ -258,6 +365,14 @@ class TDWaveGapIntegral : public TMCIntegrator {
     double FuncAtX(double *) const;
 };
 
+/**
+ * <p>Calculate the function value---actual implementation of the function
+ *
+ * <p><b>return:</b>
+ * - function value
+ *
+ * \param x point where the function should be evaluated
+ */
 inline double TDWaveGapIntegral::FuncAtX(double *x) const // x = {E, phi}, fPar = {T, Delta(T)}
 {
   double twokt(2.0*0.08617384436*fPar[0]); // kB in meV/K
@@ -265,8 +380,11 @@ inline double TDWaveGapIntegral::FuncAtX(double *x) const // x = {E, phi}, fPar
   return -1.0/(2.0*twokt*TMath::CosH(TMath::Sqrt(x[0]*x[0]+deltasq)/twokt)*TMath::CosH(TMath::Sqrt(x[0]*x[0]+deltasq)/twokt));
 }
 
-// To be integrated: anisotropic s wave gap integral
+/**
-
+ * <p>Class for the 2D Monte-Carlo integration for the calculation of the superfluid density within the semi-classical model
+ *    assuming a cylindrical Fermi surface and an "anisotropic s-wave" symmetry of the superconducting order parameter.
+ *    The integration uses the GSL integration routines.
+ */
 class TAnSWaveGapIntegral : public TMCIntegrator {
   public:
     TAnSWaveGapIntegral() {}
@@ -274,6 +392,14 @@ class TAnSWaveGapIntegral : public TMCIntegrator {
     double FuncAtX(double *) const;
 };
 
+/**
+ * <p>Calculate the function value---actual implementation of the function
+ *
+ * <p><b>return:</b>
+ * - function value
+ *
+ * \param x point where the function should be evaluated
+ */
 inline double TAnSWaveGapIntegral::FuncAtX(double *x) const // x = {E, phi}, fPar = {T, Delta(T), a}
 {
   double twokt(2.0*0.08617384436*fPar[0]); // kB in meV/K
@@ -281,8 +407,10 @@ inline double TAnSWaveGapIntegral::FuncAtX(double *x) const // x = {E, phi}, fPa
   return -1.0/(2.0*twokt*TMath::CosH(TMath::Sqrt(x[0]*x[0]+deltasq)/twokt)*TMath::CosH(TMath::Sqrt(x[0]*x[0]+deltasq)/twokt));
 }
 
-// To be integrated: Bessel function times Exponential
+/**
-
+ * <p>Class for the 1D integration of j0(a*x)*exp(-b*x)
+ *    The integration uses the GSL integration routines.
+ */
 class TIntBesselJ0Exp : public TIntegrator {
   public:
     TIntBesselJ0Exp() {}
@@ -290,6 +418,14 @@ class TIntBesselJ0Exp : public TIntegrator {
     double FuncAtX(double) const;
 };
 
+/**
+ * <p>Calculate the function value---actual implementation of the function j0(a*x)*exp(-b*x)
+ *
+ * <p><b>return:</b>
+ * - function value
+ *
+ * \param x point where the function should be evaluated
+ */
 inline double TIntBesselJ0Exp::FuncAtX(double x) const
 {
   double w0t(TMath::TwoPi()*fPar[0]*x);
@@ -302,8 +438,10 @@ inline double TIntBesselJ0Exp::FuncAtX(double x) const
   return j0 * TMath::Exp(-fPar[1]*x);
 }
 
-// To be integrated: Sine times Gaussian
+/**
-
+ * <p>Class for the 1D integration of sin(a*x)*exp(-b*x*x)
+ *    The integration uses the GSL integration routines.
+ */
 class TIntSinGss : public TIntegrator {
   public:
     TIntSinGss() {}
@@ -311,13 +449,25 @@ class TIntSinGss : public TIntegrator {
     double FuncAtX(double) const;
 };
 
+/**
+ * <p>Calculate the function value---actual implementation of the function sin(a*x)*exp(-b*x*x)
+ *
+ * <p><b>return:</b>
+ * - function value
+ *
+ * \param x point where the function should be evaluated
+ */
 inline double TIntSinGss::FuncAtX(double x) const
 {
   return TMath::Sin(TMath::TwoPi()*fPar[0]*x) * TMath::Exp(-0.5*fPar[1]*fPar[1]*x*x);
 }
 
-// To be integrated: DeRenzi Spin Glass Interpolation Integrand
+/**
-
+ * <p>Class for the 1D integration of the "DeRenzi Spin Glass Interpolation Integrand"
+ *    See Eq. (5) of R. De Renzi and S. Fanesi, Physica B 289-290, 209-212 (2000).
+ *    doi:10.1016/S0921-4526(00)00368-9
+ *    The integration uses the GSL integration routines.
+ */
 class TIntSGInterpolation : public TIntegrator {
   public:
     TIntSGInterpolation() {}
@@ -325,6 +475,14 @@ class TIntSGInterpolation : public TIntegrator {
     double FuncAtX(double) const;
 };
 
+/**
+ * <p>Calculate the function value---actual implementation of the function
+ *
+ * <p><b>return:</b>
+ * - function value
+ *
+ * \param x point where the function should be evaluated
+ */
 inline double TIntSGInterpolation::FuncAtX(double x) const
 {
   // Parameters: nu_L [MHz], a [1/us], lambda [1/us], beta [1], t [us]
@@ -333,17 +491,26 @@ inline double TIntSGInterpolation::FuncAtX(double x) const
   return (wt*TMath::Cos(wt)-TMath::Sin(wt))/(wt*wt)*TMath::Exp(-TMath::Power(expo,fPar[3]))/TMath::Power(expo,(1.0-fPar[3]));
 }
 
-
+/**
-// To be integrated: df/dE * E / sqrt(E^2 - Delta^2)
+ * <p>Class for the 1D integration for the calculation of the superfluid density within the semi-classical model
-
+ *    assuming a cylindrical Fermi surface and an isotropic s-wave symmetry of the superconducting order parameter.
+ *    The integration uses the GSL integration routines.
+ */
 class TGapIntegral : public TIntegrator {
   public:
     TGapIntegral() {}
     ~TGapIntegral() {}
     double FuncAtX(double) const; // variable: E
-
 };
 
+/**
+ * <p>Calculate the function value---actual implementation of the function df/dE * E / sqrt(E^2 - Delta^2)
+ *
+ * <p><b>return:</b>
+ * - function value
+ *
+ * \param x point where the function should be evaluated
+ */
 inline double TGapIntegral::FuncAtX(double e) const
 {
   return 1.0/(TMath::Power(TMath::CosH(TMath::Sqrt(e*e+fPar[1]*fPar[1])/fPar[0]),2.0));

src/msr2data.cpp

@@ -90,7 +90,7 @@ void msr2data_syntax()
   cout << endl;
   cout << endl << "       <run>, <run1>, <run2>, ... <runN> : run numbers";
   cout << endl << "       <extension> : msr-file extension, e.g. _tf_h13 for the file name 8472_tf_h13.msr";
-  cout << endl << "       -o<outputfile> : specify the name of the DB or column data output file; default: out.db/out.dat";
+  cout << endl << "       -o<outputfile> : specify the name of the DB or column-data output file; default: out.db/out.dat";
   cout << endl << "                        if the option '-o none' is used, no output file will be written.";
   cout << endl << "       new : before writing a new output file, delete the contents of any existing file with the same name";
   cout << endl << "       data : instead of to a DB file the data are written to a simple column structure";
@@ -101,25 +101,25 @@ void msr2data_syntax()
   cout << endl << "       nosummary : no additional data from the run data file is written to the output file";
   cout << endl << "       fit : invoke musrfit to fit the specified runs";
   cout << endl << "              All msr input files are assumed to be present, none is newly generated!";
-  cout << endl << "       fit-<template>! : generate msr-files for the runs to be processed from the <template>-run";
+  cout << endl << "       fit-<template>! : generate msr files for the runs to be processed from the <template> run";
   cout << endl << "              and call musrfit for fitting these runs";
-  cout << endl << "       fit-<template> : same as above, but the <template>-run is only used for the first file creation - ";
+  cout << endl << "       fit-<template> : same as above, but the <template> run is only used for the first file creation---";
-  cout << endl << "              the succeding files are generated using the musrfit-output from the preceding runs";
+  cout << endl << "              the successive files are generated using the musrfit output from the preceding runs";
   cout << endl << "       msr-<template> : same as above without calling musrfit";
-  cout << endl << "              In case any fitting-option is present, this option is ignored!";
+  cout << endl << "              In case any fitting option is present, this option is ignored!";
   cout << endl << "       -k : if fitting is used, pass the option --keep-mn2-output to musrfit";
   cout << endl << "       -t : if fitting is used, pass the option --title-from-data-file to musrfit";
   cout << endl;
   cout << endl << "       global : switch on the global-fit mode";
-  cout << endl << "              Within that mode all specified runs will be united in a single msr-file!";
+  cout << endl << "              Within that mode all specified runs will be united in a single msr file!";
-  cout << endl << "              The fit parameters can be either run-specific or common to all runs.";
+  cout << endl << "              The fit parameters can be either run specific or common to all runs.";
   cout << endl << "              For a complete description of this feature please refer to the manual.";
   cout << endl;
   cout << endl << "       global+[!] : operate in the global-fit mode, however, in case a global input file is created";
   cout << endl << "              all specified runs are pre-analyzed first one by one using the given template.";
   cout << endl << "              For the generation of the global input file, the run-specific parameter values are taken";
-  cout << endl << "              from this pre-analysis for each run - not just copied from the template.";
+  cout << endl << "              from this pre-analysis for each run---they are not just copied from the template.";
-  cout << endl << "              The specification of ! determines which fit-mode (see above) is used for this pre-analysis.";
+  cout << endl << "              The specification of '!' determines which fit mode (see above) is used for this pre-analysis.";
   cout << endl;
   cout << endl << "    For further information please refer to";
   cout << endl << "    https://intranet.psi.ch/MUSR/Msr2Data";
@@ -624,11 +624,14 @@ int main(int argc, char *argv[])
     // delete old db/data file if the "new" option is given
     if (!msr2data_useOption(arg, "new")) {
       fstream fileOutput;
-      fileOutput.open(outputFile.c_str(), ios::in | ios::out | ios::trunc);
+      fileOutput.open(outputFile.c_str(), ios::in);
       if (fileOutput.is_open()) {
-        cout << endl << ">> msr2data: **INFO** Deleting output file " << outputFile << " if it existed" << endl;
+        cout << endl << ">> msr2data: **INFO** Deleting output file " << outputFile << endl;
-        //fileOutput << endl;
         fileOutput.close();
+        fileOutput.open(outputFile.c_str(), ios::out | ios::trunc);
+        fileOutput.close();
+      } else {
+        cout << endl << ">> msr2data: **INFO** Ignoring the 'new' option since " << outputFile << " does not exist yet." << endl;
       }
       if (writeHeader == 2) {
         writeHeader = 1;