1 files changed, 44 insertions, 2 deletions
diff --git a/fourier.hpp b/fourier.hpp
index 86335f7..22a57e7 100644
--- a/fourier.hpp
+++ b/fourier.hpp
@@ -3,6 +3,12 @@
 #include <fftw3.h>
 #include <complex>
 #include <vector>
+#include <tuple>
+#include <gsl/gsl_sf_gamma.h>
+
+#ifndef FFTW_THREADS
+#define FFTW_THREADS 1
+#endif
 
 using Real = double;
 using Complex = std::complex<Real>;
@@ -14,19 +20,55 @@ Real ddf(Real λ, unsigned p, unsigned s, Real q);
 
 class FourierTransform {
 private:
-  std::vector<Real> a;
-  std::vector<Complex> â;
+  Real* a;
+  Complex* â;
   fftw_plan plan_r2c;
   fftw_plan plan_c2r;
+  unsigned n;
   Real Δω;
   Real Δτ;
 public:
   FourierTransform(unsigned n, Real Δω, Real Δτ, unsigned flags = 0);
   ~FourierTransform();
   std::vector<Complex> fourier(const std::vector<Real>& c);
+  std::vector<Complex> fourier();
   std::vector<Real> inverse(const std::vector<Complex>& ĉ);
+  void writeToA(unsigned i, Real ai);
+  std::vector<Real> convolve(const std::vector<Real>& Γh, const std::vector<Real>& R);
 };
 
+class LogarithmicFourierTransform {
+private:
+  Complex* a;
+  Complex* â;
+  fftw_plan a_to_â;
+  fftw_plan â_to_a;
+  unsigned N;
+  Real k;
+  Real Δτ;
+  Real Δω;
+  Real Δs;
+  Real τₛ;
+  Real ωₛ;
+  Real sₛ;
+public:
+  LogarithmicFourierTransform(unsigned N, Real k, Real Δω, Real Δτ, Real Δs);
+  ~LogarithmicFourierTransform();
+  Real τ(unsigned n) const;
+  Real ω(unsigned n) const;
+  Real t(unsigned n) const;
+  Real ν(unsigned n) const;
+  Real s(unsigned n) const;
+  std::vector<Complex> fourier(const std::vector<Real>& c, bool symmetric);
+  std::vector<Complex> fourier();
+  std::vector<Real> inverse(const std::vector<Complex>& ĉ);
+  void writeToA(unsigned i, Real ai);
+  std::vector<Real> convolve(const std::vector<Real>& Γh, const std::vector<Real>& R);
+};
+
+Complex gamma(Complex z);
+
 std::string fourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real y, unsigned log2n, Real τₘₐₓ);
 Real energy(const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real y, Real Δτ);
 std::tuple<std::vector<Complex>, std::vector<Complex>> RddfCtdfCt(FourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ);
+Real estimateZ(FourierTransform& fft, const std::vector<Real>& C, const std::vector<Complex>& Ct, const std::vector<Real>& R, const std::vector<Complex>& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real y);