renaming and adding Keren SG

This commit is contained in:
nemu 2009-01-24 13:46:22 +00:00
parent 3726cc68e8
commit f06ab1b700
8 changed files with 41 additions and 459 deletions

View File

@ -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

View File

@ -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];
}
}

View File

@ -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;
}
}
}

View File

@ -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];
}
}

View File

@ -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;
}

View File

@ -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;
}

View File

@ -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;
}

View File

@ -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