quantumvortex2d.cpp

Go to the documentation of this file.

/**
 * @file quantumvortex2d.cpp
 * @brief Implements the quantumvortex2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/quantumvortex2d.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 "quantumvortex2d.h"
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <limits>
#include <stdexcept>
 
namespace qpe {
namespace {
 
const double kPi = 3.14159265358979323846;
 
/**
 * @brief Performs finite.
 * @param v The v value used by the operation.
 */
bool finite(double v) noexcept
{
    return std::isfinite(v);
}
 
/**
 * @brief Sets reason.
 * @param reason The reason value used by the operation.
 * @param text The text value used by the operation.
 */
void setReason(std::string *reason, const char *text)
{
    if (reason) *reason = text;
}
 
/**
 * @brief Performs coordinate.
 * @param i The i value used by the operation.
 * @param n The n value used by the operation.
 * @param spacing The spacing value used by the operation.
 */
double coordinate(std::size_t i, std::size_t n, double spacing) noexcept
{
    if (n == 0u) return 0.0;
    return (static_cast<double>(i) - 0.5 * static_cast<double>(n - 1u)) * spacing;
}
 
/**
 * @brief Performs wrap phase.
 * @param phase The phase value used by the operation.
 */
double wrapPhase(double phase) noexcept
{
    const double twoPi = 2.0 * kPi;
    while (phase <= -kPi) phase += twoPi;
    while (phase > kPi) phase -= twoPi;
    return phase;
}
 
} // namespace
 
bool QuantumVortex2D::validateConfig(const QuantumVortexInitialState2DConfig &config,
                                     std::string *reason)
{
    if (!finite(config.envelopeCenterX) || !finite(config.envelopeCenterY) ||
        !finite(config.envelopeSigma) || !finite(config.coreX) || !finite(config.coreY)) {
        setReason(reason, "QuantumVortex2D: configuration values must be finite");
        return false;
    }
    if (!(config.envelopeSigma > 0.0) || config.envelopeSigma > 1000.0) {
        setReason(reason, "QuantumVortex2D: envelope sigma must be positive");
        return false;
    }
    if (config.charge != 1 && config.charge != -1) {
        setReason(reason, "QuantumVortex2D: charge must be +1 or -1");
        return false;
    }
    return true;
}
 
bool QuantumVortex2D::initialize(ComplexField2D &field,
                                 const QuantumVortexInitialState2DConfig &config,
                                 std::string *reason)
{
    if (field.empty()) {
        setReason(reason, "QuantumVortex2D: field must not be empty");
        return false;
    }
    if (!validateConfig(config, reason)) return false;
 
    const double sigma = config.envelopeSigma;
    const double inv4Sigma2 = 1.0 / (4.0 * sigma * sigma);
    const double dxGrid = static_cast<double>(field.dx());
    const double dyGrid = static_cast<double>(field.dy());
 
    for (std::size_t y = 0; y < field.height(); ++y) {
        const double py = coordinate(y, field.height(), dyGrid);
        const double ey = py - config.envelopeCenterY;
        const double vy = py - config.coreY;
        for (std::size_t x = 0; x < field.width(); ++x) {
            const double px = coordinate(x, field.width(), dxGrid);
            const double ex = px - config.envelopeCenterX;
            const double vx = px - config.coreX;
            const double envelope = std::exp(-(ex * ex + ey * ey) * inv4Sigma2);
            field(x, y) = ComplexField2D::Value(
                static_cast<float>(vx * envelope),
                static_cast<float>(static_cast<double>(config.charge) * vy * envelope));
        }
    }
 
    if (!field.normalize()) {
        setReason(reason, "QuantumVortex2D: normalization failed");
        return false;
    }
    return true;
}
 
