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

#include <cmath>
#include <algorithm>
#include "ResolutionShells.h"
#include "JFJochException.h"

ResolutionShells::ResolutionShells(float d_min, float d_max, int32_t nshells)
    : d_min(d_min),
      d_max(d_max),
      one_over_dmin2(1 / (d_min * d_min)),
      one_over_dmax2(1 / (d_max * d_max)),
      nshells(nshells) {
    if (d_min >= d_max || d_min < 0)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Wrong resolution range");

    if (nshells <= 0)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "Number of shells cannot be zero or negative");
}

std::optional<int32_t> ResolutionShells::GetShell(float d) const {
    if (d <= d_min || d > d_max)
        return {};

    if (d == d_max)
        return 0;

    const float one_over_d2 = 1 / (d * d);
    const float shell_fp = (one_over_d2 - one_over_dmax2) / (one_over_dmin2 - one_over_dmax2) * static_cast<float>(nshells);
    return std::clamp<int32_t>(static_cast<int32_t>(shell_fp), 0, nshells - 1);
}

std::vector<float> ResolutionShells::GetShellMeanOneOverResSq() const {
    std::vector<float> ret;
    const float x = (one_over_dmin2 - one_over_dmax2) / static_cast<float>(nshells);
    for (int i = 0; i < nshells; i++) {
        const float one_over_d2 = one_over_dmax2 + (static_cast<float>(i) + 0.5f) * x;
        ret.push_back(one_over_d2);
    }
    return ret;
}

std::vector<float> ResolutionShells::GetShellMinRes() const {
    std::vector<float> ret;
    const float x = (one_over_dmin2 - one_over_dmax2) / static_cast<float>(nshells);
    for (int i = 1; i < nshells; i++) {
        const float one_over_d2 = one_over_dmax2 + static_cast<float>(i) * x;
        const float d = 1 / std::sqrt(one_over_d2);
        ret.push_back(d);
    }
    ret.push_back(d_min);
    return ret;
}