momentumspace2d.cpp

Go to the documentation of this file.

/**
 * @file momentumspace2d.cpp
 * @brief Implements the momentumspace2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/momentumspace2d.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 "momentumspace2d.h"
 
#include "core/math/fft2d.h"
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <stdexcept>
 
namespace qpe {
namespace {
 
const double kPi = 3.141592653589793238462643383279502884;
 
/**
 * @brief Performs signed fft index.
 * @param index Zero-based index of the selected item.
 * @param count The count value used by the operation.
 */
long long signedFftIndex(std::size_t index, std::size_t count)
{
    const std::size_t half = count / 2u;
    if (index < half) {
        return static_cast<long long>(index);
    }
    return static_cast<long long>(index) - static_cast<long long>(count);
}
 
std::size_t shiftedSourceIndex(std::size_t shiftedIndex, std::size_t count)
{
    return (shiftedIndex + count / 2u) % count;
}
 
std::size_t nearestCoordinateIndex(const std::vector<double> &coordinates,
                                   double value)
{
    if (coordinates.empty()) {
        return 0u;
    }
 
    const std::vector<double>::const_iterator it =
        std::lower_bound(coordinates.begin(), coordinates.end(), value);
    if (it == coordinates.begin()) {
        return 0u;
    }
    if (it == coordinates.end()) {
        return coordinates.size() - 1u;
    }
 
    const std::size_t right = static_cast<std::size_t>(it - coordinates.begin());
    const std::size_t left = right - 1u;
    return std::fabs(value - coordinates[left]) <=
                   std::fabs(coordinates[right] - value)
               ? left
               : right;
}
 
/**
 * @brief Performs append sample.
 * @param value Value to assign or evaluate.
 * @param coordinate The coordinate value used by the operation.
 * @param out The out value used by the operation.
 */
void appendSample(const ComplexField2D::Value &value,
                  double coordinate,
                  WaveFunctionCrossSection2DData &out)
{
    const double re = static_cast<double>(value.real());
    const double im = static_cast<double>(value.imag());
    out.coordinate.push_back(coordinate);
    out.realPart.push_back(re);
    out.imaginaryPart.push_back(im);
    out.probability.push_back(re * re + im * im);
    out.phase.push_back(std::atan2(im, re));
}
 
} // namespace
 
