add proper p-wave (line,point) superfluid density calculation -- adopted to DKS.

src/external/BMWtools/BMWIntegrator.cpp

@@ -35,6 +35,182 @@
 #define SEED 0
 #define STATEFILE NULL
 
+//-----------------------------------------------------------------------------
+
+std::vector<double> TPointPWaveGapIntegralCuhre::fPar;
+
+//-----------------------------------------------------------------------------
+/**
+ * <p>Integrate the function using the Cuhre interface
+ *
+ * <p><b>return:</b>
+ * - value of the integral
+ */
+double TPointPWaveGapIntegralCuhre::IntegrateFunc(int tag)
+{
+  const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
+  const double EPSREL (1e-4);
+  const double EPSABS (1e-6);
+  const unsigned int VERBOSE (0);
+  const unsigned int LAST (4);
+  const unsigned int MINEVAL (0);
+  const unsigned int MAXEVAL (50000);
+
+  const unsigned int KEY (13);
+
+  int nregions, neval, fail;
+  double integral[NCOMP], error[NCOMP], prob[NCOMP];
+
+  if (tag == 0)
+    Cuhre(fNDim, NCOMP, Integrand_aa, USERDATA, NVEC,
+      EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
+      KEY, STATEFILE, SPIN,
+      &nregions, &neval, &fail, integral, error, prob);
+  else
+    Cuhre(fNDim, NCOMP, Integrand_cc, USERDATA, NVEC,
+      EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
+      KEY, STATEFILE, SPIN,
+      &nregions, &neval, &fail, integral, error, prob);
+
+  return integral[0];
+}
+
+//-----------------------------------------------------------------------------
+/**
+ * <p>Calculate the function value for the use with Cuhre---actual implementation of the function
+ * for p-wave point, aa==bb component
+ *
+ * <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 TPointPWaveGapIntegralCuhre::Integrand_aa(const int *ndim, const double x[],
+                      const int *ncomp, double f[], void *userdata) // x = {E, z}, fPar = {twokBT, Delta(T), Ec, zc}
+{
+  double z = x[1]*fPar[3];
+  double deltasq(pow(sqrt(1.0-z*z)*fPar[1],2.0));
+  f[0] = (1.0-z*z)/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0);
+  return 0;
+}
+
+//-----------------------------------------------------------------------------
+/**
+ * <p>Calculate the function value for the use with Cuhre---actual implementation of the function
+ * for p-wave point, cc component
+ *
+ * <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 TPointPWaveGapIntegralCuhre::Integrand_cc(const int *ndim, const double x[],
+                      const int *ncomp, double f[], void *userdata) // x = {E, z}, fPar = {twokBT, Delta(T), Ec, zc}
+{
+  double z = x[1]*fPar[3];
+  double deltasq(pow(sqrt(1.0-z*z)*fPar[1],2.0));
+  f[0] = (z*z)/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0);
+  return 0;
+}
+
+//-----------------------------------------------------------------------------
+
+std::vector<double> TLinePWaveGapIntegralCuhre::fPar;
+
+//-----------------------------------------------------------------------------
+/**
+ * <p>Integrate the function using the Cuhre interface
+ *
+ * <p><b>return:</b>
+ * - value of the integral
+ */
+double TLinePWaveGapIntegralCuhre::IntegrateFunc(int tag)
+{
+  const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
+  const double EPSREL (1e-4);
+  const double EPSABS (1e-6);
+  const unsigned int VERBOSE (0);
+  const unsigned int LAST (4);
+  const unsigned int MINEVAL (0);
+  const unsigned int MAXEVAL (50000);
+
+  const unsigned int KEY (13);
+
+  int nregions, neval, fail;
+  double integral[NCOMP], error[NCOMP], prob[NCOMP];
+
+  if (tag == 0)
+    Cuhre(fNDim, NCOMP, Integrand_aa, USERDATA, NVEC,
+      EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
+      KEY, STATEFILE, SPIN,
+      &nregions, &neval, &fail, integral, error, prob);
+  else
+    Cuhre(fNDim, NCOMP, Integrand_cc, USERDATA, NVEC,
+      EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
+      KEY, STATEFILE, SPIN,
+      &nregions, &neval, &fail, integral, error, prob);
+
+  return integral[0];
+}
+
+//-----------------------------------------------------------------------------
+/**
+ * <p>Calculate the function value for the use with Cuhre---actual implementation of the function
+ * for p-wave line, aa==bb component
+ *
+ * <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 TLinePWaveGapIntegralCuhre::Integrand_aa(const int *ndim, const double x[],
+                      const int *ncomp, double f[], void *userdata) // x = {E, z}, fPar = {twokBT, Delta(T), Ec, zc}
+{
+  double z = x[1]*fPar[3];
+  double deltasq(pow(z*fPar[1],2.0));
+  f[0] = (1.0-z*z)/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0);
+  return 0;
+}
+
+//-----------------------------------------------------------------------------
+/**
+ * <p>Calculate the function value for the use with Cuhre---actual implementation of the function
+ * for p-wave line, cc component
+ *
+ * <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 TLinePWaveGapIntegralCuhre::Integrand_cc(const int *ndim, const double x[],
+                      const int *ncomp, double f[], void *userdata) // x = {E, z}, fPar = {twokBT, Delta(T), Ec, zc}
+{
+  double z = x[1]*fPar[3];
+  double deltasq(pow(z*fPar[1],2.0));
+  f[0] = (z*z)/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0);
+  return 0;
+}
+
+//-----------------------------------------------------------------------------
+
 std::vector<double> TDWaveGapIntegralCuhre::fPar;
 
 /**

src/external/BMWtools/BMWIntegrator.h

@@ -249,6 +249,47 @@ inline double TMCIntegrator::IntegrateFunc(size_t dim, double *x1, double *x2)
   return fMCIntegrator->Integral(fFunc, dim, x1, x2, (this));
 }
 
+//-----------------------------------------------------------------------------
+/**
+ * <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 point p symmetry of the superconducting order parameter.
+ *    The integration uses the Cuhre algorithm of the Cuba library.
+ */
+class TPointPWaveGapIntegralCuhre {
+  public:
+    TPointPWaveGapIntegralCuhre() : fNDim(2) {}
+    ~TPointPWaveGapIntegralCuhre() { fPar.clear(); }
+    void SetParameters(const std::vector<double> &par) { fPar=par; }
+    static int Integrand_aa(const int*, const double[], const int*, double[], void*);
+    static int Integrand_cc(const int*, const double[], const int*, double[], void*);
+    double IntegrateFunc(int tag);
+
+  protected:
+    static std::vector<double> fPar; ///< parameters of the integrand
+    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 a line p symmetry of the superconducting order parameter.
+ *    The integration uses the Cuhre algorithm of the Cuba library.
+ */
+class TLinePWaveGapIntegralCuhre {
+  public:
+    TLinePWaveGapIntegralCuhre() : fNDim(2) {}
+    ~TLinePWaveGapIntegralCuhre() { fPar.clear(); }
+    void SetParameters(const std::vector<double> &par) { fPar=par; }
+    static int Integrand_aa(const int*, const double[], const int*, double[], void*);
+    static int Integrand_cc(const int*, const double[], const int*, double[], void*);
+    double IntegrateFunc(int tag);
+
+  protected:
+    static std::vector<double> fPar; ///< parameters of the integrand
+    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 a d_{x^2-y^2} symmetry of the superconducting order parameter.

src/external/libGapIntegrals/GapIntegrals.tex

@@ -74,16 +74,17 @@ E-Mail:  & \verb?andreas.suter@psi.ch? &&
 \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 developed. The aim of this implementation is the efficient calculation of integrals of the form
-\begin{equation}\label{int}
+\begin{equation}\label{eq:int_phi}
  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
+For gap symmetries which involve not only a $E$- and $\varphi$-dependence but also a $\theta$-dependence, see the special section towards the end of the memo.
+In order not to do too many unnecessary function calls during the final numerical evaluation we simplify the integral (\ref{eq:int_phi}) 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
+Using (\ref{derivative}) and doing the substitution $E'^2 = E^2-\Delta^2(\varphi,T)$, equation (\ref{eq:int_phi}) can be written as
-\begin{equation}
+\begin{equation}\label{eq:bmw_2d}
 \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\,.
@@ -167,6 +168,22 @@ The \gapint plug-in calculates $\tilde{I}(T)$ for the following $\Delta(\varphi)
     Parameters: $T_{\mathrm c}~(\mathrm{K})$, $\Delta_0~(\mathrm{meV})$, $a~(1)$, $[c_0~(1),~a_{\rm G}~(1)]$. 
     If $c_0$ and $a_{\rm G}$ are provided, the temperature dependence according to Eq.(\ref{eq:gapT_Prozorov}) will be used, 
     otherwise Eq.(\ref{eq:gapT_Manzano}) will be utilized.
+\item[\textit{p}-wave (point) \cite{Pang2015}:]
+   \begin{equation}
+     \Delta(\theta, T) = \Delta(T) \sin(\theta) = \Delta(T) \cdot \sqrt{1-z^2}
+   \end{equation}
+   \musrfit theory line: \verb?userFcn libGapIntegrals TGapPointPWave 1 2 [3 [4 5]]?\\[1.5ex]
+    Parameters: $T_{\mathrm c}~(\mathrm{K})$, $\Delta_0~(\mathrm{meV})$, [ \verb?orientation_tag?, $[c_0~(1),~a_{\rm G}~(1)]$]. If $c_0$ and $a_{\rm G}$ are provided, 
+    the temperature dependence according to Eq.(\ref{eq:gapT_Prozorov}) will be used, otherwise Eq.(\ref{eq:gapT_Manzano}) will be utilized. \\
+    \verb?orientation_tag?: $0=\{aa,bb\}$, $1=cc$, and the default $2=$ average (see Eq.\ (\ref{eq:n_avg}))
+\item[\textit{p}-wave (line) \cite{Ozaki1986}:]
+   \begin{equation}
+     \Delta(\theta, T) = \Delta(T) \cos(\theta) = \Delta(T) \cdot z
+   \end{equation}
+   \musrfit theory line: \verb?userFcn libGapIntegrals TGapLinePWave 1 2 [3 [4 5]]?\\[1.5ex]
+    Parameters: $T_{\mathrm c}~(\mathrm{K})$, $\Delta_0~(\mathrm{meV})$, [ \verb?orientation_tag?, $[c_0~(1),~a_{\rm G}~(1)]$]. If $c_0$ and $a_{\rm G}$ are provided, 
+    the temperature dependence according to Eq.(\ref{eq:gapT_Prozorov}) will be used, otherwise Eq.(\ref{eq:gapT_Manzano}) will be utilized. \\
+    \verb?orientation_tag?: $0=\{aa,bb\}$, $1=cc$, and the default $2=$ average (see Eq.\ (\ref{eq:n_avg}))
 \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. 
@@ -193,6 +210,69 @@ Obviously for this no integration is needed.
 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*{Gap Integrals for $\bm{\theta}$-, and $\bm{(\theta, \varphi)}$-dependent Gaps}%
+
+First some general formulae as found in Ref.\,\cite{Prozorov}. It assumes an anisotropic response which can be classified in 3 directions ($a$, $b$, and $c$).
+
+\noindent For the case of a 2D Fermi surface (cylindrical symmetry):
+
+\begin{equation}\label{eq:n_anisotrope_2D}
+n_{aa \atop bb}(T) = 1 - \frac{1}{2\pi k_{\rm B} T} \int_0^{2\pi} \mathrm{d}\varphi\, {\cos^2(\varphi) \atop \sin^2(\varphi)} \underbrace{\int_0^\infty \mathrm{d}\varepsilon\, \left\{ \cosh\left[\frac{\sqrt{\varepsilon^2 + \Delta^2}}{2 k_{\rm B}T}\right]\right\}^{-2}}_{= G(\Delta(\varphi), T)}
+\end{equation}
+
+\noindent For the case of a 3D Fermi surface:
+
+\begin{eqnarray}
+ n_{aa \atop bb}(T) &=& 1 - \frac{3}{4\pi k_{\rm B} T} \int_0^1 \mathrm{d}z\, (1-z^2) \int_0^{2\pi} \mathrm{d}\varphi\, {\cos^2(\varphi) \atop \sin^2(\varphi)} \cdot G(\Delta(z,\varphi), T) \label{eq:n_anisotrope_3D_aabb} \\
+ n_{cc}(T) &=& 1 - \frac{3}{2\pi k_{\rm B} T} \int_0^1 \mathrm{d}z\, z^2 \int_0^{2\pi} \mathrm{d}\varphi\, \cos^2(\varphi) \cdot G(\Delta(z,\varphi), T) \label{eq:n_anisotrope_3D_cc}
+\end{eqnarray}
+
+\noindent The ``powder averaged'' superfluid density is then defined as
+
+\begin{equation}\label{eq:n_avg}
+ n_{\rm S} = \frac{1}{3}\cdot \left[ \sqrt{n_{aa} n_{bb}} + \sqrt{n_{aa} n_{cc}} + \sqrt{n_{bb} n_{cc}} \right]
+\end{equation}
+
+\subsubsection*{Isotropic s-Wave Gap}
+
+\noindent For the 2D/3D case this means that $\Delta$ is just a constant. 
+
+\noindent For the 2D case it follows
+
+\begin{equation}
+  n_{aa \atop bb}(T) = 1 - \frac{1}{2 k_{\rm B} T} \cdot G(\Delta, T) = n_{\rm S}(T).
+\end{equation}
+
+\noindent This is the same as Eq.(\ref{eq:bmw_2d}), assuming a $\Delta \neq f(\varphi)$.
+
+\vspace{5mm}
+
+\noindent The 3D case for $\Delta \neq f(\theta, \varphi)$:
+
+\noindent The variable transformation $z = \cos(\theta)$ leads to $\mathrm{d}z = -\sin(\theta)\,\mathrm{d}\theta$, $z=0 \to \theta=\pi/2$, $z=1 \to \theta=0$, and hence to
+
+\begin{eqnarray*}
+ n_{aa \atop bb}(T) &=& 1 - \frac{3}{4\pi k_{\rm B} T} \underbrace{\int_0^{\pi/2} \mathrm{d}\theta \, \sin^3(\theta)}_{= 2/3} \, \underbrace{\int_0^{2\pi} \mathrm{d}\varphi {\cos^2(\varphi) \atop \sin^2(\varphi)}}_{=\pi} \cdot G(\Delta, T) \\
+ &=& 1 - \frac{1}{2 k_{\rm B} T} \cdot G(\Delta, T). \\
+ n_{cc}(T) &=& 1 - \frac{3}{2\pi k_{\rm B} T} \underbrace{\int_0^{\pi/2} \mathrm{d}\theta \, \cos^2(\theta)\sin(\theta)}_{=1/3} \, \underbrace{\int_0^{2\pi} \mathrm{d}\varphi \cos^2(\varphi)}_{=\pi} \cdot G(\Delta, T) \\
+ &=&  1 - \frac{1}{2 k_{\rm B} T} \cdot G(\Delta, T). 
+\end{eqnarray*}
+
+\noindent And hence
+
+$$ n_{\rm S}(T) = 1- \frac{1}{2 k_{\rm B} T} \cdot G(\Delta, T). $$
+
+\subsubsection*{3D Fermi Surface Gap $\mathbf{\Delta \neq f(\bm\varphi)}$}
+
+For this case the superfluid density integrals reduce to ($z=\cos(\theta)$)
+
+\begin{eqnarray}
+  n_{aa \atop bb}(T) &=& 1 - \frac{3}{4 k_{\rm B} T} \int_0^1 \mathrm{d}z\, (1-z^2) \cdot G(\Delta(z, T),T) \\
+  n_{cc}(T) &=& 1 - \frac{3}{2 k_{\rm B} T} \int_0^1 \mathrm{d}z\, z^2 \cdot G(\Delta(z, T),T)
+\end{eqnarray}
+
+
+
 \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.
 
@@ -208,6 +288,8 @@ The \gapint library is free software; you can redistribute it and/or modify it u
 \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{Pang2015} G.M.~Pang, \emph{et al.}, Phys.~Rev.~B~\textbf{91}~(2015)~220502(R), and references in there.
+\bibitem{Ozaki1986} M.~Ozaki, \emph{et al.}, Prog.~Theor.~Phys.~\textbf{75}~(1986)~442.
 \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
 

src/external/libGapIntegrals/TGapIntegrals.cpp

@@ -36,6 +36,8 @@
 #define TWOPI 6.28318530717958647692
 
 ClassImp(TGapSWave)
+ClassImp(TGapPointPWave)
+ClassImp(TGapLinePWave)
 ClassImp(TGapDWave)
 ClassImp(TGapCosSqDWave)
 ClassImp(TGapSinSqDWave)
@@ -46,6 +48,8 @@ ClassImp(TGapPowerLaw)
 ClassImp(TGapDirtySWave)
 
 ClassImp(TLambdaSWave)
+ClassImp(TLambdaPointPWave)
+ClassImp(TLambdaLinePWave)
 ClassImp(TLambdaDWave)
 ClassImp(TLambdaAnSWave)
 ClassImp(TLambdaNonMonDWave1)
@@ -53,6 +57,8 @@ ClassImp(TLambdaNonMonDWave2)
 ClassImp(TLambdaPowerLaw)
 
 ClassImp(TLambdaInvSWave)
+ClassImp(TLambdaInvPointPWave)
+ClassImp(TLambdaInvLinePWave)
 ClassImp(TLambdaInvDWave)
 ClassImp(TLambdaInvAnSWave)
 ClassImp(TLambdaInvNonMonDWave1)
@@ -77,6 +83,38 @@ TGapSWave::TGapSWave() {
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p> point p wave  gap integral
+ */
+TGapPointPWave::TGapPointPWave() {
+  TPointPWaveGapIntegralCuhre *gapint = new TPointPWaveGapIntegralCuhre();
+  fGapIntegral = gapint;
+  gapint = nullptr;
+
+  fTemp.clear();
+  fTempIter = fTemp.end();
+  fIntegralValues.clear();
+  fCalcNeeded.clear();
+  fPar.clear();
+}
+
+//--------------------------------------------------------------------
+/**
+ * <p> line p wave  gap integral
+ */
+TGapLinePWave::TGapLinePWave() {
+  TLinePWaveGapIntegralCuhre *gapint = new TLinePWaveGapIntegralCuhre();
+  fGapIntegral = gapint;
+  gapint = nullptr;
+
+  fTemp.clear();
+  fTempIter = fTemp.end();
+  fIntegralValues.clear();
+  fCalcNeeded.clear();
+  fPar.clear();
+}
+
 //--------------------------------------------------------------------
 /**
  * <p>
@@ -181,6 +219,22 @@ TLambdaSWave::TLambdaSWave() {
   fLambdaInvSq = new TGapSWave();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+TLambdaPointPWave::TLambdaPointPWave() {
+  fLambdaInvSq = new TGapPointPWave();
+}
+
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+TLambdaLinePWave::TLambdaLinePWave() {
+  fLambdaInvSq = new TGapLinePWave();
+}
+
 //--------------------------------------------------------------------
 /**
  * <p>
@@ -221,6 +275,22 @@ TLambdaInvSWave::TLambdaInvSWave() {
   fLambdaInvSq = new TGapSWave();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+TLambdaInvPointPWave::TLambdaInvPointPWave() {
+  fLambdaInvSq = new TGapPointPWave();
+}
+
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+TLambdaInvLinePWave::TLambdaInvLinePWave() {
+  fLambdaInvSq = new TGapLinePWave();
+}
+
 //--------------------------------------------------------------------
 /**
  * <p>
@@ -268,6 +338,36 @@ TGapSWave::~TGapSWave() {
   fPar.clear();
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+TGapPointPWave::~TGapPointPWave() {
+  delete fGapIntegral;
+  fGapIntegral = nullptr;
+
+  fTemp.clear();
+  fTempIter = fTemp.end();
+  fIntegralValues.clear();
+  fCalcNeeded.clear();
+  fPar.clear();
+}
+
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+TGapLinePWave::~TGapLinePWave() {
+  delete fGapIntegral;
+  fGapIntegral = nullptr;
+
+  fTemp.clear();
+  fTempIter = fTemp.end();
+  fIntegralValues.clear();
+  fCalcNeeded.clear();
+  fPar.clear();
+}
+
 //--------------------------------------------------------------------
 /**
  * <p>
@@ -367,6 +467,24 @@ TLambdaSWave::~TLambdaSWave() {
   fLambdaInvSq = nullptr;
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+TLambdaPointPWave::~TLambdaPointPWave() {
+  delete fLambdaInvSq;
+  fLambdaInvSq = nullptr;
+}
+
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+TLambdaLinePWave::~TLambdaLinePWave() {
+  delete fLambdaInvSq;
+  fLambdaInvSq = nullptr;
+}
+
 //--------------------------------------------------------------------
 /**
  * <p>
@@ -412,6 +530,24 @@ TLambdaInvSWave::~TLambdaInvSWave() {
   fLambdaInvSq = nullptr;
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+TLambdaInvPointPWave::~TLambdaInvPointPWave() {
+  delete fLambdaInvSq;
+  fLambdaInvSq = nullptr;
+}
+
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+TLambdaInvLinePWave::~TLambdaInvLinePWave() {
+  delete fLambdaInvSq;
+  fLambdaInvSq = nullptr;
+}
+
 //--------------------------------------------------------------------
 /**
  * <p>
@@ -528,6 +664,191 @@ double TGapSWave::operator()(double t, const std::vector<double> &par) const {
 
 }
 
+//--------------------------------------------------------------------
+/**
+ * <p>prepare the needed parameters for the integration carried out in TPointPWaveGapIntegralCuhre.
+ * For details see also the Memo GapIntegrals.pdf, , especially Eq.(19) and (20).
+ */
+double TGapPointPWave::operator()(double t, const std::vector<double> &par) const {
+
+  // parameters: [0] Tc (K), [1] Delta0 (meV), [[2] orientation tag, [[3] c0 (1), [4] aG (1)]]
+
+  assert((par.size() >= 2) && (par.size() <= 5)); // 2 parameters: see A.~Carrington and F.~Manzano, Physica~C~\textbf{385}~(2003)~205
+                                                  // 4 or 5 parameters: see R. Prozorov and R. Giannetta, Supercond. Sci. Technol. 19 (2006) R41-R67
+                                                  // and Erratum Supercond. Sci. Technol. 21 (2008) 082003
+                                                  // c0 in the original context is c0 = (pi kB Tc) / Delta0
+                                                  // orientation tag: 0=aa,bb; 1=cc; 2=(sqrt[aa bb] + sqrt[aa cc] + sqrt[bb cc])/3 (default)
+  if (t <= 0.0)
+    return 1.0;
+  else if (t >= par[0])
+    return 0.0;
+
+  // check if orientation tag is given
+  int orientation_tag(2);
+  if ((par.size()==3) || (par.size()==5))
+    orientation_tag = static_cast<int>(par[2]);
+
+  bool integralParChanged(false);
+
+  if (fPar.empty()) { // first time calling this routine
+    fPar = par;
+    integralParChanged = true;
+  } else { // check if Tc or Delta0 have changed
+    for (unsigned int i(0); i<par.size(); i++) {
+      if (par[i] != fPar[i]) {
+        fPar[i] = par[i];
+        integralParChanged = true;
+      }
+    }
+  }
+
+  bool newTemp(false);
+  unsigned int vectorIndex;
+
+  if (integralParChanged) {
+    fCalcNeeded.clear();
+    fCalcNeeded.resize(fTemp.size(), true);
+  }
+
+  fTempIter = find(fTemp.begin(), fTemp.end(), t);
+  if(fTempIter == fTemp.end()) {
+    fTemp.push_back(t);
+    vectorIndex = fTemp.size() - 1;
+    fCalcNeeded.push_back(true);
+    newTemp = true;
+  } else {
+    vectorIndex = fTempIter - fTemp.begin();
+  }
+
+  if (fCalcNeeded[vectorIndex]) {
+    double ds, ds1;
+    std::vector<double> intPar; // parameters for the integral, T & Delta(T)
+    intPar.push_back(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K
+    if ((par.size() == 2) || (par.size() == 3)) { // Carrington/Manzano
+      intPar.push_back(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51)));
+    } else { // Prozorov/Giannetta
+      intPar.push_back(par[1]*tanh(par[2]*sqrt(par[3]*(par[0]/t-1.0)))); // Delta0*tanh(c0*sqrt(aG*(Tc/T-1)))
+    }
+    intPar.push_back(4.0*(t+intPar[1])); // upper limit of energy-integration: cutoff energy
+    intPar.push_back(1.0); // upper limit of theta-integration
+
+    fGapIntegral->SetParameters(intPar);
+    if (orientation_tag == 0) // aa,bb
+      ds = 1.0-(intPar[2]*3.0)/(2.0*intPar[0])*fGapIntegral->IntegrateFunc(0); // integral prefactor is by 2 lower [Eqs.(19,20)] since intPar[0]==2kB T!
+    else if (orientation_tag == 1) // cc
+      ds = 1.0-(intPar[2]*3.0)/(intPar[0])*fGapIntegral->IntegrateFunc(1); // integral prefactor is by 2 lower [Eqs.(19,20)] since intPar[0]==2kB T!
+    else { // average
+      ds = 1.0-(intPar[2]*3.0)/(2.0*intPar[0])*fGapIntegral->IntegrateFunc(0); // integral prefactor is by 2 lower [Eqs.(19,20)] since intPar[0]==2kB T!
+      ds1 = 1.0-(intPar[2]*3.0)/(intPar[0])*fGapIntegral->IntegrateFunc(1); // integral prefactor is by 2 lower [Eqs.(19,20)] since intPar[0]==2kB T!
+      ds = (ds + 2.0 * sqrt(ds*ds1))/3.0; // since aa==bb the avg looks like this
+    }
+
+    intPar.clear();
+
+    if (newTemp)
+      fIntegralValues.push_back(ds);
+    else
+      fIntegralValues[vectorIndex] = ds;
+
+    fCalcNeeded[vectorIndex] = false;
+  }
+
+  return fIntegralValues[vectorIndex];
+}
+
+//--------------------------------------------------------------------
+/**
+ * <p>prepare the needed parameters for the integration carried out in TLinePWaveGapIntegralCuhre.
+ * For details see also the Memo GapIntegrals.pdf, especially Eq.(19) and (20).
+ */
+double TGapLinePWave::operator()(double t, const std::vector<double> &par) const {
+
+  // parameters: [0] Tc (K), [1] Delta0 (meV), [[2] orientation tag, [[3] c0 (1), [4] aG (1)]]
+
+  assert((par.size() >= 2) && (par.size() <= 5)); // 2 parameters: see A.~Carrington and F.~Manzano, Physica~C~\textbf{385}~(2003)~205
+                                                  // 4 or 5 parameters: see R. Prozorov and R. Giannetta, Supercond. Sci. Technol. 19 (2006) R41-R67
+                                                  // and Erratum Supercond. Sci. Technol. 21 (2008) 082003
+                                                  // c0 in the original context is c0 = (pi kB Tc) / Delta0
+                                                  // orientation tag: 0=aa,bb; 1=cc; 2=(sqrt[aa bb] + sqrt[aa cc] + sqrt[bb cc])/3 (default)
+  if (t <= 0.0)
+    return 1.0;
+  else if (t >= par[0])
+    return 0.0;
+
+
+  // check if orientation tag is given
+  int orientation_tag(2);
+  if ((par.size()==3) || (par.size()==5))
+    orientation_tag = static_cast<int>(par[2]);
+
+  bool integralParChanged(false);  
+
+  if (fPar.empty()) { // first time calling this routine
+    fPar = par;
+    integralParChanged = true;
+  } else { // check if parameter have changed
+    for (unsigned int i(0); i<par.size(); i++) {
+      if (par[i] != fPar[i]) {
+        fPar[i] = par[i];
+        integralParChanged = true;
+      }
+    }
+  }
+
+  bool newTemp(false);
+  unsigned int vectorIndex;
+
+  if (integralParChanged) {
+    fCalcNeeded.clear();
+    fCalcNeeded.resize(fTemp.size(), true);
+  }
+
+  fTempIter = find(fTemp.begin(), fTemp.end(), t);
+  if(fTempIter == fTemp.end()) {
+    fTemp.push_back(t);
+    vectorIndex = fTemp.size() - 1;
+    fCalcNeeded.push_back(true);
+    newTemp = true;
+  } else {
+    vectorIndex = fTempIter - fTemp.begin();
+  }
+
+  if (fCalcNeeded[vectorIndex]) {
+    double ds, ds1;
+    std::vector<double> intPar; // parameters for the integral, T & Delta(T)
+    intPar.push_back(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K
+    if ((par.size() == 2) || (par.size() == 3)) { // Carrington/Manzano
+      intPar.push_back(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51)));
+    } else { // Prozorov/Giannetta
+      intPar.push_back(par[1]*tanh(par[3]*sqrt(par[4]*(par[0]/t-1.0)))); // Delta0*tanh(c0*sqrt(aG*(Tc/T-1)))
+    }
+    intPar.push_back(4.0*(t+intPar[1])); // upper limit of energy-integration: cutoff energy
+    intPar.push_back(1.0); // upper limit of z-integration
+
+    fGapIntegral->SetParameters(intPar);
+    if (orientation_tag == 0) // aa,bb
+      ds = 1.0-(intPar[2]*3.0)/(2.0*intPar[0])*fGapIntegral->IntegrateFunc(0); // integral prefactor is by 2 lower [Eqs.(19,20)] since intPar[0]==2kB T!
+    else if (orientation_tag == 1) // cc
+      ds = 1.0-(intPar[2]*3.0)/(intPar[0])*fGapIntegral->IntegrateFunc(1); // integral prefactor is by 2 lower [Eqs.(19,20)] since intPar[0]==2kB T!
+    else { // average
+      ds = 1.0-(intPar[2]*3.0)/(2.0*intPar[0])*fGapIntegral->IntegrateFunc(0); // integral prefactor is by 2 lower [Eqs.(19,20)] since intPar[0]==2kB T!
+      ds1 = 1.0-(intPar[2]*3.0)/(intPar[0])*fGapIntegral->IntegrateFunc(1); // integral prefactor is by 2 lower [Eqs.(19,20)] since intPar[0]==2kB T!
+      ds = (ds + 2.0 * sqrt(ds*ds1))/3.0; // since aa==bb the avg looks like this
+    }
+
+    intPar.clear();
+
+    if (newTemp)
+      fIntegralValues.push_back(ds);
+    else
+      fIntegralValues[vectorIndex] = ds;
+
+    fCalcNeeded[vectorIndex] = false;
+  }
+
+  return fIntegralValues[vectorIndex];
+}
+
 //--------------------------------------------------------------------
 /**
  * <p>prepare the needed parameters for the integration carried out in TDWaveGapIntegralCuhre.
@@ -609,7 +930,6 @@ double TGapDWave::operator()(double t, const std::vector<double> &par) const {
   }
 
   return fIntegralValues[vectorIndex];
-
 }
 
 //--------------------------------------------------------------------
@@ -1103,7 +1423,40 @@ double TLambdaSWave::operator()(double t, const std::vector<double> &par) const
     return 1.0;
 
   return 1.0/sqrt((*fLambdaInvSq)(t, par));
+}
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+double TLambdaPointPWave::operator()(double t, const std::vector<double> &par) const
+{
+  assert(par.size() == 2); // two parameters: Tc, Delta0
+
+  if (t >= par[0])
+    return -1.0;
+
+  if (t <= 0.0)
+    return 1.0;
+
+  return 1.0/sqrt((*fLambdaInvSq)(t, par));
+}
+
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+double TLambdaLinePWave::operator()(double t, const std::vector<double> &par) const
+{
+  assert(par.size() == 2); // two parameters: Tc, Delta0
+
+  if (t >= par[0])
+    return -1.0;
+
+  if (t <= 0.0)
+    return 1.0;
+
+  return 1.0/sqrt((*fLambdaInvSq)(t, par));
 }
 
 //--------------------------------------------------------------------
@@ -1121,7 +1474,6 @@ double TLambdaDWave::operator()(double t, const std::vector<double> &par) const
     return 1.0;
 
   return 1.0/sqrt((*fLambdaInvSq)(t, par));
-
 }
 
 //--------------------------------------------------------------------
@@ -1139,7 +1491,6 @@ double TLambdaAnSWave::operator()(double t, const std::vector<double> &par) cons
     return 1.0;
 
   return 1.0/sqrt((*fLambdaInvSq)(t, par));
-
 }
 
 //--------------------------------------------------------------------
@@ -1157,7 +1508,6 @@ double TLambdaNonMonDWave1::operator()(double t, const std::vector<double> &par)
     return 1.0;
 
   return 1.0/sqrt((*fLambdaInvSq)(t, par));
-
 }
 
 //--------------------------------------------------------------------
@@ -1175,7 +1525,6 @@ double TLambdaNonMonDWave2::operator()(double t, const std::vector<double> &par)
     return 1.0;
 
   return 1.0/sqrt((*fLambdaInvSq)(t, par));
-
 }
 
 //--------------------------------------------------------------------
@@ -1192,7 +1541,6 @@ double TLambdaPowerLaw::operator()(double t, const std::vector<double> &par) con
     return -1.0;
 
   return 1.0/sqrt(1.0 - pow(t/par[0], par[1]));
-
 }
 
 
@@ -1211,7 +1559,40 @@ double TLambdaInvSWave::operator()(double t, const std::vector<double> &par) con
     return 1.0;
 
   return sqrt((*fLambdaInvSq)(t, par));
+}
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+double TLambdaInvPointPWave::operator()(double t, const std::vector<double> &par) const
+{
+  assert(par.size() == 2); // two parameters: Tc, Delta0
+
+  if (t >= par[0])
+    return 0.0;
+
+  if (t <= 0.0)
+    return 1.0;
+
+  return sqrt((*fLambdaInvSq)(t, par));
+}
+
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+double TLambdaInvLinePWave::operator()(double t, const std::vector<double> &par) const
+{
+  assert(par.size() == 2); // two parameters: Tc, Delta0
+
+  if (t >= par[0])
+    return 0.0;
+
+  if (t <= 0.0)
+    return 1.0;
+
+  return sqrt((*fLambdaInvSq)(t, par));
 }
 
 //--------------------------------------------------------------------
@@ -1229,7 +1610,6 @@ double TLambdaInvDWave::operator()(double t, const std::vector<double> &par) con
     return 1.0;
 
   return sqrt((*fLambdaInvSq)(t, par));
-
 }
 
 //--------------------------------------------------------------------
@@ -1247,7 +1627,6 @@ double TLambdaInvAnSWave::operator()(double t, const std::vector<double> &par) c
     return 1.0;
 
   return sqrt((*fLambdaInvSq)(t, par));
-
 }
 
 //--------------------------------------------------------------------
@@ -1265,7 +1644,6 @@ double TLambdaInvNonMonDWave1::operator()(double t, const std::vector<double> &p
     return 1.0;
 
   return sqrt((*fLambdaInvSq)(t, par));
-
 }
 
 //--------------------------------------------------------------------
@@ -1283,7 +1661,6 @@ double TLambdaInvNonMonDWave2::operator()(double t, const std::vector<double> &p
     return 1.0;
 
   return sqrt((*fLambdaInvSq)(t, par));
-
 }
 
 //--------------------------------------------------------------------
@@ -1300,7 +1677,6 @@ double TLambdaInvPowerLaw::operator()(double t, const std::vector<double> &par)
     return 0.0;
 
   return sqrt(1.0 - pow(t/par[0], par[1]));
-
 }
 
 //--------------------------------------------------------------------

src/external/libGapIntegrals/TGapIntegrals.h

@@ -63,6 +63,62 @@ private:
   ClassDef(TGapSWave,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+class TGapPointPWave : public PUserFcnBase {
+
+public:
+  TGapPointPWave();
+  virtual ~TGapPointPWave();
+
+  virtual Bool_t NeedGlobalPart() const { return false; }
+  virtual void SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) { }
+  virtual Bool_t GlobalPartIsValid() const { return true; }
+
+  double operator()(double, const std::vector<double>&) const;
+
+private:
+  TPointPWaveGapIntegralCuhre *fGapIntegral;
+  mutable std::vector<double> fTemp;
+  mutable std::vector<double>::const_iterator fTempIter;
+  mutable std::vector<double> fIntegralValues;
+  mutable std::vector<bool> fCalcNeeded;
+
+  mutable std::vector<double> fPar;
+
+  ClassDef(TGapPointPWave,1)
+};
+
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+class TGapLinePWave : public PUserFcnBase {
+
+public:
+  TGapLinePWave();
+  virtual ~TGapLinePWave();
+
+  virtual Bool_t NeedGlobalPart() const { return false; }
+  virtual void SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) { }
+  virtual Bool_t GlobalPartIsValid() const { return true; }
+
+  double operator()(double, const std::vector<double>&) const;
+
+private:
+  TLinePWaveGapIntegralCuhre *fGapIntegral;
+  mutable std::vector<double> fTemp;
+  mutable std::vector<double>::const_iterator fTempIter;
+  mutable std::vector<double> fIntegralValues;
+  mutable std::vector<bool> fCalcNeeded;
+
+  mutable std::vector<double> fPar;
+
+  ClassDef(TGapLinePWave,1)
+};
+
 //--------------------------------------------------------------------
 /**
  * <p>
@@ -297,6 +353,50 @@ private:
   ClassDef(TLambdaSWave,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+class TLambdaPointPWave : public PUserFcnBase {
+
+public:
+  TLambdaPointPWave();
+  virtual ~TLambdaPointPWave();
+
+  virtual Bool_t NeedGlobalPart() const { return false; }
+  virtual void SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) { }
+  virtual Bool_t GlobalPartIsValid() const { return true; }
+
+  double operator()(double, const std::vector<double>&) const;
+
+private:
+  TGapPointPWave *fLambdaInvSq;
+
+  ClassDef(TLambdaPointPWave,1)
+};
+
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+class TLambdaLinePWave : public PUserFcnBase {
+
+public:
+  TLambdaLinePWave();
+  virtual ~TLambdaLinePWave();
+
+  virtual Bool_t NeedGlobalPart() const { return false; }
+  virtual void SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) { }
+  virtual Bool_t GlobalPartIsValid() const { return true; }
+
+  double operator()(double, const std::vector<double>&) const;
+
+private:
+  TGapLinePWave *fLambdaInvSq;
+
+  ClassDef(TLambdaLinePWave,1)
+};
+
 //--------------------------------------------------------------------
 /**
  * <p>
@@ -428,6 +528,50 @@ private:
   ClassDef(TLambdaInvSWave,1)
 };
 
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+class TLambdaInvPointPWave : public PUserFcnBase {
+
+public:
+  TLambdaInvPointPWave();
+  virtual ~TLambdaInvPointPWave();
+
+  virtual Bool_t NeedGlobalPart() const { return false; }
+  virtual void SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) { }
+  virtual Bool_t GlobalPartIsValid() const { return true; }
+
+  double operator()(double, const std::vector<double>&) const;
+
+private:
+  TGapPointPWave *fLambdaInvSq;
+
+  ClassDef(TLambdaInvPointPWave,1)
+};
+
+//--------------------------------------------------------------------
+/**
+ * <p>
+ */
+class TLambdaInvLinePWave : public PUserFcnBase {
+
+public:
+  TLambdaInvLinePWave();
+  virtual ~TLambdaInvLinePWave();
+
+  virtual Bool_t NeedGlobalPart() const { return false; }
+  virtual void SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) { }
+  virtual Bool_t GlobalPartIsValid() const { return true; }
+
+  double operator()(double, const std::vector<double>&) const;
+
+private:
+  TGapLinePWave *fLambdaInvSq;
+
+  ClassDef(TLambdaInvLinePWave,1)
+};
+
 //--------------------------------------------------------------------
 /**
  * <p>

src/external/libGapIntegrals/TGapIntegralsLinkDef.h

@@ -35,6 +35,8 @@
 #pragma link off all functions;
 
 #pragma link C++ class TGapSWave+;
+#pragma link C++ class TGapPointPWave+;
+#pragma link C++ class TGapLinePWave+;
 #pragma link C++ class TGapDWave+;
 #pragma link C++ class TGapCosSqDWave+;
 #pragma link C++ class TGapSinSqDWave+;
@@ -44,12 +46,16 @@
 #pragma link C++ class TGapPowerLaw+;
 #pragma link C++ class TGapDirtySWave+;
 #pragma link C++ class TLambdaSWave+;
+#pragma link C++ class TLambdaPointPWave+;
+#pragma link C++ class TLambdaLinePWave+;
 #pragma link C++ class TLambdaDWave+;
 #pragma link C++ class TLambdaAnSWave+;
 #pragma link C++ class TLambdaNonMonDWave1+;
 #pragma link C++ class TLambdaNonMonDWave2+;
 #pragma link C++ class TLambdaPowerLaw+;
 #pragma link C++ class TLambdaInvSWave+;
+#pragma link C++ class TLambdaInvPointPWave+;
+#pragma link C++ class TLambdaInvLinePWave+;
 #pragma link C++ class TLambdaInvDWave+;
 #pragma link C++ class TLambdaInvAnSWave+;
 #pragma link C++ class TLambdaInvNonMonDWave1+;