probabilitycurrent2d.cpp

Go to the documentation of this file.

/**
 * @file probabilitycurrent2d.cpp
 * @brief Implements the probabilitycurrent2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/probabilitycurrent2d.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 "probabilitycurrent2d.h"
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <stdexcept>
 
namespace qpe {
 
bool ProbabilityCurrentField2D::empty() const noexcept
{
    const std::size_t count = width * height;
    return width == 0 || height == 0 ||
           jx.size() != count || jy.size() != count;
}
 
ProbabilityCurrentField2D ProbabilityCurrent2D::compute(const ComplexField2D &field,
                                                         float hbar,
                                                         float mass)
{
    if (field.empty()) {
        return ProbabilityCurrentField2D();
    }
    if (!(hbar > 0.0f) || !std::isfinite(static_cast<double>(hbar))) {
        throw std::invalid_argument("ProbabilityCurrent2D: hbar must be finite and positive");
    }
    if (!(mass > 0.0f) || !std::isfinite(static_cast<double>(mass))) {
        throw std::invalid_argument("ProbabilityCurrent2D: mass must be finite and positive");
    }
 
    ProbabilityCurrentField2D out;
    out.width = field.width();
    out.height = field.height();
    out.dx = field.dx();
    out.dy = field.dy();
    out.jx.assign(field.size(), 0.0f);
    out.jy.assign(field.size(), 0.0f);
 
    if (field.width() < 3 || field.height() < 3) {
        return out;
    }
 
    const double inv2Dx = 1.0 / (2.0 * static_cast<double>(field.dx()));
    const double inv2Dy = 1.0 / (2.0 * static_cast<double>(field.dy()));
    const double prefactor = static_cast<double>(hbar) / static_cast<double>(mass);
    const double cellArea = static_cast<double>(field.dx()) * static_cast<double>(field.dy());
 
    for (std::size_t y = 1; y + 1 < field.height(); ++y) {
        for (std::size_t x = 1; x + 1 < field.width(); ++x) {
            const std::complex<double> psi(
                static_cast<double>(field(x, y).real()),
                static_cast<double>(field(x, y).imag()));
 
            const std::complex<double> psiXp(
                static_cast<double>(field(x + 1, y).real()),
                static_cast<double>(field(x + 1, y).imag()));
            const std::complex<double> psiXm(
                static_cast<double>(field(x - 1, y).real()),
                static_cast<double>(field(x - 1, y).imag()));
            const std::complex<double> psiYp(
                static_cast<double>(field(x, y + 1).real()),
                static_cast<double>(field(x, y + 1).imag()));
            const std::complex<double> psiYm(
                static_cast<double>(field(x, y - 1).real()),
                static_cast<double>(field(x, y - 1).imag()));
 
            const std::complex<double> dPsiDx = (psiXp - psiXm) * inv2Dx;
            const std::complex<double> dPsiDy = (psiYp - psiYm) * inv2Dy;
 
            const double currentX = prefactor * std::imag(std::conj(psi) * dPsiDx);
            const double currentY = prefactor * std::imag(std::conj(psi) * dPsiDy);
            const std::size_t index = y * field.width() + x;
 
            out.jx[index] = static_cast<float>(currentX);
            out.jy[index] = static_cast<float>(currentY);
 
            const float magnitude = static_cast<float>(std::sqrt(currentX * currentX + currentY * currentY));
            out.maxMagnitude = std::max(out.maxMagnitude, magnitude);
            out.integratedX += currentX * cellArea;
            out.integratedY += currentY * cellArea;
        }
    }
 
    return out;
}
 
ProbabilityCurrentField2D ProbabilityCurrent2D::computeGaugeCovariant(const ComplexField2D &field,
                                                                       const U1GaugeField2D &gauge,
                                                                       float hbar,
                                                                       float mass)
{
    if (field.empty()) return ProbabilityCurrentField2D();
    if (!(hbar > 0.0f) || !std::isfinite(static_cast<double>(hbar)))
        throw std::invalid_argument("ProbabilityCurrent2D: hbar must be finite and positive");
    if (!(mass > 0.0f) || !std::isfinite(static_cast<double>(mass)))
        throw std::invalid_argument("ProbabilityCurrent2D: mass must be finite and positive");
    if (!gauge.matchesGrid(field))
        throw std::invalid_argument("ProbabilityCurrent2D: gauge field does not match wavefunction grid");
 
    ProbabilityCurrentField2D out;
    out.width = field.width(); out.height = field.height();
    out.dx = field.dx(); out.dy = field.dy();
    out.jx.assign(field.size(), 0.0f); out.jy.assign(field.size(), 0.0f);
    if (field.width() < 3 || field.height() < 3) return out;
 
    const double inv2Dx = 1.0 / (2.0 * static_cast<double>(field.dx()));
    const double inv2Dy = 1.0 / (2.0 * static_cast<double>(field.dy()));
    const double prefactor = static_cast<double>(hbar) / static_cast<double>(mass);
    const double cellArea = static_cast<double>(field.dx()) * static_cast<double>(field.dy());
 
    for (std::size_t y = 1; y + 1 < field.height(); ++y) {
        for (std::size_t x = 1; x + 1 < field.width(); ++x) {
            const std::complex<double> psi(field(x, y));
            const std::complex<double> right = std::complex<double>(gauge.linkX(x, y)) *
                                               std::complex<double>(field(x + 1u, y));
            const std::complex<double> left = std::conj(std::complex<double>(gauge.linkX(x - 1u, y))) *
                                              std::complex<double>(field(x - 1u, y));
            const std::complex<double> up = std::complex<double>(gauge.linkY(x, y)) *
                                            std::complex<double>(field(x, y + 1u));
            const std::complex<double> down = std::conj(std::complex<double>(gauge.linkY(x, y - 1u))) *
                                              std::complex<double>(field(x, y - 1u));
            const std::complex<double> dPsiDx = (right - left) * inv2Dx;
            const std::complex<double> dPsiDy = (up - down) * inv2Dy;
            const double currentX = prefactor * std::imag(std::conj(psi) * dPsiDx);
            const double currentY = prefactor * std::imag(std::conj(psi) * dPsiDy);
            const std::size_t index = y * field.width() + x;
            out.jx[index] = static_cast<float>(currentX);
            out.jy[index] = static_cast<float>(currentY);
            out.maxMagnitude = std::max(out.maxMagnitude,
                static_cast<float>(std::sqrt(currentX * currentX + currentY * currentY)));
            out.integratedX += currentX * cellArea;
            out.integratedY += currentY * cellArea;
        }
    }
    return out;
}
 
} // namespace qpe