Jungfraujoch/tests/BrukerReaderTest.cpp

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

// The two forms Bruker frames reach us in: a full imgCIF file (what APEX writes when it converts a
// frame, CBFlib's own form - geometry in CIF tables, no "# " lines) and the native .sfrm. Every
// fixture is synthetic and written byte for byte here.

#include <catch2/catch_all.hpp>

#include <cstring>
#include <fstream>
#include <iomanip>
#include <sstream>

#include "../common/JFJochMath.h"
#include "../reader/JFJochCBFReader.h"
#include "../reader/JFJochSFRMReader.h"
#include "../reader/MiniCBF.h"
#include "../reader/SFRM.h"

namespace {

// x-CBF_BYTE_OFFSET for small images: int8 deltas, escaping to int16.
std::string ByteOffset(const std::vector<int32_t> &pixels) {
    std::string out;
    int32_t prev = 0;
    for (const int32_t v: pixels) {
        const int32_t d = v - prev;
        if (d > -128 && d < 128) {
            out.push_back(static_cast<char>(static_cast<int8_t>(d)));
        } else {
            out.push_back(static_cast<char>(-128));
            const auto d16 = static_cast<int16_t>(d);
            out.append(reinterpret_cast<const char *>(&d16), 2);
        }
        prev = v;
    }
    return out;
}

struct ImgCIFFrame {
    double omega = 30, kappa = 60, phi = 10;
    double omega_inc = 0, kappa_inc = 0, phi_inc = 0.5;
    double two_theta = 0;
    double distance = 50;
    std::string polarization;    // extra _diffrn_radiation_wavelength lines, if any
};

// A full imgCIF frame of a kappa head scanning phi, with a detector on a 2theta arm. The image's
// fast index runs along laboratory +Y and its slow one along +X, both increasing, with the image
// centred on the arm's axis: pixel (cols-1)/2, (rows-1)/2 is where the unswung beam meets it.
void WriteImgCIF(const std::string &name, const ImgCIFFrame &f, int64_t nx, int64_t ny,
                 const std::vector<int32_t> &pixels) {
    const double pixel = 0.12;
    std::ostringstream h;
    h << "###CBF: VERSION 1.7.11\n# CBF file written by CBFlib v0.9.6\n\ndata_image\n\n"
      << "_diffrn_scan_frame.integration_time 5\n"
      << "_diffrn_radiation_wavelength.id SOURCE\n_diffrn_radiation_wavelength.wavelength 0.71073\n"
      << f.polarization << "\n"
      << "loop_\n_axis.id\n_axis.depends_on\n_axis.equipment\n_axis.type\n"
      << "_axis.vector[1]\n_axis.vector[2]\n_axis.vector[3]\n_axis.offset[1]\n_axis.offset[2]\n_axis.offset[3]\n"
      << " OMEGA . goniometer rotation 1 0 0 0 0 0\n"
      << " KAPPA OMEGA goniometer rotation 0.642788 0 -0.766044 0 0 0\n"
      << " PHI KAPPA goniometer rotation 1 0 0 0 0 0\n"
      << " TWOTHETA . detector rotation 1 0 0 0 0 0\n"
      << " DX TWOTHETA detector translation 0 0 -1 0 0 0\n"
      << " ELEMENT_X DX detector translation 0 1 0 " << -ny * pixel / 2 << " " << -nx * pixel / 2 << " 0\n"
      << " ELEMENT_Y ELEMENT_X detector translation 1 0 0 0 0 0\n\n"
      << "loop_\n_diffrn_scan_frame_axis.axis_id\n_diffrn_scan_frame_axis.displacement\n"
      << "_diffrn_scan_frame_axis.angle\n"
      << " DX " << f.distance << " ?\n TWOTHETA ? " << f.two_theta << "\n OMEGA ? " << f.omega
      << "\n PHI ? " << f.phi << "\n KAPPA ? " << f.kappa << "\n\n"
      << "loop_\n_diffrn_scan_axis.axis_id\n_diffrn_scan_axis.displacement_start\n"
      << "_diffrn_scan_axis.displacement_increment\n_diffrn_scan_axis.displacement_range\n"
      << "_diffrn_scan_axis.angle_start\n_diffrn_scan_axis.angle_increment\n_diffrn_scan_axis.angle_range\n"
      << " DX " << f.distance << " 0 0 ? ? ?\n TWOTHETA ? ? ? " << f.two_theta << " 0 0\n"
      << " OMEGA ? ? ? " << f.omega << " " << f.omega_inc << " " << f.omega_inc << "\n"
      << " PHI ? ? ? " << f.phi << " " << f.phi_inc << " " << f.phi_inc << "\n"
      << " KAPPA ? ? ? " << f.kappa << " " << f.kappa_inc << " " << f.kappa_inc << "\n\n"
      << "_diffrn_measurement_axis.measurement_id GONIOMETER\n_diffrn_measurement_axis.axis_id PHI\n\n"
      << "_diffrn_detector.id DETECTOR\n_diffrn_detector.detector CCD\n_diffrn_detector.type APEX2\n\n"
      << "loop_\n_array_structure_list_axis.axis_set_id\n_array_structure_list_axis.axis_id\n"
      << "_array_structure_list_axis.displacement\n_array_structure_list_axis.displacement_increment\n"
      << " ELEMENT_X ELEMENT_X " << pixel / 2 << " " << pixel << "\n"
      << " ELEMENT_Y ELEMENT_Y " << pixel / 2 << " " << pixel << "\n\n"
      << "loop_\n_array_structure_list.array_id\n_array_structure_list.index\n_array_structure_list.dimension\n"
      << "_array_structure_list.precedence\n_array_structure_list.direction\n_array_structure_list.axis_set_id\n"
      << " IMAGE 1 " << nx << " 1 increasing ELEMENT_X\n IMAGE 2 " << ny << " 2 increasing ELEMENT_Y\n\n"
      << "_array_intensities.array_id IMAGE\n_array_intensities.gain 15.7\n_array_intensities.overload 600000\n\n"
      << "_array_data.array_id IMAGE\n_array_data.binary_id 1\n_array_data.data\n;\n"
      << "--CIF-BINARY-FORMAT-SECTION--\nContent-Type: application/octet-stream;\n"
      << "     conversions=\"x-CBF_BYTE_OFFSET\"\nContent-Transfer-Encoding: BINARY\n"
      << "X-Binary-Element-Type: \"signed 32-bit integer\"\n"
      << "X-Binary-Number-of-Elements: " << nx * ny << "\n"
      << "X-Binary-Size-Fastest-Dimension: " << nx << "\n"
      << "X-Binary-Size-Second-Dimension: " << ny << "\n\n";
    std::ofstream out(name, std::ios::binary);
    const std::string text = h.str();
    out.write(text.data(), static_cast<std::streamsize>(text.size()));
    out.write(reinterpret_cast<const char *>(minicbf::BINARY_SEPARATOR), sizeof(minicbf::BINARY_SEPARATOR));
    const std::string data = ByteOffset(pixels);
    out.write(data.data(), static_cast<std::streamsize>(data.size()));
}

std::string FrameName(const std::string &prefix, int frame, const std::string &ext) {
    std::ostringstream name;
    name << prefix << "_" << std::setfill('0') << std::setw(4) << frame << ext;
    return name.str();
}

Coord Internal(double x, double y, double z) {
    return {static_cast<float>(x), static_cast<float>(-y), static_cast<float>(-z)};
}

}  // namespace

