1 files changed, 85 insertions, 0 deletions
diff --git a/integrator.cpp b/integrator.cpp
new file mode 100644
index 0000000..bb22a36
--- /dev/null
+++ b/integrator.cpp
@@ -0,0 +1,85 @@
+#include <getopt.h>
+#include <vector>
+#include <array>
+#include <cmath>
+#include <iostream>
+
+using Real = double;
+
+class Point : public std::array<Real, 2> {
+public:
+  Real τ() const {
+    return front();
+  }
+  Real C() const {
+    return back();
+  }
+};
+
+Real f(Real q) {
+  return 0.5 * pow(q, 2);
+}
+
+Real df(Real q) {
+  return q;
+}
+
+Real ddf(Real q) {
+  return 1;
+}
+
+Real integrate(const std::vector<Point>& Cₜ) {
+  Real I = 0;
+  for (unsigned i = 0; i < Cₜ.size() - 1; i++) {
+    Real Δτ = Cₜ[i + 1].τ() - Cₜ[i].τ();
+    Real C = (Cₜ[i + 1].C() + Cₜ[i].C()) / 2;
+    Real dC = (Cₜ[i + 1].C() - Cₜ[i].C()) / Δτ;
+    I += Δτ * df(C) * dC;
+  }
+  return I;
+}
+
+int main(int argc, char* argv[]) {
+  Real Δτ = 1e-3;
+  Real τₘₐₓ = 1e3;
+  Real τ₀ = 0;
+  Real y = 0.5;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "d:T:t:y:")) != -1) {
+    switch (opt) {
+    case 'd':
+      Δτ = atof(optarg);
+      break;
+    case 'T':
+      τₘₐₓ = atof(optarg);
+      break;
+    case 't':
+      τ₀ = atof(optarg);
+      break;
+    case 'y':
+      y = atof(optarg);
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  std::vector<Point> Cₜ;
+  Cₜ.reserve(τₘₐₓ / Δτ);
+
+  Cₜ.push_back({0, 1});
+  Cₜ.push_back({Δτ, 1 - Δτ});
+
+  while (Cₜ.back().τ() < τₘₐₓ) {
+    Real dC = -Cₜ.back().C() - 2 * pow(y, 2) * integrate(Cₜ);
+    Cₜ.push_back({Cₜ.back().τ() + Δτ, Cₜ.back().C() + Δτ * dC});
+  }
+
+  for (const Point& p : Cₜ) {
+    std::cout << p.τ() << " " << p.C() << std::endl;
+  }
+
+  return 0;
+}