diff --git a/src/external/BMWtools/BMWIntegrator.cpp b/src/external/BMWtools/BMWIntegrator.cpp index fb8ead3a..dfbc63b2 100644 --- a/src/external/BMWtools/BMWIntegrator.cpp +++ b/src/external/BMWtools/BMWIntegrator.cpp @@ -47,7 +47,7 @@ std::vector TPointPWaveGapIntegralCuhre::fPar; *

return: * - 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, - EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, - KEY, STATEFILE, SPIN, - &nregions, &neval, &fail, integral, error, prob); + 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]; } @@ -74,6 +80,7 @@ double TPointPWaveGapIntegralCuhre::IntegrateFunc() //----------------------------------------------------------------------------- /** *

Calculate the function value for the use with Cuhre---actual implementation of the function + * for p-wave point, aa==bb component * *

return: * - 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[], - const int *ncomp, double f[], void *userdata) // x = {E, theta}, fPar = {twokBT, Delta(T), Ec, thetac} +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 sinTheta = TMath::Sin(x[1]*fPar[3]); - double deltasq(TMath::Power(fPar[1]*sinTheta,2.0)); - f[0] = sinTheta/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0); + 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; +} + +//----------------------------------------------------------------------------- +/** + *

Calculate the function value for the use with Cuhre---actual implementation of the function + * for p-wave point, cc component + * + *

return: + * - 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 TLinePWaveGapIntegralCuhre::fPar; *

return: * - 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, - EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, - KEY, STATEFILE, SPIN, - &nregions, &neval, &fail, integral, error, prob); + 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]; } @@ -131,6 +167,7 @@ double TLinePWaveGapIntegralCuhre::IntegrateFunc() //----------------------------------------------------------------------------- /** *

Calculate the function value for the use with Cuhre---actual implementation of the function + * for p-wave line, aa==bb component * *

return: * - 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[], - const int *ncomp, double f[], void *userdata) // x = {E, theta}, fPar = {twokBT, Delta(T), Ec, thetac} +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 sinTheta = TMath::Sin(x[1]*fPar[3]); - double cosTheta = TMath::Cos(x[1]*fPar[3]); - double deltasq(TMath::Power(fPar[1]*cosTheta,2.0)); - f[0] = sinTheta/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0); + 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; +} + +//----------------------------------------------------------------------------- +/** + *

Calculate the function value for the use with Cuhre---actual implementation of the function + * for p-wave line, cc component + * + *

return: + * - 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; } diff --git a/src/external/BMWtools/BMWIntegrator.h b/src/external/BMWtools/BMWIntegrator.h index 5ecb7406..067e7d35 100644 --- a/src/external/BMWtools/BMWIntegrator.h +++ b/src/external/BMWtools/BMWIntegrator.h @@ -259,8 +259,9 @@ class TPointPWaveGapIntegralCuhre { TPointPWaveGapIntegralCuhre() : fNDim(2) {} ~TPointPWaveGapIntegralCuhre() { fPar.clear(); } void SetParameters(const std::vector &par) { fPar=par; } - static int Integrand(const int*, const double[], const int*, double[], void*); - double IntegrateFunc(); + 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 fPar; ///< parameters of the integrand @@ -278,8 +279,9 @@ class TLinePWaveGapIntegralCuhre { TLinePWaveGapIntegralCuhre() : fNDim(2) {} ~TLinePWaveGapIntegralCuhre() { fPar.clear(); } void SetParameters(const std::vector &par) { fPar=par; } - static int Integrand(const int*, const double[], const int*, double[], void*); - double IntegrateFunc(); + 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 fPar; ///< parameters of the integrand diff --git a/src/external/libGapIntegrals/GapIntegrals.pdf b/src/external/libGapIntegrals/GapIntegrals.pdf index 5a13b18c..37334a47 100644 Binary files a/src/external/libGapIntegrals/GapIntegrals.pdf and b/src/external/libGapIntegrals/GapIntegrals.pdf differ diff --git a/src/external/libGapIntegrals/GapIntegrals.tex b/src/external/libGapIntegrals/GapIntegrals.tex index 2dafd92b..c319c67d 100644 --- a/src/external/libGapIntegrals/GapIntegrals.tex +++ b/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] - 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, - the temperature dependence according to Eq.(\ref{eq:gapT_Prozorov}) will be used, otherwise Eq.(\ref{eq:gapT_Manzano}) will be utilized. + \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(\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] - 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, - the temperature dependence according to Eq.(\ref{eq:gapT_Prozorov}) will be used, otherwise Eq.(\ref{eq:gapT_Manzano}) will be utilized. + \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. @@ -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. diff --git a/src/external/libGapIntegrals/TGapIntegrals.cpp b/src/external/libGapIntegrals/TGapIntegrals.cpp index b90ee836..36c06122 100644 --- a/src/external/libGapIntegrals/TGapIntegrals.cpp +++ b/src/external/libGapIntegrals/TGapIntegrals.cpp @@ -671,17 +671,23 @@ double TGapSWave::operator()(double t, const std::vector &par) const { */ double TGapPointPWave::operator()(double t, const std::vector &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 - // 4 parameters: see R. Prozorov and R. Giannetta, Supercond. Sci. Technol. 19 (2006) R41-R67 + 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(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 &par) cons } if (fCalcNeeded[vectorIndex]) { - double ds; + double ds, ds1; std::vector 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 &par) cons //-------------------------------------------------------------------- /** *

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 &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 - // 4 parameters: see R. Prozorov and R. Giannetta, Supercond. Sci. Technol. 19 (2006) R41-R67 + 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; - bool integralParChanged(false); + + // check if orientation tag is given + int orientation_tag(2); + if ((par.size()==3) || (par.size()==5)) + orientation_tag = static_cast(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) const } if (fCalcNeeded[vectorIndex]) { - double ds; + double ds, ds1; std::vector 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();