490 lines
18 KiB
C++
490 lines
18 KiB
C++
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <catch2/catch_all.hpp>
|
|
#include <iostream>
|
|
|
|
#include "../image_analysis/geom_refinement/XtalOptimizer.h"
|
|
#include "../image_analysis/bragg_integration/BraggPrediction.h"
|
|
|
|
TEST_CASE("XtalOptimizer") {
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
CrystalLattice latt_i(40,40,80,90,90,90);
|
|
|
|
auto pred = CalcBraggPredictions(exp_i, latt_i,1.5, 0.001);
|
|
|
|
std::vector<SpotToSave> spots;
|
|
|
|
for (const auto &[ang, refl]: pred) {
|
|
spots.push_back(SpotToSave{
|
|
refl.predicted_x,
|
|
refl.predicted_y
|
|
});
|
|
}
|
|
|
|
XtalOptimizerData xtal_opt;
|
|
xtal_opt.latt = CrystalLattice(40.2,39.4,80.2, 90,91, 89);
|
|
xtal_opt.geom.BeamX_pxl(1010).BeamY_pxl(995).DetectorDistance_mm(200)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Monoclinic;
|
|
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, spots));
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
auto uc_i = latt_i.GetUnitCell();
|
|
auto uc_o = xtal_opt.latt.GetUnitCell();
|
|
|
|
std::cout << "Beam center: " << xtal_opt.geom.GetBeamX_pxl() << " " << xtal_opt.geom.GetBeamY_pxl() << std::endl;
|
|
std::cout << "Unit cell: " << uc_o.a << " " << uc_o.b << " " << uc_o.c << std::endl;
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.05);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.05);
|
|
|
|
CHECK(fabsf(uc_i.a - uc_o.a) < 0.1);
|
|
CHECK(fabsf(uc_i.b - uc_o.b) < 0.1);
|
|
CHECK(fabsf(uc_i.c - uc_o.c) < 0.2);
|
|
CHECK(fabsf(uc_i.alpha - uc_o.alpha) < 0.1);
|
|
CHECK(fabsf(uc_i.beta - uc_o.beta) < 0.1);
|
|
CHECK(fabsf(uc_i.gamma - uc_o.gamma) < 0.1);
|
|
}
|
|
|
|
|
|
TEST_CASE("XtalOptimizer_NoBeamCenter") {
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
CrystalLattice latt_i(40,50,80,90,95,90);
|
|
|
|
auto pred = CalcBraggPredictions(exp_i, latt_i,1.5, 0.001);
|
|
|
|
std::vector<SpotToSave> spots;
|
|
|
|
for (const auto &[ang, refl]: pred) {
|
|
spots.push_back(SpotToSave{
|
|
refl.predicted_x,
|
|
refl.predicted_y
|
|
});
|
|
}
|
|
|
|
XtalOptimizerData xtal_opt;
|
|
xtal_opt.latt = CrystalLattice(40.2,49.4,80.2, 90,94, 89);
|
|
xtal_opt.geom.BeamX_pxl(999.8).BeamY_pxl(1000.2).DetectorDistance_mm(200)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Monoclinic;
|
|
xtal_opt.refine_beam_center = false;
|
|
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, spots));
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
auto uc_i = latt_i.GetUnitCell();
|
|
auto uc_o = xtal_opt.latt.GetUnitCell();
|
|
|
|
std::cout << "Beam center: " << xtal_opt.geom.GetBeamX_pxl() << " " << xtal_opt.geom.GetBeamY_pxl() << std::endl;
|
|
std::cout << "Unit cell: " << uc_o.a << " " << uc_o.b << " " << uc_o.c << std::endl;
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - 999.8) < 0.01);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - 1000.2) < 0.01);
|
|
|
|
CHECK(fabsf(uc_i.a - uc_o.a) < 0.1);
|
|
CHECK(fabsf(uc_i.b - uc_o.b) < 0.1);
|
|
CHECK(fabsf(uc_i.c - uc_o.c) < 0.2);
|
|
CHECK(fabsf(uc_i.alpha - uc_o.alpha) < 0.1);
|
|
CHECK(fabsf(uc_i.beta - uc_o.beta) < 0.1);
|
|
CHECK(fabsf(uc_i.gamma - uc_o.gamma) < 0.1);
|
|
}
|
|
|
|
TEST_CASE("XtalOptimizer_orthorombic") {
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
CrystalLattice latt_i(40,50,80,90,90,90);
|
|
|
|
auto pred = CalcBraggPredictions(exp_i, latt_i,1.5, 0.001);
|
|
|
|
std::vector<SpotToSave> spots;
|
|
|
|
for (const auto &[ang, refl]: pred) {
|
|
spots.push_back(SpotToSave{
|
|
refl.predicted_x,
|
|
refl.predicted_y
|
|
});
|
|
}
|
|
|
|
XtalOptimizerData xtal_opt;
|
|
xtal_opt.latt = CrystalLattice(40.2,49.6,80.3, 90,91, 89);
|
|
xtal_opt.geom.BeamX_pxl(1005).BeamY_pxl(997).DetectorDistance_mm(200)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Orthorhombic;
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, spots));
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
auto uc_i = latt_i.GetUnitCell();
|
|
auto uc_o = xtal_opt.latt.GetUnitCell();
|
|
|
|
std::cout << "Beam center: " << xtal_opt.geom.GetBeamX_pxl() << " " << xtal_opt.geom.GetBeamY_pxl() << std::endl;
|
|
std::cout << "Unit cell: " << uc_o.a << " " << uc_o.b << " " << uc_o.c << std::endl;
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.1);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.1);
|
|
|
|
CHECK(fabsf(uc_i.a - uc_o.a) < 0.1);
|
|
CHECK(fabsf(uc_i.b - uc_o.b) < 0.1);
|
|
CHECK(fabsf(uc_i.c - uc_o.c) < 0.2);
|
|
CHECK(fabs(uc_o.alpha - 90) < 0.02);
|
|
CHECK(fabs(uc_o.beta - 90) < 0.02);
|
|
CHECK(fabs(uc_o.gamma - 90) < 0.02);
|
|
}
|
|
|
|
TEST_CASE("XtalOptimizer_tetragonal") {
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
CrystalLattice latt_i(40,40,80,90,90,90);
|
|
|
|
auto pred = CalcBraggPredictions(exp_i, latt_i,1.5, 0.001);
|
|
|
|
std::vector<SpotToSave> spots;
|
|
|
|
for (const auto &[ang, refl]: pred) {
|
|
spots.push_back(SpotToSave{
|
|
refl.predicted_x,
|
|
refl.predicted_y
|
|
});
|
|
}
|
|
|
|
XtalOptimizerData xtal_opt;
|
|
xtal_opt.latt = CrystalLattice(40.6,39.3,80.5, 90,91, 89);
|
|
xtal_opt.geom.BeamX_pxl(1010).BeamY_pxl(995).DetectorDistance_mm(200)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Tetragonal;
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, spots));
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
auto uc_i = latt_i.GetUnitCell();
|
|
auto uc_o = xtal_opt.latt.GetUnitCell();
|
|
|
|
std::cout << "Beam center: " << xtal_opt.geom.GetBeamX_pxl() << " " << xtal_opt.geom.GetBeamY_pxl() << std::endl;
|
|
std::cout << "Unit cell: " << uc_o.a << " " << uc_o.b << " " << uc_o.c << std::endl;
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.1);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.1);
|
|
|
|
CHECK(fabsf(uc_i.a - uc_o.a) < 0.1);
|
|
CHECK(fabsf(uc_i.b - uc_o.b) < 0.1);
|
|
CHECK(fabsf(uc_i.c - uc_o.c) < 0.2);
|
|
CHECK(fabs(uc_o.alpha - 90) < 0.02);
|
|
CHECK(fabs(uc_o.beta - 90) < 0.02);
|
|
CHECK(fabs(uc_o.gamma - 90) < 0.02);
|
|
}
|
|
|
|
TEST_CASE("XtalOptimizer_hexagonal") {
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
CrystalLattice latt_i(40,40,70,90,90,120);
|
|
|
|
auto pred = CalcBraggPredictions(exp_i, latt_i,1.5, 0.001);
|
|
|
|
std::vector<SpotToSave> spots;
|
|
|
|
for (const auto &[ang, refl]: pred) {
|
|
spots.push_back(SpotToSave{
|
|
refl.predicted_x,
|
|
refl.predicted_y
|
|
});
|
|
}
|
|
|
|
XtalOptimizerData xtal_opt;
|
|
xtal_opt.latt = CrystalLattice(39.5,39.8,70.1, 90,90, 119.5);
|
|
xtal_opt.geom.BeamX_pxl(1007).BeamY_pxl(990).DetectorDistance_mm(200)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Hexagonal;
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, spots));
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
auto uc_i = latt_i.GetUnitCell();
|
|
auto uc_o = xtal_opt.latt.GetUnitCell();
|
|
|
|
std::cout << "Beam center: " << xtal_opt.geom.GetBeamX_pxl() << " " << xtal_opt.geom.GetBeamY_pxl() << std::endl;
|
|
std::cout << "Unit cell: " << uc_o.a << " " << uc_o.b << " " << uc_o.c << " " << uc_o.alpha << " " << uc_o.beta
|
|
<< " " << uc_o.gamma << std::endl;
|
|
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.1);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.1);
|
|
|
|
CHECK(fabsf(uc_i.a - uc_o.a) < 0.1);
|
|
CHECK(fabsf(uc_i.b - uc_o.b) < 0.1);
|
|
CHECK(fabsf(uc_i.c - uc_o.c) < 0.2);
|
|
CHECK(fabs(uc_o.alpha - 90) < 0.02);
|
|
CHECK(fabs(uc_o.beta - 90) < 0.01);
|
|
CHECK(fabs(uc_o.gamma - 120) < 0.01);
|
|
}
|
|
|
|
TEST_CASE("XtalOptimizer_hexagonal_unconstrained") {
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
CrystalLattice latt_i(40,40,70,90,90,120);
|
|
|
|
auto pred = CalcBraggPredictions(exp_i, latt_i,1.5, 0.001);
|
|
|
|
std::vector<SpotToSave> spots;
|
|
|
|
for (const auto &[ang, refl]: pred) {
|
|
spots.push_back(SpotToSave{
|
|
refl.predicted_x,
|
|
refl.predicted_y
|
|
});
|
|
}
|
|
|
|
XtalOptimizerData xtal_opt;
|
|
xtal_opt.latt = CrystalLattice(39.9,39.8,70.1, 90,90, 120);
|
|
xtal_opt.geom.BeamX_pxl(1002).BeamY_pxl(998).DetectorDistance_mm(200)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Triclinic;
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, spots));
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
auto uc_i = latt_i.GetUnitCell();
|
|
auto uc_o = xtal_opt.latt.GetUnitCell();
|
|
|
|
std::cout << "Beam center: " << xtal_opt.geom.GetBeamX_pxl() << " " << xtal_opt.geom.GetBeamY_pxl() << std::endl;
|
|
std::cout << "Unit cell: " << uc_o.a << " " << uc_o.b << " " << uc_o.c << " " << uc_o.alpha << " " << uc_o.beta
|
|
<< " " << uc_o.gamma << std::endl;
|
|
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.1);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.1);
|
|
|
|
CHECK(fabsf(uc_i.a - uc_o.a) < 0.1);
|
|
CHECK(fabsf(uc_i.b - uc_o.b) < 0.1);
|
|
CHECK(fabsf(uc_i.c - uc_o.c) < 0.2);
|
|
CHECK(fabs(uc_o.alpha - 90) < 0.02);
|
|
CHECK(fabs(uc_o.beta - 90) < 0.01);
|
|
CHECK(fabs(uc_o.gamma - 120) < 0.01);
|
|
}
|
|
|
|
TEST_CASE("XtalOptimizer_cubic") {
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
CrystalLattice latt_i(Coord(40,0,0),
|
|
Coord(0, 40 / sqrt(2), -40 / sqrt(2)),
|
|
Coord(0, 40 / sqrt(2), 40 / sqrt(2)));
|
|
auto uc_i = latt_i.GetUnitCell();
|
|
|
|
auto pred = CalcBraggPredictions(exp_i, latt_i,1.5, 0.001);
|
|
|
|
std::vector<SpotToSave> spots;
|
|
|
|
for (const auto &[ang, refl]: pred) {
|
|
spots.push_back(SpotToSave{
|
|
refl.predicted_x,
|
|
refl.predicted_y
|
|
});
|
|
}
|
|
|
|
XtalOptimizerData xtal_opt;
|
|
xtal_opt.latt = CrystalLattice(Coord(39,0,0),
|
|
Coord(0, 39.5 / sqrt(2), -40.5 / sqrt(2)),
|
|
Coord(0, 39.2 / sqrt(2), 39.7 / sqrt(2)));
|
|
xtal_opt.geom.BeamX_pxl(1007).BeamY_pxl(990).DetectorDistance_mm(200)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Cubic;
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, spots));
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
auto uc_o = xtal_opt.latt.GetUnitCell();
|
|
|
|
std::cout << "Beam center: " << xtal_opt.geom.GetBeamX_pxl() << " " << xtal_opt.geom.GetBeamY_pxl() << std::endl;
|
|
std::cout << "Unit cell: " << uc_o.a << " " << uc_o.b << " " << uc_o.c << std::endl;
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.1);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.1);
|
|
|
|
CHECK(fabsf(uc_i.a - uc_o.a) < 0.1);
|
|
CHECK(fabsf(uc_i.b - uc_o.b) < 0.1);
|
|
CHECK(fabsf(uc_i.c - uc_o.c) < 0.2);
|
|
CHECK(fabs(uc_o.alpha - 90) < 0.02);
|
|
CHECK(fabs(uc_o.beta - 90) < 0.02);
|
|
CHECK(fabs(uc_o.gamma - 90) < 0.02);
|
|
}
|
|
|
|
TEST_CASE("XtalOptimizer_cubic_I") {
|
|
// Conventional cubic (take as I-conventional basis): a=b=c=50, orthogonal
|
|
const float a0 = 50.f;
|
|
const float expected_len = a0 * std::sqrt(3.f) * 0.5f; // a*sqrt(3)/2
|
|
const float expected_angle = 109.4712206f; // arccos(-1/3) in degrees
|
|
|
|
CrystalLattice latt_i(expected_len, expected_len, expected_len,
|
|
expected_angle, expected_angle, expected_angle);
|
|
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
auto uc_i = latt_i.GetUnitCell();
|
|
|
|
auto pred = CalcBraggPredictions(exp_i, latt_i, 1.5, 0.001);
|
|
|
|
std::vector<SpotToSave> spots;
|
|
|
|
for (const auto &[ang, refl]: pred) {
|
|
spots.push_back(SpotToSave{
|
|
refl.predicted_x,
|
|
refl.predicted_y
|
|
});
|
|
}
|
|
|
|
XtalOptimizerData xtal_opt;
|
|
xtal_opt.latt = latt_i;
|
|
xtal_opt.geom.BeamX_pxl(1007).BeamY_pxl(990).DetectorDistance_mm(200)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Cubic;
|
|
xtal_opt.centering = 'I';
|
|
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, spots));
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
auto uc_o = xtal_opt.latt.GetUnitCell();
|
|
|
|
std::cout << "Beam center: " << xtal_opt.geom.GetBeamX_pxl() << " " << xtal_opt.geom.GetBeamY_pxl() << std::endl;
|
|
std::cout << "Unit cell: " << uc_o.a << " " << uc_o.b << " " << uc_o.c << std::endl;
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.1);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.1);
|
|
|
|
CHECK(fabsf(uc_i.a - uc_o.a) < 0.1);
|
|
CHECK(fabsf(uc_i.b - uc_o.b) < 0.1);
|
|
CHECK(fabsf(uc_i.c - uc_o.c) < 0.2);
|
|
CHECK(fabs(uc_o.alpha - expected_angle) < 0.02);
|
|
CHECK(fabs(uc_o.beta - expected_angle) < 0.02);
|
|
CHECK(fabs(uc_o.gamma - expected_angle) < 0.02);
|
|
}
|
|
|
|
TEST_CASE("LatticeToRodrigues") {
|
|
double rod[3];
|
|
double lengths[3];
|
|
|
|
CrystalLattice latt_i(40,50,80,90,90,90);
|
|
|
|
LatticeToRodriguesAndLengths_GS(latt_i, rod, lengths);
|
|
CHECK(lengths[0] == Catch::Approx(40.0));
|
|
CHECK(lengths[1] == Catch::Approx(50.0));
|
|
CHECK(lengths[2] == Catch::Approx(80.0));
|
|
CHECK(fabs(rod[0]) < 0.001);
|
|
CHECK(fabs(rod[1]) < 0.001);
|
|
CHECK(fabs(rod[2]) < 0.001);
|
|
|
|
auto latt_o = AngleAxisAndLengthsToLattice(rod, lengths, false);
|
|
|
|
CHECK(latt_o.Vec0().Length() == Catch::Approx(40.0));
|
|
CHECK(latt_o.Vec1().Length() == Catch::Approx(50.0));
|
|
CHECK(latt_o.Vec2().Length() == Catch::Approx(80.0));
|
|
}
|
|
|
|
TEST_CASE("LatticeToRodrigues_irregular") {
|
|
double rod[3];
|
|
double lengths[3];
|
|
|
|
CrystalLattice latt_i(Coord(40,0,0),
|
|
Coord(0, 50 / sqrt(2), -50 / sqrt(2)),
|
|
Coord(0, 80 / sqrt(2), 80 / sqrt(2)));
|
|
|
|
LatticeToRodriguesAndLengths_GS(latt_i, rod, lengths);
|
|
CHECK(lengths[0] == Catch::Approx(40.0));
|
|
CHECK(lengths[1] == Catch::Approx(50.0));
|
|
CHECK(lengths[2] == Catch::Approx(80.0));
|
|
|
|
|
|
auto latt_o = AngleAxisAndLengthsToLattice(rod, lengths, false);
|
|
|
|
CHECK(latt_o.Vec0().Length() == Catch::Approx(40.0));
|
|
CHECK(latt_o.Vec1().Length() == Catch::Approx(50.0));
|
|
CHECK(latt_o.Vec2().Length() == Catch::Approx(80.0));
|
|
}
|
|
|
|
TEST_CASE("LatticeToRodrigues_Hex") {
|
|
double rod[3];
|
|
double lengths[3];
|
|
|
|
Coord a = Coord(40,0,0);
|
|
Coord b = Coord(40 / 2, 40 * sqrt(3)/ 2.0, 0);
|
|
Coord c = Coord(0, 0, 70);
|
|
|
|
RotMatrix R(1.0, Coord(0,1,1));
|
|
CrystalLattice latt_i(R*a,R*b,R*c);
|
|
|
|
LatticeToRodriguesAndLengths_Hex(latt_i, rod, lengths);
|
|
CHECK(lengths[0] == Catch::Approx(40.0));
|
|
CHECK(lengths[1] == Catch::Approx(40.0));
|
|
CHECK(lengths[2] == Catch::Approx(70.0));
|
|
|
|
auto latt_o = AngleAxisAndLengthsToLattice(rod, lengths, true);
|
|
auto uc_o = latt_o.GetUnitCell();
|
|
CHECK(uc_o.a == Catch::Approx(40.0));
|
|
CHECK(uc_o.b == Catch::Approx(40.0));
|
|
CHECK(uc_o.c == Catch::Approx(70.0));
|
|
CHECK(uc_o.alpha == Catch::Approx(90.0));
|
|
CHECK(uc_o.beta == Catch::Approx(90.0));
|
|
CHECK(uc_o.gamma == Catch::Approx(120.0));
|
|
}
|