potentialfield2d.cpp

Go to the documentation of this file.

/**
 * @file potentialfield2d.cpp
 * @brief Implements the potentialfield2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/potentialfield2d.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 "potentialfield2d.h"
 
#include "core/math/complexfield2d.h"
 
#include <algorithm>
#include <cmath>
#include <limits>
#include <stdexcept>
 
namespace qpe {
namespace {
 
/**
 * @brief Performs finite.
 * @param value Value to assign or evaluate.
 */
bool finite(float value) noexcept
{
    return std::isfinite(static_cast<double>(value));
}
 
/**
 * @brief Performs finite positive.
 * @param value Value to assign or evaluate.
 */
bool finitePositive(float value) noexcept
{
    return finite(value) && value > 0.0f;
}
 
/**
 * @brief Performs require finite.
 * @param value Value to assign or evaluate.
 * @param message The message value used by the operation.
 */
void requireFinite(float value, const char *message)
{
    if (!finite(value)) {
        throw std::invalid_argument(message);
    }
}
 
/**
 * @brief Performs point in polygon.
 * @param x X coordinate/value.
 * @param y Y coordinate/value.
 * @param vertices The vertices value used by the operation.
 */
bool pointInPolygon(double x,
                    double y,
                    const std::vector<PotentialPoint2D> &vertices) noexcept
{
    bool inside = false;
    const std::size_t n = vertices.size();
 
    for (std::size_t i = 0, j = n - 1; i < n; j = i++) {
        const double xi = static_cast<double>(vertices[i].x);
        const double yi = static_cast<double>(vertices[i].y);
        const double xj = static_cast<double>(vertices[j].x);
        const double yj = static_cast<double>(vertices[j].y);
 
        const bool crosses = ((yi > y) != (yj > y));
        if (crosses) {
            const double xCross = (xj - xi) * (y - yi) / (yj - yi) + xi;
            if (x < xCross) {
                inside = !inside;
            }
        }
    }
 
    return inside;
}
 
} // namespace
 
PotentialField2D::PotentialField2D(std::size_t width,
                                   std::size_t height,
                                   float dx,
                                   float dy)
{
    resize(width, height, dx, dy);
}
 
void PotentialField2D::resize(std::size_t width,
                              std::size_t height,
                              float dx,
                              float dy)
{
    if (!finitePositive(dx) || !finitePositive(dy)) {
        throw std::invalid_argument("PotentialField2D: dx and dy must be finite and positive");
    }
 
    width_ = width;
    height_ = height;
    dx_ = dx;
    dy_ = dy;
    data_.assign(width_ * height_, 0.0f);
    touch();
}
 
void PotentialField2D::clear() noexcept
{
    fill(0.0f);
}
 
void PotentialField2D::fill(Value value) noexcept
{
    std::fill(data_.begin(), data_.end(), value);
    touch();
}
 
bool PotentialField2D::empty() const noexcept
{
    return data_.empty();
}
 
std::size_t PotentialField2D::width() const noexcept
{
    return width_;
}
 
std::size_t PotentialField2D::height() const noexcept
{
    return height_;
}
 
std::size_t PotentialField2D::size() const noexcept
{
    return data_.size();
}
 
float PotentialField2D::dx() const noexcept
{
    return dx_;
}
 
float PotentialField2D::dy() const noexcept
{
    return dy_;
}
 
PotentialField2D::Value PotentialField2D::operator()(std::size_t x,
                                                      std::size_t y) const noexcept
{
    return data_[index(x, y)];
}
 
PotentialField2D::Value PotentialField2D::at(std::size_t x, std::size_t y) const
{
    return data_[checkedIndex(x, y)];
}
 
void PotentialField2D::set(std::size_t x, std::size_t y, Value value) noexcept
{
    data_[index(x, y)] = value;
    touch();
}
 
void PotentialField2D::setAt(std::size_t x, std::size_t y, Value value)
{
    data_[checkedIndex(x, y)] = value;
    touch();
}
 
void PotentialField2D::accumulate(const PotentialField2D &other, Value scale)
{
    requireFinite(scale, "PotentialField2D::accumulate: scale must be finite");
    if (width_ != other.width_ || height_ != other.height_ ||
        dx_ != other.dx_ || dy_ != other.dy_) {
        throw std::invalid_argument("PotentialField2D::accumulate: grid mismatch");
    }
 
    for (std::size_t i = 0; i < data_.size(); ++i) {
        data_[i] += other.data_[i] * scale;
    }
    touch();
}
 
