diff --git a/src/external/TFitPofB-lib/classes/Makefile.TFitPofB b/src/external/TFitPofB-lib/classes/Makefile.TFitPofB
index b927d1c2..4904adcb 100644
--- a/src/external/TFitPofB-lib/classes/Makefile.TFitPofB
+++ b/src/external/TFitPofB-lib/classes/Makefile.TFitPofB
@@ -12,7 +12,7 @@ OBJS += TTrimSPDataHandler.o
 OBJS += TBofZCalc.o
 OBJS += TPofBCalc.o
 OBJS += TPofTCalc.o
-OBJS += TUserLondon.o
+OBJS += TLondon1D.o
 
 SHLIB = libTFitPofB.so
 
diff --git a/src/external/TFitPofB-lib/classes/TLondon1D.cpp b/src/external/TFitPofB-lib/classes/TLondon1D.cpp
new file mode 100644
index 00000000..4ff1d0fa
--- /dev/null
+++ b/src/external/TFitPofB-lib/classes/TLondon1D.cpp
@@ -0,0 +1,426 @@
+/***************************************************************************
+
+  TLondon1D.cpp
+
+  Author: Bastian M. Wojek
+  e-mail: bastian.wojek@psi.ch
+
+  2008/05/27
+
+***************************************************************************/
+
+#include "TLondon1D.h"
+#include <iostream>
+using namespace std;
+
+//------------------
+// Destructor of the TLondon1D class -- cleaning up
+//------------------
+
+TLondon1D::~TLondon1D() {
+    fPar.clear();
+    delete fImpProfile;
+    fImpProfile = 0;
+    delete fPofT;
+    fPofT = 0;
+}
+
+//------------------
+// Constructor of the TLondon1D1L class -- reading available implantation profiles and
+// creates (a pointer to) the TPofTCalc object (with the FFT plan)
+//------------------
+
+TLondon1D1L::TLondon1D1L(const vector<unsigned int> &parNo, const vector<double> &par) {
+
+  for(unsigned int i(0); i<parNo.size(); i++) {
+    fPar.push_back(par[parNo[i]-1]);
+  }
+
+//  string rge_path("/home/l_wojek/nt/wojek/g/Bastian/ImplantationDepth/YBCO_PBCO-");
+//  string energy_arr[] = {"02_1", "02_5", "03_5", "05_0", "07_5", "10_0", "12_5", "15_0", "17_5", "19_0", "20_0", "22_5", "25_0"};
+
+  string rge_path("/home/l_wojek/TrimSP/AuYBCO_2005/AuYBCO-500000-");
+  string energy_arr[] = {"04_6", "09_6", "14_6", "18_6", "21_6", "24_6", "28_7"};
+
+  vector<string> energy_vec(energy_arr, energy_arr+(sizeof(energy_arr)/sizeof(energy_arr[0])));
+
+  vector<double> par_for_PofT;
+
+  for (unsigned int i(0); i<3; i++)
+    par_for_PofT.push_back(fPar[i]);
+
+  TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
+  fImpProfile = x;
+  x = 0;
+  delete x;
+
+  TPofTCalc *y = new TPofTCalc(par_for_PofT);
+  fPofT = y;
+  y = 0;
+  delete y;
+
+}
+
+//------------------
+// TLondon1D1L-Method that calls the procedures to create B(z), p(B) and P(t)
+// It finally returns P(t) for a given t.
+// Parameters: all the parameters for the function to be fitted through TLondon1D1L
+//------------------
+
+double TLondon1D1L::Eval(double t, const vector<double> &par) const {
+
+  // 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 == 1 || i == 2 || i == 3) {
+        cout << "You are varying dt, dB or E! These parameters have to be fixed! Quitting..." << endl; 
+        exit(-1);
+      } else if (i == 0) {
+        only_phase_changed = true;
+      } else {
+        only_phase_changed = false;
+      }
+    }
+  }
+
+  if (par_changed)
+    fCalcNeeded = true;
+
+/* DEBUGGING CODE COMMENTED -- quite a mess... sorry*/
+
+  // if model parameters have changed, recalculate B(z), P(B) and P(t)
+
+  if (fCalcNeeded) {
+
+    vector<double> par_for_PofT;
+
+//    cout << "par_for_PofT: ";
+
+    for (unsigned int i(0); i<3; i++) {
+      par_for_PofT.push_back(par[i]);
+//      cout << par[i] << " ";
+    }
+//    cout << endl;
+
+    if(!only_phase_changed) {
+
+      cout << " Parameters have changed, (re-)calculating p(B) and P(t) now..." << endl;
+
+      vector<double> par_for_BofZ;
+      vector<double> par_for_PofB;
+
+//      cout << "par_for_BofZ: ";
+
+      for (unsigned int i(4); i<par.size(); i++) {
+        par_for_BofZ.push_back(par[i]);
+//        cout << par[i] << " ";
+      }
+//      cout << endl;
+
+//      cout << "par_for_PofB: ";
+
+      for (unsigned int i(1); i<4; i++) {
+        par_for_PofB.push_back(par[i]);
+//        cout << par[i] << " ";
+      }
+//      cout << endl;
+
+      TLondon1D_1L BofZ1(par_for_BofZ);
+      TPofBCalc PofB1(BofZ1, *fImpProfile, par_for_PofB);
+      fPofT->DoFFT(PofB1);
+
+    } else {
+      cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
+    }
+
+    fPofT->CalcPol(par_for_PofT);
+
+    fCalcNeeded = false;
+  }
+
+  return fPofT->Eval(t);
+
+}
+
+//------------------
+// Constructor of the TLondon1D2L class -- reading available implantation profiles and
+// creates (a pointer to) the TPofTCalc object (with the FFT plan)
+//------------------
+
+TLondon1D2L::TLondon1D2L(const vector<unsigned int> &parNo, const vector<double> &par) : fLastTwoChanged(true) {
+
+  for(unsigned int i(0); i<parNo.size(); i++) {
+    fPar.push_back(par[parNo[i]-1]);
+  }
+
+//  string rge_path("/home/l_wojek/nt/wojek/g/Bastian/ImplantationDepth/YBCO_PBCO-");
+//  string energy_arr[] = {"02_1", "02_5", "03_5", "05_0", "07_5", "10_0", "12_5", "15_0", "17_5", "19_0", "20_0", "22_5", "25_0"};
+
+  string rge_path("/home/l_wojek/TrimSP/AuYBCO_2005/AuYBCO-500000-");
+  string energy_arr[] = {"04_6", "09_6", "14_6", "18_6", "21_6", "24_6", "28_7"};
+
+  vector<string> energy_vec(energy_arr, energy_arr+(sizeof(energy_arr)/sizeof(energy_arr[0])));
+
+  vector<double> par_for_PofT;
+
+  for (unsigned int i(0); i<3; i++)
+    par_for_PofT.push_back(fPar[i]);
+
+  TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
+  fImpProfile = x;
+  x = 0;
+  delete x;
+
+  TPofTCalc *y = new TPofTCalc(par_for_PofT);
+  fPofT = y;
+  y = 0;
+  delete y;
+
+}
+
+//------------------
+// TLondon1D2L-Method that calls the procedures to create B(z), p(B) and P(t)
+// It finally returns P(t) for a given t.
+// Parameters: all the parameters for the function to be fitted through TLondon1D1L
+//------------------
+
+double TLondon1D2L::Eval(double t, const vector<double> &par) const {
+
+  // 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 == 1 || i == 2 || i == 3) {
+        cout << "You are varying dt, dB or E! These parameters have to be fixed! Quitting..." << endl; 
+        exit(-1);
+      } else if (i == 0) {
+        only_phase_changed = true;
+      } else {
+        only_phase_changed = false;
+      }
+      if (i == fPar.size()-2 || i == fPar.size()-1)
+        fLastTwoChanged = true;
+    }
+  }
+
+  if (par_changed)
+    fCalcNeeded = true;
+
+/* DEBUGGING CODE COMMENTED -- quite a mess... sorry*/
+
+  // if model parameters have changed, recalculate B(z), P(B) and P(t)
+
+  if (fCalcNeeded) {
+
+    vector<double> par_for_PofT;
+
+//    cout << "par_for_PofT: ";
+
+    for (unsigned int i(0); i<3; i++) {
+      par_for_PofT.push_back(par[i]);
+//      cout << par[i] << " ";
+    }
+//    cout << endl;
+
+    if(!only_phase_changed) {
+
+      cout << " Parameters have changed, (re-)calculating p(B) and P(t) now..." << endl;
+
+      vector<double> par_for_BofZ;
+      vector<double> par_for_PofB;
+
+//      cout << "par_for_BofZ: ";
+
+      for (unsigned int i(4); i<par.size()-2; i++) {
+        par_for_BofZ.push_back(par[i]);
+//        cout << par[i] << " ";
+      }
+//      cout << endl;
+
+//      cout << "par_for_PofB: ";
+
+      for (unsigned int i(1); i<4; i++) {
+        par_for_PofB.push_back(par[i]);
+//        cout << par[i] << " ";
+      }
+//      cout << endl;
+
+      if(fLastTwoChanged) {
+        vector<double> interfaces;
+        interfaces.push_back(par[5]+par[6]);
+
+        vector<double> weights;
+        for(unsigned int i(par.size()-2); i<par.size(); i++)
+          weights.push_back(par[i]);
+
+        cout << "Weighting has changed, re-calculating n(z) now..." << endl;
+        fImpProfile->WeightLayers(par[3], interfaces, weights);
+      }
+
+      TLondon1D_2L BofZ2(par_for_BofZ);
+      TPofBCalc PofB2(BofZ2, *fImpProfile, par_for_PofB);
+      fPofT->DoFFT(PofB2);
+
+    } else {
+      cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
+    }
+
+    fPofT->CalcPol(par_for_PofT);
+
+    fCalcNeeded = false;
+    fLastTwoChanged = false;
+  }
+
+  return fPofT->Eval(t);
+
+}
+
+//------------------
+// Constructor of the TLondon1D3LS class -- reading available implantation profiles and
+// creates (a pointer to) the TPofTCalc object (with the FFT plan)
+//------------------
+
+TLondon1D3LS::TLondon1D3LS(const vector<unsigned int> &parNo, const vector<double> &par) : fLastThreeChanged(true) {
+
+  for(unsigned int i(0); i<parNo.size(); i++) {
+    fPar.push_back(par[parNo[i]-1]);
+  }
+
+//  string rge_path("/home/l_wojek/nt/wojek/g/Bastian/ImplantationDepth/YBCO_PBCO-");
+//  string energy_arr[] = {"02_1", "02_5", "03_5", "05_0", "07_5", "10_0", "12_5", "15_0", "17_5", "19_0", "20_0", "22_5", "25_0"};
+
+  string rge_path("/home/l_wojek/TrimSP/AuYBCO_2005/AuYBCO-500000-");
+  string energy_arr[] = {"04_6", "09_6", "14_6", "18_6", "21_6", "24_6", "28_7"};
+
+  vector<string> energy_vec(energy_arr, energy_arr+(sizeof(energy_arr)/sizeof(energy_arr[0])));
+
+  vector<double> par_for_PofT;
+
+  for (unsigned int i(0); i<3; i++)
+    par_for_PofT.push_back(fPar[i]);
+
+  TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
+  fImpProfile = x;
+  x = 0;
+  delete x;
+
+  TPofTCalc *y = new TPofTCalc(par_for_PofT);
+  fPofT = y;
+  y = 0;
+  delete y;
+
+}
+
+//------------------
+// TLondon1D3LS-Method that calls the procedures to create B(z), p(B) and P(t)
+// It finally returns P(t) for a given t.
+// Parameters: all the parameters for the function to be fitted through TLondon1D1L
+//------------------
+
+double TLondon1D3LS::Eval(double t, const vector<double> &par) const {
+
+  // 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 == 1 || i == 2 || i == 3) {
+        cout << "You are varying dt, dB or E! These parameters have to be fixed! Quitting..." << endl; 
+        exit(-1);
+      } else if (i == 0) {
+        only_phase_changed = true;
+      } else {
+        only_phase_changed = false;
+      }
+      if (i == fPar.size()-3 || i == fPar.size()-2 || i == fPar.size()-1)
+        fLastThreeChanged = true;
+    }
+  }
+
+  if (par_changed)
+    fCalcNeeded = true;
+
+/* DEBUGGING CODE COMMENTED -- quite a mess... sorry*/
+
+  // if model parameters have changed, recalculate B(z), P(B) and P(t)
+
+  if (fCalcNeeded) {
+
+    vector<double> par_for_PofT;
+
+//    cout << "par_for_PofT: ";
+
+    for (unsigned int i(0); i<3; i++) {
+      par_for_PofT.push_back(par[i]);
+//      cout << par[i] << " ";
+    }
+//    cout << endl;
+
+    if(!only_phase_changed) {
+
+      cout << " Parameters have changed, (re-)calculating p(B) and P(t) now..." << endl;
+
+      vector<double> par_for_BofZ;
+      vector<double> par_for_PofB;
+
+//      cout << "par_for_BofZ: ";
+
+      for (unsigned int i(4); i<par.size()-2; i++) {
+        par_for_BofZ.push_back(par[i]);
+//        cout << par[i] << " ";
+      }
+//      cout << endl;
+
+//      cout << "par_for_PofB: ";
+
+      for (unsigned int i(1); i<4; i++) {
+        par_for_PofB.push_back(par[i]);
+//        cout << par[i] << " ";
+      }
+//      cout << endl;
+
+      if(fLastThreeChanged) {
+        vector<double> interfaces;
+        interfaces.push_back(par[5]+par[6]);
+        interfaces.push_back(par[5]+par[6]+par[7]);
+
+        vector<double> weights;
+        for(unsigned int i(par.size()-3); i<par.size(); i++)
+          weights.push_back(par[i]);
+
+        cout << "Weighting has changed, re-calculating n(z) now..." << endl;
+        fImpProfile->WeightLayers(par[3], interfaces, weights);
+      }
+
+      TLondon1D_3L BofZ3(par_for_BofZ);
+      TPofBCalc PofB3(BofZ3, *fImpProfile, par_for_PofB);
+      fPofT->DoFFT(PofB3);
+
+    } else {
+      cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
+    }
+
+    fPofT->CalcPol(par_for_PofT);
+
+    fCalcNeeded = false;
+    fLastThreeChanged = false;
+  }
+
+  return fPofT->Eval(t);
+
+}
diff --git a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp
index de2aac13..3f2702bc 100644
--- a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp
+++ b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp
@@ -5,7 +5,7 @@
   Author: Bastian M. Wojek
   e-mail: bastian.wojek@psi.ch
 