TEST_CASE("JFJochCBFReader_FullImgCIF_KappaChain", "[HDF5][Full][portable]") {
    const int64_t nx = 24, ny = 16;
    std::vector<int32_t> pixels(nx * ny, 5);
    pixels[3] = -32;            // outside the active area
    pixels[17] = 1000;          // escapes to int16

    SECTION("phi scan on a kappa head, square-on detector") {
        for (int frame = 1; frame <= 3; frame++) {
            ImgCIFFrame f;
            f.phi = 10 + 0.5 * (frame - 1);
            WriteImgCIF(FrameName("imgcif_kappa", frame, ".cbf"), f, nx, ny, pixels);
        }
        REQUIRE(JFJochCBFReader::CanRead("imgcif_kappa_0002.cbf"));
        JFJochCBFReader reader;
        REQUIRE_NOTHROW(reader.ReadFiles("imgcif_kappa_0002.cbf"));
        CHECK(reader.GetNumberOfImages() == 3);
        CHECK(reader.CountsInADU());
        const auto x = reader.GetDataset()->experiment;

        // phi turns about its own axis as omega (30) and kappa (60) hold it; GoniometerAxis takes
        // minus the right-handed direction
        const Coord kappa_axis = Internal(0.642788, 0, -0.766044);
        const Coord expected = -(RotMatrix(static_cast<float>(30 * PI / 180), Coord(1, 0, 0))
                                 * (RotMatrix(static_cast<float>(60 * PI / 180), kappa_axis) * Coord(1, 0, 0)));
        REQUIRE(x.GetGoniometer().has_value());
        CHECK((x.GetGoniometer()->GetAxis() - expected).Length() < 1e-5f);
        CHECK(x.GetGoniometer()->GetStart_deg() == Catch::Approx(10.0));
        CHECK(x.GetGoniometer()->GetIncrement_deg() == Catch::Approx(0.5));
        CHECK(x.GetGoniometer()->GetName() == "phi");

        CHECK(x.GetDetectorDistance_mm() == Catch::Approx(50.0));
        CHECK(x.GetBeamX_pxl() == Catch::Approx((nx - 1) / 2.0));
        CHECK(x.GetBeamY_pxl() == Catch::Approx((ny - 1) / 2.0));
        CHECK(x.GetPixelSize_mm() == Catch::Approx(0.12));
        CHECK(x.GetWavelength_A() == Catch::Approx(0.71073));
        CHECK(x.GetSaturationLimit() == 600000);
        CHECK_FALSE(x.GetPolarizationFactor().has_value());   // the header states none
        // fast along laboratory +Y (internal -y), slow along +X
        const auto geom = x.GetDiffractionGeometry();
        CHECK((geom.GetFastAxis() - Coord(0, -1, 0)).Length() < 1e-5f);
        CHECK((geom.GetSlowAxis() - Coord(1, 0, 0)).Length() < 1e-5f);
        CHECK((geom.GetNormalAxis() - Coord(0, 0, 1)).Length() < 1e-5f);

        // ReadRawImage hands the byte-offset stream on undecoded; LoadImage decodes it
        const auto image = reader.LoadImage(1);
        REQUIRE(image);
        CHECK(image->Image()[0] == 5);
        CHECK(image->Image()[17] == 1000);
        CHECK(reader.GetDataset()->pixel_mask->GetMask()[3] != 0);
        CHECK(reader.GetDataset()->pixel_mask->GetMask()[0] == 0);
        reader.Close();
    }

    SECTION("a detector swung out on the 2theta arm") {
        for (int frame = 1; frame <= 2; frame++) {
            ImgCIFFrame f;
            f.phi = 10 + 0.5 * (frame - 1);
            f.two_theta = -20;
            WriteImgCIF(FrameName("imgcif_swing", frame, ".cbf"), f, nx, ny, pixels);
        }
        JFJochCBFReader reader;
        REQUIRE_NOTHROW(reader.ReadFiles("imgcif_swing_0001.cbf"));
        const auto x = reader.GetDataset()->experiment;
        // The arm carries the square-on geometry with it: distance and foot of the normal stay.
        CHECK(x.GetDetectorDistance_mm() == Catch::Approx(50.0));
        CHECK(x.GetBeamX_pxl() == Catch::Approx((nx - 1) / 2.0).margin(1e-3));
        CHECK(x.GetBeamY_pxl() == Catch::Approx((ny - 1) / 2.0).margin(1e-3));
        const Coord normal = RotMatrix(static_cast<float>(-20 * PI / 180), Coord(1, 0, 0)) * Coord(0, 0, 1);
        CHECK((x.GetDiffractionGeometry().GetNormalAxis() - normal).Length() < 1e-5f);
        reader.Close();
    }

    for (int frame = 1; frame <= 3; frame++) {
        remove(FrameName("imgcif_kappa", frame, ".cbf").c_str());
        remove(FrameName("imgcif_swing", frame, ".cbf").c_str());
    }
}