bool QuantumVortex2D::identifyTrap(const PotentialScene2D &scene,
                                   QuantumVortexTrap2D &trap,
                                   std::string *reason)
{
    trap = QuantumVortexTrap2D();
    const PotentialObject2D *harmonic = nullptr;
    for (const PotentialObject2D &object : scene.objects()) {
        if (!object.enabled) continue;
        if (object.type == PotentialObjectType2D::Harmonic) {
            if (harmonic != nullptr) {
                setReason(reason, "QuantumVortex2D: expected exactly one enabled harmonic trap");
                return false;
            }
            harmonic = &object;
        } else {
            setReason(reason, "QuantumVortex2D: vortex guided analysis expects only one harmonic object");
            return false;
        }
    }
    if (!harmonic) {
        setReason(reason, "QuantumVortex2D: harmonic trap not found");
        return false;
    }
    const double kx = static_cast<double>(harmonic->curvatureX);
    const double ky = static_cast<double>(harmonic->curvatureY);
    if (!(kx > 0.0) || !(ky > 0.0) || !finite(kx) || !finite(ky)) {
        setReason(reason, "QuantumVortex2D: trap curvature must be positive");
        return false;
    }
    const double scale = std::max(1.0, std::max(std::abs(kx), std::abs(ky)));
    if (std::abs(kx - ky) > 1.0e-6 * scale) {
        setReason(reason, "QuantumVortex2D: guided vortex requires an isotropic harmonic trap");
        return false;
    }
 
    trap.valid = true;
    trap.centerX = static_cast<double>(harmonic->centerX);
    trap.centerY = static_cast<double>(harmonic->centerY);
    trap.curvature = 0.5 * (kx + ky);
    return true;
}
 
double QuantumVortex2D::trapAngularFrequency(const QuantumVortexTrap2D &trap,
                                             double mass)
{
    if (!trap.valid || !(trap.curvature > 0.0) || !(mass > 0.0) || !finite(mass)) return 0.0;
    return std::sqrt(trap.curvature / mass);
}
 
double QuantumVortex2D::trapPeriod(const QuantumVortexTrap2D &trap,
                                   double mass)
{
    const double omega = trapAngularFrequency(trap, mass);
    return omega > 0.0 ? 2.0 * kPi / omega : 0.0;
}
 
QuantumVortexCore2D QuantumVortex2D::detectCore(
    const ComplexField2D &field,
    const QuantumVortexInitialState2DConfig &config,
    const QuantumVortexTrap2D &trap)
{
    QuantumVortexCore2D best;
    if (field.width() < 3u || field.height() < 3u || !trap.valid) return best;
 
    const double maxRadius = std::max(2.5, 4.0 * config.envelopeSigma);
    const double maxRadius2 = maxRadius * maxRadius;
    double bestDensity = std::numeric_limits<double>::infinity();
    const int expectedSign = config.charge >= 0 ? 1 : -1;
 
    for (std::size_t y = 0; y + 1u < field.height(); ++y) {
        for (std::size_t x = 0; x + 1u < field.width(); ++x) {
            const double cx = 0.5 * (coordinate(x, field.width(), field.dx()) +
                                     coordinate(x + 1u, field.width(), field.dx()));
            const double cy = 0.5 * (coordinate(y, field.height(), field.dy()) +
                                     coordinate(y + 1u, field.height(), field.dy()));
            const double tx = cx - trap.centerX;
            const double ty = cy - trap.centerY;
            if (tx * tx + ty * ty > maxRadius2) continue;
 
            const ComplexField2D::Value &p00 = field(x, y);
            const ComplexField2D::Value &p10 = field(x + 1u, y);
            const ComplexField2D::Value &p11 = field(x + 1u, y + 1u);
            const ComplexField2D::Value &p01 = field(x, y + 1u);
            const double windingPhase =
                wrapPhase(std::arg(p10) - std::arg(p00)) +
                wrapPhase(std::arg(p11) - std::arg(p10)) +
                wrapPhase(std::arg(p01) - std::arg(p11)) +
                wrapPhase(std::arg(p00) - std::arg(p01));
            const int winding = static_cast<int>(std::lround(windingPhase / (2.0 * kPi)));
            if (winding != expectedSign) continue;
 
            const double density = 0.25 * (
                static_cast<double>(std::norm(p00)) + static_cast<double>(std::norm(p10)) +
                static_cast<double>(std::norm(p11)) + static_cast<double>(std::norm(p01)));
            if (density < bestDensity) {
                bestDensity = density;
                best.valid = true;
                best.x = cx;
                best.y = cy;
                best.winding = winding;
                best.localDensity = density;
            }
        }
    }
    return best;
}
 
