changed temperature depends of the gap functions, added comments and amended the documentation

src/external/libGapIntegrals/GapIntegrals.tex

@@ -0,0 +1,198 @@
+\documentclass[twoside]{article}
+
+\usepackage[english]{babel}
+%\usepackage{a4}
+\usepackage{amssymb,amsmath,bm}
+\usepackage{graphicx,tabularx}
+\usepackage{fancyhdr}
+\usepackage{array}
+\usepackage{float}
+\usepackage{hyperref}
+\usepackage{xspace}
+\usepackage{rotating}
+\usepackage{dcolumn}
+\usepackage{geometry}
+
+\geometry{a4paper,left=20mm,right=20mm,top=20mm,bottom=20mm}
+
+% \setlength{\topmargin}{10mm}
+% \setlength{\topmargin}{-13mm}
+% % \setlength{\oddsidemargin}{0.5cm}
+% % \setlength{\evensidemargin}{0cm}
+% \setlength{\oddsidemargin}{1cm}
+% \setlength{\evensidemargin}{1cm}
+% \setlength{\textwidth}{15cm}
+\setlength{\textheight}{23.8cm}
+
+\pagestyle{fancyplain}
+\addtolength{\headwidth}{0.6cm}
+\fancyhead{}%
+\fancyhead[RE,LO]{\bf \textsc{GapIntegrals}}%
+\fancyhead[LE,RO]{\thepage}
+\cfoot{--- B.M.~Wojek / A.~Suter -- \today~ ---}
+\rfoot{\includegraphics[width=2cm]{PSI-Logo_narrow.jpg}}
+
+\DeclareMathAlphabet{\bi}{OML}{cmm}{b}{it}
+
+\newcommand{\mean}[1]{\langle #1 \rangle}
+\newcommand{\lem}{LE-$\mu$SR\xspace}
+\newcommand{\lemhead}{LE-$\bm{\mu}$SR\xspace}
+\newcommand{\musr}{$\mu$SR\xspace}
+\newcommand{\musrhead}{$\bm{\mu}$SR\xspace}
+\newcommand{\trimsp}{\textsc{trim.sp}\xspace}
+\newcommand{\musrfithead}{MUSRFIT\xspace}
+\newcommand{\musrfit}{\textsc{musrfit}\xspace}
+\newcommand{\gapint}{\textsc{GapIntegrals}\xspace}
+\newcommand{\YBCO}{YBa$_{2}$Cu$_{3}$O$_{7-\delta}$\xspace}
+\newcommand{\YBCOhead}{YBa$_{\bm{2}}$Cu$_{\bm{3}}$O$_{\bm{7-\delta}}$\xspace}
+
+\newcolumntype{d}[1]{D{.}{.}{#1}}
+\newcolumntype{C}[1]{>{\centering\arraybackslash}p{#1}}
+
+\begin{document}
+% Header info --------------------------------------------------
+\thispagestyle{empty}
+\noindent
+\begin{tabular}{@{\hspace{-0.2cm}}l@{\hspace{6cm}}r}
+\noindent\includegraphics[width=3.4cm]{PSI-Logo_narrow.jpg} &
+  {\Huge\sf Memorandum}
+\end{tabular}
+%
+\vskip 1cm
+%
+\begin{tabular}{@{\hspace{-0.5cm}}ll@{\hspace{4cm}}ll}
+Date:    & September 19, 2009 / August 19, 2014       &     & \\[3ex]
+From:    & B.M.~Wojek / Modified by A. Suter & \\
+E-Mail:  & \verb?bastian.wojek@psi.ch? / \verb?andreas.suter@psi.ch? &&
+\end{tabular}
+%
+\vskip 0.3cm
+\noindent\hrulefill
+\vskip 1cm
+%
+\section*{\musrfithead plug-in for the calculation of the temperature dependence of $\bm{1/\lambda^2}$ for various gap symmetries}%
+
+This memo is intended to give a short summary of the background on which the \gapint plug-in for \musrfit \cite{musrfit} has been developped. The aim of this implementation is the efficient calculation of integrals of the form
+\begin{equation}\label{int}
+ I(T) = 1 + \frac{1}{\pi}\int_0^{2\pi}\int_{\Delta(\varphi,T)}^{\infty}\left(\frac{\partial f}{\partial E}\right) \frac{E}{\sqrt{E^2-\Delta^2(\varphi,T)}}\mathrm{d}E\mathrm{d}\varphi\,,
+\end{equation}
+where $f = (1+\exp(E/k_{\mathrm B}T))^{-1}$, like they appear e.g. in the theoretical temperature dependence of $1/\lambda^2$~\cite{Manzano}.
+In order not to do too many unnecessary function calls during the final numerical evaluation we simplify the integral (\ref{int}) as far as possible analytically. The derivative of $f$ is given by
+\begin{equation}\label{derivative}
+\frac{\partial f}{\partial E} = -\frac{1}{k_{\mathrm B}T}\frac{\exp(E/k_{\mathrm B}T)}{\left(1+\exp(E/k_{\mathrm B}T)\right)^2} = -\frac{1}{4k_{\mathrm B}T} \frac{1}{\cosh^2\left(E/2k_{\mathrm B}T\right)}.
+\end{equation}
+Using (\ref{derivative}) and doing the substitution $E'^2 = E^2-\Delta^2(\varphi,T)$, equation (\ref{int}) can be written as
+\begin{equation}
+\begin{split}
+I(T) & = 1 - \frac{1}{4\pi k_{\mathrm B}T}\int_0^{2\pi}\int_{\Delta(\varphi,T)}^{\infty}\frac{1}{\cosh^2\left(E/2k_{\mathrm B}T\right)}\frac{E}{\sqrt{E^2-\Delta^2(\varphi,T)}}\mathrm{d}E\mathrm{d}\varphi \\
+& = 1 - \frac{1}{4\pi k_{\mathrm B}T}\int_0^{2\pi}\int_{0}^{\infty}\frac{1}{\cosh^2\left(\sqrt{E'^2+\Delta^2(\varphi,T)}/2k_{\mathrm B}T\right)}\mathrm{d}E'\mathrm{d}\varphi\,.
+\end{split}
+\end{equation}
+Since a numerical integration should be performed and the function to be integrated is exponentially approaching zero for $E'\rightarrow\infty$ the infinite $E'$ integration limit can be replaced by a cutoff energy $E_{\mathrm c}$ which has to be chosen big enough:
+\begin{equation}
+I(T) \simeq \tilde{I}(T) \equiv 1 - \frac{1}{4\pi k_{\mathrm B}T}\int_0^{2\pi}\int_{0}^{E_{\mathrm c}}\frac{1}{\cosh^2\left(\sqrt{E'^2+\Delta^2(\varphi,T)}/2k_{\mathrm B}T\right)}\mathrm{d}E'\mathrm{d}\varphi\,.
+\end{equation}
+In the case that $\Delta^2(\varphi,T)$ is periodic in $\varphi$ with a period of $\pi/2$ (valid for all gap symmetries implemented at the moment), it is enough to limit the $\varphi$-integration to one period and to multiply the result by $4$:
+\begin{equation}
+\tilde{I}(T) = 1 - \frac{1}{\pi k_{\mathrm B}T}\int_0^{\pi/2}\int_{0}^{E_{\mathrm c}}\frac{1}{\cosh^2\left(\sqrt{E'^2+\Delta^2(\varphi,T)}/2k_{\mathrm B}T\right)}\mathrm{d}E'\mathrm{d}\varphi\,.
+\end{equation}
+For the numerical integration we use algorithms of the \textsc{Cuba} library \cite{cuba} which require to have a Riemann integral over the unit square. Therefore, we have to scale the integrand by the upper limits of the integrations. Note that $E_{\mathrm c}$ and $\pi/2$ (or in general the upper limit of the $\varphi$ integration) are now treated as dimensionless scaling factors.
+\begin{equation}
+\tilde{I}(T) =  1 - \frac{E_{\mathrm c}}{2k_{\mathrm B}T}\int_0^{1_{\varphi}}\int_{0}^{1_{E}}\frac{1}{\cosh^2\left(\sqrt{(E_{\mathrm c}E)^2+\Delta^2\left(\frac{\pi}{2}\varphi,T\right)}/2k_{\mathrm B}T\right)}\mathrm{d}E\mathrm{d}\varphi
+\end{equation}
+
+\subsection*{Implemented gap functions and function calls from MUSRFIT}
+At the moment the calculation of $\tilde{I}(T)$ is implemented for various gap functions all using the approximate \textsc{BCS} temperature dependence~\cite{Prozorov}
+\begin{equation}\label{eq:gapT}
+ \Delta(\varphi,T) \simeq \Delta(\varphi,0)\,\tanh\left[\frac{\pi k_{\rm B} T_{\rm c}}{\Delta_0}\sqrt{a_{\rm G} \left(\frac{T_{\rm c}}{T}-1\right)}\right]
+\end{equation}
+\noindent with $\Delta_0$ as given below, and $a_{\rm G}$ depends on the pairing state:
+
+\begin{description}
+ \item [\textit{s}-wave:] $a_{\rm G}=1$ \qquad with $\Delta_0 = 1.76\, k_{\rm B} T_c$
+ \item [\textit{d}-wave:] $a_{\rm G}=4/3$ \quad with $\Delta_0 = 2.14\, k_{\rm B} T_c$
+\end{description}
+
+\noindent Eq.(\ref{eq:gapT}) replaces the \emph{previously} used approximation~\cite{Manzano}:
+\begin{equation}
+\Delta(\varphi,T) \simeq \Delta(\varphi)\tanh\left(1.82\left(1.018\left(\frac{T_{\mathrm c}}{T}-1\right)\right)^{0.51}\right)\,.
+\end{equation}
+The \gapint plug-in calculates $\tilde{I}(T)$ for the following $\Delta(\varphi)$:
+
+\begin{description}
+ \item[\textit{s}-wave gap:]
+   \begin{equation}
+     \Delta(\varphi) = \Delta_0
+   \end{equation}
+   \musrfit theory line: \verb?userFcn libGapIntegrals TGapSWave 1 2 [3]?\\
+   (Parameters: $T_{\mathrm c}~(\mathrm{K})$, $\Delta_0~(\mathrm{meV})$, $[a_{\rm G}~(1)]$, if $a_{\rm G}$ is not given, $a_{\rm G} = 1$)
+ \item[\textit{d}-wave gap \cite{Deutscher}:] 
+   \begin{equation}
+     \Delta(\varphi) = \Delta_0\cos\left(2\varphi\right)
+   \end{equation}
+   \musrfit theory line: \verb?userFcn libGapIntegrals TGapDWave 1 2 [3]?\\
+   (Parameters: $T_{\mathrm c}~(\mathrm{K})$, $\Delta_0~(\mathrm{meV})$, $[a_{\rm G}~(1)]$, if $a_{\rm G}$ is not given, $a_{\rm G} = 4/3$)
+ \item[non-monotonic \textit{d}-wave gap \cite{Matsui}:] 
+   \begin{equation}
+     \Delta(\varphi) = \Delta_0\left[a \cos\left(2\varphi\right) + (1-a)\cos\left(6\varphi\right)\right]
+   \end{equation}
+    \musrfit theory line: \verb?userFcn libGapIntegrals TGapNonMonDWave1 1 2 3 [4]?\\
+    (Parameters: $T_{\mathrm c}~(\mathrm{K})$, $\Delta_0~(\mathrm{meV})$, $a~(1)$, $[a_{\rm G}~(1)]$, if $a_{\rm G}$ is not given, $a_{\rm G} = 4/3$)
+ \item[non-monotonic \textit{d}-wave gap \cite{Eremin}:]
+    \begin{equation}
+      \Delta(\varphi) = \Delta_0\left[\frac{2}{3} \sqrt{\frac{a}{3}}\cos\left(2\varphi\right) / \left( 1 + a\cos^2\left(2\varphi\right)\right)^{\frac{3}{2}}\right],\,a>1/2
+    \end{equation}
+    \musrfit theory line: \verb?userFcn libGapIntegrals TGapNonMonDWave2 1 2 3 [4]?\\
+    (Parameters: $T_{\mathrm c}~(\mathrm{K})$, $\Delta_0~(\mathrm{meV})$, $a~(1)$, $a~(1)$, $[a_{\rm G}~(1)]$, if $a_{\rm G}$ is not given, $a_{\rm G} = 4/3$)
+ \item[anisotropic \textit{s}-wave gap \cite{AnisotropicSWave}:] 
+    \begin{equation}
+      \Delta(\varphi) = \Delta_0\left[1+a\cos\left(4\varphi\right)\right]\,,\,0\leqslant a\leqslant1
+    \end{equation}
+    \musrfit theory line: \verb?userFcn libGapIntegrals TGapAnSWave 1 2 3 [4]?\\
+    (Parameters: $T_{\mathrm c}~(\mathrm{K})$, $\Delta_0~(\mathrm{meV})$, $a~(1)$, $[a_{\rm G}~(1)]$, if $a_{\rm G}$ is not given, $a_{\rm G} = 1$)
+\end{description}
+
+\noindent It is also possible to calculate a power law temperature dependence (in the two fluid approximation $n=4$) and the dirty \textit{s}-wave expression. 
+Obviously for this no integration is needed.
+\begin{description}
+ \item[Power law return function:] 
+   \begin{equation}
+     \frac{\lambda(0)^2}{\lambda(T)^2} = 1-\left(\frac{T}{T_{\mathrm c}}\right)^n
+   \end{equation}
+   \musrfit theory line: \verb?userFcn libGapIntegrals TGapPowerLaw 1 2?\\
+   (Parameters: $T_{\mathrm c}~(\mathrm{K})$, $n~(1)$)
+ \item[dirty \textit{s}-wave \cite{Tinkham}:]
+    \begin{equation}
+      \frac{\lambda(0)^2}{\lambda(T)^2} = \frac{\Delta(T)}{\Delta_0}\,\tanh\left[\frac{\Delta(T)}{2 k_{\rm B} T}\right]
+    \end{equation}
+    with $\Delta(T)$ given by Eq.(\ref{eq:gapT}).\\ 
+    \musrfit theory line: \verb?userFcn libGapIntegrals TGapDirtySWave 1 2 [3]?\\
+    (Parameters: $T_{\mathrm c}~(\mathrm{K})$, $\Delta_0~(\mathrm{meV})$, $[a_{\rm G}~(1)]$, if $a_{\rm G}$ is not given, $a_{\rm G} = 1$)
+\end{description}
+
+\noindent Currently there are two gap functions to be found in the code which are \emph{not} documented here: 
+\verb?TGapCosSqDWave? and \verb?TGapSinSqDWave?. For details for these gap functions (superfluid density along the \textit{a}/\textit{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 \textit{a}-\textit{b}-anisotropy)) see the source code.
+
+\subsection*{License}
+The \gapint library 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 \cite{GPL}; either version 2 of the License, or (at your option) any later version.
+
+\bibliographystyle{plain}
+\begin{thebibliography}{1}
+
+\bibitem{musrfit} A.~Suter, \textit{\musrfit User Manual}, https://wiki.intranet.psi.ch/MUSR/MusrFit
+\bibitem{cuba} T.~Hahn, \textit{Cuba -- a library for multidimensional numerical integration}, Comput.~Phys.~Commun.~\textbf{168}~(2005)~78-95, http://www.feynarts.de/cuba/
+\bibitem{Prozorov} R.~Prozorov and R.W.~Giannetta, \textit{Magnetic penetration depth in unconventional superconductors}, Supercond.\ Sci.\ Technol.\ \textbf{19}~(2006)~R41-R67, and Erratum in Supercond.\ Sci.\ Technol.\ \textbf{21}~(2008)~082003.
+\bibitem{Manzano} A.~Carrington and F.~Manzano, Physica~C~\textbf{385}~(2003)~205
+\bibitem{Deutscher} G.~Deutscher, \textit{Andreev-Saint-James reflections: A probe of cuprate superconductors}, Rev.~Mod.~Phys.~\textbf{77}~(2005)~109-135
+\bibitem{Matsui} H.~Matsui~\textit{et al.}, \textit{Direct Observation of a Nonmonotonic $d_{x^2-y^2}$-Wave Superconducting Gap in the Electron-Doped High-T$_{\mathrm c}$ Superconductor Pr$_{0.89}$LaCe$_{0.11}$CuO$_4$}, Phys.~Rev.~Lett.~\textbf{95}~(2005)~017003
+\bibitem{Eremin} I.~Eremin, E.~Tsoncheva, and A.V.~Chubukov, \textit{Signature of the nonmonotonic $d$-wave gap in electron-doped cuprates}, Phys.~Rev.~B~\textbf{77}~(2008)~024508
+\bibitem{AnisotropicSWave} ??
+\bibitem{Tinkham} M.~Tinkham, \textit{Introduction to Superconductivity} $2^{\rm nd}$ ed. (Dover Publications, New York, 2004).
+\bibitem{GPL} http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
+
+\end{thebibliography}
+
+
+\end{document}

src/external/libGapIntegrals/TGapIntegrals.cpp

@@ -2,13 +2,13 @@
 
   TGapIntegrals.cpp
 
-  Author: Bastian M. Wojek
+  Author: Bastian M. Wojek / Andreas Suter
 
 ***************************************************************************/
 
 /***************************************************************************
  *   Copyright (C) 2009 by Bastian M. Wojek                                *
- *   bastian.wojek@psi.ch                                                  *
+ *   bastian.wojek@psi.ch / andreas.suter@psi.ch                           *
  *                                                                         *
  *   This program is free software; you can redistribute it and/or modify  *
  *   it under the terms of the GNU General Public License as published by  *
@@ -64,8 +64,10 @@ ClassImp(TLambdaInvPowerLaw)
 
 ClassImp(TFilmMagnetizationDWave)
 
-// s wave  gap integral
+//--------------------------------------------------------------------
-
+/**
+ * <p> s wave  gap integral
+ */
 TGapSWave::TGapSWave() {
   TGapIntegral *gapint = new TGapIntegral();
   fGapIntegral = gapint;
@@ -78,6 +80,10 @@ TGapSWave::TGapSWave() {
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TGapDWave::TGapDWave() {
   TDWaveGapIntegralCuhre *gapint = new TDWaveGapIntegralCuhre();
   fGapIntegral = gapint;
@@ -90,6 +96,10 @@ TGapDWave::TGapDWave() {
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TGapCosSqDWave::TGapCosSqDWave() {
   TCosSqDWaveGapIntegralCuhre *gapint = new TCosSqDWaveGapIntegralCuhre();
   fGapIntegral = gapint;
@@ -102,6 +112,10 @@ TGapCosSqDWave::TGapCosSqDWave() {
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TGapSinSqDWave::TGapSinSqDWave() {
   TSinSqDWaveGapIntegralCuhre *gapint = new TSinSqDWaveGapIntegralCuhre();
   fGapIntegral = gapint;
@@ -114,6 +128,10 @@ TGapSinSqDWave::TGapSinSqDWave() {
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TGapAnSWave::TGapAnSWave() {
   TAnSWaveGapIntegralCuhre *gapint = new TAnSWaveGapIntegralCuhre();
   fGapIntegral = gapint;
@@ -126,6 +144,10 @@ TGapAnSWave::TGapAnSWave() {
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TGapNonMonDWave1::TGapNonMonDWave1() {
   TNonMonDWave1GapIntegralCuhre *gapint = new TNonMonDWave1GapIntegralCuhre();
   fGapIntegral = gapint;
@@ -138,6 +160,10 @@ TGapNonMonDWave1::TGapNonMonDWave1() {
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TGapNonMonDWave2::TGapNonMonDWave2() {
   TNonMonDWave2GapIntegralCuhre *gapint = new TNonMonDWave2GapIntegralCuhre();
   fGapIntegral = gapint;
@@ -150,46 +176,90 @@ TGapNonMonDWave2::TGapNonMonDWave2() {
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaSWave::TLambdaSWave() {
   fLambdaInvSq = new TGapSWave();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaDWave::TLambdaDWave() {
   fLambdaInvSq = new TGapDWave();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaAnSWave::TLambdaAnSWave() {
   fLambdaInvSq = new TGapAnSWave();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaNonMonDWave1::TLambdaNonMonDWave1() {
   fLambdaInvSq = new TGapNonMonDWave1();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaNonMonDWave2::TLambdaNonMonDWave2() {
   fLambdaInvSq = new TGapNonMonDWave2();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaInvSWave::TLambdaInvSWave() {
   fLambdaInvSq = new TGapSWave();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaInvDWave::TLambdaInvDWave() {
   fLambdaInvSq = new TGapDWave();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaInvAnSWave::TLambdaInvAnSWave() {
   fLambdaInvSq = new TGapAnSWave();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaInvNonMonDWave1::TLambdaInvNonMonDWave1() {
   fLambdaInvSq = new TGapNonMonDWave1();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaInvNonMonDWave2::TLambdaInvNonMonDWave2() {
   fLambdaInvSq = new TGapNonMonDWave2();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TGapSWave::~TGapSWave() {
   delete fGapIntegral;
   fGapIntegral = 0;
@@ -201,6 +271,10 @@ TGapSWave::~TGapSWave() {
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TGapDWave::~TGapDWave() {
   delete fGapIntegral;
   fGapIntegral = 0;
@@ -212,6 +286,10 @@ TGapDWave::~TGapDWave() {
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TGapCosSqDWave::~TGapCosSqDWave() {
   delete fGapIntegral;
   fGapIntegral = 0;
@@ -223,6 +301,10 @@ TGapCosSqDWave::~TGapCosSqDWave() {
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TGapSinSqDWave::~TGapSinSqDWave() {
   delete fGapIntegral;
   fGapIntegral = 0;
@@ -234,6 +316,10 @@ TGapSinSqDWave::~TGapSinSqDWave() {
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TGapAnSWave::~TGapAnSWave() {
   delete fGapIntegral;
   fGapIntegral = 0;
@@ -245,6 +331,10 @@ TGapAnSWave::~TGapAnSWave() {
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TGapNonMonDWave1::~TGapNonMonDWave1() {
   delete fGapIntegral;
   fGapIntegral = 0;
@@ -256,6 +346,10 @@ TGapNonMonDWave1::~TGapNonMonDWave1() {
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TGapNonMonDWave2::~TGapNonMonDWave2() {
   delete fGapIntegral;
   fGapIntegral = 0;
@@ -267,60 +361,110 @@ TGapNonMonDWave2::~TGapNonMonDWave2() {
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaSWave::~TLambdaSWave() {
   delete fLambdaInvSq;
   fLambdaInvSq = 0;
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaDWave::~TLambdaDWave() {
   delete fLambdaInvSq;
   fLambdaInvSq = 0;
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaAnSWave::~TLambdaAnSWave() {
   delete fLambdaInvSq;
   fLambdaInvSq = 0;
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaNonMonDWave1::~TLambdaNonMonDWave1() {
   delete fLambdaInvSq;
   fLambdaInvSq = 0;
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaNonMonDWave2::~TLambdaNonMonDWave2() {
   delete fLambdaInvSq;
   fLambdaInvSq = 0;
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaInvSWave::~TLambdaInvSWave() {
   delete fLambdaInvSq;
   fLambdaInvSq = 0;
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaInvDWave::~TLambdaInvDWave() {
   delete fLambdaInvSq;
   fLambdaInvSq = 0;
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaInvAnSWave::~TLambdaInvAnSWave() {
   delete fLambdaInvSq;
   fLambdaInvSq = 0;
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaInvNonMonDWave1::~TLambdaInvNonMonDWave1() {
   delete fLambdaInvSq;
   fLambdaInvSq = 0;
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TLambdaInvNonMonDWave2::~TLambdaInvNonMonDWave2() {
   delete fLambdaInvSq;
   fLambdaInvSq = 0;
 }
 
 
+//--------------------------------------------------------------------
+/**
+ * <p>prepare the needed parameters for the integration carried out in TGapIntegral.
+ * For details see also the Memo GapIntegrals.pdf, especially Eq.(7) and (9).
+ */
 double TGapSWave::operator()(double t, const vector<double> &par) const {
 
-  assert(par.size() == 2); // two parameters: Tc, Delta(0)
+  assert((par.size() == 2) || (par.size() == 3)); // two or three parameters: Tc (K), Delta(0) (meV), [a (1)]
+                                                  // see R. Prozorov and R. Giannetta, Supercond. Sci. Technol. 19 (2006) R41-R67
+                                                  // and Erratum Supercond. Sci. Technol. 21 (2008) 082003
+  double aa = 1.0;
+  if (par.size() == 3)
+    aa = par[2];
 
   if (t<=0.0)
     return 1.0;
@@ -364,8 +508,8 @@ double TGapSWave::operator()(double t, const vector<double> &par) const {
     double ds;
 
     vector<double> intPar; // parameters for the integral, T & Delta(T)
-    intPar.push_back(2.0*0.08617384436*t); // 2 kB T, kB in meV/K
+    intPar.push_back(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K
-    intPar.push_back(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51)));
+    intPar.push_back(par[1]*tanh(0.2707214816*par[0]/par[1]*sqrt(aa*(par[0]/t-1.0)))); // tanh(pi kB Tc / Delta(0) * sqrt()), pi kB = 0.2707214816 meV/K
 
     fGapIntegral->SetParameters(intPar);
     ds = 1.0-1.0/intPar[0]*fGapIntegral->IntegrateFunc(0.0, 2.0*(t+intPar[1]));
@@ -384,9 +528,19 @@ double TGapSWave::operator()(double t, const vector<double> &par) const {
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>prepare the needed parameters for the integration carried out in TDWaveGapIntegralCuhre.
+ * For details see also the Memo GapIntegrals.pdf, especially Eq.(7) and (10).
+ */
 double TGapDWave::operator()(double t, const vector<double> &par) const {
 
-  assert(par.size() == 2); // two parameters: Tc, Delta(0)
+  assert((par.size() == 2) || (par.size() == 3)); // two or three parameters: Tc (K), Delta(0) (meV), [a (1)]
+                                                  // see R. Prozorov and R. Giannetta, Supercond. Sci. Technol. 19 (2006) R41-R67
+                                                  // and Erratum Supercond. Sci. Technol. 21 (2008) 082003
+  double aa = 4.0/3.0;
+  if (par.size() == 3)
+    aa = par[2];
 
   if (t<=0.0)
     return 1.0;
@@ -429,9 +583,8 @@ double TGapDWave::operator()(double t, const vector<double> &par) const {
 
     double ds;
     vector<double> intPar; // parameters for the integral, T & Delta(T)
-    intPar.push_back(2.0*0.08617384436*t); // 2 kB T, kB in meV/K
+    intPar.push_back(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K
-    //intPar.push_back(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51)));
+    intPar.push_back(par[1]*tanh(0.2707214816*par[0]/par[1]*sqrt(aa*(par[0]/t-1.0)))); // tanh(pi kB Tc / Delta(0) * sqrt()), pi kB = 0.2707214816 meV/K
-    intPar.push_back(par[1]*tanh(TMath::Pi()*0.08617384436*par[0]/par[1]*sqrt(4.0/3.0*(par[0]/t-1.0))));
     intPar.push_back(4.0*(t+intPar[1])); // upper limit of energy-integration: cutoff energy
     intPar.push_back(TMath::PiOver2()); // upper limit of phi-integration
 
@@ -456,9 +609,21 @@ double TGapDWave::operator()(double t, const vector<double> &par) const {
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>prepare the needed parameters for the integration carried out in TCosSqDWaveGapIntegralCuhre.
+ * For details see also the Memo GapIntegrals.pdf, especially Eq.(7) and (??).
+ */
 double TGapCosSqDWave::operator()(double t, const vector<double> &par) const {
 
-  assert(par.size() == 3); // three parameters: Tc, DeltaD(0), DeltaS(0)
+  assert((par.size() == 3) || (par.size() == 5)); // three or five parameters: Tc (K), DeltaD(0) (meV), DeltaS(0) (meV), [aD (1), aS (1)]
+
+  double aD = 4.0/3.0;
+  double aS = 1.0;
+  if (par.size() == 5) {
+    aD = par[3];
+    aS = par[4];
+  }
 
   if (t<=0.0)
     return 1.0;
@@ -501,12 +666,11 @@ double TGapCosSqDWave::operator()(double t, const vector<double> &par) const {
 
     double ds;
     vector<double> intPar; // parameters for the integral, T & Delta(T)
-    intPar.push_back(2.0*0.08617384436*t); // 2 kB T, kB in meV/K
+    intPar.push_back(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K
-//    intPar.push_back(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51)));
+    intPar.push_back(par[1]*tanh(0.2707214816*par[0]/par[1]*sqrt(aD*(par[0]/t-1.0)))); // DeltaD(T) : tanh(pi kB Tc / Delta(0) * sqrt()), pi kB = 0.2707214816 meV/K
-    intPar.push_back(par[1]*tanh(TMath::Pi()*0.08617384436*par[0]/par[1]*sqrt(4.0/3.0*(par[0]/t-1.0)))); // DeltaD(T)
     intPar.push_back(1.0*(t+intPar[1])); // upper limit of energy-integration: cutoff energy
     intPar.push_back(TMath::Pi()); // upper limit of phi-integration
-    intPar.push_back(par[2]*tanh(TMath::Pi()*0.08617384436*par[0]/par[2]*sqrt((par[0]/t-1.0)))); // DeltaS(T)
+    intPar.push_back(par[2]*tanh(0.2707214816*par[0]/par[2]*sqrt(aS*(par[0]/t-1.0)))); // DeltaS(T) : tanh(pi kB Tc / Delta(0) * sqrt()), pi kB = 0.2707214816 meV/K
 
 //    double xl[] = {0.0, 0.0}; // lower bound E, phi
 //    double xu[] = {4.0*(t+intPar[1]), 0.5*PI}; // upper bound E, phi
@@ -529,9 +693,21 @@ double TGapCosSqDWave::operator()(double t, const vector<double> &par) const {
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>prepare the needed parameters for the integration carried out in TSinSqDWaveGapIntegralCuhre.
+ * For details see also the Memo GapIntegrals.pdf, especially Eq.(7) and (??).
+ */
 double TGapSinSqDWave::operator()(double t, const vector<double> &par) const {
 
-  assert(par.size() == 3); // two parameters: Tc, Delta(0), s-weight
+  assert((par.size() == 3) || (par.size() == 5));  // three or five parameters: Tc (K), DeltaD(0) (meV), DeltaS(0) (meV), [aD (1), aS (1)]
+
+  double aD = 4.0/3.0;
+  double aS = 1.0;
+  if (par.size() == 5) {
+    aD = par[3];
+    aS = par[4];
+  }
 
   if (t<=0.0)
     return 1.0;
@@ -574,12 +750,11 @@ double TGapSinSqDWave::operator()(double t, const vector<double> &par) const {
 
     double ds;
     vector<double> intPar; // parameters for the integral, T & Delta(T)
-    intPar.push_back(2.0*0.08617384436*t); // 2 kB T, kB in meV/K
+    intPar.push_back(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K
-//    intPar.push_back(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51)));
+    intPar.push_back(par[1]*tanh(0.2707214816*par[0]/par[1]*sqrt(aD*(par[0]/t-1.0)))); // DeltaD(T) : tanh(pi kB Tc / Delta(0) * sqrt()), pi kB = 0.2707214816 meV/K
-    intPar.push_back(par[1]*tanh(TMath::Pi()*0.08617384436*par[0]/par[1]*sqrt(4.0/3.0*(par[0]/t-1.0))));
     intPar.push_back(1.0*(t+intPar[1])); // upper limit of energy-integration: cutoff energy
     intPar.push_back(TMath::Pi()); // upper limit of phi-integration
-    intPar.push_back(par[2]*tanh(TMath::Pi()*0.08617384436*par[0]/par[2]*sqrt((par[0]/t-1.0)))); // DeltaS(T)
+    intPar.push_back(par[2]*tanh(0.2707214816*par[0]/par[2]*sqrt(aS*(par[0]/t-1.0)))); // DeltaS(T) : tanh(pi kB Tc / Delta(0) * sqrt()), pi kB = 0.2707214816 meV/K
 
 //    double xl[] = {0.0, 0.0}; // lower bound E, phi
 //    double xu[] = {4.0*(t+intPar[1]), 0.5*PI}; // upper bound E, phi
@@ -602,9 +777,19 @@ double TGapSinSqDWave::operator()(double t, const vector<double> &par) const {
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>prepare the needed parameters for the integration carried out in TAnSWaveGapIntegralCuhre (anisotropic s-wave).
+ * For details see also the Memo GapIntegrals.pdf, especially Eq.(7) and (13).
+ */
 double TGapAnSWave::operator()(double t, const vector<double> &par) const {
 
-  assert(par.size() == 3); // three parameters: Tc, Delta(0), a
+  assert((par.size() == 3) || (par.size() == 4)); // three or four parameters: Tc (K), Delta(0) (meV), a (1), [aS_Gap (1)]
+
+  double aS = 1.0;
+  if (par.size() == 4) {
+    aS = par[3];
+  }
 
   if (t<=0.0)
     return 1.0;
@@ -647,8 +832,8 @@ double TGapAnSWave::operator()(double t, const vector<double> &par) const {
 
     double ds;
     vector<double> intPar; // parameters for the integral, T & Delta(T)
-    intPar.push_back(2.0*0.08617384436*t); // 2 kB T, kB in meV/K
+    intPar.push_back(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K
-    intPar.push_back(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51)));
+    intPar.push_back(par[1]*tanh(0.2707214816*par[0]/par[1]*sqrt(aS*(par[0]/t-1.0)))); // DeltaS(T) : tanh(pi kB Tc / Delta(0) * sqrt()), pi kB = 0.2707214816 meV/K
     intPar.push_back(par[2]);
     intPar.push_back(4.0*(t+(1.0+par[2])*intPar[1])); // upper limit of energy-integration: cutoff energy
     intPar.push_back(TMath::PiOver2()); // upper limit of phi-integration
@@ -674,9 +859,19 @@ double TGapAnSWave::operator()(double t, const vector<double> &par) const {
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>prepare the needed parameters for the integration carried out in TNonMonDWave1GapIntegralCuhre.
+ * For details see also the Memo GapIntegrals.pdf, especially Eq.(7) and (11).
+ */
 double TGapNonMonDWave1::operator()(double t, const vector<double> &par) const {
 
-  assert(par.size() == 3); // two parameters: Tc, Delta(0), a
+  assert((par.size() == 3) || (par.size() == 4)); // three or four parameters: Tc (K), Delta(0) (meV), a (1), [aD_Gap (1)]
+
+  double aD = 4.0/3.0;
+  if (par.size() == 4) {
+    aD = par[3];
+  }
 
   if (t<=0.0)
     return 1.0;
@@ -719,8 +914,8 @@ double TGapNonMonDWave1::operator()(double t, const vector<double> &par) const {
 
     double ds;
     vector<double> intPar; // parameters for the integral: 2 k_B T, Delta(T), a, E_c, phi_c
-    intPar.push_back(2.0*0.08617384436*t); // 2 kB T, kB in meV/K
+    intPar.push_back(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K
-    intPar.push_back(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51)));
+    intPar.push_back(par[1]*tanh(0.2707214816*par[0]/par[1]*sqrt(aD*(par[0]/t-1.0)))); // DeltaD(T) : tanh(pi kB Tc / Delta(0) * sqrt()), pi kB = 0.2707214816 meV/K
     intPar.push_back(par[2]);
     intPar.push_back(4.0*(t+intPar[1])); // upper limit of energy-integration: cutoff energy
     intPar.push_back(TMath::PiOver2()); // upper limit of phi-integration
@@ -743,9 +938,19 @@ double TGapNonMonDWave1::operator()(double t, const vector<double> &par) const {
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>prepare the needed parameters for the integration carried out in TNonMonDWave2GapIntegralCuhre.
+ * For details see also the Memo GapIntegrals.pdf, especially Eq.(7) and (11).
+ */
 double TGapNonMonDWave2::operator()(double t, const vector<double> &par) const {
 
-  assert(par.size() == 3); // two parameters: Tc, Delta(0), a
+  assert((par.size() == 3) || (par.size() == 4)); // three parameters: Tc (K), Delta(0) (meV), a (1), [aD_Gap (1)]
+
+  double aD = 4.0/3.0;
+  if (par.size() == 4) {
+    aD = par[3];
+  }
 
   if (t<=0.0)
     return 1.0;
@@ -788,8 +993,8 @@ double TGapNonMonDWave2::operator()(double t, const vector<double> &par) const {
 
     double ds;
     vector<double> intPar; // parameters for the integral: 2 k_B T, Delta(T), a, E_c, phi_c
-    intPar.push_back(2.0*0.08617384436*t); // 2 kB T, kB in meV/K
+    intPar.push_back(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K
-    intPar.push_back(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51)));
+    intPar.push_back(par[1]*tanh(0.2707214816*par[0]/par[1]*sqrt(aD*(par[0]/t-1.0)))); // DeltaD(T) : tanh(pi kB Tc / Delta(0) * sqrt()), pi kB = 0.2707214816 meV/K
     intPar.push_back(par[2]);
     intPar.push_back(4.0*(t+intPar[1])); // upper limit of energy-integration: cutoff energy
     intPar.push_back(TMath::PiOver2()); // upper limit of phi-integration
@@ -812,9 +1017,14 @@ double TGapNonMonDWave2::operator()(double t, const vector<double> &par) const {
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>Superfluid density in the ``two-fluid'' approximation.
+ * For details see also the Memo GapIntegrals.pdf, especially Eq.(7) and (14).
+ */
 double TGapPowerLaw::operator()(double t, const vector<double> &par) const {
 
-  assert(par.size() == 2); // two parameters: Tc, N
+  assert(par.size() == 2); // two parameters: Tc, n
 
   if(t<=0.0)
     return 1.0;
@@ -826,22 +1036,35 @@ double TGapPowerLaw::operator()(double t, const vector<double> &par) const {
 }
 
 
+//--------------------------------------------------------------------
+/**
+ * <p>Superfluid density for a dirty s-wave superconductor.
+ * For details see also the Memo GapIntegrals.pdf, especially Eq.(7) and (15).
+ */
 double TGapDirtySWave::operator()(double t, const vector<double> &par) const {
 
-  assert(par.size() == 2); // two parameters: Tc, Delta(0)
+  assert((par.size() == 2) || (par.size() == 3)); // two or three parameters: Tc (K), Delta(0) (meV), [a (1)]
+
+  double a = 1.0;
+  if (par.size() == 3) {
+    a = par[2];
+  }
 
   if (t<=0.0)
     return 1.0;
   else if (t >= par[0])
     return 0.0;
 
-  double deltaT(tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51)));
+  double deltaT(tanh(0.2707214816*par[0]/par[1]*sqrt(a*(par[0]/t-1.0)))); // DeltaS(T) : tanh(pi kB Tc / Delta(0) * sqrt()), pi kB = 0.2707214816 meV/K
-
-  return deltaT*tanh(par[1]*deltaT/(2.0*0.08617384436*t)); // Delta(T)/Delta(0)*tanh(Delta(T)/2 kB T), kB in meV/K
 
+  return deltaT*tanh(par[1]*deltaT/(0.172346648*t)); // Delta(T)/Delta(0)*tanh(Delta(T)/2 kB T), kB in meV/K
 }
 
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 double TLambdaSWave::operator()(double t, const vector<double> &par) const
 {
   assert(par.size() == 2); // two parameters: Tc, Delta0
@@ -856,6 +1079,10 @@ double TLambdaSWave::operator()(double t, const vector<double> &par) const
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 double TLambdaDWave::operator()(double t, const vector<double> &par) const
 {
   assert(par.size() == 2); // two parameters: Tc, Delta0
@@ -870,6 +1097,10 @@ double TLambdaDWave::operator()(double t, const vector<double> &par) const
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 double TLambdaAnSWave::operator()(double t, const vector<double> &par) const
 {
   assert(par.size() == 3); // three parameters: Tc, Delta0, a
@@ -884,6 +1115,10 @@ double TLambdaAnSWave::operator()(double t, const vector<double> &par) const
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 double TLambdaNonMonDWave1::operator()(double t, const vector<double> &par) const
 {
   assert(par.size() == 3); // three parameters: Tc, Delta0, a
@@ -898,6 +1133,10 @@ double TLambdaNonMonDWave1::operator()(double t, const vector<double> &par) cons
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 double TLambdaNonMonDWave2::operator()(double t, const vector<double> &par) const
 {
   assert(par.size() == 3); // three parameters: Tc, Delta0, a
@@ -912,6 +1151,10 @@ double TLambdaNonMonDWave2::operator()(double t, const vector<double> &par) cons
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 double TLambdaPowerLaw::operator()(double t, const vector<double> &par) const {
 
   assert(par.size() == 2); // two parameters: Tc, N
@@ -926,6 +1169,10 @@ double TLambdaPowerLaw::operator()(double t, const vector<double> &par) const {
 }
 
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 double TLambdaInvSWave::operator()(double t, const vector<double> &par) const
 {
   assert(par.size() == 2); // two parameters: Tc, Delta0
@@ -940,6 +1187,10 @@ double TLambdaInvSWave::operator()(double t, const vector<double> &par) const
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 double TLambdaInvDWave::operator()(double t, const vector<double> &par) const
 {
   assert(par.size() == 2); // two parameters: Tc, Delta0
@@ -954,6 +1205,10 @@ double TLambdaInvDWave::operator()(double t, const vector<double> &par) const
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 double TLambdaInvAnSWave::operator()(double t, const vector<double> &par) const
 {
   assert(par.size() == 3); // three parameters: Tc, Delta0, a
@@ -968,6 +1223,10 @@ double TLambdaInvAnSWave::operator()(double t, const vector<double> &par) const
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 double TLambdaInvNonMonDWave1::operator()(double t, const vector<double> &par) const
 {
   assert(par.size() == 3); // three parameters: Tc, Delta0, a
@@ -982,6 +1241,10 @@ double TLambdaInvNonMonDWave1::operator()(double t, const vector<double> &par) c
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 double TLambdaInvNonMonDWave2::operator()(double t, const vector<double> &par) const
 {
   assert(par.size() == 3); // three parameters: Tc, Delta0, a
@@ -996,6 +1259,10 @@ double TLambdaInvNonMonDWave2::operator()(double t, const vector<double> &par) c
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 double TLambdaInvPowerLaw::operator()(double t, const vector<double> &par) const {
 
   assert(par.size() == 2); // two parameters: Tc, N
@@ -1009,12 +1276,20 @@ double TLambdaInvPowerLaw::operator()(double t, const vector<double> &par) const
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TFilmMagnetizationDWave::TFilmMagnetizationDWave()
 {
   fLambdaInvSq = new TGapDWave();
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 TFilmMagnetizationDWave::~TFilmMagnetizationDWave()
 {
   delete fLambdaInvSq;
@@ -1022,6 +1297,10 @@ TFilmMagnetizationDWave::~TFilmMagnetizationDWave()
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 double TFilmMagnetizationDWave::operator()(double t, const vector<double> &par) const
 {
   assert(par.size() == 4); // four parameters: Tc, Delta0, lambda0, film-thickness

src/external/libGapIntegrals/TGapIntegrals.h

@@ -37,6 +37,10 @@ using namespace std;
 #include "PUserFcnBase.h"
 #include "BMWIntegrator.h"
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TGapSWave : public PUserFcnBase {
 
 public:
@@ -61,6 +65,10 @@ private:
   ClassDef(TGapSWave,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TGapDWave : public PUserFcnBase {
 
 public:
@@ -85,6 +93,10 @@ private:
   ClassDef(TGapDWave,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TGapCosSqDWave : public PUserFcnBase {
 
 public:
@@ -109,6 +121,10 @@ private:
   ClassDef(TGapCosSqDWave,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TGapSinSqDWave : public PUserFcnBase {
 
 public:
@@ -133,7 +149,10 @@ private:
   ClassDef(TGapSinSqDWave,1)
 };
 
-
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TGapAnSWave : public PUserFcnBase {
 
 public:
@@ -158,6 +177,10 @@ private:
   ClassDef(TGapAnSWave,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TGapNonMonDWave1 : public PUserFcnBase {
 
 public:
@@ -182,6 +205,10 @@ private:
   ClassDef(TGapNonMonDWave1,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TGapNonMonDWave2 : public PUserFcnBase {
 
 public:
@@ -207,6 +234,10 @@ private:
 };
 
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TGapPowerLaw : public PUserFcnBase {
 
 public:
@@ -224,6 +255,10 @@ private:
   ClassDef(TGapPowerLaw,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TGapDirtySWave : public PUserFcnBase {
 
 public:
@@ -242,6 +277,10 @@ private:
 };
 
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TLambdaSWave : public PUserFcnBase {
 
 public:
@@ -260,6 +299,10 @@ private:
   ClassDef(TLambdaSWave,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TLambdaDWave : public PUserFcnBase {
 
 public:
@@ -278,6 +321,10 @@ private:
   ClassDef(TLambdaDWave,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TLambdaAnSWave : public PUserFcnBase {
 
 public:
@@ -296,6 +343,10 @@ private:
   ClassDef(TLambdaAnSWave,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TLambdaNonMonDWave1 : public PUserFcnBase {
 
 public:
@@ -314,6 +365,10 @@ private:
   ClassDef(TLambdaNonMonDWave1,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TLambdaNonMonDWave2 : public PUserFcnBase {
 
 public:
@@ -332,7 +387,10 @@ private:
   ClassDef(TLambdaNonMonDWave2,1)
 };
 
-
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TLambdaPowerLaw : public PUserFcnBase {
 
 public:
@@ -350,6 +408,10 @@ private:
   ClassDef(TLambdaPowerLaw,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TLambdaInvSWave : public PUserFcnBase {
 
 public:
@@ -368,6 +430,10 @@ private:
   ClassDef(TLambdaInvSWave,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TLambdaInvDWave : public PUserFcnBase {
 
 public:
@@ -386,6 +452,10 @@ private:
   ClassDef(TLambdaInvDWave,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TLambdaInvAnSWave : public PUserFcnBase {
 
 public:
@@ -404,6 +474,10 @@ private:
   ClassDef(TLambdaInvAnSWave,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TLambdaInvNonMonDWave1 : public PUserFcnBase {
 
 public:
@@ -422,6 +496,10 @@ private:
   ClassDef(TLambdaInvNonMonDWave1,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TLambdaInvNonMonDWave2 : public PUserFcnBase {
 
 public:
@@ -440,7 +518,10 @@ private:
   ClassDef(TLambdaInvNonMonDWave2,1)
 };
 
-
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TLambdaInvPowerLaw : public PUserFcnBase {
 
 public:
@@ -458,6 +539,10 @@ private:
   ClassDef(TLambdaInvPowerLaw,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
 class TFilmMagnetizationDWave : public PUserFcnBase {
 
 public: