quantumbackflow2d.cpp

Go to the documentation of this file.

/**
 * @file quantumbackflow2d.cpp
 * @brief Implements the quantumbackflow2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/quantumbackflow2d.cpp. This file is part of Quantum Phenomena Explorer (QPE),
 * maintained by Officina Turini. See LICENSE_IT.txt / LICENSE.txt for the
 * QPE freeware source-available terms and THIRD_PARTY_NOTICES.txt for
 * third-party components.
 */
 
#include "quantumbackflow2d.h"
 
#include "core/math/fft2d.h"
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <limits>
 
namespace qpe {
namespace {
 
constexpr double kPi = 3.141592653589793238462643383279502884;
constexpr double kTwoPi = 2.0 * kPi;
 
/**
 * @brief Performs finite.
 * @param v The v value used by the operation.
 */
bool finite(double v) noexcept { return std::isfinite(v); }
 
/**
 * @brief Performs fft wave number.
 * @param index Zero-based index of the selected item.
 * @param count The count value used by the operation.
 * @param spacing The spacing value used by the operation.
 */
double fftWaveNumber(std::size_t index, std::size_t count, double spacing) noexcept
{
    const long long n = (index <= count / 2u)
        ? static_cast<long long>(index)
        : static_cast<long long>(index) - static_cast<long long>(count);
    return kTwoPi * static_cast<double>(n) /
           (static_cast<double>(count) * spacing);
}
 
/**
 * @brief Performs coordinate.
 * @param index Zero-based index of the selected item.
 * @param count The count value used by the operation.
 * @param spacing The spacing value used by the operation.
 */
double coordinate(std::size_t index, std::size_t count, double spacing) noexcept
{
    return (static_cast<double>(index) - 0.5 * static_cast<double>(count)) * spacing;
}
 
std::size_t nearestColumn(const ComplexField2D &field, double x) noexcept
{
    const double f = x / static_cast<double>(field.dx()) +
                     0.5 * static_cast<double>(field.width());
    long long i = static_cast<long long>(std::llround(f));
    if (i < 1) i = 1;
    if (i > static_cast<long long>(field.width()) - 2)
        i = static_cast<long long>(field.width()) - 2;
    return static_cast<std::size_t>(i);
}
 
} // namespace
 
bool QuantumBackflow2D::validateConfig(const QuantumBackflow2DConfig &c,
                                        std::string *reason) noexcept
{
    const bool ok = finite(c.detectorX) && finite(c.slowWaveNumber) &&
                    finite(c.fastWaveNumber) && finite(c.momentumSigma) &&
                    finite(c.fastAmplitudeRatio) && finite(c.relativePhaseDegrees) &&
                    finite(c.initialCenterX) && finite(c.transverseSigma) &&
                    finite(c.analysisEndTime) &&
                    c.slowWaveNumber > 0.0 && c.slowWaveNumber <= 8.0 &&
                    c.fastWaveNumber > c.slowWaveNumber && c.fastWaveNumber <= 12.0 &&
                    c.momentumSigma >= 0.03 && c.momentumSigma <= 3.0 &&
                    c.fastAmplitudeRatio > 0.0 && c.fastAmplitudeRatio <= 2.0 &&
                    c.transverseSigma >= 0.2 && c.transverseSigma <= 12.0 &&
                    c.analysisEndTime >= 0.1 && c.analysisEndTime <= 100.0 &&
                    std::abs(c.detectorX) <= 1000.0 &&
                    std::abs(c.initialCenterX) <= 1000.0;
    if (!ok && reason)
        *reason = "QuantumBackflow2D: invalid configuration";
    return ok;
}
 
bool QuantumBackflow2D::prepareState(ComplexField2D &field,
                                      const QuantumBackflow2DConfig &c,
                                      double hbar,
                                      std::string *reason)
{
    if (field.empty() || !FFT2D::supportsSize(field.width(), field.height()) ||
        !finite(hbar) || hbar <= 0.0 || !validateConfig(c, reason)) {
        if (reason && reason->empty()) *reason = "QuantumBackflow2D: invalid field or hbar";
        return false;
    }
 
    ComplexField2D spectrum(field.width(), field.height(), field.dx(), field.dy());
    spectrum.fill(ComplexField2D::Value(0.0f, 0.0f));
 
    const double lx = static_cast<double>(field.width()) * static_cast<double>(field.dx());
    const double ly = static_cast<double>(field.height()) * static_cast<double>(field.dy());
    const double x0IndexCoordinates = c.initialCenterX + 0.5 * lx;
    const double y0IndexCoordinates = 0.5 * ly;
    const double phase = c.relativePhaseDegrees * kPi / 180.0;
    const std::complex<double> relative(std::cos(phase), std::sin(phase));
    const double invTwoSigma2 = 0.5 / (c.momentumSigma * c.momentumSigma);
    const double transverseKsigma = 1.0 / c.transverseSigma;
    const double invTwoKySigma2 = 0.5 / (transverseKsigma * transverseKsigma);
 
    for (std::size_t y = 0; y < field.height(); ++y) {
        const double ky = fftWaveNumber(y, field.height(), static_cast<double>(field.dy()));
        const double gy = std::exp(-ky * ky * invTwoKySigma2);
        const std::complex<double> shiftY = std::exp(std::complex<double>(0.0, -ky * y0IndexCoordinates));
        for (std::size_t x = 0; x < field.width(); ++x) {
            const double kx = fftWaveNumber(x, field.width(), static_cast<double>(field.dx()));
            if (kx <= 0.0) continue; // exact positive-momentum support by construction
 
            const double d1 = kx - c.slowWaveNumber;
            const double d2 = kx - c.fastWaveNumber;
            const double g1 = std::exp(-d1 * d1 * invTwoSigma2);
            const double g2 = std::exp(-d2 * d2 * invTwoSigma2);
            const std::complex<double> ax = g1 + c.fastAmplitudeRatio * relative * g2;
            const std::complex<double> shiftX = std::exp(std::complex<double>(0.0, -kx * x0IndexCoordinates));
            const std::complex<double> a = ax * gy * shiftX * shiftY;
            spectrum(x, y) = ComplexField2D::Value(static_cast<float>(a.real()), static_cast<float>(a.imag()));
        }
    }
 
    FFT2D::inverse(spectrum);
    field = spectrum;
    if (!field.normalize()) {
        if (reason) *reason = "QuantumBackflow2D: state normalization failed";
        return false;
    }
    return true;
}
 
double QuantumBackflow2D::positiveMomentumFraction(const ComplexField2D &field) noexcept
{
    if (field.empty() || !FFT2D::supportsSize(field.width(), field.height())) return 0.0;
    try {
        ComplexField2D spectrum = field;
        FFT2D::forward(spectrum);
        double positive = 0.0;
        double total = 0.0;
        for (std::size_t y = 0; y < field.height(); ++y) {
            for (std::size_t x = 0; x < field.width(); ++x) {
                const double w = std::norm(spectrum(x, y));
                total += w;
                if (fftWaveNumber(x, field.width(), static_cast<double>(field.dx())) > 0.0)
                    positive += w;
            }
        }
        return total > 0.0 ? positive / total : 0.0;
    } catch (...) {
        return 0.0;
    }
}
 
QuantumBackflowSample2D QuantumBackflow2D::measure(const ComplexField2D &field,
                                                     const QuantumBackflow2DConfig &c,
                                                     double mass,
                                                     double hbar,
                                                     double time) noexcept
{
    QuantumBackflowSample2D out;
    out.time = time;
    if (field.empty() || !finite(mass) || mass <= 0.0 || !finite(hbar) || hbar <= 0.0)
        return out;
 
    const std::size_t ix = nearestColumn(field, c.detectorX);
    const double dx = static_cast<double>(field.dx());
    const double dy = static_cast<double>(field.dy());
    const double currentScale = hbar / mass;
 
    double left = 0.0;
    double right = 0.0;
    double current = 0.0;
    for (std::size_t y = 0; y < field.height(); ++y) {
        for (std::size_t x = 0; x < field.width(); ++x) {
            const double p = std::norm(field(x, y)) * dx * dy;
            if (coordinate(x, field.width(), dx) < c.detectorX) left += p;
            else right += p;
        }
        const std::complex<double> psi(field(ix, y).real(), field(ix, y).imag());
        const std::complex<double> plus(field(ix + 1u, y).real(), field(ix + 1u, y).imag());
        const std::complex<double> minus(field(ix - 1u, y).real(), field(ix - 1u, y).imag());
        const std::complex<double> dpsi = (plus - minus) / (2.0 * dx);
        current += currentScale * std::imag(std::conj(psi) * dpsi) * dy;
    }
    out.detectorCurrent = current;
    out.probabilityLeft = left;
    out.probabilityRight = right;
    return out;
}
 
void QuantumBackflowHistory2D::clear()
{
    configured_ = false;
    config_ = QuantumBackflow2DConfig();
    mass_ = 1.0;
    hbar_ = 1.0;
    positiveMomentumFraction_ = 0.0;
    minimumCurrent_ = 0.0;
    maximumBackwardFlux_ = 0.0;
    samples_.clear();
}
 
bool QuantumBackflowHistory2D::configure(const ComplexField2D &field,
                                          const QuantumBackflow2DConfig &config,
                                          double mass,
                                          double hbar)
{
    clear();
    if (field.empty() || !QuantumBackflow2D::validateConfig(config, nullptr) ||
        !finite(mass) || mass <= 0.0 || !finite(hbar) || hbar <= 0.0)
        return false;
    config_ = config;
    mass_ = mass;
    hbar_ = hbar;
    positiveMomentumFraction_ = QuantumBackflow2D::positiveMomentumFraction(field);
    configured_ = true;
    return true;
}
 
void QuantumBackflowHistory2D::append(double time, const ComplexField2D &field)
{
    if (!configured_ || field.empty()) return;
    QuantumBackflowSample2D s = QuantumBackflow2D::measure(field, config_, mass_, hbar_, time);
    if (!samples_.empty()) {
        const QuantumBackflowSample2D &p = samples_.back();
        const double dt = std::max(0.0, s.time - p.time);
        const double backward0 = std::max(0.0, -p.detectorCurrent);
        const double backward1 = std::max(0.0, -s.detectorCurrent);
        s.cumulativeBackwardFlux = p.cumulativeBackwardFlux + 0.5 * dt * (backward0 + backward1);
    }
    minimumCurrent_ = samples_.empty() ? s.detectorCurrent : std::min(minimumCurrent_, s.detectorCurrent);
    maximumBackwardFlux_ = std::max(maximumBackwardFlux_, s.cumulativeBackwardFlux);
    samples_.push_back(s);
    trim();
}
 
void QuantumBackflowHistory2D::trim()
{
    if (maximumSamples_ == 0u) { samples_.clear(); return; }
    if (samples_.size() > maximumSamples_)
        samples_.erase(samples_.begin(), samples_.begin() + static_cast<std::ptrdiff_t>(samples_.size() - maximumSamples_));
}
 
} // namespace qpe