quantumrevivals2d.cpp

Go to the documentation of this file.

/**
 * @file quantumrevivals2d.cpp
 * @brief Implements the quantumrevivals2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/quantumrevivals2d.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 "quantumrevivals2d.h"
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <stdexcept>
#include <vector>
 
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 near.
 * @param a The a value used by the operation.
 * @param b The b value used by the operation.
 * @param scale The scale value used by the operation.
 */
bool near(double a, double b, double scale = 1.0) noexcept
{
    return std::fabs(a - b) <= 1.0e-4 * std::max(1.0, scale);
}
 
} // namespace
 
bool QuantumRevivals2D::identifyBox(const PotentialScene2D &scene,
                                    QuantumBoxGeometry2D &geometry,
                                    std::string *reason) noexcept
{
    geometry = QuantumBoxGeometry2D();
    if (reason) reason->clear();
 
    std::vector<const PotentialObject2D *> vertical;
    std::vector<const PotentialObject2D *> horizontal;
 
    for (const PotentialObject2D &object : scene.objects()) {
        if (!object.enabled) continue;
        if (object.type != PotentialObjectType2D::Rectangle || !(object.value > 0.0f)) {
            fail(reason, "Quantum revival box expects four enabled positive Rectangle walls.");
            return false;
        }
        if (object.height > object.width) vertical.push_back(&object);
        else if (object.width > object.height) horizontal.push_back(&object);
        else {
            fail(reason, "Quantum revival box walls must be clearly vertical or horizontal.");
            return false;
        }
    }
 
    if (vertical.size() != 2u || horizontal.size() != 2u) {
        fail(reason, "Quantum revival box needs exactly two vertical and two horizontal walls.");
        return false;
    }
 
    std::sort(vertical.begin(), vertical.end(),
              [](const PotentialObject2D *a, const PotentialObject2D *b) {
                  return a->centerX < b->centerX;
              });
    std::sort(horizontal.begin(), horizontal.end(),
              [](const PotentialObject2D *a, const PotentialObject2D *b) {
                  return a->centerY < b->centerY;
              });
 
    const PotentialObject2D &left = *vertical[0];
    const PotentialObject2D &right = *vertical[1];
    const PotentialObject2D &bottom = *horizontal[0];
    const PotentialObject2D &top = *horizontal[1];
 
    const double leftInner = static_cast<double>(left.centerX) + 0.5 * left.width;
    const double rightInner = static_cast<double>(right.centerX) - 0.5 * right.width;
    const double bottomInner = static_cast<double>(bottom.centerY) + 0.5 * bottom.height;
    const double topInner = static_cast<double>(top.centerY) - 0.5 * top.height;
    const double width = rightInner - leftInner;
    const double height = topInner - bottomInner;
 
    if (!(width > 0.0) || !(height > 0.0)) {
        fail(reason, "Quantum revival box walls overlap or leave no finite interior.");
        return false;
    }
 
    const double wallValue = left.value;
    if (!near(right.value, wallValue, wallValue) ||
        !near(bottom.value, wallValue, wallValue) ||
        !near(top.value, wallValue, wallValue)) {
        fail(reason, "Quantum revival box walls must use the same potential height.");
        return false;
    }
    if (!near(left.width, right.width, left.width) ||
        !near(bottom.height, top.height, bottom.height)) {
        fail(reason, "Opposite quantum-box walls must use equal thickness.");
        return false;
    }
 
    const double centerX = 0.5 * (leftInner + rightInner);
    const double centerY = 0.5 * (bottomInner + topInner);
    if (!near(left.centerY, centerY, height) || !near(right.centerY, centerY, height) ||
        !near(bottom.centerX, centerX, width) || !near(top.centerX, centerX, width)) {
        fail(reason, "Quantum revival box walls are not centred on one rectangular interior.");
        return false;
    }
 
    geometry.valid = true;
    geometry.centerX = centerX;
    geometry.centerY = centerY;
    geometry.innerWidth = width;
    geometry.innerHeight = height;
    geometry.wallThicknessX = left.width;
    geometry.wallThicknessY = bottom.height;
    geometry.wallPotential = wallValue;
    return true;
}
 
