From 8551665524e51a2a122656486f7e850d035e5359 Mon Sep 17 00:00:00 2001 From: "Bastian M. Wojek" Date: Wed, 12 May 2010 14:53:35 +0000 Subject: [PATCH] Added dirty limit s-wave expression for the superfluid density --- src/external/BMWIntegrator/BMWIntegrator.cpp | 67 ++++++ src/external/BMWIntegrator/BMWIntegrator.h | 28 +++ .../libGapIntegrals/TGapIntegrals.cpp | 214 +++++++++++++++++- src/external/libGapIntegrals/TGapIntegrals.h | 55 +++++ .../libGapIntegrals/TGapIntegralsLinkDef.h | 3 + 5 files changed, 366 insertions(+), 1 deletion(-) diff --git a/src/external/BMWIntegrator/BMWIntegrator.cpp b/src/external/BMWIntegrator/BMWIntegrator.cpp index 7701eea2..220070ae 100644 --- a/src/external/BMWIntegrator/BMWIntegrator.cpp +++ b/src/external/BMWIntegrator/BMWIntegrator.cpp @@ -66,6 +66,73 @@ void TDWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[], return; } +std::vector TCosSqDWaveGapIntegralCuhre::fPar; + +double TCosSqDWaveGapIntegralCuhre::IntegrateFunc() +{ + const unsigned int NCOMP(1); + const double EPSREL (1e-8); + const double EPSABS (1e-6); + const unsigned int VERBOSE (0); + const unsigned int LAST (4); + const unsigned int MINEVAL (0); + const unsigned int MAXEVAL (500000); + + const unsigned int KEY (13); + + int nregions, neval, fail; + double integral[NCOMP], error[NCOMP], prob[NCOMP]; + + Cuhre(fNDim, NCOMP, Integrand, + EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, + KEY, + &nregions, &neval, &fail, integral, error, prob); + + return integral[0]; +} + +void TCosSqDWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[], + const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, DeltaD(T), Ec, phic, DeltaS(T)} +{ + double deltasq(TMath::Power(fPar[1]*TMath::Cos(2.0*x[1]*fPar[3]) + fPar[4], 2.0)); + f[0] = TMath::Power(TMath::Cos(x[1]*fPar[3])/TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0); + return; +} + +std::vector TSinSqDWaveGapIntegralCuhre::fPar; + +double TSinSqDWaveGapIntegralCuhre::IntegrateFunc() +{ + const unsigned int NCOMP(1); + const double EPSREL (1e-8); + const double EPSABS (1e-10); + const unsigned int VERBOSE (0); + const unsigned int LAST (4); + const unsigned int MINEVAL (0); + const unsigned int MAXEVAL (500000); + + const unsigned int KEY (13); + + int nregions, neval, fail; + double integral[NCOMP], error[NCOMP], prob[NCOMP]; + + Cuhre(fNDim, NCOMP, Integrand, + EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, + KEY, + &nregions, &neval, &fail, integral, error, prob); + + return integral[0]; +} + +void TSinSqDWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[], + const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, DeltaD(T), Ec, phic, DeltaS(T)} +{ + double deltasq(TMath::Power(fPar[1]*TMath::Cos(2.0*x[1]*fPar[3]) + fPar[4],2.0)); + f[0] = TMath::Power(TMath::Sin(x[1]*fPar[3]),2.0)/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0); + return; +} + + std::vector TAnSWaveGapIntegralCuhre::fPar; double TAnSWaveGapIntegralCuhre::IntegrateFunc() diff --git a/src/external/BMWIntegrator/BMWIntegrator.h b/src/external/BMWIntegrator/BMWIntegrator.h index 804fcd2f..5b6b250c 100644 --- a/src/external/BMWIntegrator/BMWIntegrator.h +++ b/src/external/BMWIntegrator/BMWIntegrator.h @@ -143,6 +143,34 @@ class TDWaveGapIntegralCuhre { }; +class TCosSqDWaveGapIntegralCuhre { + public: + TCosSqDWaveGapIntegralCuhre() : fNDim(2) {} + ~TCosSqDWaveGapIntegralCuhre() { fPar.clear(); } + void SetParameters(const std::vector &par) { fPar=par; } + static void Integrand(const int*, const double[], const int*, double[]); + double IntegrateFunc(); + + protected: + static vector fPar; + unsigned int fNDim; + +}; + +class TSinSqDWaveGapIntegralCuhre { + public: + TSinSqDWaveGapIntegralCuhre() : fNDim(2) {} + ~TSinSqDWaveGapIntegralCuhre() { fPar.clear(); } + void SetParameters(const std::vector &par) { fPar=par; } + static void Integrand(const int*, const double[], const int*, double[]); + double IntegrateFunc(); + + protected: + static vector fPar; + unsigned int fNDim; + +}; + class TAnSWaveGapIntegralCuhre { public: TAnSWaveGapIntegralCuhre() : fNDim(2) {} diff --git a/src/external/libGapIntegrals/TGapIntegrals.cpp b/src/external/libGapIntegrals/TGapIntegrals.cpp index 092047a2..f7414956 100644 --- a/src/external/libGapIntegrals/TGapIntegrals.cpp +++ b/src/external/libGapIntegrals/TGapIntegrals.cpp @@ -43,10 +43,13 @@ using namespace std; ClassImp(TGapSWave) ClassImp(TGapDWave) +ClassImp(TGapCosSqDWave) +ClassImp(TGapSinSqDWave) ClassImp(TGapAnSWave) ClassImp(TGapNonMonDWave1) ClassImp(TGapNonMonDWave2) ClassImp(TGapPowerLaw) +ClassImp(TGapDirtySWave) ClassImp(TLambdaSWave) ClassImp(TLambdaDWave) @@ -90,6 +93,30 @@ TGapDWave::TGapDWave() { fPar.clear(); } +TGapCosSqDWave::TGapCosSqDWave() { + TCosSqDWaveGapIntegralCuhre *gapint = new TCosSqDWaveGapIntegralCuhre(); + fGapIntegral = gapint; + gapint = 0; + + fTemp.clear(); + fTempIter = fTemp.end(); + fIntegralValues.clear(); + fCalcNeeded.clear(); + fPar.clear(); +} + +TGapSinSqDWave::TGapSinSqDWave() { + TSinSqDWaveGapIntegralCuhre *gapint = new TSinSqDWaveGapIntegralCuhre(); + fGapIntegral = gapint; + gapint = 0; + + fTemp.clear(); + fTempIter = fTemp.end(); + fIntegralValues.clear(); + fCalcNeeded.clear(); + fPar.clear(); +} + TGapAnSWave::TGapAnSWave() { TAnSWaveGapIntegralCuhre *gapint = new TAnSWaveGapIntegralCuhre(); fGapIntegral = gapint; @@ -188,6 +215,28 @@ TGapDWave::~TGapDWave() { fPar.clear(); } +TGapCosSqDWave::~TGapCosSqDWave() { + delete fGapIntegral; + fGapIntegral = 0; + + fTemp.clear(); + fTempIter = fTemp.end(); + fIntegralValues.clear(); + fCalcNeeded.clear(); + fPar.clear(); +} + +TGapSinSqDWave::~TGapSinSqDWave() { + delete fGapIntegral; + fGapIntegral = 0; + + fTemp.clear(); + fTempIter = fTemp.end(); + fIntegralValues.clear(); + fCalcNeeded.clear(); + fPar.clear(); +} + TGapAnSWave::~TGapAnSWave() { delete fGapIntegral; fGapIntegral = 0; @@ -384,7 +433,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(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(4.0*(t+intPar[1])); // upper limit of energy-integration: cutoff energy intPar.push_back(TMath::PiOver2()); // upper limit of phi-integration @@ -409,6 +459,152 @@ double TGapDWave::operator()(double t, const vector &par) const { } +double TGapCosSqDWave::operator()(double t, const vector &par) const { + + assert(par.size() == 3); // three parameters: Tc, DeltaD(0), DeltaS(0) + + if (t<=0.0) + return 1.0; + else if (t >= par[0]) + return 0.0; + + bool integralParChanged(false); + + if (fPar.empty()) { // first time calling this routine + fPar = par; + integralParChanged = true; + } else { // check if Tc or Delta0 have changed + for (unsigned int i(0); i 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(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) + +// double xl[] = {0.0, 0.0}; // lower bound E, phi +// double xu[] = {4.0*(t+intPar[1]), 0.5*PI}; // upper bound E, phi + + fGapIntegral->SetParameters(intPar); +// ds = 1.0+4.0/PI*fGapIntegral->IntegrateFunc(2, xl, xu); + ds = 1.0-2.0*intPar[2]/intPar[0]*fGapIntegral->IntegrateFunc(); + + intPar.clear(); + + if (newTemp) + fIntegralValues.push_back(ds); + else + fIntegralValues[vectorIndex] = ds; + + fCalcNeeded[vectorIndex] = false; + } + + return fIntegralValues[vectorIndex]; + +} + +double TGapSinSqDWave::operator()(double t, const vector &par) const { + + assert(par.size() == 3); // two parameters: Tc, Delta(0), s-weight + + if (t<=0.0) + return 1.0; + else if (t >= par[0]) + return 0.0; + + bool integralParChanged(false); + + if (fPar.empty()) { // first time calling this routine + fPar = par; + integralParChanged = true; + } else { // check if Tc or Delta0 have changed + for (unsigned int i(0); i 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(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) + +// double xl[] = {0.0, 0.0}; // lower bound E, phi +// double xu[] = {4.0*(t+intPar[1]), 0.5*PI}; // upper bound E, phi + + fGapIntegral->SetParameters(intPar); +// ds = 1.0+4.0/PI*fGapIntegral->IntegrateFunc(2, xl, xu); + ds = 1.0-2.0*intPar[2]/intPar[0]*fGapIntegral->IntegrateFunc(); + + intPar.clear(); + + if (newTemp) + fIntegralValues.push_back(ds); + else + fIntegralValues[vectorIndex] = ds; + + fCalcNeeded[vectorIndex] = false; + } + + return fIntegralValues[vectorIndex]; + +} + double TGapAnSWave::operator()(double t, const vector &par) const { assert(par.size() == 3); // three parameters: Tc, Delta(0), a @@ -633,6 +829,22 @@ double TGapPowerLaw::operator()(double t, const vector &par) const { } +double TGapDirtySWave::operator()(double t, const vector &par) const { + + assert(par.size() == 2); // two parameters: Tc, Delta(0) + + 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 TLambdaSWave::operator()(double t, const vector &par) const { assert(par.size() == 2); // two parameters: Tc, Delta0 diff --git a/src/external/libGapIntegrals/TGapIntegrals.h b/src/external/libGapIntegrals/TGapIntegrals.h index 3c4c54e9..ae5abfe9 100644 --- a/src/external/libGapIntegrals/TGapIntegrals.h +++ b/src/external/libGapIntegrals/TGapIntegrals.h @@ -80,6 +80,47 @@ private: ClassDef(TGapDWave,1) }; +class TGapCosSqDWave : public PUserFcnBase { + +public: + TGapCosSqDWave(); + ~TGapCosSqDWave(); + + double operator()(double, const vector&) const; + +private: + TCosSqDWaveGapIntegralCuhre *fGapIntegral; + mutable vector fTemp; + mutable vector::const_iterator fTempIter; + mutable vector fIntegralValues; + mutable vector fCalcNeeded; + + mutable vector fPar; + + ClassDef(TGapCosSqDWave,1) +}; + +class TGapSinSqDWave : public PUserFcnBase { + +public: + TGapSinSqDWave(); + ~TGapSinSqDWave(); + + double operator()(double, const vector&) const; + +private: + TSinSqDWaveGapIntegralCuhre *fGapIntegral; + mutable vector fTemp; + mutable vector::const_iterator fTempIter; + mutable vector fIntegralValues; + mutable vector fCalcNeeded; + + mutable vector fPar; + + ClassDef(TGapSinSqDWave,1) +}; + + class TGapAnSWave : public PUserFcnBase { public: @@ -154,6 +195,20 @@ private: ClassDef(TGapPowerLaw,1) }; +class TGapDirtySWave : public PUserFcnBase { + +public: + TGapDirtySWave() {} + ~TGapDirtySWave() {} + + double operator()(double, const vector&) const; + +private: + + ClassDef(TGapDirtySWave,1) +}; + + class TLambdaSWave : public PUserFcnBase { public: diff --git a/src/external/libGapIntegrals/TGapIntegralsLinkDef.h b/src/external/libGapIntegrals/TGapIntegralsLinkDef.h index 005c8876..65a5c5e7 100644 --- a/src/external/libGapIntegrals/TGapIntegralsLinkDef.h +++ b/src/external/libGapIntegrals/TGapIntegralsLinkDef.h @@ -39,10 +39,13 @@ #pragma link C++ class TGapSWave+; #pragma link C++ class TGapDWave+; +#pragma link C++ class TGapCosSqDWave+; +#pragma link C++ class TGapSinSqDWave+; #pragma link C++ class TGapAnSWave+; #pragma link C++ class TGapNonMonDWave1+; #pragma link C++ class TGapNonMonDWave2+; #pragma link C++ class TGapPowerLaw+; +#pragma link C++ class TGapDirtySWave+; #pragma link C++ class TLambdaSWave+; #pragma link C++ class TLambdaDWave+; #pragma link C++ class TLambdaAnSWave+;