Added p flag and calculating energy

1 files changed, 30 insertions, 22 deletions

diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
index 9c97f1b..40a958e 100644
--- a/fourier_integrator.cpp
+++ b/fourier_integrator.cpp
@@ -9,17 +9,15 @@ using Real = double;
 using Complex = std::complex<Real>;
 using namespace std::complex_literals;
 
-unsigned p = 2;
-
-Real f(Real q) {
+inline Real f(unsigned p, Real q) {
   return 0.5 * pow(q, p);
 }
 
-Real df(Real q) {
+inline Real df(unsigned p, Real q) {
   return 0.5 * p * pow(q, p - 1);
 }
 
-Real ddf(Real q) {
+inline Real ddf(unsigned p, Real q) {
   return 0.5 * p * (p - 1) * pow(q, p - 2);
 }
 
@@ -65,21 +63,27 @@ public:
 };
 
 int main(int argc, char* argv[]) {
-  Real Δω = 1e-3;
-  Real Δτ = 1e-3;
+  unsigned p = 2;
   Real τ₀ = 0;
   Real y = 0.5;
+
+  unsigned log2n = 8;
+  Real τₘₐₓ = 20 / M_PI;
+
   unsigned iterations = 10;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "d:T:t:y:I:")) != -1) {
+  while ((opt = getopt(argc, argv, "p:2:T:t:y:I:")) != -1) {
     switch (opt) {
-    case 'd':
-      Δω = atof(optarg);
+    case 'p':
+      p = atoi(optarg);
+      break;
+    case '2':
+      log2n = atoi(optarg);
       break;
     case 'T':
-      Δτ = atof(optarg);
+      τₘₐₓ = atof(optarg) / M_PI;
       break;
     case 't':
       τ₀ = atof(optarg);
@@ -95,12 +99,14 @@ int main(int argc, char* argv[]) {
     }
   }
 
+  unsigned n = pow(2, log2n);
+
+  Real Δτ = τₘₐₓ / n;
+  Real Δω = 1 / τₘₐₓ;
+
   Real z = (-1+sqrt(1+2*τ₀)) / (2 * τ₀);
   Real Γ₀ = 1;
 
-  Real ωₘₐₓ = 1 / Δτ;
-  unsigned n = ωₘₐₓ / Δω;
-
   std::vector<Real> C(2 * n);
   std::vector<Real> R(2 * n);
 
@@ -121,13 +127,13 @@ int main(int argc, char* argv[]) {
   for (unsigned it = 0; it < iterations; it++) {
     std::vector<Real> RddfC(C.size());
     for (unsigned i = 0; i < C.size(); i++) {
-      RddfC[i] = R[i] * ddf(C[i]);
+      RddfC[i] = R[i] * ddf(p, C[i]);
     }
     std::vector<Complex> RddfCt = fft.fourier(RddfC);
 
     std::vector<Real> dfC(C.size());
     for (unsigned i = 0; i < C.size(); i++) {
-      dfC[i] = df(C[i]);
+      dfC[i] = df(p, C[i]);
     }
     std::vector<Complex> dfCt = fft.fourier(dfC);
 
@@ -164,15 +170,17 @@ int main(int argc, char* argv[]) {
     R = Rnew;
   }
 
-  for (unsigned i = 0; i < C.size(); i++) {
-    std::cout << i * Δτ * M_PI << " " << C[i] << std::endl;
+  Real energy = 0;
+
+  for (unsigned i = 0; i < n; i++) {
+    energy += y * R[i] * df(p, C[i]) * M_PI * Δτ;
   }
 
-  /*
-  for (unsigned i = 0; i < Ĉ.size(); i++) {
-    std::cout << i * Δω << " " << Ĉ[i].real() * Δω / (2 * M_PI) << std::endl;
+  std::cerr << energy << std::endl;
+
+  for (unsigned i = 0; i < C.size(); i++) {
+    std::cout << i * Δτ * M_PI << " " << C[i] << std::endl;
   }
 
-  */
   return 0;
 }