double QuantumRevivals2D::revivalTime(double boxLength,
                                      double mass,
                                      double hbar) noexcept
{
    if (!(boxLength > 0.0) || !(mass > 0.0) || !(hbar > 0.0) ||
        !std::isfinite(boxLength) || !std::isfinite(mass) || !std::isfinite(hbar)) {
        return 0.0;
    }
    return 4.0 * mass * boxLength * boxLength / (kPi * hbar);
}
 
double QuantumRevivals2D::revivalTimeX(const QuantumBoxGeometry2D &geometry,
                                       double mass,
                                       double hbar) noexcept
{
    return geometry.valid ? revivalTime(geometry.innerWidth, mass, hbar) : 0.0;
}
 
double QuantumRevivals2D::revivalTimeY(const QuantumBoxGeometry2D &geometry,
                                       double mass,
                                       double hbar) noexcept
{
    return geometry.valid ? revivalTime(geometry.innerHeight, mass, hbar) : 0.0;
}
 
std::vector<float> QuantumRevivals2D::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);
            const double re = value.real();
            const double im = value.imag();
            sum += re * re + im * im;
        }
        result[x] = static_cast<float>(sum * dy);
    }
    return result;
}
 
double QuantumRevivals2D::fidelity(const ComplexField2D &reference,
                                   const ComplexField2D &field) noexcept
{
    if (reference.empty() || field.empty() ||
        reference.width() != field.width() || reference.height() != field.height() ||
        std::fabs(reference.dx() - field.dx()) > 1.0e-9f ||
        std::fabs(reference.dy() - field.dy()) > 1.0e-9f) {
        return 0.0;
    }
 
    std::complex<double> overlap(0.0, 0.0);
    double normA = 0.0;
    double normB = 0.0;
    for (std::size_t i = 0u; 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;
        normA += std::norm(a);
        normB += std::norm(b);
    }
    if (!(normA > 0.0) || !(normB > 0.0)) return 0.0;
    const double value = std::norm(overlap) / (normA * normB);
    return std::max(0.0, std::min(1.0, value));
}
 
QuantumRevivalHistory2D::QuantumRevivalHistory2D(std::size_t maximumSamples)
    : maximumSamples_(maximumSamples > 0u ? maximumSamples : 1u)
{
    samples_.reserve(std::min<std::size_t>(maximumSamples_, 1200u));
}
 
void QuantumRevivalHistory2D::clear() noexcept
{
    reference_ = ComplexField2D();
    samples_.clear();
}
 
bool QuantumRevivalHistory2D::captureReference(const ComplexField2D &field)
{
    if (field.empty()) {
        clear();
        return false;
    }
    reference_ = field;
    samples_.clear();
    return true;
}
 
void QuantumRevivalHistory2D::setMaximumSamples(std::size_t maximumSamples)
{
    if (maximumSamples == 0u) {
        throw std::invalid_argument("QuantumRevivalHistory2D maximumSamples must be > 0");
    }
    maximumSamples_ = maximumSamples;
    trimToLimit();
}
 
bool QuantumRevivalHistory2D::append(double time, const ComplexField2D &field)
{
    if (!hasReference() || field.empty() || field.width() != reference_.width() ||
        field.height() != reference_.height()) {
        return false;
    }
 
    QuantumRevivalSample2D sample;
    sample.time = time;
    sample.fidelity = QuantumRevivals2D::fidelity(reference_, field);
    sample.marginalX = QuantumRevivals2D::marginalProbabilityX(field);
    samples_.push_back(sample);
    trimToLimit();
    return true;
}
 
void QuantumRevivalHistory2D::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