| Rev. | ed56057f569a00172c163f91d7f5eab90ecaee64 |
|---|---|
| Größe | 28,638 Bytes |
| Zeit | 2006-12-10 23:58:13 |
| Autor | iselllo |
| Log Message | initial import |
Fork#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_integration.h>
#include <iostream>
#include <cmath>
#include <vector>
#include <algorithm>
#include <complex>
#include <fstream>
#include <blitz/array.h>
#define REAL(z,i) ((z)[2*(i)])
#define IMAG(z,i) ((z)[2*(i)+1])
using namespace std;
using namespace blitz;
/*
void rhs_eval (double x, Array<double,1> y, Array<double,1>& dydx)
{
double k=4;
dydx(0)=y(1);
dydx(1)=-k*y(0);
}
*/
//This function is a simple implementation of a 4th order Runge-Kutta and it can be used to
// solve a simple problem using a fixed time-step.
// the rhs_eval function has to be provided by the user.
void rk4_fixed (double& x, Array<double,1>& y,
void (*rhs_eval)(double, Array<double,1>, Array<double,1>&),
double h)
{
// Array y assumed to be of extent n, where n is no. of coupled
// equations
int n = y.extent(0);
// Declare local arrays
Array<double,1> k1(n), k2(n), k3(n), k4(n), f(n), dydx(n);
// Zeroth intermediate step
(*rhs_eval) (x, y, dydx);
for (int j = 0; j < n; j++)
{
k1(j) = h * dydx(j);
f(j) = y(j) + k1(j) / 2.;
}
// First intermediate step
(*rhs_eval) (x + h / 2., f, dydx);
for (int j = 0; j < n; j++)
{
k2(j) = h * dydx(j);
f(j) = y(j) + k2(j) / 2.;
}
// Second intermediate step
(*rhs_eval) (x + h / 2., f, dydx);
for (int j = 0; j < n; j++)
{
k3(j) = h * dydx(j);
f(j) = y(j) + k3(j);
}
// Third intermediate step
(*rhs_eval) (x + h, f, dydx);
for (int j = 0; j < n; j++)
{
k4(j) = h * dydx(j);
}
// Actual step
for (int j = 0; j < n; j++)
{
y(j) += k1(j) / 6. + k2(j) / 3. + k3(j) / 3. + k4(j) / 6.;
}
x += h;
return;
}
// This function implements an adaptive step 4-th order Runge-Kutta. It makes use of the previous function
// and again the rhs_eval function has to be provided by the user.
void rk4_adaptive (double& x, Array<double,1>& y,
void (*rhs_eval)(double, Array<double,1>, Array<double,1>&),
double& h, double& t_err, double acc,
double S, int& rept, int maxrept,
double h_min, double h_max, int flag)
{
// Array y assumed to be of extent n, where n is no. of coupled
// equations
int n = y.extent(0);
// Declare local arrays
Array<double,1> y0(n), y1(n);
// Declare repetition counter
static int count = 0;
// Save initial data
double x0 = x;
y0 = y;
// Take full step
rk4_fixed (x, y, rhs_eval, h);
// Save data
y1 = y;
// Restore initial data
x = x0;
y = y0;
// Take two half-steps
rk4_fixed (x, y, rhs_eval, h/2.);
rk4_fixed (x, y, rhs_eval, h/2.);
// Calculate truncation error
t_err = 0.;
double err, err1, err2;
if (flag == 0)
{
// Use absolute truncation error
for (int i = 0; i < n; i++)
{
err = fabs (y(i) - y1(i));
t_err = (err > t_err) ? err : t_err;
}
}
else if (flag == 1)
{
// Use relative truncation error
for (int i = 0; i < n; i++)
{
err = fabs ((y(i) - y1(i)) / y(i));
t_err = (err > t_err) ? err : t_err;
}
}
else
{
// Use mixed truncation error
for (int i = 0; i < n; i++)
{
err1 = fabs ((y(i) - y1(i)) / y(i));
err2 = fabs (y(i) - y1(i));
err = (err1 < err2) ? err1 : err2;
t_err = (err > t_err) ? err : t_err;
}
}
// Prevent small truncation error from rounding to zero
if (t_err == 0.) t_err = 1.e-15;
// Calculate new step-length
double h_est = h * pow (fabs (acc / t_err), 0.2);
// Prevent step-length from changing by more than factor S
if (h_est / h > S)
h *= S;
else if (h_est / h < 1. / S)
h /= S;
else
h = h_est;
// Prevent step-length from exceeding h_max
h = (fabs(h) > h_max) ? h_max * h / fabs(h) : h;
// Abort if step-length falls below h_min
if (fabs(h) < h_min)
{
printf ("Error - |h| < hmin\n");
exit (1);
}
// If truncation error acceptable take step
if ((t_err <= acc) || (count >= maxrept))
{
rept = count;
count = 0;
}
// If truncation error unacceptable repeat step
else
{
count++;
x = x0;
y = y0;
rk4_adaptive (x, y, rhs_eval, h, t_err, acc,
S, rept, maxrept, h_min, h_max, flag);
}
return;
}
// now I start providing the functions I need in order to solve the DNLS in the homogeneous case
// I do not need to change the names of the functions, because their arguments are different and the C++
//compiler is smart enough to understand which function is being called, so the function names can remain
// the same.
//this function provides the RHS for the homogenous DNLS
void rhs_eval (double x, Array<complex<double>,1> y,
Array<complex<double>,1>& dydx,double U, double Jt,
int m,Array<complex<double>,1> state,int periodic)
{
int end=(state.extent(0)-1);
typedef complex<double> dcomp;
dcomp invi=dcomp(0,-1);
if (m>end)
{
cout<<"m is too large"<<endl;
}
if (m<0)
{
cout<<"m is too small"<<endl;
}
if (m==0)
{
if (periodic==0)
{
// cout<<"Periodic is zero"<<endl;
dydx(0)=invi*(Jt*(state(m+1))+
U*abs(y(0))*abs(y(0))*y(0));
}
else if (periodic==1)
{
dydx(0)=invi*(Jt*(state(m+1)+state(end))+
U*abs(y(0))*abs(y(0))*y(0));
}
else
{
cout<<"mistake in chosing the boundary conditions "<<endl;
}
}
// In these two cases I impose periodic boundary conditions
else if (m==end)
{
if (periodic == 0)
{
dydx(0)=invi*(Jt*(state(m-1))+
U*abs(y(0))*abs(y(0))*y(0));
}
else if (periodic ==1)
{
dydx(0)=invi*(Jt*(state(m-1)+state(0))+
U*abs(y(0))*abs(y(0))*y(0));
}
else
{
cout<<"mistake in chosing the boundary conditions for end "<<endl;
}
//cout<<"dydx(0) for m=end is "<<dydx(0)<<endl;
}
else
{
dydx(0)=invi*(Jt*(state(m-1)+state(m+1))+
U*abs(y(0))*abs(y(0))*y(0));
}
}
// this function now solves the homogeneous DNLS in the lattice using the previous function
// providing the RHS.
void rk4_fixed (double& x, Array<complex<double>,1>& y,
void (*rhs_eval)(double, Array<complex<double>,1>,
Array<complex<double>,1>& ,double,
double,int ,Array<complex<double>,1>, int),
double h, int m, Array<complex<double>,1> state,
double U, double Jt,int periodic)
{
// Array y assumed to be of extent n, where n is no. of coupled
// equations
int n = y.extent(0);
int end=(state.extent(0)-1); //final well in the initial array
// Declare local arrays
// Array<double,1> k1(n), k2(n), k3(n), k4(n), f(n), dydx(n);
Array<complex<double>,1> k1(n), k2(n), k3(n), k4(n), f(n), dydx(n);
// Zeroth intermediate step
(*rhs_eval) (x, y, dydx,U,Jt,m,state,periodic);
for (int j = 0; j < n; j++)
{
k1(j) = h * dydx(j);
f(j) = y(j) + k1(j) / 2.;
}
// First intermediate step
(*rhs_eval) (x + h / 2., f, dydx,U,Jt,m,state,periodic);
for (int j = 0; j < n; j++)
{
k2(j) = h * dydx(j);
f(j) = y(j) + k2(j) / 2.;
}
// Second intermediate step
(*rhs_eval) (x + h / 2., f, dydx,U,Jt,m,state,periodic);
for (int j = 0; j < n; j++)
{
k3(j) = h * dydx(j);
f(j) = y(j) + k3(j);
}
// Third intermediate step
(*rhs_eval) (x + h, f, dydx,U,Jt,m,state,periodic);
for (int j = 0; j < n; j++)
{
k4(j) = h * dydx(j);
}
// Actual step
for (int j = 0; j < n; j++)
{
y(j) += k1(j) / 6. + k2(j) / 3. + k3(j) / 3. + k4(j) / 6.;
}
if (m==end)
{
x += h; //I update time only once all the initial state has been evolved
}
return;
}
// I explicitely add the function t work out the RHS of the equation only in the case of the homogenous
// DNLS, where it is non-trivial
// now I add the function to use the adaptive time step
// this function, instead, solves the homogeneous DNLS in the lattice using some adaptive step
// 4th order RK.
void rk4_adaptive (double& x, Array<complex<double>,1>& y,
void (*rhs_eval)(double, Array<complex<double>,1>,
Array<complex<double>,1>& ,double,
double,int ,Array<complex<double>,1>, int),
double& h, double& t_err, double acc,
double S, int& rept, int maxrept,
double h_min, double h_max, int flag,
int m, Array<complex<double>,1> state,
double U, double Jt, int periodic )
{
// Array y assumed to be of extent n, where n is no. of coupled
// equations
int n = y.extent(0);
// Declare local arrays
// Array<double,1> y0(n), y1(n);
Array<complex<double>,1> y0(n), y1(n);
// Declare repetition counter
static int count = 0;
// Save initial data
double x0 = x;
y0 = y;
// Take full step
// rk4_fixed (x, y, rhs_eval, h);
rk4_fixed (x, y, rhs_eval, h, m, state, U, Jt,periodic);
// Save data
y1 = y;
// Restore initial data
x = x0;
y = y0;
// Take two half-steps
// rk4_fixed (x, y, rhs_eval, h/2.);
// rk4_fixed (x, y, rhs_eval, h/2.);
rk4_fixed (x, y, rhs_eval,
h/2, m, state,
U, Jt,periodic);
rk4_fixed (x, y, rhs_eval,
h/2, m, state,
U, Jt,periodic);
// Calculate truncation error
t_err = 0.;
double err, err1, err2;
if (flag == 0)
{
// Use absolute truncation error
for (int i = 0; i < n; i++)
{
// err = fabs (y(i) - y1(i));
err = abs (y(i) - y1(i));
t_err = (err > t_err) ? err : t_err;
}
}
else if (flag == 1)
{
// Use relative truncation error
for (int i = 0; i < n; i++)
{
// err = fabs ((y(i) - y1(i)) / y(i));
err = abs (y(i) - y1(i));
t_err = (err > t_err) ? err : t_err;
}
}
else
{
// Use mixed truncation error
for (int i = 0; i < n; i++)
{
// err1 = fabs ((y(i) - y1(i)) / y(i));
// err2 = fabs (y(i) - y1(i));
err1 = abs ((y(i) - y1(i)) / y(i));
err2 = abs (y(i) - y1(i));
err = (err1 < err2) ? err1 : err2;
t_err = (err > t_err) ? err : t_err;
}
}
// Prevent small truncation error from rounding to zero
if (t_err == 0.) t_err = 1.e-15;
// Calculate new step-length
double h_est = h * pow (fabs (acc / t_err), 0.2);
// Prevent step-length from changing by more than factor S
if (h_est / h > S)
h *= S;
else if (h_est / h < 1. / S)
h /= S;
else
h = h_est;
// Prevent step-length from exceeding h_max
h = (fabs(h) > h_max) ? h_max * h / fabs(h) : h;
// Abort if step-length falls below h_min
if (fabs(h) < h_min)
{
printf ("Error - |h| < hmin\n");
exit (1);
}
// If truncation error acceptable take step
if ((t_err <= acc) || (count >= maxrept))
{
rept = count;
count = 0;
}
// If truncation error unacceptable repeat step
else
{
count++;
x = x0;
y = y0;
rk4_adaptive (x, y, rhs_eval, h, t_err, acc,
S, rept, maxrept, h_min, h_max, flag,m,state,U,Jt,periodic);
}
return;
}
//this function provides the RHS for the homogenous DNLS with an
// extra term in the tunnelling part
void rhs_eval2 (double x, Array<complex<double>,1> y,
Array<complex<double>,1>& dydx,double U, double Jt,
int m,Array<complex<double>,1> state,int periodic)
{
int end=(state.extent(0)-1);
typedef complex<double> dcomp;
dcomp invi=dcomp(0,-1);
if (m>end)
{
cout<<"m is too large"<<endl;
}
if (m<0)
{
cout<<"m is too small"<<endl;
}
if (m==0)
{
if (periodic==0)
{
// cout<<"Periodic is zero"<<endl;
dydx(0)=invi*(Jt*(state(m+1)-2.0*y(0))+
U*abs(y(0))*abs(y(0))*y(0));
}
else if (periodic==1)
{
dydx(0)=invi*(Jt*(state(m+1)+state(end)-2.0*y(0))+
U*abs(y(0))*abs(y(0))*y(0));
}
else
{
cout<<"mistake in chosing the boundary conditions "<<endl;
}
}
// In these two cases I impose periodic boundary conditions
else if (m==end)
{
if (periodic == 0)
{
dydx(0)=invi*(Jt*(state(m-1)-2.0*y(0))+
U*abs(y(0))*abs(y(0))*y(0));
}
else if (periodic ==1)
{
dydx(0)=invi*(Jt*(state(m-1)+state(0)-2.0*y(0))+
U*abs(y(0))*abs(y(0))*y(0));
}
else
{
cout<<"mistake in chosing the boundary conditions for end "<<endl;
}
//cout<<"dydx(0) for m=end is "<<dydx(0)<<endl;
}
else
{
dydx(0)=invi*(Jt*(state(m-1)+state(m+1)-2.0*y(0))+
U*abs(y(0))*abs(y(0))*y(0));
}
}
// Once I introduce rhs_eval2, then I do not have to modify the appropriate RK
// functions (I just have to call them with the new rhs_eval2
/*
// this function now solves the homogeneous DNLS in the lattice using the previous function
// providing the RHS.
void rk4_fixed2 (double& x, Array<complex<double>,1>& y,
void (*rhs_eval2)(double, Array<complex<double>,1>,
Array<complex<double>,1>& ,double,
double,int ,Array<complex<double>,1>, int),
double h, int m, Array<complex<double>,1> state,
double U, double Jt,int periodic)
{
// Array y assumed to be of extent n, where n is no. of coupled
// equations
int n = y.extent(0);
int end=(state.extent(0)-1); //final well in the initial array
// Declare local arrays
// Array<double,1> k1(n), k2(n), k3(n), k4(n), f(n), dydx(n);
Array<complex<double>,1> k1(n), k2(n), k3(n), k4(n), f(n), dydx(n);
// Zeroth intermediate step
(*rhs_eval2) (x, y, dydx,U,Jt,m,state,periodic);
for (int j = 0; j < n; j++)
{
k1(j) = h * dydx(j);
f(j) = y(j) + k1(j) / 2.;
}
// First intermediate step
(*rhs_eval2) (x + h / 2., f, dydx,U,Jt,m,state,periodic);
for (int j = 0; j < n; j++)
{
k2(j) = h * dydx(j);
f(j) = y(j) + k2(j) / 2.;
}
// Second intermediate step
(*rhs_eval2) (x + h / 2., f, dydx,U,Jt,m,state,periodic);
for (int j = 0; j < n; j++)
{
k3(j) = h * dydx(j);
f(j) = y(j) + k3(j);
}
// Third intermediate step
(*rhs_eval2) (x + h, f, dydx,U,Jt,m,state,periodic);
for (int j = 0; j < n; j++)
{
k4(j) = h * dydx(j);
}
// Actual step
for (int j = 0; j < n; j++)
{
y(j) += k1(j) / 6. + k2(j) / 3. + k3(j) / 3. + k4(j) / 6.;
}
if (m==end)
{
x += h; //I update time only once all the initial state has been evolved
}
return;
}
void rk4_adaptive2 (double& x, Array<complex<double>,1>& y,
void (*rhs_eval2)(double, Array<complex<double>,1>,
Array<complex<double>,1>& ,double,
double,int ,Array<complex<double>,1>, int),
double& h, double& t_err, double acc,
double S, int& rept, int maxrept,
double h_min, double h_max, int flag,
int m, Array<complex<double>,1> state,
double U, double Jt, int periodic )
{
// Array y assumed to be of extent n, where n is no. of coupled
// equations
int n = y.extent(0);
// Declare local arrays
// Array<double,1> y0(n), y1(n);
Array<complex<double>,1> y0(n), y1(n);
// Declare repetition counter
static int count = 0;
// Save initial data
double x0 = x;
y0 = y;
// Take full step
// rk4_fixed (x, y, rhs_eval, h);
rk4_fixed (x, y, rhs_eval2, h, m, state, U, Jt,periodic);
// Save data
y1 = y;
// Restore initial data
x = x0;
y = y0;
// Take two half-steps
// rk4_fixed (x, y, rhs_eval, h/2.);
// rk4_fixed (x, y, rhs_eval, h/2.);
rk4_fixed (x, y, rhs_eval2,
h/2, m, state,
U, Jt,periodic);
rk4_fixed (x, y, rhs_eval2,
h/2, m, state,
U, Jt,periodic);
// Calculate truncation error
t_err = 0.;
double err, err1, err2;
if (flag == 0)
{
// Use absolute truncation error
for (int i = 0; i < n; i++)
{
// err = fabs (y(i) - y1(i));
err = abs (y(i) - y1(i));
t_err = (err > t_err) ? err : t_err;
}
}
else if (flag == 1)
{
// Use relative truncation error
for (int i = 0; i < n; i++)
{
// err = fabs ((y(i) - y1(i)) / y(i));
err = abs (y(i) - y1(i));
t_err = (err > t_err) ? err : t_err;
}
}
else
{
// Use mixed truncation error
for (int i = 0; i < n; i++)
{
// err1 = fabs ((y(i) - y1(i)) / y(i));
// err2 = fabs (y(i) - y1(i));
err1 = abs ((y(i) - y1(i)) / y(i));
err2 = abs (y(i) - y1(i));
err = (err1 < err2) ? err1 : err2;
t_err = (err > t_err) ? err : t_err;
}
}
// Prevent small truncation error from rounding to zero
if (t_err == 0.) t_err = 1.e-15;
// Calculate new step-length
double h_est = h * pow (fabs (acc / t_err), 0.2);
// Prevent step-length from changing by more than factor S
if (h_est / h > S)
h *= S;
else if (h_est / h < 1. / S)
h /= S;
else
h = h_est;
// Prevent step-length from exceeding h_max
h = (fabs(h) > h_max) ? h_max * h / fabs(h) : h;
// Abort if step-length falls below h_min
if (fabs(h) < h_min)
{
printf ("Error - |h| < hmin\n");
exit (1);
}
// If truncation error acceptable take step
if ((t_err <= acc) || (count >= maxrept))
{
rept = count;
count = 0;
}
// If truncation error unacceptable repeat step
else
{
count++;
x = x0;
y = y0;
rk4_adaptive (x, y, rhs_eval2, h, t_err, acc,
S, rept, maxrept, h_min, h_max, flag,m,state,U,Jt,periodic);
}
return;
}
*/
double real_integration(double f[],int NN,double length)
{
double integral=0;
double h=length/(NN-1);
for (int i=0;i<(NN-4);i=i+4)
{
double S2=h/45.0*(14.0*f[i]+
64.0*f[i+1] +
24.0*f[i+2]+
64.0*f[i+3]+
14.0*f[i+4]);
integral = integral + S2;
}
return integral;
}
double complex_integration (double f[],int NN,double length)
{
double integral=0;
double h=length/(NN-1);
for (int i=0;i<(NN-4);i=i+4)
{
double S2=h/45.0*(14.0*REAL(f,i)+
64.0*REAL(f,(i+1)) +
24.0*REAL(f,(i+2))+
64.0*REAL(f,(i+3))+
14.0*REAL(f,(i+4)));
double S1=h/45.0*(14.0*IMAG(f,i)+
64.0*IMAG(f,(i+1)) +
24.0*IMAG(f,(i+2))+
64.0*IMAG(f,(i+3))+
14.0*IMAG(f,(i+4)));
integral = integral + S2+S1;
}
return integral;
}
void my_save(char fileout[80],int NN, double arr[])
{
ofstream my_output;
my_output.open(fileout);
for(int i=0; i<NN; i++)
{
my_output<<arr[i]<<endl;
}
my_output.close( );
}
void my_save_vec(char fileout[80],int NN, vector<double> arr)
{
ofstream my_output;
my_output.open(fileout);
for(int i=0; i<NN; i++)
{
my_output<<arr[i]<<endl;
}
my_output.close( );
}
void my_density(double arr[],int NN,double dens[])
{
for (int i=0;i<NN;i++)
{
dens[i]=REAL(arr,i)*REAL(arr,i)+IMAG(arr,i)*IMAG(arr,i);
}
}
void calc_cent_site(double Jt, double U, Array<complex<double>,1> inistate, vector<double> & Cinteraction,vector<double> & Ckin)
{int nsites=inistate.extent(0);
double onehalf=0.5;
int middle=(nsites-1)/2; // remember that the first element in an array is arr(0) and not arr(1)
double number = onehalf*U*abs(inistate(middle))*abs(inistate(middle))
*abs(inistate(middle))*abs(inistate(middle));
Cinteraction.push_back(number);
number=-Jt*abs(inistate(middle)-inistate(middle-1))
*abs(inistate(middle)-inistate(middle-1))
-Jt*abs(inistate(middle+1)-inistate(middle))
*abs(inistate(middle+1)-inistate(middle))
;
Ckin.push_back(number);
}
void calc_lattice(double Jt, double U, Array<complex<double>,1> inistate,vector<double> & Tinteraction,vector<double> & Tkin)
{
int nsites=inistate.extent(0);
double onehalf=0.5;
double number1 =0;
double number2 =0 ;
for (int i=0;i<nsites;i++)
{
number1=number1+onehalf*U*abs(inistate(i))*abs(inistate(i))*abs(inistate(i))*abs(inistate(i));
// the previous expression accounts for the total interaction
if (i == 0)
{
number2=number2-Jt*abs(inistate(i)-inistate(nsites))*abs(inistate(i)-inistate(nsites));
}
else
{
number2=number2-Jt*abs(inistate(i)-inistate(i-1))*abs(inistate(i)-inistate(i-1));
}
}
Tinteraction.push_back(number1);
Tkin.push_back(number2);
}
void FFT(double data[],int NN,double V[],double delta_t,double length )
{
//first I'll do the evolution in real space for delta_t/2
for (int i=0;i<NN;i++)
{
double alpha=V[i]*delta_t/2.0;
double tempr=cos(alpha)*REAL(data,i)+sin(alpha)*IMAG(data,i);
double tempi=cos(alpha)*IMAG(data,i)-sin(alpha)*REAL(data,i);
REAL(data,i)=tempr;
IMAG(data,i)=tempi;
}
// this way I have made the first part of the time evolution
// Now it is time for the FFT
gsl_fft_complex_radix2_forward (data, 1, NN);
double pi=acos(-1.0);
double dk=2.0*pi/length;
double mom[NN];
for (int i=0;i<NN;i++)
{
mom[i]=dk*i;
if (i>(NN/2))
{
mom[i]=dk*(i-NN);
}
}
//evolution in momentum space
for (int i=0;i<NN;i++)
{
double alpha=mom[i]*mom[i]*delta_t/2;
double tempr=cos(alpha)*REAL(data,i)+sin(alpha)*IMAG(data,i);
double tempi=cos(alpha)*IMAG(data,i)-sin(alpha)*REAL(data,i);
REAL(data,i)=tempr;
IMAG(data,i)=tempi;
}
// now FFT back to the real space
gsl_fft_complex_radix2_inverse (data, 1, NN);
// last part of the evolution under the effect of the potential
//in real space
for (int i=0;i<NN;i++)
{
double alpha=V[i]*delta_t/2.0;
double tempr=cos(alpha)*REAL(data,i)+sin(alpha)*IMAG(data,i);
double tempi=cos(alpha)*IMAG(data,i)-sin(alpha)*REAL(data,i);
REAL(data,i)=tempr;
IMAG(data,i)=tempi;
}
}
void nonlinFFT(double data[],int NN,double V[],double delta_t,double length,double G )
{
//first I'll do the evolution in real space for delta_t/2
for (int i=0;i<NN;i++)
{
double alpha=V[i]*delta_t/2.0+G*delta_t/2.0*(REAL(data,i)*REAL(data,i)+IMAG(data,i)*IMAG(data,i));
double tempr=cos(alpha)*REAL(data,i)+sin(alpha)*IMAG(data,i);
double tempi=cos(alpha)*IMAG(data,i)-sin(alpha)*REAL(data,i);
REAL(data,i)=tempr;
IMAG(data,i)=tempi;
}
// this way I have made the first part of the time evolution
// Now it is time for the FFT
gsl_fft_complex_radix2_forward (data, 1, NN);
double pi=acos(-1.0);
double dk=2.0*pi/length;
double mom[NN];
for (int i=0;i<NN;i++)
{
mom[i]=dk*i;
if (i>(NN/2))
{
mom[i]=dk*(i-NN);
}
}
//evolution in momentum space
for (int i=0;i<NN;i++)
{
double alpha=mom[i]*mom[i]*delta_t/2;
double tempr=cos(alpha)*REAL(data,i)+sin(alpha)*IMAG(data,i);
double tempi=cos(alpha)*IMAG(data,i)-sin(alpha)*REAL(data,i);
REAL(data,i)=tempr;
IMAG(data,i)=tempi;
}
// now FFT back to the real space
gsl_fft_complex_radix2_inverse (data, 1, NN);
// last part of the evolution under the effect of the potential
//in real space
for (int i=0;i<NN;i++)
{
double alpha=V[i]*delta_t/2.0+G*delta_t/2.0*(REAL(data,i)*REAL(data,i)+IMAG(data,i)*IMAG(data,i));
double tempr=cos(alpha)*REAL(data,i)+sin(alpha)*IMAG(data,i);
double tempi=cos(alpha)*IMAG(data,i)-sin(alpha)*REAL(data,i);
REAL(data,i)=tempr;
IMAG(data,i)=tempi;
}
}
void imaginaryFFT(double data[],int NN,double V[],double delta_t,double length,double G )
{
//first I'll do the evolution in real space for delta_t/2
for (int i=0;i<NN;i++)
{
double alpha=V[i]*delta_t/2.0+G*delta_t/2.0*(REAL(data,i)*REAL(data,i)+IMAG(data,i)*IMAG(data,i));
double tempr=exp(-alpha)*REAL(data,i);
double tempi=exp(-alpha)*IMAG(data,i);
REAL(data,i)=tempr;
IMAG(data,i)=tempi;
}
// this way I have made the first part of the time evolution
// Now it is time for the FFT
gsl_fft_complex_radix2_forward (data, 1, NN);
double pi=acos(-1.0);
double dk=2.0*pi/length;
double mom[NN];
for (int i=0;i<NN;i++)
{
mom[i]=dk*i;
if (i>(NN/2))
{
mom[i]=dk*(i-NN);
}
}
//evolution in momentum space
for (int i=0;i<NN;i++)
{
double alpha=mom[i]*mom[i]*delta_t/2;
double tempr=exp(-alpha)*REAL(data,i);
double tempi=exp(-alpha)*IMAG(data,i);
REAL(data,i)=tempr;
IMAG(data,i)=tempi;
}
// now FFT back to the real space
gsl_fft_complex_radix2_inverse (data, 1, NN);
// last part of the evolution under the effect of the potential
//in real space
for (int i=0;i<NN;i++)
{
double alpha=V[i]*delta_t/2.0+G*delta_t/2.0*(REAL(data,i)*REAL(data,i)+IMAG(data,i)*IMAG(data,i));
double tempr=exp(-alpha)*REAL(data,i);
double tempi=exp(-alpha)*IMAG(data,i);
REAL(data,i)=tempr;
IMAG(data,i)=tempi;
}
}
void renormalization(double g[],int NN,double length)
{
double f[NN];
for (int i=0;i<NN; i++)
{
f[i]=REAL(g,i)*REAL(g,i)+IMAG(g,i)*IMAG(g,i);
}
double integral=0;
double h=length/(NN-1);
for (int i=0;i<(NN-4);i=i+4)
{
double S2=h/45.0*(14.0*f[i]+
64.0*f[i+1] +
24.0*f[i+2]+
64.0*f[i+3]+
14.0*f[i+4]);
integral = integral + S2;
}
for (int i=0;i<NN;i++)
{
REAL(g,i)=REAL(g,i)/sqrt(integral);
IMAG(g,i)=IMAG(g,i)/sqrt(integral);
}
}
void renorm_imaginaryFFT(double data[],int NN,double V[],double delta_t,double length,double G )
{
//first I'll do the evolution in real space for delta_t/2
for (int i=0;i<NN;i++)
{
double alpha=V[i]*delta_t/2.0+G*delta_t/2.0*(REAL(data,i)*REAL(data,i)+IMAG(data,i)*IMAG(data,i));
double tempr=exp(-alpha)*REAL(data,i);
double tempi=exp(-alpha)*IMAG(data,i);
REAL(data,i)=tempr;
IMAG(data,i)=tempi;
}
// this way I have made the first part of the time evolution
// Now it is time for the FFT
gsl_fft_complex_radix2_forward (data, 1, NN);
double pi=acos(-1.0);
double dk=2.0*pi/length;
double mom[NN];
for (int i=0;i<NN;i++)
{
mom[i]=dk*i;
if (i>(NN/2))
{
mom[i]=dk*(i-NN);
}
}
//evolution in momentum space
for (int i=0;i<NN;i++)
{
double alpha=mom[i]*mom[i]*delta_t/2;
double tempr=exp(-alpha)*REAL(data,i);
double tempi=exp(-alpha)*IMAG(data,i);
REAL(data,i)=tempr;
IMAG(data,i)=tempi;
}
// now FFT back to the real space
gsl_fft_complex_radix2_inverse (data, 1, NN);
// last part of the evolution under the effect of the potential
//in real space
for (int i=0;i<NN;i++)
{
double alpha=V[i]*delta_t/2.0+G*delta_t/2.0*(REAL(data,i)*REAL(data,i)+IMAG(data,i)*IMAG(data,i));
double tempr=exp(-alpha)*REAL(data,i);
double tempi=exp(-alpha)*IMAG(data,i);
REAL(data,i)=tempr;
IMAG(data,i)=tempi;
}
renormalization(data, NN, length);
}
void imaginary_evol_FFT(double data[],int NN,double V[],double delta_t,double length,double G,double time )
{
// this routine is OK to evolve a system in imaginary time (it already has the iteration
// process in-built )
// Now I find out how many times I want to iterate the FFT routine for time evolution
/*
static cast example below
double rate =
static_cast<double>(total)/days;
*/
int Niter = static_cast<int>(time/delta_t);
for (int i=0;i<Niter;i++)
{
renorm_imaginaryFFT(data, NN,V, delta_t,length, G );
}
}
// now a similar function for real-time propagation and a static potential
void real_evol_FFT(double data[],int NN,double V[],double delta_t,double length,double G,double time )
{
// this routine is OK to evolve a system in imaginary time (it already has the iteration
// process in-built )
// Now I find out how many times I want to iterate the FFT routine for time evolution
/*
static cast example below
double rate =
static_cast<double>(total)/days;
*/
int Niter = static_cast<int>(time/delta_t);
for (int i=0;i<Niter;i++)
{
nonlinFFT(data, NN,V, delta_t,length, G );
}
}