-  2008/05/25
+  2008/05/27
 
 ***************************************************************************/
 
@@ -30,7 +30,7 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
   fBmax = BofZ.GetBmax();
 
   double BB, BBnext;
-  double zm, zp, zNextm, zNextp, dz;
+  double zm, zp, zNext, dz;
 
   // fill not used Bs before Bmin with 0.0
 
@@ -54,7 +54,7 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
   seconds = time (NULL);
 
   char debugfile[50];
-  int n = sprintf (debugfile, "test_Bz_%f_%ld.dat", fBmin, seconds);
+  int n = sprintf (debugfile, "test_Bz_%ld_%f.dat", seconds, fBmin);
 
   if (n > 0) {
     ofstream of(debugfile);
@@ -66,7 +66,7 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
   }
 
   char debugfile1[50];
-  int n1 = sprintf (debugfile1, "test_NZ_%f_%ld.dat", para[2], seconds);
+  int n1 = sprintf (debugfile1, "test_NZ_%ld_%f.dat", seconds, para[2]);
 
   if (n1 > 0) {
     ofstream of1(debugfile1);
@@ -80,6 +80,8 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
 ---------------------------------------------------------*/
 
   double nn;
+  bool zNextFound(false);
+
   for ( ; BB <= fBmax ; BB += para[1]) {
     BBnext = BB + para[1];
     fB.push_back(BB);
@@ -92,16 +94,21 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
         for (unsigned int k(0); k < j; k++) {
           if ( ( bofzBZ[j-k] <= BBnext && bofzBZ[j-k-1] >= BBnext ) ) {
 //           cout << "1 " << j << " " << k << endl;
-            zNextm = (BBnext-bofzBZ[j-k-1])*ddZ/(bofzBZ[j-k]-bofzBZ[j-k-1]) + bofzZ[j-k-1];
+            zNext = (BBnext-bofzBZ[j-k-1])*ddZ/(bofzBZ[j-k]-bofzBZ[j-k-1]) + bofzZ[j-k-1];
+            zNextFound = true;
             break;
           }
         }
 
-        dz = zNextm-zm;
-        nn = dataTrimSP.GetNofZ(zm, para[2]);
-        if (nn != -1.0) {
+        if(zNextFound) {
+          zNextFound = false;
+
+          dz = zNext-zm;
+          nn = dataTrimSP.GetNofZ(zm, para[2]);
+          if (nn != -1.0) {
 //          cout << "zNext = " << zNextm << ", zm = " << zm << ", dz = " << dz << endl;
-          *(fPB.end()-1) += nn*fabs(dz/para[1]);
+            *(fPB.end()-1) += nn*fabs(dz/para[1]);
+          }
         }
 
       } else if (bofzBZ[j] <= BB && bofzBZ[j+1] >= BB ) {
@@ -110,16 +117,21 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
         for (unsigned int k(0); k < bofzZ.size() - j - 1; k++) {
           if ( ( bofzBZ[j+k] <= BBnext && bofzBZ[j+k+1] >= BBnext ) ) {
 //            cout << "2 " << j << " " << k << endl;
-            zNextp = (BBnext-bofzBZ[j+k])*ddZ/(bofzBZ[j+k+1]-bofzBZ[j+k]) + bofzZ[j+k];
+            zNext = (BBnext-bofzBZ[j+k])*ddZ/(bofzBZ[j+k+1]-bofzBZ[j+k]) + bofzZ[j+k];
+            zNextFound = true;
             break;
           }
         }
 
-        dz = zNextp-zp;
-        nn = dataTrimSP.GetNofZ(zp, para[2]);
-        if (nn != -1.0) {
-//          cout << "zNext = " << zNextp << ", zp = " << zp << ", dz = " << dz << endl;
-          *(fPB.end()-1) += nn*fabs(dz/para[1]);
+        if(zNextFound) {
+          zNextFound = false;
+
+          dz = zNext-zp;
+          nn = dataTrimSP.GetNofZ(zp, para[2]);
+          if (nn != -1.0) {
+//            cout << "zNext = " << zNextp << ", zp = " << zp << ", dz = " << dz << endl;
+            *(fPB.end()-1) += nn*fabs(dz/para[1]);
+          }
         }
      }
     }
diff --git a/src/external/TFitPofB-lib/classes/TUserLondon.cpp b/src/external/TFitPofB-lib/classes/TUserLondon.cpp
deleted file mode 100644
index db6ca04a..00000000
--- a/src/external/TFitPofB-lib/classes/TUserLondon.cpp
+++ /dev/null
@@ -1,209 +0,0 @@
-/***************************************************************************
-
-  TUserLondon.cpp
-
-  Author: Bastian M. Wojek
-  e-mail: bastian.wojek@psi.ch
-
-  2008/05/26
-
-***************************************************************************/
-
-#include "TUserLondon.h"
-#include <iostream>
-using namespace std;
-
-//------------------
-// Constructor of the TUserLondon class -- reading available implantation profiles and
-// creates (a pointer to) the TPofTCalc object (with the FFT plan)
-//------------------
-
-TUserLondon::TUserLondon(const vector<unsigned int> &parNo, const vector<double> &par) 
-  : fCalcNeeded(true), fLastTwoChanged(true), fLastThreeChanged(true) {
-
-  for(unsigned int i(0); i<parNo.size(); i++) {
-    fPar.push_back(par[parNo[i]-1]);
-  }
-
-//  string rge_path("/home/l_wojek/nt/wojek/g/Bastian/ImplantationDepth/YBCO_PBCO-");
-//  string energy_arr[] = {"02_1", "02_5", "03_5", "05_0", "07_5", "10_0", "12_5", "15_0", "17_5", "19_0", "20_0", "22_5", "25_0"};
-
-  string rge_path("/home/l_wojek/TrimSP/AuYBCO_2005/AuYBCO-500000-");
-  string energy_arr[] = {"04_6", "09_6", "14_6", "18_6", "21_6", "24_6", "28_7"};
-
-  vector<string> energy_vec(energy_arr, energy_arr+(sizeof(energy_arr)/sizeof(energy_arr[0])));
-
-  vector<double> par_for_PofT;
-
-  for (unsigned int i(1); i<4; i++)
-    par_for_PofT.push_back(fPar[i]);
-
-  TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
-  fImpProfile = x;
-  x = 0;
-  delete x;
-
-  TPofTCalc *y = new TPofTCalc(par_for_PofT);
-  fPofT = y;
-  y = 0;
-  delete y;
-
-}
-
-//------------------
-// Destructor of the TUserLondon class -- cleaning up
-//------------------
-
-TUserLondon::~TUserLondon() {
-    fPar.clear();
-    delete fImpProfile;
-    fImpProfile = 0;
-    delete fPofT;
-    fPofT = 0;
-}
-
-//------------------
-// Method that calls the procedures to create B(z), p(B) and P(t)
-// It finally returns P(t) for a given t.
-// Parameters: all the parameters for the function to be fitted through TUserLondon
-//------------------
-
-double TUserLondon::Eval(double t, const vector<double> &par) const {
-
-  // 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 || i == 2 || i == 3 || i == 4) {
-        cout << "You are varying the model-flag, dt, dB or E! These parameters have to be fixed! Quitting..." << endl; 
-        exit(-1);
-      } else if (i == 1) {
-        only_phase_changed = true;
-      } else {
-        only_phase_changed = false;
-      }
-      if(i == fPar.size()-3) {
-        fLastThreeChanged = true; 
-      } else if (i == fPar.size()-2 || i == fPar.size()-1) {
-        fLastTwoChanged = true;
-        fLastThreeChanged = true;
-      }
-    }
-  }
-
-  if (par_changed)
-    fCalcNeeded = true;
-
-
-/* DEBUGGING CODE COMMENTED -- quite a mess... sorry*/
-
-  // if model parameters have changed, recalculate B(z), P(B) and P(t)
-
-  if (fCalcNeeded) {
-
-    vector<double> par_for_PofT;
-
-//    cout << "par_for_PofT: ";
-
-    for (unsigned int i(1); i<4; i++) {
-      par_for_PofT.push_back(par[i]);
-//      cout << par[i] << " ";
-    }
-//    cout << endl;
-
-    if(!only_phase_changed) {
-
-      cout << " Parameters have changed, (re-)calculating p(B) and P(t) now..." << endl;
-
-      vector<double> par_for_BofZ;
-      vector<double> par_for_PofB;
-
-//      cout << "par_for_BofZ: ";
-
-      for (unsigned int i(5); i<par.size(); i++) {
-        par_for_BofZ.push_back(par[i]);
-//        cout << par[i] << " ";
-      }
-//      cout << endl;
-
-//      cout << "par_for_PofB: ";
-
-      for (unsigned int i(2); i<5; i++) {
-        par_for_PofB.push_back(par[i]);
-//        cout << par[i] << " ";
-      }
-//      cout << endl;
-
-      switch (int(par[0])) {
-        case 1:
-        {
-//          cout << "Found the 1D-London model." << endl;
-          TLondon1D_1L BofZ1(par_for_BofZ);
-          TPofBCalc PofB1(BofZ1, *fImpProfile, par_for_PofB);
-          fPofT->DoFFT(PofB1);
-        }
-          break;
-        case 2:
-        {
-//          cout << "Found the 1D-London model.2L" << endl;
-          if(fLastTwoChanged) {
-            vector<double> interfaces;
-            interfaces.push_back(par[6]+par[7]);
-
-            vector<double> weights;
-            for(unsigned int i(11); i<par.size(); i++)
-              weights.push_back(par[i]);
-
-            cout << "Weighting has changed, re-calculating n(z) now..." << endl;
-            fImpProfile->WeightLayers(par[4], interfaces, weights);
-          }
-          TLondon1D_2L BofZ2(par_for_BofZ);
-          TPofBCalc PofB2(BofZ2, *fImpProfile, par_for_PofB);
-          fPofT->DoFFT(PofB2);
-        }
-          break;
-        case 3:
-        {
-//          cout << "Found the 1D-London model.3L" << endl;
-          if(fLastThreeChanged) {
-            vector<double> interfaces;
-            interfaces.push_back(par[6]+par[7]);
-            interfaces.push_back(par[6]+par[7]+par[8]);
-
-            vector<double> weights;
-            for(unsigned int i(12); i<par.size(); i++)
-              weights.push_back(par[i]);
-
-            cout << "Weighting has changed, re-calculating n(z) now..." << endl;
-            fImpProfile->WeightLayers(par[4], interfaces, weights);
-          }
-          TLondon1D_3L BofZ3(par_for_BofZ);
-          TPofBCalc PofB3(BofZ3, *fImpProfile, par_for_PofB);
-          fPofT->DoFFT(PofB3);
-        }
-          break;
-        default:
-          cout << "The user function you specified with the first parameter, does not exist. Quitting..." << endl;
-          exit(-1);
-      }
-    } else {
-      cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
-    }
-
-    fPofT->CalcPol(par_for_PofT);
-
-    fCalcNeeded = false;
-    fLastTwoChanged = false;
-    fLastThreeChanged = false;
-
-  }
-
-  return fPofT->Eval(t);
-
-}
-
diff --git a/src/external/TFitPofB-lib/include/TLondon1D.h b/src/external/TFitPofB-lib/include/TLondon1D.h
new file mode 100644
index 00000000..3cf2709d
--- /dev/null
+++ b/src/external/TFitPofB-lib/include/TLondon1D.h
@@ -0,0 +1,68 @@
+/***************************************************************************
+
+  TLondon1D.h
+
+  Author: Bastian M. Wojek
+  e-mail: bastian.wojek@psi.ch
+
+  2008/05/27
+
+***************************************************************************/
+
+#ifndef _TLondon1D_H_
+#define _TLondon1D_H_
+
+#include "TPofTCalc.h"
+
+class TLondon1D {
+
+public:
+  TLondon1D() : fCalcNeeded(true) {}
+  virtual ~TLondon1D();
+
+  virtual double Eval(double, const vector<double>&) const = 0;
+
+protected:
+  mutable vector<double> fPar;
+  TTrimSPData *fImpProfile;
+  TPofTCalc *fPofT;
+  mutable bool fCalcNeeded;
+};
+
+class TLondon1D1L : public TLondon1D {
+
+public:
+  TLondon1D1L(const vector<unsigned int>& , const vector<double>&);
+  ~TLondon1D1L() {}
+
+  double Eval(double, const vector<double>&) const;
+
+};
+
+class TLondon1D2L : public TLondon1D {
+
+public:
+  TLondon1D2L(const vector<unsigned int>& , const vector<double>&);
+  ~TLondon1D2L() {}
+
+  double Eval(double, const vector<double>&) const;
+
+private:
+  mutable bool fLastTwoChanged;
+
+};
+
+class TLondon1D3LS : public TLondon1D {
+
+public:
+  TLondon1D3LS(const vector<unsigned int>& , const vector<double>&);
+  ~TLondon1D3LS() {}
+
+  double Eval(double, const vector<double>&) const;
+
+private:
+  mutable bool fLastThreeChanged;
+
+};
+
+#endif //_TUserLondon1D_H_
diff --git a/src/external/TFitPofB-lib/include/TUserLondon.h b/src/external/TFitPofB-lib/include/TUserLondon.h
deleted file mode 100644
index 1a73ca03..00000000
--- a/src/external/TFitPofB-lib/include/TUserLondon.h
+++ /dev/null
@@ -1,35 +0,0 @@
-/***************************************************************************
-
-  TUserLondon.h
-
-  Author: Bastian M. Wojek
-  e-mail: bastian.wojek@psi.ch
-
-  2008/05/25
-
-***************************************************************************/
-
-#ifndef _TUserLondon_H_
-#define _TUserLondon_H_
-
-#include "TPofTCalc.h"
-
-class TUserLondon {
-
-public:
-  TUserLondon(const vector<unsigned int>& , const vector<double>&);
-  ~TUserLondon();
-
-  double Eval(double, const vector<double>&) const;
-
-private:
-  mutable vector<double> fPar;
-  TTrimSPData *fImpProfile;
-  TPofTCalc *fPofT;
-  mutable bool fCalcNeeded;
-  mutable bool fLastTwoChanged;
-  mutable bool fLastThreeChanged;
-
-};
-
-#endif //_TUserLondon_H_
diff --git a/src/external/TFitPofB-lib/test/test.cpp b/src/external/TFitPofB-lib/test/test.cpp
index 3590058c..79bba95b 100644
--- a/src/external/TFitPofB-lib/test/test.cpp
+++ b/src/external/TFitPofB-lib/test/test.cpp
@@ -1,4 +1,4 @@
-#include "TUserLondon.h"
+#include "TLondon1D.h"
 #include <iostream>
 #include <fstream>
 