const PotentialField2D::Value *PotentialField2D::data() const noexcept
{
    return data_.data();
}
 
bool PotentialField2D::matchesGrid(const ComplexField2D &field) const noexcept
{
    return width_ == field.width() &&
           height_ == field.height() &&
           dx_ == field.dx() &&
           dy_ == field.dy();
}
 
PotentialField2D::Value PotentialField2D::minimum() const noexcept
{
    if (data_.empty()) {
        return 0.0f;
    }
    return *std::min_element(data_.begin(), data_.end());
}
 
PotentialField2D::Value PotentialField2D::maximum() const noexcept
{
    if (data_.empty()) {
        return 0.0f;
    }
    return *std::max_element(data_.begin(), data_.end());
}
 
std::uint64_t PotentialField2D::revision() const noexcept
{
    return revision_;
}
 
double PotentialField2D::xCoordinate(std::size_t x) const noexcept
{
    if (width_ == 0) {
        return 0.0;
    }
    const double centre = 0.5 * static_cast<double>(width_ - 1u);
    return (static_cast<double>(x) - centre) * static_cast<double>(dx_);
}
 
double PotentialField2D::yCoordinate(std::size_t y) const noexcept
{
    if (height_ == 0) {
        return 0.0;
    }
    const double centre = 0.5 * static_cast<double>(height_ - 1u);
    return (static_cast<double>(y) - centre) * static_cast<double>(dy_);
}
 
void PotentialField2D::addRectangle(float centerX,
                                    float centerY,
                                    float width,
                                    float height,
                                    Value value)
{
    requireFinite(centerX, "PotentialField2D::addRectangle: centerX must be finite");
    requireFinite(centerY, "PotentialField2D::addRectangle: centerY must be finite");
    requireFinite(value, "PotentialField2D::addRectangle: value must be finite");
    if (!finitePositive(width) || !finitePositive(height)) {
        throw std::invalid_argument("PotentialField2D::addRectangle: width and height must be finite and positive");
    }
 
    const double halfW = 0.5 * static_cast<double>(width);
    const double halfH = 0.5 * static_cast<double>(height);
 
    for (std::size_t y = 0; y < height_; ++y) {
        const double py = yCoordinate(y);
        if (std::abs(py - static_cast<double>(centerY)) > halfH) {
            continue;
        }
        for (std::size_t x = 0; x < width_; ++x) {
            const double px = xCoordinate(x);
            if (std::abs(px - static_cast<double>(centerX)) <= halfW) {
                data_[index(x, y)] = value;
            }
        }
    }
    touch();
}
 
void PotentialField2D::addCircle(float centerX,
                                 float centerY,
                                 float radius,
                                 Value value)
{
    requireFinite(centerX, "PotentialField2D::addCircle: centerX must be finite");
    requireFinite(centerY, "PotentialField2D::addCircle: centerY must be finite");
    requireFinite(value, "PotentialField2D::addCircle: value must be finite");
    if (!finitePositive(radius)) {
        throw std::invalid_argument("PotentialField2D::addCircle: radius must be finite and positive");
    }
 
    const double r2 = static_cast<double>(radius) * static_cast<double>(radius);
    for (std::size_t y = 0; y < height_; ++y) {
        const double dy = yCoordinate(y) - static_cast<double>(centerY);
        for (std::size_t x = 0; x < width_; ++x) {
            const double dx = xCoordinate(x) - static_cast<double>(centerX);
            if (dx * dx + dy * dy <= r2) {
                data_[index(x, y)] = value;
            }
        }
    }
    touch();
}
 
void PotentialField2D::addBarrier(float x1,
                                  float y1,
                                  float x2,
                                  float y2,
                                  float thickness,
                                  Value value)
{
    requireFinite(x1, "PotentialField2D::addBarrier: x1 must be finite");
    requireFinite(y1, "PotentialField2D::addBarrier: y1 must be finite");
    requireFinite(x2, "PotentialField2D::addBarrier: x2 must be finite");
    requireFinite(y2, "PotentialField2D::addBarrier: y2 must be finite");
    requireFinite(value, "PotentialField2D::addBarrier: value must be finite");
    if (!finitePositive(thickness)) {
        throw std::invalid_argument("PotentialField2D::addBarrier: thickness must be finite and positive");
    }
 
    const double ax = static_cast<double>(x1);
    const double ay = static_cast<double>(y1);
    const double bx = static_cast<double>(x2);
    const double by = static_cast<double>(y2);
    const double vx = bx - ax;
    const double vy = by - ay;
    const double length2 = vx * vx + vy * vy;
    const double halfT = 0.5 * static_cast<double>(thickness);
    const double limit2 = halfT * halfT;
 
    for (std::size_t y = 0; y < height_; ++y) {
        const double py = yCoordinate(y);
        for (std::size_t x = 0; x < width_; ++x) {
            const double px = xCoordinate(x);
 
            double cx = ax;
            double cy = ay;
            if (length2 > std::numeric_limits<double>::epsilon()) {
                double t = ((px - ax) * vx + (py - ay) * vy) / length2;
                t = std::max(0.0, std::min(1.0, t));
                cx = ax + t * vx;
                cy = ay + t * vy;
            }
 
            const double dx = px - cx;
            const double dy = py - cy;
            if (dx * dx + dy * dy <= limit2) {
                data_[index(x, y)] = value;
            }
        }
    }
    touch();
}
 
