diff --git a/src/tests/dynGaussKT_LF/dynGaussKT_LF.pro b/src/tests/dynGaussKT_LF/dynGaussKT_LF.pro
deleted file mode 100644
index b4e57901..00000000
--- a/src/tests/dynGaussKT_LF/dynGaussKT_LF.pro
+++ /dev/null
@@ -1,20 +0,0 @@
-#------------------------------------------------------
-# dynGaussKT_LF.pro
-# qmake file for dynGaussKT_LF
-#
-# Andreas Suter, 2009/01/13
-#
-# $Id$
-#
-#------------------------------------------------------
-
-MAKEFILE = Makefile
-
-CONFIG -= qt
-CONFIG += warn_on console debug
-
-SOURCES = dynGaussKT_LF.cpp
-
-TARGET=dynGaussKT_LF
-
-unix:LIBS += -lgsl -lgslcblas -lm
diff --git a/src/tests/dynGaussKT_LF/nr/lubksb.cpp b/src/tests/dynGaussKT_LF/nr/lubksb.cpp
deleted file mode 100644
index c16c4914..00000000
--- a/src/tests/dynGaussKT_LF/nr/lubksb.cpp
+++ /dev/null
@@ -1,24 +0,0 @@
-#include "nr.h"
-
-void NR::lubksb(Mat_I_DP &a, Vec_I_INT &indx, Vec_IO_DP &b)
-{
-	int i,ii=0,ip,j;
-	DP sum;
-
-	int n=a.nrows();
-	for (i=0;i<n;i++) {
-		ip=indx[i];
-		sum=b[ip];
-		b[ip]=b[i];
-		if (ii != 0)
-			for (j=ii-1;j<i;j++) sum -= a[i][j]*b[j];
-		else if (sum != 0.0)
-			ii=i+1;
-		b[i]=sum;
-	}
-	for (i=n-1;i>=0;i--) {
-		sum=b[i];
-		for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
-		b[i]=sum/a[i][i];
-	}
-}
diff --git a/src/tests/dynGaussKT_LF/nr/ludcmp.cpp b/src/tests/dynGaussKT_LF/nr/ludcmp.cpp
deleted file mode 100644
index 127c57b2..00000000
--- a/src/tests/dynGaussKT_LF/nr/ludcmp.cpp
+++ /dev/null
@@ -1,53 +0,0 @@
-#include <cmath>
-#include "nr.h"
-using namespace std;
-
-void NR::ludcmp(Mat_IO_DP &a, Vec_O_INT &indx, DP &d)
-{
-	const DP TINY=1.0e-20;
-	int i,imax,j,k;
-	DP big,dum,sum,temp;
-
-	int n=a.nrows();
-	Vec_DP vv(n);
-	d=1.0;
-	for (i=0;i<n;i++) {
-		big=0.0;
-		for (j=0;j<n;j++)
-			if ((temp=fabs(a[i][j])) > big) big=temp;
-		if (big == 0.0) nrerror("Singular matrix in routine ludcmp");
-		vv[i]=1.0/big;
-	}
-	for (j=0;j<n;j++) {
-		for (i=0;i<j;i++) {
-			sum=a[i][j];
-			for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
-			a[i][j]=sum;
-		}
-		big=0.0;
-		for (i=j;i<n;i++) {
-			sum=a[i][j];
-			for (k=0;k<j;k++) sum -= a[i][k]*a[k][j];
-			a[i][j]=sum;
-			if ((dum=vv[i]*fabs(sum)) >= big) {
-				big=dum;
-				imax=i;
-			}
-		}
-		if (j != imax) {
-			for (k=0;k<n;k++) {
-				dum=a[imax][k];
-				a[imax][k]=a[j][k];
-				a[j][k]=dum;
-			}
-			d = -d;
-			vv[imax]=vv[j];
-		}
-		indx[j]=imax;
-		if (a[j][j] == 0.0) a[j][j]=TINY;
-		if (j != n-1) {
-			dum=1.0/(a[j][j]);
-			for (i=j+1;i<n;i++) a[i][j] *= dum;
-		}
-	}
-}
diff --git a/src/tests/dynGaussKT_LF/nr/voltra.cpp b/src/tests/dynGaussKT_LF/nr/voltra.cpp
deleted file mode 100644
index c593a4fa..00000000
--- a/src/tests/dynGaussKT_LF/nr/voltra.cpp
+++ /dev/null
@@ -1,36 +0,0 @@
-#include "nr.h"
-
-void NR::voltra(const DP t0, const DP h, Vec_O_DP &t, Mat_O_DP &f,
-	DP g(const int, const DP),
-	DP ak(const int, const int, const DP, const DP))
-{
-	int i,j,k,l;
-	DP d,sum;
-
-	int m=f.nrows();
-	int n=f.ncols();
-	Vec_INT indx(m);
-	Vec_DP b(m);
-	Mat_DP a(m,m);
-	t[0]=t0;
-	for (k=0;k<m;k++) f[k][0]=g(k,t[0]);
-	for (i=1;i<n;i++) {
-		t[i]=t[i-1]+h;
-		for (k=0;k<m;k++) {
-			sum=g(k,t[i]);
-			for (l=0;l<m;l++) {
-				sum += 0.5*h*ak(k,l,t[i],t[0])*f[l][0];
-				for (j=1;j<i;j++)
-					sum += h*ak(k,l,t[i],t[j])*f[l][j];
-				if (k == l)
-					a[k][l]=1.0-0.5*h*ak(k,l,t[i],t[i]);
-				else
-					a[k][l] = -0.5*h*ak(k,l,t[i],t[i]);
-			}
-			b[k]=sum;
-		}
-		ludcmp(a,indx,d);
-		lubksb(a,indx,b);
-		for (k=0;k<m;k++) f[k][i]=b[k];
-	}
-}
diff --git a/src/tests/dynGaussKT_LF/nr/xvoltra.cpp b/src/tests/dynGaussKT_LF/nr/xvoltra.cpp
deleted file mode 100644
index 887aedae..00000000
--- a/src/tests/dynGaussKT_LF/nr/xvoltra.cpp
+++ /dev/null
@@ -1,44 +0,0 @@
-#include <iostream>
-#include <iomanip>
-#include <cmath>
-#include "nr.h"
-using namespace std;
-
-// Driver for routine voltra
-
-DP g(const int k, const DP t)
-{
-        return (k == 0 ? cosh(t)+t*sin(t) :
-          2.0*sin(t)+t*(SQR(sin(t))+exp(t)));
-}
-
-DP ak(const int k, const int l, const DP t, const DP s)
-{
-        return ((k == 0) ?
-          (l == 0 ? -exp(t-s) : -cos(t-s)) :
-          (l == 0 ? -exp(t+s) : -t*cos(s)));
-}
-
-int main(void)
-{
-        const int N=10,M=2;
-        const DP H=0.05;
-        int nn;
-        DP t0=0.0;
-        Vec_DP t(N);
-        Mat_DP f(M,N);
-
-        NR::voltra(t0,H,t,f,g,ak);
-        // exact soln is f[1]=exp(-t), f[2]=2sin(t)
-        cout << "    abscissa     voltra       real";
-        cout << "         voltra       real" << endl;
-        cout << "                 answer1      answer1";
-        cout << "      answer2      answer2" << endl << endl;
-        cout << fixed << setprecision(6);
-        for (nn=0;nn<N;nn++) {
-          cout << setw(12) << t[nn] << setw(13) << f[0][nn];
-          cout << setw(13) << exp(-t[nn]) << setw(13) << f[1][nn];
-          cout << setw(13) << 2.0*sin(t[nn]) << endl;
-        }
-        return 0;
-}
diff --git a/src/tests/dynKT_LF/dynGaussKT_LF.cpp b/src/tests/dynKT_LF/dynGaussKT_LF.cpp
deleted file mode 100644
index 91fb8075..00000000
--- a/src/tests/dynKT_LF/dynGaussKT_LF.cpp
+++ /dev/null
@@ -1,250 +0,0 @@
-#include <sys/time.h>
-#include <unistd.h>
-
-#include <iostream>
-#include <iomanip>
-#include <vector>
-#include <cmath>
-using namespace std;
-
-#include <gsl/gsl_sf_bessel.h>
-
-typedef vector<double> PDoubleVector;
-
-#define PI 3.14159265359
-
-// NR p.797
-void voltra(const double h, PDoubleVector &t, PDoubleVector &param, PDoubleVector &f,
-            double g(const int, const PDoubleVector &, const PDoubleVector &, const PDoubleVector &), PDoubleVector &gi)
-{
-  int i,j;
-  double sum;
-
-  int n = f.size();
-  double a;
-  t[0] = 0.0;
-  f[0]=g(0,t,param,gi);
-
-  for (i=1; i<n; i++) {
-    sum = g(i,t,param,gi);
-    sum += 0.5*h*param[2]*g(i,t,param,gi)*f[0];
-    for (j=1; j<i; j++) {
-      sum += h*param[2]*g(i-j,t,param,gi)*f[j];
-    }
-    a = 1.0-0.5*h*param[2]*g(0,t,param,gi);
-
-    f[i]=sum/a;
-  }
-}
-
-double g_gauss(const int i, const PDoubleVector &tvec, const PDoubleVector &param, const PDoubleVector &gi)
-{
-  // param: 0=w0, 1=sigma, 2=nu
-  double result;
-  double t = tvec[i];
-  double Dt2 = pow(param[1]*t, 2.0)/2.0;
-
-  if (param[0] == 0.0) {
-    result = 0.333333333333333333333333 + 0.66666666666666666666 * (1.0-2.0*Dt2) * exp(-Dt2);
-  } else {
-    result = 1.0 - 2.0*param[1]*param[1]/(param[0]*param[0])*(1.0-exp(-Dt2)*cos(param[0]*t)) +
-             2.0*pow(param[1],4.0)/pow(param[0],3.0)*gi[i];
-  }
-
-  result *= exp(-param[2]*t);
-
-  return result;
-}
-
-double g_lorentz(const int i, const PDoubleVector &tvec, const PDoubleVector &param, const PDoubleVector &gi)
-{
-  // param: 0=w0, 1=lambda, 2=nu
-  double result;
-  double t = tvec[i];
-  double at = param[1]*t;
-
-  if (param[0] == 0.0) {
-    result = 0.333333333333333333333333 + 0.66666666666666666666 * (1.0-at) * exp(-at);
-  } else {
-    double awL = param[1]/param[0];
-    double wLt = param[0]*t;
-    double expat = exp(-param[1]*t);
-    result = 1.0 - awL*gsl_sf_bessel_j1(wLt)*expat - awL*awL*(gsl_sf_bessel_j0(wLt)*expat-1.0) -
-             (1.0+awL*awL)*param[1]*gi[i];
-  }
-
-  result *= exp(-param[2]*t);
-
-  return result;
-}
-
-void calc_gi_gauss(const double h, PDoubleVector &t, PDoubleVector &param, PDoubleVector &f)
-{
-  // if w0=0 nothing to be done
-  if (param[0] == 0.0)
-    return;
-
-  double dtHalf = h/2.0;
-  PDoubleVector hh(t.size());
-
-  hh[0] = 0.0;
-  f[0] = 0.0;
-  for (unsigned int i=1; i<t.size(); i++) {
-    hh[i] = exp(-0.5*pow(param[1]*t[i],2.0))*sin(param[0]*t[i]);
-    f[i] = f[i-1] + dtHalf*(hh[i]+hh[i-1]);
-  }
-}
-
-void calc_gi_lorentz(const double h, PDoubleVector &t, PDoubleVector &param, PDoubleVector &f)
-{
-  // if w0=0 nothing to be done
-  if (param[0] == 0.0)
-    return;
-
-  double dtHalf = h/2.0;
-  PDoubleVector hh(t.size());
-
-  hh[0] = 0.0;
-  f[0] = 0.0;
-  for (unsigned int i=1; i<t.size(); i++) {
-    hh[i] = gsl_sf_bessel_j0(param[0]*t[i]) * exp(-param[1]*t[i]);
-    f[i] = f[i-1] + dtHalf*(hh[i]+hh[i-1]);
-  }
-}
-
-int main(int argc, char *argv[])
-{
-
-  bool useKeren = false;
-
-  if (argc < 4) {
-    cout << endl << "usage: dynGaussKT_LF w0 width nu [G|L N]";
-    cout << endl << "          w0:    external field in Mc/s";
-    cout << endl << "          width: static field width in Mc/s";
-    cout << endl << "          nu:    hopping rate in Mc/s";
-    cout << endl << "          G/L:   G=Gaussian field distribution; L=Lorentzain field distribution";
-    cout << endl << "          if G/L not given, G is set as default";
-    cout << endl << "          N:     number of sampling points";
-    cout << endl << "          if N is not given, calc N internally";
-    cout << endl << endl;
-    return 0;
-  }
-
-  PDoubleVector param(3);
-  const double Tmax = 15.0;
-  unsigned int N;
-  bool gaussian = true;
-
-  // feed parameter vector
-  param[0] = atof(argv[1]); // w0
-  param[1] = atof(argv[2]); // width
-  param[2] = atof(argv[3]); // nu
-  if (argc == 6) {
-    N = atoi(argv[5]); // N
-  } else {
-    // w0 criteria, i.e. w0 T = 2 pi, ts = T/16, N = Tmax/ts, if N < 300, N == 300
-    double val = 8.0/PI*Tmax*param[0];
-    if (val < 250)
-      N = 250;
-    else
-      N = static_cast<unsigned int>(val);
-
-    // nu/Delta criteria
-    if (param[1] != 0.0) { // Delta != 0
-      val = param[2]/param[1]; // nu/Delta
-      if (val > 5.0) {
-        useKeren = true;
-        N = 3000;
-      }
-    }
-  }
-
-  if (argc > 4) {
-    if (*argv[4] == 'L') {
-      gaussian = false;
-    }
-  }
-
-  const double H = Tmax/N;
-
-  PDoubleVector t(N);
-  PDoubleVector f(N);
-  PDoubleVector gi(N);
-  PDoubleVector keren(N);
-
-  // get start time
-  struct timeval tv_start, tv_stop;
-  double t1, t2;
-  gettimeofday(&tv_start, 0);
-
-  // feed time vector
-  for (unsigned int i=0; i<N; i++) {
-    t[i] = H*i;
-  }
-
-  if (gaussian) {
-    calc_gi_gauss(H, t, param, gi);
-  } else {
-    calc_gi_lorentz(H, t, param, gi);
-  }
-  gettimeofday(&tv_stop, 0);
-  t1 = (tv_stop.tv_sec - tv_start.tv_sec)*1000.0 + (tv_stop.tv_usec - tv_start.tv_usec)/1000.0;
-
-  gettimeofday(&tv_start, 0);
-  if (gaussian) {
-    voltra(H, t, param, f, g_gauss, gi);
-  } else {
-    voltra(H, t, param, f, g_lorentz, gi);
-  }
-
-  // get stop time
-  gettimeofday(&tv_stop, 0);
-  t2 = (tv_stop.tv_sec - tv_start.tv_sec)*1000.0 + (tv_stop.tv_usec - tv_start.tv_usec)/1000.0;
-
-  // calculate keren LF
-  double w02, nu2, Delta2, Gamma_t;
-  for (unsigned int i=0; i<t.size(); i++) {
-    w02 = pow(param[0], 2.0);
-    Delta2 = pow(param[1], 2.0);
-    nu2 = pow(param[2], 2.0);
-    Gamma_t = 2.0*Delta2/pow(w02+nu2,2.0)*((w02+nu2)*param[2]*t[i]+(w02-nu2)*(1.0-exp(-param[2]*t[i])*cos(param[0]*t[i]))-2.0*param[2]*param[0]*exp(-param[2]*t[i])*sin(param[0]*t[i]));
-    keren[i] = exp(-Gamma_t);
-  }
-
-  if (gaussian) {
-    if (useKeren)
-      cout << "# use Keren = true" << endl;
-    else
-      cout << "# use Keren = false" << endl;
-  }
-
-  cout << "# N = " << N << endl;
-  if (gaussian) {
-    cout << "# Gaussian field distribution" << endl;
-    cout << "# w0 = " << param[0] << ", sigma = " << param[1] << ", nu = " << param[2] << endl;
-  } else {
-    cout << "# Lorentzian field distribution" << endl;
-    cout << "# w0 = " << param[0] << ", lambda = " << param[1] << ", nu = " << param[2] << endl;
-  }
-  cout << "# calculation time: t1 = " << t1 << " (ms), t2 = " << t2 << " (ms)" << endl;
-
-  if (gaussian) {
-    cout <<  "#       time" <<  setw(13) << "Pz_dyn_LF" << setw(13) << "g" << setw(13) << "gi" << setw(13) << "keren" << endl;
-  } else {
-    cout <<  "#       time" <<  setw(13) << "Pz_dyn_LF" << setw(13) << "g" << setw(13) << "gi" << endl;
-  }
-  cout << fixed << setprecision(6);
-  if (gaussian) {
-    for (unsigned int nn=0; nn<N; nn++) {
-      cout << setw(12) << t[nn] << setw(13) << f[nn] << setw(13) << g_gauss(nn,t,param,gi) << setw(13) << gi[nn] << setw(13) << keren[nn];
-      cout << endl;
-     }
-   } else {
-    for (unsigned int nn=0; nn<N; nn++) {
-      cout << setw(12) << t[nn] << setw(13) << f[nn] << setw(13) << g_lorentz(nn,t,param,gi) << setw(13) << gi[nn];
-      cout << endl;
-     }
-   }
-
-  return 0;
-}
diff --git a/src/tests/dynGaussKT_LF/dynGaussKT_LF.cpp b/src/tests/dynKT_LF/dynKT_LF.cpp
similarity index 79%
rename from src/tests/dynGaussKT_LF/dynGaussKT_LF.cpp
rename to src/tests/dynKT_LF/dynKT_LF.cpp
index 91fb8075..8cbed311 100644
--- a/src/tests/dynGaussKT_LF/dynGaussKT_LF.cpp
+++ b/src/tests/dynKT_LF/dynKT_LF.cpp
@@ -5,6 +5,7 @@
 #include <iomanip>
 #include <vector>
 #include <cmath>
