quantumtalbot2d.cpp

Go to the documentation of this file.

/**
 * @file quantumtalbot2d.cpp
 * @brief Implements the quantumtalbot2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/quantumtalbot2d.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 "quantumtalbot2d.h"
 
#include "core/quantum/gaussianpacket2d.h"
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <limits>
#include <stdexcept>
 
namespace qpe {
namespace {
 
constexpr double kPi = 3.141592653589793238462643383279502884;
 
/**
 * @brief Performs finite.
 * @param v The v value used by the operation.
 */
bool finite(double v) noexcept
{
    return std::isfinite(v);
}
 
/**
 * @brief Performs fail.
 * @param reason The reason value used by the operation.
 * @param message The message value used by the operation.
 */
void fail(std::string *reason, const char *message) noexcept
{
    if (reason) *reason = message;
}
 
/**
 * @brief Performs wrapped distance.
 * @param x X coordinate/value.
 * @param period The period value used by the operation.
 */
double wrappedDistance(double x, double period) noexcept
{
    if (!(period > 0.0)) return 0.0;
    return x - period * std::floor(x / period + 0.5);
}
 
} // namespace
 
bool QuantumTalbot2D::validateConfig(const QuantumTalbotInitialState2DConfig &config,
                                     std::string *reason) noexcept
{
    if (reason) reason->clear();
    if (!(config.period > 0.0) || !finite(config.period)) {
        fail(reason, "Talbot period must be finite and positive.");
        return false;
    }
    if (!(config.toothSigma > 0.0) || !finite(config.toothSigma) ||
        config.toothSigma >= 0.5 * config.period) {
        fail(reason, "Talbot tooth sigma must be positive and smaller than half the period.");
        return false;
    }
    if (!(config.sigmaY > 0.0) || !finite(config.sigmaY)) {
        fail(reason, "Talbot longitudinal sigma must be finite and positive.");
        return false;
    }
    if (!finite(config.offsetX) || !finite(config.centerY) || !finite(config.momentumY)) {
        fail(reason, "Talbot initial-state position and momentum values must be finite.");
        return false;
    }
    return true;
}
 
bool QuantumTalbot2D::initialize(ComplexField2D &field,
                                 const QuantumTalbotInitialState2DConfig &config,
                                 double hbar,
                                 std::string *reason)
{
    if (reason) reason->clear();
    if (field.empty()) {
        fail(reason, "Talbot initial state requires a non-empty field.");
        return false;
    }
    if (!validateConfig(config, reason)) return false;
    if (!(hbar > 0.0) || !finite(hbar)) {
        fail(reason, "Talbot initial state requires positive finite hbar.");
        return false;
    }
 
    const double physicalWidth = static_cast<double>(field.width()) * field.dx();
    const double periods = physicalWidth / config.period;
    const double nearest = std::floor(periods + 0.5);
    if (nearest < 2.0 || std::fabs(periods - nearest) > 1.0e-5) {
        fail(reason, "Talbot comb period must divide the FFT domain width exactly.");
        return false;
    }
 
    const double periodCells = config.period / static_cast<double>(field.dx());
    const double nearestCells = std::floor(periodCells + 0.5);
    if (nearestCells < 4.0 || std::fabs(periodCells - nearestCells) > 1.0e-5 ||
        (static_cast<long long>(nearestCells) % 2LL) != 0LL) {
        fail(reason, "Talbot period must span an even integer number of grid cells.");
        return false;
    }
 
    const double inv4SigmaX2 = 1.0 / (4.0 * config.toothSigma * config.toothSigma);
    const double inv4SigmaY2 = 1.0 / (4.0 * config.sigmaY * config.sigmaY);
    const double invHbar = 1.0 / hbar;
 
    for (std::size_t y = 0u; y < field.height(); ++y) {
        const double py = GaussianPacket2D::yCoordinate(field, y);
        const double dy = py - config.centerY;
        const double envelopeY = std::exp(-(dy * dy) * inv4SigmaY2);
        const double phase = config.momentumY * py * invHbar;
        const std::complex<float> phaseY(
            static_cast<float>(std::cos(phase)),
            static_cast<float>(std::sin(phase)));
 
        for (std::size_t x = 0u; x < field.width(); ++x) {
            const double px = GaussianPacket2D::xCoordinate(field, x);
            const double d = wrappedDistance(px - config.offsetX, config.period);
            const double comb = std::exp(-(d * d) * inv4SigmaX2);
            field(x, y) = static_cast<float>(comb * envelopeY) * phaseY;
        }
    }
 
    if (!field.normalize()) {
        fail(reason, "Talbot initial-state normalization failed.");
        return false;
    }
    return true;
}
 
double QuantumTalbot2D::talbotTime(double period,
                                   double mass,
                                   double hbar) noexcept
{
    if (!(period > 0.0) || !(mass > 0.0) || !(hbar > 0.0) ||
        !finite(period) || !finite(mass) || !finite(hbar)) {
        return 0.0;
    }
    return mass * period * period / (kPi * hbar);
}
 
