squarebarrier1d.cpp

Go to the documentation of this file.

/**
 * @file squarebarrier1d.cpp
 * @brief Implements the squarebarrier1d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/squarebarrier1d.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 "squarebarrier1d.h"
 
#include <algorithm>
#include <cmath>
#include <stdexcept>
 
namespace qpe {
namespace {
 
/**
 * @brief Performs finite positive.
 * @param v The v value used by the operation.
 */
bool finitePositive(double v) noexcept
{
    return std::isfinite(v) && v > 0.0;
}
 
/**
 * @brief Performs finite non negative.
 * @param v The v value used by the operation.
 */
bool finiteNonNegative(double v) noexcept
{
    return std::isfinite(v) && v >= 0.0;
}
 
} // namespace
 
SquareBarrier1DResult SquareBarrier1D::evaluate(double energy,
                                                 double barrierHeight,
                                                 double barrierWidth,
                                                 double mass,
                                                 double hbar)
{
    if (!finiteNonNegative(energy) ||
        !finiteNonNegative(barrierHeight) ||
        !finiteNonNegative(barrierWidth) ||
        !finitePositive(mass) ||
        !finitePositive(hbar)) {
        throw std::invalid_argument("SquareBarrier1D: invalid parameters");
    }
 
    SquareBarrier1DResult result;
    result.energy = energy;
    result.barrierHeight = barrierHeight;
    result.barrierWidth = barrierWidth;
    result.tunnelingRegime = energy < barrierHeight;
 
    // Degenerate cases: no barrier, zero thickness, or a non-zero-energy
    // particle encountering V0=0 all transmit completely.
    if (barrierHeight == 0.0 || barrierWidth == 0.0) {
        result.transmission = 1.0;
        result.reflection = 0.0;
        return result;
    }
 
    // At exactly zero incident energy the stationary plane-wave flux tends to
    // zero and the finite positive barrier is fully reflecting in this limit.
    if (energy == 0.0) {
        result.transmission = 0.0;
        result.reflection = 1.0;
        return result;
    }
 
    const double hbar2 = hbar * hbar;
    double denominatorExtra = 0.0;
 
    const double scale = std::max({1.0, std::abs(energy), std::abs(barrierHeight)});
    const bool atTop = std::abs(energy - barrierHeight) <= 1.0e-12 * scale;
 
    if (atTop) {
        // Limit E -> V0 of either branch.
        denominatorExtra = mass * barrierHeight * barrierWidth * barrierWidth /
                           (2.0 * hbar2);
    } else if (energy < barrierHeight) {
        const double delta = barrierHeight - energy;
        const double kappa = std::sqrt(2.0 * mass * delta) / hbar;
        const double s = std::sinh(kappa * barrierWidth);
        denominatorExtra = barrierHeight * barrierHeight * s * s /
                           (4.0 * energy * delta);
    } else {
        const double delta = energy - barrierHeight;
        const double q = std::sqrt(2.0 * mass * delta) / hbar;
        const double s = std::sin(q * barrierWidth);
        denominatorExtra = barrierHeight * barrierHeight * s * s /
                           (4.0 * energy * delta);
    }
 
    const double t = 1.0 / (1.0 + denominatorExtra);
    result.transmission = std::max(0.0, std::min(1.0, t));
    result.reflection = 1.0 - result.transmission;
    return result;
}
 
} // namespace qpe