+#include <fstream>
 using namespace std;
 
 #include <gsl/gsl_sf_bessel.h>
@@ -165,6 +166,12 @@ int main(int argc, char *argv[])
     }
   }
 
+  char fln[128];
+  if (gaussian)
+    sprintf(fln, "dynKT_LF_w0_%1.1f_width%1.1f_nu%1.1f_N%d_G.dat", param[0], param[1], param[2], N);
+  else 
+    sprintf(fln, "dynKT_LF_w0_%1.1f_width%1.1f_nu%1.1f_N%d_L.dat", param[0], param[1], param[2], N);
+
   const double H = Tmax/N;
 
   PDoubleVector t(N);
@@ -202,49 +209,51 @@ int main(int argc, char *argv[])
   t2 = (tv_stop.tv_sec - tv_start.tv_sec)*1000.0 + (tv_stop.tv_usec - tv_start.tv_usec)/1000.0;
 
   // calculate keren LF
-  double w02, nu2, Delta2, Gamma_t;
+  double w02, nu2, width2, Gamma_t;
   for (unsigned int i=0; i<t.size(); i++) {
     w02 = pow(param[0], 2.0);
-    Delta2 = pow(param[1], 2.0);
+    width2 = pow(param[1], 2.0);
     nu2 = pow(param[2], 2.0);
-    Gamma_t = 2.0*Delta2/pow(w02+nu2,2.0)*((w02+nu2)*param[2]*t[i]+(w02-nu2)*(1.0-exp(-param[2]*t[i])*cos(param[0]*t[i]))-2.0*param[2]*param[0]*exp(-param[2]*t[i])*sin(param[0]*t[i]));
-    keren[i] = exp(-Gamma_t);
+    Gamma_t = 2.0*width2/pow(w02+nu2,2.0)*((w02+nu2)*param[2]*t[i]+(w02-nu2)*(1.0-exp(-param[2]*t[i])*cos(param[0]*t[i]))-2.0*param[2]*param[0]*exp(-param[2]*t[i])*sin(param[0]*t[i]));
+    if (gaussian)
+      keren[i] = exp(-Gamma_t);
+    else
+      keren[i] = exp(-sqrt(2.0*Gamma_t));
+  }
+
+  ofstream fout;
+
+  // open mlog-file
+  fout.open(fln, iostream::out);
+  if (!fout.is_open()) {
+    return -1;
   }
 
   if (gaussian) {
     if (useKeren)
-      cout << "# use Keren = true" << endl;
+      fout << "# use Keren = true" << endl;
     else
-      cout << "# use Keren = false" << endl;
+      fout << "# use Keren = false" << endl;
   }
 
