quantumbouncer2d.cpp

Go to the documentation of this file.

/**
 * @file quantumbouncer2d.cpp
 * @brief Implements the quantumbouncer2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/quantumbouncer2d.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 "quantumbouncer2d.h"
 
#include "core/quantum/gaussianpacket2d.h"
#include "core/quantum/wavepacketobservables2d.h"
 
#include <algorithm>
#include <cmath>
#include <stdexcept>
 
namespace qpe {
namespace {
 
/**
 * @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;
}
 
const double kAiryZeros[] = {
    2.338107410459767,
    4.087949444130970,
    5.520559828095551,
    6.786708090071759,
    7.944133587120853,
    9.022650853340981,
    10.04017434155809,
    11.00852430373326
};
 
} // namespace
 
bool QuantumBouncer2D::identify(const PotentialScene2D &scene,
                                QuantumBouncerGeometry2D &geometry,
                                std::string *reason) noexcept
{
    geometry = QuantumBouncerGeometry2D();
    if (reason) reason->clear();
 
    const PotentialObject2D *floor = nullptr;
    const PotentialObject2D *linear = nullptr;
    for (const PotentialObject2D &object : scene.objects()) {
        if (!object.enabled) continue;
        if (object.type == PotentialObjectType2D::Rectangle && object.value > 0.0f &&
            object.width > 4.0f * object.height) {
            if (floor != nullptr) {
                fail(reason, "Quantum bouncer expects exactly one horizontal positive Rectangle floor.");
                return false;
            }
            floor = &object;
        } else if (object.type == PotentialObjectType2D::Linear &&
                   std::fabs(object.slopeX) < 1.0e-7f && object.slopeY > 0.0f) {
            if (linear != nullptr) {
                fail(reason, "Quantum bouncer expects exactly one downward Linear-Y force term.");
                return false;
            }
            linear = &object;
        } else {
            fail(reason, "Quantum bouncer scene must contain only one floor Rectangle and one Linear-Y term.");
            return false;
        }
    }
 
    if (floor == nullptr || linear == nullptr) {
        fail(reason, "Quantum bouncer requires one floor Rectangle and one positive Linear-Y slope.");
        return false;
    }
 
    geometry.valid = true;
    geometry.floorY = static_cast<double>(floor->centerY) + 0.5 * floor->height;
    geometry.floorPotential = floor->value;
    geometry.force = linear->slopeY;
    return true;
}
 
double QuantumBouncer2D::airyLength(const QuantumBouncerGeometry2D &geometry,
                                    double mass,
                                    double hbar) noexcept
{
    if (!geometry.valid || !(geometry.force > 0.0) || !(mass > 0.0) || !(hbar > 0.0) ||
        !std::isfinite(mass) || !std::isfinite(hbar)) {
        return 0.0;
    }
    return std::cbrt((hbar * hbar) / (2.0 * mass * geometry.force));
}
 
double QuantumBouncer2D::airyEnergy(int stateIndex,
                                    const QuantumBouncerGeometry2D &geometry,
                                    double mass,
                                    double hbar) noexcept
{
    if (stateIndex < 0 || stateIndex >= static_cast<int>(sizeof(kAiryZeros) / sizeof(kAiryZeros[0]))) {
        return 0.0;
    }
    const double ell = airyLength(geometry, mass, hbar);
    return ell > 0.0 ? geometry.force * ell * kAiryZeros[stateIndex] : 0.0;
}
 
double QuantumBouncer2D::airyTurningHeight(int stateIndex,
                                           const QuantumBouncerGeometry2D &geometry,
                                           double mass,
                                           double hbar) noexcept
{
    const double ell = airyLength(geometry, mass, hbar);
    if (!(ell > 0.0) || stateIndex < 0 ||
        stateIndex >= static_cast<int>(sizeof(kAiryZeros) / sizeof(kAiryZeros[0]))) {
        return geometry.floorY;
    }
    return geometry.floorY + ell * kAiryZeros[stateIndex];
}
 
std::vector<float> QuantumBouncer2D::marginalProbabilityY(const ComplexField2D &field)
{
    std::vector<float> result(field.height(), 0.0f);
    if (field.empty()) return result;
 
    const double dx = static_cast<double>(field.dx());
    for (std::size_t y = 0u; y < field.height(); ++y) {
        double sum = 0.0;
        for (std::size_t x = 0u; x < field.width(); ++x) {
            const ComplexField2D::Value value = field(x, y);
            const double re = value.real();
            const double im = value.imag();
            sum += re * re + im * im;
        }
        result[y] = static_cast<float>(sum * dx);
    }
    return result;
}
 
double QuantumBouncer2D::probabilityNearFloor(const ComplexField2D &field,
                                               const QuantumBouncerGeometry2D &geometry,
                                               double bandHeight) noexcept
{
    if (field.empty() || !geometry.valid || !(bandHeight > 0.0)) return 0.0;
    const double cellArea = static_cast<double>(field.dx()) * field.dy();
    double probability = 0.0;
    for (std::size_t y = 0u; y < field.height(); ++y) {
        const double py = GaussianPacket2D::yCoordinate(field, y);
        if (py < geometry.floorY || py > geometry.floorY + bandHeight) continue;
        for (std::size_t x = 0u; x < field.width(); ++x) {
            const ComplexField2D::Value value = field(x, y);
            probability += (static_cast<double>(value.real()) * value.real() +
                            static_cast<double>(value.imag()) * value.imag()) * cellArea;
        }
    }
    return probability;
}
 
double QuantumBouncer2D::classicalHeight(double time,
                                         const QuantumBouncerGeometry2D &geometry,
                                         double mass,
                                         double initialY,
                                         double initialMomentumY) noexcept
{
    if (!geometry.valid || !(geometry.force > 0.0) || !(mass > 0.0) || !(time >= 0.0)) {
        return initialY;
    }
 
    const double a = geometry.force / mass;
    double y = std::max(initialY, geometry.floorY);
    double v = initialMomentumY / mass;
    double remaining = time;
 
    for (int impact = 0; impact < 128 && remaining > 1.0e-12; ++impact) {
        const double height = std::max(0.0, y - geometry.floorY);
        const double discriminant = std::max(0.0, v * v + 2.0 * a * height);
        const double hit = (v + std::sqrt(discriminant)) / a;
        if (!(hit > 1.0e-12) || hit > remaining) {
            return std::max(geometry.floorY, y + v * remaining - 0.5 * a * remaining * remaining);
        }
        y = geometry.floorY;
        v = -(v - a * hit); // elastic reflection
        remaining -= hit;
    }
    return std::max(geometry.floorY, y);
}
 
QuantumBouncerHistory2D::QuantumBouncerHistory2D(std::size_t maximumSamples)
    : maximumSamples_(maximumSamples > 0u ? maximumSamples : 1u)
{
    samples_.reserve(std::min<std::size_t>(maximumSamples_, 1200u));
}
 
void QuantumBouncerHistory2D::clear() noexcept
{
    geometry_ = QuantumBouncerGeometry2D();
    initialY_ = 0.0;
    initialMomentumY_ = 0.0;
    dy_ = 0.0f;
    gridHeight_ = 0u;
    samples_.clear();
}
 
bool QuantumBouncerHistory2D::configure(const QuantumBouncerGeometry2D &geometry,
                                        double initialY,
                                        double initialMomentumY)
{
    clear();
    if (!geometry.valid || !std::isfinite(initialY) || !std::isfinite(initialMomentumY)) return false;
    geometry_ = geometry;
    initialY_ = initialY;
    initialMomentumY_ = initialMomentumY;
    return true;
}
 
void QuantumBouncerHistory2D::setMaximumSamples(std::size_t maximumSamples)
{
    if (maximumSamples == 0u) throw std::invalid_argument("QuantumBouncerHistory2D maximumSamples must be > 0");
    maximumSamples_ = maximumSamples;
    trimToLimit();
}
 
bool QuantumBouncerHistory2D::append(double time, const ComplexField2D &field,
                                     double mass, double hbar)
{
    if (!configured() || field.empty() || !(mass > 0.0) || !(hbar > 0.0)) return false;
    QuantumBouncerSample2D sample;
    sample.time = time;
    const SpatialObservables2D spatial = WavePacketObservables2D::spatial(field);
    sample.meanY = spatial.meanY;
    sample.sigmaY = spatial.sigmaY;
    sample.norm = spatial.norm2;
    const double ell = QuantumBouncer2D::airyLength(geometry_, mass, hbar);
    sample.floorProbability = QuantumBouncer2D::probabilityNearFloor(field, geometry_, 2.0 * ell);
    sample.marginalY = QuantumBouncer2D::marginalProbabilityY(field);
    dy_ = field.dy();
    gridHeight_ = field.height();
    samples_.push_back(sample);
    trimToLimit();
    return true;
}
 
void QuantumBouncerHistory2D::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