Added dirty limit s-wave expression for the superfluid density

src/external/BMWIntegrator/BMWIntegrator.cpp

@@ -66,6 +66,73 @@ void TDWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
   return;
 }
 
+std::vector<double> 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<double> 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<double> TAnSWaveGapIntegralCuhre::fPar;
 
 double TAnSWaveGapIntegralCuhre::IntegrateFunc()

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<double> &par) { fPar=par; }
+    static void Integrand(const int*, const double[], const int*, double[]);
+    double IntegrateFunc();
+
+  protected:
+    static vector<double> fPar;
+    unsigned int fNDim;
+
+};
+
+class TSinSqDWaveGapIntegralCuhre {
+  public:
+    TSinSqDWaveGapIntegralCuhre() : fNDim(2) {}
+    ~TSinSqDWaveGapIntegralCuhre() { fPar.clear(); }
+    void SetParameters(const std::vector<double> &par) { fPar=par; }
+    static void Integrand(const int*, const double[], const int*, double[]);
+    double IntegrateFunc();
+
+  protected:
+    static vector<double> fPar;
+    unsigned int fNDim;
+
+};
+
 class TAnSWaveGapIntegralCuhre {
   public:
     TAnSWaveGapIntegralCuhre() : fNDim(2) {}

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<double> &par) const {
     double ds;
     vector<double> 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<double> &par) const {
 
 }
 
+double TGapCosSqDWave::operator()(double t, const vector<double> &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<par.size(); i++) {
+      if (par[i] != fPar[i]) {
+        fPar[i] = par[i];
+        integralParChanged = true;
+      }
+    }
+  }
+
+  bool newTemp(false);
+  unsigned int vectorIndex;
+
+  if (integralParChanged) {
+    fCalcNeeded.clear();
+    fCalcNeeded.resize(fTemp.size(), true);
+  }
+
+  fTempIter = find(fTemp.begin(), fTemp.end(), t);
+  if(fTempIter == fTemp.end()) {
+    fTemp.push_back(t);
+    vectorIndex = fTemp.size() - 1;
+    fCalcNeeded.push_back(true);
+    newTemp = true;
+  } else {
+    vectorIndex = fTempIter - fTemp.begin();
+  }
+
+  if (fCalcNeeded[vectorIndex]) {
+
+    double ds;
+    vector<double> 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<double> &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<par.size(); i++) {
+      if (par[i] != fPar[i]) {
+        fPar[i] = par[i];
+        integralParChanged = true;
+      }
+    }
+  }
+
+  bool newTemp(false);
+  unsigned int vectorIndex;
+
+  if (integralParChanged) {
+    fCalcNeeded.clear();
+    fCalcNeeded.resize(fTemp.size(), true);
+  }
+
+  fTempIter = find(fTemp.begin(), fTemp.end(), t);
+  if(fTempIter == fTemp.end()) {
+    fTemp.push_back(t);
+    vectorIndex = fTemp.size() - 1;
+    fCalcNeeded.push_back(true);
+    newTemp = true;
+  } else {
+    vectorIndex = fTempIter - fTemp.begin();
+  }
+
+  if (fCalcNeeded[vectorIndex]) {
+
+    double ds;
+    vector<double> 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<double> &par) const {
 
   assert(par.size() == 3); // three parameters: Tc, Delta(0), a
@@ -633,6 +829,22 @@ double TGapPowerLaw::operator()(double t, const vector<double> &par) const {
 }
 
 
+double TGapDirtySWave::operator()(double t, const vector<double> &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<double> &par) const
 {
   assert(par.size() == 2); // two parameters: Tc, Delta0

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<double>&) const;
+
+private:
+  TCosSqDWaveGapIntegralCuhre *fGapIntegral;
+  mutable vector<double> fTemp;
+  mutable vector<double>::const_iterator fTempIter;
+  mutable vector<double> fIntegralValues;
+  mutable vector<bool> fCalcNeeded;
+
+  mutable vector<double> fPar;
+
+  ClassDef(TGapCosSqDWave,1)
+};
+
+class TGapSinSqDWave : public PUserFcnBase {
+
+public:
+  TGapSinSqDWave();
+  ~TGapSinSqDWave();
+
+  double operator()(double, const vector<double>&) const;
+
+private:
+  TSinSqDWaveGapIntegralCuhre *fGapIntegral;
+  mutable vector<double> fTemp;
+  mutable vector<double>::const_iterator fTempIter;
+  mutable vector<double> fIntegralValues;
+  mutable vector<bool> fCalcNeeded;
+
+  mutable vector<double> 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<double>&) const;
+
+private:
+
+  ClassDef(TGapDirtySWave,1)
+};
+
+
 class TLambdaSWave : public PUserFcnBase {
 
 public:

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+;