diff --git a/src/external/TFitPofB-lib/classes/Makefile.TFitPofB b/src/external/TFitPofB-lib/classes/Makefile.TFitPofB index 128e26f4..eec4dc90 100644 --- a/src/external/TFitPofB-lib/classes/Makefile.TFitPofB +++ b/src/external/TFitPofB-lib/classes/Makefile.TFitPofB @@ -6,16 +6,16 @@ ROOTCFLAGS = $(shell $(ROOTSYS)/bin/root-config --cflags) #--------------------------------------------------- OS = LINUX -CXX = g++ -CXXFLAGS = -g -Wall -Wno-trigraphs -fPIC +CXX = g++-4.2.4 +CXXFLAGS = -O3 -Wall -fopenmp -Wno-trigraphs -fPIC MUSRFITINCLUDE = ../../../include #MUSRFITINCLUDE = /home/l_wojek/rep/analysis/musrfit/src/include LOCALINCLUDE = ../include ROOTINCLUDE = $(ROOTSYS)/include INCLUDES = -I$(LOCALINCLUDE) -I$(MUSRFITINCLUDE) -I$(ROOTINCLUDE) -LD = g++ -LDFLAGS = -g -SOFLAGS = -O -shared +LD = g++-4.2.4 +LDFLAGS = -O3 +SOFLAGS = -O3 -shared -fopenmp -lfftw3_threads -lfftw3 -lm -lpthread # the output from the root-config script: CXXFLAGS += $(ROOTCFLAGS) @@ -29,6 +29,7 @@ OBJS += TPofBCalc.o OBJS += TPofTCalc.o OBJS += TFitPofBStartupHandler.o TFitPofBStartupHandlerDict.o OBJS += TLondon1D.o TLondon1DDict.o +OBJS += TSkewedGss.o TSkewedGssDict.o SHLIB = libTFitPofB.so @@ -59,10 +60,14 @@ TLondon1DDict.cpp: ../include/TLondon1D.h ../include/TLondon1DLinkDef.h @echo "Generating dictionary $@..." rootcint -f $@ -c -p -I$(MUSRFITINCLUDE) $^ +TSkewedGssDict.cpp: ../include/TSkewedGss.h ../include/TSkewedGssLinkDef.h + @echo "Generating dictionary $@..." + rootcint -f $@ -c -p -I$(MUSRFITINCLUDE) $^ + TFitPofBStartupHandlerDict.cpp: ../include/TFitPofBStartupHandler.h ../include/TFitPofBStartupHandlerLinkDef.h @echo "Generating dictionary $@..." rootcint -f $@ -c -p $^ - + install: all @echo "Installing shared lib: libTFitPofB.so" ifeq ($(OS),LINUX) diff --git a/src/external/TFitPofB-lib/classes/TBofZCalc.cpp b/src/external/TFitPofB-lib/classes/TBofZCalc.cpp index 449caa45..5b70b3a2 100644 --- a/src/external/TFitPofB-lib/classes/TBofZCalc.cpp +++ b/src/external/TFitPofB-lib/classes/TBofZCalc.cpp @@ -5,26 +5,38 @@ Author: Bastian M. Wojek e-mail: bastian.wojek@psi.ch - 2008/06/30 + 2009/04/25 ***************************************************************************/ #include "TBofZCalc.h" #include #include -#include -#include +//#include +//#include +#include -vector DataZ() const{ - if (fZ.empty) - Calculate(); - return fZ; -} +//------------------ +// Method filling the z and B(z)-vectors for the case B(z) should be used numerically +//------------------ -vector DataBZ() const{ - if (fBZ.empty) - Calculate(); - return fBZ; +void TBofZCalc::Calculate() +{ + if (!fBZ.empty()) + return; + + double ZZ; + int j; + fZ.resize(fSteps); + fBZ.resize(fSteps); + +#pragma omp parallel for default(shared) private(j,ZZ) schedule(dynamic) + for (j=0; j DataBZ() const{ // Parameters: Bext[G], deadlayer[nm], lambda[nm] //------------------ -TLondon1D_HS::TLondon1D_HS(unsigned int steps, const vector ¶m) { - +TLondon1D_HS::TLondon1D_HS(const vector ¶m, unsigned int steps) +{ + fSteps = steps; fDZ = 200.0/double(steps); - double ZZ, BBz; - - for (unsigned int j(0); j > TLondon1D_HS::GetInverseAndDerivative(double BB) const +{ + vector< pair > inv; + + if(BB <= 0.0 || BB > fParam[0]) + return inv; + + pair invAndDerivative; + + invAndDerivative.first = fParam[1] - fParam[2]*log(BB/fParam[0]); + invAndDerivative.second = -fParam[2]/BB; + + inv.push_back(invAndDerivative); + + return inv; +} + + //------------------ // Constructor of the TLondon1D_1L class // 1D-London screening in a thin superconducting film // Parameters: Bext[G], deadlayer[nm], thickness[nm], lambda[nm] //------------------ -TLondon1D_1L::TLondon1D_1L(unsigned int steps, const vector ¶m) { - - double N(cosh(param[2]/2.0/param[3])); - +TLondon1D_1L::TLondon1D_1L(const vector ¶m, unsigned int steps) +{ + fSteps = steps; fDZ = param[2]/double(steps); - double ZZ, BBz; + fParam = param; + fMinZ = -1.0; + fMinB = -1.0; - for (unsigned int j(0); j=0.); + +// lambdas have to be greater than zero + assert(param[3]!=0.); + +// Calculate the coefficients of the exponentials + double N0(param[0]/(1.0+exp(param[2]/param[3]))); + + fCoeff[0]=N0*exp((param[1]+param[2])/param[3]); + fCoeff[1]=N0*exp(-param[1]/param[3]); + +// none of the coefficients should be zero + for(unsigned int i(0); i<2; i++) + assert(fCoeff[i]); + + SetBmin(); +} + +double TLondon1D_1L::GetBofZ(double ZZ) const +{ + if(ZZ < 0. || ZZ < fParam[1] || ZZ > fParam[1]+fParam[2]) + return fParam[0]; + + return fCoeff[0]*exp(-ZZ/fParam[3])+fCoeff[1]*exp(ZZ/fParam[3]); +} + +double TLondon1D_1L::GetBmax() const +{ + // return applied field + return fParam[0]; +} + +double TLondon1D_1L::GetBmin() const +{ + // return field minimum + return fMinB; +} + +void TLondon1D_1L::SetBmin() +{ + double b_a(fCoeff[1]/fCoeff[0]); + assert (b_a>0.); + + double minZ; + // check if the minimum is in the first layer + minZ=-0.5*fParam[3]*log(b_a); + if (minZ > fParam[1] && minZ <= fParam[2]) { + fMinZ = minZ; + fMinB = GetBofZ(minZ); + return; } + assert(fMinZ > 0. && fMinB > 0.); + return; } +vector< pair > TLondon1D_1L::GetInverseAndDerivative(double BB) const +{ + vector< pair > inv; + + if(BB <= fMinB || BB > fParam[0]) + return inv; + + double inverse[2]; + pair invAndDerivative; + + inverse[0]=fParam[3]*log((BB-sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]))/(2.0*fCoeff[1])); + inverse[1]=fParam[3]*log((BB+sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]))/(2.0*fCoeff[1])); + + if(inverse[0] > fParam[1] && inverse[0] < fMinZ) { + invAndDerivative.first = inverse[0]; + invAndDerivative.second = -fParam[3]/sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]); + inv.push_back(invAndDerivative); + } + if(inverse[1] > fMinZ && inverse[1] <= fParam[1]+fParam[2]) { + invAndDerivative.first = inverse[1]; + invAndDerivative.second = +fParam[3]/sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]); + inv.push_back(invAndDerivative); + } + + return inv; +} + + //------------------ // Constructor of the TLondon1D_2L class // 1D-London screening in a thin superconducting film, two layers, two lambda // Parameters: Bext[G], deadlayer[nm], thickness1[nm], thickness2[nm], lambda1[nm], lambda2[nm] //------------------ -TLondon1D_2L::TLondon1D_2L(unsigned int steps, const vector ¶m) { - - double N1(param[5]*cosh(param[3]/param[5])*sinh(param[2]/param[4]) + param[4]*cosh(param[2]/param[4])*sinh(param[3]/param[5])); - double N2(4.0*N1); - +TLondon1D_2L::TLondon1D_2L(const vector ¶m, unsigned int steps) +{ + fSteps = steps; fDZ = (param[2]+param[3])/double(steps); - double ZZ, BBz; + fParam = param; + fMinTag = -1; + fMinZ = -1.0; + fMinB = -1.0; - for (unsigned int j(0); j=0.); + +// lambdas have to be greater than zero + for(unsigned int i(4); i<6; i++){ + assert(param[i]!=0.); + } + fInterfaces[0]=param[1]; + fInterfaces[1]=param[1]+param[2]; + fInterfaces[2]=param[1]+param[2]+param[3]; + +// Calculate the coefficients of the exponentials + double B0(param[0]), dd(param[1]), d1(param[2]), d2(param[3]), l1(param[4]), l2(param[5]); + + double N0((1.0+exp(2.0*d1/l1))*(-1.0+exp(2.0*d2/l2))*l1+(-1.0+exp(2.0*d1/l1))*(1.0+exp(2.0*d2/l2))*l2); + + double N1(4.0*(l2*cosh(d2/l2)*sinh(d1/l1)+l1*cosh(d1/l1)*sinh(d2/l2))); + + fCoeff[0]=B0*exp((d1+dd)/l1)*(-2.0*exp(d2/l2)*l2+exp(d1/l1)*(l2-l1)+exp(d1/l1+2.0*d2/l2)*(l1+l2))/N0; + + fCoeff[1]=-B0*exp(-(dd+d1)/l1-d2/l2)*(l1-exp(2.0*d2/l2)*l1+(1.0-2.0*exp(d1/l1+d2/l2)+exp(2.0*d2/l2))*l2)/N1; + + fCoeff[2]=B0*exp((d1+d2+dd)/l2)*(-(1.0+exp(2.0*d1/l1)-2.0*exp(d1/l1+d2/l2))*l1+(-1.0+exp(2.0*d1/l1))*l2)/N0; + + fCoeff[3]=B0*exp(-d1/l1-(d1+d2+dd)/l2)*(-2.0*exp(d1/l1)*l1+exp(d2/l2)*(l1-l2+exp(2.0*d1/l1)*(l1+l2)))/N1; + +// none of the coefficients should be zero + for(unsigned int i(0); i<4; i++) + assert(fCoeff[i]); + + SetBmin(); +} + +double TLondon1D_2L::GetBofZ(double ZZ) const +{ + if(ZZ < 0. || ZZ < fInterfaces[0] || ZZ > fInterfaces[2]) + return fParam[0]; + + if (ZZ < fInterfaces[1]) { + return fCoeff[0]*exp(-ZZ/fParam[4])+fCoeff[1]*exp(ZZ/fParam[4]); + } else { + return fCoeff[2]*exp(-ZZ/fParam[5])+fCoeff[3]*exp(ZZ/fParam[5]); + } +} + +double TLondon1D_2L::GetBmax() const +{ + // return applied field + return fParam[0]; +} + +double TLondon1D_2L::GetBmin() const +{ + // return field minimum + return fMinB; +} + +void TLondon1D_2L::SetBmin() +{ + double b_a(fCoeff[1]/fCoeff[0]); + assert (b_a>0.); + + double minZ; + // check if the minimum is in the first layer + minZ=-0.5*fParam[4]*log(b_a); + if (minZ > fInterfaces[0] && minZ <= fInterfaces[1]) { + fMinTag = 1; + fMinZ = minZ; + fMinB = GetBofZ(minZ); + return; } + double d_c(fCoeff[3]/fCoeff[2]); + assert (d_c>0.); + + // check if the minimum is in the second layer + minZ=-0.5*fParam[5]*log(d_c); + if (minZ > fInterfaces[1] && minZ <= fInterfaces[2]) { + fMinTag = 2; + fMinZ = minZ; + fMinB = GetBofZ(minZ); + return; + } + + assert(fMinZ > 0. && fMinB > 0.); + return; +} + +vector< pair > TLondon1D_2L::GetInverseAndDerivative(double BB) const +{ + vector< pair > inv; + + if(BB <= fMinB || BB > fParam[0]) + return inv; + + double inverse[3]; + pair invAndDerivative; + + switch(fMinTag) + { + case 1: + inverse[0]=fParam[4]*log((BB-sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]))/(2.0*fCoeff[1])); + inverse[1]=fParam[4]*log((BB+sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]))/(2.0*fCoeff[1])); + inverse[2]=fParam[5]*log((BB+sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]))/(2.0*fCoeff[3])); + + if(inverse[0] > fInterfaces[0] && inverse[0] < fMinZ) { + invAndDerivative.first = inverse[0]; + invAndDerivative.second = -fParam[4]/sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]); + inv.push_back(invAndDerivative); + } + if(inverse[1] > fMinZ && inverse[1] <= fInterfaces[1]) { + invAndDerivative.first = inverse[1]; + invAndDerivative.second = +fParam[4]/sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]); + inv.push_back(invAndDerivative); + } + if(inverse[2] > fInterfaces[1] && inverse[2] <= fInterfaces[2]) { + invAndDerivative.first = inverse[2]; + invAndDerivative.second = +fParam[5]/sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]); + inv.push_back(invAndDerivative); + } + break; + + case 2: + inverse[0]=fParam[4]*log((BB-sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]))/(2.0*fCoeff[1])); + inverse[1]=fParam[5]*log((BB-sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]))/(2.0*fCoeff[3])); + inverse[2]=fParam[5]*log((BB+sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]))/(2.0*fCoeff[3])); + + if(inverse[0] > fInterfaces[0] && inverse[0] <= fInterfaces[1]) { + invAndDerivative.first = inverse[0]; + invAndDerivative.second = -fParam[4]/sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]); + inv.push_back(invAndDerivative); + } + if(inverse[1] > fInterfaces[1] && inverse[1] < fMinZ) { + invAndDerivative.first = inverse[1]; + invAndDerivative.second = -fParam[5]/sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]); + inv.push_back(invAndDerivative); + } + if(inverse[2] > fMinZ && inverse[2] <= fInterfaces[2]) { + invAndDerivative.first = inverse[2]; + invAndDerivative.second = +fParam[5]/sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]); + inv.push_back(invAndDerivative); + } + break; + + default: + break; + } + return inv; } //------------------ @@ -149,67 +409,62 @@ TLondon1D_2L::TLondon1D_2L(unsigned int steps, const vector ¶m) { // Parameters: Bext[G], deadlayer[nm], thickness1[nm], thickness2[nm], thickness3[nm], lambda1[nm], lambda2[nm], lambda3[nm] //------------------ -TLondon1D_3L::TLondon1D_3L(unsigned int steps, const vector ¶m) - : fSteps(steps), fDZ((param[2]+param[3]+param[4])/double(steps)), fParam(param), fMinTag(-1), fMinZ(-1.0), fMinB(-1.0) +TLondon1D_3L::TLondon1D_3L(const vector ¶m, unsigned int steps) +// : fSteps(steps), fDZ((param[2]+param[3]+param[4])/double(steps)), fParam(param), fMinTag(-1), fMinZ(-1.0), fMinB(-1.0) { +// no members of TLondon1D_3L, therefore the initialization list cannot be used!! + fSteps = steps; + fDZ = (param[2]+param[3]+param[4])/double(steps); + fParam = param; + fMinTag = -1; + fMinZ = -1.0; + fMinB = -1.0; + // thicknesses have to be greater or equal to zero for(unsigned int i(1); i<5; i++) assert(param[i]>=0.); // lambdas have to be greater than zero - for(unsigned int i(5); i0.); + for(unsigned int i(5); i<8; i++) + assert(param[i]!=0.); - fInterfaces={param[1], param[1]+param[2], param[1]+param[2]+param[3], param[1]+param[2]+param[3]+param[4]}; + fInterfaces[0]=param[1]; + fInterfaces[1]=param[1]+param[2]; + fInterfaces[2]=param[1]+param[2]+param[3]; + fInterfaces[3]=param[1]+param[2]+param[3]+param[4]; // Calculate the coefficients of the exponentials - double N0(4.0*((1.0+exp(2.0*param[2]/param[5]))*param[5]*(param[7]*cosh(param[4]/param[7])*sinh(param[3]/param[6]) + \ - param[6]*cosh(param[3]/param[6])*sinh(param[4]/param[7]))+(-1.0+exp(2.0*param[2]/param[5])) * \ - param[6]*(param[7]*cosh(param[3]/param[6])*cosh(param[4]/param[7])+param[6]*sinh(param[3]/param[6])*sinh(param[4]/param[7])))); + double B0(param[0]), dd(param[1]), d1(param[2]), d2(param[3]), d3(param[4]), l1(param[5]), l2(param[6]), l3(param[7]); - double N2(2.0*param[7]*cosh(param[4]/param[7])*((-1.0+exp(2.0*param[2]/param[5]))*param[6]*cosh(param[3]/param[6]) + \ - (1.0+exp(2.0*param[2]/param[5]))*param[5]*sinh(param[3]/param[6])) + \ - 4.0*exp(param[2]/param[5])*param[6]*(param[5]*cosh(param[2]/param[5])*cosh(param[3]/param[6]) + \ - param[6]*sinh(param[2]/param[5])*sinh(param[3]/param[6]))*sinh(param[4]/param[7])); + double N0(4.0*((1.0+exp(2.0*d1/l1))*l1*(l3*cosh(d3/l3)*sinh(d2/l2) + l2*cosh(d2/l2)*sinh(d3/l3))+(-1.0+exp(2.0*d1/l1)) * \ + l2*(l3*cosh(d2/l2)*cosh(d3/l3)+l2*sinh(d2/l2)*sinh(d3/l3)))); - double N3(2.0*((1.0+exp(2.0*param[2]/param[5]))*param[5]*(param[7]*cosh(param[4]/param[7])*sinh(param[3]/param[6]) + \ - param[6]*cosh(param[3]/param[6])*sinh(param[4]/param[7])) + \ - (-1.0+exp(2.0*param[2]/param[5]))*param[6]*(param[7]*cosh(param[3]/param[6]) * \ - cosh(param[4]/param[7])+param[6]*sinh(param[3]/param[6])*sinh(param[4]/param[7])))); + double N2(2.0*l3*cosh(d3/l3)*((-1.0+exp(2.0*d1/l1))*l2*cosh(d2/l2) + (1.0+exp(2.0*d1/l1))*l1*sinh(d2/l2)) + \ + 4.0*exp(d1/l1)*l2*(l1*cosh(d1/l1)*cosh(d2/l2) + l2*sinh(d1/l1)*sinh(d2/l2))*sinh(d3/l3)); - double N5((1.0+exp(2.0*param[2]/param[5]))*param[5]*(param[7]*cosh(param[4]/param[7])*sinh(param[3]/param[6]) + \ - param[6]*cosh(param[3]/param[6])*sinh(param[4]/param[7])) + \ - (-1.0+exp(2.0*param[2]/param[5]))*param[6]*(param[7]*cosh(param[3]/param[6]) * \ - cosh(param[4]/param[7])+param[6]*sinh(param[3]/param[6])*sinh(param[4]/param[7]))); + double N3(2.0*((1.0+exp(2.0*d1/l1))*l1*(l3*cosh(d3/l3)*sinh(d2/l2) + l2*cosh(d2/l2)*sinh(d3/l3)) + \ + (-1.0+exp(2.0*d1/l1))*l2*(l3*cosh(d2/l2) * cosh(d3/l3)+l2*sinh(d2/l2)*sinh(d3/l3)))); - fCoeff[0]=(B0*exp((param[2]+param[1])/param[5]-param[3]/param[6]-param[4]/param[7]) * \ - (exp(param[2]/param[5])*(param[5]-param[6])*(-param[6]+exp(2.0*param[4]/param[7])*(param[6]-param[7])-param[7]) - \ - 4.0*exp(param[3]/param[6]+param[4]/param[7]) * param[6]*param[7] + \ - exp(param[2]/param[5]+2.0*param[3]/param[6])*(param[5]+param[6]) * \ - (-param[6]+param[7]+exp(2.0*param[4]/param[7])*(param[6]+param[7]))))/N0; + double N5((1.0+exp(2.0*d1/l1))*l1*(l3*cosh(d3/l3)*sinh(d2/l2) + l2*cosh(d2/l2)*sinh(d3/l3)) + \ + (-1.0+exp(2.0*d1/l1))*l2*(l3*cosh(d2/l2) * cosh(d3/l3)+l2*sinh(d2/l2)*sinh(d3/l3))); - fCoeff[1]=(B0*exp(-param[1]/param[5]-param[3]/param[6]-param[4]/param[7])*(-param[6]*((-1.0+exp(2.0*param[3]/param[6])) * \ - (-1.0+exp(2.0*param[4]/param[7]))*param[6]+(1.0+exp(2.0*param[3]/param[6]) - \ - 4.0*exp(param[2]/param[5]+param[3]/param[6]+param[4]/param[7]) + exp(2.0*param[4]/param[7])*(1.0+exp(2.0*param[3]/param[6])))*param[7]) + \ - 4.0*exp(param[3]/param[6]+param[4]/param[7])*param[5]*(param[7]*cosh(param[4]/param[7]) * \ - sinh(param[3]/param[6])+param[6]*cosh(param[3]/param[6])*sinh(param[4]/param[7]))))/N0; + fCoeff[0]=(B0*exp((d1+dd)/l1-d2/l2-d3/l3) * (exp(d1/l1)*(l1-l2)*(-l2+exp(2.0*d3/l3)*(l2-l3)-l3) - \ + 4.0*exp(d2/l2+d3/l3) * l2*l3 + exp(d1/l1+2.0*d2/l2)*(l1+l2) * (-l2+l3+exp(2.0*d3/l3)*(l2+l3))))/N0; - fCoeff[2]=(B0*exp((param[2]+param[1])/param[6])*(-exp(2.0*param[2]/param[5])*(param[5]-param[6])*param[7]-(param[5]+param[6])*param[7] + \ - 2.0*exp(param[2]/param[5]+param[3]/param[6])*param[5]*(param[7]*cosh(param[4]/param[7])+param[6]*sinh(param[4]/param[7]))))/N2; + fCoeff[1]=(B0*exp(-dd/l1-d2/l2-d3/l3)*(-l2*((-1.0+exp(2.0*d2/l2)) * (-1.0+exp(2.0*d3/l3))*l2+(1.0+exp(2.0*d2/l2) - \ + 4.0*exp(d1/l1+d2/l2+d3/l3) + exp(2.0*d3/l3)*(1.0+exp(2.0*d2/l2)))*l3) + 4.0*exp(d2/l2+d3/l3)*l1*(l3*cosh(d3/l3) * \ + sinh(d2/l2)+l2*cosh(d2/l2)*sinh(d3/l3))))/N0; - fCoeff[3]=(B0*exp(-(param[2]+param[3]+param[1])/param[6]-param[4]/param[7])*(exp(param[3]/param[6]+param[4]/param[7]) * \ - (param[5]-param[6]+exp(2.0*param[2]/param[5])*(param[5]+param[6]))*param[7] - \ - exp(param[2]/param[5])*param[5]*(param[6]+param[7]+exp(2.0*param[4]/param[7])*(param[7]-param[6]))))/N3; + fCoeff[2]=(B0*exp((d1+dd)/l2)*(-exp(2.0*d1/l1)*(l1-l2)*l3-(l1+l2)*l3 + 2.0*exp(d1/l1+d2/l2)*l1*(l3*cosh(d3/l3)+l2*sinh(d3/l3))))/N2; - fCoeff[4]=(0.25*B0*exp(-param[3]/param[6]+(param[2]+param[3]+param[1])/param[7]) * \ - (4.0*exp(param[2]/param[5]+param[3]/param[6]+param[4]/param[7])*param[5]*param[6] - \ - (-1.0+exp(2.0*param[2]/param[5]))*param[6]*(-param[6]+exp(2.0*param[3]/param[6])*(param[6]-param[7])-param[7]) - \ - (1.0+exp(2.0*param[2]/param[5])*param[5]*(param[6]+exp(2.0*param[3]/param[6])*(param[6]-param[7])+param[7])))/N5; + fCoeff[3]=(B0*exp(-(d1+d2+dd)/l2-d3/l3)*(exp(d2/l2+d3/l3) * (l1-l2+exp(2.0*d1/l1)*(l1+l2))*l3 - \ + exp(d1/l1)*l1*(l2+l3+exp(2.0*d3/l3)*(l3-l2))))/N3; - fCoeff[5]=(B0*exp(-(param[2]+param[3]+param[4]+param[1])/param[7]+param[2]/param[5]) * \ - (-param[5]*param[6]+exp(param[4]/param[7])*(param[6]*sinh(param[2]/param[5])*(param[7]*cosh(param[3]/param[6]) + \ - param[6]*sinh(param[3]/param[6])) + param[5]*cosh(param[2]/param[5]) * \ - (param[6]*cosh(param[3]/param[6])+param[7]*sinh(param[3]/param[6])))))/N5; + fCoeff[4]=(0.25*B0*exp(-d2/l2+(d1+d2+dd)/l3) * (4.0*exp(d1/l1+d2/l2+d3/l3)*l1*l2 - \ + (-1.0+exp(2.0*d1/l1))*l2*(-l2+exp(2.0*d2/l2)*(l2-l3)-l3) - (1.0+exp(2.0*d1/l1))*l1*(l2+exp(2.0*d2/l2)*(l2-l3)+l3)))/N5; + + fCoeff[5]=(B0*exp(-(d1+d2+d3+dd)/l3+d1/l1) * (-l1*l2+exp(d3/l3)*(l2*sinh(d1/l1)*(l3*cosh(d2/l2) + \ + l2*sinh(d2/l2)) + l1*cosh(d1/l1) * (l2*cosh(d2/l2)+l3*sinh(d2/l2)))))/N5; // none of the coefficients should be zero for(unsigned int i(0); i<6; i++) @@ -218,70 +473,17 @@ TLondon1D_3L::TLondon1D_3L(unsigned int steps, const vector ¶m) SetBmin(); } -void TLondon1D_3L::Calculate() -{ - if (!fBZ.empty()) - return; - - double ZZ; - int j; - fZ.resize(fSteps); - fBZ.resize(fSteps); - -#pragma omp parallel for default(shared) private(j,ZZ) schedule(dynamic) - for (j=0; j fInterfaces[3]) return fParam[0]; - double N1(fParam[7]*cosh(fParam[4]/fParam[7])*((exp(2.0*fParam[2]/fParam[5])-1.0)*fParam[6]*cosh(fParam[3]/fParam[6]) + \ - (1.0+exp(2.0*fParam[2]/fParam[5]))*fParam[5]*sinh(fParam[3]/fParam[6])) + \ - 2.0*exp(fParam[2]/fParam[5])*fParam[6]*(fParam[5]*cosh(fParam[2]/fParam[5])*cosh(fParam[3]/fParam[6]) + \ - fParam[6]*sinh(fParam[2]/fParam[5])*sinh(fParam[3]/fParam[6]))*sinh(fParam[4]/fParam[7])); - - double N21(2.0*(fParam[7]*cosh(fParam[4]/fParam[7])*((-1.0+exp(2.0*fParam[2]/fParam[5]))*fParam[6]*cosh(fParam[3]/fParam[6]) + \ - (1.0+exp(2.0*fParam[2]/fParam[5]))*fParam[5]*sinh(fParam[3]/fParam[6])) + \ - 2.0*exp(fParam[2]/fParam[5])*fParam[6]*(fParam[5]*cosh(fParam[2]/fParam[5])*cosh(fParam[3]/fParam[6]) + \ - fParam[6]*sinh(fParam[2]/fParam[5])*sinh(fParam[3]/fParam[6]))*sinh(fParam[4]/fParam[7]))); - - double N22(((fParam[5]+exp(2.0*fParam[2]/fParam[5])*(fParam[5]-fParam[6])+fParam[6])*(-fParam[7]*cosh(fParam[4]/fParam[7]) + \ - fParam[6]*sinh(fParam[4]/fParam[7]))+exp(2.0*fParam[3]/fParam[6])*(fParam[5]-fParam[6] + \ - exp(2.0*fParam[2]/fParam[5])*(fParam[5]+fParam[6]))*(fParam[7]*cosh(fParam[4]/fParam[7])+fParam[6]*sinh(fParam[4]/fParam[7])))); - - double N3(4.0*((1.0+exp(2.0*fParam[2]/fParam[5]))*fParam[5]*(fParam[7]*cosh(fParam[4]/fParam[7])*sinh(fParam[3]/fParam[6]) + \ - fParam[6]*cosh(fParam[3]/fParam[6])*sinh(fParam[4]/fParam[7])) + \ - (-1.0+exp(2.0*fParam[2]/fParam[5]))*fParam[6]*(fParam[7]*cosh(fParam[3]/fParam[6])*cosh(fParam[4]/fParam[7]) + \ - fParam[6]*sinh(fParam[3]/fParam[6])*sinh(fParam[4]/fParam[7])))); - if (ZZ < fInterfaces[1]) { - return (2.0*fParam[0]*exp(fParam[2]/fParam[5])*(-fParam[6]*fParam[7]*sinh((fParam[1]-ZZ)/fParam[5]) + \ - fParam[7]*cosh(fParam[4]/fParam[7])*(fParam[6]*cosh(fParam[3]/fParam[6])*sinh((fParam[2]+fParam[1]-ZZ)/fParam[5]) + \ - fParam[5]*cosh((fParam[2]+fParam[1]-ZZ)/fParam[5])*sinh(fParam[3]/fParam[6])) + \ - fParam[6]*(fParam[5]*cosh((fParam[2]+fParam[1]-ZZ)/fParam[5])*cosh(fParam[3]/fParam[6]) + \ - fParam[6]*sinh((fParam[2]+fParam[1]-ZZ)/fParam[5])*sinh(fParam[3]/fParam[6]))*sinh(fParam[4]/fParam[7])))/N1; + return fCoeff[0]*exp(-ZZ/fParam[5])+fCoeff[1]*exp(ZZ/fParam[5]); } else if (ZZ < fInterfaces[2]) { - return (fParam[0]*exp((fParam[2]+fParam[1]-ZZ)/fParam[6])*(-fParam[5]-fParam[6]+exp(2.0*fParam[2]/fParam[5])*(fParam[6]-fParam[5]))*fParam[7])/N21 + \ - (fParam[0]*exp(-(fParam[2]+fParam[1]+ZZ)/fParam[6])*(exp((fParam[3]+2.0*ZZ)/fParam[6])*(fParam[5]-fParam[6] + \ - exp(2.0*fParam[2]/fParam[5])*(fParam[5]+fParam[6]))*fParam[7] + \ - 2.0*exp(fParam[2]/fParam[5])*fParam[5]*(exp(2.0*ZZ/fParam[6])*(-fParam[7]*cosh(fParam[4]/fParam[7]) + \ - fParam[6]*sinh(fParam[4]/fParam[7]))+(cosh(2.0*(fParam[2]+fParam[3]+fParam[1])/fParam[6]) + \ - sinh(2.0*(fParam[2]+fParam[3]+fParam[1])/fParam[6]))*(fParam[7]*cosh(fParam[4]/fParam[7])+fParam[6]*sinh(fParam[4]/fParam[7])))))/N22; + return fCoeff[2]*exp(-ZZ/fParam[6])+fCoeff[3]*exp(ZZ/fParam[6]); } else { - return (fParam[0]*exp(-(fParam[2]+fParam[3]+fParam[1]+ZZ)/fParam[7])*(-exp(-fParam[3]/fParam[6] + \ - 2.0*(fParam[2]+fParam[3]+fParam[1])/fParam[7])*(-4.0*exp(fParam[2]/fParam[5]+fParam[3]/fParam[6]+fParam[4]/fParam[7])*fParam[5]*fParam[6] + \ - (-1.0+exp(2.0*fParam[2]/fParam[5]))*fParam[6]*(-fParam[6]+exp(2.0*fParam[3]/fParam[6])*(fParam[6]-fParam[7])-fParam[7]) + \ - (1.0+exp(2.0*fParam[2]/fParam[5]))*fParam[5]*(fParam[6]+exp(2.0*fParam[3]/fParam[6])*(fParam[6]-fParam[7])+fParam[7])) + \ - 4.0*exp(fParam[2]/fParam[5]-(fParam[4]-2.0*ZZ)/fParam[7])*(-fParam[5]*fParam[6]+exp(fParam[4]/fParam[7]) * \ - (fParam[6]*sinh(fParam[2]/fParam[5])*(fParam[7]*cosh(fParam[3]/fParam[6])+fParam[6]*sinh(fParam[3]/fParam[6])) + \ - fParam[5]*cosh(fParam[2]/fParam[5])*(fParam[6]*cosh(fParam[3]/fParam[6])+fParam[7]*sinh(fParam[3]/fParam[6]))))))/N3; + return fCoeff[4]*exp(-ZZ/fParam[7])+fCoeff[5]*exp(ZZ/fParam[7]); } } @@ -336,7 +538,7 @@ void TLondon1D_3L::SetBmin() return; } - assert(fminZ > 0. && fminB > 0.); + assert(fMinZ > 0. && fMinB > 0.); return; } @@ -344,7 +546,7 @@ vector< pair > TLondon1D_3L::GetInverseAndDerivative(double BB) { vector< pair > inv; - if(BB <= fminB || BB > fParam[0]) + if(BB <= fMinB || BB > fParam[0]) return inv; double inverse[4]; @@ -358,12 +560,12 @@ vector< pair > TLondon1D_3L::GetInverseAndDerivative(double BB) inverse[2]=fParam[6]*log((BB+sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]))/(2.0*fCoeff[3])); inverse[3]=fParam[7]*log((BB+sqrt(BB*BB-4.0*fCoeff[4]*fCoeff[5]))/(2.0*fCoeff[5])); - if(inverse[0] > fInterfaces[0] && inverse[0] < fminZ) { + if(inverse[0] > fInterfaces[0] && inverse[0] < fMinZ) { invAndDerivative.first = inverse[0]; invAndDerivative.second = -fParam[5]/sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]); inv.push_back(invAndDerivative); } - if(inverse[1] > fminZ && inverse[1] <= fInterfaces[1]) { + if(inverse[1] > fMinZ && inverse[1] <= fInterfaces[1]) { invAndDerivative.first = inverse[1]; invAndDerivative.second = +fParam[5]/sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]); inv.push_back(invAndDerivative); @@ -391,12 +593,12 @@ vector< pair > TLondon1D_3L::GetInverseAndDerivative(double BB) invAndDerivative.second = -fParam[5]/sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]); inv.push_back(invAndDerivative); } - if(inverse[1] > fInterfaces[1] && inverse[1] < fminZ) { + if(inverse[1] > fInterfaces[1] && inverse[1] < fMinZ) { invAndDerivative.first = inverse[1]; invAndDerivative.second = -fParam[6]/sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]); inv.push_back(invAndDerivative); } - if(inverse[2] > fminZ && inverse[2] <= fInterfaces[2]) { + if(inverse[2] > fMinZ && inverse[2] <= fInterfaces[2]) { invAndDerivative.first = inverse[2]; invAndDerivative.second = +fParam[6]/sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]); inv.push_back(invAndDerivative); @@ -425,12 +627,12 @@ vector< pair > TLondon1D_3L::GetInverseAndDerivative(double BB) invAndDerivative.second = -fParam[6]/sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]); inv.push_back(invAndDerivative); } - if(inverse[2] > fInterfaces[2] && inverse[2] < fminZ) { + if(inverse[2] > fInterfaces[2] && inverse[2] < fMinZ) { invAndDerivative.first = inverse[2]; invAndDerivative.second = -fParam[7]/sqrt(BB*BB-4.0*fCoeff[4]*fCoeff[5]); inv.push_back(invAndDerivative); } - if(inverse[3] > fminZ && inverse[3] <= fInterfaces[3]) { + if(inverse[3] > fMinZ && inverse[3] <= fInterfaces[3]) { invAndDerivative.first = inverse[3]; invAndDerivative.second = +fParam[7]/sqrt(BB*BB-4.0*fCoeff[4]*fCoeff[5]); inv.push_back(invAndDerivative); @@ -450,74 +652,274 @@ vector< pair > TLondon1D_3L::GetInverseAndDerivative(double BB) // Parameters: Bext[G], deadlayer[nm], thickness1[nm], thickness2[nm], thickness3[nm], lambda1[nm], lambda2[nm] //------------------ -TLondon1D_3LS::TLondon1D_3LS(unsigned int steps, const vector ¶m) { - - double N1(8.0*(param[5]*param[6]*cosh(param[3]/param[6])*sinh((param[2]+param[4])/param[5]) + ((param[5]*param[5]*cosh(param[2]/param[5])*cosh(param[4]/param[5])) + (param[6]*param[6]*sinh(param[2]/param[5])*sinh(param[4]/param[5])))*sinh(param[3]/param[6]))); - - double N2(2.0*param[5]*param[6]*cosh(param[3]/param[6])*sinh((param[2]+param[4])/param[5]) + 2.0*(param[5]*param[5]*cosh(param[2]/param[5])*cosh(param[4]/param[5]) + param[6]*param[6]*sinh(param[2]/param[5])*sinh(param[4]/param[5]))*sinh(param[3]/param[6])); - - double N3(8.0*(param[5]*param[6]*cosh(param[3]/param[6])*sinh((param[2]+param[4])/param[5]) + (param[5]*param[5]*cosh(param[2]/param[5])*cosh(param[4]/param[5]) + param[6]*param[6]*sinh(param[2]/param[5])*sinh(param[4]/param[5]))*sinh(param[3]/param[6]))); - +TLondon1D_3LS::TLondon1D_3LS(const vector ¶m, unsigned int steps) +{ + fSteps = steps; fDZ = (param[2]+param[3]+param[4])/double(steps); - double ZZ, BBz; + fParam = param; + fMinTag = -1; + fMinZ = -1.0; + fMinB = -1.0; - for (unsigned int j(0); j=0.); + +// lambdas have to be greater than zero + for(unsigned int i(5); i<7; i++){ + assert(param[i]!=0.); + } + fInterfaces[0]=param[1]; + fInterfaces[1]=param[1]+param[2]; + fInterfaces[2]=param[1]+param[2]+param[3]; + fInterfaces[3]=param[1]+param[2]+param[3]+param[4]; + +// Calculate the coefficients of the exponentials + double B0(param[0]), dd(param[1]), d1(param[2]), d2(param[3]), d3(param[4]), l1(param[5]), l2(param[6]); + + double N0(8.0*(l1*l2*cosh(d2/l2)*sinh((d1+d3)/l1)+(l1*l1*cosh(d1/l1)*cosh(d3/l1)+l2*l2*sinh(d1/l1)*sinh(d3/l1))*sinh(d2/l2))); + + fCoeff[0]=B0*exp((dd-d3)/l1-d2/l2)*(-4.0*exp(d3/l1+d2/l2)*l1*l2-exp(d1/l1)*(l1-l2)*(l1+exp(2.0*d3/l1)*(l1-l2)+l2) + \ + exp(d1/l1+2.0*d2/l2)*(l1+l2)*(l1-l2+exp(2.0*d3/l1)*(l1+l2)))/N0; + + fCoeff[1]=B0*exp(-(d1+d3+dd)/l1-d2/l2)*((1.0+exp(2.0*d3/l1))*(-1.0+exp(2.0*d2/l2))*l1*l1-2.0*(1.0-2.0*exp((d1+d3)/l1+d2/l2) + \ + exp(2.0*d2/l2))*l1*l2-(-1.0+exp(2.0*d3/l1))*(-1.0+exp(2.0*d2/l2))*l2*l2)/N0; + + fCoeff[2]=2.0*B0*exp(-(d1+d3)/l1+(d1+dd)/l2)*l1*(-exp(d3/l1)*(l1+exp(2.0*d1/l1)*(l1-l2)+l2) + \ + exp(d1/l1+d2/l2)*(l1-l2+exp(2.0*d3/l1)*(l1+l2)))/N0; + + fCoeff[3]=2.0*B0*exp(-(d1+d3)/l1-(d1+d2+dd)/l2)*l1*(-exp(d1/l1)*(l1+exp(2.0*d3/l1)*(l1-l2)+l2) + \ + exp(d3/l1+d2/l2)*(l1-l2+exp(2.0*d1/l1)*(l1+l2)))/N0; + + fCoeff[4]=B0*exp((d2+dd)/l1-d2/l2)*((1.0+exp(2.0*d1/l1))*(-1.0+exp(2.0*d2/l2))*l1*l1-2.0*(1.0-2.0*exp((d1+d3)/l1+d2/l2) + \ + exp(2.0*d2/l2))*l1*l2-(-1.0+exp(2.0*d1/l1))*(-1.0+exp(2.0*d2/l2))*l2*l2)/N0; + + fCoeff[5]=B0*exp(-(2.0*d1+d2+d3+dd)/l1-d2/l2)*(-4.0*exp(d1/l1+d2/l2)*l1*l2+exp(d3/l1)*(-l1*l1-exp(2.0*d1/l1)*(l1-l2)*(l1-l2)+l2*l2) + \ + exp(d3/l1+2.0*d2/l2)*(l1*l1-l2*l2+exp(2.0*d1/l1)*(l1+l2)*(l1+l2)))/N0; + +// none of the coefficients should be zero + for(unsigned int i(0); i<6; i++) + assert(fCoeff[i]); + + SetBmin(); +} + +double TLondon1D_3LS::GetBofZ(double ZZ) const +{ + if(ZZ < 0. || ZZ < fInterfaces[0] || ZZ > fInterfaces[3]) + return fParam[0]; + + if (ZZ < fInterfaces[1]) { + return fCoeff[0]*exp(-ZZ/fParam[5])+fCoeff[1]*exp(ZZ/fParam[5]); + } else if (ZZ < fInterfaces[2]) { + return fCoeff[2]*exp(-ZZ/fParam[6])+fCoeff[3]*exp(ZZ/fParam[6]); + } else { + return fCoeff[4]*exp(-ZZ/fParam[5])+fCoeff[5]*exp(ZZ/fParam[5]); } } -//------------------ -// Constructor of the TLondon1D_4L class -// 1D-London screening in a thin superconducting film, four layers, four lambdas -// Parameters: Bext[G], deadlayer[nm], thickness1[nm], thickness2[nm], thickness3[nm], thickness4[nm], -// lambda1[nm], lambda2[nm], lambda3[nm], lambda4[nm] -//------------------ - -TLondon1D_4L::TLondon1D_4L(unsigned int steps, const vector ¶m) { - - double N1((param[6]+exp(2.0*param[2]/param[6])*(param[6]-param[7])+param[7])*(-1.0*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])-param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))+exp(2.0*param[3]/param[7])*(param[6]-param[7]+exp(2.0*param[2]/param[6])*(param[6]+param[7]))*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))); - - double N11(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])*(param[7]*cosh(param[3]/param[7])-param[6]*sinh(param[3]/param[7]))+param[7]*(-1.0*param[6]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])*(-1.0*param[6]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))+param[8]*(param[7]*cosh(param[3]/param[7])-param[6]*sinh(param[3]/param[7]))*sinh(param[4]/param[8]))*sinh(param[5]/param[9])); - - double N12(exp(2.0*(param[2]+param[1])/param[6])*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])*(param[7]*cosh(param[3]/param[7])+param[6]*sinh(param[3]/param[7]))+param[7]*(param[6]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])*(param[6]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))+param[8]*(param[7]*cosh(param[3]/param[7])+param[6]*sinh(param[3]/param[7]))*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))); - - double N2((param[6]+param[7]+(param[6]-param[7])*exp(2.0*param[2]/param[6]))*(-1.0*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])-param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))+exp(2.0*param[3]/param[7])*(param[6]-param[7]+exp(2.0*param[2]/param[6])*(param[6]+param[7]))*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))); - - double N21(exp(param[3]/param[7])*param[8]*param[9]*(param[6]*cosh(param[2]/param[6])+param[7]*sinh(param[2]/param[6]))+param[6]*param[9]*cosh(param[5]/param[9])*(-1.0*param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[6]*param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9])); - - double N22(exp((2.0*param[2]+param[3]+2.0*param[1])/param[7])*(-1.0*param[6]*param[8]*param[9]*cosh(param[2]/param[6])+param[7]*param[8]*param[9]*sinh(param[2]/param[6])+exp(param[3]/param[7])*param[6]*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+exp(param[3]/param[7])*param[6]*param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))); - - double N3(2.0*((param[6]+exp(2.0*param[2]/param[6])*(param[6]-param[7])+param[7])*(-1.0*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])-param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))+exp(2.0*param[3]/param[7])*(param[6]-param[7]+exp(2.0*param[2]/param[6])*(param[6]+param[7]))*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9])))); - - double N4(4.0*((param[6]+exp(2.0*param[2]/param[6])*(param[6]-param[7])+param[7])*(-1.0*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])-param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))+exp(2.0*param[3]/param[7])*(param[6]-param[7]+exp(2.0*param[2]/param[6])*(param[6]+param[7]))*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9])))); - - double N42(-1.0*param[6]*param[7]*param[8]+exp(param[5]/param[9])*(param[6]*cosh(param[2]/param[6])*(param[8]*sinh(param[3]/param[7])*(param[9]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))+param[7]*cosh(param[3]/param[7])*(param[8]*cosh(param[4]/param[8])+param[9]*sinh(param[4]/param[8])))+param[7]*sinh(param[2]/param[6])*(param[8]*cosh(param[3]/param[7])*(param[9]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))+param[7]*sinh(param[3]/param[7])*(param[8]*cosh(param[4]/param[8])+param[9]*sinh(param[4]/param[8]))))); - - fDZ = (param[2]+param[3]+param[4]+param[5])/double(steps); - double ZZ, BBz; - - for (unsigned int j(0); j0.); + + double minZ; + // check if the minimum is in the first layer + minZ=-0.5*fParam[5]*log(b_a); + if (minZ > fInterfaces[0] && minZ <= fInterfaces[1]) { + fMinTag = 1; + fMinZ = minZ; + fMinB = GetBofZ(minZ); + return; + } + + double d_c(fCoeff[3]/fCoeff[2]); + assert (d_c>0.); + + // check if the minimum is in the second layer + minZ=-0.5*fParam[6]*log(d_c); + if (minZ > fInterfaces[1] && minZ <= fInterfaces[2]) { + fMinTag = 2; + fMinZ = minZ; + fMinB = GetBofZ(minZ); + return; + } + + double f_e(fCoeff[5]/fCoeff[4]); + assert (f_e>0.); + + // check if the minimum is in the third layer + minZ=-0.5*fParam[5]*log(f_e); + if (minZ > fInterfaces[2] && minZ <= fInterfaces[3]) { + fMinTag = 3; + fMinZ = minZ; + fMinB = GetBofZ(minZ); + return; + } + + assert(fMinZ > 0. && fMinB > 0.); + return; +} + +vector< pair > TLondon1D_3LS::GetInverseAndDerivative(double BB) const +{ + vector< pair > inv; + + if(BB <= fMinB || BB > fParam[0]) + return inv; + + double inverse[4]; + pair invAndDerivative; + + switch(fMinTag) + { + case 1: + inverse[0]=fParam[5]*log((BB-sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]))/(2.0*fCoeff[1])); + inverse[1]=fParam[5]*log((BB+sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]))/(2.0*fCoeff[1])); + inverse[2]=fParam[6]*log((BB+sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]))/(2.0*fCoeff[3])); + inverse[3]=fParam[5]*log((BB+sqrt(BB*BB-4.0*fCoeff[4]*fCoeff[5]))/(2.0*fCoeff[5])); + + if(inverse[0] > fInterfaces[0] && inverse[0] < fMinZ) { + invAndDerivative.first = inverse[0]; + invAndDerivative.second = -fParam[5]/sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]); + inv.push_back(invAndDerivative); + } + if(inverse[1] > fMinZ && inverse[1] <= fInterfaces[1]) { + invAndDerivative.first = inverse[1]; + invAndDerivative.second = +fParam[5]/sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]); + inv.push_back(invAndDerivative); + } + if(inverse[2] > fInterfaces[1] && inverse[2] <= fInterfaces[2]) { + invAndDerivative.first = inverse[2]; + invAndDerivative.second = +fParam[6]/sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]); + inv.push_back(invAndDerivative); + } + if(inverse[3] > fInterfaces[2] && inverse[3] <= fInterfaces[3]) { + invAndDerivative.first = inverse[3]; + invAndDerivative.second = +fParam[5]/sqrt(BB*BB-4.0*fCoeff[4]*fCoeff[5]); + inv.push_back(invAndDerivative); + } + break; + + case 2: + inverse[0]=fParam[5]*log((BB-sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]))/(2.0*fCoeff[1])); + inverse[1]=fParam[6]*log((BB-sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]))/(2.0*fCoeff[3])); + inverse[2]=fParam[6]*log((BB+sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]))/(2.0*fCoeff[3])); + inverse[3]=fParam[5]*log((BB+sqrt(BB*BB-4.0*fCoeff[4]*fCoeff[5]))/(2.0*fCoeff[5])); + + if(inverse[0] > fInterfaces[0] && inverse[0] <= fInterfaces[1]) { + invAndDerivative.first = inverse[0]; + invAndDerivative.second = -fParam[5]/sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]); + inv.push_back(invAndDerivative); + } + if(inverse[1] > fInterfaces[1] && inverse[1] < fMinZ) { + invAndDerivative.first = inverse[1]; + invAndDerivative.second = -fParam[6]/sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]); + inv.push_back(invAndDerivative); + } + if(inverse[2] > fMinZ && inverse[2] <= fInterfaces[2]) { + invAndDerivative.first = inverse[2]; + invAndDerivative.second = +fParam[6]/sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]); + inv.push_back(invAndDerivative); + } + if(inverse[3] > fInterfaces[2] && inverse[3] <= fInterfaces[3]) { + invAndDerivative.first = inverse[3]; + invAndDerivative.second = +fParam[5]/sqrt(BB*BB-4.0*fCoeff[4]*fCoeff[5]); + inv.push_back(invAndDerivative); + } + + break; + + case 3: + inverse[0]=fParam[5]*log((BB-sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]))/(2.0*fCoeff[1])); + inverse[1]=fParam[6]*log((BB-sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]))/(2.0*fCoeff[3])); + inverse[2]=fParam[5]*log((BB-sqrt(BB*BB-4.0*fCoeff[4]*fCoeff[5]))/(2.0*fCoeff[5])); + inverse[3]=fParam[5]*log((BB+sqrt(BB*BB-4.0*fCoeff[4]*fCoeff[5]))/(2.0*fCoeff[5])); + + if(inverse[0] > fInterfaces[0] && inverse[0] <= fInterfaces[1]) { + invAndDerivative.first = inverse[0]; + invAndDerivative.second = -fParam[5]/sqrt(BB*BB-4.0*fCoeff[0]*fCoeff[1]); + inv.push_back(invAndDerivative); + } + if(inverse[1] > fInterfaces[1] && inverse[1] <= fInterfaces[2]) { + invAndDerivative.first = inverse[1]; + invAndDerivative.second = -fParam[6]/sqrt(BB*BB-4.0*fCoeff[2]*fCoeff[3]); + inv.push_back(invAndDerivative); + } + if(inverse[2] > fInterfaces[2] && inverse[2] < fMinZ) { + invAndDerivative.first = inverse[2]; + invAndDerivative.second = -fParam[5]/sqrt(BB*BB-4.0*fCoeff[4]*fCoeff[5]); + inv.push_back(invAndDerivative); + } + if(inverse[3] > fMinZ && inverse[3] <= fInterfaces[3]) { + invAndDerivative.first = inverse[3]; + invAndDerivative.second = +fParam[5]/sqrt(BB*BB-4.0*fCoeff[4]*fCoeff[5]); + inv.push_back(invAndDerivative); + } + + break; + + default: + break; + } + return inv; +} + + +// //------------------ +// // Constructor of the TLondon1D_4L class +// // 1D-London screening in a thin superconducting film, four layers, four lambdas +// // Parameters: Bext[G], deadlayer[nm], thickness1[nm], thickness2[nm], thickness3[nm], thickness4[nm], +// // lambda1[nm], lambda2[nm], lambda3[nm], lambda4[nm] +// //------------------ +// +// TLondon1D_4L::TLondon1D_4L(unsigned int steps, const vector ¶m) { +// +// double N1((param[6]+exp(2.0*param[2]/param[6])*(param[6]-param[7])+param[7])*(-1.0*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])-param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))+exp(2.0*param[3]/param[7])*(param[6]-param[7]+exp(2.0*param[2]/param[6])*(param[6]+param[7]))*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))); +// +// double N11(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])*(param[7]*cosh(param[3]/param[7])-param[6]*sinh(param[3]/param[7]))+param[7]*(-1.0*param[6]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])*(-1.0*param[6]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))+param[8]*(param[7]*cosh(param[3]/param[7])-param[6]*sinh(param[3]/param[7]))*sinh(param[4]/param[8]))*sinh(param[5]/param[9])); +// +// double N12(exp(2.0*(param[2]+param[1])/param[6])*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])*(param[7]*cosh(param[3]/param[7])+param[6]*sinh(param[3]/param[7]))+param[7]*(param[6]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])*(param[6]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))+param[8]*(param[7]*cosh(param[3]/param[7])+param[6]*sinh(param[3]/param[7]))*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))); +// +// double N2((param[6]+param[7]+(param[6]-param[7])*exp(2.0*param[2]/param[6]))*(-1.0*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])-param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))+exp(2.0*param[3]/param[7])*(param[6]-param[7]+exp(2.0*param[2]/param[6])*(param[6]+param[7]))*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))); +// +// double N21(exp(param[3]/param[7])*param[8]*param[9]*(param[6]*cosh(param[2]/param[6])+param[7]*sinh(param[2]/param[6]))+param[6]*param[9]*cosh(param[5]/param[9])*(-1.0*param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[6]*param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9])); +// +// double N22(exp((2.0*param[2]+param[3]+2.0*param[1])/param[7])*(-1.0*param[6]*param[8]*param[9]*cosh(param[2]/param[6])+param[7]*param[8]*param[9]*sinh(param[2]/param[6])+exp(param[3]/param[7])*param[6]*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+exp(param[3]/param[7])*param[6]*param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))); +// +// double N3(2.0*((param[6]+exp(2.0*param[2]/param[6])*(param[6]-param[7])+param[7])*(-1.0*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])-param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))+exp(2.0*param[3]/param[7])*(param[6]-param[7]+exp(2.0*param[2]/param[6])*(param[6]+param[7]))*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9])))); +// +// double N4(4.0*((param[6]+exp(2.0*param[2]/param[6])*(param[6]-param[7])+param[7])*(-1.0*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])-param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))+exp(2.0*param[3]/param[7])*(param[6]-param[7]+exp(2.0*param[2]/param[6])*(param[6]+param[7]))*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9])))); +// +// double N42(-1.0*param[6]*param[7]*param[8]+exp(param[5]/param[9])*(param[6]*cosh(param[2]/param[6])*(param[8]*sinh(param[3]/param[7])*(param[9]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))+param[7]*cosh(param[3]/param[7])*(param[8]*cosh(param[4]/param[8])+param[9]*sinh(param[4]/param[8])))+param[7]*sinh(param[2]/param[6])*(param[8]*cosh(param[3]/param[7])*(param[9]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))+param[7]*sinh(param[3]/param[7])*(param[8]*cosh(param[4]/param[8])+param[9]*sinh(param[4]/param[8]))))); +// +// fDZ = (param[2]+param[3]+param[4]+param[5])/double(steps); +// double ZZ, BBz; +// +// for (unsigned int j(0); j &par) const { fParForPofB[2] = par[1]; // energy // fParForPofB[3] = par[3]; // deadlayer for convolution of field profile - TLondon1D_HS BofZ(fNSteps, fParForBofZ); + TLondon1D_HS BofZ(fParForBofZ); TPofBCalc PofB(BofZ, *fImpProfile, fParForPofB); fPofT->DoFFT(PofB); @@ -367,7 +370,7 @@ double TLondon1D1L::operator()(double t, const vector &par) const { fParForPofB[2] = par[1]; // energy - TLondon1D_1L BofZ1(fNSteps, fParForBofZ); + TLondon1D_1L BofZ1(fParForBofZ); TPofBCalc PofB1(BofZ1, *fImpProfile, fParForPofB); fPofT->DoFFT(PofB1); @@ -524,7 +527,7 @@ double TLondon1D2L::operator()(double t, const vector &par) const { fImpProfile->WeightLayers(par[1], interfaces, weights); } - TLondon1D_2L BofZ2(fNSteps, fParForBofZ); + TLondon1D_2L BofZ2(fParForBofZ); TPofBCalc PofB2(BofZ2, *fImpProfile, fParForPofB); fPofT->DoFFT(PofB2); @@ -687,7 +690,7 @@ double TLondon1D3L::operator()(double t, const vector &par) const { fImpProfile->WeightLayers(par[1], interfaces, weights); } - TLondon1D_3L BofZ3(fNSteps, fParForBofZ); + TLondon1D_3L BofZ3(fParForBofZ); TPofBCalc PofB3(BofZ3, *fImpProfile, fParForPofB); fPofT->DoFFT(PofB3); @@ -836,7 +839,7 @@ double TLondon1D3LS::operator()(double t, const vector &par) const { fImpProfile->WeightLayers(par[1], interfaces, weights); } - TLondon1D_3LS BofZ3S(fNSteps, fParForBofZ); + TLondon1D_3LS BofZ3S(fParForBofZ); TPofBCalc PofB3S(BofZ3S, *fImpProfile, fParForPofB); fPofT->DoFFT(PofB3S); @@ -854,169 +857,169 @@ double TLondon1D3LS::operator()(double t, const vector &par) const { } -//------------------ -// Constructor of the TLondon1D4L class -- reading available implantation profiles and -// creates (a pointer to) the TPofTCalc object (with the FFT plan) -//------------------ - -TLondon1D4L::TLondon1D4L() : fCalcNeeded(true), fFirstCall(true), fLastFourChanged(true) { - // read startup file - string startup_path_name("TFitPofB_startup.xml"); - - TSAXParser *saxParser = new TSAXParser(); - TFitPofBStartupHandler *startupHandler = new TFitPofBStartupHandler(); - saxParser->ConnectToHandler("TFitPofBStartupHandler", startupHandler); - int status (saxParser->ParseFile(startup_path_name.c_str())); - // check for parse errors - if (status) { // error - cout << endl << "**WARNING** reading/parsing TFitPofB_startup.xml failed." << endl; - } - - fNSteps = startupHandler->GetNSteps(); - fWisdom = startupHandler->GetWisdomFile(); - string rge_path(startupHandler->GetDataPath()); - vector energy_vec(startupHandler->GetEnergyList()); - - fParForPofT.push_back(0.0); - fParForPofT.push_back(startupHandler->GetDeltat()); - fParForPofT.push_back(startupHandler->GetDeltaB()); - - fParForPofB.push_back(startupHandler->GetDeltat()); - fParForPofB.push_back(startupHandler->GetDeltaB()); - fParForPofB.push_back(0.0); - - TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); - fImpProfile = x; - x = 0; - delete x; - - TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT); - fPofT = y; - y = 0; - delete y; - - // clean up - if (saxParser) { - delete saxParser; - saxParser = 0; - } - if (startupHandler) { - delete startupHandler; - startupHandler = 0; - } -} - -//------------------ -// TLondon1D4L-Method that calls the procedures to create B(z), p(B) and P(t) -// It finally returns P(t) for a given t. -// Parameters: all the parameters for the function to be fitted through TLondon1D4L -//------------------ - -double TLondon1D4L::operator()(double t, const vector &par) const { - - assert(par.size() == 16); - - if(t<0.0) - return 0.0; - - // check if the function is called the first time and if yes, read in parameters - - if(fFirstCall){ - fPar = par; - -/* for (unsigned int i(0); iWeightLayers(par[1], interfaces, weights); - } - - TLondon1D_4L BofZ4(fNSteps, fParForBofZ); - TPofBCalc PofB4(BofZ4, *fImpProfile, fParForPofB); - fPofT->DoFFT(PofB4); - - }/* else { - cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; - }*/ - - fPofT->CalcPol(fParForPofT); - - fCalcNeeded = false; - fLastFourChanged = false; - } - - return fPofT->Eval(t); - -} +// //------------------ +// // Constructor of the TLondon1D4L class -- reading available implantation profiles and +// // creates (a pointer to) the TPofTCalc object (with the FFT plan) +// //------------------ +// +// TLondon1D4L::TLondon1D4L() : fCalcNeeded(true), fFirstCall(true), fLastFourChanged(true) { +// // read startup file +// string startup_path_name("TFitPofB_startup.xml"); +// +// TSAXParser *saxParser = new TSAXParser(); +// TFitPofBStartupHandler *startupHandler = new TFitPofBStartupHandler(); +// saxParser->ConnectToHandler("TFitPofBStartupHandler", startupHandler); +// int status (saxParser->ParseFile(startup_path_name.c_str())); +// // check for parse errors +// if (status) { // error +// cout << endl << "**WARNING** reading/parsing TFitPofB_startup.xml failed." << endl; +// } +// +// fNSteps = startupHandler->GetNSteps(); +// fWisdom = startupHandler->GetWisdomFile(); +// string rge_path(startupHandler->GetDataPath()); +// vector energy_vec(startupHandler->GetEnergyList()); +// +// fParForPofT.push_back(0.0); +// fParForPofT.push_back(startupHandler->GetDeltat()); +// fParForPofT.push_back(startupHandler->GetDeltaB()); +// +// fParForPofB.push_back(startupHandler->GetDeltat()); +// fParForPofB.push_back(startupHandler->GetDeltaB()); +// fParForPofB.push_back(0.0); +// +// TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); +// fImpProfile = x; +// x = 0; +// delete x; +// +// TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT); +// fPofT = y; +// y = 0; +// delete y; +// +// // clean up +// if (saxParser) { +// delete saxParser; +// saxParser = 0; +// } +// if (startupHandler) { +// delete startupHandler; +// startupHandler = 0; +// } +// } +// +// //------------------ +// // TLondon1D4L-Method that calls the procedures to create B(z), p(B) and P(t) +// // It finally returns P(t) for a given t. +// // Parameters: all the parameters for the function to be fitted through TLondon1D4L +// //------------------ +// +// double TLondon1D4L::operator()(double t, const vector &par) const { +// +// assert(par.size() == 16); +// +// if(t<0.0) +// return 0.0; +// +// // check if the function is called the first time and if yes, read in parameters +// +// if(fFirstCall){ +// fPar = par; +// +// /* for (unsigned int i(0); iWeightLayers(par[1], interfaces, weights); +// } +// +// TLondon1D_4L BofZ4(fNSteps, fParForBofZ); +// TPofBCalc PofB4(BofZ4, *fImpProfile, fParForPofB); +// fPofT->DoFFT(PofB4); +// +// }/* else { +// cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; +// }*/ +// +// fPofT->CalcPol(fParForPofT); +// +// fCalcNeeded = false; +// fLastFourChanged = false; +// } +// +// return fPofT->Eval(t); +// +// } //------------------ // Constructor of the TLondon1D3LSub class -- reading available implantation profiles and @@ -1024,6 +1027,10 @@ double TLondon1D4L::operator()(double t, const vector &par) const { //------------------ TLondon1D3LSub::TLondon1D3LSub() : fCalcNeeded(true), fFirstCall(true), fWeightsChanged(true) { +// omp_set_nested(1); +// omp_set_dynamic(1); +// omp_set_num_threads(4); + // read startup file string startup_path_name("TFitPofB_startup.xml"); @@ -1085,7 +1092,10 @@ double TLondon1D3LSub::operator()(double t, const vector &par) const { // check if the function is called the first time and if yes, read in parameters +//#pragma omp critical +//{ if(fFirstCall){ + fPar = par; /* for (unsigned int i(0); i &par) const { fFirstCall=false; } + // check if any parameter has changed bool par_changed(false); @@ -1164,7 +1175,7 @@ double TLondon1D3LSub::operator()(double t, const vector &par) const { fImpProfile->WeightLayers(par[1], interfaces, weights); } - TLondon1D_3L BofZ3(fNSteps, fParForBofZ); + TLondon1D_3L BofZ3(fParForBofZ); TPofBCalc PofB3(BofZ3, *fImpProfile, fParForPofB); // Add background contribution from the substrate @@ -1181,9 +1192,19 @@ double TLondon1D3LSub::operator()(double t, const vector &par) const { fCalcNeeded = false; fWeightsChanged = false; - } - + } +//} return fPofT->Eval(t); } + + +double TLondon1D3Lestimate::operator()(double z, const vector& par) const { + + TLondon1D_3L BofZ(par); + + return BofZ.GetBofZ(z); + +} + diff --git a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp index 1739fdd4..34618f83 100644 --- a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp +++ b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp @@ -13,18 +13,189 @@ #include #include #include +#include /* USED FOR DEBUGGING----------------------------------- #include #include /-------------------------------------------------------*/ +// Do not actually calculate P(B) but take it from a B and a P(B) vector of the same size + +TPofBCalc::TPofBCalc( const vector& b, const vector& pb, double dt) { + assert(b.size() == pb.size() && b.size() >= 2); + + fB = b; + fPB = pb; + + vector::const_iterator iter, iterB; + iterB = b.begin(); + + for(iter = pb.begin(); iter != pb.end(); iter++){ + if(*iter != 0.0) { + fBmin = *iterB; + cout << fBmin << endl; + break; + } + iterB++; + } + + for( ; iter != b.end(); iter++){ + if(*iter == 0.0) { + fBmax = *(iterB-1); + cout << fBmax << endl; + break; + } + iterB++; + } + + fDT = dt; // needed if a convolution should be done + fDB = b[1]-b[0]; + + cout << fDB << endl; +} + + +TPofBCalc::TPofBCalc( const string &type, const vector ¶ ) : fDT(para[0]), fDB(para[1]) { + +if (type == "skg"){ // skewed Gaussian + + fBmin = 0.0; + fBmax = para[2]/gBar+10.0*fabs(para[4])/(2.0*pi*gBar); + + double BB; + double expominus(para[3]*para[3]/(2.0*pi*pi*gBar*gBar)); + double expoplus(para[4]*para[4]/(2.0*pi*pi*gBar*gBar)); + double B0(para[2]/gBar); + double BmaxFFT(1.0/gBar/fDT); + + for ( BB = 0.0 ; BB < B0 ; BB += fDB ) { + fB.push_back(BB); + fPB.push_back(exp(-(BB-B0)*(BB-B0)/expominus)); + } + + for ( ; BB < fBmax ; BB += fDB ) { + fB.push_back(BB); + fPB.push_back(exp(-(BB-B0)*(BB-B0)/expoplus)); + } + + unsigned int lastZerosStart(fB.size()); + + for ( ; BB <= BmaxFFT ; BB += fDB ) { + fB.push_back(BB); + fPB.push_back(0.0); + } + + // make sure that we have an even number of elements in p(B) for FFT, so we do not have to care later + + if (fB.size() % 2) { + fB.push_back(BB); + fPB.push_back(0.0); + } + + // normalize p(B) + + double pBsum = 0.0; + for (unsigned int i(0); i ¶ ) : fDT(para[0]), fDB(para[1]) +{ + + fBmin = BofZ.GetBmin(); + fBmax = BofZ.GetBmax(); + + double BB; + + // fill not used Bs before Bmin with 0.0 + + for (BB = 0.0; BB < fBmin; BB += fDB) + fB.push_back(BB); + fPB.resize(fB.size(),0.0); + + unsigned int firstZerosEnd(fB.size()); + + // calculate p(B) from the inverse of B(z) + + for ( ; BB <= fBmax ; BB += fDB) { + fB.push_back(BB); + fPB.push_back(0.0); + + vector< pair > inv; + inv = BofZ.GetInverseAndDerivative(BB); + + for (unsigned int i(0); i ¶ ) : fDT(para[0]), fDB(para[1]) { +TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, const vector ¶, unsigned int zonk ) : fDT(para[0]), fDB(para[1]) { fBmin = BofZ.GetBmin(); fBmax = BofZ.GetBmax(); @@ -34,7 +205,7 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons // fill not used Bs before Bmin with 0.0 - for ( BB = 0.0 ; BB < fBmin ; BB += fDB ) { + for (BB = 0.0; BB < fBmin; BB += fDB) { fB.push_back(BB); fPB.push_back(0.0); } @@ -45,7 +216,7 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons vector bofzZ(BofZ.DataZ()); vector bofzBZ(BofZ.DataBZ()); - double ddZ(BofZ.GetdZ()); + double ddZ(BofZ.GetDZ()); /* USED FOR DEBUGGING----------------------------------- cout << "Bmin = " << fBmin << ", Bmax = " << fBmax << endl; diff --git a/src/external/TFitPofB-lib/classes/TPofTCalc.cpp b/src/external/TFitPofB-lib/classes/TPofTCalc.cpp index 978c1cf8..72ada3ef 100644 --- a/src/external/TFitPofB-lib/classes/TPofTCalc.cpp +++ b/src/external/TFitPofB-lib/classes/TPofTCalc.cpp @@ -41,6 +41,8 @@ #include #include +#include + #include #include #include @@ -62,9 +64,25 @@ TPofTCalc::TPofTCalc (const string &wisdom, const vector &par) : fWisdom(wisdom) { + int init_threads(fftw_init_threads()); + if (!init_threads) + cout << "TPofTCalc::TPofTCalc: Couldn't initialize multiple FFTW-threads ..." << endl; + else + fftw_plan_with_nthreads(2); + fNFFT = ( int(1.0/gBar/par[1]/par[2]+1.0) % 2 ) ? int(1.0/gBar/par[1]/par[2]+2.0) : int(1.0/gBar/par[1]/par[2]+1.0); fTBin = 1.0/gBar/double(fNFFT-1)/par[2]; + fT.resize(fNFFT/2+1); + fPT.resize(fNFFT/2+1); + + int i; + +#pragma omp parallel for default(shared) private(i) schedule(dynamic) + for (i=0; i &par) { double sinph(sin(par[0]*PI/180.0)), cosph(cos(par[0]*PI/180.0)); + int i; - double polTemp(0.0); - fT.clear(); - fPT.clear(); - - for (unsigned int i(0); i //--------------------- double TPofTCalc::Eval(double t) const { - for (unsigned int i(0); i -#include using namespace std; //-------------------- // Base class for any kind of theory function B(z) -// In principle only constructors for the different models have to be implemented //-------------------- class TBofZCalc { @@ -30,33 +28,52 @@ public: virtual ~TBofZCalc() { fZ.clear(); fBZ.clear(); + fParam.clear(); } - virtual vector DataZ() const; - virtual vector DataBZ() const; - virtual void Calculate() = 0; + virtual vector DataZ() const {return fZ;} + virtual vector DataBZ() const {return fBZ;} + virtual void Calculate(); virtual double GetBofZ(double) const = 0; virtual double GetBmin() const = 0; virtual double GetBmax() const = 0; - double GetdZ() const {return fDZ;} + double GetDZ() const {return fDZ;} protected: - vector fZ; - vector fBZ; + int fSteps; double fDZ; vector fParam; - unsigned int fSteps; + vector fZ; + vector fBZ; +}; + +//-------------------- +// Base class for any kind of theory function B(z) where the inverse and its derivative are given analytically +//-------------------- + +class TBofZCalcInverse : public TBofZCalc { + +public: + + TBofZCalcInverse() {} + virtual ~TBofZCalcInverse() {} + + virtual vector< pair > GetInverseAndDerivative(double) const = 0; }; //-------------------- // Class "for Meissner screening" in a superconducting half-space //-------------------- -class TLondon1D_HS : public TBofZCalc { +class TLondon1D_HS : public TBofZCalcInverse { public: - TLondon1D_HS(unsigned int, const vector& ); + TLondon1D_HS(const vector&, unsigned int steps = 3000); + double GetBofZ(double) const; + double GetBmin() const; + double GetBmax() const; + vector< pair > GetInverseAndDerivative(double) const; }; @@ -64,11 +81,22 @@ public: // Class "for Meissner screening" in a thin superconducting film //-------------------- -class TLondon1D_1L : public TBofZCalc { +class TLondon1D_1L : public TBofZCalcInverse { public: - TLondon1D_1L(unsigned int, const vector& ); + TLondon1D_1L(const vector&, unsigned int steps = 3000); + double GetBofZ(double) const; + double GetBmin() const; + double GetBmax() const; + vector< pair > GetInverseAndDerivative(double) const; + +private: + void SetBmin(); + + double fMinZ; + double fMinB; + double fCoeff[2]; }; @@ -76,28 +104,15 @@ public: // Class "for Meissner screening" in a thin superconducting film - bilayer with two different lambdas //-------------------- -class TLondon1D_2L : public TBofZCalc { +class TLondon1D_2L : public TBofZCalcInverse { public: - TLondon1D_2L(unsigned int, const vector& ); - -}; - -//-------------------- -// Class "for Meissner screening" in a thin superconducting film - tri-layer with three different lambdas -//-------------------- - -class TLondon1D_3L : public TBofZCalc { - -public: - - TLondon1D_3L(unsigned int, const vector& ); - void Calculate(); + TLondon1D_2L(const vector&, unsigned int steps = 3000); double GetBofZ(double) const; - vector< pair > GetInverseAndDerivative(double) const; double GetBmin() const; double GetBmax() const; + vector< pair > GetInverseAndDerivative(double) const; private: void SetBmin(); @@ -105,33 +120,68 @@ private: int fMinTag; double fMinZ; double fMinB; - double fCoeff[6]; - double fInterfaces[4]; + double fInterfaces[3]; + double fCoeff[4]; +}; +//-------------------- +// Class "for Meissner screening" in a thin superconducting film - tri-layer with three different lambdas +//-------------------- + +class TLondon1D_3L : public TBofZCalcInverse { + +public: + + TLondon1D_3L(const vector&, unsigned int steps = 3000); + double GetBofZ(double) const; + double GetBmin() const; + double GetBmax() const; + vector< pair > GetInverseAndDerivative(double) const; + +private: + void SetBmin(); + + int fMinTag; + double fMinZ; + double fMinB; + double fInterfaces[4]; + double fCoeff[6]; }; //-------------------- // Class "for Meissner screening" in a thin superconducting film - tri-layer with two different lambdas //-------------------- -class TLondon1D_3LS : public TBofZCalc { +class TLondon1D_3LS : public TBofZCalcInverse { public: - TLondon1D_3LS(unsigned int, const vector& ); + TLondon1D_3LS(const vector&, unsigned int steps = 3000); + double GetBofZ(double) const; + double GetBmin() const; + double GetBmax() const; + vector< pair > GetInverseAndDerivative(double) const; +private: + void SetBmin(); + + int fMinTag; + double fMinZ; + double fMinB; + double fInterfaces[4]; + double fCoeff[6]; }; -//-------------------- -// Class "for Meissner screening" in a thin superconducting film - four layers with four different lambdas -//-------------------- - -class TLondon1D_4L : public TBofZCalc { - -public: - - TLondon1D_4L(unsigned int, const vector& ); - -}; +// //-------------------- +// // Class "for Meissner screening" in a thin superconducting film - four layers with four different lambdas +// //-------------------- +// +// class TLondon1D_4L : public TBofZCalc { +// +// public: +// +// TLondon1D_4L(unsigned int, const vector& ); +// +// }; #endif // _BofZCalc_H_ diff --git a/src/external/TFitPofB-lib/include/TLondon1D.h b/src/external/TFitPofB-lib/include/TLondon1D.h index f151ee50..9faf330b 100644 --- a/src/external/TFitPofB-lib/include/TLondon1D.h +++ b/src/external/TFitPofB-lib/include/TLondon1D.h @@ -140,30 +140,30 @@ private: ClassDef(TLondon1D3LS,1) }; -class TLondon1D4L : public PUserFcnBase { - -public: - // default constructor - TLondon1D4L(); - ~TLondon1D4L(); - - double operator()(double, const vector&) const; - -private: - mutable vector fPar; - TTrimSPData *fImpProfile; - TPofTCalc *fPofT; - mutable bool fCalcNeeded; - mutable bool fFirstCall; - mutable vector fParForPofT; - mutable vector fParForBofZ; - mutable vector fParForPofB; - string fWisdom; - unsigned int fNSteps; - mutable bool fLastFourChanged; - - ClassDef(TLondon1D4L,1) -}; +// class TLondon1D4L : public PUserFcnBase { +// +// public: +// // default constructor +// TLondon1D4L(); +// ~TLondon1D4L(); +// +// double operator()(double, const vector&) const; +// +// private: +// mutable vector fPar; +// TTrimSPData *fImpProfile; +// TPofTCalc *fPofT; +// mutable bool fCalcNeeded; +// mutable bool fFirstCall; +// mutable vector fParForPofT; +// mutable vector fParForBofZ; +// mutable vector fParForPofB; +// string fWisdom; +// unsigned int fNSteps; +// mutable bool fLastFourChanged; +// +// ClassDef(TLondon1D4L,1) +// }; class TLondon1D3LSub : public PUserFcnBase { @@ -190,4 +190,18 @@ private: ClassDef(TLondon1D3LSub,1) }; +// Class for fitting directly B(z) without any P(B)-calculation + +class TLondon1D3Lestimate : public PUserFcnBase { + +public: + // default constructor + TLondon1D3Lestimate() {} + ~TLondon1D3Lestimate() {} + + double operator()(double, const vector&) const; + + ClassDef(TLondon1D3Lestimate,1) +}; + #endif //_TLondon1D_H_ diff --git a/src/external/TFitPofB-lib/include/TLondon1DLinkDef.h b/src/external/TFitPofB-lib/include/TLondon1DLinkDef.h index 2097638d..8b723dac 100644 --- a/src/external/TFitPofB-lib/include/TLondon1DLinkDef.h +++ b/src/external/TFitPofB-lib/include/TLondon1DLinkDef.h @@ -21,9 +21,11 @@ #pragma link C++ class TLondon1D2L+; #pragma link C++ class TLondon1D3L+; #pragma link C++ class TLondon1D3LS+; -#pragma link C++ class TLondon1D4L+; +//#pragma link C++ class TLondon1D4L+; #pragma link C++ class TLondon1D3LSub+; +#pragma link C++ class TLondon1D3Lestimate; + #endif //__CINT__ // root dictionary stuff -------------------------------------------------- diff --git a/src/external/TFitPofB-lib/include/TPofBCalc.h b/src/external/TFitPofB-lib/include/TPofBCalc.h index 5a273cef..101b5ccf 100644 --- a/src/external/TFitPofB-lib/include/TPofBCalc.h +++ b/src/external/TFitPofB-lib/include/TPofBCalc.h @@ -22,7 +22,10 @@ class TPofBCalc { public: - TPofBCalc( const TBofZCalc&, const TTrimSPData&, const vector& ); + TPofBCalc( const string&, const vector& ); + TPofBCalc( const TBofZCalc&, const TTrimSPData&, const vector&, unsigned int ); + TPofBCalc( const TBofZCalcInverse&, const TTrimSPData&, const vector& ); + TPofBCalc( const vector&, const vector& , double dt = 0.01 ); ~TPofBCalc() { fB.clear(); fPB.clear(); diff --git a/src/external/TFitPofB-lib/include/TPofTCalc.h b/src/external/TFitPofB-lib/include/TPofTCalc.h index 83cf8bb8..a81d86b8 100644 --- a/src/external/TFitPofB-lib/include/TPofTCalc.h +++ b/src/external/TFitPofB-lib/include/TPofTCalc.h @@ -40,7 +40,7 @@ private: vector fT; vector fPT; double fTBin; - unsigned int fNFFT; + int fNFFT; const string fWisdom; };