gaussianpacket2d_tests.cpp

Go to the documentation of this file.

/**
 * @file gaussianpacket2d_tests.cpp
 * @brief Implements regression and invariant tests for gaussianpacket2d tests.
 * @ingroup qpe_tests
 *
 * Source: @c tests/gaussianpacket2d_tests.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_tests.h"
 
#include "core/math/complexfield2d.h"
#include "core/quantum/gaussianpacket2d.h"
 
#include <cmath>
#include <complex>
#include <stdexcept>
 
namespace qpe {
namespace tests {
namespace {
 
/**
 * @brief Performs nearly equal.
 * @param a The a value used by the operation.
 * @param b The b value used by the operation.
 * @param tolerance The tolerance value used by the operation.
 */
bool nearlyEqual(double a, double b, double tolerance)
{
    return std::abs(a - b) <= tolerance;
}
 
/**
 * @brief Performs fail.
 * @param errorMessage The error message value used by the operation.
 * @param message The message value used by the operation.
 */
bool fail(std::string &errorMessage, const char *message)
{
    errorMessage = message;
    return false;
}
 
/**
 * @brief Performs wrapped angle difference.
 * @param a The a value used by the operation.
 * @param b The b value used by the operation.
 */
double wrappedAngleDifference(double a, double b)
{
    const double pi = 3.14159265358979323846;
    double d = a - b;
    while (d > pi) d -= 2.0 * pi;
    while (d < -pi) d += 2.0 * pi;
    return d;
}
 
} // namespace
 
/**
 * @brief Runs gaussian packet2 d tests.
 * @param errorMessage The error message value used by the operation.
 */
bool runGaussianPacket2DTests(std::string &errorMessage)
{
    errorMessage.clear();
 
    try {
        // Odd dimensions put an actual sample at x=y=0 and make the
        // expectation-value checks straightforward.
        ComplexField2D field(257, 257, 0.05f, 0.05f);
 
        GaussianPacket2DParams p;
        p.centerX = -1.25f;
        p.centerY = 0.70f;
        p.sigma = 0.60f;
        p.momentumX = 2.10f;
        p.momentumY = -0.80f;
        p.hbar = 1.0f;
 
        GaussianPacket2D::initialize(field, p);
 
        if (!nearlyEqual(field.norm2(), 1.0, 5.0e-6)) {
            return fail(errorMessage, "Gaussian packet is not normalized");
        }
 
        double meanX = 0.0;
        double meanY = 0.0;
        double meanX2 = 0.0;
        double meanY2 = 0.0;
        const double cellArea = static_cast<double>(field.dx()) *
                                static_cast<double>(field.dy());
 
        for (std::size_t y = 0; y < field.height(); ++y) {
            const double yy = GaussianPacket2D::yCoordinate(field, y);
            for (std::size_t x = 0; x < field.width(); ++x) {
                const double xx = GaussianPacket2D::xCoordinate(field, x);
                const double probability = static_cast<double>(std::norm(field(x, y)));
                const double w = probability * cellArea;
 
                meanX += xx * w;
                meanY += yy * w;
                meanX2 += xx * xx * w;
                meanY2 += yy * yy * w;
            }
        }
 
        if (!nearlyEqual(meanX, p.centerX, 2.0e-3) ||
            !nearlyEqual(meanY, p.centerY, 2.0e-3)) {
            return fail(errorMessage, "Gaussian packet centre is incorrect");
        }
 
        const double varX = meanX2 - meanX * meanX;
        const double varY = meanY2 - meanY * meanY;
        const double expectedVariance = static_cast<double>(p.sigma) * p.sigma;
 
        if (!nearlyEqual(varX, expectedVariance, 3.0e-3) ||
            !nearlyEqual(varY, expectedVariance, 3.0e-3)) {
            return fail(errorMessage, "Gaussian packet sigma is incorrect");
        }
 
        // Check that the phase increment along X is p_x * dx / hbar.
        const std::size_t cx = field.width() / 2;
        const std::size_t cy = field.height() / 2;
        const std::complex<float> a = field(cx, cy);
        const std::complex<float> b = field(cx + 1, cy);
        const double measuredPhaseStep = std::arg(b * std::conj(a));
        const double expectedPhaseStep =
            static_cast<double>(p.momentumX) * field.dx() / p.hbar;
 
        if (std::abs(wrappedAngleDifference(measuredPhaseStep, expectedPhaseStep)) > 2.0e-5) {
            return fail(errorMessage, "Gaussian packet X phase gradient is incorrect");
        }
 
        // Step 7.2 also exposes py, so verify the independent Y phase gradient.
        const std::complex<float> c = field(cx, cy + 1);
        const double measuredPhaseStepY = std::arg(c * std::conj(a));
        const double expectedPhaseStepY =
            static_cast<double>(p.momentumY) * field.dy() / p.hbar;
 
        if (std::abs(wrappedAngleDifference(measuredPhaseStepY, expectedPhaseStepY)) > 2.0e-5) {
            return fail(errorMessage, "Gaussian packet Y phase gradient is incorrect");
        }
 
        bool invalidSigmaRejected = false;
        try {
            GaussianPacket2DParams bad = p;
            bad.sigma = 0.0f;
            GaussianPacket2D::initialize(field, bad);
        } catch (const std::invalid_argument &) {
            invalidSigmaRejected = true;
        }
 
        if (!invalidSigmaRejected) {
            return fail(errorMessage, "invalid sigma was not rejected");
        }
    } catch (const std::exception &e) {
        errorMessage = e.what();
        return false;
    }
 
    return true;
}
 
} // namespace tests
} // namespace qpe