// What the full imgCIF path must never do: read a series it cannot describe as one rotation about
// one axis as if it could.
TEST_CASE("JFJochCBFReader_FullImgCIF_Polarization", "[HDF5][Full][portable]") {
    const int64_t nx = 24, ny = 16;
    const std::vector<int32_t> pixels(nx * ny, 5);
    auto read = [&](const std::string &prefix, const std::string &lines) {
        for (int frame = 1; frame <= 2; frame++) {
            ImgCIFFrame f;
            f.phi = 10 + 0.5 * (frame - 1);
            f.polarization = lines;
            WriteImgCIF(FrameName(prefix, frame, ".cbf"), f, nx, ny, pixels);
        }
        JFJochCBFReader reader;
        reader.ReadFiles(FrameName(prefix, 1, ".cbf"));
        return reader.GetDataset()->experiment.GetPolarizationFactor();
    };
    // the degree as stated, in the horizontal plane when the plane's normal is the laboratory Y axis
    const auto horizontal = read("imgcif_pol_h", "_diffrn_radiation_wavelength.polarizn_source_ratio 0.0444\n"
                                                 "_diffrn_radiation_wavelength.polarizn_source_norm 0\n");
    REQUIRE(horizontal.has_value());
    CHECK(*horizontal == Catch::Approx(0.0444));
    // a normal turned by 90 degrees puts the plane vertical, which the correction takes as the opposite sign
    const auto vertical = read("imgcif_pol_v", "_diffrn_radiation_wavelength.polarizn_source_ratio 0.5\n"
                                               "_diffrn_radiation_wavelength.polarizn_source_norm 90\n");
    REQUIRE(vertical.has_value());
    CHECK(*vertical == Catch::Approx(-0.5));
}

