LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 4 Activity

Added some more simple triangular Vortex lattice models to libTFitPofB - UNTESTEDklog wojek NO GUARANTEE, esp for AGL

Browse Source
This commit is contained in:
Bastian M. Wojek
2009-11-16 21:58:35 +00:00
parent 08efd371ac
commit ad673f7d36
8 changed files with 841 additions and 171 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
src/external/TFitPofB-lib/classes/Makefile.libTFitPofB vendored
View File

@ -38,7 +38,7 @@ endif
# -- Linux # -- Linux
ifeq ($(OS),LINUX) ifeq ($(OS),LINUX)
CXX = g++ CXX = g++
CXXFLAGS = -O3 -Wall -Wno-trigraphs -fPIC CXXFLAGS = -g -O3 -Wall -Wno-trigraphs -fPIC
PMUSRPATH = ../../../include PMUSRPATH = ../../../include
MNPATH = $(ROOTSYS)/include MNPATH = $(ROOTSYS)/include
LOCALPATH = ../include LOCALPATH = ../include
@ -109,7 +109,7 @@ endif
# some definitions: headers (used to generate *Dict* stuff), sources, objects,... # some definitions: headers (used to generate *Dict* stuff), sources, objects,...
OBJS = OBJS =
OBJS += TTrimSPDataHandler.o OBJS += TTrimSPDataHandler.o
OBJS += TBulkVortexFieldCalc.o OBJS += TBulkTriVortexFieldCalc.o
OBJS += TBofZCalc.o OBJS += TBofZCalc.o
OBJS += TPofBCalc.o OBJS += TPofBCalc.o
OBJS += TPofTCalc.o OBJS += TPofTCalc.o
@ -120,7 +120,7 @@ OBJS += TSkewedGss.o TSkewedGssDict.o
INST_HEADER = INST_HEADER =
INST_HEADER += ../include/TBofZCalc.h INST_HEADER += ../include/TBofZCalc.h
INST_HEADER += ../include/TBulkVortexFieldCalc.h INST_HEADER += ../include/TBulkTriVortexFieldCalc.h
INST_HEADER += ../include/TFitPofBStartupHandler.h INST_HEADER += ../include/TFitPofBStartupHandler.h
INST_HEADER += ../include/TLondon1D.h INST_HEADER += ../include/TLondon1D.h
INST_HEADER += ../include/TPofBCalc.h INST_HEADER += ../include/TPofBCalc.h

591
src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp → src/external/TFitPofB-lib/classes/TBulkTriVortexFieldCalc.cpp vendored
View File