void PotentialField2D::addPolygon(const std::vector<PotentialPoint2D> &vertices,
                                  Value value)
{
    requireFinite(value, "PotentialField2D::addPolygon: value must be finite");
    if (vertices.size() < 3) {
        throw std::invalid_argument("PotentialField2D::addPolygon: at least three vertices are required");
    }
    for (const PotentialPoint2D &p : vertices) {
        if (!finite(p.x) || !finite(p.y)) {
            throw std::invalid_argument("PotentialField2D::addPolygon: vertices must be finite");
        }
    }
 
    for (std::size_t y = 0; y < height_; ++y) {
        const double py = yCoordinate(y);
        for (std::size_t x = 0; x < width_; ++x) {
            const double px = xCoordinate(x);
            if (pointInPolygon(px, py, vertices)) {
                data_[index(x, y)] = value;
            }
        }
    }
    touch();
}
 
void PotentialField2D::addSingleSlit(float centerX,
                                     float barrierWidth,
                                     float slitCenterY,
                                     float slitHeight,
                                     Value value)
{
    requireFinite(centerX, "PotentialField2D::addSingleSlit: centerX must be finite");
    requireFinite(slitCenterY, "PotentialField2D::addSingleSlit: slitCenterY must be finite");
    requireFinite(value, "PotentialField2D::addSingleSlit: value must be finite");
    if (!finitePositive(barrierWidth) || !finitePositive(slitHeight)) {
        throw std::invalid_argument("PotentialField2D::addSingleSlit: dimensions must be finite and positive");
    }
 
    const double halfW = 0.5 * static_cast<double>(barrierWidth);
    const double halfSlit = 0.5 * static_cast<double>(slitHeight);
    for (std::size_t y = 0; y < height_; ++y) {
        const double py = yCoordinate(y);
        const bool inOpening = std::abs(py - static_cast<double>(slitCenterY)) <= halfSlit;
        if (inOpening) {
            continue;
        }
        for (std::size_t x = 0; x < width_; ++x) {
            const double px = xCoordinate(x);
            if (std::abs(px - static_cast<double>(centerX)) <= halfW) {
                data_[index(x, y)] = value;
            }
        }
    }
    touch();
}
 