TEST_CASE("JFJochCBFReader_FullImgCIF_refuses_what_it_cannot_describe", "[HDF5][Full][portable]") {
    const int64_t nx = 8, ny = 8;
    const std::vector<int32_t> pixels(nx * ny, 1);

    SECTION("two axes moving at once") {
        for (int frame = 1; frame <= 2; frame++) {
            ImgCIFFrame f;
            f.omega_inc = 0.5;
            WriteImgCIF(FrameName("imgcif_two", frame, ".cbf"), f, nx, ny, pixels);
        }
        CHECK(JFJochCBFReader::CanRead("imgcif_two_0001.cbf"));
        JFJochCBFReader reader;
        CHECK_THROWS_AS(reader.ReadFiles("imgcif_two_0001.cbf"), JFJochException);
    }

    SECTION("kappa moved between frames of one name template") {
        for (int frame = 1; frame <= 2; frame++) {
            ImgCIFFrame f;
            f.phi = 10 + 0.5 * (frame - 1);
            f.kappa = (frame == 1) ? 60 : 0;
            WriteImgCIF(FrameName("imgcif_reset", frame, ".cbf"), f, nx, ny, pixels);
        }
        JFJochCBFReader reader;
        CHECK_THROWS_AS(reader.ReadFiles("imgcif_reset_0001.cbf"), JFJochException);
    }

    SECTION("nothing turns") {
        for (int frame = 1; frame <= 2; frame++) {
            ImgCIFFrame f;
            f.phi_inc = 0;
            WriteImgCIF(FrameName("imgcif_still", frame, ".cbf"), f, nx, ny, pixels);
        }
        JFJochCBFReader reader;
        CHECK_THROWS_AS(reader.ReadFiles("imgcif_still_0001.cbf"), JFJochException);
    }

    for (int frame = 1; frame <= 2; frame++)
        for (const char *p: {"imgcif_two", "imgcif_reset", "imgcif_still"})
            remove(FrameName(p, frame, ".cbf").c_str());
}