MomentumSpace2DResult MomentumSpace2D::compute(const ComplexField2D &field,
                                               double hbar)
{
    if (field.empty()) {
        throw std::invalid_argument("MomentumSpace2D: field must not be empty");
    }
    if (!(hbar > 0.0) || !std::isfinite(hbar)) {
        throw std::invalid_argument("MomentumSpace2D: hbar must be finite and positive");
    }
    if (!FFT2D::supportsSize(field.width(), field.height())) {
        throw std::invalid_argument("MomentumSpace2D: FFT grid dimensions must be powers of two");
    }
 
    const std::size_t width = field.width();
    const std::size_t height = field.height();
    const double dx = static_cast<double>(field.dx());
    const double dy = static_cast<double>(field.dy());
    const double dpX = 2.0 * kPi * hbar / (static_cast<double>(width) * dx);
    const double dpY = 2.0 * kPi * hbar / (static_cast<double>(height) * dy);
    const double amplitudeScale = dx * dy / (2.0 * kPi * hbar);
 
    ComplexField2D transformed = field;
    FFT2D::forward(transformed);
 
    MomentumSpace2DResult out;
    out.hbar = hbar;
    out.spectrum.resize(width, height,
                        static_cast<float>(dpX),
                        static_cast<float>(dpY));
    out.pxCoordinates.resize(width);
    out.pyCoordinates.resize(height);
 
    for (std::size_t x = 0u; x < width; ++x) {
        const std::size_t sourceX = shiftedSourceIndex(x, width);
        out.pxCoordinates[x] =
            static_cast<double>(signedFftIndex(sourceX, width)) * dpX;
    }
    for (std::size_t y = 0u; y < height; ++y) {
        const std::size_t sourceY = shiftedSourceIndex(y, height);
        out.pyCoordinates[y] =
            static_cast<double>(signedFftIndex(sourceY, height)) * dpY;
    }
 
    const double xOriginIndex = 0.5 * static_cast<double>(width - 1u);
    const double yOriginIndex = 0.5 * static_cast<double>(height - 1u);
    double rawNorm = 0.0;
    double rawPx = 0.0;
    double rawPy = 0.0;
    double rawPx2 = 0.0;
    double rawPy2 = 0.0;
    double peakProbability = -1.0;
 
    for (std::size_t y = 0u; y < height; ++y) {
        const std::size_t sourceY = shiftedSourceIndex(y, height);
        const double py = out.pyCoordinates[y];
        for (std::size_t x = 0u; x < width; ++x) {
            const std::size_t sourceX = shiftedSourceIndex(x, width);
            const double px = out.pxCoordinates[x];
 
            // FFT2D assumes sample coordinate j*dx. The simulated spatial grid
            // is centred around zero, so restore the corresponding Fourier
            // phase. It does not change |psi_tilde|^2, but it makes the stored
            // complex spectrum physically consistent for future extensions.
            const double phase = 2.0 * kPi *
                (static_cast<double>(sourceX) * xOriginIndex / static_cast<double>(width) +
                 static_cast<double>(sourceY) * yOriginIndex / static_cast<double>(height));
            const std::complex<double> originPhase(std::cos(phase), std::sin(phase));
 
            const ComplexField2D::Value &src = transformed(sourceX, sourceY);
            const std::complex<double> value(
                static_cast<double>(src.real()),
                static_cast<double>(src.imag()));
            const std::complex<double> scaled = value * originPhase * amplitudeScale;
            out.spectrum(x, y) = ComplexField2D::Value(
                static_cast<float>(scaled.real()),
                static_cast<float>(scaled.imag()));
 
            const double probability = std::norm(scaled);
            const double weight = probability * dpX * dpY;
            rawNorm += weight;
            rawPx += px * weight;
            rawPy += py * weight;
            rawPx2 += px * px * weight;
            rawPy2 += py * py * weight;
 
            if (probability > peakProbability) {
                peakProbability = probability;
                out.peakPx = px;
                out.peakPy = py;
            }
        }
    }
 
    out.norm2 = rawNorm;
    if (rawNorm > 0.0 && std::isfinite(rawNorm)) {
        out.meanPx = rawPx / rawNorm;
        out.meanPy = rawPy / rawNorm;
        const double variancePx = std::max(0.0, rawPx2 / rawNorm - out.meanPx * out.meanPx);
        const double variancePy = std::max(0.0, rawPy2 / rawNorm - out.meanPy * out.meanPy);
        out.sigmaPx = std::sqrt(variancePx);
        out.sigmaPy = std::sqrt(variancePy);
    }
 
    return out;
}
 
WaveFunctionCrossSection2DData MomentumSpace2D::sample(
    const MomentumSpace2DResult &momentum,
    CrossSectionAxis axis,
    double position)
{
    WaveFunctionCrossSection2DData out;
    out.axis = axis;
    out.requestedPosition = position;
 
    if (momentum.empty() ||
        momentum.pxCoordinates.size() != momentum.spectrum.width() ||
        momentum.pyCoordinates.size() != momentum.spectrum.height()) {
        return out;
    }
 
    if (axis == CrossSectionAxis::HorizontalY) {
        const std::size_t y = nearestCoordinateIndex(momentum.pyCoordinates, position);
        out.actualPosition = momentum.pyCoordinates[y];
        out.coordinate.reserve(momentum.spectrum.width());
        out.probability.reserve(momentum.spectrum.width());
        out.realPart.reserve(momentum.spectrum.width());
        out.imaginaryPart.reserve(momentum.spectrum.width());
        out.phase.reserve(momentum.spectrum.width());
        for (std::size_t x = 0u; x < momentum.spectrum.width(); ++x) {
            appendSample(momentum.spectrum(x, y), momentum.pxCoordinates[x], out);
        }
    } else {
        const std::size_t x = nearestCoordinateIndex(momentum.pxCoordinates, position);
        out.actualPosition = momentum.pxCoordinates[x];
        out.coordinate.reserve(momentum.spectrum.height());
        out.probability.reserve(momentum.spectrum.height());
        out.realPart.reserve(momentum.spectrum.height());
        out.imaginaryPart.reserve(momentum.spectrum.height());
        out.phase.reserve(momentum.spectrum.height());
        for (std::size_t y = 0u; y < momentum.spectrum.height(); ++y) {
            appendSample(momentum.spectrum(x, y), momentum.pyCoordinates[y], out);
        }
    }
 
    return out;
}
 
} // namespace qpe