@@ -95,14 +95,13 @@ int main(){
     of8 << i << " " << poft(i) << endl;
   }
   of8.close();
-  
-  
-  
-  
-  
+
 */
-  unsigned int parNo_arr[] = {1, 3, 5, 7, 9, 11, 12, 13, 14, 15, 16, 17, 18};
-  double par_arr[] = {2.0, 999.0, 0.0, 999.0, 0.01, 999.0, 0.01, 999.0, 21.6, 999.0, 100.0, 5.0, 70.0, 75.0, 180.0, 500.0, 1.0, 0.3};
+
+/**************** Test TLondon1D1L *********************************
+
+  unsigned int parNo_arr[] = {1, 3, 5, 7, 9, 10, 11, 12};
+  double par_arr[] = {0.0, 999.0, 0.01, 999.0, 0.01, 999.0, 21.6, 999.0, 100.0, 5.0, 190.0, 180.0};
   
   vector<unsigned int> parNo_vec(parNo_arr, parNo_arr+(sizeof(parNo_arr)/sizeof(parNo_arr[0])));
   vector<double> par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0])));
@@ -113,7 +112,45 @@ int main(){
     par_vec_sub.push_back(par_vec[parNo_vec[i]-1]);
   }
   
-  TUserLondon fitter(parNo_vec, par_vec);
+  TLondon1D1L fitter(parNo_vec, par_vec);
+
+************************************************************************/
+
+/**************** Test TLondon1D2L **********************************
+
+  unsigned int parNo_arr[] = {1, 3, 5, 7, 9, 10, 11, 12, 13, 14, 15, 16};
+  double par_arr[] = {0.0, 999.0, 0.01, 999.0, 0.01, 999.0, 21.6, 999.0, 100.0, 5.0, 70.0, 70.0, 180.0, 500.0, 1.0, 0.3};
+  
+  vector<unsigned int> parNo_vec(parNo_arr, parNo_arr+(sizeof(parNo_arr)/sizeof(parNo_arr[0])));
+  vector<double> par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0])));
+  
+  vector<double> par_vec_sub;
+  
+  for(unsigned int i(0); i<parNo_vec.size(); i++) {
+    par_vec_sub.push_back(par_vec[parNo_vec[i]-1]);
+  }
+  
+  TLondon1D2L fitter(parNo_vec, par_vec);
+
+************************************************************************/
+
+/**************** Test TLondon1D3LS *********************************/
+
+  unsigned int parNo_arr[] = {1, 3, 5, 7, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18};
+  double par_arr[] = {0.0, 999.0, 0.01, 999.0, 0.01, 999.0, 21.6, 999.0, 100.0, 5.0, 70.0, 50.0, 70.0, 180.0, 500.0, 1.0, 0.3, 1.0};
+  
+  vector<unsigned int> parNo_vec(parNo_arr, parNo_arr+(sizeof(parNo_arr)/sizeof(parNo_arr[0])));
+  vector<double> par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0])));
+  
+  vector<double> par_vec_sub;
+  
+  for(unsigned int i(0); i<parNo_vec.size(); i++) {
+    par_vec_sub.push_back(par_vec[parNo_vec[i]-1]);
+  }
+  
+  TLondon1D3LS fitter(parNo_vec, par_vec);
+
+/************************************************************************/
 
   ofstream of01("test_fitter01.dat");
   ofstream of02("test_fitter02.dat");
