Jungfraujoch/image_analysis/WriteMmcif.cpp

// SPDX-FileCopyrightText: 2025 Paul Scherrer Institute
// SPDX-License-Identifier: GPL-3.0-only

#include "WriteMmcif.h"

#include <cmath>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <stdexcept>
#include <ctime>
#include <chrono>

namespace {

/// Current date in ISO-8601 (YYYY-MM-DD) for the _audit block.
std::string CurrentDateISO() {
    auto now = std::chrono::system_clock::now();
    auto t   = std::chrono::system_clock::to_time_t(now);
    std::tm tm{};
#ifdef _WIN32
    gmtime_s(&tm, &t);
#else
    gmtime_r(&t, &tm);
#endif
    char buf[32];
    std::strftime(buf, sizeof(buf), "%Y-%m-%d", &tm);
    return buf;
}

/// Format a double with given decimal places; returns "?" for non-finite.
std::string Fmt(double val, int decimals = 4) {
    if (!std::isfinite(val))
        return "?";
    std::ostringstream ss;
    ss << std::fixed << std::setprecision(decimals) << val;
    return ss.str();
}

/// Quote a CIF string value; returns "?" for empty.
std::string CifStr(const std::string& s) {
    if (s.empty())
        return "?";
    // If it contains spaces or special chars, single-quote it
    if (s.find(' ') != std::string::npos ||
        s.find('\'') != std::string::npos ||
        s.find('#') != std::string::npos)
        return "'" + s + "'";
    return s;
}

} // namespace

void WriteMmcifReflections(std::ostream& out,
                           const std::vector<FrenchWilsonReflection>& reflections,
                           const MmcifMetadata& meta) {

    out << std::fixed;

    // ---------- data block ----------
    out << "data_" << meta.data_block_name << "\n";
    out << "#\n";

    // ---------- _audit ----------
    out << "_audit.revision_id     1\n";
    out << "_audit.creation_date   " << CurrentDateISO() << "\n";
    out << "_audit.update_record   'Initial release'\n";
    out << "#\n";

    // ---------- _software ----------
    out << "_software.name         'Jungfraujoch'\n";
    if (meta.software_version.has_value())
        out << "_software.version      " << CifStr(meta.software_version.value()) << "\n";
    out << "_software.classification reduction\n";
    out << "#\n";

    // ---------- _cell ----------
    out << "_cell.length_a     " << Fmt(meta.unit_cell.a, 3) << "\n";
    out << "_cell.length_b     " << Fmt(meta.unit_cell.b, 3) << "\n";
    out << "_cell.length_c     " << Fmt(meta.unit_cell.c, 3) << "\n";
    out << "_cell.angle_alpha  " << Fmt(meta.unit_cell.alpha, 2) << "\n";
    out << "_cell.angle_beta   " << Fmt(meta.unit_cell.beta, 2) << "\n";
    out << "_cell.angle_gamma  " << Fmt(meta.unit_cell.gamma, 2) << "\n";

    // ---------- _symmetry ----------
    out << "_symmetry.space_group_name_H-M  " << CifStr(meta.space_group_name) << "\n";
    out << "_symmetry.Int_Tables_number      " << meta.space_group_number << "\n";
    out << "#\n";

    // ---------- _diffrn_source / _diffrn_detector ----------
    if (meta.source.has_value()) {
        out << "_diffrn_source.pdbx_synchrotron_site  "
            << CifStr(meta.source.value()) << "\n";
    }

    if (meta.beamline.has_value()) {
        out << "_diffrn_source.pdbx_synchrotron_beamline  "
            << CifStr(meta.beamline.value()) << "\n";
    }

    if (meta.wavelength_A.has_value()) {
        out << "_diffrn_radiation_wavelength.wavelength  "
            << Fmt(meta.wavelength_A.value(), 5) << "\n";
    }

    if (!meta.detector_name.empty()) {
        out << "_diffrn_detector.detector  " << CifStr(meta.detector_name) << "\n";
    }

    if (meta.detector_distance_mm.has_value()) {
        out << "_diffrn_detector.distance  "
            << Fmt(meta.detector_distance_mm.value(), 2) << "\n";
    }

    if (meta.sample_temperature_K.has_value()) {
        out << "_diffrn.ambient_temp  "
            << Fmt(meta.sample_temperature_K.value(), 1) << "\n";
    }

    if (meta.sample_name.has_value() && !meta.sample_name->empty()) {
        out << "_entity.id    1\n";
        out << "_entity.type  polymer\n";
        out << "_entity.pdbx_description  " << CifStr(meta.sample_name.value()) << "\n";
    }

    // ---------- _refln loop ----------
    out << "loop_\n";
    out << "_refln.index_h\n";
    out << "_refln.index_k\n";
    out << "_refln.index_l\n";
    out << "_refln.F_meas_au\n";
    out << "_refln.F_meas_sigma_au\n";
    out << "_refln.intensity_meas\n";
    out << "_refln.intensity_sigma\n";
    out << "_refln.status\n";

    for (const auto& r : reflections) {
        out << std::setw(5) << r.h << " "
            << std::setw(5) << r.k << " "
            << std::setw(5) << r.l << " "
            << std::setw(12) << Fmt(r.F, 4) << " "
            << std::setw(12) << Fmt(r.sigmaF, 4) << " "
            << std::setw(14) << Fmt(r.I, 4) << " "
            << std::setw(14) << Fmt(r.sigmaI, 4) << " "
            << "o"  // 'o' = observed
            << "\n";
    }

    out << "#\n";
    out << "# End of reflections\n";
}

void WriteMmcifReflections(const std::string& path,
                           const std::vector<FrenchWilsonReflection>& reflections,
                           const MmcifMetadata& meta) {
    std::ofstream file(path);
    if (!file)
        throw std::runtime_error("WriteMmcifReflections: cannot open " + path);
    WriteMmcifReflections(file, reflections, meta);
    file.close();
    if (!file)
        throw std::runtime_error("WriteMmcifReflections: I/O error writing " + path);
}