void PotentialField2D::addDoubleSlit(float centerX,
                                     float barrierWidth,
                                     float slitCenterY,
                                     float slitSeparation,
                                     float slitHeight,
                                     Value value)
{
    requireFinite(centerX, "PotentialField2D::addDoubleSlit: centerX must be finite");
    requireFinite(slitCenterY, "PotentialField2D::addDoubleSlit: slitCenterY must be finite");
    requireFinite(value, "PotentialField2D::addDoubleSlit: value must be finite");
    if (!finitePositive(barrierWidth) || !finitePositive(slitSeparation) || !finitePositive(slitHeight)) {
        throw std::invalid_argument("PotentialField2D::addDoubleSlit: dimensions must be finite and positive");
    }
 
    const double halfW = 0.5 * static_cast<double>(barrierWidth);
    const double halfSlit = 0.5 * static_cast<double>(slitHeight);
    const double halfSep = 0.5 * static_cast<double>(slitSeparation);
    const double slit1 = static_cast<double>(slitCenterY) - halfSep;
    const double slit2 = static_cast<double>(slitCenterY) + halfSep;
 
    for (std::size_t y = 0; y < height_; ++y) {
        const double py = yCoordinate(y);
        const bool inOpening = std::abs(py - slit1) <= halfSlit ||
                               std::abs(py - slit2) <= halfSlit;
        if (inOpening) {
            continue;
        }
        for (std::size_t x = 0; x < width_; ++x) {
            const double px = xCoordinate(x);
            if (std::abs(px - static_cast<double>(centerX)) <= halfW) {
                data_[index(x, y)] = value;
            }
        }
    }
    touch();
}
 
 
void PotentialField2D::addGaussian(float centerX,
                                   float centerY,
                                   float sigmaX,
                                   float sigmaY,
                                   Value amplitude)
{
    requireFinite(centerX, "PotentialField2D::addGaussian: centerX must be finite");
    requireFinite(centerY, "PotentialField2D::addGaussian: centerY must be finite");
    requireFinite(amplitude, "PotentialField2D::addGaussian: amplitude must be finite");
    if (!finitePositive(sigmaX) || !finitePositive(sigmaY)) {
        throw std::invalid_argument("PotentialField2D::addGaussian: sigmas must be finite and positive");
    }
 
    const double inv2sx2 = 1.0 / (2.0 * static_cast<double>(sigmaX) * static_cast<double>(sigmaX));
    const double inv2sy2 = 1.0 / (2.0 * static_cast<double>(sigmaY) * static_cast<double>(sigmaY));
    for (std::size_t y = 0; y < height_; ++y) {
        const double dy = yCoordinate(y) - static_cast<double>(centerY);
        for (std::size_t x = 0; x < width_; ++x) {
            const double dx = xCoordinate(x) - static_cast<double>(centerX);
            data_[index(x, y)] += static_cast<Value>(
                static_cast<double>(amplitude) * std::exp(-(dx * dx * inv2sx2 + dy * dy * inv2sy2)));
        }
    }
    touch();
}
 
void PotentialField2D::addHarmonic(float centerX,
                                   float centerY,
                                   Value offset,
                                   float curvatureX,
                                   float curvatureY)
{
    requireFinite(centerX, "PotentialField2D::addHarmonic: centerX must be finite");
    requireFinite(centerY, "PotentialField2D::addHarmonic: centerY must be finite");
    requireFinite(offset, "PotentialField2D::addHarmonic: offset must be finite");
    requireFinite(curvatureX, "PotentialField2D::addHarmonic: curvatureX must be finite");
    requireFinite(curvatureY, "PotentialField2D::addHarmonic: curvatureY must be finite");
 
    for (std::size_t y = 0; y < height_; ++y) {
        const double dy = yCoordinate(y) - static_cast<double>(centerY);
        for (std::size_t x = 0; x < width_; ++x) {
            const double dx = xCoordinate(x) - static_cast<double>(centerX);
            data_[index(x, y)] += static_cast<Value>(
                static_cast<double>(offset) +
                0.5 * static_cast<double>(curvatureX) * dx * dx +
                0.5 * static_cast<double>(curvatureY) * dy * dy);
        }
    }
    touch();
}
 
void PotentialField2D::addQuadrupole(float centerX,
                                      float centerY,
                                      Value q1,
                                      Value q2)
{
    requireFinite(centerX, "PotentialField2D::addQuadrupole: centerX must be finite");
    requireFinite(centerY, "PotentialField2D::addQuadrupole: centerY must be finite");
    requireFinite(q1, "PotentialField2D::addQuadrupole: q1 must be finite");
    requireFinite(q2, "PotentialField2D::addQuadrupole: q2 must be finite");
 
    for (std::size_t y = 0; y < height_; ++y) {
        const double dy = yCoordinate(y) - static_cast<double>(centerY);
        for (std::size_t x = 0; x < width_; ++x) {
            const double dx = xCoordinate(x) - static_cast<double>(centerX);
            data_[index(x, y)] += static_cast<Value>(
                static_cast<double>(q1) * (dx * dx - dy * dy) +
                2.0 * static_cast<double>(q2) * dx * dy);
        }
    }
    touch();
}
 
void PotentialField2D::addStandingWaveX(float period,
                                           Value depth,
                                           float phaseDegrees)
{
    if (!finitePositive(period)) {
        throw std::invalid_argument("PotentialField2D::addStandingWaveX: period must be finite and positive");
    }
    requireFinite(depth, "PotentialField2D::addStandingWaveX: depth must be finite");
    requireFinite(phaseDegrees, "PotentialField2D::addStandingWaveX: phase must be finite");
 
    const double pi = 3.141592653589793238462643383279502884;
    const double halfPhase = 0.5 * static_cast<double>(phaseDegrees) * pi / 180.0;
    const double wave = pi / static_cast<double>(period);
    for (std::size_t x = 0; x < width_; ++x) {
        const double argument = wave * xCoordinate(x) + halfPhase;
        const double c = std::cos(argument);
        const Value v = static_cast<Value>(static_cast<double>(depth) * c * c);
        for (std::size_t y = 0; y < height_; ++y) {
            data_[index(x, y)] += v;
        }
    }
    touch();
}
 
