diff --git a/src/external/TFitPofB-lib/classes/Makefile.libTFitPofB b/src/external/TFitPofB-lib/classes/Makefile.libTFitPofB
index 9a815f09..ac87df73 100644
--- a/src/external/TFitPofB-lib/classes/Makefile.libTFitPofB
+++ b/src/external/TFitPofB-lib/classes/Makefile.libTFitPofB
@@ -37,16 +37,16 @@ endif
 
 # -- Linux
 ifeq ($(OS),LINUX)
-CXX           = g++-4.4.2
-CXXFLAGS      = -O3 -fopenmp -Wall -Wno-trigraphs -fPIC
+CXX           = g++
+CXXFLAGS      = -O3 -Wall -Wno-trigraphs -fPIC
 PMUSRPATH     = ../../../include
 MNPATH        = $(ROOTSYS)/include
 LOCALPATH     = ../include
 FFTW3PATH     = /usr/include
-INCLUDES      = -I$(PMUSRPATH) -I$(MNPATH) -I$(LOCALPATH) -I$(FFTW3PATH)
-LD            = g++-4.4.2
+INCLUDES      = -I$(LOCALPATH) -I$(PMUSRPATH) -I$(MNPATH) -I$(FFTW3PATH)
+LD            = g++
 LDFLAGS       = -O
-SOFLAGS       = -shared -fopenmp
+SOFLAGS       = -shared
 SHLIB         = libTFitPofB.so
 endif
 
@@ -58,7 +58,7 @@ PMUSRPATH     = ../../../include
 MNPATH        = $(ROOTSYS)/include
 LOCALPATH     = ../include
 FFTW3PATH     = /usr/include
-INCLUDES      = -I$(PMUSRPATH) -I$(MNPATH) -I$(LOCALPATH) -I$(FFTW3PATH)
+INCLUDES      = -I$(LOCALPATH) -I$(PMUSRPATH) -I$(MNPATH) -I$(FFTW3PATH)
 LD            = g++
 LDFLAGS       = -O -Wl,--enable-auto-import -Wl,--enable-runtime-pseudo-reloc
 SOFLAGS       = -shared -Wl,--export-all-symbols
@@ -74,7 +74,7 @@ MNPATH        = $(ROOTSYS)/include
 FINKPATH      = /sw/include
 LOCALPATH     = ../include
 FFTW3PATH     = /sw/include
-INCLUDES      = -I$(PMUSRPATH) -I$(MNPATH) -I$(LOCALPATH) -I$(FINKPATH) -I$(FFTW3PATH)
+INCLUDES      = -I$(LOCALPATH) -I$(PMUSRPATH) -I$(MNPATH) -I$(FINKPATH) -I$(FFTW3PATH)
 LD            = g++
 LDFLAGS       = -O -Xlinker -bind_at_load
 SOFLAGS       = -dynamiclib -flat_namespace -undefined suppress -Wl,-x
