energyobservables2d.cpp

Go to the documentation of this file.

/**
 * @file energyobservables2d.cpp
 * @brief Implements the energyobservables2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/energyobservables2d.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 "energyobservables2d.h"
 
#include <cmath>
#include <stdexcept>
 
namespace qpe {
namespace {
 
std::size_t wrapIndex(long long index, std::size_t count) noexcept
{
    const long long n = static_cast<long long>(count);
    long long wrapped = index % n;
    if (wrapped < 0) {
        wrapped += n;
    }
    return static_cast<std::size_t>(wrapped);
}
 
std::complex<double> sample(const ComplexField2D &field,
                            long long x,
                            long long y) noexcept
{
    const ComplexField2D::Value &v = field(
        wrapIndex(x, field.width()),
        wrapIndex(y, field.height()));
    return std::complex<double>(static_cast<double>(v.real()),
                                static_cast<double>(v.imag()));
}
 
std::complex<double> derivativeX4(const ComplexField2D &field,
                                  std::size_t x,
                                  std::size_t y) noexcept
{
    const long long xi = static_cast<long long>(x);
    const long long yi = static_cast<long long>(y);
    const double h = static_cast<double>(field.dx());
    return (-sample(field, xi + 2, yi)
            + 8.0 * sample(field, xi + 1, yi)
            - 8.0 * sample(field, xi - 1, yi)
            + sample(field, xi - 2, yi)) / (12.0 * h);
}
 
std::complex<double> derivativeY4(const ComplexField2D &field,
                                  std::size_t x,
                                  std::size_t y) noexcept
{
    const long long xi = static_cast<long long>(x);
    const long long yi = static_cast<long long>(y);
    const double h = static_cast<double>(field.dy());
    return (-sample(field, xi, yi + 2)
            + 8.0 * sample(field, xi, yi + 1)
            - 8.0 * sample(field, xi, yi - 1)
            + sample(field, xi, yi - 2)) / (12.0 * h);
}
 
} // namespace
 
EnergyObservables2DResult EnergyObservables2D::compute(
    const ComplexField2D &field,
    const PotentialField2D &potential,
    double hbar,
    double mass)
{
    if (field.empty()) {
        throw std::invalid_argument("EnergyObservables2D: field must not be empty");
    }
    if (!potential.matchesGrid(field)) {
        throw std::invalid_argument("EnergyObservables2D: potential grid must match psi");
    }
    if (!(hbar > 0.0) || !std::isfinite(hbar)) {
        throw std::invalid_argument("EnergyObservables2D: hbar must be finite and positive");
    }
    if (!(mass > 0.0) || !std::isfinite(mass)) {
        throw std::invalid_argument("EnergyObservables2D: mass must be finite and positive");
    }
    if (field.width() < 5u || field.height() < 5u) {
        throw std::invalid_argument("EnergyObservables2D: grid must be at least 5x5");
    }
 
    EnergyObservables2DResult out;
    const double cellArea = static_cast<double>(field.dx()) *
                            static_cast<double>(field.dy());
    const double kineticScale = (hbar * hbar) / (2.0 * mass);
 
    double kineticRaw = 0.0;
    double potentialRaw = 0.0;
 
    for (std::size_t y = 0u; y < field.height(); ++y) {
        for (std::size_t x = 0u; x < field.width(); ++x) {
            const ComplexField2D::Value &vf = field(x, y);
            const std::complex<double> psi(static_cast<double>(vf.real()),
                                           static_cast<double>(vf.imag()));
            const double probabilityDensity = std::norm(psi);
            out.norm2 += probabilityDensity * cellArea;
 
            const std::complex<double> dxPsi = derivativeX4(field, x, y);
            const std::complex<double> dyPsi = derivativeY4(field, x, y);
            kineticRaw += kineticScale *
                          (std::norm(dxPsi) + std::norm(dyPsi)) * cellArea;
            potentialRaw += static_cast<double>(potential(x, y)) *
                            probabilityDensity * cellArea;
        }
    }
 
    if (!(out.norm2 > 0.0) || !std::isfinite(out.norm2)) {
        return out;
    }
 
    out.kinetic = kineticRaw / out.norm2;
    out.potential = potentialRaw / out.norm2;
    out.total = out.kinetic + out.potential;
    return out;
}
 
} // namespace qpe