leonarski_fandClaude Opus 5.5
0c1c963379
A pure move of three method files out of rugnux/: SpotWidth (spot width and bandwidth, which set the integration radius) joins JFJochBraggIntegration, which now links JFJochImagePreprocessing for the ImagePreprocessorBuffer it reads; SpindleCuspLoss and LeftoverLattices join JFJochIndexing. Only include paths and CMake lists change. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi
94 lines
4.9 KiB
C++
94 lines
4.9 KiB
C++
// 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 "../image_analysis/indexing/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());
|
|
}
|
|
}
|