TEST_CASE("ImgCIF_Item_and_Loop", "[HDF5][Full][portable]") {
    const std::string text = "\n_a.x 1\nloop_\n_b.id\n_b.v\n p 1\n q 2\n\n_c.id 'quoted value'\n";
    CHECK(imgcif::Item(text, "_a.x").value_or("") == "1");
    CHECK(imgcif::Item(text, "_c.id").value_or("") == "quoted value");
    CHECK_FALSE(imgcif::Item(text, "_b.id").has_value());     // heads a loop
    const auto rows = imgcif::ParseLoop(text, "_b.id");
    REQUIRE(rows.size() == 2);
    CHECK(rows[1].at("_b.v") == "2");
    // a tag outside any loop is not the rows of the loop before it
    CHECK(imgcif::ParseLoop(text, "_c.id").empty());
}

namespace {

// One FORMAT 100 frame: 1-byte pixels, the 16- and 32-bit overflow tables, baseline 32.
struct SFRMFrame {
    int64_t axis = 2;
    double start = 0, increment = 0.5;
    std::array<double, 4> angles{0, 0, 0, 0};    // 2theta omega phi chi
    double distance_cm = 5.0;
    double monochromator_2theta = 0;   // FILTER2
    std::vector<uint8_t> pixels;
    std::vector<uint16_t> over16;
    std::vector<int32_t> over32;
};

std::string Line(const std::string &key, const std::string &value) {
    std::string k = key;
    k.resize(7, ' ');
    std::string v = value;
    v.resize(72, ' ');
    return k + ":" + v;
}

void Pad(std::string &s, size_t to) {
    s.resize((s.size() + to - 1) / to * to, '\0');
}

void WriteSFRM(const std::string &name, const SFRMFrame &f, int64_t n) {
    std::ostringstream angles;
    for (double a: f.angles)
        angles << a << " ";
    std::string h;
    h += Line("FORMAT", "100");
    h += Line("VERSION", "11");
    h += Line("HDRBLKS", "15");
    h += Line("NOVERFL", "-1 " + std::to_string(f.over16.size()) + " " + std::to_string(f.over32.size()));
    h += Line("NPIXELB", "1 1");
    h += Line("NROWS", std::to_string(n));
    h += Line("NCOLS", std::to_string(n));
    h += Line("WORDORD", "0");
    h += Line("LONGORD", "0");
    h += Line("START", std::to_string(f.start));
    h += Line("INCREME", std::to_string(f.increment));
    h += Line("ANGLES", angles.str());
    h += Line("AXIS", std::to_string(f.axis));
    h += Line("DISTANC", std::to_string(f.distance_cm));
    h += Line("CENTER", std::to_string(n / 2.0) + " " + std::to_string(n / 2.0));
    h += Line("WAVELEN", "1.54184 1.54056 1.54439");
    h += Line("CUMULAT", "10.0");
    h += Line("DETTYPE", "CCD-TEST 64.0 0.1 0 0.254 0.05 1");
    h += Line("NEXP", "2 0 32 0 1");
    h += Line("CCDPARM", "14 10 170 0 900000");
    h += Line("DETPAR", "0 0 0 0 0 0");
    h += Line("FILTER2", std::to_string(f.monochromator_2theta) + " 0");
    h.resize(15 * 512, ' ');
    std::string body(f.pixels.begin(), f.pixels.end());
    Pad(body, 16);
    std::string t16(reinterpret_cast<const char *>(f.over16.data()), f.over16.size() * 2);
    Pad(t16, 16);
    std::string t32(reinterpret_cast<const char *>(f.over32.data()), f.over32.size() * 4);
    Pad(t32, 16);
    std::ofstream out(name, std::ios::binary);
    for (const auto *s: {&h, &body, &t16, &t32})
        out.write(s->data(), static_cast<std::streamsize>(s->size()));
}

}  // namespace

