quantumkleinanglescan2d.cpp

Go to the documentation of this file.

/**
 * @file quantumkleinanglescan2d.cpp
 * @brief Implements the quantumkleinanglescan2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/quantumkleinanglescan2d.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 "quantumkleinanglescan2d.h"
 
#include <algorithm>
#include <array>
#include <cmath>
#include <complex>
 
namespace qpe {
namespace {
 
constexpr double kPi = 3.141592653589793238462643383279502884;
using C = std::complex<double>;
using Vec2 = std::array<C, 2>;
using Matrix4 = std::array<std::array<C, 4>, 4>;
using Vector4 = std::array<C, 4>;
 
/**
 * @brief Performs finite.
 * @param x X coordinate/value.
 */
bool finite(double x) noexcept { return std::isfinite(x); }
 
/**
 * @brief Performs eigen spinor.
 * @param px The px value used by the operation.
 * @param py The py value used by the operation.
 * @param relativeEnergy The relative energy value used by the operation.
 * @param lightSpeed The light speed value used by the operation.
 * @param restEnergy The rest energy value used by the operation.
 */
Vec2 eigenSpinor(const C &px, double py, double relativeEnergy,
                 double lightSpeed, double restEnergy) noexcept
{
    const C denom(relativeEnergy + restEnergy, 0.0);
    if (std::abs(denom) > 1.0e-12) {
        return Vec2{{C(1.0, 0.0), lightSpeed * (px + C(0.0, py)) / denom}};
    }
 
    const C alternate(relativeEnergy - restEnergy, 0.0);
    if (std::abs(alternate) > 1.0e-12) {
        return Vec2{{lightSpeed * (px - C(0.0, py)) / alternate, C(1.0, 0.0)}};
    }
    return Vec2{{C(1.0, 0.0), C(0.0, 0.0)}};
}
 
/**
 * @brief Performs solve4x4.
 * @param a The a value used by the operation.
 * @param b The b value used by the operation.
 * @param x X coordinate/value.
 */
bool solve4x4(Matrix4 a, Vector4 b, Vector4 &x) noexcept
{
    for (int col = 0; col < 4; ++col) {
        int pivot = col;
        double best = std::abs(a[static_cast<std::size_t>(col)][static_cast<std::size_t>(col)]);
        for (int row = col + 1; row < 4; ++row) {
            const double candidate = std::abs(a[static_cast<std::size_t>(row)][static_cast<std::size_t>(col)]);
            if (candidate > best) {
                best = candidate;
                pivot = row;
            }
        }
        if (!(best > 1.0e-14) || !finite(best)) return false;
        if (pivot != col) {
            std::swap(a[static_cast<std::size_t>(pivot)], a[static_cast<std::size_t>(col)]);
            std::swap(b[static_cast<std::size_t>(pivot)], b[static_cast<std::size_t>(col)]);
        }
 
        const C diag = a[static_cast<std::size_t>(col)][static_cast<std::size_t>(col)];
        for (int j = col; j < 4; ++j)
            a[static_cast<std::size_t>(col)][static_cast<std::size_t>(j)] /= diag;
        b[static_cast<std::size_t>(col)] /= diag;
 
        for (int row = 0; row < 4; ++row) {
            if (row == col) continue;
            const C f = a[static_cast<std::size_t>(row)][static_cast<std::size_t>(col)];
            if (std::abs(f) == 0.0) continue;
            for (int j = col; j < 4; ++j)
                a[static_cast<std::size_t>(row)][static_cast<std::size_t>(j)] -=
                    f * a[static_cast<std::size_t>(col)][static_cast<std::size_t>(j)];
            b[static_cast<std::size_t>(row)] -= f * b[static_cast<std::size_t>(col)];
        }
    }
    x = b;
    return true;
}
 
/**
 * @brief Performs rotated config.
 * @param base The base value used by the operation.
 * @param angleDegrees The angle degrees value used by the operation.
 */
QuantumKleinTunneling2DConfig rotatedConfig(const QuantumKleinTunneling2DConfig &base,
                                             double angleDegrees) noexcept
{
    QuantumKleinTunneling2DConfig c = base;
    const double p = QuantumKleinTunneling2D::incidentMomentumMagnitude(base);
    const double theta = angleDegrees * kPi / 180.0;
    c.momentumX = p * std::cos(theta);
    c.momentumY = p * std::sin(theta);
    return c;
}
 
} // namespace
 
