-
Notifications
You must be signed in to change notification settings - Fork 0
/
BEC_Energy.cpp
executable file
·129 lines (106 loc) · 3.55 KB
/
BEC_Energy.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
#include"fftw3.h"
#include<vector>
#include<iostream>
#include"headers/split_step.h"
#include"headers/bec.h"
#include<algorithm>
#include<functional>
#include<iterator>
#include<cmath>
#include<complex>
#include"headers/BEC_Groundstate.h"
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>
using namespace std;
int main(int argc, char* argv[]){
typedef std::vector<double> Vector;
typedef std::vector<complex<double> > Cvector;
typedef boost::numeric::ublas::matrix<double> Matrix;
typedef complex<double> complex_type;
typedef fftw_complex complex_c_type;
typedef vector<complex<double> >::size_type size_t;
typedef complex<double> complex_type;
typedef fftw_complex complex_c_type;
double K = 0.1;
double g = 0.1;
double maxt = 20;
int period = 10000;
const double L = 10;
double mynorm= 0;
int N=1024;
const double deltat = 2.*pi/double(period);
const complex_type I(0,1);
double t=0;
int stepcounter=0;
//parsing command line
switch(argc){
case 1:
break;
case 2:
K=atof(argv[1]);
break;
case 3:
K=atof(argv[1]);
g=atof(argv[2]);
break;
case 4:
K=atof(argv[1]);
g=atof(argv[2]);
maxt=atof(argv[3]);
break;
case 5:
K=atof(argv[1]);
g=atof(argv[2]);
maxt=atof(argv[3]);
period = atoi(argv[4]);
break;
case 6:
K=atof(argv[1]);
g=atof(argv[2]);
maxt=atof(argv[3]);
period = atoi(argv[4]);
N=atoi(argv[5]);
break;
default:
break;
}
//finding Groundstate
Vector groundstate(N, 100.);
Cvector data(N);
Cvector mom(N);
for(int i = 0; i<N; ++i){
double xpos = L*0.5-double(i)*L/double(N);
data[i]=1./sqrt(sqrt(pi) )*exp(-xpos*xpos/2.);
}
for(int i = 0; i<20000;++i){
data=split_step::S3imag(L, deltat, g, data, bec::Potential<double, true, true>());
mynorm = L_2_norm(data.begin(), data.end(), L);
//cout<<mynorm<<endl;
transform(data.begin(), data.end(), data.begin(), bind2nd(divides<complex<double> >(), mynorm) );
}
//properly normalize the whole thing, so that we may use it with the propagator.
mynorm = L_2_norm(data.begin(), data.end(), L);
transform(data.begin(), data.end(), data.begin(), bind2nd(divides<complex<double> >(), mynorm) );
cout<<"System size"<<'\t'<<"Points"<<'\t'<<"Coupling"<<'\t'<<"deltat"<<'\t'<<"Kick strength"<<'\t'<<"Max time"<<'\t'<<"Period"<<endl;
cout<<L<<'\t'<<N<<'\t'<<g<<'\t'<<deltat<<'\t'<<K<<'\t'<<maxt<<'\t'<<period<<endl<<endl;
cout<<"t"<<'\t'<<"Energy"<<endl;
// mom=FourierTransform<fft::forward, fft::Estimate>(data);
cout<<t<<'\t'<<bec::energy(data, bec::Potential<double, true, true>(), L, g) <<endl;
vector<complex_type> kick(N);
for (size_t i = 0; i<N; ++i) {
double xpos = L*0.5-double(i)*L/double(N);
kick[i]= exp(-I*K*sin(xpos));
}
transform(data.begin(), data.end(), kick.begin(), data.begin(), multiplies<complex_type>() );
for(;t<maxt; t+=deltat, ++stepcounter){
data=split_step::S3(L, deltat, g, data, bec::Potential<double, true, true>());
mynorm = L_2_norm(data.begin(), data.end(), L);
transform(data.begin(), data.end(), data.begin(), bind2nd(divides<complex<double> >(), mynorm) );
//mom=FourierTransform<fft::forward, fft::Estimate>(data);
if(!(stepcounter%1000)){
cout<<t<<'\t'<<bec::energy(data, bec::Potential<double, true, true>(), L, g)<<endl;
}
}
return 0;
}