void PotentialField2D::addLinear(float centerX,
                                 float centerY,
                                 Value offset,
                                 float slopeX,
                                 float slopeY)
{
    requireFinite(centerX, "PotentialField2D::addLinear: centerX must be finite");
    requireFinite(centerY, "PotentialField2D::addLinear: centerY must be finite");
    requireFinite(offset, "PotentialField2D::addLinear: offset must be finite");
    requireFinite(slopeX, "PotentialField2D::addLinear: slopeX must be finite");
    requireFinite(slopeY, "PotentialField2D::addLinear: slopeY must be finite");
 
    for (std::size_t y = 0; y < height_; ++y) {
        const double dy = yCoordinate(y) - static_cast<double>(centerY);
        for (std::size_t x = 0; x < width_; ++x) {
            const double dx = xCoordinate(x) - static_cast<double>(centerX);
            data_[index(x, y)] += static_cast<Value>(
                static_cast<double>(offset) +
                static_cast<double>(slopeX) * dx +
                static_cast<double>(slopeY) * dy);
        }
    }
    touch();
}
 
void PotentialField2D::addSoftenedRadial(float centerX,
                                         float centerY,
                                         Value coupling,
                                         float softening)
{
    requireFinite(centerX, "PotentialField2D::addSoftenedRadial: centerX must be finite");
    requireFinite(centerY, "PotentialField2D::addSoftenedRadial: centerY must be finite");
    requireFinite(coupling, "PotentialField2D::addSoftenedRadial: coupling must be finite");
    if (!finitePositive(softening)) {
        throw std::invalid_argument("PotentialField2D::addSoftenedRadial: softening must be finite and positive");
    }
 
    const double eps2 = static_cast<double>(softening) * static_cast<double>(softening);
    for (std::size_t y = 0; y < height_; ++y) {
        const double dy = yCoordinate(y) - static_cast<double>(centerY);
        for (std::size_t x = 0; x < width_; ++x) {
            const double dx = xCoordinate(x) - static_cast<double>(centerX);
            data_[index(x, y)] += static_cast<Value>(
                static_cast<double>(coupling) / std::sqrt(dx * dx + dy * dy + eps2));
        }
    }
    touch();
}
 
std::size_t PotentialField2D::index(std::size_t x, std::size_t y) const noexcept
{
    return y * width_ + x;
}
 
std::size_t PotentialField2D::checkedIndex(std::size_t x, std::size_t y) const
{
    if (x >= width_ || y >= height_) {
        throw std::out_of_range("PotentialField2D: coordinates outside field");
    }
    return index(x, y);
}
 
void PotentialField2D::touch() noexcept
{
    ++revision_;
    if (revision_ == 0) {
        revision_ = 1;
    }
}
 
void PotentialField2D::addStadiumExterior(float centerX,
                                          float centerY,
                                          float straightLength,
                                          float radius,
                                          Value value)
{
    requireFinite(centerX, "PotentialField2D::addStadiumExterior: centerX must be finite");
    requireFinite(centerY, "PotentialField2D::addStadiumExterior: centerY must be finite");
    requireFinite(value, "PotentialField2D::addStadiumExterior: value must be finite");
    if (!finitePositive(straightLength) || !finitePositive(radius)) {
        throw std::invalid_argument("PotentialField2D::addStadiumExterior: dimensions must be finite and positive");
    }
 
    const double halfStraight = 0.5 * static_cast<double>(straightLength);
    const double r = static_cast<double>(radius);
    const double r2 = r * r;
 
    for (std::size_t y = 0; y < height_; ++y) {
        const double ly = yCoordinate(y) - static_cast<double>(centerY);
        for (std::size_t x = 0; x < width_; ++x) {
            const double lx = xCoordinate(x) - static_cast<double>(centerX);
            bool inside = false;
            if (std::abs(lx) <= halfStraight) {
                inside = std::abs(ly) <= r;
            } else {
                const double capX = lx > 0.0 ? halfStraight : -halfStraight;
                const double dx = lx - capX;
                inside = dx * dx + ly * ly <= r2;
            }
            if (!inside) {
                data_[index(x, y)] = value;
            }
        }
    }
    touch();
}
 
} // namespace qpe