Added mixed support

1 files changed, 24 insertions, 7 deletions

diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
index 829a074..dd1b2d3 100644
--- a/fourier_integrator.cpp
+++ b/fourier_integrator.cpp
@@ -9,18 +9,30 @@ using Real = double;
 using Complex = std::complex<Real>;
 using namespace std::complex_literals;
 
-inline Real f(unsigned p, Real q) {
+inline Real fP(unsigned p, Real q) {
   return 0.5 * pow(q, p);
 }
 
-inline Real df(unsigned p, Real q) {
+inline Real dfP(unsigned p, Real q) {
   return 0.5 * p * pow(q, p - 1);
 }
 
-inline Real ddf(unsigned p, Real q) {
+inline Real ddfP(unsigned p, Real q) {
   return 0.5 * p * (p - 1) * pow(q, p - 2);
 }
 
+inline Real f(Real λ, unsigned p, unsigned s, Real q) {
+  return (1 - λ) * fP(p, q) + λ * fP(s, q);
+}
+
+inline Real df(Real λ, unsigned p, unsigned s, Real q) {
+  return (1 - λ) * dfP(p, q) + λ * dfP(s, q);
+}
+
+inline Real ddf(Real λ, unsigned p, unsigned s, Real q) {
+  return (1 - λ) * ddfP(p, q) + λ * ddfP(s, q);
+}
+
 class FourierTransform {
 private:
   std::vector<Real> a;
@@ -64,6 +76,8 @@ public:
 
 int main(int argc, char* argv[]) {
   unsigned p = 2;
+  unsigned s = 2;
+  Real λ = 0.5;
   Real τ₀ = 0;
   Real yₘₐₓ = 0.5;
   Real Δy = 0.05;
@@ -77,11 +91,14 @@ int main(int argc, char* argv[]) {
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "p:2:T:t:y:d:I:g:")) != -1) {
+  while ((opt = getopt(argc, argv, "p:s:2:T:t:y:d:I:g:")) != -1) {
     switch (opt) {
     case 'p':
       p = atoi(optarg);
       break;
+    case 's':
+      s = atoi(optarg);
+      break;
     case '2':
       log2n = atoi(optarg);
       break;
@@ -142,13 +159,13 @@ int main(int argc, char* argv[]) {
       it++;
       std::vector<Real> RddfC(C.size());
       for (unsigned i = 0; i < C.size(); i++) {
-        RddfC[i] = R[i] * ddf(p, C[i]);
+        RddfC[i] = R[i] * ddf(λ, p, s, 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(p, C[i]);
+        dfC[i] = df(λ, p, s, C[i]);
       }
       std::vector<Complex> dfCt = fft.fourier(dfC);
 
@@ -184,7 +201,7 @@ int main(int argc, char* argv[]) {
       Real energy = 0;
 
       for (unsigned i = 0; i < n; i++) {
-        energy += y * R[i] * df(p, C[i]) * M_PI * Δτ;
+        energy += y * R[i] * df(λ, p, s, C[i]) * M_PI * Δτ;
       }
 
       std::cerr << "y " << y << " iteration " << it << ": " << sqrt(ΔC / C.size()) << " " << pow(y, 2) / z << " " << energy << std::endl;