@ -1,6 +1,6 @@
/*************************************************************************** /***************************************************************************
TBulkVortexFieldCalc.cpp TBulkTriVortexFieldCalc.cpp
Author: Bastian M. Wojek, Alexander Maisuradze Author: Bastian M. Wojek, Alexander Maisuradze
e-mail: bastian.wojek@psi.ch e-mail: bastian.wojek@psi.ch
@ -29,7 +29,7 @@
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/ ***************************************************************************/
#include "TBulkVortexFieldCalc.h" #include "TBulkTriVortexFieldCalc.h"
#include <cstdlib> #include <cstdlib>
#include <cmath> #include <cmath>
@ -37,25 +37,25 @@
#include <iostream> #include <iostream>
using namespace std; using namespace std;
#include "TMath.h"
#define PI 3.14159265358979323846 #define PI 3.14159265358979323846
#define TWOPI 6.28318530717958647692 #define TWOPI 6.28318530717958647692
const double fluxQuantum(2.067833667e7); // in this case this is Gauss times square nm const double fluxQuantum(2.067833667e7); // in this case this is Gauss times square nm
const double sqrt3(sqrt(3.0)); const double sqrt3(sqrt(3.0));
//const double pi_4sqrt3(0.25*PI/sqrt3);
const double pi_4sqrt3(0.25*PI/sqrt(3.0)); const double pi_4sqrt3(0.25*PI/sqrt(3.0));
//const double pi_4sqrt3(PI*sqrt(3.0));
double getXi(const double hc2) { // get xi given Hc2 in Gauss
double getXi(const double &hc2) { // get xi given Hc2 in Gauss
if (hc2) if (hc2)
return sqrt(fluxQuantum/(TWOPI*hc2)); return sqrt(fluxQuantum/(TWOPI*hc2));
else else
return 0.0; return 0.0;
} }
double getHc2(const double xi) { // get Hc2 given xi in nm double getHc2(const double &xi) { // get Hc2 given xi in nm
if (xi) if (xi)
return fluxQuantum/(TWOPI*xi*xi); return fluxQuantum/(TWOPI*xi*xi);
else else
@ -163,78 +163,17 @@ TBulkTriVortexLondonFieldCalc::TBulkTriVortexLondonFieldCalc(const string& wisdo
fFFTplan = fftw_plan_dft_c2r_2d(fSteps, fSteps, fFFTin, fFFTout, FFTW_ESTIMATE); fFFTplan = fftw_plan_dft_c2r_2d(fSteps, fSteps, fFFTin, fFFTout, FFTW_ESTIMATE);
} }
// double TBulkTriVortexLondonFieldCalc::GetBatPoint(double relX, double relY) const {
//
// double field(fParam[0]), lambda(fParam[1]), xi(fParam[2]);
// double xisq_2(0.5*xi*xi), lambdasq(lambda*lambda);
//
// double latConstTr(sqrt(fluxQuantum/field)*sqrt(sqrt(4.0/3.0)));
// double Hc2(getHc2(xi));
//
// double xCoord(latConstTr*relX); // in nanometers
// double yCoord(latConstTr*relY);
//
// if ((field < Hc2) && (lambda > xi/sqrt(2.0))) {
// double KLatTr(4.0*PI/(latConstTr*sqrt3));
// double fourierSum(1.0), fourierAdd(0.0), expo;
//
// // k = 0, l = 0 gives 1.0, already in fourierSum
//
// // l = 0, only integer k
// for (double k (1.0); k < static_cast<double>(fNFourierComp); k+=1.0){
// expo = 3.0*pow(KLatTr*k, 2.0);
// fourierAdd += exp(-xisq_2*expo)/(lambdasq*expo + 1.0)*cos(sqrt3*KLatTr*k*xCoord);
// }
//
// fourierSum += 2.0*fourierAdd;
// fourierAdd = 0.0;
//
// // k = 0, only integer l
// for (double l (1.0); l < static_cast<double>(fNFourierComp); l+=1.0){
// expo = pow(KLatTr*l, 2.0);
// fourierAdd += exp(-xisq_2*expo)/(lambdasq*expo + 1.0)*cos(KLatTr*l*yCoord);
// }
//
// fourierSum += 2.0*fourierAdd;
// fourierAdd = 0.0;
//
// // k != 0, l != 0 both integers
// for (double k (1.0); k < static_cast<double>(fNFourierComp); k+=1.0){
// for (double l (1.0); l < static_cast<double>(fNFourierComp); l+=1.0){
// expo = KLatTr*KLatTr*(3.0*k*k + l*l);
// fourierAdd += exp(-xisq_2*expo)/(lambdasq*expo + 1.0)*cos(sqrt3*KLatTr*k*xCoord)*cos(KLatTr*l*yCoord);
// }
// }
//
// fourierSum += 4.0*fourierAdd;
// fourierAdd = 0.0;
//
// // k != 0, l != 0 both half-integral numbers
// for (double k (0.5); k < static_cast<double>(fNFourierComp); k+=1.0){
// for (double l (0.5); l < static_cast<double>(fNFourierComp); l+=1.0){
// expo = KLatTr*KLatTr*(3.0*k*k + l*l);
// fourierAdd += exp(-xisq_2*expo)/(lambdasq*expo + 1.0)*cos(sqrt3*KLatTr*k*xCoord)*cos(KLatTr*l*yCoord);
// }
// }
//
// fourierSum += 4.0*fourierAdd;
//
// // for(int mm = -fNFourierComp; mm <= static_cast<int>(fNFourierComp); mm++) {
// // for(int nn = -fNFourierComp; nn <= static_cast<int>(fNFourierComp); nn++) {
// // fourierSum += cos(KLatTr*(xCoord*mm*(0.5*sqrt(3.0)) + yCoord*mm*0.5 + yCoord*nn))*exp(-(0.5*fParam[1]*fParam[1]*KLatTr*KLatTr)*
// // (0.75*mm*mm + (nn + 0.5*mm)*(nn + 0.5*mm)))/(1.0+(fParam[0]*KLatTr*fParam[0]*KLatTr)*(0.75*mm*mm + (nn + 0.5*mm)*(nn + 0.5*mm)));
// // }
// // }
// // cout << " " << fourierSum << ", ";
// return field*fourierSum;
// }
// else
// return field;
//
// }
void TBulkTriVortexLondonFieldCalc::CalculateGrid() const { void TBulkTriVortexLondonFieldCalc::CalculateGrid() const {
// SetParameters - method has to be called from the user before the calculation!! // SetParameters - method has to be called from the user before the calculation!!
if (fParam.size() < 3) {
cout << endl << "The SetParameters-method has to be called before B(x,y) can be calculated!" << endl;
return;
}
if (!fParam[0] || !fParam[1] || !fParam[2]) {
cout << endl << "The field, penetration depth and coherence length have to have finite values in order to calculate B(x,y)!" << endl;
return;
}
double field(fabs(fParam[0])), lambda(fabs(fParam[1])), xi(fabs(fParam[2])); double field(fabs(fParam[0])), lambda(fabs(fParam[1])), xi(fabs(fParam[2]));
double Hc2(getHc2(xi)); double Hc2(getHc2(xi));
@ -249,8 +188,8 @@ void TBulkTriVortexLondonFieldCalc::CalculateGrid() const {
// ... but first check that the field is not larger than Hc2 and that we are dealing with a type II SC // ... but first check that the field is not larger than Hc2 and that we are dealing with a type II SC
if ((field >= Hc2) || (lambda < xi/sqrt(2.0))) { if ((field >= Hc2) || (lambda < xi/sqrt(2.0))) {
int m; int m;
#pragma omp parallel for default(shared) private(m) schedule(dynamic) #pragma omp parallel for default(shared) private(m) schedule(dynamic)
for (m = 0; m < NFFTsq; m++) { for (m = 0; m < NFFTsq; m++) {
fFFTout[m] = field; fFFTout[m] = field;
} }
@ -260,90 +199,80 @@ void TBulkTriVortexLondonFieldCalc::CalculateGrid() const {
} }
// ... now fill in the Fourier components if everything was okay above // ... now fill in the Fourier components if everything was okay above
double Gsq; double Gsq, ll;
int k, l, lNFFT_2; int k, l, lNFFT_2;
// omp causes problems with the fftw_complex*... comment it out for the moment for (l = 0; l < NFFT_2; l += 2) {
//#pragma omp parallel default(shared) private(l,lNFFT_2,k,Gsq) lNFFT_2 = l*(NFFT_2 + 1);
{ ll = 3.0*static_cast<double>(l*l);
// #pragma omp sections nowait for (k = 0; k < NFFT_2; k += 2) {
{ Gsq = static_cast<double>(k*k) + ll;
// #pragma omp section fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
// #pragma omp parallel for default(shared) private(l,lNFFT_2,k) schedule(dynamic) fFFTin[lNFFT_2 + k][1] = 0.0;
for (l = 0; l < NFFT_2; l += 2) { fFFTin[lNFFT_2 + k + 1][0] = 0.0;
lNFFT_2 = l*(NFFT_2 + 1); fFFTin[lNFFT_2 + k + 1][1] = 0.0;
for (k = 0; k < NFFT_2; k += 2) { }
Gsq = 3.0*static_cast<double>(k*k) + static_cast<double>(l*l); k = NFFT_2;
fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq); Gsq = static_cast<double>(k*k) + ll;
fFFTin[lNFFT_2 + k][1] = 0.0; fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k + 1][0] = 0.0; fFFTin[lNFFT_2 + k][1] = 0.0;
fFFTin[lNFFT_2 + k + 1][1] = 0.0; }
}
k = NFFT_2;
Gsq = 3.0*static_cast<double>(k*k) + static_cast<double>(l*l);
fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][1] = 0.0;
}
// #pragma omp section
// #pragma omp parallel for default(shared) private(l,lNFFT_2,k) schedule(dynamic)
for (l = NFFT_2; l < NFFT; l += 2) {
lNFFT_2 = l*(NFFT_2 + 1);
for (k = 0; k < NFFT_2; k += 2) {
Gsq = 3.0*static_cast<double>(k*k) + static_cast<double>((NFFT-l)*(NFFT-l));
fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][1] = 0.0;
fFFTin[lNFFT_2 + k + 1][0] = 0.0;
fFFTin[lNFFT_2 + k + 1][1] = 0.0;
}
k = NFFT_2;
Gsq = 3.0*static_cast<double>(k*k) + static_cast<double>((NFFT-l)*(NFFT-l));
fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][1] = 0.0;
}
// intermediate rows for (l = NFFT_2; l < NFFT; l += 2) {
// #pragma omp section lNFFT_2 = l*(NFFT_2 + 1);
// #pragma omp parallel for default(shared) private(l,lNFFT_2,k) schedule(dynamic) ll = 3.0*static_cast<double>((NFFT-l)*(NFFT-l));
for (l = 1; l < NFFT_2; l += 2) { for (k = 0; k < NFFT_2; k += 2) {
lNFFT_2 = l*(NFFT_2 + 1); Gsq = static_cast<double>(k*k) + ll;
for (k = 0; k < NFFT_2; k += 2) { fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
Gsq = 3.0*static_cast<double>((k + 1)*(k + 1)) + static_cast<double>(l*l); fFFTin[lNFFT_2 + k][1] = 0.0;
fFFTin[lNFFT_2 + k][0] = 0.0; fFFTin[lNFFT_2 + k + 1][0] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0; fFFTin[lNFFT_2 + k + 1][1] = 0.0;
fFFTin[lNFFT_2 + k + 1][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq); }
fFFTin[lNFFT_2 + k + 1][1] = 0.0; k = NFFT_2;
} Gsq = static_cast<double>(k*k) + ll;
k = NFFT_2; fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][0] = 0.0; fFFTin[lNFFT_2 + k][1] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0; }
}
// #pragma omp section // intermediate rows
// #pragma omp parallel for default(shared) private(l,lNFFT_2,k) schedule(dynamic)
for (l = NFFT_2 + 1; l < NFFT; l += 2) {
lNFFT_2 = l*(NFFT_2 + 1);
for (k = 0; k < NFFT_2; k += 2) {
Gsq = 3.0*static_cast<double>((k+1)*(k+1)) + static_cast<double>((NFFT-l)*(NFFT-l));
fFFTin[lNFFT_2 + k][0] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0;
fFFTin[lNFFT_2 + k + 1][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k + 1][1] = 0.0;
}
k = NFFT_2;
fFFTin[lNFFT_2 + k][0] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0;
}
} /* end of sections */
} /* end of parallel section */ for (l = 1; l < NFFT_2; l += 2) {
lNFFT_2 = l*(NFFT_2 + 1);
ll = 3.0*static_cast<double>(l*l);
for (k = 0; k < NFFT_2; k += 2) {
Gsq = static_cast<double>((k + 1)*(k + 1)) + ll;
fFFTin[lNFFT_2 + k][0] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0;
fFFTin[lNFFT_2 + k + 1][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k + 1][1] = 0.0;
}
k = NFFT_2;
fFFTin[lNFFT_2 + k][0] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0;
}
for (l = NFFT_2 + 1; l < NFFT; l += 2) {
lNFFT_2 = l*(NFFT_2 + 1);
ll = 3.0*static_cast<double>((NFFT-l)*(NFFT-l));
for (k = 0; k < NFFT_2; k += 2) {
Gsq = static_cast<double>((k+1)*(k+1)) + ll;
fFFTin[lNFFT_2 + k][0] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0;
fFFTin[lNFFT_2 + k + 1][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k + 1][1] = 0.0;
}
k = NFFT_2;
fFFTin[lNFFT_2 + k][0] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0;
}
// Do the Fourier transform to get B(x,y) // Do the Fourier transform to get B(x,y)
fftw_execute(fFFTplan); fftw_execute(fFFTplan);
// Multiply by the applied field // Multiply by the applied field
#pragma omp parallel for default(shared) private(l) schedule(dynamic) #pragma omp parallel for default(shared) private(l) schedule(dynamic)
for (l = 0; l < NFFTsq; l++) { for (l = 0; l < NFFTsq; l++) {
fFFTout[l] *= field; fFFTout[l] *= field;
} }
@ -355,6 +284,353 @@ void TBulkTriVortexLondonFieldCalc::CalculateGrid() const {
} }
TBulkTriVortexMLFieldCalc::TBulkTriVortexMLFieldCalc(const string& wisdom, const unsigned int steps) {
fWisdom = wisdom;
if (steps % 2) {
fSteps = steps + 1;
} else {
fSteps = steps;
}
fParam.resize(3);
fGridExists = false;
int init_threads(fftw_init_threads());
if (init_threads)
fftw_plan_with_nthreads(2);
fFFTin = new fftw_complex[(fSteps/2 + 1) * fSteps];
fFFTout = new double[fSteps*fSteps];
// cout << "Check for the FFT plan..." << endl;
// Load wisdom from file if it exists and should be used
fUseWisdom = true;
int wisdomLoaded(0);
FILE *wordsOfWisdomR;
wordsOfWisdomR = fopen(fWisdom.c_str(), "r");
if (wordsOfWisdomR == NULL) {
fUseWisdom = false;
} else {
wisdomLoaded = fftw_import_wisdom_from_file(wordsOfWisdomR);
fclose(wordsOfWisdomR);
}
if (!wisdomLoaded) {
fUseWisdom = false;
}
// create the FFT plan
if (fUseWisdom)
fFFTplan = fftw_plan_dft_c2r_2d(fSteps, fSteps, fFFTin, fFFTout, FFTW_EXHAUSTIVE);
else
fFFTplan = fftw_plan_dft_c2r_2d(fSteps, fSteps, fFFTin, fFFTout, FFTW_ESTIMATE);
}
void TBulkTriVortexMLFieldCalc::CalculateGrid() const {
// SetParameters - method has to be called from the user before the calculation!!
if (fParam.size() < 3) {
cout << endl << "The SetParameters-method has to be called before B(x,y) can be calculated!" << endl;
return;
}
if (!fParam[0] || !fParam[1] || !fParam[2]) {
cout << endl << "The field, penetration depth and coherence length have to have finite values in order to calculate B(x,y)!" << endl;
return;
}
double field(fabs(fParam[0])), lambda(fabs(fParam[1])), xi(fabs(fParam[2]));
double Hc2(getHc2(xi));
const int NFFT(fSteps);
const int NFFT_2(fSteps/2);
const int NFFTsq(fSteps*fSteps);
// fill the field Fourier components in the matrix
// ... but first check that the field is not larger than Hc2 and that we are dealing with a type II SC
if ((field >= Hc2) || (lambda < xi/sqrt(2.0))) {
int m;
#pragma omp parallel for default(shared) private(m) schedule(dynamic)
for (m = 0; m < NFFTsq; m++) {
fFFTout[m] = field;
}
// Set the flag which shows that the calculation has been done
fGridExists = true;
return;
}
double latConstTr(sqrt(2.0*fluxQuantum/(field*sqrt3)));
double oneMb(1.0-field/Hc2);
double xisq_2_scaled(2.0/(3.0*oneMb)*pow(xi*PI/latConstTr,2.0)), lambdasq_scaled(4.0/(3.0*oneMb)*pow(lambda*PI/latConstTr,2.0));
// ... now fill in the Fourier components if everything was okay above
double Gsq, ll;
int k, l, lNFFT_2;
for (l = 0; l < NFFT_2; l += 2) {
lNFFT_2 = l*(NFFT_2 + 1);
ll = 3.0*static_cast<double>(l*l);
for (k = 0; k < NFFT_2; k += 2) {
Gsq = static_cast<double>(k*k) + ll;
fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][1] = 0.0;
fFFTin[lNFFT_2 + k + 1][0] = 0.0;
fFFTin[lNFFT_2 + k + 1][1] = 0.0;
}
k = NFFT_2;
Gsq = static_cast<double>(k*k) + ll;
fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][1] = 0.0;
}
for (l = NFFT_2; l < NFFT; l += 2) {
lNFFT_2 = l*(NFFT_2 + 1);
ll = 3.0*static_cast<double>((NFFT-l)*(NFFT-l));
for (k = 0; k < NFFT_2; k += 2) {
Gsq = static_cast<double>(k*k) + ll;
fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][1] = 0.0;
fFFTin[lNFFT_2 + k + 1][0] = 0.0;
fFFTin[lNFFT_2 + k + 1][1] = 0.0;
}
k = NFFT_2;
Gsq = static_cast<double>(k*k) + ll;
fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][1] = 0.0;
}
// intermediate rows
for (l = 1; l < NFFT_2; l += 2) {
lNFFT_2 = l*(NFFT_2 + 1);
ll = 3.0*static_cast<double>(l*l);
for (k = 0; k < NFFT_2; k += 2) {
Gsq = static_cast<double>((k + 1)*(k + 1)) + ll;
fFFTin[lNFFT_2 + k][0] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0;
fFFTin[lNFFT_2 + k + 1][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k + 1][1] = 0.0;
}
k = NFFT_2;
fFFTin[lNFFT_2 + k][0] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0;
}
for (l = NFFT_2 + 1; l < NFFT; l += 2) {
lNFFT_2 = l*(NFFT_2 + 1);
ll = 3.0*static_cast<double>((NFFT-l)*(NFFT-l));
for (k = 0; k < NFFT_2; k += 2) {
Gsq = static_cast<double>((k+1)*(k+1)) + ll;
fFFTin[lNFFT_2 + k][0] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0;
fFFTin[lNFFT_2 + k + 1][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k + 1][1] = 0.0;
}
k = NFFT_2;
fFFTin[lNFFT_2 + k][0] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0;
}
// Do the Fourier transform to get B(x,y)
fftw_execute(fFFTplan);
// Multiply by the applied field
#pragma omp parallel for default(shared) private(l) schedule(dynamic)
for (l = 0; l < NFFTsq; l++) {
fFFTout[l] *= field;
}
// Set the flag which shows that the calculation has been done
fGridExists = true;
return;
}
TBulkTriVortexAGLFieldCalc::TBulkTriVortexAGLFieldCalc(const string& wisdom, const unsigned int steps) {
fWisdom = wisdom;
if (steps % 2) {
fSteps = steps + 1;
} else {
fSteps = steps;
}
fParam.resize(3);
fGridExists = false;
int init_threads(fftw_init_threads());
if (init_threads)
fftw_plan_with_nthreads(2);
fFFTin = new fftw_complex[(fSteps/2 + 1) * fSteps];
fFFTout = new double[fSteps*fSteps];
// cout << "Check for the FFT plan..." << endl;
// Load wisdom from file if it exists and should be used
fUseWisdom = true;
int wisdomLoaded(0);
FILE *wordsOfWisdomR;
wordsOfWisdomR = fopen(fWisdom.c_str(), "r");
if (wordsOfWisdomR == NULL) {
fUseWisdom = false;
} else {
wisdomLoaded = fftw_import_wisdom_from_file(wordsOfWisdomR);
fclose(wordsOfWisdomR);
}
if (!wisdomLoaded) {
fUseWisdom = false;
}
// create the FFT plan
if (fUseWisdom)
fFFTplan = fftw_plan_dft_c2r_2d(fSteps, fSteps, fFFTin, fFFTout, FFTW_EXHAUSTIVE);
else
fFFTplan = fftw_plan_dft_c2r_2d(fSteps, fSteps, fFFTin, fFFTout, FFTW_ESTIMATE);
}
void TBulkTriVortexAGLFieldCalc::CalculateGrid() const {
// SetParameters - method has to be called from the user before the calculation!!
if (fParam.size() < 3) {
cout << endl << "The SetParameters-method has to be called before B(x,y) can be calculated!" << endl;
return;
}
if (!fParam[0] || !fParam[1] || !fParam[2]) {
cout << endl << "The field, penetration depth and coherence length have to have finite values in order to calculate B(x,y)!" << endl;
return;
}
double field(fabs(fParam[0])), lambda(fabs(fParam[1])), xi(fabs(fParam[2]));
double Hc2(getHc2(xi));
const int NFFT(fSteps);
const int NFFT_2(fSteps/2);
const int NFFTsq(fSteps*fSteps);
// fill the field Fourier components in the matrix
// ... but first check that the field is not larger than Hc2 and that we are dealing with a type II SC
if ((field >= Hc2) || (lambda < xi/sqrt(2.0))) {
int m;
#pragma omp parallel for default(shared) private(m) schedule(dynamic)
for (m = 0; m < NFFTsq; m++) {
fFFTout[m] = field;
}
// Set the flag which shows that the calculation has been done
fGridExists = true;
return;
}
double latConstTr(sqrt(2.0*fluxQuantum/(field*sqrt3)));
double b(field/Hc2);
double xisq_scaled(8.0/3.0*pow(xi*PI/latConstTr,2.0)*(1.0+pow(b,4.0))*(1.0-2.0*b*pow(1.0-b,2.0)));
double lambdasq_scaled(4.0/3.0*pow(lambda*PI/latConstTr,2.0));
// ... now fill in the Fourier components if everything was okay above
double Gsq, sqrtXiSqScGsq, ll;
int k, l, lNFFT_2;
for (l = 0; l < NFFT_2; l += 2) {
lNFFT_2 = l*(NFFT_2 + 1);
ll = 3.0*static_cast<double>(l*l);
for (k = 0; k < NFFT_2; k += 2) {
Gsq = static_cast<double>(k*k) + ll;
sqrtXiSqScGsq = ((k || l) ? sqrt(xisq_scaled*Gsq) : 1.0E-9);
fFFTin[lNFFT_2 + k][0] = sqrtXiSqScGsq*TMath::BesselK1(sqrtXiSqScGsq)/(lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][1] = 0.0;
fFFTin[lNFFT_2 + k + 1][0] = 0.0;
fFFTin[lNFFT_2 + k + 1][1] = 0.0;
}
k = NFFT_2;
Gsq = static_cast<double>(k*k) + ll;
sqrtXiSqScGsq = sqrt(xisq_scaled*Gsq);
fFFTin[lNFFT_2 + k][0] = sqrtXiSqScGsq*TMath::BesselK1(sqrtXiSqScGsq)/(lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][1] = 0.0;
}
for (l = NFFT_2; l < NFFT; l += 2) {
lNFFT_2 = l*(NFFT_2 + 1);
ll = 3.0*static_cast<double>((NFFT-l)*(NFFT-l));
for (k = 0; k < NFFT_2; k += 2) {
Gsq = static_cast<double>(k*k) + ll;
sqrtXiSqScGsq = sqrt(xisq_scaled*Gsq);
fFFTin[lNFFT_2 + k][0] = sqrtXiSqScGsq*TMath::BesselK1(sqrtXiSqScGsq)/(lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][1] = 0.0;
fFFTin[lNFFT_2 + k + 1][0] = 0.0;
fFFTin[lNFFT_2 + k + 1][1] = 0.0;
}
k = NFFT_2;
Gsq = static_cast<double>(k*k) + ll;
sqrtXiSqScGsq = sqrt(xisq_scaled*Gsq);
fFFTin[lNFFT_2 + k][0] = sqrtXiSqScGsq*TMath::BesselK1(sqrtXiSqScGsq)/(lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][1] = 0.0;
}
// intermediate rows
for (l = 1; l < NFFT_2; l += 2) {
lNFFT_2 = l*(NFFT_2 + 1);
ll = 3.0*static_cast<double>(l*l);
for (k = 0; k < NFFT_2; k += 2) {
Gsq = static_cast<double>((k + 1)*(k + 1)) + ll;
sqrtXiSqScGsq = sqrt(xisq_scaled*Gsq);
fFFTin[lNFFT_2 + k][0] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0;
fFFTin[lNFFT_2 + k + 1][0] = sqrtXiSqScGsq*TMath::BesselK1(sqrtXiSqScGsq)/(lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k + 1][1] = 0.0;
}
k = NFFT_2;
fFFTin[lNFFT_2 + k][0] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0;
}
for (l = NFFT_2 + 1; l < NFFT; l += 2) {
lNFFT_2 = l*(NFFT_2 + 1);
ll = 3.0*static_cast<double>((NFFT-l)*(NFFT-l));
for (k = 0; k < NFFT_2; k += 2) {
Gsq = static_cast<double>((k+1)*(k+1)) + ll;
sqrtXiSqScGsq = sqrt(xisq_scaled*Gsq);
fFFTin[lNFFT_2 + k][0] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0;
fFFTin[lNFFT_2 + k + 1][0] = sqrtXiSqScGsq*TMath::BesselK1(sqrtXiSqScGsq)/(lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k + 1][1] = 0.0;
}
k = NFFT_2;
fFFTin[lNFFT_2 + k][0] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0;
}
// Do the Fourier transform to get B(x,y)
fftw_execute(fFFTplan);
// Multiply by the applied field
#pragma omp parallel for default(shared) private(l) schedule(dynamic)
for (l = 0; l < NFFTsq; l++) {
fFFTout[l] *= field*(1.0-pow(b,4.0));
}
// Set the flag which shows that the calculation has been done
fGridExists = true;
return;
}
TBulkTriVortexNGLFieldCalc::TBulkTriVortexNGLFieldCalc(const string& wisdom, const unsigned int steps) TBulkTriVortexNGLFieldCalc::TBulkTriVortexNGLFieldCalc(const string& wisdom, const unsigned int steps)
: fLatticeConstant(0.0), fKappa(0.0), fSumAk(0.0), fSumOmegaSq(0.0), fSumSum(0.0) : fLatticeConstant(0.0), fKappa(0.0), fSumAk(0.0), fSumOmegaSq(0.0), fSumSum(0.0)
{ {
@ -419,13 +695,15 @@ TBulkTriVortexNGLFieldCalc::TBulkTriVortexNGLFieldCalc(const string& wisdom, con
// create the FFT plans // create the FFT plans
if (fUseWisdom) { if (fUseWisdom) {
fFFTplanAkToOmega = fftw_plan_dft_c2r_2d(fSteps, fSteps, fFFTin, fOmegaMatrix, FFTW_EXHAUSTIVE); // use the first plan from the base class here - it will be destroyed by the base class destructor
fFFTplan = fftw_plan_dft_c2r_2d(fSteps, fSteps, fFFTin, fOmegaMatrix, FFTW_EXHAUSTIVE);
fFFTplanBkToBandQ = fftw_plan_dft_2d(fSteps, fSteps, fBkMatrix, fBkMatrix, FFTW_BACKWARD, FFTW_EXHAUSTIVE); fFFTplanBkToBandQ = fftw_plan_dft_2d(fSteps, fSteps, fBkMatrix, fBkMatrix, FFTW_BACKWARD, FFTW_EXHAUSTIVE);
fFFTplanOmegaToAk = fftw_plan_dft_r2c_2d(fSteps, fSteps, fOmegaMatrix, fFFTin, FFTW_EXHAUSTIVE); fFFTplanOmegaToAk = fftw_plan_dft_r2c_2d(fSteps, fSteps, fOmegaMatrix, fFFTin, FFTW_EXHAUSTIVE);
fFFTplanOmegaToBk = fftw_plan_dft_2d(fSteps, fSteps, fBkMatrix, fBkMatrix, FFTW_FORWARD, FFTW_EXHAUSTIVE); fFFTplanOmegaToBk = fftw_plan_dft_2d(fSteps, fSteps, fBkMatrix, fBkMatrix, FFTW_FORWARD, FFTW_EXHAUSTIVE);
} }
else { else {
fFFTplanAkToOmega = fftw_plan_dft_c2r_2d(fSteps, fSteps, fFFTin, fOmegaMatrix, FFTW_ESTIMATE); // use the first plan from the base class here - it will be destroyed by the base class destructor
fFFTplan = fftw_plan_dft_c2r_2d(fSteps, fSteps, fFFTin, fOmegaMatrix, FFTW_ESTIMATE);
fFFTplanBkToBandQ = fftw_plan_dft_2d(fSteps, fSteps, fBkMatrix, fBkMatrix, FFTW_BACKWARD, FFTW_ESTIMATE); fFFTplanBkToBandQ = fftw_plan_dft_2d(fSteps, fSteps, fBkMatrix, fBkMatrix, FFTW_BACKWARD, FFTW_ESTIMATE);
fFFTplanOmegaToAk = fftw_plan_dft_r2c_2d(fSteps, fSteps, fOmegaMatrix, fFFTin, FFTW_ESTIMATE); fFFTplanOmegaToAk = fftw_plan_dft_r2c_2d(fSteps, fSteps, fOmegaMatrix, fFFTin, FFTW_ESTIMATE);
fFFTplanOmegaToBk = fftw_plan_dft_2d(fSteps, fSteps, fBkMatrix, fBkMatrix, FFTW_FORWARD, FFTW_ESTIMATE); fFFTplanOmegaToBk = fftw_plan_dft_2d(fSteps, fSteps, fBkMatrix, fBkMatrix, FFTW_FORWARD, FFTW_ESTIMATE);
@ -436,7 +714,6 @@ TBulkTriVortexNGLFieldCalc::~TBulkTriVortexNGLFieldCalc() {
// clean up // clean up
fftw_destroy_plan(fFFTplanAkToOmega);
fftw_destroy_plan(fFFTplanBkToBandQ); fftw_destroy_plan(fFFTplanBkToBandQ);
fftw_destroy_plan(fFFTplanOmegaToAk); fftw_destroy_plan(fFFTplanOmegaToAk);
fftw_destroy_plan(fFFTplanOmegaToBk); fftw_destroy_plan(fFFTplanOmegaToBk);
@ -973,7 +1250,7 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
// Do the Fourier transform to get omega(x,y) - Abrikosov // Do the Fourier transform to get omega(x,y) - Abrikosov
fftw_execute(fFFTplanAkToOmega); fftw_execute(fFFTplan);
#pragma omp parallel for default(shared) private(l) schedule(dynamic) #pragma omp parallel for default(shared) private(l) schedule(dynamic)
for (l = 0; l < NFFTsq; l++) { for (l = 0; l < NFFTsq; l++) {
@ -1059,7 +1336,7 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
// Do the Fourier transform to get omega(x,y) // Do the Fourier transform to get omega(x,y)
fftw_execute(fFFTplanAkToOmega); fftw_execute(fFFTplan);
#pragma omp parallel for default(shared) private(l) schedule(dynamic) #pragma omp parallel for default(shared) private(l) schedule(dynamic)
for (l = 0; l < NFFTsq; l++) { for (l = 0; l < NFFTsq; l++) {
@ -1125,7 +1402,7 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
// Do the Fourier transform to get omega(x,y) // Do the Fourier transform to get omega(x,y)
fftw_execute(fFFTplanAkToOmega); fftw_execute(fFFTplan);
#pragma omp parallel for default(shared) private(l) schedule(dynamic) #pragma omp parallel for default(shared) private(l) schedule(dynamic)
for (l = 0; l < NFFTsq; l++) { for (l = 0; l < NFFTsq; l++) {

316
src/external/TFitPofB-lib/classes/TVortex.cpp vendored
View File

@ -39,6 +39,8 @@ using namespace std;
#include "TFitPofBStartupHandler.h" #include "TFitPofBStartupHandler.h"
ClassImp(TBulkTriVortexLondon) ClassImp(TBulkTriVortexLondon)
ClassImp(TBulkTriVortexML)
ClassImp(TBulkTriVortexAGL)
ClassImp(TBulkTriVortexNGL) ClassImp(TBulkTriVortexNGL)
//------------------ //------------------
@ -58,6 +60,40 @@ TBulkTriVortexLondon::~TBulkTriVortexLondon() {
fParForPofT.clear(); fParForPofT.clear();
} }
//------------------
// Destructor of the TBulkTriVortexML class -- cleaning up
//------------------
TBulkTriVortexML::~TBulkTriVortexML() {
delete fPofT;
fPofT = 0;
delete fPofB;
fPofT = 0;
delete fVortex;
fVortex = 0;
fPar.clear();
fParForVortex.clear();
fParForPofB.clear();
fParForPofT.clear();
}
//------------------
// Destructor of the TBulkTriVortexAGL class -- cleaning up
//------------------
TBulkTriVortexAGL::~TBulkTriVortexAGL() {
delete fPofT;
fPofT = 0;
delete fPofB;
fPofT = 0;
delete fVortex;
fVortex = 0;
fPar.clear();
fParForVortex.clear();
fParForPofB.clear();
fParForPofT.clear();
}
//------------------ //------------------
// Destructor of the TBulkTriVortexNGL class -- cleaning up // Destructor of the TBulkTriVortexNGL class -- cleaning up
//------------------ //------------------
@ -136,7 +172,7 @@ TBulkTriVortexLondon::TBulkTriVortexLondon() : fCalcNeeded(true), fFirstCall(tru
//------------------ //------------------
// TBulkTriVortexLondon-Method that calls the procedures to create B(x,y), p(B) and P(t) // TBulkTriVortexLondon-Method that calls the procedures to create B(x,y), p(B) and P(t)
// It finally returns P(t) for a given t. // It finally returns P(t) for a given t.
// Parameters: all the parameters for the function to be fitted through TBulkTriVortexLondon (phase, appl.field, lambda, xi, [not implemented: bkg weight]) // Parameters: all the parameters for the function to be fitted through TBulkTriVortexLondon (phase, av.field, lambda, xi, [not implemented: bkg weight])
//------------------ //------------------
double TBulkTriVortexLondon::operator()(double t, const vector<double> &par) const { double TBulkTriVortexLondon::operator()(double t, const vector<double> &par) const {
@ -213,6 +249,284 @@ double TBulkTriVortexLondon::operator()(double t, const vector<double> &par) con
} }
//------------------
// Constructor of the TBulkTriVortexML class
// creates (a pointer to) the TPofTCalc object (with the FFT plan)
//------------------
TBulkTriVortexML::TBulkTriVortexML() : fCalcNeeded(true), fFirstCall(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;
}
fGridSteps = startupHandler->GetGridSteps();
fWisdom = startupHandler->GetWisdomFile();
fParForVortex.resize(3); // field, lambda, xi
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); // Bkg-Field
fParForPofB.push_back(0.005); // Bkg-width
fParForPofB.push_back(0.0); // Bkg-weight
TBulkTriVortexMLFieldCalc *x = new TBulkTriVortexMLFieldCalc(fWisdom, fGridSteps);
fVortex = x;
x = 0;
TPofBCalc *y = new TPofBCalc(fParForPofB);
fPofB = y;
y = 0;
TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT);
fPofT = z;
z = 0;
// clean up
if (saxParser) {
delete saxParser;
saxParser = 0;
}
if (startupHandler) {
delete startupHandler;
startupHandler = 0;
}
}
//------------------
// TBulkTriVortexML-Method that calls the procedures to create B(x,y), 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 TBulkTriVortexML (phase, av.field, lambda, xi, [not implemented: bkg weight])
//------------------
double TBulkTriVortexML::operator()(double t, const vector<double> &par) const {
assert(par.size() == 4 || par.size() == 5);
if(t<0.0)
return cos(par[0]*0.017453293);
// check if the function is called the first time and if yes, read in parameters
if(fFirstCall){
fPar = par;
for (unsigned int i(0); i < 3; i++) {
fParForVortex[i] = fPar[i+1];
}
fFirstCall = false;
}
// check if any parameter has changed
bool par_changed(false);
bool only_phase_changed(false);
for (unsigned int i(0); i<fPar.size(); i++) {
if( fPar[i]-par[i] ) {
fPar[i] = par[i];
par_changed = true;
if (i == 0) {
only_phase_changed = true;
} else {
only_phase_changed = false;
}
}
}
if (par_changed)
fCalcNeeded = true;
// if model parameters have changed, recalculate B(x,y), P(B) and P(t)
if (fCalcNeeded) {
fParForPofT[0] = par[0]; // phase
if(!only_phase_changed) {
// cout << " Parameters have changed, (re-)calculating p(B) and P(t) now..." << endl;
for (unsigned int i(0); i < 3; i++) {
fParForVortex[i] = par[i+1];
}
fParForPofB[2] = par[1]; // Bkg-Field
//fParForPofB[3] = 0.005; // Bkg-width (in principle zero)
fVortex->SetParameters(fParForVortex);
fVortex->CalculateGrid();
fPofB->UnsetPBExists();
fPofB->Calculate(fVortex, fParForPofB);
fPofT->DoFFT();
}/* else {
cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
}*/
fPofT->CalcPol(fParForPofT);
fCalcNeeded = false;
}
return fPofT->Eval(t);
}
//------------------
// Constructor of the TBulkTriVortexAGL class
// creates (a pointer to) the TPofTCalc object (with the FFT plan)
//------------------
TBulkTriVortexAGL::TBulkTriVortexAGL() : fCalcNeeded(true), fFirstCall(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;
}
fGridSteps = startupHandler->GetGridSteps();
fWisdom = startupHandler->GetWisdomFile();
fParForVortex.resize(3); // field, lambda, xi
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); // Bkg-Field
fParForPofB.push_back(0.005); // Bkg-width
fParForPofB.push_back(0.0); // Bkg-weight
TBulkTriVortexAGLFieldCalc *x = new TBulkTriVortexAGLFieldCalc(fWisdom, fGridSteps);
fVortex = x;
x = 0;
TPofBCalc *y = new TPofBCalc(fParForPofB);
fPofB = y;
y = 0;
TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT);
fPofT = z;
z = 0;
// clean up
if (saxParser) {
delete saxParser;
saxParser = 0;
}
if (startupHandler) {
delete startupHandler;
startupHandler = 0;
}
}
//------------------
// TBulkTriVortexAGL-Method that calls the procedures to create B(x,y), 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 TBulkTriVortexAGL (phase, av.field, lambda, xi, [not implemented: bkg weight])
//------------------
double TBulkTriVortexAGL::operator()(double t, const vector<double> &par) const {
assert(par.size() == 4 || par.size() == 5);
if(t<0.0)
return cos(par[0]*0.017453293);
// check if the function is called the first time and if yes, read in parameters
if(fFirstCall){
fPar = par;
for (unsigned int i(0); i < 3; i++) {
fParForVortex[i] = fPar[i+1];
}
fFirstCall = false;
}
// check if any parameter has changed
bool par_changed(false);
bool only_phase_changed(false);
for (unsigned int i(0); i<fPar.size(); i++) {
if( fPar[i]-par[i] ) {
fPar[i] = par[i];
par_changed = true;
if (i == 0) {
only_phase_changed = true;
} else {
only_phase_changed = false;
}
}
}
if (par_changed)
fCalcNeeded = true;
// if model parameters have changed, recalculate B(x,y), P(B) and P(t)
if (fCalcNeeded) {
fParForPofT[0] = par[0]; // phase
if(!only_phase_changed) {
// cout << " Parameters have changed, (re-)calculating p(B) and P(t) now..." << endl;
for (unsigned int i(0); i < 3; i++) {
fParForVortex[i] = par[i+1];
}
fParForPofB[2] = par[1]; // Bkg-Field
//fParForPofB[3] = 0.005; // Bkg-width (in principle zero)
fVortex->SetParameters(fParForVortex);
fVortex->CalculateGrid();
fPofB->UnsetPBExists();
fPofB->Calculate(fVortex, fParForPofB);
fPofT->DoFFT();
}/* else {
cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
}*/
fPofT->CalcPol(fParForPofT);
fCalcNeeded = false;
}
return fPofT->Eval(t);
}
//------------------ //------------------
// Constructor of the TBulkTriVortexNGL class // Constructor of the TBulkTriVortexNGL class
// creates (a pointer to) the TPofTCalc object (with the FFT plan) // creates (a pointer to) the TPofTCalc object (with the FFT plan)

35
src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h → src/external/TFitPofB-lib/include/TBulkTriVortexFieldCalc.h vendored
View File

@ -1,6 +1,6 @@
/*************************************************************************** /***************************************************************************
TBulkVortexFieldCalc.h TBulkTriVortexFieldCalc.h
Author: Bastian M. Wojek, Alexander Maisuradze Author: Bastian M. Wojek, Alexander Maisuradze
e-mail: bastian.wojek@psi.ch e-mail: bastian.wojek@psi.ch
@ -85,6 +85,36 @@ public:
}; };
//--------------------
// Class for triangular vortex lattice, modified London model
//--------------------
class TBulkTriVortexMLFieldCalc : public TBulkVortexFieldCalc {
public:
TBulkTriVortexMLFieldCalc(const string&, const unsigned int steps = 256);
~TBulkTriVortexMLFieldCalc() {}
void CalculateGrid() const;
};
//--------------------
// Class for triangular vortex lattice, analytical GL approximation
//--------------------
class TBulkTriVortexAGLFieldCalc : public TBulkVortexFieldCalc {
public:
TBulkTriVortexAGLFieldCalc(const string&, const unsigned int steps = 256);
~TBulkTriVortexAGLFieldCalc() {}
void CalculateGrid() const;
};
//-------------------- //--------------------
// Class for triangular vortex lattice, Minimisation of the GL free energy à la Brandt // Class for triangular vortex lattice, Minimisation of the GL free energy à la Brandt
//-------------------- //--------------------
@ -129,11 +159,10 @@ private:
mutable double fSumAk; mutable double fSumAk;
mutable double fSumOmegaSq; mutable double fSumOmegaSq;
mutable double fSumSum; mutable double fSumSum;
fftw_plan fFFTplanAkToOmega;
fftw_plan fFFTplanBkToBandQ; fftw_plan fFFTplanBkToBandQ;
fftw_plan fFFTplanOmegaToAk; fftw_plan fFFTplanOmegaToAk;
fftw_plan fFFTplanOmegaToBk; fftw_plan fFFTplanOmegaToBk;
}; };
#endif // _TBulkVortexFieldCalc_H_ #endif // _TBulkTriVortexFieldCalc_H_

2
src/external/TFitPofB-lib/include/TPofBCalc.h vendored
View File

@ -34,7 +34,7 @@
#include "TBofZCalc.h" #include "TBofZCalc.h"
#include "TTrimSPDataHandler.h" #include "TTrimSPDataHandler.h"
#include "TBulkVortexFieldCalc.h" #include "TBulkTriVortexFieldCalc.h"
#define gBar 0.0135538817 #define gBar 0.0135538817
#define pi 3.14159265358979323846 #define pi 3.14159265358979323846

48
src/external/TFitPofB-lib/include/TVortex.h vendored
View File

@ -59,6 +59,54 @@ private:
ClassDef(TBulkTriVortexLondon,1) ClassDef(TBulkTriVortexLondon,1)
}; };
class TBulkTriVortexML : public PUserFcnBase {
public:
TBulkTriVortexML();
~TBulkTriVortexML();
double operator()(double, const vector<double>&) const;
private:
mutable vector<double> fPar;
TBulkTriVortexMLFieldCalc *fVortex;
TPofBCalc *fPofB;
TPofTCalc *fPofT;
mutable bool fCalcNeeded;
mutable bool fFirstCall;
mutable vector<double> fParForVortex;
mutable vector<double> fParForPofB;
mutable vector<double> fParForPofT;
string fWisdom;
unsigned int fGridSteps;
ClassDef(TBulkTriVortexML,1)
};
class TBulkTriVortexAGL : public PUserFcnBase {
public:
TBulkTriVortexAGL();
~TBulkTriVortexAGL();
double operator()(double, const vector<double>&) const;
private:
mutable vector<double> fPar;
TBulkTriVortexAGLFieldCalc *fVortex;
TPofBCalc *fPofB;
TPofTCalc *fPofT;
mutable bool fCalcNeeded;
mutable bool fFirstCall;
mutable vector<double> fParForVortex;
mutable vector<double> fParForPofB;
mutable vector<double> fParForPofT;
string fWisdom;
unsigned int fGridSteps;
ClassDef(TBulkTriVortexAGL,1)
};
class TBulkTriVortexNGL : public PUserFcnBase { class TBulkTriVortexNGL : public PUserFcnBase {
public: public:

2
src/external/TFitPofB-lib/include/TVortexLinkDef.h vendored
View File

@ -37,6 +37,8 @@
#pragma link off all functions; #pragma link off all functions;
#pragma link C++ class TBulkTriVortexLondon+; #pragma link C++ class TBulkTriVortexLondon+;
#pragma link C++ class TBulkTriVortexML+;
#pragma link C++ class TBulkTriVortexAGL+;
#pragma link C++ class TBulkTriVortexNGL+; #pragma link C++ class TBulkTriVortexNGL+;
#endif //__CINT__ #endif //__CINT__

12
src/external/TFitPofB-lib/test/testVortex-v2.cpp vendored
View File

@ -19,18 +19,18 @@ int main(){
// parForVortex[2] = 4.0; //xi // parForVortex[2] = 4.0; //xi
vector<double> parForPofB; vector<double> parForPofB;
parForPofB.push_back(0.01); //dt parForPofB.push_back(0.005); //dt
parForPofB.push_back(0.1); //dB parForPofB.push_back(20.0); //dB
vector<double> parForPofT; vector<double> parForPofT;
parForPofT.push_back(0.0); //phase parForPofT.push_back(0.0); //phase
parForPofT.push_back(0.01); //dt parForPofT.push_back(0.005); //dt
parForPofT.push_back(0.1); //dB parForPofT.push_back(20.0); //dB
TBulkTriVortexLondonFieldCalc *vortexLattice = new TBulkTriVortexLondonFieldCalc("/home/l_wojek/analysis/WordsOfWisdom.dat", NFFT); TBulkTriVortexLondonFieldCalc *vortexLattice = new TBulkTriVortexLondonFieldCalc("/home/l_wojek/analysis/WordsOfWisdom.dat", NFFT);
parForVortex[0] = 10.0; //app.field parForVortex[0] = 3000.0; //app.field
parForVortex[1] = 200.0; //lambda parForVortex[1] = 100.0; //lambda
parForVortex[2] = 4.0; //xi parForVortex[2] = 4.0; //xi
vortexLattice->SetParameters(parForVortex); vortexLattice->SetParameters(parForVortex);