diff --git a/src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp b/src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp
index 1d5efd8f..d780e770 100644
--- a/src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp
+++ b/src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp
@@ -358,11 +358,23 @@ void TBulkTriVortexLondonFieldCalc::CalculateGrid() const {
 TBulkTriVortexNGLFieldCalc::TBulkTriVortexNGLFieldCalc(const string& wisdom, const unsigned int steps) : fLatticeConstant(0.0), fSumAk(0.0), fSumOmegaSq(0.0)
 {
   fWisdom = wisdom;
-  if (steps % 2) {
-    fSteps = steps + 1;
-  } else {
-    fSteps = steps;
+  switch (steps % 4) {
+    case 0:
+      fSteps = steps;
+      break;
+    case 1:
+      fSteps = steps + 3;
+      break;
+    case 2:
+      fSteps = steps + 2;
+      break;
+    case 3:
+      fSteps = steps + 1;
+      break;
+    default:
+      break;
   }
+
   fParam.resize(3);
   fGridExists = false;
 
@@ -377,9 +389,10 @@ TBulkTriVortexNGLFieldCalc::TBulkTriVortexNGLFieldCalc(const string& wisdom, con
   fAkMatrix = new fftw_complex[stepsSq];
   fBkMatrix = new fftw_complex[stepsSq];
   fRealSpaceMatrix = new fftw_complex[stepsSq];
-  fBMatrix = new double[stepsSq];
   fOmegaMatrix = new double[stepsSq];
   fOmegaSqMatrix = new double[stepsSq];
+  fOmegaDDiffXMatrix = new double[stepsSq];
+  fOmegaDDiffYMatrix = new double[stepsSq];
   fOmegaDiffXMatrix = new double[stepsSq];
   fOmegaDiffYMatrix = new double[stepsSq];
   fGMatrix = new double[stepsSq];
@@ -387,6 +400,7 @@ TBulkTriVortexNGLFieldCalc::TBulkTriVortexNGLFieldCalc(const string& wisdom, con
   fQyMatrix = new double[stepsSq];
   fQxMatrixA = new double[stepsSq];
   fQyMatrixA = new double[stepsSq];
+  fAbrikosovCheck = new double[stepsSq];
 
   fCheckAkConvergence = new double[fSteps];
   fCheckBkConvergence = new double[fSteps];
@@ -438,11 +452,12 @@ TBulkTriVortexNGLFieldCalc::~TBulkTriVortexNGLFieldCalc() {
   delete[] fAkMatrix; fAkMatrix = 0;
   delete[] fBkMatrix; fBkMatrix = 0;
   delete[] fRealSpaceMatrix; fRealSpaceMatrix = 0;
-  delete[] fBMatrix; fBMatrix = 0;
   delete[] fOmegaMatrix; fOmegaMatrix = 0;
   delete[] fOmegaSqMatrix; fOmegaSqMatrix = 0;
   delete[] fOmegaDiffXMatrix; fOmegaDiffXMatrix = 0;
   delete[] fOmegaDiffYMatrix; fOmegaDiffYMatrix = 0;
+  delete[] fOmegaDDiffXMatrix; fOmegaDiffXMatrix = 0;
+  delete[] fOmegaDDiffYMatrix; fOmegaDiffYMatrix = 0;
   delete[] fGMatrix; fGMatrix = 0;
   delete[] fQxMatrix; fQxMatrix = 0;
   delete[] fQyMatrix; fQyMatrix = 0;
@@ -450,22 +465,39 @@ TBulkTriVortexNGLFieldCalc::~TBulkTriVortexNGLFieldCalc() {
   delete[] fQyMatrixA; fQyMatrixA = 0;
   delete[] fCheckAkConvergence; fCheckAkConvergence = 0;
   delete[] fCheckBkConvergence; fCheckBkConvergence = 0;
+  delete[] fAbrikosovCheck; fAbrikosovCheck = 0;
 }
 
 void TBulkTriVortexNGLFieldCalc::CalculateIntermediateMatrices() const {
   const int NFFT(fSteps);
+  const int NFFT_2(fSteps/2);
   const int NFFTsq(fSteps*fSteps);
 
 // Derivatives of omega
 
-//  const double denomY(2.0*fLatticeConstant/static_cast<double>(NFFT));
-  const double denomY(2.0/static_cast<double>(NFFT));
+  const double denomY(2.0*fLatticeConstant/static_cast<double>(NFFT));
+//  const double denomY(2.0/static_cast<double>(NFFT));
 //  const double denomY(2.0);
+//  const double denomX(2.0);
   const double denomX(denomY*sqrt3);
   const double kappa(fParam[1]/fParam[2]);
   const double fourKappaSq(4.0*kappa*kappa);
+  
+  // Ensure that omega at the vortex-core positions is zero
+  fOmegaMatrix[0] = 0.0;
+  fOmegaMatrix[(NFFT+1)*NFFT_2] = 0.0;
+
   int i;
 
+  // Ensure that omega is NOT negative
+#pragma omp parallel for default(shared) private(i) schedule(dynamic)
+  for (i = 0; i < NFFTsq; i += 1) {
+    if (fOmegaMatrix[i] < 0.0) {
+      cout << "Omega negative for index " << i << endl;
+      fOmegaMatrix[i] = 0.0;
+    }
+  }
+
 #pragma omp parallel for default(shared) private(i) schedule(dynamic)
   for (i = 0; i < NFFTsq; i += 1) {
     if (!(i % NFFT)) { // first column
@@ -490,6 +522,12 @@ void TBulkTriVortexNGLFieldCalc::CalculateIntermediateMatrices() const {
     fOmegaDiffYMatrix[i] = (fOmegaMatrix[i+NFFT]-fOmegaMatrix[i-NFFT])/denomY;
   }
 
+  // Ensure that the derivatives at the vortex-core positions are zero
+  fOmegaDiffXMatrix[0] = 0.0;
+  fOmegaDiffXMatrix[(NFFT+1)*NFFT_2] = 0.0;
+  fOmegaDiffYMatrix[0] = 0.0;
+  fOmegaDiffYMatrix[(NFFT+1)*NFFT_2] = 0.0;
+
 // g-Matrix
 
 #pragma omp parallel for default(shared) private(i) schedule(dynamic)
@@ -501,6 +539,32 @@ void TBulkTriVortexNGLFieldCalc::CalculateIntermediateMatrices() const {
     }
   }
 
+// Laplacian for Abrikosov check
+
+#pragma omp parallel for default(shared) private(i) schedule(dynamic)
+  for (i = 0; i < NFFTsq; i += 1) {
+    if (!(i % NFFT)) { // first column
+      fOmegaDDiffXMatrix[i] = (fOmegaDiffXMatrix[i+1]-fOmegaDiffXMatrix[i+NFFT-1])/denomX;
+    } else if (!((i + 1) % NFFT)) { // last column
+      fOmegaDDiffXMatrix[i] = (fOmegaDiffXMatrix[i-NFFT+1]-fOmegaDiffXMatrix[i-1])/denomX;
+    } else {
+      fOmegaDDiffXMatrix[i] = (fOmegaDiffXMatrix[i+1]-fOmegaDiffXMatrix[i-1])/denomX;
+    }
+  }
+
+  for (i = 0; i < NFFT; i++) { // first row
+    fOmegaDDiffYMatrix[i] = (fOmegaDiffYMatrix[i+NFFT]-fOmegaDiffYMatrix[NFFTsq-NFFT+i])/denomY;
+  }
+
+  for (i = NFFTsq - NFFT; i < NFFTsq; i++) { // last row
+    fOmegaDDiffYMatrix[i] = (fOmegaDiffYMatrix[i-NFFTsq+NFFT]-fOmegaDiffYMatrix[i-NFFT])/denomY;
+  }
+
+#pragma omp parallel for default(shared) private(i) schedule(dynamic)
+  for (i = NFFT; i < NFFTsq - NFFT; i++) { // the rest
+    fOmegaDDiffYMatrix[i] = (fOmegaDiffYMatrix[i+NFFT]-fOmegaDiffYMatrix[i-NFFT])/denomY;
+  }
+
   return;
 }
 
@@ -779,8 +843,10 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
   }
 
   double field(fabs(fParam[0])), lambda(fabs(fParam[1])), xi(fabs(fParam[2]));
-  double Hc2(getHc2(xi)), kappa(lambda/xi), S(TWOPI/(kappa*field));
-  double scaledB(field/fluxQuantum*2.0*PI*xi*lambda);
+  double Hc2(getHc2(xi)), kappa(lambda/xi), Hc2_kappa(Hc2/kappa);
+  double scaledB(field/Hc2_kappa);
+
+  cout << "4pi/S = " << 2.0*kappa*scaledB << endl;
 
 //  fLatticeConstant = sqrt(2.0*fluxQuantum/(scaledB*sqrt3));
   fLatticeConstant = sqrt(2.0*TWOPI/(kappa*scaledB*sqrt3));
@@ -965,7 +1031,7 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
 
 #pragma omp parallel for default(shared) private(l) schedule(dynamic)
   for (l = 0; l < NFFT; l++) {
-    fCheckAkConvergence[l] = fAkMatrix[NFFT_2/2*NFFT + l][0];
+    fCheckAkConvergence[l] = fAkMatrix[l][0];
   }
 
   fSumAk = 0.0;
@@ -984,13 +1050,34 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
     fOmegaMatrix[l] = fSumAk - fRealSpaceMatrix[l][0];
   }
 
+  double sumOmega(0.0);
+  for (l = 0; l < NFFTsq; l++) {
+    sumOmega += fOmegaMatrix[l];
+  }
+  sumOmega /= static_cast<double>(NFFTsq);
+
+  cout << "sumOmega = " << sumOmega << endl;
+
   CalculateIntermediateMatrices();
 
+///////////// Check the Abrikosov solution
+
+  for (l = 0; l < NFFTsq; l++) {
+    if (!fOmegaMatrix[l]) {
+      fAbrikosovCheck[l] = 0.0;
+    } else {
+      fAbrikosovCheck[l] = (fOmegaDiffXMatrix[l]*fOmegaDiffXMatrix[l]+fOmegaDiffYMatrix[l]*fOmegaDiffYMatrix[l])/(fOmegaMatrix[l]*fOmegaMatrix[l]) - (fOmegaDDiffXMatrix[l] + fOmegaDDiffYMatrix[l])/fOmegaMatrix[l];
+    }
+  }
+
+
+  
+
   double denomQA;
 
 #pragma omp parallel for default(shared) private(l,denomQA) schedule(dynamic)
   for (l = 0; l < NFFTsq; l++) {
-    if (!fOmegaMatrix[l]) {
+    if (!fOmegaMatrix[l] || !l || (l == (NFFT+1)*NFFT_2)) {
       fQxMatrixA[l] = 0.0;
       fQyMatrixA[l] = 0.0;
     } else {
@@ -1006,14 +1093,18 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
 
 #pragma omp parallel for default(shared) private(l) schedule(dynamic)
   for (l = 0; l < NFFTsq; l++) {
-    fBMatrix[l] = scaledB;
+    fFFTout[l] = scaledB;
   }
 
   bool akConverged(false), bkConverged(false), akInitiallyConverged(false), firstBkCalculation(true);
   double coeff1, coeff2;
 
+  int count(0);
+
   while (!akConverged || !bkConverged) {
 
+//  if (count == 3) break;
+
 #pragma omp parallel for default(shared) private(l) schedule(dynamic)
     for (l = 0; l < NFFTsq; l++) {
       fOmegaSqMatrix[l] = fOmegaMatrix[l]*fOmegaMatrix[l];
@@ -1021,7 +1112,7 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
 
 #pragma omp parallel for default(shared) private(l) schedule(dynamic)
     for (l = 0; l < NFFTsq; l++) {
-      fRealSpaceMatrix[l][0] = fOmegaSqMatrix[l] + fOmegaMatrix[l]*(fQxMatrix[l]*fQxMatrix[l] + fQyMatrix[l]*fQyMatrix[l] - 2.0) + fGMatrix[l];
+      fRealSpaceMatrix[l][0] = fOmegaMatrix[l]*fOmegaMatrix[l] - 2.0*fOmegaMatrix[l] + fOmegaMatrix[l]*(fQxMatrix[l]*fQxMatrix[l] + fQyMatrix[l]*fQyMatrix[l]) + (!fOmegaMatrix[l] ? 0.0 : (fOmegaDiffXMatrix[l]*fOmegaDiffXMatrix[l]+fOmegaDiffYMatrix[l]*fOmegaDiffYMatrix[l])/(4.0*kappa*kappa*fOmegaMatrix[l]));
       fRealSpaceMatrix[l][1] = 0.0;
     }
 
@@ -1048,7 +1139,15 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
       fOmegaMatrix[l] = fSumAk - fRealSpaceMatrix[l][0];
     }
 
-    CalculateIntermediateMatrices();
+  CalculateIntermediateMatrices();
+
+  sumOmega = 0.0;
+  for (l = 0; l < NFFTsq; l++) {
+    sumOmega += fOmegaMatrix[l];
+  }
+  sumOmega /= static_cast<double>(NFFTsq);
+  cout << "sumOmega = " << sumOmega << endl;
+
 //
 #pragma omp parallel for default(shared) private(l) schedule(dynamic)
     for (l = 0; l < NFFTsq; l++) {
@@ -1062,8 +1161,10 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
 
     double coeffSum(0.0);
     for (l = 0; l < NFFTsq; l++) {
-      coeffSum += fOmegaMatrix[l]*(1.0 - (fQxMatrix[l]*fQxMatrix[l] + fQyMatrix[l]*fQyMatrix[l])) - fGMatrix[l];
+      coeffSum += fOmegaMatrix[l] - fOmegaMatrix[l]*(fQxMatrix[l]*fQxMatrix[l] + fQyMatrix[l]*fQyMatrix[l]) - (!fOmegaMatrix[l] ? 0.0 : (fOmegaDiffXMatrix[l]*fOmegaDiffXMatrix[l]+fOmegaDiffYMatrix[l]*fOmegaDiffYMatrix[l])/(4.0*kappa*kappa*fOmegaMatrix[l]));
     }
+    
+    cout << "coeffSum = " << coeffSum << ", fSumOmegaSq = " << fSumOmegaSq << endl;
 
     for (l = 0; l < NFFTsq; l++) {
       fAkMatrix[l][0] *= coeffSum/fSumOmegaSq;
@@ -1072,7 +1173,8 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
     akConverged = true;
 
     for (l = 0; l < NFFT; l++) {
-      if (fabs(fCheckAkConvergence[l] - fAkMatrix[NFFT_2/2*NFFT + l][0]) > 1.0E-12) {
+      if ((fabs(fCheckAkConvergence[l] - fAkMatrix[l][0])/fAkMatrix[l][0] > 0.025) && (fabs(fAkMatrix[l][0]) > 1.0E-4)) {
+        cout << "old: " << fCheckAkConvergence[l] << ", new: " << fAkMatrix[l][0] << endl;
         akConverged = false;
         break;
       }
@@ -1081,7 +1183,7 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
     if (!akConverged) {
 #pragma omp parallel for default(shared) private(l) schedule(dynamic)
       for (l = 0; l < NFFT; l++) {
-        fCheckAkConvergence[l] = fAkMatrix[NFFT_2/2*NFFT + l][0];
+        fCheckAkConvergence[l] = fAkMatrix[l][0];
       }
     } else {
       akInitiallyConverged = true;
@@ -1106,12 +1208,20 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
     }
 
     CalculateIntermediateMatrices();
+    
+  sumOmega = 0.0;
+  for (l = 0; l < NFFTsq; l++) {
+    sumOmega += fOmegaMatrix[l];
+  }
+  sumOmega /= static_cast<double>(NFFTsq);
+   cout << "sumOmega = " << sumOmega << endl;
 
+//break;
   if (akInitiallyConverged) {// break;
 
 #pragma omp parallel for default(shared) private(l) schedule(dynamic)
     for (l = 0; l < NFFTsq; l++) {
-      fRealSpaceMatrix[l][0] = (fOmegaMatrix[l]-fSumAk)*fBMatrix[l] + fSumAk*scaledB - fQxMatrix[l]*fOmegaDiffYMatrix[l] + fQyMatrix[l]*fOmegaDiffXMatrix[l];
+      fRealSpaceMatrix[l][0] = fOmegaMatrix[l]*fFFTout[l] + fSumAk*(scaledB - fFFTout[l]) - fQxMatrix[l]*fOmegaDiffYMatrix[l] + fQyMatrix[l]*fOmegaDiffXMatrix[l];
       fRealSpaceMatrix[l][1] = 0.0;
     }
 
@@ -1120,7 +1230,7 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
     if (firstBkCalculation) {
 #pragma omp parallel for default(shared) private(l) schedule(dynamic)
       for (l = 0; l < NFFT; l++) {
-        fCheckBkConvergence[l] = fBkMatrix[NFFT_2/2*NFFT + l][0];
+        fCheckBkConvergence[l] = fBkMatrix[l][0];
       }
       firstBkCalculation = false;
       akConverged = false;
@@ -1134,13 +1244,24 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
 
 #pragma omp parallel for default(shared) private(l) schedule(dynamic)
     for (l = 0; l < NFFTsq; l++) {
-      fBMatrix[l] = scaledB + fRealSpaceMatrix[l][0];
+      fFFTout[l] = scaledB + fRealSpaceMatrix[l][0];
     }
 
+  // Ensure that the field is NOT negative
+#pragma omp parallel for default(shared) private(l) schedule(dynamic)
+  for (l = 0; l < NFFTsq; l++) {
+    if (fFFTout[l] < 0.0)
+      fFFTout[l] = 0.0;
+  }
+//  break;
+
+    
+
     bkConverged = true;
 
     for (l = 0; l < NFFT; l++) {
-      if (fabs(fCheckBkConvergence[l] - fBkMatrix[NFFT_2/2*NFFT + l][0]) > 1.0E-12) {
+      if ((fabs(fCheckBkConvergence[l] - fBkMatrix[l][0])/fBkMatrix[l][0] > 0.025) && (fabs(fBkMatrix[l][0]) > 1.0E-4)) {
+        cout << "old: " << fCheckBkConvergence[l] << ", new: " << fBkMatrix[l][0] << endl;
         bkConverged = false;
         break;
       }
@@ -1151,10 +1272,13 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
     if (!bkConverged) {
 #pragma omp parallel for default(shared) private(l) schedule(dynamic)
       for (l = 0; l < NFFT; l++) {
-        fCheckBkConvergence[l] = fBkMatrix[NFFT_2/2*NFFT + l][0];
+        fCheckBkConvergence[l] = fBkMatrix[l][0];
       }
-    } else if (akConverged) {
+    } else if (count == 1) {
       break;
+    } else {
+      bkConverged = false;
+      count++;
     }
 
 // I need a copy of Bk... store it to Ak since already tooooooooooooooooooo much memory is allocated
@@ -1191,6 +1315,11 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const {
   }
   }
 
+#pragma omp parallel for default(shared) private(l) schedule(dynamic)
+  for (l = 0; l < NFFTsq; l++) {
+    fFFTout[l] *= Hc2_kappa;
+  }
+
   // Set the flag which shows that the calculation has been done
 
   fGridExists = true;
diff --git a/src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h b/src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h
index 82b3d9df..01f4f9c0 100644
--- a/src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h
+++ b/src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h
@@ -101,11 +101,12 @@ public:
   fftw_complex* GetAkMatrix() const {return fAkMatrix;}
   fftw_complex* GetBkMatrix() const {return fBkMatrix;}
   double* GetOmegaMatrix() const {return fOmegaMatrix;}
-  double* GetBMatrix() const {return fBMatrix;}
+  double* GetBMatrix() const {return fFFTout;}
   double* GetOmegaDiffXMatrix() const {return fOmegaDiffXMatrix;}
   double* GetOmegaDiffYMatrix() const {return fOmegaDiffYMatrix;}
   double* GetQxMatrix() const {return fQxMatrix;}
   double* GetQyMatrix() const {return fQyMatrix;}
+  double* GetAbrikosovCheck() const {return fAbrikosovCheck;}
 
 private:
 
@@ -117,11 +118,13 @@ private:
   fftw_complex *fAkMatrix;
   fftw_complex *fBkMatrix;
   mutable fftw_complex *fRealSpaceMatrix;
-  mutable double *fBMatrix;
+  mutable double *fAbrikosovCheck;
   mutable double *fOmegaMatrix;
   mutable double *fOmegaSqMatrix;
   mutable double *fOmegaDiffXMatrix;
   mutable double *fOmegaDiffYMatrix;
+  mutable double *fOmegaDDiffXMatrix;
+  mutable double *fOmegaDDiffYMatrix;
   mutable double *fQxMatrix;
   mutable double *fQyMatrix;
   mutable double *fQxMatrixA;