@@ -132,60 +169,60 @@ int main(){
   }
   of01.close();
   
-  par_vec_sub[1] += 10.0;
-  par_vec_sub[10] -= 20.0;
+  par_vec_sub[0] += 10.0;
+  par_vec_sub[10] -= 10.0;
   
   for (double i(0.); i<12.0; i+=0.003) {
     of02 << i << " " << fitter.Eval(i, par_vec_sub) << endl;
   }
   of02.close();
   
-    par_vec_sub[1] += 10.0;
-    par_vec_sub[10] -= 20.0;
+    par_vec_sub[0] += 10.0;
+    par_vec_sub[10] -= 10.0;
   
   for (double i(0.); i<12.0; i+=0.003) {
     of03 << i << " " << fitter.Eval(i, par_vec_sub) << endl;
   }
   of03.close();
   
-    par_vec_sub[1] += 10.0;
+    par_vec_sub[0] += 10.0;
   
   for (double i(0.); i<12.0; i+=0.003) {
     of04 << i << " " << fitter.Eval(i, par_vec_sub) << endl;
   }
   of04.close();
   
-    par_vec_sub[1] += 10.0;
-    par_vec_sub[10] -= 20.0;
+    par_vec_sub[0] += 10.0;
+    par_vec_sub[10] -= 10.0;
   
   for (double i(0.); i<12.0; i+=0.003) {
     of05 << i << " " << fitter.Eval(i, par_vec_sub) << endl;
   }
   of05.close();
   
