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

Tried to fix the ASCII export from musrview in the case of a Fourier-power-difference (some more tests needed); additionally minor changes to the BMWlibs

Browse Source
This commit is contained in:
Bastian M. Wojek
2011-04-10 16:27:36 +00:00
parent ae83e1a296
commit 2c514a881c
13 changed files with 492 additions and 39 deletions
Download Patch File Download Diff File Expand all files Collapse all files

197
src/external/libFitPofB/classes/TBulkTriVortexFieldCalc.cpp vendored
View File

@ -804,6 +804,195 @@ void TBulkTriVortexAGLFieldCalc::CalculateGrid() const {
TBulkTriVortexAGLIIFieldCalc::TBulkTriVortexAGLIIFieldCalc(const string& wisdom, const unsigned int steps) {
fWisdom = wisdom;
if (steps % 2) {
fSteps = steps + 1;
} else {
fSteps = steps;
}
fParam.resize(3);
fGridExists = false;
#ifdef HAVE_LIBFFTW3_THREADS
int init_threads(fftw_init_threads());
if (init_threads) {
#ifdef HAVE_GOMP
fftw_plan_with_nthreads(omp_get_num_procs());
#else
fftw_plan_with_nthreads(2);
#endif /* HAVE_GOMP */
}
#endif /* HAVE_LIBFFTW3_THREADS */
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 TBulkTriVortexAGLIIFieldCalc::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 vortex-core radius have to have finite values in order to calculate B(x,y)!" << endl;
return;
}
double field(fabs(fParam[0])), lambda(fabs(fParam[1])), xiV(fabs(fParam[2]));
double Hc2(getHc2(xiV)); // use the vortex-core radius for Hc2-calculation (which is wrong since xi_GL should be used instead: one would therefore need one more parameter)
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 < xiV/sqrt(2.0))) {
int m;
#ifdef HAVE_GOMP
#pragma omp parallel for default(shared) private(m) schedule(dynamic)
#endif
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 fInf(1.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, sqrtfInfSqPlusLsqGsq, 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;
sqrtfInfSqPlusLsqGsq = sqrt(fInf*fInf + lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][0] = ((!k && !l) ? 1.0 : \
fInf*TMath::BesselK1(xiV/lambda*sqrtfInfSqPlusLsqGsq)/(sqrtfInfSqPlusLsqGsq*TMath::BesselK1(xiV/lambda*fInf)));
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;
sqrtfInfSqPlusLsqGsq = sqrt(fInf*fInf + lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][0] = fInf*TMath::BesselK1(xiV/lambda*sqrtfInfSqPlusLsqGsq)/(sqrtfInfSqPlusLsqGsq*TMath::BesselK1(xiV/lambda*fInf));
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;
sqrtfInfSqPlusLsqGsq = sqrt(fInf*fInf + lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][0] = fInf*TMath::BesselK1(xiV/lambda*sqrtfInfSqPlusLsqGsq)/(sqrtfInfSqPlusLsqGsq*TMath::BesselK1(xiV/lambda*fInf));
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;
sqrtfInfSqPlusLsqGsq = sqrt(fInf*fInf + lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][0] = fInf*TMath::BesselK1(xiV/lambda*sqrtfInfSqPlusLsqGsq)/(sqrtfInfSqPlusLsqGsq*TMath::BesselK1(xiV/lambda*fInf));
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;
sqrtfInfSqPlusLsqGsq = sqrt(fInf*fInf + lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][0] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0;
fFFTin[lNFFT_2 + k + 1][0] = fInf*TMath::BesselK1(xiV/lambda*sqrtfInfSqPlusLsqGsq)/(sqrtfInfSqPlusLsqGsq*TMath::BesselK1(xiV/lambda*fInf));
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;
sqrtfInfSqPlusLsqGsq = sqrt(fInf*fInf + lambdasq_scaled*Gsq);
fFFTin[lNFFT_2 + k][0] = 0.0;
fFFTin[lNFFT_2 + k][1] = 0.0;
fFFTin[lNFFT_2 + k + 1][0] = fInf*TMath::BesselK1(xiV/lambda*sqrtfInfSqPlusLsqGsq)/(sqrtfInfSqPlusLsqGsq*TMath::BesselK1(xiV/lambda*fInf));
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
#ifdef HAVE_GOMP
#pragma omp parallel for default(shared) private(l) schedule(dynamic)
#endif
for (l = 0; l < NFFTsq; ++l) {
fFFTout[l] *= field;
}
// Set the flag which shows that the calculation has been done
fGridExists = true;
return;
}
TBulkTriVortexNGLFieldCalc::TBulkTriVortexNGLFieldCalc(const string& wisdom, const unsigned int steps)
: fLatticeConstant(0.0), fKappa(0.0), fSumAk(0.0), fSumOmegaSq(0.0), fSumSum(0.0)
{
@ -1728,11 +1917,11 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
for (l = 0; l < NFFT; l++) {
if (fFFTin[l][0]){
if (((fabs(fFFTin[l][0]) > 1.0E-6) && (fabs(fCheckAkConvergence[l] - fFFTin[l][0])/fFFTin[l][0] > 1.0E-6)) || \
if (((fabs(fFFTin[l][0]) > 1.0E-6) && (fabs(fCheckAkConvergence[l] - fFFTin[l][0])/fFFTin[l][0] > 1.0E-3)) || \
(fCheckAkConvergence[l]/fFFTin[l][0] < 0.0)) {
//cout << "old: " << fCheckAkConvergence[l] << ", new: " << fFFTin[l][0] << endl;
cout << "old: " << fCheckAkConvergence[l] << ", new: " << fFFTin[l][0] << endl;
akConverged = false;
//cout << "index = " << l << endl;
cout << "index = " << l << endl;
break;
}
}
@ -1808,7 +1997,7 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
for (l = 0; l < NFFT; l++) {
if (fBkMatrix[l][0]) {
if (((fabs(fBkMatrix[l][0]) > 1.0E-6) && (fabs(fCheckBkConvergence[l] - fBkMatrix[l][0])/fabs(fBkMatrix[l][0]) > 1.0E-6)) || \
if (((fabs(fBkMatrix[l][0]) > 1.0E-6) && (fabs(fCheckBkConvergence[l] - fBkMatrix[l][0])/fabs(fBkMatrix[l][0]) > 1.0E-3)) || \
(fCheckBkConvergence[l]/fBkMatrix[l][0] < 0.0)) {
// cout << "old: " << fCheckBkConvergence[l] << ", new: " << fBkMatrix[l][0] << endl;
bkConverged = false;