// Copyright (2019-2022) Paul Scherrer Institute
// SPDX-License-Identifier: GPL-3.0-or-later

#include <cmath>
#include "DiffractionGeometry.h"

DiffractionGeometry::DiffractionGeometry(float in_beam_x_pxl, float in_beam_y_pxl, float in_detector_distance_mm,
                                         float in_pixel_size_mm, float in_wavelength_A, Coord in_scattering_vector)
        : beam_x_pxl(in_beam_x_pxl),
        beam_y_pxl(in_beam_y_pxl),
        detector_distance_mm(in_detector_distance_mm),
        pixel_size_mm(in_pixel_size_mm),
        wavelength_A(in_wavelength_A),
        scattering_vector(in_scattering_vector) {}

Coord DiffractionGeometry::DetectorToLab(float x_pxl, float y_pxl) const {
    // Assumes planar detector, 90 deg towards beam
    return {(x_pxl - beam_x_pxl) * pixel_size_mm ,
            (y_pxl - beam_y_pxl) * pixel_size_mm ,
            detector_distance_mm};
}

Coord DiffractionGeometry::DetectorToRecip(float x_pxl, float y_pxl) const {
    return DetectorToLab(x_pxl, y_pxl).Normalize() / wavelength_A - scattering_vector;
}

float DiffractionGeometry::PxlToRes(float detector_x, float detector_y) const {
    auto lab = DetectorToLab(detector_x, detector_y);

    float beam_path = lab.Length();
    if (beam_path == detector_distance_mm) return INFINITY;

    float cos_2theta = detector_distance_mm / beam_path;
    // cos(2theta) = cos(theta)^2 - sin(theta)^2
    // cos(2theta) = 1 - 2*sin(theta)^2
    // Technically two solutions for two theta, but it makes sense only to take positive one in this case
    float sin_theta = sqrt((1-cos_2theta)/2);
    return wavelength_A / (2 * sin_theta);
}

float DiffractionGeometry::ResToPxl(float resolution) const {
    if (resolution == 0)
        return INFINITY;

    float sin_theta = wavelength_A / (2 * resolution);
    float theta = asin(sin_theta);
    float tan_2theta = tan(2 * theta);
    return tan_2theta * detector_distance_mm / pixel_size_mm;
}