gaussianpacket2d.cpp

Go to the documentation of this file.

/**
 * @file gaussianpacket2d.cpp
 * @brief Implements the gaussianpacket2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/gaussianpacket2d.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 "gaussianpacket2d.h"
 
#include <cmath>
#include <stdexcept>
 
namespace qpe {
namespace {
 
/**
 * @brief Performs finite.
 * @param v The v value used by the operation.
 */
bool finite(float v) noexcept
{
    return std::isfinite(static_cast<double>(v));
}
 
} // namespace
 
double GaussianPacket2D::xCoordinate(const ComplexField2D &field, std::size_t x) noexcept
{
    if (field.width() == 0) {
        return 0.0;
    }
    const double centre = 0.5 * static_cast<double>(field.width() - 1u);
    return (static_cast<double>(x) - centre) * static_cast<double>(field.dx());
}
 
double GaussianPacket2D::yCoordinate(const ComplexField2D &field, std::size_t y) noexcept
{
    if (field.height() == 0) {
        return 0.0;
    }
    const double centre = 0.5 * static_cast<double>(field.height() - 1u);
    return (static_cast<double>(y) - centre) * static_cast<double>(field.dy());
}
 
void GaussianPacket2D::initialize(ComplexField2D &field,
                                  const GaussianPacket2DParams &params)
{
    if (field.empty()) {
        throw std::invalid_argument("GaussianPacket2D: field must not be empty");
    }
 
    if (!(params.sigma > 0.0f) || !finite(params.sigma)) {
        throw std::invalid_argument("GaussianPacket2D: sigma must be finite and positive");
    }
 
    if (!(params.hbar > 0.0f) || !finite(params.hbar)) {
        throw std::invalid_argument("GaussianPacket2D: hbar must be finite and positive");
    }
 
    if (!finite(params.centerX) || !finite(params.centerY) ||
        !finite(params.momentumX) || !finite(params.momentumY)) {
        throw std::invalid_argument("GaussianPacket2D: all parameters must be finite");
    }
 
    const double sigma = static_cast<double>(params.sigma);
    const double inv4Sigma2 = 1.0 / (4.0 * sigma * sigma);
    const double invHbar = 1.0 / static_cast<double>(params.hbar);
 
    for (std::size_t y = 0; y < field.height(); ++y) {
        const double py = yCoordinate(field, y);
        const double dy = py - static_cast<double>(params.centerY);
 
        for (std::size_t x = 0; x < field.width(); ++x) {
            const double px = xCoordinate(field, x);
            const double dx = px - static_cast<double>(params.centerX);
 
            const double envelope = std::exp(-(dx * dx + dy * dy) * inv4Sigma2);
            const double phase =
                (static_cast<double>(params.momentumX) * px +
                 static_cast<double>(params.momentumY) * py) * invHbar;
 
            field(x, y) = ComplexField2D::Value(
                static_cast<float>(envelope * std::cos(phase)),
                static_cast<float>(envelope * std::sin(phase)));
        }
    }
 
    if (!field.normalize()) {
        throw std::runtime_error("GaussianPacket2D: normalization failed");
    }
}
 
} // namespace qpe