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

src/external/BMWtools/BMWIntegrator.cpp

@@ -47,7 +47,7 @@ std::vector<double> TPointPWaveGapIntegralCuhre::fPar;
  * <p><b>return:</b>
  * - value of the integral
  */
-double TPointPWaveGapIntegralCuhre::IntegrateFunc()
+double TPointPWaveGapIntegralCuhre::IntegrateFunc(int tag)
 {
   const unsigned int NCOMP(1);
   const unsigned int NVEC(1);
@@ -63,10 +63,16 @@ double TPointPWaveGapIntegralCuhre::IntegrateFunc()
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Cuhre(fNDim, NCOMP, Integrand, USERDATA, NVEC,
+  if (tag == 0)
-    EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
+    Cuhre(fNDim, NCOMP, Integrand_aa, USERDATA, NVEC,
-    KEY, STATEFILE, SPIN,
+      EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
-    &nregions, &neval, &fail, integral, error, prob);
+      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];
 }
@@ -74,6 +80,7 @@ double TPointPWaveGapIntegralCuhre::IntegrateFunc()
 //-----------------------------------------------------------------------------
 /**
  * <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
@@ -84,12 +91,35 @@ double TPointPWaveGapIntegralCuhre::IntegrateFunc()
  * \param f function value
  * \param userdata additional user parameters (required by the interface, NULL here)
  */
-int TPointPWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
+int TPointPWaveGapIntegralCuhre::Integrand_aa(const int *ndim, const double x[],
-                      const int *ncomp, double f[], void *userdata) // x = {E, theta}, fPar = {twokBT, Delta(T), Ec, thetac}
+                      const int *ncomp, double f[], void *userdata) // x = {E, z}, fPar = {twokBT, Delta(T), Ec, zc}
 {
-  double sinTheta = TMath::Sin(x[1]*fPar[3]);
+  double z = x[1]*fPar[3];
-  double deltasq(TMath::Power(fPar[1]*sinTheta,2.0));
+  double deltasq(pow(sqrt(1.0-z*z)*fPar[1],2.0));
-  f[0] = sinTheta/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),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;
 }
 
@@ -104,7 +134,7 @@ std::vector<double> TLinePWaveGapIntegralCuhre::fPar;
  * <p><b>return:</b>
  * - value of the integral
  */
-double TLinePWaveGapIntegralCuhre::IntegrateFunc()
+double TLinePWaveGapIntegralCuhre::IntegrateFunc(int tag)
 {
   const unsigned int NCOMP(1);
   const unsigned int NVEC(1);
@@ -120,10 +150,16 @@ double TLinePWaveGapIntegralCuhre::IntegrateFunc()
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Cuhre(fNDim, NCOMP, Integrand, USERDATA, NVEC,
+  if (tag == 0)
-    EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
+    Cuhre(fNDim, NCOMP, Integrand_aa, USERDATA, NVEC,
-    KEY, STATEFILE, SPIN,
+      EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
-    &nregions, &neval, &fail, integral, error, prob);
+      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];
 }
@@ -131,6 +167,7 @@ double TLinePWaveGapIntegralCuhre::IntegrateFunc()
 //-----------------------------------------------------------------------------
 /**
  * <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
@@ -141,13 +178,35 @@ double TLinePWaveGapIntegralCuhre::IntegrateFunc()
  * \param f function value
  * \param userdata additional user parameters (required by the interface, NULL here)
  */
-int TLinePWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
+int TLinePWaveGapIntegralCuhre::Integrand_aa(const int *ndim, const double x[],
-                      const int *ncomp, double f[], void *userdata) // x = {E, theta}, fPar = {twokBT, Delta(T), Ec, thetac}
+                      const int *ncomp, double f[], void *userdata) // x = {E, z}, fPar = {twokBT, Delta(T), Ec, zc}
 {
-  double sinTheta = TMath::Sin(x[1]*fPar[3]);
+  double z = x[1]*fPar[3];
-  double cosTheta = TMath::Cos(x[1]*fPar[3]);
+  double deltasq(pow(z*fPar[1],2.0));
-  double deltasq(TMath::Power(fPar[1]*cosTheta,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);
-  f[0] = sinTheta/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;
 }
 

src/external/BMWtools/BMWIntegrator.h

@@ -259,8 +259,9 @@ class TPointPWaveGapIntegralCuhre {
     TPointPWaveGapIntegralCuhre() : fNDim(2) {}
     ~TPointPWaveGapIntegralCuhre() { fPar.clear(); }
     void SetParameters(const std::vector<double> &par) { fPar=par; }
-    static int Integrand(const int*, const double[], const int*, double[], void*);
+    static int Integrand_aa(const int*, const double[], const int*, double[], void*);
-    double IntegrateFunc();
+    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
@@ -278,8 +279,9 @@ class TLinePWaveGapIntegralCuhre {
     TLinePWaveGapIntegralCuhre() : fNDim(2) {}
     ~TLinePWaveGapIntegralCuhre() { fPar.clear(); }
     void SetParameters(const std::vector<double> &par) { fPar=par; }
-    static int Integrand(const int*, const double[], const int*, double[], void*);
+    static int Integrand_aa(const int*, const double[], const int*, double[], void*);
-    double IntegrateFunc();
+    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

src/external/libGapIntegrals/GapIntegrals.tex

@@ -170,18 +170,20 @@ The \gapint plug-in calculates $\tilde{I}(T)$ for the following $\Delta(\varphi)
     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(\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]?\\[1.5ex]
+   \musrfit theory line: \verb?userFcn libGapIntegrals TGapPointPWave 1 2 [3 [4 5]]?\\[1.5ex]
-    Parameters: $T_{\mathrm c}~(\mathrm{K})$, $\Delta_0~(\mathrm{meV})$, $[c_0~(1),~a_{\rm G}~(1)]$. If $c_0$ and $a_{\rm G}$ are provided, 
+    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.
+    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(\theta, T) = \Delta(T) \cos(\theta) = \Delta(T) \cdot z
    \end{equation}
-   \musrfit theory line: \verb?userFcn libGapIntegrals TGapLinePWave 1 2 [3 4]?\\[1.5ex]
+   \musrfit theory line: \verb?userFcn libGapIntegrals TGapLinePWave 1 2 [3 [4 5]]?\\[1.5ex]
-    Parameters: $T_{\mathrm c}~(\mathrm{K})$, $\Delta_0~(\mathrm{meV})$, $[c_0~(1),~a_{\rm G}~(1)]$. If $c_0$ and $a_{\rm G}$ are provided, 
+    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.
+    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. 
@@ -222,12 +224,12 @@ n_{aa \atop bb}(T) = 1 - \frac{1}{2\pi k_{\rm B} T} \int_0^{2\pi} \mathrm{d}\var
 
 \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) \atop \sin^2(\varphi)} \cdot G(\Delta(z,\varphi), T) \label{eq:n_anisotrope_3D_cc}
+ 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}
+\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}
 
