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

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

#define DEG_TO_RAD static_cast<float>(M_PI/180.0)

CrystalLattice::CrystalLattice(const UnitCell &cell) {
    vec[0] = {cell.a, 0, 0};
    vec[1] = {cell.b * cosf(cell.gamma * DEG_TO_RAD), cell.b * sinf(cell.gamma * DEG_TO_RAD), 0};
    float cx = cell.c * cosf(cell.beta * DEG_TO_RAD);
    float cy = cell.c
                * (cosf(cell.alpha * DEG_TO_RAD) - cosf(cell.beta * DEG_TO_RAD) * cosf(cell.gamma * DEG_TO_RAD))
                / sinf(cell.gamma * DEG_TO_RAD);
    vec[2] = {cx, cy, sqrtf(cell.c*cell.c-cx*cx-cy*cy)};
}

Coord &CrystalLattice::Vec0() {
    return vec[0];
}

Coord &CrystalLattice::Vec1() {
    return vec[1];
}

Coord &CrystalLattice::Vec2() {
    return vec[2];
}

const Coord &CrystalLattice::Vec0() const {
    return vec[0];
}

const Coord &CrystalLattice::Vec1() const {
    return vec[1];
}

const Coord &CrystalLattice::Vec2() const {
    return vec[2];
}

UnitCell CrystalLattice::GetUnitCell() const {
    UnitCell cell{};
    cell.a = vec[0].Length();
    cell.b = vec[1].Length();
    cell.c = vec[2].Length();
    cell.alpha = angle_deg(vec[1], vec[2]);
    cell.beta = angle_deg(vec[0], vec[2]);
    cell.gamma = angle_deg(vec[0], vec[1]);
    return cell;
}

std::vector<float> CrystalLattice::GetVector() const {
    std::vector<float> output(9);
    for (int i = 0; i < 3; i++) {
        output[3 * i + 0] = vec[i].x;
        output[3 * i + 1] = vec[i].y;
        output[3 * i + 2] = vec[i].z;
    }
    return output;
}