diff --git a/src/external/libLFRelaxation/Makefile.libLFRelaxation b/src/external/libLFRelaxation/Makefile.libLFRelaxation
new file mode 100644
index 00000000..35d53b1f
--- /dev/null
+++ b/src/external/libLFRelaxation/Makefile.libLFRelaxation
@@ -0,0 +1,64 @@
+#---------------------------------------------------
+# get compilation flags from root-config
+
+ROOTCFLAGS    = $(shell $(ROOTSYS)/bin/root-config --cflags)
+# ROOTLIBS    = $(shell $(ROOTSYS)/bin/root-config --libs)
+
+#---------------------------------------------------
+
+OS            = LINUX
+CXX           = g++
+CXXFLAGS      = -g -Wall -Wno-trigraphs -fPIC
+MUSRFITINCLUDE = ../../include
+#MUSRFITINCLUDE = /home/l_wojek/rep/analysis/musrfit/src/include
+LOCALINCLUDE  = .
+ROOTINCLUDE   = $(ROOTSYS)/include/root
+INCLUDES      = -I$(LOCALINCLUDE) -I$(MUSRFITINCLUDE) -I$(ROOTINCLUDE)
+LD            = g++
+LDFLAGS       = -g
+SOFLAGS       = -O -shared
+
+# the output from the root-config script:
+CXXFLAGS      += $(ROOTCFLAGS)
+LDFLAGS       +=
+
+# some definitions: headers (used to generate *Dict* stuff), sources, objects,...
+OBJS =
+OBJS += TIntegrator.o
+OBJS += TLFRelaxation.o TLFRelaxationDict.o
+
+SHLIB = libLFRelaxation.so
+
+# make the shared lib:
+#
+all:  $(SHLIB)
+
+$(SHLIB): $(OBJS)
+	@echo "---> Building shared library $(SHLIB) ..."
+	/bin/rm -f $(SHLIB)
+	$(LD) $(LDFLAGS) $(OBJS) $(SOFLAGS) -o $(SHLIB)
+	@echo "done"
+
+# clean up: remove all object file (and core files)
+# semicolon needed to tell make there is no source
+# for this target!
+#
+clean:; 	@rm -f $(OBJS) $(SHLIB) *Dict* core* *~
+	@echo "---> removing $(OBJS) $(SHLIB)"
+
+#
+$(OBJS): %.o: %.cpp
+	$(CXX) $(INCLUDES) $(CXXFLAGS) -c $<
+
+# Generate the ROOT CINT dictionary
+
+TLFRelaxationDict.cpp: ./TLFRelaxation.h ./TLFRelaxationLinkDef.h
+	@echo "Generating dictionary $@..."
+	rootcint -f $@ -c -p -I$(MUSRFITINCLUDE) $^
+
+install: all
+	@echo "Installing shared lib: libLFRelaxation.so"
+ifeq ($(OS),LINUX)
+	cp -pv $(SHLIB) $(ROOTSYS)/lib
+	cp -pv $(LOCALINCLUDE)/*.h $(ROOTSYS)/include
+endif
diff --git a/src/external/libLFRelaxation/README b/src/external/libLFRelaxation/README
new file mode 100644
index 00000000..143deb04
--- /dev/null
+++ b/src/external/libLFRelaxation/README
@@ -0,0 +1,19 @@
+/***************************************************************************
+
+  Author: Bastian M. Wojek
+  e-mail: bastian.wojek@psi.ch
+
+  2008/12/05
+
+***************************************************************************/
+
+Implementation of a userFcn-interface to Gaussian and Lorentzian static and dynamic LF relaxation functions.
+At the moment this is localized to l_wojek@pc5405, because an absolute path had to be set. Of course this can be easily 
+changed in the code if needed.
+
+The functions are then called from within musrfit as:
+
+userFcn  libLFRelaxation.so TLFStatGssKT 1 2       (frequency rate)
+userFcn  libLFRelaxation.so TLFStatLorKT 1 2       (frequency rate)
+userFcn  libLFRelaxation.so TLFDynGssKT 1 2 3      (frequency rate fluct.rate)
+userFcn  libLFRelaxation.so TLFDynLorKT 1 2 3      (frequency rate fluct.rate)
diff --git a/src/external/libLFRelaxation/TIntegrator.cpp b/src/external/libLFRelaxation/TIntegrator.cpp
new file mode 100644
index 00000000..3a21e2e9
--- /dev/null
+++ b/src/external/libLFRelaxation/TIntegrator.cpp
@@ -0,0 +1,52 @@
+/***************************************************************************
+
+  TIntegrator.cpp
+
+  Author: Bastian M. Wojek
+  e-mail: bastian.wojek@psi.ch
+
+  2008/12/03
+
+***************************************************************************/
+
+#include "TIntegrator.h"
+#include "TMath.h"
+
+using namespace std;
+
+TIntegrator::TIntegrator() : fFunc(0) {
+  ROOT::Math::GSLIntegrator *integrator = new ROOT::Math::GSLIntegrator(ROOT::Math::Integration::kADAPTIVE,ROOT::Math::Integration::kGAUSS51);
+  fIntegrator = integrator;
+  integrator = 0;
+  delete integrator;
+}
+
+TIntegrator::~TIntegrator(){
+  delete fIntegrator;
+  fIntegrator=0;
+  fFunc=0;
+}
+
+inline double TIntegrator::FuncAtXgsl(double x, void *obj)
+{
+  return ((TIntegrator*)obj)->FuncAtX(x);
+}
+
+double TIntegrator::IntegrateFunc(double x1, double x2)
+{
+  fFunc = &TIntegrator::FuncAtXgsl;
+  return fIntegrator->Integral(fFunc, (this), x1, x2);
+}
+
+
+inline double TIntBesselJ0Exp::FuncAtX(double x) const
+{
+  return TMath::BesselJ0(TMath::TwoPi()*fPar[0]*x) * TMath::Exp(-fPar[1]*x);
+}
+
+
+inline double TIntSinGss::FuncAtX(double x) const
+{
+  return TMath::Sin(TMath::TwoPi()*fPar[0]*x) * TMath::Exp(-0.5*fPar[1]*fPar[1]*x*x);
+}
+
diff --git a/src/external/libLFRelaxation/TIntegrator.h b/src/external/libLFRelaxation/TIntegrator.h
new file mode 100644
index 00000000..108526c4
--- /dev/null
+++ b/src/external/libLFRelaxation/TIntegrator.h
@@ -0,0 +1,52 @@
+/***************************************************************************
+
+  TIntegrator.h
+
+  Author: Bastian M. Wojek
+  e-mail: bastian.wojek@psi.ch
+
+  2008/12/03
+
+***************************************************************************/
+
+#ifndef _TIntegrator_H_
+#define _TIntegrator_H_
+
+#include "Math/GSLIntegrator.h"
+
+#include<vector>
+
+using namespace std;
+
+class TIntegrator {
+  public:
+    TIntegrator();
+    virtual ~TIntegrator();
+    void SetParameters(const std::vector<double> &par) const { fPar=par; }
+    virtual double FuncAtX(double) const = 0;
+    double IntegrateFunc(double, double);
+
+  protected:
+    mutable vector<double> fPar;
+
+  private:
+    static double FuncAtXgsl(double, void *);
+    ROOT::Math::GSLIntegrator *fIntegrator;
+    mutable double (*fFunc)(double, void *);
+};
+
+class TIntBesselJ0Exp : public TIntegrator {
+  public:
+    TIntBesselJ0Exp() {}
+    ~TIntBesselJ0Exp() {}
+    double FuncAtX(double) const;
+};
+
+class TIntSinGss : public TIntegrator {
+  public:
+    TIntSinGss() {}
+    ~TIntSinGss() {}
+    double FuncAtX(double) const;
+};
+
+#endif //_TIntegrator_H_
diff --git a/src/external/libLFRelaxation/TLFRelaxation.cpp b/src/external/libLFRelaxation/TLFRelaxation.cpp
new file mode 100644
index 00000000..2c02b49c
--- /dev/null
+++ b/src/external/libLFRelaxation/TLFRelaxation.cpp
@@ -0,0 +1,375 @@
+/***************************************************************************
+
+  TLFRelaxation.cpp
+
+  Author: Bastian M. Wojek
+  e-mail: bastian.wojek@psi.ch
+
+  2008/12/04
+
+***************************************************************************/
+
+#include <cassert>
+
+#include "TIntegrator.h"
+#include "TLFRelaxation.h"
+
+using namespace std;
+
+ClassImp(TLFStatGssKT)
+ClassImp(TLFStatLorKT)
+ClassImp(TLFDynGssKT)
+ClassImp(TLFDynLorKT)
+
+// LF Static Gaussian KT
+
+TLFStatGssKT::TLFStatGssKT() {
+  TIntSinGss *intSinGss = new TIntSinGss();
+  fIntSinGss = intSinGss;
+  intSinGss = 0;
+  delete intSinGss;
+}
+
+TLFStatGssKT::~TLFStatGssKT() {
+    delete fIntSinGss;
+    fIntSinGss = 0;
+}
+
+double TLFStatGssKT::operator()(double t, const vector<double> &par) const {
+
+  assert(par.size() == 2); // two parameters nu=gbar*B,sigma
+
+  if(t<0.0)
+    return 1.0;
+
+  double sigsq(par[1]*par[1]);
+
+  if(TMath::Abs(par[0])<0.00135538817){
+    return (0.33333333333333333333+0.66666666666666666667*(1.0-sigsq*t*t)*TMath::Exp(-0.5*sigsq*t*t));
+  }
+
+  fIntSinGss->SetParameters(par);
+
+  double omegaL(TMath::TwoPi()*par[0]);
+  double coeff1(2.0*sigsq/(omegaL*omegaL));
+  double coeff2(2.0*sigsq*sigsq/(omegaL*omegaL*omegaL));
+
+  return (1.0-(coeff1*(1.0-TMath::Exp(-0.5*sigsq*t*t)*TMath::Cos(omegaL*t)))+(coeff2*fIntSinGss->IntegrateFunc(0.,t)));
+}
+
+// LF Static Lorentzian KT
+
+TLFStatLorKT::TLFStatLorKT() {
+  TIntBesselJ0Exp *intBesselJ0Exp = new TIntBesselJ0Exp();
+  fIntBesselJ0Exp = intBesselJ0Exp;
+  intBesselJ0Exp = 0;
+  delete intBesselJ0Exp;
+}
+
+TLFStatLorKT::~TLFStatLorKT() {
+    delete fIntBesselJ0Exp;
+    fIntBesselJ0Exp = 0;
+}
+
+double TLFStatLorKT::operator()(double t, const vector<double> &par) const {
+
+  assert(par.size() == 2); // two parameters nu=gbar*B,rate
+
+  if(t<0.0)
+    return 1.0;
+
+  if(TMath::Abs(par[0])<0.00135538817){
+    return (0.33333333333333333333+0.66666666666666666667*(1.0-par[1]*t)*TMath::Exp(-par[1]*t));
+  }
+
+  fIntBesselJ0Exp->SetParameters(par);
+
+  double omegaL(TMath::TwoPi()*par[0]);
+  double coeff1(par[1]/omegaL);
+  double coeff2(coeff1*coeff1);
+  double coeff3((1.0+coeff2)*par[1]);
+
+  return (1.0-(coeff1*TMath::Exp(-par[1]*t)*TMath::BesselJ1(omegaL*t))-(coeff2*(TMath::BesselJ0(omegaL*t)*TMath::Exp(-par[1]*t)-1.0))-coeff3*fIntBesselJ0Exp->IntegrateFunc(0.,t));
+}
+
+// LF Dynamic Gaussian KT
+
+TLFDynGssKT::TLFDynGssKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("/home/l_wojek/analysis/WordsOfWisdom.dat"), fNSteps(524288), fDt(0.000040), fCounter(0) {
+  // Calculate d_omega and C for given NFFT and dt
+  fDw = TMath::Pi()/fNSteps/fDt;
+  fC = 2.0*TMath::Log(double(fNSteps))/(double(fNSteps-1)*fDt);
+
+  // Load FFTW Wisdom
+  int wisdomLoaded(0);
+
+  FILE *wordsOfWisdomR;
+  wordsOfWisdomR = fopen(fWisdom.c_str(), "r");
+  if (wordsOfWisdomR == NULL) {
+    cout << "TLFDynGssKT::TLFDynGssKT: Couldn't open wisdom file ..." << endl;
+  } else {
+    wisdomLoaded = fftw_import_wisdom_from_file(wordsOfWisdomR);
+    fclose(wordsOfWisdomR);
+  }
+
+  if (!wisdomLoaded) {
+    cout << "TLFDynGssKT::TLFDynGssKT: No wisdom is imported..." << endl;
+  }
+  // END of WisdomLoading
+
+  // allocating memory for the FFtransform pairs and create the FFT plans
+
+  fFFTtime = (double *)malloc(sizeof(double) * fNSteps);
+  fFFTfreq = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (fNSteps/2+1));
+  fFFTplanFORW = fftw_plan_dft_r2c_1d(fNSteps, fFFTtime, fFFTfreq, FFTW_EXHAUSTIVE);
+  fFFTplanBACK = fftw_plan_dft_c2r_1d(fNSteps, fFFTfreq, fFFTtime, FFTW_EXHAUSTIVE);
+}
+
+TLFDynGssKT::~TLFDynGssKT() {
+  // export FFTW Wisdom so it has not to be checked for the FFT-plan next time
+  FILE *wordsOfWisdomW;
+  wordsOfWisdomW = fopen(fWisdom.c_str(), "w");
+  if (wordsOfWisdomW == NULL) {
+    cout << "TLFDynGssKT::~TLFDynGssKT: Could not open file ... No wisdom is exported..." << endl;
+  } else {
+    fftw_export_wisdom_to_file(wordsOfWisdomW);
+    fclose(wordsOfWisdomW);
+  }
+  // END of Wisdom Export
+
+  // clean up
+  fftw_destroy_plan(fFFTplanFORW);
+  fftw_destroy_plan(fFFTplanBACK);
+  free(fFFTtime);
+  fftw_free(fFFTfreq);
+  cout << "TLFDynGssKT::~TLFDynGssKT(): " << fCounter << " full FFT cycles needed..." << endl;
+}
+
+double TLFDynGssKT::operator()(double t, const vector<double> &par) const {
+
+  assert(par.size() == 3); // three parameters nuL=gbar*B,sigma,fluct.rate nu
+
+  if(t<0.0)
+    return 1.0;
+
+  if(t>20.0)
+    return 0.0;
+
+  if(fFirstCall){
+    fPar = par;
+    fFirstCall=false;
+  }
+
+  for (unsigned int i(0); i<par.size(); i++) {
+    if( fPar[i]-par[i] ) {
+      fPar[i] = par[i];
+      fCalcNeeded=true;
+    }
+  }
+
+  double sigsq(par[1]*par[1]); // sigma^2
+  double omegaL(TMath::TwoPi()*par[0]); // Larmor frequency
+  double nusq(par[2]*par[2]); // nu^2
+
+  if(par[1]){
+    if(par[2]/par[1] > 5. || omegaL > 20.0*par[1]){
+      if(TMath::Abs(par[0])<0.00135538817){
+        return TMath::Exp(-2.0*sigsq/nusq*(TMath::Exp(-par[2]*t)-1.0+par[2]*t)); // ZF Abragam Delta^2->2*Delta^2
+      }
+      double omegaLnusqp(omegaL*omegaL+nusq);
+      double omegaLnusqm(omegaL*omegaL-nusq);
+      return TMath::Exp(-2.0*sigsq/(omegaLnusqp*omegaLnusqp)*(omegaLnusqp*par[2]*t+omegaLnusqm*(1.0-TMath::Exp(-par[2]*t)*TMath::Cos(omegaL*t))-2.0*par[2]*omegaL*TMath::Exp(-par[2]*t)*TMath::Sin(omegaL*t))); // Keren
+    }
+  }
+
+  if(fCalcNeeded){
+
+    double tt(0.);
+
+    if(TMath::Abs(par[0])<0.00135538817){
+      for(unsigned int i(0); i<fNSteps; i++) {
+        tt=double(i)*fDt;
+        fFFTtime[i]=(0.33333333333333333333+0.66666666666666666667*(1.0-sigsq*tt*tt)*TMath::Exp(-0.5*sigsq*tt*tt))*TMath::Exp(-(fC+par[2])*tt)*fDt;
+      }
+    } else {
+
+      double coeff1(2.0*sigsq/(omegaL*omegaL));
+      double coeff2(2.0*sigsq*sigsq/(omegaL*omegaL*omegaL));
+
+      fFFTtime[0] = 1.0*fDt;
+
+      for(unsigned int i(1); i<fNSteps; i++) {
+        tt=(double(i)-0.5)*fDt;
+        fFFTtime[i]=(TMath::Sin(omegaL*tt) * TMath::Exp(-0.5*sigsq*tt*tt))*fDt;
+      }
+
+      double totoIntegrale(0.);
+
+      for(unsigned int i(1); i<fNSteps; i++) {
+        tt=double(i)*fDt;
+        totoIntegrale+=fFFTtime[i];
+        fFFTtime[i]=(1.0-(coeff1*(1.0-TMath::Exp(-0.5*sigsq*tt*tt)*TMath::Cos(omegaL*tt)))+(coeff2*totoIntegrale))*TMath::Exp(-(fC+par[2])*tt)*fDt;
+      }
+    }
+
+    // Transform to frequency domain
+
+    fftw_execute(fFFTplanFORW);
+
+    // calculate F(s)
+
+    double denom(1.0);
+
+    for (unsigned int i(0); i<fNSteps/2+1; i++) {
+      denom=(1.0-(par[2]*fFFTfreq[i][0]))*(1.0-(par[2]*fFFTfreq[i][0])) + (nusq*fFFTfreq[i][1]*fFFTfreq[i][1]);
+      fFFTfreq[i][0] = (fFFTfreq[i][0]-(par[2]*fFFTfreq[i][0]*fFFTfreq[i][0])-(par[2]*fFFTfreq[i][1]*fFFTfreq[i][1]))/denom;
+      fFFTfreq[i][1] = fFFTfreq[i][1]/denom;
+    }
+
+    // Transform back to time domain
+
+    fftw_execute(fFFTplanBACK);
+
+//    for (unsigned int i(0); i<fNSteps; i++) {
+//      fFFTtime[i]=(fDw*TMath::Exp(fC*i*fDt)/TMath::Pi()*fFFTtime[i]);
+//    }
+    fCalcNeeded=false;
+    fCounter++;
+  }
+//  return fFFTtime[int(t/fDt)];
+  return fDw*TMath::Exp(fC*t)/TMath::Pi()*fFFTtime[int(t/fDt)];
+}
+
+// LF Dynamic Lorentz KT
+
+TLFDynLorKT::TLFDynLorKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("/home/l_wojek/analysis/WordsOfWisdom.dat"), fNSteps(524288), fDt(0.000040), fCounter(0) {
+  // Calculate d_omega and C for given NFFT and dt
+  fDw = TMath::Pi()/fNSteps/fDt;
+  fC = 2.0*TMath::Log(double(fNSteps))/(double(fNSteps-1)*fDt);
+
+  // Load FFTW Wisdom
+  int wisdomLoaded(0);
+
+  FILE *wordsOfWisdomR;
+  wordsOfWisdomR = fopen(fWisdom.c_str(), "r");
+  if (wordsOfWisdomR == NULL) {
+    cout << "TLFDynGssKT::TLFDynGssKT: Couldn't open wisdom file ..." << endl;
+  } else {
+    wisdomLoaded = fftw_import_wisdom_from_file(wordsOfWisdomR);
+    fclose(wordsOfWisdomR);
+  }
+
+  if (!wisdomLoaded) {
+    cout << "TLFDynGssKT::TLFDynGssKT: No wisdom is imported..." << endl;
+  }
+  // END of WisdomLoading
+
+  // allocating memory for the FFtransform pairs and create the FFT plans
+
+  fFFTtime = (double *)malloc(sizeof(double) * fNSteps);
+  fFFTfreq = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (fNSteps/2+1));
+  fFFTplanFORW = fftw_plan_dft_r2c_1d(fNSteps, fFFTtime, fFFTfreq, FFTW_EXHAUSTIVE);
+  fFFTplanBACK = fftw_plan_dft_c2r_1d(fNSteps, fFFTfreq, fFFTtime, FFTW_EXHAUSTIVE);
+}
+
+TLFDynLorKT::~TLFDynLorKT() {
+  // export FFTW Wisdom so it has not to be checked for the FFT-plan next time
+  FILE *wordsOfWisdomW;
+  wordsOfWisdomW = fopen(fWisdom.c_str(), "w");
+  if (wordsOfWisdomW == NULL) {
+    cout << "TLFDynGssKT::~TLFDynGssKT: Could not open file ... No wisdom is exported..." << endl;
+  } else {
+    fftw_export_wisdom_to_file(wordsOfWisdomW);
+    fclose(wordsOfWisdomW);
+  }
+  // END of Wisdom Export
+
+  // clean up
+  fftw_destroy_plan(fFFTplanFORW);
+  fftw_destroy_plan(fFFTplanBACK);
+  free(fFFTtime);
+  fftw_free(fFFTfreq);
+  cout << "TLFDynLorKT::~TLFDynLorKT(): " << fCounter << " full FFT cyles needed..." << endl;
+}
+
+double TLFDynLorKT::operator()(double t, const vector<double> &par) const {
+
+  assert(par.size() == 3); // three parameters nuL=gbar*B,sigma,fluct.rate nu
+
+  if(t<0.0)
+    return 1.0;
+
+  if(t>20.0)
+    return 0.0;
+
+  if(fFirstCall){
+    fPar = par;
+    fFirstCall=false;
+  }
+
+  for (unsigned int i(0); i<par.size(); i++) {
+    if( fPar[i]-par[i] ) {
+      fPar[i] = par[i];
+      fCalcNeeded=true;
+    }
+  }
+
+  if(fCalcNeeded){
+
+    double tt(0.);
+
+    if(TMath::Abs(par[0])<0.00135538817){
+      for(unsigned int i(0); i<fNSteps; i++) {
+        tt=double(i)*fDt;
+        fFFTtime[i]=(0.33333333333333333333+0.66666666666666666667*(1.0-par[1]*tt)*TMath::Exp(-par[1]*tt))*TMath::Exp(-(fC+par[2])*tt)*fDt;
+      }
+    } else {
+
+      double omegaL(TMath::TwoPi()*par[0]); // Larmor frequency
+      double coeff1(par[1]/omegaL);
+      double coeff2(coeff1*coeff1);
+      double coeff3((1.0+coeff2)*par[1]);
+
+      fFFTtime[0] = 1.0*fDt;
+
+      for(unsigned int i(1); i<fNSteps; i++) {
+        tt=(double(i)-0.5)*fDt;
+        fFFTtime[i]=TMath::BesselJ0(omegaL*tt)*TMath::Exp(-par[1]*tt)*fDt;
+      }
+
+      double totoIntegrale(0.);
+
+      for(unsigned int i(1); i<fNSteps; i++) {
+        tt=double(i)*fDt;
+        totoIntegrale+=fFFTtime[i];
+        fFFTtime[i]=(1.0-(coeff1*TMath::Exp(-par[1]*tt)*TMath::BesselJ1(omegaL*tt))-(coeff2*(TMath::BesselJ0(omegaL*tt)*TMath::Exp(-par[1]*tt)-1.0))-coeff3*totoIntegrale)*TMath::Exp(-(fC+par[2])*tt)*fDt;
+      }
+    }
+
+    // Transform to frequency domain
+
+    fftw_execute(fFFTplanFORW);
+
+    // calculate F(s)
+
+    double denom(1.0);
+    double nusq(par[2]*par[2]); // nu^2
+
+    for (unsigned int i(0); i<fNSteps/2+1; i++) {
+      denom=(1.0-(par[2]*fFFTfreq[i][0]))*(1.0-(par[2]*fFFTfreq[i][0])) + (nusq*fFFTfreq[i][1]*fFFTfreq[i][1]);
+      fFFTfreq[i][0] = (fFFTfreq[i][0]-(par[2]*fFFTfreq[i][0]*fFFTfreq[i][0])-(par[2]*fFFTfreq[i][1]*fFFTfreq[i][1]))/denom;
+      fFFTfreq[i][1] = fFFTfreq[i][1]/denom;
+    }
+
+    // Transform back to time domain
+
+    fftw_execute(fFFTplanBACK);
+
+//    for (unsigned int i(0); i<fNSteps; i++) {
+//      fFFTtime[i]=(fDw*TMath::Exp(fC*i*fDt)/TMath::Pi()*fFFTtime[i]);
+//    }
+    fCalcNeeded=false;
+    fCounter++;
+  }
+//  return fFFTtime[int(t/fDt)];
+  return fDw*TMath::Exp(fC*t)/TMath::Pi()*fFFTtime[int(t/fDt)];
+}
diff --git a/src/external/libLFRelaxation/TLFRelaxation.h b/src/external/libLFRelaxation/TLFRelaxation.h
new file mode 100644
index 00000000..b8964316
--- /dev/null
+++ b/src/external/libLFRelaxation/TLFRelaxation.h
@@ -0,0 +1,105 @@
+/***************************************************************************
+
+  TLFRelaxation.h
+
+  Author: Bastian M. Wojek
+  e-mail: bastian.wojek@psi.ch
+
+  2008/12/04
+
+***************************************************************************/
+
+#ifndef _TLFRelaxation_H_
+#define _TLFRelaxation_H_
+
+#include<vector>
+#include<cstdio>
+
+using namespace std;
+
+#include "TMath.h"
+#include "PUserFcnBase.h"
+#include "fftw3.h"
+#include "TIntegrator.h"
+
+class TLFStatGssKT : public PUserFcnBase {
+
+public:
+  TLFStatGssKT();
+  ~TLFStatGssKT();
+
+  double operator()(double, const vector<double>&) const;
+
+private:
+  TIntSinGss *fIntSinGss;
+
+  ClassDef(TLFStatGssKT,1)
+};
+
+class TLFStatLorKT : public PUserFcnBase {
+
+public:
+  TLFStatLorKT();
+  ~TLFStatLorKT();
+
+  double operator()(double, const vector<double>&) const;
+
+private:
+  TIntBesselJ0Exp *fIntBesselJ0Exp;
+
+  ClassDef(TLFStatLorKT,1)
+};
+
+class TLFDynGssKT : public PUserFcnBase {
+
+public:
+  TLFDynGssKT();
+  ~TLFDynGssKT();
+
+  double operator()(double, const vector<double>&) const;
+
+private:
+  mutable vector<double> fPar;
+  mutable bool fCalcNeeded;
+  mutable bool fFirstCall;
+  string fWisdom;
+  unsigned int fNSteps;
+  double fDt;
+  double fDw;
+  double fC;
+  fftw_plan fFFTplanFORW;
+  fftw_plan fFFTplanBACK;
+  double *fFFTtime;
+  fftw_complex *fFFTfreq;
+  mutable unsigned int fCounter;
+
+  ClassDef(TLFDynGssKT,1)
+};
+
+class TLFDynLorKT : public PUserFcnBase {
+
+public:
+  TLFDynLorKT();
+  ~TLFDynLorKT();
+
+  double operator()(double, const vector<double>&) const;
+
+private:
+  mutable vector<double> fPar;
+  mutable bool fCalcNeeded;
+  mutable bool fFirstCall;
+  string fWisdom;
+  unsigned int fNSteps;
+  double fDt;
+  double fDw;
+  double fC;
+  fftw_plan fFFTplanFORW;
+  fftw_plan fFFTplanBACK;
+  double *fFFTtime;
+  fftw_complex *fFFTfreq;
+  mutable unsigned int fCounter;
+
+  ClassDef(TLFDynLorKT,1)
+};
+
+#endif //_LFRelaxation_H_
diff --git a/src/external/libLFRelaxation/TLFRelaxationLinkDef.h b/src/external/libLFRelaxation/TLFRelaxationLinkDef.h
new file mode 100644
index 00000000..0c21ff24
--- /dev/null
+++ b/src/external/libLFRelaxation/TLFRelaxationLinkDef.h
@@ -0,0 +1,26 @@
+/***************************************************************************
+
+  TLFRelaxationLinkDef.h
+
+  Author: Bastian M. Wojek
+  e-mail: bastian.wojek@psi.ch
+
+  2008/12/04
+
+***************************************************************************/
+
+// root dictionary stuff --------------------------------------------------
+#ifdef __CINT__
+
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ class TLFStatGssKT+;
+#pragma link C++ class TLFStatLorKT+;
+#pragma link C++ class TLFDynGssKT+;
+#pragma link C++ class TLFDynLorKT+;
+
+#endif //__CINT__
+// root dictionary stuff --------------------------------------------------
+