diff --git a/doc/html/user/MUSR/BmwLibs.html b/doc/html/user/MUSR/BmwLibs.html index da9601fe..ffbae1c8 100644 --- a/doc/html/user/MUSR/BmwLibs.html +++ b/doc/html/user/MUSR/BmwLibs.html @@ -1,6 +1,6 @@ - + @@ -8,7 +8,7 @@ - + @@ -105,7 +105,7 @@ pre {
-
+

@@ -144,7 +144,7 @@ pre {
Topic revision: r5 - 10 Jul 2011, wojek
@@ -195,7 +195,7 @@ pre {
- +

diff --git a/doc/html/user/MUSR/LibFitPofB.html b/doc/html/user/MUSR/LibFitPofB.html index a37cfcf8..a1fe1a17 100644 --- a/doc/html/user/MUSR/LibFitPofB.html +++ b/doc/html/user/MUSR/LibFitPofB.html @@ -1,6 +1,6 @@ - + @@ -8,7 +8,7 @@ - + @@ -105,7 +105,7 @@ pre {
-
+

@@ -420,7 +420,7 @@ An example XML file looks as follows:
Topic revision: r16 - 10 Jul 2011, wojek
@@ -471,7 +471,7 @@ An example XML file looks as follows:
- +

diff --git a/doc/html/user/MUSR/LibZFRelaxation.html b/doc/html/user/MUSR/LibZFRelaxation.html index bffa062e..7faed20a 100644 --- a/doc/html/user/MUSR/LibZFRelaxation.html +++ b/doc/html/user/MUSR/LibZFRelaxation.html @@ -1,6 +1,6 @@ - + @@ -8,7 +8,7 @@ - + @@ -105,7 +105,7 @@ pre {
-
+

@@ -225,7 +225,7 @@ The parameters are:
Topic revision: r2 - 10 Jul 2011, wojek
@@ -276,7 +276,7 @@ The parameters are:
- +

diff --git a/doc/html/user/MUSR/Msr2Data.html b/doc/html/user/MUSR/Msr2Data.html index 99cb6f86..967c12f3 100644 --- a/doc/html/user/MUSR/Msr2Data.html +++ b/doc/html/user/MUSR/Msr2Data.html @@ -1,6 +1,6 @@ - + @@ -8,7 +8,7 @@ - + @@ -105,7 +105,7 @@ pre {
-
+

@@ -349,7 +349,7 @@ For reporting bugs or requesting new features and improvements please use the - @@ -422,7 +422,7 @@ For reporting bugs or requesting new features and improvements please use the
- +

diff --git a/doc/html/user/MUSR/MusrFit.html b/doc/html/user/MUSR/MusrFit.html index 6b8d44c7..a3851131 100644 --- a/doc/html/user/MUSR/MusrFit.html +++ b/doc/html/user/MUSR/MusrFit.html @@ -1,6 +1,6 @@ - + @@ -8,7 +8,7 @@ - + @@ -105,7 +105,7 @@ pre {
-
+

@@ -787,20 +787,20 @@ The RUN block is used to collect the data needed for a particular run to be fitt The tokens following the RUN statement are used to identify the run, the potential location where the run might be found, and the file format in which the run data has been saved. In order to understand the meaning of all the above tokens, a short digression is needed.

-Where is musrfit looking for data files? There is a specific order how this is done:
    +Where is musrfit looking for data files? There is a specific order how this is done:
    1. Check if the file is found in the current directory
    2. Check if the path (or multiple paths) was (were) given in the XML startup file.
    3. Check if there is a system variable MUSRFULLDATAPATH. This system variable can contain multiple search paths separated by colons, e.g.
       export MUSRFULLDATAPATH=/mnt/data/nemu/wkm/:/mnt/data/nemu/his/:/afs/psi.ch/user/s/smith/
       
    4. Construct the search path from the RUN-block information in the following way: Based on the RUN line in the RUN block, default paths will be generated, e.g. for
      -RUN lem07_his_2018 MUE4 PSI ROOT-NPP
      +RUN 2007/lem07_his_2018 MUE4 PSI ROOT-NPP
       
      the generated search path will look like
      -musrFullDataPathToken/DATA/Facility/Beamline/Year/runName.ext
      +musrFullDataPathToken/DATA/Facility/Beamline/runName.ext
       
      where musrFullDataPathToken is extracted from the MUSRFULLDATAPATH token by token, for the above example this might lead to the path
       /afs/psi.ch/user/s/smith/DATA/PSI/MUE4/2007/lem07_his_2018.root
       
      -
+

Here are some valid examples for the first line of a RUN block:
 RUN 2007/lem07_his_2018 MUE4 PSI ROOT-NPP
@@ -1453,7 +1453,7 @@ For reporting bugs or requesting new features and improvements please use the 
 
 
-
Topic revision: r114 - 22 Aug 2014, AndreasSuter
@@ -1526,7 +1526,7 @@ For reporting bugs or requesting new features and improvements please use the
- +

diff --git a/doc/html/user/MUSR/MusrFitAcknowledgements.html b/doc/html/user/MUSR/MusrFitAcknowledgements.html index 9b92479b..dad31e35 100644 --- a/doc/html/user/MUSR/MusrFitAcknowledgements.html +++ b/doc/html/user/MUSR/MusrFitAcknowledgements.html @@ -1,6 +1,6 @@ - + @@ -8,7 +8,7 @@ - + @@ -105,7 +105,7 @@ pre {
-
+

@@ -142,7 +142,7 @@ pre {
Topic revision: r4 - 10 Jul 2011, wojek
@@ -193,7 +193,7 @@ pre {
- +

diff --git a/doc/html/user/MUSR/MusrFitSetup.html b/doc/html/user/MUSR/MusrFitSetup.html index 1bac6475..42a91445 100644 --- a/doc/html/user/MUSR/MusrFitSetup.html +++ b/doc/html/user/MUSR/MusrFitSetup.html @@ -1,6 +1,6 @@ - + @@ -8,7 +8,7 @@ - + @@ -105,7 +105,7 @@ pre {
-
+

@@ -830,7 +830,7 @@ musrview test-histo-ROOT-NPP.msr
Topic revision: r56 - 09 May 2014, AndreasSuter
@@ -881,7 +881,7 @@ musrview test-histo-ROOT-NPP.msr
- +

diff --git a/doc/html/user/MUSR/MusrGui.html b/doc/html/user/MUSR/MusrGui.html index 2b24697b..4a533721 100644 --- a/doc/html/user/MUSR/MusrGui.html +++ b/doc/html/user/MUSR/MusrGui.html @@ -1,6 +1,6 @@ - + @@ -8,7 +8,7 @@ - + @@ -105,7 +105,7 @@ pre {
-
+

@@ -324,7 +324,7 @@ For reporting bugs or requesting new features and improvements please use the - @@ -397,7 +397,7 @@ For reporting bugs or requesting new features and improvements please use the
- +

diff --git a/doc/html/user/MUSR/QuickStart.html b/doc/html/user/MUSR/QuickStart.html index 418feb0d..7ac2c6f0 100644 --- a/doc/html/user/MUSR/QuickStart.html +++ b/doc/html/user/MUSR/QuickStart.html @@ -1,6 +1,6 @@ - + @@ -8,7 +8,7 @@ - + @@ -105,7 +105,7 @@ pre {
-
+

@@ -280,7 +280,7 @@ RUN 2008/lem08_his_8472 MUE4 PSI ROOT-NPP (name beamline institute dat
Topic revision: r7 - 10 Jul 2011, wojek
@@ -331,7 +331,7 @@ RUN 2008/lem08_his_8472 MUE4 PSI ROOT-NPP (name beamline institute dat
- +

diff --git a/doc/html/user/MUSR/TutorialSingleHisto.html b/doc/html/user/MUSR/TutorialSingleHisto.html index 83df529c..e354fc6d 100644 --- a/doc/html/user/MUSR/TutorialSingleHisto.html +++ b/doc/html/user/MUSR/TutorialSingleHisto.html @@ -1,6 +1,6 @@ - + @@ -8,7 +8,7 @@ - + @@ -105,7 +105,7 @@ pre {
-
+

@@ -281,7 +281,7 @@ This page only summarizes the very basic features and options of the programs co -
+

Attachments (6)

 
@@ -303,7 +303,7 @@ This page only summarizes the very basic features and options of the programs co
Topic revision: r9 - 02 Sep 2011, wojek
@@ -354,7 +354,7 @@ This page only summarizes the very basic features and options of the programs co
- +

diff --git a/doc/html/user/MUSR/WebHome.html b/doc/html/user/MUSR/WebHome.html index d4929216..3ae90265 100644 --- a/doc/html/user/MUSR/WebHome.html +++ b/doc/html/user/MUSR/WebHome.html @@ -1,6 +1,6 @@ - + @@ -8,7 +8,7 @@ - + @@ -105,7 +105,7 @@ pre {
-
+

@@ -156,7 +156,7 @@ pre {
Topic revision: r43 - 09 May 2014, AndreasSuter
@@ -207,7 +207,7 @@ pre {
- +

diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index 42dfd9f0..66dba079 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -579,7 +579,8 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) sstr.Remove(TString::kLeading, ' '); if (str.BeginsWith("fun")) { if (FilterNumber(sstr, "fun", 0, number)) { - sstr = fFuncHandler->GetFuncString(number-1); + idx = GetFuncIndex(number); // get index of the function number + sstr = fFuncHandler->GetFuncString(idx); sstr.ToLower(); fout << sstr.Data() << endl; } @@ -1228,6 +1229,8 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) fout << str.Data() << endl; if (messages) cout << endl << str.Data() << endl; + + // check if per run block maxLH needs to be written } } else { fout << "*** FIT DID NOT CONVERGE ***" << endl; diff --git a/src/classes/PRunDataHandler.cpp b/src/classes/PRunDataHandler.cpp index bfae557f..08a7a4be 100644 --- a/src/classes/PRunDataHandler.cpp +++ b/src/classes/PRunDataHandler.cpp @@ -1057,15 +1057,11 @@ Bool_t PRunDataHandler::FileExistsCheck(PMsrRunBlock &runInfo, const UInt_t idx) } pstr->ToUpper(); runInfo.SetBeamline(*pstr, idx); - TDatime datetime; - TString dt; - dt += datetime.GetYear(); for (Int_t i=0; iGetEntries(); i++) { ostr = dynamic_cast(tokens->At(i)); str = ostr->GetString() + TString("/DATA/") + *runInfo.GetInstitute(idx) + TString("/") + *runInfo.GetBeamline(idx) + TString("/") + - dt + TString("/") + *runInfo.GetRunName(idx); TestFileName(str, ext); if (!str.IsNull()) { // found @@ -5167,7 +5163,15 @@ Bool_t PRunDataHandler::WritePsiBinFile(TString fln) vector svec; TString str, date; TDatime dt; - dt.Set(fData[0].GetStartDateTime()); + int year, month, day; + + // 28-Aug-2014, TP: the following line does not work, it generates the wrong date + //dt.Set(fData[0].GetStartDateTime()); + // the following generates the correct date entry + date.Append(*fData[0].GetStartDate()); + sscanf((const char*)date.Data(),"%04d-%02d-%02d", &year, &month, &day); + dt.Set(year, month, day, 0, 0, 0); + date.Form("%02d-", dt.GetDay()); date.Append(GetMonth(dt.GetMonth())); date.Append("-"); @@ -5183,7 +5187,14 @@ Bool_t PRunDataHandler::WritePsiBinFile(TString fln) svec.clear(); // run stop date - dt.Set(fData[0].GetStopDateTime()); + // 28-Aug-2014, TP: the following line does not work, it generates the wrong date + //dt.Set(fData[0].GetStopDateTime()); + // the following generates the correct date entry + date.Clear(); + date.Append(*fData[0].GetStopDate()); + sscanf((const char*)date.Data(),"%04d-%02d-%02d", &year, &month, &day); + dt.Set(year, month, day, 0, 0, 0); + date.Form("%02d-", dt.GetDay()); date.Append(GetMonth(dt.GetMonth())); date.Append("-"); diff --git a/src/external/libGapIntegrals/GapIntegrals.pdf b/src/external/libGapIntegrals/GapIntegrals.pdf index 86112612..80214db8 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 new file mode 100644 index 00000000..57f8bf4d --- /dev/null +++ b/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} diff --git a/src/external/libGapIntegrals/PSI-Logo_narrow.jpg b/src/external/libGapIntegrals/PSI-Logo_narrow.jpg new file mode 100644 index 00000000..9bda1610 Binary files /dev/null and b/src/external/libGapIntegrals/PSI-Logo_narrow.jpg differ diff --git a/src/external/libGapIntegrals/TGapIntegrals.cpp b/src/external/libGapIntegrals/TGapIntegrals.cpp index 448556d3..6a3545f8 100644 --- a/src/external/libGapIntegrals/TGapIntegrals.cpp +++ b/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 - +//-------------------------------------------------------------------- +/** + *

s wave gap integral + */ TGapSWave::TGapSWave() { TGapIntegral *gapint = new TGapIntegral(); fGapIntegral = gapint; @@ -78,6 +80,10 @@ TGapSWave::TGapSWave() { fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ TGapDWave::TGapDWave() { TDWaveGapIntegralCuhre *gapint = new TDWaveGapIntegralCuhre(); fGapIntegral = gapint; @@ -90,6 +96,10 @@ TGapDWave::TGapDWave() { fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ TGapCosSqDWave::TGapCosSqDWave() { TCosSqDWaveGapIntegralCuhre *gapint = new TCosSqDWaveGapIntegralCuhre(); fGapIntegral = gapint; @@ -102,6 +112,10 @@ TGapCosSqDWave::TGapCosSqDWave() { fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ TGapSinSqDWave::TGapSinSqDWave() { TSinSqDWaveGapIntegralCuhre *gapint = new TSinSqDWaveGapIntegralCuhre(); fGapIntegral = gapint; @@ -114,6 +128,10 @@ TGapSinSqDWave::TGapSinSqDWave() { fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ TGapAnSWave::TGapAnSWave() { TAnSWaveGapIntegralCuhre *gapint = new TAnSWaveGapIntegralCuhre(); fGapIntegral = gapint; @@ -126,6 +144,10 @@ TGapAnSWave::TGapAnSWave() { fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ TGapNonMonDWave1::TGapNonMonDWave1() { TNonMonDWave1GapIntegralCuhre *gapint = new TNonMonDWave1GapIntegralCuhre(); fGapIntegral = gapint; @@ -138,6 +160,10 @@ TGapNonMonDWave1::TGapNonMonDWave1() { fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ TGapNonMonDWave2::TGapNonMonDWave2() { TNonMonDWave2GapIntegralCuhre *gapint = new TNonMonDWave2GapIntegralCuhre(); fGapIntegral = gapint; @@ -150,46 +176,90 @@ TGapNonMonDWave2::TGapNonMonDWave2() { fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaSWave::TLambdaSWave() { fLambdaInvSq = new TGapSWave(); } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaDWave::TLambdaDWave() { fLambdaInvSq = new TGapDWave(); } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaAnSWave::TLambdaAnSWave() { fLambdaInvSq = new TGapAnSWave(); } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaNonMonDWave1::TLambdaNonMonDWave1() { fLambdaInvSq = new TGapNonMonDWave1(); } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaNonMonDWave2::TLambdaNonMonDWave2() { fLambdaInvSq = new TGapNonMonDWave2(); } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaInvSWave::TLambdaInvSWave() { fLambdaInvSq = new TGapSWave(); } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaInvDWave::TLambdaInvDWave() { fLambdaInvSq = new TGapDWave(); } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaInvAnSWave::TLambdaInvAnSWave() { fLambdaInvSq = new TGapAnSWave(); } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaInvNonMonDWave1::TLambdaInvNonMonDWave1() { fLambdaInvSq = new TGapNonMonDWave1(); } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaInvNonMonDWave2::TLambdaInvNonMonDWave2() { fLambdaInvSq = new TGapNonMonDWave2(); } +//-------------------------------------------------------------------- +/** + *

+ */ TGapSWave::~TGapSWave() { delete fGapIntegral; fGapIntegral = 0; @@ -201,6 +271,10 @@ TGapSWave::~TGapSWave() { fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ TGapDWave::~TGapDWave() { delete fGapIntegral; fGapIntegral = 0; @@ -212,6 +286,10 @@ TGapDWave::~TGapDWave() { fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ TGapCosSqDWave::~TGapCosSqDWave() { delete fGapIntegral; fGapIntegral = 0; @@ -223,6 +301,10 @@ TGapCosSqDWave::~TGapCosSqDWave() { fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ TGapSinSqDWave::~TGapSinSqDWave() { delete fGapIntegral; fGapIntegral = 0; @@ -234,6 +316,10 @@ TGapSinSqDWave::~TGapSinSqDWave() { fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ TGapAnSWave::~TGapAnSWave() { delete fGapIntegral; fGapIntegral = 0; @@ -245,6 +331,10 @@ TGapAnSWave::~TGapAnSWave() { fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ TGapNonMonDWave1::~TGapNonMonDWave1() { delete fGapIntegral; fGapIntegral = 0; @@ -256,6 +346,10 @@ TGapNonMonDWave1::~TGapNonMonDWave1() { fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ TGapNonMonDWave2::~TGapNonMonDWave2() { delete fGapIntegral; fGapIntegral = 0; @@ -267,60 +361,110 @@ TGapNonMonDWave2::~TGapNonMonDWave2() { fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaSWave::~TLambdaSWave() { delete fLambdaInvSq; fLambdaInvSq = 0; } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaDWave::~TLambdaDWave() { delete fLambdaInvSq; fLambdaInvSq = 0; } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaAnSWave::~TLambdaAnSWave() { delete fLambdaInvSq; fLambdaInvSq = 0; } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaNonMonDWave1::~TLambdaNonMonDWave1() { delete fLambdaInvSq; fLambdaInvSq = 0; } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaNonMonDWave2::~TLambdaNonMonDWave2() { delete fLambdaInvSq; fLambdaInvSq = 0; } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaInvSWave::~TLambdaInvSWave() { delete fLambdaInvSq; fLambdaInvSq = 0; } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaInvDWave::~TLambdaInvDWave() { delete fLambdaInvSq; fLambdaInvSq = 0; } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaInvAnSWave::~TLambdaInvAnSWave() { delete fLambdaInvSq; fLambdaInvSq = 0; } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaInvNonMonDWave1::~TLambdaInvNonMonDWave1() { delete fLambdaInvSq; fLambdaInvSq = 0; } +//-------------------------------------------------------------------- +/** + *

+ */ TLambdaInvNonMonDWave2::~TLambdaInvNonMonDWave2() { delete fLambdaInvSq; fLambdaInvSq = 0; } +//-------------------------------------------------------------------- +/** + *

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 &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 &par) const { double ds; vector 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(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51))); + intPar.push_back(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K + 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 &par) const { } +//-------------------------------------------------------------------- +/** + *

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 &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 &par) const { double ds; vector 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(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51))); - 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(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K + 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(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 &par) const { } +//-------------------------------------------------------------------- +/** + *

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 &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 &par) const { double ds; vector 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(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51))); - 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(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K + 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(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 &par) const { } +//-------------------------------------------------------------------- +/** + *

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 &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 &par) const { double ds; vector 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(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51))); - 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(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K + 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(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 &par) const { } +//-------------------------------------------------------------------- +/** + *

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 &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 &par) const { double ds; vector 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(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51))); + intPar.push_back(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K + 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 &par) const { } +//-------------------------------------------------------------------- +/** + *

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 &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 &par) const { double ds; vector 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(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51))); + intPar.push_back(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K + 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 &par) const { } +//-------------------------------------------------------------------- +/** + *

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 &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 &par) const { double ds; vector 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(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51))); + intPar.push_back(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K + 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 &par) const { } +//-------------------------------------------------------------------- +/** + *

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 &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 &par) const { } +//-------------------------------------------------------------------- +/** + *

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 &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))); - - 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 + 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/(0.172346648*t)); // Delta(T)/Delta(0)*tanh(Delta(T)/2 kB T), kB in meV/K } +//-------------------------------------------------------------------- +/** + *

+ */ double TLambdaSWave::operator()(double t, const vector &par) const { assert(par.size() == 2); // two parameters: Tc, Delta0 @@ -856,6 +1079,10 @@ double TLambdaSWave::operator()(double t, const vector &par) const } +//-------------------------------------------------------------------- +/** + *

+ */ double TLambdaDWave::operator()(double t, const vector &par) const { assert(par.size() == 2); // two parameters: Tc, Delta0 @@ -870,6 +1097,10 @@ double TLambdaDWave::operator()(double t, const vector &par) const } +//-------------------------------------------------------------------- +/** + *

+ */ double TLambdaAnSWave::operator()(double t, const vector &par) const { assert(par.size() == 3); // three parameters: Tc, Delta0, a @@ -884,6 +1115,10 @@ double TLambdaAnSWave::operator()(double t, const vector &par) const } +//-------------------------------------------------------------------- +/** + *

+ */ double TLambdaNonMonDWave1::operator()(double t, const vector &par) const { assert(par.size() == 3); // three parameters: Tc, Delta0, a @@ -898,6 +1133,10 @@ double TLambdaNonMonDWave1::operator()(double t, const vector &par) cons } +//-------------------------------------------------------------------- +/** + *

+ */ double TLambdaNonMonDWave2::operator()(double t, const vector &par) const { assert(par.size() == 3); // three parameters: Tc, Delta0, a @@ -912,6 +1151,10 @@ double TLambdaNonMonDWave2::operator()(double t, const vector &par) cons } +//-------------------------------------------------------------------- +/** + *

+ */ double TLambdaPowerLaw::operator()(double t, const vector &par) const { assert(par.size() == 2); // two parameters: Tc, N @@ -926,6 +1169,10 @@ double TLambdaPowerLaw::operator()(double t, const vector &par) const { } +//-------------------------------------------------------------------- +/** + *

+ */ double TLambdaInvSWave::operator()(double t, const vector &par) const { assert(par.size() == 2); // two parameters: Tc, Delta0 @@ -940,6 +1187,10 @@ double TLambdaInvSWave::operator()(double t, const vector &par) const } +//-------------------------------------------------------------------- +/** + *

+ */ double TLambdaInvDWave::operator()(double t, const vector &par) const { assert(par.size() == 2); // two parameters: Tc, Delta0 @@ -954,6 +1205,10 @@ double TLambdaInvDWave::operator()(double t, const vector &par) const } +//-------------------------------------------------------------------- +/** + *

+ */ double TLambdaInvAnSWave::operator()(double t, const vector &par) const { assert(par.size() == 3); // three parameters: Tc, Delta0, a @@ -968,6 +1223,10 @@ double TLambdaInvAnSWave::operator()(double t, const vector &par) const } +//-------------------------------------------------------------------- +/** + *

+ */ double TLambdaInvNonMonDWave1::operator()(double t, const vector &par) const { assert(par.size() == 3); // three parameters: Tc, Delta0, a @@ -982,6 +1241,10 @@ double TLambdaInvNonMonDWave1::operator()(double t, const vector &par) c } +//-------------------------------------------------------------------- +/** + *

+ */ double TLambdaInvNonMonDWave2::operator()(double t, const vector &par) const { assert(par.size() == 3); // three parameters: Tc, Delta0, a @@ -996,6 +1259,10 @@ double TLambdaInvNonMonDWave2::operator()(double t, const vector &par) c } +//-------------------------------------------------------------------- +/** + *

+ */ double TLambdaInvPowerLaw::operator()(double t, const vector &par) const { assert(par.size() == 2); // two parameters: Tc, N @@ -1009,12 +1276,20 @@ double TLambdaInvPowerLaw::operator()(double t, const vector &par) const } +//-------------------------------------------------------------------- +/** + *

+ */ TFilmMagnetizationDWave::TFilmMagnetizationDWave() { fLambdaInvSq = new TGapDWave(); fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ TFilmMagnetizationDWave::~TFilmMagnetizationDWave() { delete fLambdaInvSq; @@ -1022,6 +1297,10 @@ TFilmMagnetizationDWave::~TFilmMagnetizationDWave() fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ double TFilmMagnetizationDWave::operator()(double t, const vector &par) const { assert(par.size() == 4); // four parameters: Tc, Delta0, lambda0, film-thickness diff --git a/src/external/libGapIntegrals/TGapIntegrals.h b/src/external/libGapIntegrals/TGapIntegrals.h index 53835148..603b08a9 100644 --- a/src/external/libGapIntegrals/TGapIntegrals.h +++ b/src/external/libGapIntegrals/TGapIntegrals.h @@ -37,6 +37,10 @@ using namespace std; #include "PUserFcnBase.h" #include "BMWIntegrator.h" +//-------------------------------------------------------------------- +/** + *

+ */ class TGapSWave : public PUserFcnBase { public: @@ -61,6 +65,10 @@ private: ClassDef(TGapSWave,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ class TGapDWave : public PUserFcnBase { public: @@ -85,6 +93,10 @@ private: ClassDef(TGapDWave,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ class TGapCosSqDWave : public PUserFcnBase { public: @@ -109,6 +121,10 @@ private: ClassDef(TGapCosSqDWave,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ class TGapSinSqDWave : public PUserFcnBase { public: @@ -133,7 +149,10 @@ private: ClassDef(TGapSinSqDWave,1) }; - +//-------------------------------------------------------------------- +/** + *

+ */ class TGapAnSWave : public PUserFcnBase { public: @@ -158,6 +177,10 @@ private: ClassDef(TGapAnSWave,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ class TGapNonMonDWave1 : public PUserFcnBase { public: @@ -182,6 +205,10 @@ private: ClassDef(TGapNonMonDWave1,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ class TGapNonMonDWave2 : public PUserFcnBase { public: @@ -207,6 +234,10 @@ private: }; +//-------------------------------------------------------------------- +/** + *

+ */ class TGapPowerLaw : public PUserFcnBase { public: @@ -224,6 +255,10 @@ private: ClassDef(TGapPowerLaw,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ class TGapDirtySWave : public PUserFcnBase { public: @@ -242,6 +277,10 @@ private: }; +//-------------------------------------------------------------------- +/** + *

+ */ class TLambdaSWave : public PUserFcnBase { public: @@ -260,6 +299,10 @@ private: ClassDef(TLambdaSWave,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ class TLambdaDWave : public PUserFcnBase { public: @@ -278,6 +321,10 @@ private: ClassDef(TLambdaDWave,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ class TLambdaAnSWave : public PUserFcnBase { public: @@ -296,6 +343,10 @@ private: ClassDef(TLambdaAnSWave,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ class TLambdaNonMonDWave1 : public PUserFcnBase { public: @@ -314,6 +365,10 @@ private: ClassDef(TLambdaNonMonDWave1,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ class TLambdaNonMonDWave2 : public PUserFcnBase { public: @@ -332,7 +387,10 @@ private: ClassDef(TLambdaNonMonDWave2,1) }; - +//-------------------------------------------------------------------- +/** + *

+ */ class TLambdaPowerLaw : public PUserFcnBase { public: @@ -350,6 +408,10 @@ private: ClassDef(TLambdaPowerLaw,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ class TLambdaInvSWave : public PUserFcnBase { public: @@ -368,6 +430,10 @@ private: ClassDef(TLambdaInvSWave,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ class TLambdaInvDWave : public PUserFcnBase { public: @@ -386,6 +452,10 @@ private: ClassDef(TLambdaInvDWave,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ class TLambdaInvAnSWave : public PUserFcnBase { public: @@ -404,6 +474,10 @@ private: ClassDef(TLambdaInvAnSWave,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ class TLambdaInvNonMonDWave1 : public PUserFcnBase { public: @@ -422,6 +496,10 @@ private: ClassDef(TLambdaInvNonMonDWave1,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ class TLambdaInvNonMonDWave2 : public PUserFcnBase { public: @@ -440,7 +518,10 @@ private: ClassDef(TLambdaInvNonMonDWave2,1) }; - +//-------------------------------------------------------------------- +/** + *

+ */ class TLambdaInvPowerLaw : public PUserFcnBase { public: @@ -458,6 +539,10 @@ private: ClassDef(TLambdaInvPowerLaw,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ class TFilmMagnetizationDWave : public PUserFcnBase { public: