quantumscar2d.cpp

Go to the documentation of this file.

/**
 * @file quantumscar2d.cpp
 * @brief Implements the quantumscar2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/quantumscar2d.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 "quantumscar2d.h"
 
#include "core/quantum/gaussianpacket2d.h"
 
#include <algorithm>
#include <cmath>
#include <complex>
 
namespace qpe {
namespace {
 
constexpr double kPi = 3.141592653589793238462643383279502884;
 
/**
 * @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 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() &&
           std::fabs(static_cast<double>(a.dx()) - static_cast<double>(b.dx())) <= 1.0e-9 &&
           std::fabs(static_cast<double>(a.dy()) - static_cast<double>(b.dy())) <= 1.0e-9;
}
 
/**
 * @brief Performs enhancement from weights.
 * @param bandWeight The band weight value used by the operation.
 * @param totalWeight The total weight value used by the operation.
 * @param geometry The geometry value used by the operation.
 * @param halfWidth The half width value used by the operation.
 */
double enhancementFromWeights(double bandWeight,
                              double totalWeight,
                              const QuantumScarGeometry2D &geometry,
                              double halfWidth) noexcept
{
    if (!(totalWeight > 0.0)) return 0.0;
    const double area = QuantumScars2D::accessibleArea(geometry);
    const double bandArea = QuantumScars2D::orbitBandArea(geometry, halfWidth);
    if (!(area > 0.0) || !(bandArea > 0.0)) return 0.0;
    const double geometricFraction = std::min(1.0, bandArea / area);
    if (!(geometricFraction > 0.0)) return 0.0;
    return (bandWeight / totalWeight) / geometricFraction;
}
 
} // namespace
 
bool QuantumScars2D::validateConfig(const QuantumScarConfig2D &config,
                                     std::string *reason) noexcept
{
    if (reason) reason->clear();
    if (!(config.orbitBandHalfWidth > 0.0) || !std::isfinite(config.orbitBandHalfWidth) ||
        config.orbitBandHalfWidth > 20.0) {
        fail(reason, "Scar orbit-band half-width must be finite and positive.");
        return false;
    }
    if (!(config.referenceEnergy >= 0.0) || !std::isfinite(config.referenceEnergy) ||
        config.referenceEnergy > 10000.0) {
        fail(reason, "Scar reference energy must be finite and non-negative.");
        return false;
    }
    return true;
}
 
bool QuantumScars2D::identify(const PotentialScene2D &scene,
                              QuantumScarGeometry2D &geometry,
                              std::string *reason) noexcept
{
    geometry = QuantumScarGeometry2D();
    if (reason) reason->clear();
 
    const PotentialObject2D *stadium = nullptr;
    for (const PotentialObject2D &object : scene.objects()) {
        if (!object.enabled) continue;
        if (object.type != PotentialObjectType2D::StadiumExterior) {
            fail(reason, "Quantum Scar analysis expects one Stadium Exterior potential object.");
            return false;
        }
        if (stadium != nullptr) {
            fail(reason, "Quantum Scar analysis expects exactly one Stadium Exterior object.");
            return false;
        }
        stadium = &object;
    }
 
    if (!stadium || stadium->value <= 0.0f || stadium->width <= 0.0f || stadium->radius <= 0.0f) {
        fail(reason, "Stadium Exterior must have positive wall height, straight length and radius.");
        return false;
    }
 
    geometry.valid = true;
    geometry.centerX = stadium->centerX;
    geometry.centerY = stadium->centerY;
    geometry.straightLength = stadium->width;
    geometry.radius = stadium->radius;
    geometry.wallPotential = stadium->value;
    return true;
}
 
bool QuantumScars2D::inside(const QuantumScarGeometry2D &geometry,
                            double x,
                            double y) noexcept
{
    if (!geometry.valid || !(geometry.straightLength > 0.0) || !(geometry.radius > 0.0)) {
        return false;
    }
    const double lx = x - geometry.centerX;
    const double ly = y - geometry.centerY;
    const double halfStraight = 0.5 * geometry.straightLength;
    if (std::fabs(lx) <= halfStraight) return std::fabs(ly) <= geometry.radius;
    const double capX = lx > 0.0 ? halfStraight : -halfStraight;
    const double dx = lx - capX;
    return dx * dx + ly * ly <= geometry.radius * geometry.radius;
}
 
double QuantumScars2D::accessibleArea(const QuantumScarGeometry2D &geometry) noexcept
{
    if (!geometry.valid || !(geometry.straightLength > 0.0) || !(geometry.radius > 0.0)) return 0.0;
    return 2.0 * geometry.radius * geometry.straightLength +
           kPi * geometry.radius * geometry.radius;
}
 
double QuantumScars2D::orbitBandArea(const QuantumScarGeometry2D &geometry,
                                     double halfWidth) noexcept
{
    if (!geometry.valid || !(halfWidth > 0.0) || !(geometry.radius > 0.0)) return 0.0;
    const double w = std::min(geometry.radius, halfWidth);
    const double r = geometry.radius;
    const double central = 2.0 * w * geometry.straightLength;
    const double root = std::sqrt(std::max(0.0, r * r - w * w));
    // The two semicircular caps together make one complete circle. This is the
    // area of the horizontal |y|<=w band through that circle.
    const double caps = 2.0 * (w * root + r * r * std::asin(std::min(1.0, w / r)));
    return std::min(accessibleArea(geometry), central + caps);
}
 
double QuantumScars2D::axialOrbitRoundTripLength(const QuantumScarGeometry2D &geometry) noexcept
{
    if (!geometry.valid) return 0.0;
    const double oneWay = geometry.straightLength + 2.0 * geometry.radius;
    return 2.0 * oneWay;
}
 
double QuantumScars2D::axialOrbitPeriod(const QuantumScarGeometry2D &geometry,
                                        double speed) noexcept
{
    if (!(speed > 0.0)) return 0.0;
    return axialOrbitRoundTripLength(geometry) / speed;
}
 
double QuantumScars2D::fidelity(const ComplexField2D &reference,
                                const ComplexField2D &field) noexcept
{
    if (!sameGrid(reference, field)) return 0.0;
    const double cellArea = static_cast<double>(field.dx()) * static_cast<double>(field.dy());
    std::complex<double> overlap(0.0, 0.0);
    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 += std::conj(a) * b * cellArea;
    }
    return std::norm(overlap);
}
 
double QuantumScars2D::orbitEnhancement(const ComplexField2D &field,
                                        const QuantumScarGeometry2D &geometry,
                                        double halfWidth) noexcept
{
    if (field.empty() || !geometry.valid || !(halfWidth > 0.0)) return 0.0;
    double total = 0.0;
    double band = 0.0;
    const double cellArea = static_cast<double>(field.dx()) * static_cast<double>(field.dy());
    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);
            if (!inside(geometry, px, py)) continue;
            const double rho = static_cast<double>(std::norm(field(x, y))) * cellArea;
            total += rho;
            if (std::fabs(py - geometry.centerY) <= halfWidth) band += rho;
        }
    }
    return enhancementFromWeights(band, total, geometry, halfWidth);
}
 
QuantumScarHistory2D::QuantumScarHistory2D(std::size_t maximumSamples)
    : maximumSamples_(std::max<std::size_t>(1u, maximumSamples))
{
}
 
void QuantumScarHistory2D::clear() noexcept
{
    configured_ = false;
    reference_ = ComplexField2D();
    geometry_ = QuantumScarGeometry2D();
    config_ = QuantumScarConfig2D();
    referenceEnergy_ = 0.0;
    hbar_ = 1.0;
    samples_.clear();
    meanDensitySum_.clear();
    spectralAccumulator_.clear();
    meanDensity_.clear();
    filteredDensity_.clear();
    accumulatedFrames_ = 0u;
    lastAccumulationTime_ = 0.0;
    totalTimeWeight_ = 0.0;
    hasAccumulationTime_ = false;
}
 