TEST_CASE("SFRM_Decode_overflow_tables", "[HDF5][Full][portable]") {
    const int64_t n = 8;
    SFRMFrame f;
    f.pixels.assign(n * n, 40);         // 8 counts above the baseline of 32
    f.pixels[0] = 0;                    // outside the active area: below the baseline
    f.pixels[5] = 255;                  // -> 16-bit table: 1000
    f.pixels[9] = 255;                  // -> 16-bit 65535 -> 32-bit 100000
    f.pixels[20] = 255;                 // -> 16-bit 65535 -> 32-bit -14 (signed)
    f.pixels[30] = 255;                 // -> 16-bit 300
    f.over16 = {1000, 65535, 65535, 300};
    f.over32 = {100000, -14};
    WriteSFRM("sfrm_tables_0001.sfrm", f, n);

    REQUIRE(sfrm::CanRead("sfrm_tables_0001.sfrm"));
    std::vector<int32_t> v;
    const auto h = sfrm::Read("sfrm_tables_0001.sfrm", v);
    CHECK(h.format == 100);
    CHECK(h.baseline == 32);
    REQUIRE(v.size() == static_cast<size_t>(n * n));
    CHECK(v[1] == 8);
    CHECK(v[0] == -32);
    CHECK(v[5] == 1000 - 32);
    CHECK(v[9] == 100000 - 32);
    CHECK(v[20] == -14 - 32);
    CHECK(v[30] == 300 - 32);
    remove("sfrm_tables_0001.sfrm");

    SECTION("a table shorter than the pixels that point into it is an error") {
        SFRMFrame g = f;
        g.over16 = {1000};
        g.over32 = {};
        WriteSFRM("sfrm_short_0001.sfrm", g, n);
        CHECK_THROWS_AS(sfrm::Read("sfrm_short_0001.sfrm", v), JFJochException);
        remove("sfrm_short_0001.sfrm");
    }
}

TEST_CASE("SFRM_Decode_format86_ascii_overflows", "[HDF5][Full][portable]") {
    const int64_t n = 4;
    std::string h;
    h += Line("FORMAT", "86");
    h += Line("VERSION", "9");
    h += Line("HDRBLKS", "5");
    h += Line("NOVERFL", "2");
    h += Line("NPIXELB", "1");
    h += Line("NROWS", std::to_string(n));
    h += Line("NCOLS", std::to_string(n));
    h.resize(5 * 512, ' ');
    std::string body(n * n, static_cast<char>(7));
    body[2] = static_cast<char>(255);
    body[11] = static_cast<char>(255);
    std::string table = "    70000      2   123456     11";
    Pad(table, 512);
    std::ofstream("sfrm_old_0001.sfrm", std::ios::binary) << h << body << table;

    std::vector<int32_t> v;
    const auto hdr = sfrm::Read("sfrm_old_0001.sfrm", v);
    CHECK(hdr.format == 86);
    CHECK(v[0] == 7);
    CHECK(v[2] == 70000);
    CHECK(v[11] == 123456);
    remove("sfrm_old_0001.sfrm");
}