@@ -252,7 +254,7 @@ n_{aa \atop bb}(T) = 1 - \frac{1}{2\pi k_{\rm B} T} \int_0^{2\pi} \mathrm{d}\var
 \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) \atop \sin^2(\varphi)}}_{=\pi} \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*}
 
@@ -260,6 +262,17 @@ n_{aa \atop bb}(T) = 1 - \frac{1}{2\pi k_{\rm B} T} \int_0^{2\pi} \mathrm{d}\var
 
 $$ 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.
 

src/external/libGapIntegrals/TGapIntegrals.cpp

@@ -671,17 +671,23 @@ double TGapSWave::operator()(double t, const std::vector<double> &par) const {
  */
 double TGapPointPWave::operator()(double t, const std::vector<double> &par) const {
 
-  // parameters: [0] Tc (K), [1] Delta0 (meV), [[2] c0 (1), [3] aG (1)]
+  // parameters: [0] Tc (K), [1] Delta0 (meV), [[2] orientation tag, [[3] c0 (1), [4] aG (1)]]
 
-  assert((par.size() == 2) || (par.size() == 4)); // 2 parameters: see A.~Carrington and F.~Manzano, Physica~C~\textbf{385}~(2003)~205
+  assert((par.size() >= 2) && (par.size() <= 5)); // 2 parameters: see A.~Carrington and F.~Manzano, Physica~C~\textbf{385}~(2003)~205
-                                                  // 4 parameters: see R. Prozorov and R. Giannetta, Supercond. Sci. Technol. 19 (2006) R41-R67
+                                                  // 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
@@ -715,19 +721,27 @@ double TGapPointPWave::operator()(double t, const std::vector<double> &par) cons
   }
 
   if (fCalcNeeded[vectorIndex]) {
-    double ds;
+    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) { // Carrington/Manzano
+    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(TMath::Pi()); // upper limit of theta-integration
+    intPar.push_back(1.0); // upper limit of theta-integration
 
     fGapIntegral->SetParameters(intPar);
-    ds = 1.0-(intPar[2]*PI)/(2.0*intPar[0])*fGapIntegral->IntegrateFunc();
+    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();
 
@@ -745,27 +759,34 @@ double TGapPointPWave::operator()(double t, const std::vector<double> &par) cons
 //--------------------------------------------------------------------
 /**
  * <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).
+ * 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] c0 (1), [3] aG (1)]
+  // parameters: [0] Tc (K), [1] Delta0 (meV), [[2] orientation tag, [[3] c0 (1), [4] aG (1)]]
 
-  assert((par.size() == 2) || (par.size() == 4)); // 2 parameters: see A.~Carrington and F.~Manzano, Physica~C~\textbf{385}~(2003)~205
+  assert((par.size() >= 2) && (par.size() <= 5)); // 2 parameters: see A.~Carrington and F.~Manzano, Physica~C~\textbf{385}~(2003)~205
-                                                  // 4 parameters: see R. Prozorov and R. Giannetta, Supercond. Sci. Technol. 19 (2006) R41-R67
+                                                  // 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
+  } else { // check if parameter have changed
     for (unsigned int i(0); i<par.size(); i++) {
       if (par[i] != fPar[i]) {
         fPar[i] = par[i];
@@ -793,19 +814,27 @@ double TGapLinePWave::operator()(double t, const std::vector<double> &par) const
   }
 
   if (fCalcNeeded[vectorIndex]) {
-    double ds;
+    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) { // Carrington/Manzano
+    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(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(TMath::Pi()); // upper limit of theta-integration
+    intPar.push_back(1.0); // upper limit of z-integration
 
     fGapIntegral->SetParameters(intPar);
-    ds = 1.0-(intPar[2]*PI)/(2.0*intPar[0])*fGapIntegral->IntegrateFunc();
+    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();