-  cout << "# N = " << N << endl;
+  fout << "# N = " << N << endl;
   if (gaussian) {
-    cout << "# Gaussian field distribution" << endl;
-    cout << "# w0 = " << param[0] << ", sigma = " << param[1] << ", nu = " << param[2] << endl;
+    fout << "# Gaussian field distribution" << endl;
+    fout << "# w0 = " << param[0] << ", sigma = " << param[1] << ", nu = " << param[2] << endl;
   } else {
-    cout << "# Lorentzian field distribution" << endl;
-    cout << "# w0 = " << param[0] << ", lambda = " << param[1] << ", nu = " << param[2] << endl;
+    fout << "# Lorentzian field distribution" << endl;
+    fout << "# w0 = " << param[0] << ", lambda = " << param[1] << ", nu = " << param[2] << endl;
   }
-  cout << "# calculation time: t1 = " << t1 << " (ms), t2 = " << t2 << " (ms)" << endl;
+  fout << "# calculation time: t1 = " << t1 << " (ms), t2 = " << t2 << " (ms)" << endl;
 
-  if (gaussian) {
-    cout <<  "#       time" <<  setw(13) << "Pz_dyn_LF" << setw(13) << "g" << setw(13) << "gi" << setw(13) << "keren" << endl;
-  } else {
-    cout <<  "#       time" <<  setw(13) << "Pz_dyn_LF" << setw(13) << "g" << setw(13) << "gi" << endl;
+  fout <<  "#       time" <<  setw(13) << "Pz_dyn_LF" << setw(13) << "g" << setw(13) << "gi" << setw(13) << "keren" << endl;
+  fout << fixed << setprecision(6);
+  for (unsigned int nn=0; nn<N; nn++) {
+    fout << setw(12) << t[nn] << setw(13) << f[nn] << setw(13) << g_gauss(nn,t,param,gi) << setw(13) << gi[nn] << setw(13) << keren[nn];
+    fout << endl;
   }
-  cout << fixed << setprecision(6);
-  if (gaussian) {
-    for (unsigned int nn=0; nn<N; nn++) {
-      cout << setw(12) << t[nn] << setw(13) << f[nn] << setw(13) << g_gauss(nn,t,param,gi) << setw(13) << gi[nn] << setw(13) << keren[nn];
-      cout << endl;
-     }
-   } else {
-    for (unsigned int nn=0; nn<N; nn++) {
-      cout << setw(12) << t[nn] << setw(13) << f[nn] << setw(13) << g_lorentz(nn,t,param,gi) << setw(13) << gi[nn];
-      cout << endl;
-     }
-   }
+
+  fout.close();
 
   return 0;
 }
diff --git a/src/tests/dynKT_LF/dynGaussKT_LF.pro b/src/tests/dynKT_LF/dynKT_LF.pro
similarity index 71%
rename from src/tests/dynKT_LF/dynGaussKT_LF.pro
rename to src/tests/dynKT_LF/dynKT_LF.pro
index b4e57901..fd1178ca 100644
--- a/src/tests/dynKT_LF/dynGaussKT_LF.pro
+++ b/src/tests/dynKT_LF/dynKT_LF.pro
@@ -1,6 +1,6 @@
 #------------------------------------------------------
-# dynGaussKT_LF.pro
-# qmake file for dynGaussKT_LF
+# dynKT_LF.pro
+# qmake file for dynKT_LF
 #
 # Andreas Suter, 2009/01/13
 #
@@ -13,8 +13,8 @@ MAKEFILE = Makefile
 CONFIG -= qt
 CONFIG += warn_on console debug
 
-SOURCES = dynGaussKT_LF.cpp
+SOURCES = dynKT_LF.cpp
 
-TARGET=dynGaussKT_LF
+TARGET=dynKT_LF
 
 unix:LIBS += -lgsl -lgslcblas -lm