Write and read from binary files

1 files changed, 46 insertions, 20 deletions

diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
index db5e07e..3c33946 100644
--- a/fourier_integrator.cpp
+++ b/fourier_integrator.cpp
@@ -4,6 +4,7 @@
 #include <iostream>
 #include <fftw3.h>
 #include <complex>
+#include <fstream>
 
 using Real = double;
 using Complex = std::complex<Real>;
@@ -79,19 +80,22 @@ int main(int argc, char* argv[]) {
   unsigned s = 2;
   Real λ = 0.5;
   Real τ₀ = 0;
+  Real y₀ = 0;
   Real yₘₐₓ = 0.5;
   Real Δy = 0.05;
 
   unsigned log2n = 8;
-  Real τₘₐₓ = 20 / M_PI;
+  Real τₘₐₓ = 20;
 
   unsigned maxIterations = 1000;
   Real ε = 1e-14;
   Real γ = 0;
 
+  bool loadData = false;
+
   int opt;
 
-  while ((opt = getopt(argc, argv, "p:s:2:T:t:y:d:I:g:")) != -1) {
+  while ((opt = getopt(argc, argv, "p:s:2:T:t:0:y:d:I:g:l")) != -1) {
     switch (opt) {
     case 'p':
       p = atoi(optarg);
@@ -103,11 +107,14 @@ int main(int argc, char* argv[]) {
       log2n = atoi(optarg);
       break;
     case 'T':
-      τₘₐₓ = atof(optarg) / M_PI;
+      τₘₐₓ = atof(optarg);
       break;
     case 't':
       τ₀ = atof(optarg);
       break;
+    case '0':
+      y₀ = atof(optarg);
+      break;
     case 'y':
       yₘₐₓ = atof(optarg);
       break;
@@ -120,6 +127,9 @@ int main(int argc, char* argv[]) {
     case 'g':
       γ = atof(optarg);
       break;
+    case 'l':
+      loadData = true;
+      break;
     default:
       exit(1);
     }
@@ -127,8 +137,8 @@ int main(int argc, char* argv[]) {
 
   unsigned n = pow(2, log2n);
 
-  Real Δτ = τₘₐₓ / n;
-  Real Δω = 1 / τₘₐₓ;
+  Real Δτ = τₘₐₓ / M_PI / n;
+  Real Δω = M_PI / τₘₐₓ;
 
   Real z = (-1+sqrt(1+2*τ₀)) / (2 * τ₀);
   Real Γ₀ = 1;
@@ -138,19 +148,30 @@ int main(int argc, char* argv[]) {
 
   FourierTransform fft(n, Δω, Δτ);
 
-  for (unsigned i = 0; i < n; i++) {
-    Real τ = i * Δτ * M_PI;
-    C[i] = Γ₀ / 2 * (exp(-z * τ) - z * τ₀ * exp(-τ / τ₀)) / (z - pow(z, 3) * pow(τ₀, 2));
-    if (i > 0) {
-      C[2 * n - i] = C[i];
+  if (!loadData) {
+    // start from the exact solution for τ₀ = 0
+    for (unsigned i = 0; i < n; i++) {
+      Real τ = i * Δτ * M_PI;
+      C[i] = Γ₀ / 2 * (exp(-z * τ) - z * τ₀ * exp(-τ / τ₀)) / (z - pow(z, 3) * pow(τ₀, 2));
+      if (i > 0) {
+        C[2 * n - i] = C[i];
+      }
+      R[i] = exp(-z * τ);
     }
-    R[i] = exp(-z * τ);
+  } else {
+    std::string file_end = std::to_string(p) + "_" + std::to_string(s) + "_" + std::to_string(λ) + "_" + std::to_string(τ₀) + "_" + std::to_string(y₀) + "_" + std::to_string(log2n) + "_" + std::to_string(τₘₐₓ) + ".dat";
+    std::ifstream cfile("C_"+file_end, std::ios::binary);
+    cfile.read((char*)(C.data()), C.size() * sizeof(Real));
+    cfile.close();
+    std::ifstream rfile("R_"+file_end, std::ios::binary);
+    rfile.read((char*)(R.data()), R.size() * sizeof(Real));
+    rfile.close();
   }
 
   std::vector<Complex> Ct = fft.fourier(C);
   std::vector<Complex> Rt = fft.fourier(R);
 
-  Real y = 0;
+  Real y = y₀;
 
   while (y += Δy, y <= yₘₐₓ) {
     Real ΔC = 1;;
@@ -181,6 +202,9 @@ int main(int argc, char* argv[]) {
 
       std::vector<Real> Cnew = fft.inverse(Ct);
       std::vector<Real> Rnew = fft.inverse(Rt);
+      for (unsigned i = n; i < 2 * n; i++) {
+        Rnew[i] = 0;
+      }
 
       ΔC = 0;
       for (unsigned i = 0; i < Cnew.size(); i++) {
@@ -198,6 +222,7 @@ int main(int argc, char* argv[]) {
 
       z *= Cnew[0];
 
+      std::cerr << it << " " << ΔC << std::endl;
 
       if (it > maxIterations) {
         it = 0;
@@ -213,14 +238,15 @@ int main(int argc, char* argv[]) {
 
     std::cerr << "y " << y << " " << energy << std::endl;
 
-    for (unsigned i = 0; i < C.size(); i++) {
-      std::cout << i * Δτ * M_PI << " " << C[i] << " ";
-    }
-    std::cout << std::endl;
-    for (unsigned i = 0; i < R.size(); i++) {
-      std::cout << i * Δτ * M_PI << " " << R[i] << " ";
-    }
-    std::cout << std::endl;
+    std::string file_end = std::to_string(p) + "_" + std::to_string(s) + "_" + std::to_string(λ) + "_" + std::to_string(τ₀) + "_" + std::to_string(y) + "_" + std::to_string(log2n) + "_" + std::to_string(τₘₐₓ) + ".dat";
+
+    std::ofstream outfile("C_" + file_end, std::ios::out | std::ios::binary);
+    outfile.write((const char*)(C.data()), C.size() * sizeof(Real));
+    outfile.close();
+
+    std::ofstream outfileR("R_" + file_end, std::ios::out | std::ios::binary);
+    outfileR.write((const char*)(R.data()), R.size() * sizeof(Real));
+    outfileR.close();
   }
 
   return 0;