double QuantumVortex2D::angularMomentumZ(const ComplexField2D &field,
                                         double hbar,
                                         double originX,
                                         double originY)
{
    if (field.width() < 3u || field.height() < 3u || !(hbar > 0.0) || !finite(hbar)) return 0.0;
    const double area = static_cast<double>(field.dx()) * static_cast<double>(field.dy());
    const double inv2Dx = 1.0 / (2.0 * static_cast<double>(field.dx()));
    const double inv2Dy = 1.0 / (2.0 * static_cast<double>(field.dy()));
    double value = 0.0;
    double norm = 0.0;
 
    for (std::size_t y = 1u; y + 1u < field.height(); ++y) {
        const double py = coordinate(y, field.height(), field.dy()) - originY;
        for (std::size_t x = 1u; x + 1u < field.width(); ++x) {
            const double px = coordinate(x, field.width(), field.dx()) - originX;
            const std::complex<double> psi(field(x, y).real(), field(x, y).imag());
            const std::complex<double> dPsiDx =
                (std::complex<double>(field(x + 1u, y).real(), field(x + 1u, y).imag()) -
                 std::complex<double>(field(x - 1u, y).real(), field(x - 1u, y).imag())) * inv2Dx;
            const std::complex<double> dPsiDy =
                (std::complex<double>(field(x, y + 1u).real(), field(x, y + 1u).imag()) -
                 std::complex<double>(field(x, y - 1u).real(), field(x, y - 1u).imag())) * inv2Dy;
            const std::complex<double> lzPsi = std::complex<double>(0.0, -hbar) *
                (px * dPsiDy - py * dPsiDx);
            value += std::real(std::conj(psi) * lzPsi) * area;
            norm += std::norm(psi) * area;
        }
    }
    return norm > 1.0e-30 ? value / norm : 0.0;
}
 
void QuantumVortex2D::idealCorePosition(double time,
                                        const QuantumVortexInitialState2DConfig &config,
                                        const QuantumVortexTrap2D &trap,
                                        double mass,
                                        double &x,
                                        double &y)
{
    x = config.coreX;
    y = config.coreY;
    const double omega = trapAngularFrequency(trap, mass);
    if (!(omega > 0.0) || !finite(time)) return;
    const double dx = config.coreX - trap.centerX;
    const double dy = config.coreY - trap.centerY;
    const double theta = static_cast<double>(config.charge) * omega * time;
    const double c = std::cos(theta);
    const double s = std::sin(theta);
    x = trap.centerX + dx * c - dy * s;
    y = trap.centerY + dx * s + dy * c;
}
 
QuantumVortexHistory2D::QuantumVortexHistory2D(std::size_t maxSamples)
    : maxSamples_(std::max<std::size_t>(32u, maxSamples))
{
}
 
void QuantumVortexHistory2D::clear()
{
    configured_ = false;
    config_ = QuantumVortexInitialState2DConfig();
    trap_ = QuantumVortexTrap2D();
    samples_.clear();
}
 
void QuantumVortexHistory2D::configure(const QuantumVortexInitialState2DConfig &config,
                                       const QuantumVortexTrap2D &trap)
{
    clear();
    if (!trap.valid || !QuantumVortex2D::validateConfig(config, nullptr)) return;
    config_ = config;
    trap_ = trap;
    configured_ = true;
}
 
bool QuantumVortexHistory2D::configured() const noexcept { return configured_; }
bool QuantumVortexHistory2D::empty() const noexcept { return samples_.empty(); }
 
void QuantumVortexHistory2D::append(double time,
                                    const ComplexField2D &field,
                                    double hbar)
{
    if (!configured_ || field.empty() || !finite(time)) return;
    QuantumVortexSample2D sample;
    sample.time = time;
    sample.core = QuantumVortex2D::detectCore(field, config_, trap_);
    sample.angularMomentumZ = QuantumVortex2D::angularMomentumZ(field, hbar,
                                                               trap_.centerX, trap_.centerY);
    sample.norm = field.norm2();
    samples_.push_back(sample);
    trim();
}
 
const QuantumVortexInitialState2DConfig &QuantumVortexHistory2D::config() const noexcept { return config_; }
const QuantumVortexTrap2D &QuantumVortexHistory2D::trap() const noexcept { return trap_; }
const std::vector<QuantumVortexSample2D> &QuantumVortexHistory2D::samples() const noexcept { return samples_; }
 
void QuantumVortexHistory2D::trim()
{
    if (samples_.size() <= maxSamples_) return;
    const std::size_t removeCount = samples_.size() - maxSamples_;
    samples_.erase(samples_.begin(), samples_.begin() + static_cast<std::ptrdiff_t>(removeCount));
}
 
} // namespace qpe