QuantumKleinAnglePoint2D QuantumKleinAngleScan2D::evaluate(
    const QuantumKleinTunneling2DConfig &baseConfig, double angleDegrees, double hbar) noexcept
{
    QuantumKleinAnglePoint2D out;
    out.angleDegrees = angleDegrees;
 
    if (!QuantumKleinTunneling2D::validateConfig(baseConfig, nullptr) ||
        !finite(angleDegrees) || angleDegrees < 0.0 || angleDegrees >= 89.9 ||
        !finite(hbar) || hbar <= 0.0) {
        return out;
    }
 
    const QuantumKleinTunneling2DConfig c = rotatedConfig(baseConfig, angleDegrees);
    out.interiorLongitudinalMomentumSquared =
        QuantumKleinTunneling2D::interiorLongitudinalMomentumSquared(c);
    out.regime = QuantumKleinTunneling2D::classifyBarrierRegime(c);
 
    const double p = QuantumKleinTunneling2D::incidentMomentumMagnitude(c);
    const double theta = angleDegrees * kPi / 180.0;
    const double px = p * std::cos(theta);
    const double py = p * std::sin(theta);
    if (!(px > 1.0e-10)) return out;
 
    const double cc = c.lightSpeed;
    const double rest = c.restMass * cc * cc;
    const double e = QuantumKleinTunneling2D::incidentEnergy(c);
    const double er = e - c.barrierHeight;
    const C qx = std::sqrt(C(out.interiorLongitudinalMomentumSquared, 0.0));
 
    const Vec2 incident = eigenSpinor(C(px, 0.0), py, e, cc, rest);
    const Vec2 reflected = eigenSpinor(C(-px, 0.0), py, e, cc, rest);
    const Vec2 barrierForward = eigenSpinor(qx, py, er, cc, rest);
    const Vec2 barrierBackward = eigenSpinor(-qx, py, er, cc, rest);
    const Vec2 transmitted = incident;
 
    const C phaseBarrierForward = std::exp(C(0.0, 1.0) * qx * c.barrierWidth / hbar);
    const C phaseBarrierBackward = std::exp(C(0.0, -1.0) * qx * c.barrierWidth / hbar);
    const C phaseTransmitted = std::exp(C(0.0, px * c.barrierWidth / hbar));
 
    Matrix4 a{};
    Vector4 b{};
    for (int component = 0; component < 2; ++component) {
        const std::size_t k = static_cast<std::size_t>(component);
        a[k][0] = reflected[k];
        a[k][1] = -barrierForward[k];
        a[k][2] = -barrierBackward[k];
        a[k][3] = C(0.0, 0.0);
        b[k] = -incident[k];
 
        a[k + 2u][0] = C(0.0, 0.0);
        a[k + 2u][1] = barrierForward[k] * phaseBarrierForward;
        a[k + 2u][2] = barrierBackward[k] * phaseBarrierBackward;
        a[k + 2u][3] = -transmitted[k] * phaseTransmitted;
        b[k + 2u] = C(0.0, 0.0);
    }
 
    Vector4 solution{};
    if (!solve4x4(a, b, solution)) return out;
 
    const double r = std::norm(solution[0]);
    const double t = std::norm(solution[3]);
    if (!finite(r) || !finite(t)) return out;
 
    // Left and right asymptotic media are identical, so the Dirac flux ratio
    // reduces to |amplitude|^2 for the reflected/transmitted coefficients.
    out.reflection = std::max(0.0, std::min(1.0, r));
    out.transmission = std::max(0.0, std::min(1.0, t));
    out.valid = true;
    return out;
}
 
QuantumKleinAngleScan2DResult QuantumKleinAngleScan2D::scan(
    const QuantumKleinTunneling2DConfig &baseConfig,
    double minimumAngleDegrees, double maximumAngleDegrees,
    int samples, double hbar) noexcept
{
    QuantumKleinAngleScan2DResult out;
    if (!QuantumKleinTunneling2D::validateConfig(baseConfig, nullptr) ||
        !finite(minimumAngleDegrees) || !finite(maximumAngleDegrees) ||
        minimumAngleDegrees < 0.0 || maximumAngleDegrees <= minimumAngleDegrees ||
        maximumAngleDegrees >= 89.9 || samples < 2 || samples > 721 ||
        !finite(hbar) || hbar <= 0.0) {
        return out;
    }
 
    out.hbar = hbar;
    out.momentumMagnitude = QuantumKleinTunneling2D::incidentMomentumMagnitude(baseConfig);
    out.incidentEnergy = QuantumKleinTunneling2D::incidentEnergy(baseConfig);
    out.barrierHeight = baseConfig.barrierHeight;
    out.barrierWidth = baseConfig.barrierWidth;
    out.minimumAngleDegrees = minimumAngleDegrees;
    out.maximumAngleDegrees = maximumAngleDegrees;
    out.points.reserve(static_cast<std::size_t>(samples));
 
    bool allValid = true;
    for (int i = 0; i < samples; ++i) {
        const double f = static_cast<double>(i) / static_cast<double>(samples - 1);
        const double angle = minimumAngleDegrees + f * (maximumAngleDegrees - minimumAngleDegrees);
        const QuantumKleinAnglePoint2D point = evaluate(baseConfig, angle, hbar);
        allValid = allValid && point.valid;
        out.points.push_back(point);
    }
    out.valid = allValid && out.points.size() == static_cast<std::size_t>(samples);
    return out;
}
 
} // namespace qpe