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
|
#include "fourier.hpp"
#include <getopt.h>
#include <iostream>
#include <fstream>
int main(int argc, char* argv[]) {
unsigned p = 2;
unsigned s = 2;
Real λ = 0.5;
Real τ₀ = 0;
Real y₀ = 0;
Real yₘₐₓ = 0.5;
Real Δy = 0.05;
unsigned log2n = 8;
Real τₘₐₓ = 20;
int opt;
while ((opt = getopt(argc, argv, "p:s:2:T:t:0:y:d:")) != -1) {
switch (opt) {
case 'p':
p = atoi(optarg);
break;
case 's':
s = atoi(optarg);
break;
case '2':
log2n = atoi(optarg);
break;
case 'T':
τₘₐₓ = atof(optarg);
break;
case 't':
τ₀ = atof(optarg);
break;
case '0':
y₀ = atof(optarg);
break;
case 'y':
yₘₐₓ = atof(optarg);
break;
case 'd':
Δy = atof(optarg);
break;
default:
exit(1);
}
}
unsigned n = pow(2, log2n);
Real Δτ = τₘₐₓ / M_PI / n;
Real Δω = M_PI / τₘₐₓ;
Real y = y₀;
FourierTransform fft(n, Δω, Δτ, FFTW_ESTIMATE);
while (y += Δy, y <= yₘₐₓ) {
std::vector<Real> C(2 * n);
std::vector<Real> R(2 * n);
std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
if (cfile.is_open()) {
cfile.read((char*)(C.data()), C.size() * sizeof(Real));
cfile.close();
std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
rfile.read((char*)(R.data()), R.size() * sizeof(Real));
rfile.close();
Real e = energy(C, R, p, s, λ, y, Δτ);
std::vector<Complex> Ct = fft.fourier(C);
std::vector<Complex> Rt = fft.fourier(R);
auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
Real z = ((std::conj(Rt[0]) + pow(y, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real();
std::cout << y << " " << e << " " << Ct[0].real() << " " << z << std::endl;
}
}
return 0;
}
|
