#include <vector>
#include <complex>
#include <cmath>
#include <algorithm>
#include <cassert>
#include <fftw3.h>

// Function declaration
std::vector<std::vector<double>> ModS2T(const std::vector<std::vector<double>>& data, double dt, double du, double ERMv, double tau0, const std::vector<std::vector<double>>& Filter, int R_x = 1) {
    // Check R_x argument
    R_x = std::max(1, R_x);
    assert(R_x > 0 && "R_x must be positive");

    std::vector<std::vector<double>> processedData = data;

    // Applies linear interpolation (if necessary)
    if (R_x > 1) {
        // Implementation of linear interpolation
        // Note: This part needs a more detailed implementation
    }

    // Zero-padding
    int nt0 = processedData.size();
    int nu0 = processedData[0].size();
    int nt = std::max(static_cast<int>(std::pow(2, std::ceil(std::log2(nt0)) + 1)), static_cast<int>(Filter.size()));
    int nu = std::max(2 * nu0, static_cast<int>(Filter[0].size()));

    // Data and image grids
    std::vector<std::vector<double>> f(nt, std::vector<double>(nu));
    std::vector<std::vector<double>> ku(nt, std::vector<double>(nu));
    std::vector<std::vector<double>> fkz(nt, std::vector<double>(nu));

    for (int i = 0; i < nt; ++i) {
        for (int j = 0; j < nu; ++j) {
            f[i][j] = (i - nt/2) / (dt * nt);
            ku[i][j] = (j - nu/2) / (du * nu);
            fkz[i][j] = ERMv * std::copysign(1.0, f[i][j]) * std::sqrt(std::pow(ku[i][j], 2) + std::pow(f[i][j] / ERMv, 2));
        }
    }

    // Converting data to frequency domain
    std::vector<std::vector<std::complex<double>>> ftdata(nt, std::vector<std::complex<double>>(nu));
    // Note: FFTW library should be used here for efficient FFT computation

    if (tau0 != 0) {
        for (int i = 0; i < nt; ++i) {
            for (int j = 0; j < nu; ++j) {
                ftdata[i][j] *= std::exp(std::complex<double>(0, -2 * M_PI * (fkz[i][j] - f[i][j]) * tau0));
            }
        }
    }

    // Applies filter (if it exists)
    if (!Filter.empty()) {
        // Implementation of filter application
        // Note: This part needs a more detailed implementation
    }

    // Generating image spectrum (by linear interpolation)
    std::vector<std::vector<std::complex<double>>> ftimage(nt, std::vector<std::complex<double>>(nu));
    // Note: Linear interpolation should be implemented here

    // Converting image to space domain
    std::vector<std::vector<double>> image(nt0, std::vector<double>(nu0));
    // Note: Inverse FFT should be applied here using FFTW library

    for (int i = 0; i < nt0; ++i) {
        for (int j = 0; j < nu0; ++j) {
            image[i][j] *= R_x;
        }
    }

    // Evaluates the gain
    double enImage = 0, enData = 0;
    for (const auto& row : image) {
        for (double val : row) {
            enImage += val * val;
        }
    }
    for (const auto& row : data) {
        for (double val : row) {
            enData += val * val;
        }
    }
    double gain = std::sqrt(enImage / enData);

    return image;
}