-    par_vec_sub[1] += 10.0;
-    par_vec_sub[10] -= 20.0;
+    par_vec_sub[0] += 10.0;
+    par_vec_sub[10] -= 10.0;
   
   for (double i(0.); i<12.0; i+=0.003) {
     of06 << i << " " << fitter.Eval(i, par_vec_sub) << endl;
   }
   of06.close();
   
-    par_vec_sub[1] += 10.0;
+    par_vec_sub[0] += 10.0;
   
   for (double i(0.); i<12.0; i+=0.003) {
     of07 << i << " " << fitter.Eval(i, par_vec_sub) << endl;
   }
   of07.close();
   
-    par_vec_sub[1] += 10.0;
+    par_vec_sub[0] += 10.0;
   
   for (double i(0.); i<12.0; i+=0.003) {
     of08 << i << " " << fitter.Eval(i, par_vec_sub) << endl;
   }
   of08.close();
   
-    par_vec_sub[1] = 0.0;
+    par_vec_sub[0] = 0.0;
     par_vec_sub[10] = 1000.0;
   
   for (double i(0.); i<12.0; i+=0.003) {
@@ -193,7 +230,7 @@ int main(){
   }
   of09.close();
   
-    par_vec_sub[10] = 500.0;
+    par_vec_sub[0] = 0.0;
     par_vec_sub[12] = 1.0;
   
   for (double i(0.); i<12.0; i+=0.003) {
@@ -201,63 +238,63 @@ int main(){
   }
   of10.close();
 
-  
-/*    vector<double> data01, data02, data03, data04, data05, data06, data07, data08, data09, data10;
+/*
+    vector<double> data01, data02, data03, data04, data05, data06, data07, data08, data09, data10;
     
     for (double i(0.); i<12.0; i+=0.003)
     data01.push_back(fitter.Eval(i, par_vec_sub));
 
   
-  par_vec_sub[1] += 10.0;
-  par_vec_sub[8] -= 10.0;
+  par_vec_sub[0] += 10.0;
+  par_vec_sub[10] -= 10.0;
   
   for (double i(0.); i<12.0; i+=0.003)
     data02.push_back(fitter.Eval(i, par_vec_sub));
 
   
-    par_vec_sub[1] += 10.0;
-    par_vec_sub[8] -= 10.0;
+    par_vec_sub[0] += 10.0;
+    par_vec_sub[10] -= 10.0;
   
   for (double i(0.); i<12.0; i+=0.003)
   data03.push_back(fitter.Eval(i, par_vec_sub));
   
-    par_vec_sub[1] += 10.0;
+    par_vec_sub[0] += 10.0;
   
   for (double i(0.); i<12.0; i+=0.003)
    data04.push_back(fitter.Eval(i, par_vec_sub));
    
 
   
-    par_vec_sub[1] += 10.0;
-    par_vec_sub[8] -= 10.0;
+    par_vec_sub[0] += 10.0;
+    par_vec_sub[10] -= 10.0;
   
   for (double i(0.); i<12.0; i+=0.003)
     data05.push_back(fitter.Eval(i, par_vec_sub));
   
-    par_vec_sub[1] += 10.0;
-    par_vec_sub[8] -= 10.0;
+    par_vec_sub[0] += 10.0;
+    par_vec_sub[10] -= 10.0;
   
   for (double i(0.); i<12.0; i+=0.003) 
     data06.push_back(fitter.Eval(i, par_vec_sub));
   
-    par_vec_sub[1] += 10.0;
+    par_vec_sub[0] += 10.0;
   
   for (double i(0.); i<12.0; i+=0.003) 
     data07.push_back(fitter.Eval(i, par_vec_sub));
   
-    par_vec_sub[1] += 10.0;
-    par_vec_sub[7] = 190.0;
+    par_vec_sub[0] += 10.0;
+    par_vec_sub[10] = 190.0;
   
   for (double i(0.); i<12.0; i+=0.003) 
     data08.push_back(fitter.Eval(i, par_vec_sub));
   
-    par_vec_sub[1] += 10.0;
-    par_vec_sub[8] -= 10.0;
+    par_vec_sub[0] += 10.0;
+    par_vec_sub[10] -= 10.0;
   
   for (double i(0.); i<12.0; i+=0.003) 
     data09.push_back(fitter.Eval(i, par_vec_sub));
   
-    par_vec_sub[1] += 10.0;
+    par_vec_sub[0] += 10.0;
   
   for (double i(0.); i<12.0; i+=0.003) 
     data10.push_back(fitter.Eval(i, par_vec_sub));