quantumloschmidt2d.cpp

Go to the documentation of this file.

/**
 * @file quantumloschmidt2d.cpp
 * @brief Implements the quantumloschmidt2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/quantumloschmidt2d.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 "quantumloschmidt2d.h"
 
#include "core/quantum/gaussianpacket2d.h"
 
#include <algorithm>
#include <cmath>
#include <complex>
 
namespace qpe {
namespace {
 
/**
 * @brief Performs finite.
 * @param value Value to assign or evaluate.
 */
bool finite(double value) noexcept
{
    return std::isfinite(value);
}
 
/**
 * @brief Performs fail.
 * @param reason The reason value used by the operation.
 * @param message The message value used by the operation.
 */
void fail(const char **reason, const char *message) noexcept
{
    if (reason) *reason = message;
}
 
/**
 * @brief Performs same grid.
 * @param a The a value used by the operation.
 * @param b The b value used by the operation.
 */
bool sameGrid(const ComplexField2D &a, const ComplexField2D &b) noexcept
{
    return !a.empty() && !b.empty() &&
           a.width() == b.width() && a.height() == b.height() &&
           a.dx() == b.dx() && a.dy() == b.dy();
}
 
} // namespace
 
bool QuantumLoschmidtEcho2D::validateConfig(const QuantumLoschmidt2DConfig &config,
                                             const char **reason) noexcept
{
    if (!config.enabled) return true;
    if (!(config.reversalTime > 0.0) || !finite(config.reversalTime)) {
        fail(reason, "Loschmidt reversal time must be finite and positive.");
        return false;
    }
    if (!finite(config.momentumErrorX) || !finite(config.momentumErrorY)) {
        fail(reason, "Loschmidt momentum-error values must be finite.");
        return false;
    }
    if (std::fabs(config.momentumErrorX) > 1000.0 ||
        std::fabs(config.momentumErrorY) > 1000.0) {
        fail(reason, "Loschmidt momentum-error values are unreasonably large.");
        return false;
    }
    return true;
}
 
bool QuantumLoschmidtEcho2D::applyTimeReversal(ComplexField2D &field,
                                                const QuantumLoschmidt2DConfig &config,
                                                double hbar) noexcept
{
    if (field.empty() || !(hbar > 0.0) || !finite(hbar) ||
        !validateConfig(config, nullptr)) {
        return false;
    }
 
    const double invHbar = 1.0 / hbar;
    for (std::size_t y = 0; y < field.height(); ++y) {
        const double py = GaussianPacket2D::yCoordinate(field, y);
        for (std::size_t x = 0; x < field.width(); ++x) {
            const double px = GaussianPacket2D::xCoordinate(field, x);
            const std::complex<float> reversed = std::conj(field(x, y));
            const double phase =
                (config.momentumErrorX * px + config.momentumErrorY * py) * invHbar;
            const std::complex<float> kick(
                static_cast<float>(std::cos(phase)),
                static_cast<float>(std::sin(phase)));
            field(x, y) = reversed * kick;
        }
    }
    return field.normalize();
}
 
double QuantumLoschmidtEcho2D::stateEchoFidelity(const ComplexField2D &reference,
                                                  const ComplexField2D &field) noexcept
{
    if (!sameGrid(reference, field)) return 0.0;
 
    // Target = conj(reference).  Therefore conj(target) = reference in the
    // inner product <target|field>.
    std::complex<double> overlap(0.0, 0.0);
    double normReference = 0.0;
    double normField = 0.0;
    const double cell = static_cast<double>(field.dx()) * static_cast<double>(field.dy());
 
    for (std::size_t i = 0; i < field.size(); ++i) {
        const std::complex<double> a(reference.data()[i].real(), reference.data()[i].imag());
        const std::complex<double> b(field.data()[i].real(), field.data()[i].imag());
        overlap += a * b;
        normReference += std::norm(a);
        normField += std::norm(b);
    }
    normReference *= cell;
    normField *= cell;
    overlap *= cell;
    const double denom = normReference * normField;
    if (!(denom > 0.0) || !finite(denom)) return 0.0;
    const double f = std::norm(overlap) / denom;
    return std::max(0.0, std::min(1.0, f));
}
 