std::vector<float> QuantumTalbot2D::marginalProbabilityX(const ComplexField2D &field)
{
    std::vector<float> result(field.width(), 0.0f);
    if (field.empty()) return result;
 
    const double dy = static_cast<double>(field.dy());
    for (std::size_t x = 0u; x < field.width(); ++x) {
        double sum = 0.0;
        for (std::size_t y = 0u; y < field.height(); ++y) {
            const ComplexField2D::Value value = field(x, y);
            sum += static_cast<double>(value.real()) * value.real() +
                   static_cast<double>(value.imag()) * value.imag();
        }
        result[x] = static_cast<float>(sum * dy);
    }
    return result;
}
 
double QuantumTalbot2D::profileCorrelation(const std::vector<float> &a,
                                           const std::vector<float> &b,
                                           double dx) noexcept
{
    if (a.empty() || a.size() != b.size() || !(dx > 0.0) || !finite(dx)) return 0.0;
 
    double numerator = 0.0;
    double normA = 0.0;
    double normB = 0.0;
    for (std::size_t i = 0u; i < a.size(); ++i) {
        const double av = std::max(0.0, static_cast<double>(a[i]));
        const double bv = std::max(0.0, static_cast<double>(b[i]));
        numerator += av * bv;
        normA += av * av;
        normB += bv * bv;
    }
    if (!(normA > 0.0) || !(normB > 0.0)) return 0.0;
    const double value = numerator / std::sqrt(normA * normB);
    return std::max(0.0, std::min(1.0, value));
}
 
std::vector<float> QuantumTalbot2D::circularShiftHalfPeriod(const std::vector<float> &profile,
                                                            double period,
                                                            double dx)
{
    if (profile.empty() || !(period > 0.0) || !(dx > 0.0)) return std::vector<float>();
    const double shiftExact = 0.5 * period / dx;
    const long long shift = static_cast<long long>(std::floor(shiftExact + 0.5));
    if (std::fabs(shiftExact - static_cast<double>(shift)) > 1.0e-5) {
        return std::vector<float>();
    }
 
    std::vector<float> result(profile.size(), 0.0f);
    const long long n = static_cast<long long>(profile.size());
    for (long long i = 0; i < n; ++i) {
        long long source = (i - shift) % n;
        if (source < 0) source += n;
        result[static_cast<std::size_t>(i)] = profile[static_cast<std::size_t>(source)];
    }
    return result;
}
 
QuantumTalbotHistory2D::QuantumTalbotHistory2D(std::size_t maximumSamples)
    : maximumSamples_(std::max<std::size_t>(1u, maximumSamples))
{
}
 
void QuantumTalbotHistory2D::clear() noexcept
{
    gridWidth_ = 0u;
    gridHeight_ = 0u;
    dx_ = 0.0f;
    dy_ = 0.0f;
    config_ = QuantumTalbotInitialState2DConfig();
    referenceMarginal_.clear();
    halfShiftReference_.clear();
    samples_.clear();
}
 
bool QuantumTalbotHistory2D::captureReference(const ComplexField2D &field,
                                              const QuantumTalbotInitialState2DConfig &config)
{
    clear();
    std::string reason;
    if (field.empty() || !QuantumTalbot2D::validateConfig(config, &reason)) return false;
    referenceMarginal_ = QuantumTalbot2D::marginalProbabilityX(field);
    halfShiftReference_ = QuantumTalbot2D::circularShiftHalfPeriod(
        referenceMarginal_, config.period, static_cast<double>(field.dx()));
    if (referenceMarginal_.empty() || halfShiftReference_.empty()) {
        clear();
        return false;
    }
    gridWidth_ = field.width();
    gridHeight_ = field.height();
    dx_ = field.dx();
    dy_ = field.dy();
    config_ = config;
    return true;
}
 
void QuantumTalbotHistory2D::setMaximumSamples(std::size_t maximumSamples)
{
    maximumSamples_ = std::max<std::size_t>(1u, maximumSamples);
    trimToLimit();
}
 
bool QuantumTalbotHistory2D::append(double time, const ComplexField2D &field)
{
    if (!hasReference() || field.empty() || field.width() != gridWidth_ ||
        field.height() != gridHeight_ || std::fabs(field.dx() - dx_) > 1.0e-9f ||
        std::fabs(field.dy() - dy_) > 1.0e-9f || !finite(time)) {
        return false;
    }
 
    QuantumTalbotSample2D sample;
    sample.time = time;
    sample.marginalX = QuantumTalbot2D::marginalProbabilityX(field);
    sample.directCorrelation = QuantumTalbot2D::profileCorrelation(
        sample.marginalX, referenceMarginal_, static_cast<double>(dx_));
    sample.halfShiftCorrelation = QuantumTalbot2D::profileCorrelation(
        sample.marginalX, halfShiftReference_, static_cast<double>(dx_));
    samples_.push_back(sample);
    trimToLimit();
    return true;
}
 
void QuantumTalbotHistory2D::trimToLimit()
{
    if (samples_.size() <= maximumSamples_) return;
    const std::size_t removeCount = samples_.size() - maximumSamples_;
    samples_.erase(samples_.begin(), samples_.begin() + static_cast<std::ptrdiff_t>(removeCount));
}
 
} // namespace qpe