Jungfraujoch/tests/SpindleCuspLossTest.cpp

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>

#include "../rugnux/SpindleCuspLoss.h"

using Catch::Matchers::WithinAbs;

namespace {
    // lambda = 1 A against d_min = 1 / (2 sin 15 deg) puts theta_max at exactly 15 deg.
    constexpr double WVL = 1.0;
    const double D_MIN = 1.0 / (2.0 * std::sin(15.0 * M_PI / 180.0));

    Coord Perpendicular(const Coord &v) {
        const Coord seed = std::fabs(v.Normalize().x) < 0.9f ? Coord(1, 0, 0) : Coord(0, 1, 0);
        return (v % seed).Normalize();
    }
}

TEST_CASE("SpindleCuspLoss_Orbit", "[Rugnux][Spindle]") {
    SECTION("P1 repairs nothing, anywhere") {
        const CrystalLattice latt(UnitCell{40, 50, 60, 83, 95, 102});
        const auto r = SpindleUnrepairedFraction(*gemmi::find_spacegroup_by_name("P 1"), latt,
                                                 Coord(0.3f, -0.5f, 0.8f), WVL, D_MIN);
        REQUIRE(r.has_value());
        CHECK_THAT(r->cone_fraction, WithinAbs(1.0, 1e-9));
        // Radially weighted double-cone content at theta_max = 15 deg (verified by Monte Carlo).
        CHECK_THAT(r->lost_unique_fraction, WithinAbs(0.0203, 0.0005));
        CHECK_THAT(r->theta_max_deg, WithinAbs(15.0, 1e-6));
    }

    SECTION("a lone diad on the spindle, and one perpendicular to it, both lose the whole cone") {
        const CrystalLattice latt(UnitCell{40, 50, 60, 90, 100, 90});   // P2, unique axis b
        const auto &sg = *gemmi::find_spacegroup_by_name("P 2");
        const Coord diad = latt.Vec1();   // b is perpendicular to a and c, so it IS the axis
        const auto on = SpindleUnrepairedFraction(sg, latt, diad, WVL, D_MIN);
        REQUIRE(on.has_value());
        CHECK_THAT(on->cone_fraction, WithinAbs(1.0, 1e-4));
        const auto perp = SpindleUnrepairedFraction(sg, latt, Perpendicular(diad), WVL, D_MIN);
        REQUIRE(perp.has_value());
        CHECK_THAT(perp->cone_fraction, WithinAbs(1.0, 1e-4));
    }

    SECTION("the same diad at 60 deg from the spindle loses nothing") {
        const CrystalLattice latt(UnitCell{40, 50, 60, 90, 100, 90});
        const Coord diad = latt.Vec1().Normalize();
        const Coord spindle = diad * 0.5f + Perpendicular(diad) * static_cast<float>(std::sqrt(0.75));
        const auto r = SpindleUnrepairedFraction(*gemmi::find_spacegroup_by_name("P 2"), latt,
                                                 spindle, WVL, D_MIN);
        REQUIRE(r.has_value());
        CHECK_THAT(r->lost_unique_fraction, WithinAbs(0.0, 1e-9));
    }

    SECTION("an axis of order >= 3 perpendicular to the spindle repairs the cone completely") {
        // The situation the per-image worst-case bound must flag but the known group clears.
        const CrystalLattice latt(UnitCell{50, 50, 60, 90, 90, 120});
        const Coord three_fold = latt.Vec2();   // c, the 3-fold of P3
        const auto r = SpindleUnrepairedFraction(*gemmi::find_spacegroup_by_name("P 3"), latt,
                                                 Perpendicular(three_fold), WVL, D_MIN);
        REQUIRE(r.has_value());
        CHECK_THAT(r->lost_unique_fraction, WithinAbs(0.0, 1e-9));
    }

    SECTION("a dihedral group repairs an in-plane diad the warning heuristic used to flag") {
        // 622 with an in-plane 2-fold exactly on the spindle: the old any-axis-within-15-deg rule
        // warned here, but the principal 6-fold maps the cone off itself and nothing is lost.
        const CrystalLattice latt(UnitCell{50, 50, 60, 90, 90, 120});
        const auto r = SpindleUnrepairedFraction(*gemmi::find_spacegroup_by_name("P 6 2 2"), latt,
                                                 latt.Vec0(), WVL, D_MIN);
        REQUIRE(r.has_value());
        CHECK_THAT(r->lost_unique_fraction, WithinAbs(0.0, 1e-9));
    }

    SECTION("a cubic group is never severe, whatever the mounting") {
        const CrystalLattice latt(UnitCell{60, 60, 60, 90, 90, 90});
        const auto &sg = *gemmi::find_spacegroup_by_name("P 4 3 2");
        for (const auto &spindle : {Coord(1, 0, 0), Coord(0.3f, -0.5f, 0.8f), Coord(1, 1, 1)}) {
            const auto r = SpindleUnrepairedFraction(sg, latt, spindle, WVL, D_MIN);
            REQUIRE(r.has_value());
            CHECK_THAT(r->lost_unique_fraction, WithinAbs(0.0, 1e-9));
        }
    }

    SECTION("no spindle, no wavelength or no resolution gives no answer") {
        const CrystalLattice latt(UnitCell{40, 50, 60, 90, 90, 90});
        const auto &sg = *gemmi::find_spacegroup_by_name("P 1");
        CHECK_FALSE(SpindleUnrepairedFraction(sg, latt, Coord(0, 0, 0), WVL, D_MIN).has_value());
        CHECK_FALSE(SpindleUnrepairedFraction(sg, latt, Coord(0, 0, 1), 0.0, D_MIN).has_value());
        CHECK_FALSE(SpindleUnrepairedFraction(sg, latt, Coord(0, 0, 1), WVL, 0.0).has_value());
    }
}