double QuantumLoschmidtEcho2D::densityEchoFidelity(const ComplexField2D &reference,
                                                    const ComplexField2D &field) noexcept
{
    if (!sameGrid(reference, field)) return 0.0;
 
    double bc = 0.0;
    double normReference = 0.0;
    double normField = 0.0;
    const double cell = static_cast<double>(field.dx()) * static_cast<double>(field.dy());
    for (std::size_t i = 0; i < field.size(); ++i) {
        const double pa = static_cast<double>(std::norm(reference.data()[i]));
        const double pb = static_cast<double>(std::norm(field.data()[i]));
        bc += std::sqrt(std::max(0.0, pa * pb));
        normReference += pa;
        normField += pb;
    }
    bc *= cell;
    normReference *= cell;
    normField *= cell;
    const double denom = normReference * normField;
    if (!(denom > 0.0) || !finite(denom)) return 0.0;
    const double f = (bc * bc) / denom;
    return std::max(0.0, std::min(1.0, f));
}
 
bool QuantumLoschmidtEcho2D::centroid(const ComplexField2D &field,
                                       double &meanX,
                                       double &meanY) noexcept
{
    meanX = 0.0;
    meanY = 0.0;
    if (field.empty()) return false;
 
    double weight = 0.0;
    for (std::size_t y = 0; y < field.height(); ++y) {
        const double py = GaussianPacket2D::yCoordinate(field, y);
        for (std::size_t x = 0; x < field.width(); ++x) {
            const double px = GaussianPacket2D::xCoordinate(field, x);
            const double p = static_cast<double>(std::norm(field(x, y)));
            weight += p;
            meanX += p * px;
            meanY += p * py;
        }
    }
    if (!(weight > 0.0) || !finite(weight)) {
        meanX = 0.0;
        meanY = 0.0;
        return false;
    }
    meanX /= weight;
    meanY /= weight;
    return finite(meanX) && finite(meanY);
}
 
QuantumLoschmidtHistory2D::QuantumLoschmidtHistory2D(std::size_t maximumSamples)
    : maximumSamples_(std::max<std::size_t>(2u, maximumSamples))
{
}
 
void QuantumLoschmidtHistory2D::clear() noexcept
{
    configured_ = false;
    reference_ = ComplexField2D();
    config_ = QuantumLoschmidt2DConfig();
    samples_.clear();
}
 
bool QuantumLoschmidtHistory2D::configure(const ComplexField2D &reference,
                                           const QuantumLoschmidt2DConfig &config)
{
    clear();
    if (reference.empty() || !QuantumLoschmidtEcho2D::validateConfig(config, nullptr)) {
        return false;
    }
    reference_ = reference;
    config_ = config;
    configured_ = config.enabled;
    return configured_;
}
 
bool QuantumLoschmidtHistory2D::append(double time,
                                        const ComplexField2D &field,
                                        bool reversalApplied)
{
    if (!configured_ || !sameGrid(reference_, field) || !finite(time)) return false;
 
    QuantumLoschmidtSample2D sample;
    sample.time = time;
    sample.densityEcho = QuantumLoschmidtEcho2D::densityEchoFidelity(reference_, field);
    sample.stateEcho = QuantumLoschmidtEcho2D::stateEchoFidelity(reference_, field);
    sample.norm = field.norm2();
    sample.reversalApplied = reversalApplied;
    QuantumLoschmidtEcho2D::centroid(field, sample.meanX, sample.meanY);
 
    // Keep the plot speed-independent: store on a simulation-time cadence
    // derived from the nominal 0..2T echo window, plus the reversal event.
    const double totalWindow = std::max(1.0e-9, 2.0 * config_.reversalTime);
    const double storeInterval = totalWindow /
        static_cast<double>(std::max<std::size_t>(1u, maximumSamples_ - 1u));
    const bool first = samples_.empty();
    const bool reversalEdge = reversalApplied &&
        (samples_.empty() || !samples_.back().reversalApplied);
    const bool intervalElapsed = samples_.empty() ||
        time - samples_.back().time >= std::max(1.0e-9, storeInterval) - 1.0e-12;
    if (first || reversalEdge || intervalElapsed) {
        samples_.push_back(sample);
        trimToLimit();
    }
    return true;
}
 
void QuantumLoschmidtHistory2D::trimToLimit()
{
    if (samples_.size() <= maximumSamples_) return;
    const std::size_t excess = samples_.size() - maximumSamples_;
    samples_.erase(samples_.begin(), samples_.begin() + static_cast<std::ptrdiff_t>(excess));
}
 
} // namespace qpe