blochoscillations2d.cpp

Go to the documentation of this file.

/**
 * @file blochoscillations2d.cpp
 * @brief Implements the blochoscillations2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/blochoscillations2d.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 "blochoscillations2d.h"
 
#include <algorithm>
#include <cmath>
#include <vector>
 
namespace qpe {
namespace {
 
constexpr double kPi = 3.141592653589793238462643383279502884;
 
/**
 * @brief Performs near.
 * @param a The a value used by the operation.
 * @param b The b value used by the operation.
 * @param scale The scale value used by the operation.
 */
bool near(double a, double b, double scale = 1.0) noexcept
{
    const double tolerance = 1.0e-4 * std::max(1.0, scale);
    return std::fabs(a - b) <= tolerance;
}
 
/**
 * @brief Performs fail.
 * @param reason The reason value used by the operation.
 * @param message The message value used by the operation.
 */
void fail(std::string *reason, const char *message) noexcept
{
    if (reason) {
        *reason = message;
    }
}
 
} // namespace
 
bool BlochOscillations2D::identify(const PotentialScene2D &scene,
                                   BlochLatticeGeometry2D &geometry,
                                   std::string *reason) noexcept
{
    geometry = BlochLatticeGeometry2D();
    if (reason) reason->clear();
 
    std::vector<const PotentialObject2D *> wells;
    const PotentialObject2D *linear = nullptr;
 
    for (const PotentialObject2D &object : scene.objects()) {
        if (!object.enabled) continue;
 
        if (object.type == PotentialObjectType2D::Rectangle && object.value < 0.0f) {
            wells.push_back(&object);
            continue;
        }
 
        if (object.type == PotentialObjectType2D::Linear) {
            if (linear != nullptr) {
                fail(reason, "Bloch lattice expects exactly one enabled Linear force term.");
                return false;
            }
            linear = &object;
            continue;
        }
 
        fail(reason, "Bloch lattice accepts only negative Rectangle wells plus one Linear X term.");
        return false;
    }
 
    if (wells.size() < 5u) {
        fail(reason, "Bloch lattice needs at least five enabled periodic wells.");
        return false;
    }
    if (linear == nullptr || std::fabs(static_cast<double>(linear->slopeX)) <= 1.0e-12 ||
        std::fabs(static_cast<double>(linear->slopeY)) > 1.0e-9) {
        fail(reason, "Bloch lattice needs one non-zero Linear slopeX and zero slopeY.");
        return false;
    }
 
    std::sort(wells.begin(), wells.end(),
              [](const PotentialObject2D *a, const PotentialObject2D *b) {
                  return a->centerX < b->centerX;
              });
 
    const PotentialObject2D &first = *wells.front();
    const double spacing = static_cast<double>(wells[1]->centerX - wells[0]->centerX);
    if (!(spacing > 0.0) || !std::isfinite(spacing)) {
        fail(reason, "Bloch lattice spacing is invalid.");
        return false;
    }
 
    for (std::size_t i = 0u; i < wells.size(); ++i) {
        const PotentialObject2D &well = *wells[i];
        if (!near(well.value, first.value, std::fabs(first.value)) ||
            !near(well.width, first.width, first.width) ||
            !near(well.height, first.height, first.height) ||
            !near(well.centerY, first.centerY, std::max(1.0f, std::fabs(first.centerY)))) {
            fail(reason, "Bloch lattice wells must have equal depth, width, height and centerY.");
            return false;
        }
        if (i > 0u) {
            const double localSpacing = static_cast<double>(wells[i]->centerX - wells[i - 1u]->centerX);
            if (!near(localSpacing, spacing, spacing)) {
                fail(reason, "Bloch lattice wells are not equally spaced.");
                return false;
            }
        }
    }
 
    geometry.valid = true;
    geometry.spacing = spacing;
    geometry.wellWidth = first.width;
    geometry.wellDepth = -static_cast<double>(first.value);
    geometry.slopeX = linear->slopeX;
    geometry.forceX = -static_cast<double>(linear->slopeX);
    geometry.firstCenterX = wells.front()->centerX;
    geometry.lastCenterX = wells.back()->centerX;
    geometry.wellCount = static_cast<int>(wells.size());
    return true;
}
 
double BlochOscillations2D::blochPeriod(const BlochLatticeGeometry2D &geometry,
                                        double hbar) noexcept
{
    if (!geometry.valid || !(geometry.spacing > 0.0) || !(hbar > 0.0) ||
        !std::isfinite(hbar) || !std::isfinite(geometry.forceX) ||
        std::fabs(geometry.forceX) <= 1.0e-15) {
        return 0.0;
    }
    return 2.0 * kPi * hbar /
           (std::fabs(geometry.forceX) * geometry.spacing);
}
 
double BlochOscillations2D::zoneEdgeMomentum(const BlochLatticeGeometry2D &geometry,
                                             double hbar) noexcept
{
    if (!geometry.valid || !(geometry.spacing > 0.0) || !(hbar > 0.0) ||
        !std::isfinite(hbar)) {
        return 0.0;
    }
    return kPi * hbar / geometry.spacing;
}
 
double BlochOscillations2D::reducedZoneMomentum(const BlochLatticeGeometry2D &geometry,
                                                double initialMomentum,
                                                double time,
                                                double hbar) noexcept
{
    const double edge = zoneEdgeMomentum(geometry, hbar);
    if (!(edge > 0.0) || !std::isfinite(initialMomentum) || !std::isfinite(time)) {
        return 0.0;
    }
 
    // dp/dt = F. One reciprocal-lattice period in physical momentum is 2*edge.
    const double raw = initialMomentum + geometry.forceX * time;
    const double width = 2.0 * edge;
    double wrapped = std::fmod(raw + edge, width);
    if (wrapped < 0.0) wrapped += width;
    return wrapped - edge;
}
 
} // namespace qpe