From 3726cc68e806d93e39de680204e47303a56c5e38 Mon Sep 17 00:00:00 2001
From: nemu <nemu>
Date: Sat, 24 Jan 2009 11:12:37 +0000
Subject: [PATCH] renaming train

---
 src/tests/dynKT_LF/dynGaussKT_LF.cpp | 250 +++++++++++++++++++++++++++
 src/tests/dynKT_LF/dynGaussKT_LF.pro |  20 +++
 src/tests/dynKT_LF/nr/lubksb.cpp     |  24 +++
 src/tests/dynKT_LF/nr/ludcmp.cpp     |  53 ++++++
 src/tests/dynKT_LF/nr/voltra.cpp     |  36 ++++
 src/tests/dynKT_LF/nr/xvoltra.cpp    |  44 +++++
 6 files changed, 427 insertions(+)
 create mode 100644 src/tests/dynKT_LF/dynGaussKT_LF.cpp
 create mode 100644 src/tests/dynKT_LF/dynGaussKT_LF.pro
 create mode 100644 src/tests/dynKT_LF/nr/lubksb.cpp
 create mode 100644 src/tests/dynKT_LF/nr/ludcmp.cpp
 create mode 100644 src/tests/dynKT_LF/nr/voltra.cpp
 create mode 100644 src/tests/dynKT_LF/nr/xvoltra.cpp

diff --git a/src/tests/dynKT_LF/dynGaussKT_LF.cpp b/src/tests/dynKT_LF/dynGaussKT_LF.cpp
new file mode 100644
index 00000000..91fb8075
--- /dev/null
+++ b/src/tests/dynKT_LF/dynGaussKT_LF.cpp
@@ -0,0 +1,250 @@
+#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/dynKT_LF/dynGaussKT_LF.pro b/src/tests/dynKT_LF/dynGaussKT_LF.pro
new file mode 100644
index 00000000..b4e57901
--- /dev/null
+++ b/src/tests/dynKT_LF/dynGaussKT_LF.pro
@@ -0,0 +1,20 @@
+#------------------------------------------------------
+# 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/dynKT_LF/nr/lubksb.cpp b/src/tests/dynKT_LF/nr/lubksb.cpp
new file mode 100644
index 00000000..c16c4914
--- /dev/null
+++ b/src/tests/dynKT_LF/nr/lubksb.cpp
@@ -0,0 +1,24 @@
+#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/dynKT_LF/nr/ludcmp.cpp b/src/tests/dynKT_LF/nr/ludcmp.cpp
new file mode 100644
index 00000000..127c57b2
--- /dev/null
+++ b/src/tests/dynKT_LF/nr/ludcmp.cpp
@@ -0,0 +1,53 @@
+#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/dynKT_LF/nr/voltra.cpp b/src/tests/dynKT_LF/nr/voltra.cpp
new file mode 100644
index 00000000..c593a4fa
--- /dev/null
+++ b/src/tests/dynKT_LF/nr/voltra.cpp
@@ -0,0 +1,36 @@
+#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/dynKT_LF/nr/xvoltra.cpp b/src/tests/dynKT_LF/nr/xvoltra.cpp
new file mode 100644
index 00000000..887aedae
--- /dev/null
+++ b/src/tests/dynKT_LF/nr/xvoltra.cpp
@@ -0,0 +1,44 @@
+#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;
+}