bool QuantumScarHistory2D::configure(const ComplexField2D &reference,
                                     const QuantumScarGeometry2D &geometry,
                                     const QuantumScarConfig2D &config,
                                     double referenceEnergy,
                                     double hbar)
{
    clear();
    if (reference.empty() || !geometry.valid || !(config.orbitBandHalfWidth > 0.0) ||
        !(referenceEnergy > 0.0) || !(hbar > 0.0) ||
        !std::isfinite(referenceEnergy) || !std::isfinite(hbar)) {
        return false;
    }
    reference_ = reference;
    geometry_ = geometry;
    config_ = config;
    referenceEnergy_ = referenceEnergy;
    hbar_ = hbar;
    meanDensitySum_.assign(reference.size(), 0.0);
    spectralAccumulator_.assign(reference.size(), std::complex<double>(0.0, 0.0));
    meanDensity_.assign(reference.size(), 0.0f);
    filteredDensity_.assign(reference.size(), 0.0f);
    configured_ = true;
    return true;
}
 
void QuantumScarHistory2D::setMaximumSamples(std::size_t maximumSamples)
{
    maximumSamples_ = std::max<std::size_t>(1u, maximumSamples);
    trimToLimit();
}
 
bool QuantumScarHistory2D::append(double time, const ComplexField2D &field)
{
    if (!configured_ || !sameGrid(reference_, field) || !std::isfinite(time)) return false;
 
    // Use elapsed simulation time as the quadrature weight. This keeps the
    // finite-time projection invariant when the GUI speed changes the number
    // of solver steps performed between display refreshes. The first frame gets
    // a tiny seed weight so the map is immediately defined; it becomes
    // negligible as soon as physical elapsed time accumulates.
    double weight = hasAccumulationTime_ ? 0.0 : 1.0e-9;
    if (hasAccumulationTime_) {
        const double elapsed = time - lastAccumulationTime_;
        if (elapsed > 0.0 && std::isfinite(elapsed)) weight = elapsed;
    }
    lastAccumulationTime_ = time;
    hasAccumulationTime_ = true;
 
    const double phaseAngle = referenceEnergy_ * time / hbar_;
    const std::complex<double> phase(std::cos(phaseAngle), std::sin(phaseAngle));
    const std::size_t n = field.size();
    for (std::size_t i = 0; i < n; ++i) {
        const std::complex<double> psi(field.data()[i].real(), field.data()[i].imag());
        const double rho = std::norm(psi);
        meanDensitySum_[i] += rho * weight;
        spectralAccumulator_[i] += psi * phase * weight;
    }
    totalTimeWeight_ += weight;
    ++accumulatedFrames_;
 
    const double invWeight = totalTimeWeight_ > 0.0 ? 1.0 / totalTimeWeight_ : 0.0;
    double filteredMax = 0.0;
    for (std::size_t i = 0; i < n; ++i) {
        meanDensity_[i] = static_cast<float>(meanDensitySum_[i] * invWeight);
        const double v = std::norm(spectralAccumulator_[i]);
        filteredMax = std::max(filteredMax, v);
    }
    const double invFilteredMax = filteredMax > 0.0 ? 1.0 / filteredMax : 0.0;
    for (std::size_t i = 0; i < n; ++i) {
        filteredDensity_[i] = static_cast<float>(std::norm(spectralAccumulator_[i]) * invFilteredMax);
    }
 
    QuantumScarSample2D sample;
    sample.time = time;
    sample.fidelity = QuantumScars2D::fidelity(reference_, field);
    sample.instantaneousOrbitEnhancement =
        QuantumScars2D::orbitEnhancement(field, geometry_, config_.orbitBandHalfWidth);
    sample.filteredOrbitEnhancement = filteredEnhancement();
    samples_.push_back(sample);
    trimToLimit();
    return true;
}
 
double QuantumScarHistory2D::filteredEnhancement() const noexcept
{
    if (!configured_ || spectralAccumulator_.empty()) return 0.0;
    double total = 0.0;
    double band = 0.0;
    for (std::size_t y = 0; y < reference_.height(); ++y) {
        const double py = GaussianPacket2D::yCoordinate(reference_, y);
        for (std::size_t x = 0; x < reference_.width(); ++x) {
            const double px = GaussianPacket2D::xCoordinate(reference_, x);
            if (!QuantumScars2D::inside(geometry_, px, py)) continue;
            const std::size_t i = y * reference_.width() + x;
            const double w = std::norm(spectralAccumulator_[i]);
            total += w;
            if (std::fabs(py - geometry_.centerY) <= config_.orbitBandHalfWidth) band += w;
        }
    }
    return enhancementFromWeights(band, total, geometry_, config_.orbitBandHalfWidth);
}
 
void QuantumScarHistory2D::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