Added a few more functions to libFitPofB

src/external/TFitPofB-lib/classes/TBulkTriVortexFieldCalc.cpp

@@ -1985,17 +1985,11 @@ void TBulkAnisotropicTriVortexLondonFieldCalc::CalculateGrid() const {
   double field(fabs(fParam[0])), lambdaX(fabs(fParam[1])), lambdaY(fabs(fParam[2])), xiX(fabs(fParam[3])), xiY(fabs(fParam[4]));
   double Hc2(fluxQuantum/(TWOPI*xiX*xiY));
 
-//  double latConstTr(sqrt(2.0*fluxQuantum/(field*sqrt3)));
+  double coeff(PI*PI*field/(sqrt3*fluxQuantum));
-  double xiXsq_2_scaled(field/fluxQuantum*pow(xiX*PI,2.0)/(sqrt3*lambdaY)*lambdaX);
+  double xiXsq_2_scaled(3.0*coeff*xiX*xiX*lambdaX/lambdaY);
-  double xiYsq_2_scaled(field/fluxQuantum*pow(xiY*PI,2.0)*lambdaY/lambdaX*sqrt3);
+  double xiYsq_2_scaled(coeff*xiY*xiY*lambdaY/lambdaX);
-  double lambdaXsq_scaled(2.0*field/fluxQuantum*PI*PI*sqrt3*lambdaY*lambdaX);
+  double lambdaXsq_scaled(2.0*coeff*lambdaX*lambdaY);
-  double lambdaYsq_scaled(2.0*field/fluxQuantum*PI*PI/sqrt3*lambdaX*lambdaY);
+  double lambdaYsq_scaled(3.0*lambdaXsq_scaled);
-/*
-  double xiXsq_2_scaled(field/fluxQuantum*pow(xiX*PI,2.0)*sqrt3*sqrt(lambdaY/lambdaX));
-  double xiYsq_2_scaled(field/fluxQuantum*pow(xiY*PI,2.0)*sqrt(lambdaX/(lambdaY*3.0)));
-  double lambdaXsq_scaled(2.0*field/fluxQuantum*PI*PI/sqrt3*pow(lambdaX,2.5)/sqrt(lambdaY));
-  double lambdaYsq_scaled(2.0*field/fluxQuantum*PI*PI*sqrt3/sqrt(lambdaX)*pow(lambdaY,2.5));
-*/
 
   const int NFFT(fSteps);
   const int NFFT_2(fSteps/2);
@@ -2026,14 +2020,14 @@ void TBulkAnisotropicTriVortexLondonFieldCalc::CalculateGrid() const {
     ll = static_cast<double>(l*l);
     for (k = 0; k < NFFT_2; k += 2) {
       kk = static_cast<double>(k*k);
-      fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*ll + xiYsq_2_scaled*kk))/(1.0+lambdaXsq_scaled*kk+lambdaYsq_scaled*ll);
+      fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
       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;
     kk = static_cast<double>(k*k);
-    fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*ll + xiYsq_2_scaled*kk))/(1.0+lambdaXsq_scaled*kk+lambdaYsq_scaled*ll);
+    fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
     fFFTin[lNFFT_2 + k][1] = 0.0;
   }
 
@@ -2043,14 +2037,14 @@ void TBulkAnisotropicTriVortexLondonFieldCalc::CalculateGrid() const {
     ll = static_cast<double>((NFFT-l)*(NFFT-l));
     for (k = 0; k < NFFT_2; k += 2) {
       kk = static_cast<double>(k*k);
-      fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*ll + xiYsq_2_scaled*kk))/(1.0+lambdaXsq_scaled*kk+lambdaYsq_scaled*ll);
+      fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
       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;
     kk = static_cast<double>(k*k);
-    fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*ll + xiYsq_2_scaled*kk))/(1.0+lambdaXsq_scaled*kk+lambdaYsq_scaled*ll);
+    fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
     fFFTin[lNFFT_2 + k][1] = 0.0;
   }
 