TEST_CASE("JFJochSFRMReader_Geometry", "[HDF5][Full][portable]") {
    const int64_t n = 8;
    SFRMFrame f;
    f.pixels.assign(n * n, 40);
    f.pixels[0] = 0;

    SECTION("omega scan") {
        for (int frame = 1; frame <= 3; frame++) {
            f.start = 0.5 * (frame - 1);
            f.angles = {0, f.start, 0, 54.74};
            WriteSFRM(FrameName("sfrm_omega", frame, ".sfrm"), f, n);
        }
        REQUIRE(JFJochSFRMReader::CanRead("sfrm_omega_0001.sfrm"));
        JFJochSFRMReader reader;
        REQUIRE_NOTHROW(reader.ReadFiles("sfrm_omega_0001.sfrm"));
        CHECK(reader.GetNumberOfImages() == 3);
        const auto x = reader.GetDataset()->experiment;
        REQUIRE(x.GetGoniometer().has_value());
        // omega is the outermost axis, so a fixed chi does not tilt it
        CHECK((x.GetGoniometer()->GetAxis() - Coord(1, 0, 0)).Length() < 1e-5f);
        CHECK(x.GetGoniometer()->GetIncrement_deg() == Catch::Approx(0.5));
        // DETTYPE: 64 pixels/cm on a 512 frame, so 1/64 cm * 512/8 per pixel; distance with the
        // 0.1 cm face-to-grid offset
        CHECK(x.GetPixelSize_mm() == Catch::Approx(10.0 / 64.0 * 512.0 / n));
        CHECK(x.GetDetectorDistance_mm() == Catch::Approx(51.0));
        CHECK(x.GetBeamX_pxl() == Catch::Approx((n - 1) / 2.0).margin(1e-3));
        CHECK(x.GetBeamY_pxl() == Catch::Approx((n - 1) / 2.0).margin(1e-3));
        CHECK(x.GetSaturationLimit() == 900000);
        CHECK(x.GetWavelength_A() == Catch::Approx(1.54184));
        // laboratory frames with no monochromator: unpolarized, not the synchrotron default
        REQUIRE(x.GetPolarizationFactor().has_value());
        CHECK(*x.GetPolarizationFactor() == Catch::Approx(0.0).margin(1e-6));
        CHECK(reader.GetDataset()->pixel_mask->GetMask()[0] != 0);
        CHECK(reader.GetDataset()->pixel_mask->GetMask()[1] == 0);

        JFJochReaderRawImage raw;
        REQUIRE(reader.ReadRawImage(2, raw));
        CHECK(reinterpret_cast<const int32_t *>(raw.image.GetCompressed())[1] == 8);
        reader.Close();
    }

    SECTION("a sweep that turns backwards is stated forwards about the opposite axis") {
        for (int frame = 1; frame <= 3; frame++) {
            f.start = 10.0 - 0.5 * (frame - 1);
            f.increment = -0.5;
            f.angles = {0, f.start, 0, 54.74};
            f.monochromator_2theta = 12.17;
            WriteSFRM(FrameName("sfrm_back", frame, ".sfrm"), f, n);
        }
        JFJochSFRMReader reader;
        REQUIRE_NOTHROW(reader.ReadFiles("sfrm_back_0001.sfrm"));
        CHECK(reader.GetNumberOfImages() == 3);
        const auto x = reader.GetDataset()->experiment;
        REQUIRE(x.GetGoniometer().has_value());
        CHECK((x.GetGoniometer()->GetAxis() - Coord(-1, 0, 0)).Length() < 1e-5f);
        CHECK(x.GetGoniometer()->GetIncrement_deg() == Catch::Approx(0.5));
        CHECK(x.GetGoniometer()->GetStart_deg() == Catch::Approx(-10.0));
        // a graphite monochromator: the 1 - cos^2(2theta_m) APEX's conversion states, horizontal
        REQUIRE(x.GetPolarizationFactor().has_value());
        const double c = std::cos(12.17 * PI / 180.0);
        CHECK(*x.GetPolarizationFactor() == Catch::Approx(1.0 - c * c));
        reader.Close();
    }

    SECTION("phi scan on a fixed-chi head turns about the inclined phi axis") {
        for (int frame = 1; frame <= 2; frame++) {
            f.axis = 3;
            f.start = 0.5 * (frame - 1);
            f.angles = {0, 0, f.start, 54.74};
            WriteSFRM(FrameName("sfrm_phi", frame, ".sfrm"), f, n);
        }
        JFJochSFRMReader reader;
        REQUIRE_NOTHROW(reader.ReadFiles("sfrm_phi_0001.sfrm"));
        const auto x = reader.GetDataset()->experiment;
        // chi turns about the beam: imgCIF +Z, internal -z
        const Coord expected = -(RotMatrix(static_cast<float>(54.74 * PI / 180), Coord(0, 0, -1)) * Coord(1, 0, 0));
        CHECK((x.GetGoniometer()->GetAxis() - expected).Length() < 1e-5f);
        reader.Close();
    }

    SECTION("chi moved between frames of one run is refused") {
        for (int frame = 1; frame <= 2; frame++) {
            f.axis = 3;
            f.start = 0.5 * (frame - 1);
            f.angles = {0, 0, f.start, frame == 1 ? 54.74 : 0.0};
            WriteSFRM(FrameName("sfrm_reset", frame, ".sfrm"), f, n);
        }
        JFJochSFRMReader reader;
        CHECK_THROWS_AS(reader.ReadFiles("sfrm_reset_0001.sfrm"), JFJochException);
    }

    for (int frame = 1; frame <= 3; frame++)
        for (const char *p: {"sfrm_omega", "sfrm_phi", "sfrm_reset"})
            remove(FrameName(p, frame, ".sfrm").c_str());
}