@@ -2063,7 +2057,7 @@ void TBulkAnisotropicTriVortexLondonFieldCalc::CalculateGrid() const {
       kk = static_cast<double>((k + 1)*(k + 1));
       fFFTin[lNFFT_2 + k][0] = 0.0;
       fFFTin[lNFFT_2 + k][1] = 0.0;
-      fFFTin[lNFFT_2 + k + 1][0] = exp(-(xiXsq_2_scaled*ll + xiYsq_2_scaled*kk))/(1.0+lambdaXsq_scaled*kk+lambdaYsq_scaled*ll);
+      fFFTin[lNFFT_2 + k + 1][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
       fFFTin[lNFFT_2 + k + 1][1] = 0.0;
     }
     k = NFFT_2;
@@ -2078,7 +2072,7 @@ void TBulkAnisotropicTriVortexLondonFieldCalc::CalculateGrid() const {
       kk = static_cast<double>((k+1)*(k+1));
       fFFTin[lNFFT_2 + k][0] = 0.0;
       fFFTin[lNFFT_2 + k][1] = 0.0;
-      fFFTin[lNFFT_2 + k + 1][0] = exp(-(xiXsq_2_scaled*ll + xiYsq_2_scaled*kk))/(1.0+lambdaXsq_scaled*kk+lambdaYsq_scaled*ll);
+      fFFTin[lNFFT_2 + k + 1][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
       fFFTin[lNFFT_2 + k + 1][1] = 0.0;
     }
     k = NFFT_2;
@@ -2105,3 +2099,375 @@ void TBulkAnisotropicTriVortexLondonFieldCalc::CalculateGrid() const {
 
 }
 
+TBulkAnisotropicTriVortexMLFieldCalc::TBulkAnisotropicTriVortexMLFieldCalc(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 TBulkAnisotropicTriVortexMLFieldCalc::CalculateGrid() const {
+  // SetParameters - method has to be called from the user before the calculation!!
+  if (fParam.size() < 5) {
+    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] || !fParam[3] || !fParam[4]) {
+    cout << endl << "The field, penetration depths and coherence lengths have to have finite values in order to calculate B(x,y)!" \
+         << endl;
+    return;
+  }
+
+  double field(fabs(fParam[0])), lambdaX(fabs(fParam[1])), lambdaY(fabs(fParam[2])), xiX(fabs(fParam[3])), xiY(fabs(fParam[4]));
+  double Hc2(fluxQuantum/(TWOPI*xiX*xiY));
+  double oneMb(1.0-field/Hc2);
+
+  double coeff(PI*PI*field/(sqrt3*fluxQuantum*oneMb));
+  double xiXsq_2_scaled(3.0*coeff*xiX*xiX*lambdaX/lambdaY);
+  double xiYsq_2_scaled(coeff*xiY*xiY*lambdaY/lambdaX);
+  double lambdaXsq_scaled(2.0*coeff*lambdaX*lambdaY);
+  double lambdaYsq_scaled(3.0*lambdaXsq_scaled);
+
+  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) || (sqrt(lambdaX*lambdaY) < sqrt(xiX*xiY)/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;
+  }
+
+  // ... now fill in the Fourier components if everything was okay above
+  double kk, ll;
+  int k, l, lNFFT_2;
+
+  for (l = 0; l < NFFT_2; l += 2) {
+    lNFFT_2 = l*(NFFT_2 + 1);
+    ll = static_cast<double>(l*l);
+    for (k = 0; k < NFFT_2; k += 2) {
+      kk = static_cast<double>(k*k);
+      fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+      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;
+    kk = static_cast<double>(k*k);
+    fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+    fFFTin[lNFFT_2 + k][1] = 0.0;
+  }
+
+
+  for (l = NFFT_2; l < NFFT; l += 2) {
+    lNFFT_2 = l*(NFFT_2 + 1);
+    ll = static_cast<double>((NFFT-l)*(NFFT-l));
+    for (k = 0; k < NFFT_2; k += 2) {
+      kk = static_cast<double>(k*k);
+      fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+      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;
+    kk = static_cast<double>(k*k);
+    fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+    fFFTin[lNFFT_2 + k][1] = 0.0;
+  }
+
+  // intermediate rows
+
+  for (l = 1; l < NFFT_2; l += 2) {
+    lNFFT_2 = l*(NFFT_2 + 1);
+    ll = static_cast<double>(l*l);
+    for (k = 0; k < NFFT_2; k += 2) {
+      kk = static_cast<double>((k + 1)*(k + 1));
+      fFFTin[lNFFT_2 + k][0] = 0.0;
+      fFFTin[lNFFT_2 + k][1] = 0.0;
+      fFFTin[lNFFT_2 + k + 1][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+      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 = static_cast<double>((NFFT-l)*(NFFT-l));
+    for (k = 0; k < NFFT_2; k += 2) {
+      kk = static_cast<double>((k+1)*(k+1));
+      fFFTin[lNFFT_2 + k][0] = 0.0;
+      fFFTin[lNFFT_2 + k][1] = 0.0;
+      fFFTin[lNFFT_2 + k + 1][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+      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;
+
+}
+
+TBulkAnisotropicTriVortexAGLFieldCalc::TBulkAnisotropicTriVortexAGLFieldCalc(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 TBulkAnisotropicTriVortexAGLFieldCalc::CalculateGrid() const {
+  // SetParameters - method has to be called from the user before the calculation!!
+  if (fParam.size() < 5) {
+    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] || !fParam[3] || !fParam[4]) {
+    cout << endl << "The field, penetration depths and coherence lengths have to have finite values in order to calculate B(x,y)!" \
+         << endl;
+    return;
+  }
+
+  double field(fabs(fParam[0])), lambdaX(fabs(fParam[1])), lambdaY(fabs(fParam[2])), xiX(fabs(fParam[3])), xiY(fabs(fParam[4]));
+  double Hc2(fluxQuantum/(TWOPI*xiX*xiY));
+
+  double b(field/Hc2);
+  double coeff1(1.0 - pow(b,4.0));
+  double coeff2(2.0*(1.0 + pow(b,4.0))*(1.0 - 2.0*b*pow(1.0-b,2.0)));
+
+  double coeff(2.0*PI*PI*field/(sqrt3*fluxQuantum));
+  double xiXsq_scaled(3.0*coeff*xiX*xiX*lambdaX/lambdaY*coeff2);
+  double xiYsq_scaled(coeff*xiY*xiY*lambdaY/lambdaX*coeff2);
+  double lambdaXsq_scaled(coeff*lambdaX*lambdaY);
+  double lambdaYsq_scaled(3.0*lambdaXsq_scaled);
+
+  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) || (sqrt(lambdaX*lambdaY) < sqrt(xiX*xiY)/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;
+  }
+
+  // ... now fill in the Fourier components if everything was okay above
+  double kk, ll, u;
+  int k, l, lNFFT_2;
+
+  for (l = 0; l < NFFT_2; l += 2) {
+    lNFFT_2 = l*(NFFT_2 + 1);
+    ll = static_cast<double>(l*l);
+    for (k = 0; k < NFFT_2; k += 2) {
+      kk = static_cast<double>(k*k);
+      u = sqrt(xiXsq_scaled*kk + xiYsq_scaled*ll);
+      fFFTin[lNFFT_2 + k][0] = ((!k && !l) ? 1.0 : coeff1*u*TMath::BesselK1(u)/(lambdaXsq_scaled*ll+lambdaYsq_scaled*kk));
+      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;
+    kk = static_cast<double>(k*k);
+    u = sqrt(xiXsq_scaled*kk + xiYsq_scaled*ll);
+    fFFTin[lNFFT_2 + k][0] = coeff1*u*TMath::BesselK1(u)/(lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+    fFFTin[lNFFT_2 + k][1] = 0.0;
+  }
+
+
+  for (l = NFFT_2; l < NFFT; l += 2) {
+    lNFFT_2 = l*(NFFT_2 + 1);
+    ll = static_cast<double>((NFFT-l)*(NFFT-l));
+    for (k = 0; k < NFFT_2; k += 2) {
+      kk = static_cast<double>(k*k);
+      u = sqrt(xiXsq_scaled*kk + xiYsq_scaled*ll);
+      fFFTin[lNFFT_2 + k][0] = coeff1*u*TMath::BesselK1(u)/(lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+      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;
+    kk = static_cast<double>(k*k);
+    u = sqrt(xiXsq_scaled*kk + xiYsq_scaled*ll);
+    fFFTin[lNFFT_2 + k][0] = coeff1*u*TMath::BesselK1(u)/(lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+    fFFTin[lNFFT_2 + k][1] = 0.0;
+  }
+
+  // intermediate rows
+
+  for (l = 1; l < NFFT_2; l += 2) {
+    lNFFT_2 = l*(NFFT_2 + 1);
+    ll = static_cast<double>(l*l);
+    for (k = 0; k < NFFT_2; k += 2) {
+      kk = static_cast<double>((k + 1)*(k + 1));
+      u = sqrt(xiXsq_scaled*kk + xiYsq_scaled*ll);
+      fFFTin[lNFFT_2 + k][0] = 0.0;
+      fFFTin[lNFFT_2 + k][1] = 0.0;
+      fFFTin[lNFFT_2 + k + 1][0] = coeff1*u*TMath::BesselK1(u)/(lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+      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 = static_cast<double>((NFFT-l)*(NFFT-l));
+    for (k = 0; k < NFFT_2; k += 2) {
+      kk = static_cast<double>((k+1)*(k+1));
+      u = sqrt(xiXsq_scaled*kk + xiYsq_scaled*ll);
+      fFFTin[lNFFT_2 + k][0] = 0.0;
+      fFFTin[lNFFT_2 + k][1] = 0.0;
+      fFFTin[lNFFT_2 + k + 1][0] = coeff1*u*TMath::BesselK1(u)/(lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+      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;
+
+}

src/external/TFitPofB-lib/classes/TVortex.cpp

@@ -901,6 +901,7 @@ void TBulkAnisotropicTriVortexLondonGlobal::Calc(const vector<double> &par) cons
   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){
@@ -929,8 +930,10 @@ void TBulkAnisotropicTriVortexLondonGlobal::Calc(const vector<double> &par) cons
     }
   }
 
-  if (par_changed)
+  if (par_changed) {
     fCalcNeeded = true;
+    fValid = false;
+  }
 
   // if model parameters have changed, recalculate B(x,y), P(B) and P(t)
 
@@ -962,11 +965,11 @@ void TBulkAnisotropicTriVortexLondonGlobal::Calc(const vector<double> &par) cons
     fPofT->CalcPol(fParForPofT);
 
     fCalcNeeded = false;
+    fValid = true;
   }
 /*
   return fPofT->Eval(t);
 */
-  fValid = true;
   return;
 }
 
@@ -1057,3 +1060,481 @@ double TBulkAnisotropicTriVortexLondon::operator()(double t, const std::vector<d
 
   return fGlobalUserFcn->GetPolarizationPointer()->Eval(t);
 }
+
+//------------------
+// Destructor of the TBulkAnisotropicTriVortexMLGlobal class -- cleaning up
+//------------------
+
+TBulkAnisotropicTriVortexMLGlobal::~TBulkAnisotropicTriVortexMLGlobal() {
+    delete fPofT;
+    fPofT = 0;
+    delete fPofB;
+    fPofT = 0;
+    delete fVortex;
+    fVortex = 0;
+    fPar.clear();
+    fParForVortex.clear();
+    fParForPofB.clear();
+    fParForPofT.clear();
+}
+
+//------------------
+// Constructor of the TBulkAnisotropicTriVortexMLGlobal class
+// creates (a pointer to) the TPofTCalc object (with the FFT plan)
+//------------------
+
+TBulkAnisotropicTriVortexMLGlobal::TBulkAnisotropicTriVortexMLGlobal() : 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
+      cerr << endl << "**ERROR** reading/parsing TFitPofB_startup.xml failed." \
+           << endl << "**ERROR** Please make sure that the file exists in the local directory and it is set up correctly!" \
+           << endl;
+      assert(false);
+    }
+
+    fGridSteps = startupHandler->GetGridSteps();
+    fWisdom = startupHandler->GetWisdomFile();
+
+    fParForVortex.resize(5); // field, lambdaX, lambdaY, xiX, xiY
+
+    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
+
+    fVortex = new TBulkAnisotropicTriVortexMLFieldCalc(fWisdom, fGridSteps);
+
+    fPofB = new TPofBCalc(fParForPofB);
+
+    fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT);
+
+    // clean up
+    if (saxParser) {
+      delete saxParser;
+      saxParser = 0;
+    }
+    if (startupHandler) {
+      delete startupHandler;
+      startupHandler = 0;
+    }
+}
+
+void TBulkAnisotropicTriVortexMLGlobal::Calc(const vector<double> &par) const {
+
+  assert(par.size() == 6);
+/*
+  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 < 5; 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;
+    fValid = false;
+  }
+
+  // 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 < 5; 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;
+    fValid = true;
+  }
+/*
+  return fPofT->Eval(t);
+*/
+  return;
+}
+
+//------------------------------------------------------------------------------------
+/**
+ * <p> Constructor.
+ */
+TBulkAnisotropicTriVortexML::TBulkAnisotropicTriVortexML()
+{
+  fValid = false;
+  fInvokedGlobal = false;
+  fIdxGlobal = -1;
+  fGlobalUserFcn = 0;
+}
+
+//------------------------------------------------------------------------------------
+/**
+ * <p> Destructor.
+ */
+TBulkAnisotropicTriVortexML::~TBulkAnisotropicTriVortexML()
+{
+  if ((fGlobalUserFcn != 0) && fInvokedGlobal) {
+    delete fGlobalUserFcn;
+    fGlobalUserFcn = 0;
+  }
+}
+
+//------------------------------------------------------------------------------------
+/**
+ * <p> Used to invoke/retrieve the proper global user function
+ *
+ * \param globalPart reference to the global user function object vector
+ * \param idx global user function index within the theory tree
+ */
+void TBulkAnisotropicTriVortexML::SetGlobalPart(vector<void *> &globalPart, UInt_t idx)
+{
+  fIdxGlobal = static_cast<Int_t>(idx);
+
+  if ((Int_t)globalPart.size() <= fIdxGlobal) { // global user function not present, invoke it
+    fGlobalUserFcn = new TBulkAnisotropicTriVortexMLGlobal();
+    if (fGlobalUserFcn == 0) { // global user function object couldn't be invoked -> error
+      fValid = false;
+      cerr << endl << ">> TBulkAnisotropicTriVortexML::SetGlobalPart(): **ERROR** Couldn't invoke global user function object, sorry ..." << endl;
+    } else {  // global user function object could be invoked -> resize to global user function vector and keep the pointer to the corresponding object
+      globalPart.resize(fIdxGlobal+1);
+      globalPart[fIdxGlobal] = dynamic_cast<TBulkAnisotropicTriVortexMLGlobal*>(fGlobalUserFcn);
+      fValid = true;
+      fInvokedGlobal = true;
+    }
+  } else { // global user function already present hence just retrieve a pointer to it
+    fValid = true;
+    fGlobalUserFcn = (TBulkAnisotropicTriVortexMLGlobal*)globalPart[fIdxGlobal];
+  }
+}
+
+//------------------------------------------------------------------------------------
+/**
+ * <p> Used to check if the user function is OK.
+ *
+ * <b>return:</b> true if both, the user function and the global user function object are valid
+ */
+bool TBulkAnisotropicTriVortexML::GlobalPartIsValid() const
+{
+  return (fValid && fGlobalUserFcn->IsValid());
+}
+
+//------------------------------------------------------------------------------------
+/**
+ * <p> function operator
+ *
+ * <b>return:</b> function value
+ *
+ * \param t time in \f$(\mu\mathrm{s})\f$, or x-axis value for non-muSR fit
+ * \param param parameter vector
+ */
+double TBulkAnisotropicTriVortexML::operator()(double t, const std::vector<double> &param) const
+{
+  // expected parameters: phase, field, lambdaX, lambdaY, xiX, xiY
+
+  assert(param.size() == 6);
+  assert(fGlobalUserFcn);
+
+  if(t<0.0)
+    return cos(param[0]*0.017453293);
+
+  // call the global user function object
+  fGlobalUserFcn->Calc(param);
+
+  return fGlobalUserFcn->GetPolarizationPointer()->Eval(t);
+}
+
+//------------------
+// Destructor of the TBulkAnisotropicTriVortexAGLGlobal class -- cleaning up
+//------------------
+
+TBulkAnisotropicTriVortexAGLGlobal::~TBulkAnisotropicTriVortexAGLGlobal() {
+    delete fPofT;
+    fPofT = 0;
+    delete fPofB;
+    fPofT = 0;
+    delete fVortex;
+    fVortex = 0;
+    fPar.clear();
+    fParForVortex.clear();
+    fParForPofB.clear();
+    fParForPofT.clear();
+}
+
+//------------------
+// Constructor of the TBulkAnisotropicTriVortexAGLGlobal class
+// creates (a pointer to) the TPofTCalc object (with the FFT plan)
+//------------------
+
+TBulkAnisotropicTriVortexAGLGlobal::TBulkAnisotropicTriVortexAGLGlobal() : 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
+      cerr << endl << "**ERROR** reading/parsing TFitPofB_startup.xml failed." \
+           << endl << "**ERROR** Please make sure that the file exists in the local directory and it is set up correctly!" \
+           << endl;
+      assert(false);
+    }
+
+    fGridSteps = startupHandler->GetGridSteps();
+    fWisdom = startupHandler->GetWisdomFile();
+
+    fParForVortex.resize(5); // field, lambdaX, lambdaY, xiX, xiY
+
+    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
+
+    fVortex = new TBulkAnisotropicTriVortexAGLFieldCalc(fWisdom, fGridSteps);
+
+    fPofB = new TPofBCalc(fParForPofB);
+
+    fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT);
+
+    // clean up
+    if (saxParser) {
+      delete saxParser;
+      saxParser = 0;
+    }
+    if (startupHandler) {
+      delete startupHandler;
+      startupHandler = 0;
+    }
+}
+
+void TBulkAnisotropicTriVortexAGLGlobal::Calc(const vector<double> &par) const {
+
+  assert(par.size() == 6);
+/*
+  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 < 5; 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;
+    fValid = false;
+  }
+
+  // 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 < 5; 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;
+    fValid = true;
+  }
+/*
+  return fPofT->Eval(t);
+*/
+  return;
+}
+
+//------------------------------------------------------------------------------------
+/**
+ * <p> Constructor.
+ */
+TBulkAnisotropicTriVortexAGL::TBulkAnisotropicTriVortexAGL()
+{
+  fValid = false;
+  fInvokedGlobal = false;
+  fIdxGlobal = -1;
+  fGlobalUserFcn = 0;
+}
+
+//------------------------------------------------------------------------------------
+/**
+ * <p> Destructor.
+ */
+TBulkAnisotropicTriVortexAGL::~TBulkAnisotropicTriVortexAGL()
+{
+  if ((fGlobalUserFcn != 0) && fInvokedGlobal) {
+    delete fGlobalUserFcn;
+    fGlobalUserFcn = 0;
+  }
+}
+
+//------------------------------------------------------------------------------------
+/**
+ * <p> Used to invoke/retrieve the proper global user function
+ *
+ * \param globalPart reference to the global user function object vector
+ * \param idx global user function index within the theory tree
+ */
+void TBulkAnisotropicTriVortexAGL::SetGlobalPart(vector<void *> &globalPart, UInt_t idx)
+{
+  fIdxGlobal = static_cast<Int_t>(idx);
+
+  if ((Int_t)globalPart.size() <= fIdxGlobal) { // global user function not present, invoke it
+    fGlobalUserFcn = new TBulkAnisotropicTriVortexAGLGlobal();
+    if (fGlobalUserFcn == 0) { // global user function object couldn't be invoked -> error
+      fValid = false;
+      cerr << endl << ">> TBulkAnisotropicTriVortexAGL::SetGlobalPart(): **ERROR** Couldn't invoke global user function object, sorry ..." << endl;
+    } else {  // global user function object could be invoked -> resize to global user function vector and keep the pointer to the corresponding object
+      globalPart.resize(fIdxGlobal+1);
+      globalPart[fIdxGlobal] = dynamic_cast<TBulkAnisotropicTriVortexAGLGlobal*>(fGlobalUserFcn);
+      fValid = true;
+      fInvokedGlobal = true;
+    }
+  } else { // global user function already present hence just retrieve a pointer to it
+    fValid = true;
+    fGlobalUserFcn = (TBulkAnisotropicTriVortexAGLGlobal*)globalPart[fIdxGlobal];
+  }
+}
+
+//------------------------------------------------------------------------------------
+/**
+ * <p> Used to check if the user function is OK.
+ *
+ * <b>return:</b> true if both, the user function and the global user function object are valid
+ */
+bool TBulkAnisotropicTriVortexAGL::GlobalPartIsValid() const
+{
+  return (fValid && fGlobalUserFcn->IsValid());
+}
+
+//------------------------------------------------------------------------------------
+/**
+ * <p> function operator
+ *
+ * <b>return:</b> function value
+ *
+ * \param t time in \f$(\mu\mathrm{s})\f$, or x-axis value for non-muSR fit
+ * \param param parameter vector
+ */
+double TBulkAnisotropicTriVortexAGL::operator()(double t, const std::vector<double> &param) const
+{
+  // expected parameters: phase, field, lambdaX, lambdaY, xiX, xiY
+
+  assert(param.size() == 6);
+  assert(fGlobalUserFcn);
+
+  if(t<0.0)
+    return cos(param[0]*0.017453293);
+
+  // call the global user function object
+  fGlobalUserFcn->Calc(param);
+
+  return fGlobalUserFcn->GetPolarizationPointer()->Eval(t);
+}

src/external/TFitPofB-lib/include/TBulkTriVortexFieldCalc.h

@@ -87,7 +87,7 @@ public:
 };
 
 /**
- * <p>Class for the calculation of the spatial field distribution B(x,y) within a 2D "triangular vortex lattice"
+ * <p>Class for the calculation of the spatial field distribution B(x,y) within a 2D triangular vortex lattice
  * using the London model with a Gaussian cutoff for an anisotropic superconductor with the field applied along one of the principal axes 
  */
 class TBulkAnisotropicTriVortexLondonFieldCalc : public TBulkVortexFieldCalc {
@@ -134,6 +134,23 @@ public:
 
 };
 
+/**
+ * <p>Class for the calculation of the spatial field distribution B(x,y) within a 2D triangular vortex lattice
+ * using the modified London model with a Gaussian cutoff for an anisotropic superconductor with the field applied 
+ * along one of the principal axes 
+ */
+class TBulkAnisotropicTriVortexMLFieldCalc : public TBulkVortexFieldCalc {
+
+public:
+
+  TBulkAnisotropicTriVortexMLFieldCalc(const string&, const unsigned int steps = 256);
+  ~TBulkAnisotropicTriVortexMLFieldCalc() {}
+
+  void CalculateGrid() const;
+  bool IsTriangular() const {return true;}
+
+};
+
 /**
  * <p>Class for the calculation of the spatial field distribution B(x,y) within a 2D triangular vortex lattice
  * using the analytical Ginzburg-Landau approximation
@@ -150,6 +167,24 @@ public:
 
 };
 
+/**
+ * <p>Class for the calculation of the spatial field distribution B(x,y) within a 2D triangular vortex lattice
+ * using the analytical Ginzburg-Landau approximation for an anisotropic superconductor with the field applied 
+ * along one of the principal axes
+ */
+class TBulkAnisotropicTriVortexAGLFieldCalc : public TBulkVortexFieldCalc {
+
+public:
+
+  TBulkAnisotropicTriVortexAGLFieldCalc(const string&, const unsigned int steps = 256);
+  ~TBulkAnisotropicTriVortexAGLFieldCalc() {}
+
+  void CalculateGrid() const;
+  bool IsTriangular() const {return true;}
+
+};
+
+
 /**
  * <p>Class for the calculation of the spatial field distribution B(x,y) within a 2D triangular vortex lattice
  * using an iterative minimization of the Ginzburg-Landau free energy after E.H. Brandt

src/external/TFitPofB-lib/include/TVortex.h

@@ -202,7 +202,7 @@ private:
 
 /**
  * <p>Class implementing the global part of the musrfit user function interface for calculating muon spin depolarization functions
- * originating from the spatial field distribution B(x,y) within a 2D "triangular vortex lattice"
+ * originating from the spatial field distribution B(x,y) within a 2D triangular vortex lattice
  * calculated using the London model with a Gaussian cutoff for an anisotropic superconductor with the field applied along
  * one of the principal axes
  */
@@ -215,17 +215,13 @@ public:
 
   bool IsValid() { return fValid; }
 
-  // the following function will check if something needs to be calculated, which 
-  // is the case if param != fPrevParam
   void Calc(const vector<double>&) const;
 
-  // this routine will return the calculated values, e.g. B(z,E) for TMyFunction::operator()()
-  // (...) here sketches only that some parameters are likley to be fed
   const TPofTCalc* GetPolarizationPointer() const { return fPofT; }
 
 private:
-  mutable bool fValid;
+  mutable bool fValid; ///< tag indicating if the global part has been calculated
-  mutable vector<double> fPar;
+  mutable vector<double> fPar; ///< parameter vector
   TBulkAnisotropicTriVortexLondonFieldCalc *fVortex; ///< spatial field distribution B(x,y) within the vortex lattice
   TPofBCalc *fPofB; ///< static field distribution P(B)
   TPofTCalc *fPofT; ///< muon spin polarization p(t)
@@ -268,4 +264,132 @@ private:
   ClassDef(TBulkAnisotropicTriVortexLondon,1)
 };
 
+/**
+ * <p>Class implementing the global part of the musrfit user function interface for calculating muon spin depolarization functions
+ * originating from the spatial field distribution B(x,y) within a 2D triangular vortex lattice
+ * calculated using the modified London model with a Gaussian cutoff for an anisotropic superconductor with the field applied along
+ * one of the principal axes
+ */
+class TBulkAnisotropicTriVortexMLGlobal {
+
+public:
+  // default constructor and destructor
+  TBulkAnisotropicTriVortexMLGlobal();
+  ~TBulkAnisotropicTriVortexMLGlobal();
+
+  bool IsValid() { return fValid; }
+
+  void Calc(const vector<double>&) const;
+
+  const TPofTCalc* GetPolarizationPointer() const { return fPofT; }
+
+private:
+  mutable bool fValid; ///< tag indicating if the global part has been calculated
+  mutable vector<double> fPar; ///< parameter vector
+  TBulkAnisotropicTriVortexMLFieldCalc *fVortex; ///< spatial field distribution B(x,y) within the vortex lattice
+  TPofBCalc *fPofB; ///< static field distribution P(B)
+  TPofTCalc *fPofT; ///< muon spin polarization p(t)
+  mutable bool fCalcNeeded; ///< tag needed to avoid unnecessary calculations if the core parameters were unchanged
+  mutable bool fFirstCall; ///< tag for checking if the function operator is called the first time
+  mutable vector<double> fParForVortex; ///< parameters for the calculation of B(x,y)
+  mutable vector<double> fParForPofB; ///< parameters for the calculation of P(B)
+  mutable vector<double> fParForPofT; ///< parameters for the calculation of p(t)
+  string fWisdom; ///< file name of the FFTW wisdom file
+  unsigned int fGridSteps; ///< number of points in x- and y-direction for which B(x,y) is calculated
+
+  // definition of the class for the ROOT-dictionary
+  ClassDef(TBulkAnisotropicTriVortexMLGlobal,1)
+};
+
+/**
+ * <p>Class implementing the musrfit user function interface for calculating muon spin depolarization functions
+ * originating from the spatial field distribution B(x,y) within a 2D "triangular vortex lattice"
+ * calculated using the modified London model with a Gaussian cutoff for an anisotropic superconductor with the field applied along
+ * one of the principal axes
+ */
+class TBulkAnisotropicTriVortexML : public PUserFcnBase {
+
+public:
+  TBulkAnisotropicTriVortexML();
+  ~TBulkAnisotropicTriVortexML();
+
+  virtual bool NeedGlobalPart() const { return true; }
+  virtual void SetGlobalPart(vector<void *> &globalPart, unsigned int idx);
+  virtual bool GlobalPartIsValid() const;
+
+  double operator()(double, const vector<double>&) const;
+
+private:
+  bool fValid;
+  bool fInvokedGlobal;
+  int  fIdxGlobal;
+  TBulkAnisotropicTriVortexMLGlobal *fGlobalUserFcn;
+
+  ClassDef(TBulkAnisotropicTriVortexML,1)
+};
+
+/**
+ * <p>Class implementing the global part of the musrfit user function interface for calculating muon spin depolarization functions
+ * originating from the spatial field distribution B(x,y) within a 2D triangular vortex lattice
+ * calculated using the analytical Ginzburg-Landau model for an anisotropic superconductor with the field applied along
+ * one of the principal axes
+ */
+class TBulkAnisotropicTriVortexAGLGlobal {
+
+public:
+  // default constructor and destructor
+  TBulkAnisotropicTriVortexAGLGlobal();
+  ~TBulkAnisotropicTriVortexAGLGlobal();
+
+  bool IsValid() { return fValid; }
+
+  void Calc(const vector<double>&) const;
+
+  const TPofTCalc* GetPolarizationPointer() const { return fPofT; }
+
+private:
+  mutable bool fValid; ///< tag indicating if the global part has been calculated
+  mutable vector<double> fPar; ///< parameter vector
+  TBulkAnisotropicTriVortexAGLFieldCalc *fVortex; ///< spatial field distribution B(x,y) within the vortex lattice
+  TPofBCalc *fPofB; ///< static field distribution P(B)
+  TPofTCalc *fPofT; ///< muon spin polarization p(t)
+  mutable bool fCalcNeeded; ///< tag needed to avoid unnecessary calculations if the core parameters were unchanged
+  mutable bool fFirstCall; ///< tag for checking if the function operator is called the first time
+  mutable vector<double> fParForVortex; ///< parameters for the calculation of B(x,y)
+  mutable vector<double> fParForPofB; ///< parameters for the calculation of P(B)
+  mutable vector<double> fParForPofT; ///< parameters for the calculation of p(t)
+  string fWisdom; ///< file name of the FFTW wisdom file
+  unsigned int fGridSteps; ///< number of points in x- and y-direction for which B(x,y) is calculated
+
+  // definition of the class for the ROOT-dictionary
+  ClassDef(TBulkAnisotropicTriVortexAGLGlobal,1)
+};
+
+/**
+ * <p>Class implementing the musrfit user function interface for calculating muon spin depolarization functions
+ * originating from the spatial field distribution B(x,y) within a 2D "triangular vortex lattice"
+ * calculated using the analytical Ginzburg-Landau model for an anisotropic superconductor with the field applied along
+ * one of the principal axes
+ */
+class TBulkAnisotropicTriVortexAGL : public PUserFcnBase {
+
+public:
+  TBulkAnisotropicTriVortexAGL();
+  ~TBulkAnisotropicTriVortexAGL();
+
+  virtual bool NeedGlobalPart() const { return true; }
+  virtual void SetGlobalPart(vector<void *> &globalPart, unsigned int idx);
+  virtual bool GlobalPartIsValid() const;
+
+  double operator()(double, const vector<double>&) const;
+
+private:
+  bool fValid;
+  bool fInvokedGlobal;
+  int  fIdxGlobal;
+  TBulkAnisotropicTriVortexAGLGlobal *fGlobalUserFcn;
+
+  ClassDef(TBulkAnisotropicTriVortexAGL,1)
+};
+
 #endif //_TVortex_H_

src/external/TFitPofB-lib/include/TVortexLinkDef.h

@@ -43,6 +43,10 @@
 #pragma link C++ class TBulkTriVortexNGL+;
 #pragma link C++ class TBulkAnisotropicTriVortexLondon+;
 #pragma link C++ class TBulkAnisotropicTriVortexLondonGlobal+;
+#pragma link C++ class TBulkAnisotropicTriVortexML+;
+#pragma link C++ class TBulkAnisotropicTriVortexMLGlobal+;
+#pragma link C++ class TBulkAnisotropicTriVortexAGL+;
+#pragma link C++ class TBulkAnisotropicTriVortexAGLGlobal+;
 
 #endif //__CINT__
 // root dictionary stuff --------------------------------------------------