remove debug statement

CMakeLists.txt

@@ -205,6 +205,8 @@ if(NOT PCRE2_FOUND)
     add_subdirectory(pcre2-simple)
 endif()
 
+find_package(glm REQUIRED)
+
 # Using Eigen3 is a bit of a thing. We don't want to build it completely since
 # we only need a couple of header files. Nothing special. But often, eigen3 is
 # already installed and then we prefer that.

docs/symmetry.rst

@@ -8,7 +8,7 @@ Points
 
 The most basic type in use is :cpp:type:`cif::point`. It can be thought of as a point in space with three coordinates, but it is also often used as a vector in 3d space. To keep the interface simple there's no separate vector type.
 
-Many functions are available in :ref:`file_cif++_point.hpp` that work on points. There are functions to calculate the :cpp:func:`cif::distance` between two points and also function to calculate dot products, cross products and dihedral angles between sets of points.
+Many functions are available in :ref:`file_cif++_point.hpp` that work on points. There are functions to calculate the :cpp:func:`glm::distance` between two points and also function to calculate dot products, cross products and dihedral angles between sets of points.
 
 Quaternions
 -----------
@@ -91,7 +91,7 @@ To give an idea how this works, here's a piece of code copied from one of the un
         auto sa2 = c.symmetry_copy(p2, cif::sym_op(symm2));
 
         // The distance between these symmetry atoms should be equal to the distance in the struct_conn record
-        assert(cif::distance(sa1, sa2) == dist);
+        assert(glm::distance(sa1, sa2) == dist);
 
         // And to show how you can obtain the closest symmetry copy of an atom near another one:
         // here we request the symmetry copy of p2 that lies closest to p1

include/cif++/cif++.hpp

@@ -40,7 +40,6 @@
 #include "cif++/gzio.hpp"
 #include "cif++/item.hpp"
 #include "cif++/iterator.hpp"
-#include "cif++/matrix.hpp"
 #include "cif++/model.hpp"
 #include "cif++/parser.hpp"
 #include "cif++/pdb.hpp"

include/cif++/matrix.hpp

@@ -26,751 +26,4 @@
 
 #pragma once
 
-#include <array>
+#warning "Using this file is deprecated"
-#include <cassert>
-#include <cmath>
-#include <ostream>
-#include <type_traits>
-#include <vector>
-
-/**
- * @file matrix.hpp
- *
- * Some basic matrix operations and classes to hold matrices.
- *
- * We're using expression templates for optimal performance.
- *
- */
-
-namespace cif
-{
-// --------------------------------------------------------------------
-// We're using expression templates here
-
-/**
- * @brief Base for the matrix expression templates
- * This all uses the Curiously recurring template pattern
- *
- * @tparam M The type of the derived class
- */
-template <typename M>
-class matrix_expression // NOLINT(bugprone-crtp-constructor-accessibility)
-{
-  public:
-	[[nodiscard]] constexpr std::size_t dim_m() const { return static_cast<const M &>(*this).dim_m(); } ///< Return the size (dimension) in direction m
-	[[nodiscard]] constexpr std::size_t dim_n() const { return static_cast<const M &>(*this).dim_n(); } ///< Return the size (dimension) in direction n
-
-	[[nodiscard]] constexpr bool empty() const { return dim_m() == 0 or dim_n() == 0; } ///< Convenient way to test for empty matrices
-
-	/** Return a reference to element [ @a i, @a j ] */
-	[[nodiscard]] constexpr auto &operator()(std::size_t i, std::size_t j)
-	{
-		return static_cast<M &>(*this).operator()(i, j);
-	}
-
-	/** Return the value of element [ @a i, @a j ] */
-	[[nodiscard]] constexpr auto operator()(std::size_t i, std::size_t j) const
-	{
-		return static_cast<const M &>(*this).operator()(i, j);
-	}
-
-	/** Swap the contents of rows @a r1 and @a r2 */
-	void swap_row(std::size_t r1, std::size_t r2)
-	{
-		for (std::size_t c = 0; c < dim_m(); ++c)
-		{
-			auto v = operator()(r1, c);
-			operator()(r1, c) = operator()(r2, c);
-			operator()(r2, c) = v;
-		}
-	}
-
-	/** Swap the contents of columns @a c1 and @a c2 */
-	void swap_col(std::size_t c1, std::size_t c2)
-	{
-		for (std::size_t r = 0; r < dim_n(); ++r)
-		{
-			auto &a = operator()(r, c1);
-			auto &b = operator()(r, c2);
-			std::swap(a, b);
-		}
-	}
-
-	/** write the matrix @a m to std::ostream @a os */
-	friend std::ostream &operator<<(std::ostream &os, const matrix_expression &m)
-	{
-		os << '[';
-
-		for (std::size_t i = 0; i < m.dim_m(); ++i)
-		{
-			os << '[';
-
-			for (std::size_t j = 0; j < m.dim_n(); ++j)
-			{
-				os << m(i, j);
-				if (j + 1 < m.dim_n())
-					os << ", ";
-			}
-
-			if (i + 1 < m.dim_m())
-				os << ", ";
-
-			os << ']';
-		}
-
-		os << ']';
-
-		return os;
-	}
-
-	/// compare two matrices
-	template <typename M2>
-	constexpr bool operator==(const matrix_expression<M2> &m) const
-	{
-		bool same = false;
-		if (dim_m() == m.dim_m() and dim_n() == m.dim_n())
-		{
-			same = true;
-			for (std::size_t i = 0; same and i < m.dim_m(); ++i)
-			{
-				for (std::size_t j = 0; same and j < m.dim_n(); ++j)
-					same = operator()(i, j) == m(i, j);
-			}
-		}
-
-		return same;
-	}
-};
-
-// --------------------------------------------------------------------
-
-/**
- * @brief Storage class implementation of matrix_expression.
- *
- * @tparam F The type of the stored values
- *
- * matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
- * element m i,j is mapped to [i * n + j] and thus storage is row major
- */
-
-template <typename F = float>
-class matrix : public matrix_expression<matrix<F>>
-{
-  public:
-	/** The value type */
-	using value_type = F;
-
-	/**
-	 * @brief Copy construct a new matrix object using @a m
-	 *
-	 * @tparam M2 Type of @a m
-	 * @param m The matrix expression to copy values from
-	 */
-	template <typename M2>
-	matrix(const matrix_expression<M2> &m)
-		: m_m(m.dim_m())
-		, m_n(m.dim_n())
-		, m_data(m_m * m_n)
-	{
-		for (std::size_t i = 0; i < m_m; ++i)
-		{
-			for (std::size_t j = 0; j < m_n; ++j)
-				operator()(i, j) = m(i, j);
-		}
-	}
-
-	/**
-	 * @brief Construct a new matrix object with dimension @a m and @a n
-	 * setting the values to @a v
-	 *
-	 * @param m Requested dimension M
-	 * @param n Requested dimension N
-	 * @param v Value to store in each element
-	 */
-	matrix(std::size_t m, std::size_t n, value_type v = 0)
-		: m_m(m)
-		, m_n(n)
-		, m_data(m_m * m_n)
-	{
-		std::fill(m_data.begin(), m_data.end(), v);
-	}
-
-	/** @cond */
-	matrix() = default;
-	matrix(matrix &&m) = default;
-	matrix(const matrix &m) = default;
-	matrix &operator=(matrix &&m) = default;
-	matrix &operator=(const matrix &m) = default;
-	/** @endcond */
-
-	[[nodiscard]] constexpr std::size_t dim_m() const { return m_m; } ///< Return dimension m
-	[[nodiscard]] constexpr std::size_t dim_n() const { return m_n; } ///< Return dimension n
-
-	/** Return the value of element [ @a i, @a j ] */
-	[[nodiscard]] constexpr value_type operator()(std::size_t i, std::size_t j) const
-	{
-		assert(i < m_m);
-		assert(j < m_n);
-		return m_data[i * m_n + j];
-	}
-
-	/** Return a reference to element [ @a i, @a j ] */
-	[[nodiscard]] constexpr value_type &operator()(std::size_t i, std::size_t j)
-	{
-		assert(i < m_m);
-		assert(j < m_n);
-		return m_data[i * m_n + j];
-	}
-
-  private:
-	std::size_t m_m = 0, m_n = 0;
-	std::vector<value_type> m_data;
-};
-
-// --------------------------------------------------------------------
-// special case, 3x3 matrix
-
-/**
- * @brief Storage class implementation of matrix_expression
- * with compile time fixed size.
- *
- * @tparam F The type of the stored values
- *
- * matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
- * element m i,j is mapped to [i * n + j] and thus storage is row major
- */
-
-template <typename F, std::size_t M, std::size_t N>
-class matrix_fixed : public matrix_expression<matrix_fixed<F, M, N>>
-{
-  public:
-	/** The value type */
-	using value_type = F;
-
-	/** The storage size */
-	static constexpr std::size_t kSize = M * N;
-
-	/** Copy constructor */
-	template <typename M2>
-	matrix_fixed(const M2 &m)
-	{
-		assert(M == m.dim_m() and N == m.dim_n());
-		for (std::size_t i = 0; i < M; ++i)
-		{
-			for (std::size_t j = 0; j < N; ++j)
-				operator()(i, j) = m(i, j);
-		}
-	}
-
-	/** default constructor */
-	matrix_fixed(value_type v = 0)
-	{
-		m_data.fill(v);
-	}
-
-	/** Alternate constructor taking an array of values to store */
-	matrix_fixed(const F (&v)[kSize])
-	{
-		fill(v, std::make_index_sequence<kSize>{});
-	}
-
-	/** @cond */
-	matrix_fixed(matrix_fixed &&m) = default;
-	matrix_fixed(const matrix_fixed &m) = default;
-	matrix_fixed &operator=(matrix_fixed &&m) = default;
-	matrix_fixed &operator=(const matrix_fixed &m) = default;
-	/** @endcond */
-
-	/** Store the values in @a a in the matrix */
-	template <std::size_t... Ixs>
-	matrix_fixed &fill(const F (&a)[kSize], std::index_sequence<Ixs...>)
-	{
-		m_data = { a[Ixs]... };
-		return *this;
-	}
-
-	[[nodiscard]] constexpr std::size_t dim_m() const { return M; } ///< Return dimension m
-	[[nodiscard]] constexpr std::size_t dim_n() const { return N; } ///< Return dimension n
-
-	/** Return the value of element [ @a i, @a j ] */
-	[[nodiscard]] constexpr value_type operator()(std::size_t i, std::size_t j) const
-	{
-		assert(i < M);
-		assert(j < N);
-		return m_data[i * N + j];
-	}
-
-	/** Return a reference to element [ @a i, @a j ] */
-	[[nodiscard]] constexpr value_type &operator()(std::size_t i, std::size_t j)
-	{
-		assert(i < M);
-		assert(j < N);
-		return m_data[i * N + j];
-	}
-
-  private:
-	std::array<value_type, M * N> m_data;
-};
-
-/** typedef of a fixed matrix of size 3x3 */
-template <typename F>
-using matrix3x3 = matrix_fixed<F, 3, 3>;
-
-/** typedef of a fixed matrix of size 4x4 */
-template <typename F>
-using matrix4x4 = matrix_fixed<F, 4, 4>;
-
-// --------------------------------------------------------------------
-
-/**
- * @brief Storage class implementation of symmetric matrix_expression
- *
- * @tparam F The type of the stored values
- *
- * matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
- * element m i,j is mapped to [i * n + j] and thus storage is row major
- */
-template <typename F = float>
-class symmetric_matrix : public matrix_expression<symmetric_matrix<F>>
-{
-  public:
-	/** The value type */
-	using value_type = F;
-
-	/** constructor for a matrix of size @a n x @a n elements with value @a v */
-	symmetric_matrix(std::size_t n, value_type v = 0)
-		: m_n(n)
-		, m_data((m_n * (m_n + 1)) / 2)
-	{
-		std::fill(m_data.begin(), m_data.end(), v);
-	}
-
-	/** @cond */
-	symmetric_matrix() = default;
-	symmetric_matrix(symmetric_matrix &&m) = default;
-	symmetric_matrix(const symmetric_matrix &m) = default;
-	symmetric_matrix &operator=(symmetric_matrix &&m) = default;
-	symmetric_matrix &operator=(const symmetric_matrix &m) = default;
-	/** @endcond */
-
-	[[nodiscard]] constexpr std::size_t dim_m() const { return m_n; } ///< Return dimension m
-	[[nodiscard]] constexpr std::size_t dim_n() const { return m_n; } ///< Return dimension n
-
-	/** Return the value of element [ @a i, @a j ] */
-	[[nodiscard]] constexpr value_type operator()(std::size_t i, std::size_t j) const
-	{
-		return i < j
-		           ? m_data[(j * (j + 1)) / 2 + i]
-		           : m_data[(i * (i + 1)) / 2 + j];
-	}
-
-	/** Return a reference to element [ @a i, @a j ] */
-	[[nodiscard]] constexpr value_type &operator()(std::size_t i, std::size_t j)
-	{
-		if (i > j)
-			std::swap(i, j);
-		assert(j < m_n);
-		return m_data[(j * (j + 1)) / 2 + i];
-	}
-
-  private:
-	std::size_t m_n;
-	std::vector<value_type> m_data;
-};
-
-// --------------------------------------------------------------------
-
-/**
- * @brief Storage class implementation of symmetric matrix_expression
- * with compile time fixed size.
- *
- * @tparam F The type of the stored values
- *
- * matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
- * element m i,j is mapped to [i * n + j] and thus storage is row major
- */
-template <typename F, std::size_t M>
-class symmetric_matrix_fixed : public matrix_expression<symmetric_matrix_fixed<F, M>>
-{
-  public:
-	/** The value type */
-	using value_type = F;
-
-	/** constructor with all elements set to value @a v */
-	symmetric_matrix_fixed(value_type v = 0)
-	{
-		std::fill(m_data.begin(), m_data.end(), v);
-	}
-
-	/** @cond */
-	symmetric_matrix_fixed(symmetric_matrix_fixed &&m) = default;
-	symmetric_matrix_fixed(const symmetric_matrix_fixed &m) = default;
-	symmetric_matrix_fixed &operator=(symmetric_matrix_fixed &&m) = default;
-	symmetric_matrix_fixed &operator=(const symmetric_matrix_fixed &m) = default;
-	/** @endcond */
-
-	[[nodiscard]] constexpr std::size_t dim_m() const { return M; } ///< Return dimension m
-	[[nodiscard]] constexpr std::size_t dim_n() const { return M; } ///< Return dimension n
-
-	/** Return the value of element [ @a i, @a j ] */
-	[[nodiscard]] constexpr value_type operator()(std::size_t i, std::size_t j) const
-	{
-		return i < j
-		           ? m_data[(j * (j + 1)) / 2 + i]
-		           : m_data[(i * (i + 1)) / 2 + j];
-	}
-
-	/** Return a reference to element [ @a i, @a j ] */
-	[[nodiscard]] constexpr value_type &operator()(std::size_t i, std::size_t j)
-	{
-		if (i > j)
-			std::swap(i, j);
-		assert(j < M);
-		return m_data[(j * (j + 1)) / 2 + i];
-	}
-
-  private:
-	std::array<value_type, (M * (M + 1)) / 2> m_data;
-};
-
-/** typedef of a fixed symmetric matrix of size 3x3 */
-template <typename F>
-using symmetric_matrix3x3 = symmetric_matrix_fixed<F, 3>;
-
-/** typedef of a fixed symmetric matrix of size 4x4 */
-template <typename F>
-using symmetric_matrix4x4 = symmetric_matrix_fixed<F, 4>;
-
-// --------------------------------------------------------------------
-
-/// A transposed matrix view
-
-template <typename M>
-class transposed_matrix : public cif::matrix_expression<transposed_matrix<M>>
-{
-  public:
-	transposed_matrix(const M &m)
-		: m_m(m)
-	{
-	}
-
-	[[nodiscard]] constexpr std::size_t dim_m() const { return m_m.dim_n(); } ///< Return dimension m
-	[[nodiscard]] constexpr std::size_t dim_n() const { return m_m.dim_m(); } ///< Return dimension n
-
-	/** Access to the value of element [ @a i, @a j ] */
-	[[nodiscard]] constexpr auto operator()(std::size_t i, std::size_t j) const
-	{
-		return m_m(j, i);
-	}
-
-  private:
-	const M &m_m;
-};
-
-// --------------------------------------------------------------------
-
-/**
- * @brief implementation of symmetric matrix_expression with a value
- * of 1 for the diagonal values and 0 for all the others.
- *
- * @tparam F The type of the stored values
- *
- * matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
- * element m i,j is mapped to [i * n + j] and thus storage is row major
- */
-template <typename F = float>
-class identity_matrix : public matrix_expression<identity_matrix<F>>
-{
-  public:
-	/** the value type */
-	using value_type = F;
-
-	/** constructor taking a dimension @a n */
-	identity_matrix(std::size_t n)
-		: m_n(n)
-	{
-	}
-
-	[[nodiscard]] constexpr std::size_t dim_m() const { return m_n; } ///< Return dimension m
-	[[nodiscard]] constexpr std::size_t dim_n() const { return m_n; } ///< Return dimension n
-
-	/** Return the value of element [ @a i, @a j ] */
-	[[nodiscard]] constexpr value_type operator()(std::size_t i, std::size_t j) const
-	{
-		return static_cast<value_type>(i == j ? 1 : 0);
-	}
-
-  private:
-	std::size_t m_n;
-};
-
-// --------------------------------------------------------------------
-// matrix functions, implemented as expression templates
-
-/**
- * @brief Implementation of a substraction operation as a matrix expression
- *
- * @tparam M1 Type of matrix 1
- * @tparam M2 Type of matrix 2
- */
-template <typename M1, typename M2>
-class matrix_subtraction : public matrix_expression<matrix_subtraction<M1, M2>>
-{
-  public:
-	/** constructor */
-	matrix_subtraction(const M1 &m1, const M2 &m2)
-		: m_m1(m1)
-		, m_m2(m2)
-	{
-		assert(m_m1.dim_m() == m_m2.dim_m());
-		assert(m_m1.dim_n() == m_m2.dim_n());
-	}
-
-	[[nodiscard]] constexpr std::size_t dim_m() const { return m_m1.dim_m(); } ///< Return dimension m
-	[[nodiscard]] constexpr std::size_t dim_n() const { return m_m1.dim_n(); } ///< Return dimension n
-
-	/** Access to the value of element [ @a i, @a j ] */
-	[[nodiscard]] constexpr auto operator()(std::size_t i, std::size_t j) const
-	{
-		return m_m1(i, j) - m_m2(i, j);
-	}
-
-  private:
-	const M1 &m_m1;
-	const M2 &m_m2;
-};
-
-/** operator to subtract two matrices and return a matrix expression */
-template <typename M1, typename M2>
-auto operator-(const matrix_expression<M1> &m1, const matrix_expression<M2> &m2)
-{
-	return matrix_subtraction(m1, m2);
-}
-
-/**
- * @brief Implementation of a multiplication operation as a matrix expression
- *
- * @tparam M1 Type of matrix 1
- * @tparam M2 Type of matrix 2
- */
-template <typename M1, typename M2>
-class matrix_matrix_multiplication : public matrix_expression<matrix_matrix_multiplication<M1, M2>>
-{
-  public:
-	/** constructor */
-	matrix_matrix_multiplication(const M1 &m1, const M2 &m2)
-		: m_m1(m1)
-		, m_m2(m2)
-	{
-		assert(m1.dim_n() == m2.dim_m());
-	}
-
-	[[nodiscard]] constexpr std::size_t dim_m() const { return m_m1.dim_m(); } ///< Return dimension m
-	[[nodiscard]] constexpr std::size_t dim_n() const { return m_m1.dim_n(); } ///< Return dimension n
-
-	/** Access to the value of element [ @a i, @a j ] */
-	[[nodiscard]] constexpr auto operator()(std::size_t i, std::size_t j) const
-	{
-		using value_type = decltype(m_m1(0, 0));
-
-		value_type result = {};
-
-		for (std::size_t k = 0; k < m_m1.dim_n(); ++k)
-			result += m_m1(i, k) * m_m2(k, j);
-
-		return result;
-	}
-
-  private:
-	const M1 &m_m1;
-	const M2 &m_m2;
-};
-
-/**
- * @brief Implementation of a multiplication operation of a matrix and a scalar value as a matrix expression
- *
- * @tparam M1 Type of matrix
- * @tparam M2 Type of scalar value
- */
-template <typename M, typename T>
-class matrix_scalar_multiplication : public matrix_expression<matrix_scalar_multiplication<M, T>>
-{
-  public:
-	/** value type */
-	using value_type = T;
-
-	/** constructor */
-	matrix_scalar_multiplication(const M &m, value_type v)
-		: m_m(m)
-		, m_v(v)
-	{
-	}
-
-	[[nodiscard]] constexpr std::size_t dim_m() const { return m_m.dim_m(); } ///< Return dimension m
-	[[nodiscard]] constexpr std::size_t dim_n() const { return m_m.dim_n(); } ///< Return dimension n
-
-	/** Access to the value of element [ @a i, @a j ] */
-	[[nodiscard]] constexpr auto operator()(std::size_t i, std::size_t j) const
-	{
-		return m_m(i, j) * m_v;
-	}
-
-  private:
-	const M &m_m;
-	value_type m_v;
-};
-
-/** First implementation of operator*, enabled if the second parameter is a scalar */
-template <typename M1, typename T>
-auto operator*(const matrix_expression<M1> &m, T v)
-	requires(std::is_floating_point_v<T>)
-{
-	return matrix_scalar_multiplication(m, v);
-}
-
-/** First implementation of operator*, enabled if the second parameter is not a scalar and thus must be a matrix, right? */
-template <typename M1, typename M2>
-auto operator*(const matrix_expression<M1> &m1, const matrix_expression<M2> &m2)
-	requires(not std::is_floating_point_v<M2>)
-{
-	return matrix_matrix_multiplication(m1, m2);
-}
-
-// --------------------------------------------------------------------
-
-/// A sub-view on a matrix
-template <typename M2>
-class sub_matrix : public matrix_expression<sub_matrix<M2>>
-{
-  public:
-	/// Constructor
-	sub_matrix(const M2 &m, int i, int j)
-		: m_m(m)
-		, m_i(i)
-		, m_j(j)
-	{
-	}
-
-	[[nodiscard]] constexpr std::size_t dim_m() const { return m_m.dim_m() - 1; } ///< Return dimension m
-	[[nodiscard]] constexpr std::size_t dim_n() const { return m_m.dim_n() - 1; } ///< Return dimension n
-
-	/** Access to the value of element [ @a i, @a j ] */
-	[[nodiscard]] constexpr auto operator()(std::size_t i, std::size_t j) const
-	{
-		return m_m(
-			i >= m_i ? i + 1 : i,
-			j >= m_j ? j + 1 : j);
-	}
-
-  private:
-	const M2 &m_m;
-	std::size_t m_i, m_j;
-};
-
-// --------------------------------------------------------------------
-
-/** Generic routine to calculate the determinant of a matrix
- *
- * @note This is currently only implemented for fixed matrices of size 3x3
- */
-template <typename M>
-auto determinant(const M &m);
-
-/** Implementation of the determinant function for fixed size matrices of size 3x3 */
-template <typename F = float>
-auto determinant(const matrix3x3<F> &m)
-{
-	return (m(0, 0) * ((m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1))) +
-			m(0, 1) * ((m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2))) +
-			m(0, 2) * ((m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0))));
-}
-
-/** Implementation of the determinant function for fixed size matrices of size 4x4 */
-template <typename F = float>
-F determinant(const matrix4x4<F> &m)
-{
-	return m(0, 0) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 0))) -
-	       m(0, 1) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 1))) +
-	       m(0, 2) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 2))) -
-	       m(0, 3) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 3)));
-}
-
-// --------------------------------------------------------------------
-
-/** Generic routine to calculate the inverse of a matrix
- *
- * @note This is currently only implemented for fixed matrices of size 3x3
- */
-template <typename M>
-M inverse(const M &m);
-
-/** Implementation of the inverse function for fixed size matrices of size 3x3 */
-template <typename F = float>
-matrix3x3<F> inverse(const matrix3x3<F> &m)
-{
-	F det = determinant(m);
-
-	matrix3x3<F> result;
-
-	result(0, 0) = (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1)) / det;
-	result(1, 0) = (m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2)) / det;
-	result(2, 0) = (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0)) / det;
-	result(0, 1) = (m(2, 1) * m(0, 2) - m(2, 2) * m(0, 1)) / det;
-	result(1, 1) = (m(2, 2) * m(0, 0) - m(2, 0) * m(0, 2)) / det;
-	result(2, 1) = (m(2, 0) * m(0, 1) - m(2, 1) * m(0, 0)) / det;
-	result(0, 2) = (m(0, 1) * m(1, 2) - m(0, 2) * m(1, 1)) / det;
-	result(1, 2) = (m(0, 2) * m(1, 0) - m(0, 0) * m(1, 2)) / det;
-	result(2, 2) = (m(0, 0) * m(1, 1) - m(0, 1) * m(1, 0)) / det;
-
-	return result;
-}
-
-// --------------------------------------------------------------------
-
-/**
- * @brief Implementation of a cofactor calculation as a matrix expression
- *
- * @tparam M Type of matrix
- */
-template <typename M>
-class matrix_cofactors : public matrix_expression<matrix_cofactors<M>>
-{
-  public:
-	/** constructor */
-	matrix_cofactors(const M &m)
-		: m_m(m)
-	{
-	}
-
-	[[nodiscard]] constexpr std::size_t dim_m() const { return m_m.dim_m(); } ///< Return dimension m
-	[[nodiscard]] constexpr std::size_t dim_n() const { return m_m.dim_n(); } ///< Return dimension n
-
-	/** Access to the value of element [ @a i, @a j ] */
-	[[nodiscard]] constexpr auto operator()(std::size_t i, std::size_t j) const
-	{
-		const std::size_t ixs[4][3] = {
-			{ 1, 2, 3 },
-			{ 0, 2, 3 },
-			{ 0, 1, 3 },
-			{ 0, 1, 2 }
-		};
-
-		const std::size_t *ix = ixs[i];
-		const std::size_t *iy = ixs[j];
-
-		auto result =
-			m_m(ix[0], iy[0]) * m_m(ix[1], iy[1]) * m_m(ix[2], iy[2]) +
-			m_m(ix[0], iy[1]) * m_m(ix[1], iy[2]) * m_m(ix[2], iy[0]) +
-			m_m(ix[0], iy[2]) * m_m(ix[1], iy[0]) * m_m(ix[2], iy[1]) -
-			m_m(ix[0], iy[2]) * m_m(ix[1], iy[1]) * m_m(ix[2], iy[0]) -
-			m_m(ix[0], iy[1]) * m_m(ix[1], iy[0]) * m_m(ix[2], iy[2]) -
-			m_m(ix[0], iy[0]) * m_m(ix[1], iy[2]) * m_m(ix[2], iy[1]);
-
-		return (i + j) % 2 == 1 ? -result : result;
-	}
-
-  private:
-	const M &m_m;
-};
-
-} // namespace cif

include/cif++/model.hpp

@@ -99,7 +99,7 @@ class atom
 		{
 			auto r = row();
 			if (r)
-				std::tie(m_location.m_x, m_location.m_y, m_location.m_z) = r.get<float, float, float>("Cartn_x", "Cartn_y", "Cartn_z");
+				std::tie(m_location.x, m_location.y, m_location.z) = r.get<float, float, float>("Cartn_x", "Cartn_y", "Cartn_z");
 		}
 
 		// constructor for a symmetry copy of an atom
@@ -298,16 +298,16 @@ class atom
 	/// \brief Rotate the position of this atom by \a q
 	void rotate(quaternion q)
 	{
-		auto loc = get_location();
+		set_location(q * get_location());
-		loc.rotate(q);
-		set_location(loc);
 	}
 
 	/// \brief rotate the coordinates of this atom by \a q around point \a p
 	void rotate(quaternion q, point p)
 	{
 		auto loc = get_location();
-		loc.rotate(q, p);
+		loc -= p;
+		loc = q * loc;
+		loc += p;
 		set_location(loc);
 	}
 
@@ -316,7 +316,7 @@ class atom
 	{
 		auto loc = get_location();
 		loc += t;
-		loc.rotate(q);
+		loc = q * loc;
 		set_location(loc);
 	}
 
@@ -325,7 +325,7 @@ class atom
 	{
 		auto loc = get_location();
 		loc += t1;
-		loc.rotate(q);
+		loc = q * loc;
 		loc += t2;
 		set_location(loc);
 	}
@@ -453,16 +453,6 @@ inline float distance(const atom &a, const atom &b)
 	return distance(a.get_location(), b.get_location());
 }
 
-/** Calculate the square of the distance between atoms @a and @a b in ångström
- *
- * @note Use this whenever possible instead of simply using distance since
- * this function does not have to calculate a square root which is expensive.
- */
-inline float distance_squared(const atom &a, const atom &b)
-{
-	return distance_squared(a.get_location(), b.get_location());
-}
-
 // --------------------------------------------------------------------
 
 /**

include/cif++/point.hpp

@@ -26,21 +26,19 @@
 
 #pragma once
 
-#include <array>
 #include <cmath>
-#include <complex>
-#include <cstdint>
 #include <cstdlib>
 #include <format>
-#include <functional>
+#include <glm/ext/quaternion_common.hpp>
+#include <glm/ext/quaternion_geometric.hpp>
+#include <glm/ext/quaternion_trigonometric.hpp>
+#include <glm/glm.hpp>
+#include <glm/gtc/quaternion.hpp>
+#include <glm/trigonometric.hpp>
 #include <limits>
 #include <numbers>
 #include <optional>
-#include <ostream>
 #include <tuple>
-#include <type_traits>
-#include <utility>
-#include <valarray>
 #include <vector>
 
 #if __has_include(<clipper/core/coords.h>)
@@ -54,659 +52,36 @@
 
 /** \file point.hpp
  *
- * This file contains the definition for *cif::point* as well as
+ * This file contains the definition for *point* as well as
  * lots of routines and classes that can manipulate points.
  */
 
 namespace cif
 {
 
-// --------------------------------------------------------------------
+// Using glm now
-/**
- * @brief A stripped down quaternion implementation, based on boost::math::quaternion
- *
- * We use quaternions to do rotations in 3d space. Quaternions are faster than
- * matrix calculations and they also suffer less from drift caused by rounding
- * errors.
- *
- * Like complex number, quaternions do have a meaningful notion of "real part",
- * but unlike them there is no meaningful notion of "imaginary part".
- * Instead there is an "unreal part" which itself is a quaternion, and usually
- * nothing simpler (as opposed to the complex number case).
- * However, for practicality, there are accessors for the other components
- * (these are necessary for the templated copy constructor, for instance).
- *
- * @note Quaternion multiplication is *NOT* commutative;
- * symbolically, "q *= rhs;" means "q = q * rhs;"
- * and "q /= rhs;" means "q = q * inverse_of(rhs);"
- */
 
 template <typename T>
-class quaternion_type
+using quaternion_type = glm::qua<T>;
-{
-  public:
-	/// \brief the value type of the elements, usually this is float
-	using value_type = T;
 
-	/// \brief constructor with the four members
-	constexpr explicit quaternion_type(value_type const &value_a = {}, value_type const &value_b = {}, value_type const &value_c = {}, value_type const &value_d = {})
-		: a(value_a)
-		, b(value_b)
-		, c(value_c)
-		, d(value_d)
-	{
-	}
-
-	/// \brief constructor taking two complex values as input
-	constexpr explicit quaternion_type(std::complex<value_type> const &z0, std::complex<value_type> const &z1 = std::complex<value_type>())
-		: a(z0.real())
-		, b(z0.imag())
-		, c(z1.real())
-		, d(z1.imag())
-	{
-	}
-
-	constexpr quaternion_type(quaternion_type const &) = default; ///< Copy constructor
-	constexpr quaternion_type(quaternion_type &&) = default;      ///< Copy constructor
-
-	/// \brief Copy constructor accepting a quaternion with a different value_type
-	template <typename X>
-	constexpr explicit quaternion_type(quaternion_type<X> const &rhs)
-		: a(static_cast<value_type>(rhs.a))
-		, b(static_cast<value_type>(rhs.b))
-		, c(static_cast<value_type>(rhs.c))
-		, d(static_cast<value_type>(rhs.d))
-	{
-	}
-
-	// accessors
-
-	/// \brief See class description, return the *real* part of the quaternion
-	[[nodiscard]] constexpr value_type real() const
-	{
-		return a;
-	}
-
-	/// \brief See class description, return the *unreal* part of the quaternion
-	[[nodiscard]] constexpr quaternion_type unreal() const
-	{
-		return { 0, b, c, d };
-	}
-
-	/// \brief swap
-	constexpr void swap(quaternion_type &o)
-	{
-		std::swap(a, o.a);
-		std::swap(b, o.b);
-		std::swap(c, o.c);
-		std::swap(d, o.d);
-	}
-
-	// assignment operators
-
-	/// \brief Assignment operator accepting a quaternion with optionally another value_type
-	template <typename X>
-	constexpr quaternion_type &operator=(quaternion_type<X> const &rhs)
-	{
-		a = static_cast<value_type>(rhs.a);
-		b = static_cast<value_type>(rhs.b);
-		c = static_cast<value_type>(rhs.c);
-		d = static_cast<value_type>(rhs.d);
-
-		return *this;
-	}
-
-	/// \brief Assignment operator
-	constexpr quaternion_type &operator=(quaternion_type const &rhs) = default;
-
-	/// \brief Assignment operator that sets the *real* part to @a rhs and the *unreal* parts to zero
-	constexpr quaternion_type &operator=(value_type const &rhs)
-	{
-		a = rhs;
-
-		b = c = d = static_cast<value_type>(0);
-
-		return *this;
-	}
-
-	/// \brief Assignment operator that sets the *real* part to the real part of @a rhs
-	/// and the first *unreal* part to the imaginary part of of @a rhs. The other *unreal*
-	// parts are set to zero.
-	constexpr quaternion_type &operator=(std::complex<value_type> const &rhs)
-	{
-		a = rhs.real();
-		b = rhs.imag();
-
-		c = d = static_cast<value_type>(0);
-
-		return *this;
-	}
-
-	// other assignment-related operators
-
-	/// \brief operator += adding value @a rhs to the *real* part
-	constexpr quaternion_type &operator+=(value_type const &rhs)
-	{
-		a += rhs;
-		return *this;
-	}
-
-	/// \brief operator += adding the real part of @a rhs to the *real* part
-	/// and the imaginary part of @a rhs to the first *unreal* part
-	constexpr quaternion_type &operator+=(std::complex<value_type> const &rhs)
-	{
-		a += std::real(rhs);
-		b += std::imag(rhs);
-		return *this;
-	}
-
-	/// \brief operator += adding the parts of @a rhs to the equivalent part of this
-	template <class X>
-	constexpr quaternion_type &operator+=(quaternion_type<X> const &rhs)
-	{
-		a += rhs.a;
-		b += rhs.b;
-		c += rhs.c;
-		d += rhs.d;
-		return *this;
-	}
-
-	/// \brief operator -= subtracting value @a rhs from the *real* part
-	constexpr quaternion_type &operator-=(value_type const &rhs)
-	{
-		a -= rhs;
-		return *this;
-	}
-
-	/// \brief operator -= subtracting the real part of @a rhs from the *real* part
-	/// and the imaginary part of @a rhs from the first *unreal* part
-	constexpr quaternion_type &operator-=(std::complex<value_type> const &rhs)
-	{
-		a -= std::real(rhs);
-		b -= std::imag(rhs);
-		return *this;
-	}
-
-	/// \brief operator -= subtracting the parts of @a rhs from the equivalent part of this
-	template <class X>
-	constexpr quaternion_type &operator-=(quaternion_type<X> const &rhs)
-	{
-		a -= rhs.a;
-		b -= rhs.b;
-		c -= rhs.c;
-		d -= rhs.d;
-		return *this;
-	}
-
-	/// \brief multiply all parts with value @a rhs
-	constexpr quaternion_type &operator*=(value_type const &rhs)
-	{
-		a *= rhs;
-		b *= rhs;
-		c *= rhs;
-		d *= rhs;
-		return *this;
-	}
-
-	/// \brief multiply with complex number @a rhs
-	constexpr quaternion_type &operator*=(std::complex<value_type> const &rhs)
-	{
-		value_type ar = rhs.real();
-		value_type br = rhs.imag();
-		quaternion_type result(a * ar - b * br, a * br + b * ar, c * ar + d * br, -c * br + d * ar);
-		swap(result);
-		return *this;
-	}
-
-	/// \brief multiply @a a with @a b and return the result
-	friend constexpr quaternion_type operator*(const quaternion_type &a, const quaternion_type &b)
-	{
-		auto result = a;
-		result *= b;
-		return result;
-	}
-
-	/// \brief multiply with quaternion @a rhs
-	template <typename X>
-	constexpr quaternion_type &operator*=(quaternion_type<X> const &rhs)
-	{
-		auto ar = static_cast<value_type>(rhs.a);
-		auto br = static_cast<value_type>(rhs.b);
-		auto cr = static_cast<value_type>(rhs.c);
-		auto dr = static_cast<value_type>(rhs.d);
-
-		quaternion_type result(a * ar - b * br - c * cr - d * dr, a * br + b * ar + c * dr - d * cr, a * cr - b * dr + c * ar + d * br, a * dr + b * cr - c * br + d * ar);
-		swap(result);
-		return *this;
-	}
-
-	/// \brief divide all parts by @a rhs
-	constexpr quaternion_type &operator/=(value_type const &rhs)
-	{
-		a /= rhs;
-		b /= rhs;
-		c /= rhs;
-		d /= rhs;
-		return *this;
-	}
-
-	/// \brief divide by complex number @a rhs
-	constexpr quaternion_type &operator/=(std::complex<value_type> const &rhs)
-	{
-		value_type ar = rhs.real();
-		value_type br = rhs.imag();
-		value_type denominator = ar * ar + br * br;
-		quaternion_type result((+a * ar + b * br) / denominator, (-a * br + b * ar) / denominator, (+c * ar - d * br) / denominator, (+c * br + d * ar) / denominator);
-		swap(result);
-		return *this;
-	}
-
-	/// \brief divide by quaternion @a rhs
-	template <typename X>
-	constexpr quaternion_type &operator/=(quaternion_type<X> const &rhs)
-	{
-		auto ar = static_cast<value_type>(rhs.a);
-		auto br = static_cast<value_type>(rhs.b);
-		auto cr = static_cast<value_type>(rhs.c);
-		auto dr = static_cast<value_type>(rhs.d);
-
-		value_type denominator = ar * ar + br * br + cr * cr + dr * dr;
-		quaternion_type result((+a * ar + b * br + c * cr + d * dr) / denominator, (-a * br + b * ar - c * dr + d * cr) / denominator, (-a * cr + b * dr + c * ar - d * br) / denominator, (-a * dr - b * cr + c * br + d * ar) / denominator);
-		swap(result);
-		return *this;
-	}
-
-	/// \brief normalise the values so that the length of the result is exactly 1
-	friend constexpr quaternion_type normalize(quaternion_type q)
-	{
-		std::valarray<value_type> t(4);
-
-		t[0] = q.a;
-		t[1] = q.b;
-		t[2] = q.c;
-		t[3] = q.d;
-
-		t *= t;
-
-		value_type length = std::sqrt(t.sum());
-
-		if (length > 0.001)
-			q /= static_cast<value_type>(length);
-		else
-			q = quaternion_type(1, 0, 0, 0);
-
-		return q;
-	}
-
-	/// \brief return the conjugate of this
-	friend constexpr quaternion_type conj(quaternion_type q)
-	{
-		return quaternion_type{ +q.a, -q.b, -q.c, -q.d };
-	}
-
-	[[nodiscard]] constexpr value_type get_a() const { return a; } ///< Return part a
-	[[nodiscard]] constexpr value_type get_b() const { return b; } ///< Return part b
-	[[nodiscard]] constexpr value_type get_c() const { return c; } ///< Return part c
-	[[nodiscard]] constexpr value_type get_d() const { return d; } ///< Return part d
-
-	/// \brief compare with @a rhs
-	constexpr bool operator==(const quaternion_type &rhs) const
-	{
-		return a == rhs.a and b == rhs.b and c == rhs.c and d == rhs.d;
-	}
-
-	/// \brief compare with @a rhs
-	constexpr bool operator!=(const quaternion_type &rhs) const
-	{
-		return a != rhs.a or b != rhs.b or c != rhs.c or d != rhs.d;
-	}
-
-	/// \brief test for all zero values
-	constexpr explicit operator bool() const
-	{
-		return a != 0 or b != 0 or c != 0 or d != 0;
-	}
-
-	/// \brief for debugging e.g.
-	friend std::ostream &operator<<(std::ostream &os, const quaternion_type &rhs)
-	{
-		os << std::format("{{ a: {}, b: {}, c: {}, d: {} }}", rhs.a, rhs.b, rhs.c, rhs.d);
-		return os;
-	}
-
-  private:
-	value_type a, b, c, d;
-};
-
-/**
- * @brief This code is similar to the code in boost so I copy the documentation as well:
- *
- * > spherical is a simple transposition of polar, it takes as inputs a (positive)
- * > magnitude and a point on the hypersphere, given by three angles. The first of
- * > these, theta has a natural range of -pi to +pi, and the other two have natural
- * > ranges of -pi/2 to +pi/2 (as is the case with the usual spherical coordinates in
- * > **R**<sup>3</sup>). Due to the many symmetries and periodicities, nothing untoward happens if
- * > the magnitude is negative or the angles are outside their natural ranges. The
- * > expected degeneracies (a magnitude of zero ignores the angles settings...) do
- * > happen however.
- */
-
-template <typename T>
-inline quaternion_type<T> spherical(T const &rho, T const &theta, T const &phi1, T const &phi2)
-{
-	T cos_phi1 = std::cos(phi1);
-	T cos_phi2 = std::cos(phi2);
-
-	T a = std::cos(theta) * cos_phi1 * cos_phi2;
-	T b = std::sin(theta) * cos_phi1 * cos_phi2;
-	T c = std::sin(phi1) * cos_phi2;
-	T d = std::sin(phi2);
-
-	quaternion_type result(a, b, c, d);
-	result *= rho;
-
-	return result;
-}
-
-/// \brief By default we use the float version of a quaternion
 using quaternion = quaternion_type<float>;
 
+template <typename T>
+using point_type = glm::vec<3, T>;
+
+using point = point_type<float>;
+
 // --------------------------------------------------------------------
 
-/**
+constexpr point rotate(point pt, quaternion q)
- * @brief 3D point: a location with x, y and z coordinates as floating point.
- *
- * Note that you can simply use structured binding to get access to the
- * individual parts like so:
- *
- * @code{.cpp}
- * float x, y, z;
- * tie(x, y, z) = atom.get_location();
- * @endcode
- */
-
-template <typename F>
-struct point_type
 {
-	/// \brief the value type of the x, y and z members
+	return q * pt;
-	using value_type = F;
+}
 
-	value_type m_x, ///< The x part of the location
+constexpr point rotate(point pt, quaternion q, point a)
-		m_y,        ///< The y part of the location
+{
-		m_z;        ///< The z part of the location
+	return q * (pt - a) + a;
-
+}
-	/// \brief default constructor, initialises the values to zero
-	constexpr point_type()
-		: m_x(0)
-		, m_y(0)
-		, m_z(0)
-	{
-	}
-
-	/// \brief constructor taking three values
-	constexpr point_type(value_type x, value_type y, value_type z)
-		: m_x(x)
-		, m_y(y)
-		, m_z(z)
-	{
-	}
-
-	/// \brief Copy constructor
-	template <typename PF>
-	constexpr point_type(const point_type<PF> &pt)
-		: m_x(static_cast<F>(pt.m_x))
-		, m_y(static_cast<F>(pt.m_y))
-		, m_z(static_cast<F>(pt.m_z))
-	{
-	}
-
-	/// \brief constructor taking a tuple of three values
-	constexpr point_type(const std::tuple<value_type, value_type, value_type> &pt)
-		: point_type(std::get<0>(pt), std::get<1>(pt), std::get<2>(pt))
-	{
-	}
-
-#if HAVE_LIBCLIPPER
-	/// \brief Construct a point using the values in clipper coordinate @a pt
-	constexpr point_type(const clipper::Coord_orth &pt)
-		: m_x(pt[0])
-		, m_y(pt[1])
-		, m_z(pt[2])
-	{
-	}
-
-	/// \brief Assign a point using the values in clipper coordinate @a rhs
-	constexpr point_type &operator=(const clipper::Coord_orth &rhs)
-	{
-		m_x = rhs[0];
-		m_y = rhs[1];
-		m_z = rhs[2];
-		return *this;
-	}
-#endif
-
-	/// \brief Assignment operator
-	template <typename PF>
-	constexpr point_type &operator=(const point_type<PF> &rhs)
-	{
-		m_x = static_cast<F>(rhs.m_x);
-		m_y = static_cast<F>(rhs.m_y);
-		m_z = static_cast<F>(rhs.m_z);
-		return *this;
-	}
-
-	[[nodiscard]] constexpr value_type &get_x() { return m_x; }      ///< Get a reference to x
-	[[nodiscard]] constexpr value_type get_x() const { return m_x; } ///< Get the value of x
-	constexpr void set_x(value_type x) { m_x = x; }                  ///< Set the value of x to @a x
-
-	[[nodiscard]] constexpr value_type &get_y() { return m_y; }      ///< Get a reference to y
-	[[nodiscard]] constexpr value_type get_y() const { return m_y; } ///< Get the value of y
-	constexpr void set_y(value_type y) { m_y = y; }                  ///< Set the value of y to @a y
-
-	[[nodiscard]] constexpr value_type &get_z() { return m_z; }      ///< Get a reference to z
-	[[nodiscard]] constexpr value_type get_z() const { return m_z; } ///< Get the value of z
-	constexpr void set_z(value_type z) { m_z = z; }                  ///< Set the value of z to @a z
-
-	/// \brief add @a rhs
-	constexpr point_type &operator+=(const point_type &rhs)
-	{
-		m_x += rhs.m_x;
-		m_y += rhs.m_y;
-		m_z += rhs.m_z;
-
-		return *this;
-	}
-
-	/// \brief add @a d to all members
-	constexpr point_type &operator+=(value_type d)
-	{
-		m_x += d;
-		m_y += d;
-		m_z += d;
-
-		return *this;
-	}
-
-	/// \brief Add the points @a lhs and @a rhs and return the result
-	template <typename F2>
-	friend constexpr auto operator+(const point_type &lhs, const point_type<F2> &rhs)
-	{
-		return point_type<std::common_type_t<value_type, F2>>(lhs.m_x + rhs.m_x, lhs.m_y + rhs.m_y, lhs.m_z + rhs.m_z);
-	}
-
-	/// \brief subtract @a rhs
-	constexpr point_type &operator-=(const point_type &rhs)
-	{
-		m_x -= rhs.m_x;
-		m_y -= rhs.m_y;
-		m_z -= rhs.m_z;
-
-		return *this;
-	}
-
-	/// \brief subtract @a d from all members
-	constexpr point_type &operator-=(value_type d)
-	{
-		m_x -= d;
-		m_y -= d;
-		m_z -= d;
-
-		return *this;
-	}
-
-	/// \brief Subtract the points @a lhs and @a rhs and return the result
-	template <typename F2>
-	friend constexpr auto operator-(const point_type &lhs, const point_type<F2> &rhs)
-	{
-		return point_type<std::common_type_t<value_type, F2>>(lhs.m_x - rhs.m_x, lhs.m_y - rhs.m_y, lhs.m_z - rhs.m_z);
-	}
-
-	/// \brief Return the negative copy of @a pt
-	friend constexpr point_type operator-(const point_type &pt)
-	{
-		return point_type(-pt.m_x, -pt.m_y, -pt.m_z);
-	}
-
-	/// \brief multiply all members with @a rhs
-	constexpr point_type &operator*=(value_type rhs)
-	{
-		m_x *= rhs;
-		m_y *= rhs;
-		m_z *= rhs;
-		return *this;
-	}
-
-	/// \brief multiply point @a pt with value @a f and return the result
-	template <typename F2>
-	friend constexpr auto operator*(const point_type &pt, F2 f)
-	{
-		return point_type<std::common_type_t<value_type, F2>>(pt.m_x * f, pt.m_y * f, pt.m_z * f);
-	}
-
-	/// \brief multiply point @a pt with value @a f and return the result
-	template <typename F2>
-	friend constexpr auto operator*(F2 f, const point_type &pt)
-	{
-		return point_type<std::common_type_t<value_type, F2>>(pt.m_x * f, pt.m_y * f, pt.m_z * f);
-	}
-
-	/// \brief divide all members by @a rhs
-	constexpr point_type &operator/=(value_type rhs)
-	{
-		m_x /= rhs;
-		m_y /= rhs;
-		m_z /= rhs;
-		return *this;
-	}
-
-	/// \brief divide point @a pt by value @a f and return the result
-	template <typename F2>
-	friend constexpr auto operator/(const point_type &pt, F2 f)
-	{
-		return point_type<std::common_type_t<value_type, F2>>(pt.m_x / f, pt.m_y / f, pt.m_z / f);
-	}
-
-	/**
-	 * @brief looking at this point as a vector, normalise it which
-	 * means dividing all members by the length making the length
-	 * effectively 1.
-	 *
-	 * @return The previous length of this vector
-	 */
-	constexpr value_type normalize()
-	{
-		auto length = m_x * m_x + m_y * m_y + m_z * m_z;
-		if (length > 0)
-		{
-			length = std::sqrt(length);
-			operator/=(length);
-		}
-		return length;
-	}
-
-	/// \brief Rotate this point using the quaterion @a q
-	constexpr void rotate(const quaternion &q)
-	{
-		quaternion_type<value_type> p(0, m_x, m_y, m_z);
-
-		p = q * p * conj(q);
-
-		m_x = p.get_b();
-		m_y = p.get_c();
-		m_z = p.get_d();
-	}
-
-	/// \brief Rotate this point using the quaterion @a q by first
-	/// moving the point to @a pivot and after rotating moving it
-	/// back
-	constexpr void rotate(const quaternion &q, point_type pivot)
-	{
-		operator-=(pivot);
-		rotate(q);
-		operator+=(pivot);
-	}
-
-#if HAVE_LIBCLIPPER
-	/// \brief Make it possible to pass a point to clipper functions expecting a clipper coordinate
-	operator clipper::Coord_orth() const
-	{
-		return clipper::Coord_orth(m_x, m_y, m_z);
-	}
-#endif
-
-	/// \brief Allow access to this point as if it is a tuple of three const value_type's
-	constexpr operator std::tuple<const value_type &, const value_type &, const value_type &>() const
-	{
-		return std::make_tuple(std::ref(m_x), std::ref(m_y), std::ref(m_z));
-	}
-
-	/// \brief Allow access to this point as if it is a tuple of three value_type's
-	constexpr operator std::tuple<value_type &, value_type &, value_type &>()
-	{
-		return std::make_tuple(std::ref(m_x), std::ref(m_y), std::ref(m_z));
-	}
-
-#if defined(__cpp_impl_three_way_comparison)
-	/// \brief a default spaceship operator
-	constexpr auto operator<=>(const point_type &rhs) const = default;
-#else
-	/// \brief a default equals operator
-	constexpr bool operator==(const point_type &rhs) const
-	{
-		return m_x == rhs.m_x and m_y == rhs.m_y and m_z == rhs.m_z;
-	}
-
-	/// \brief a default not-equals operator
-	constexpr bool operator!=(const point_type &rhs) const
-	{
-		return not operator==(rhs);
-	}
-#endif
-
-	// consider point as a vector... perhaps I should rename point?
-
-	/// \brief looking at the point as if it is a vector, return the squared length
-	[[nodiscard]] constexpr value_type length_sq() const
-	{
-		return m_x * m_x + m_y * m_y + m_z * m_z;
-	}
-
-	/// \brief looking at the point as if it is a vector, return the length
-	[[nodiscard]] constexpr value_type length() const
-	{
-		return std::sqrt(length_sq());
-	}
-
-	/// \brief Print out the point @a pt to @a os
-	friend std::ostream &operator<<(std::ostream &os, const point_type &pt)
-	{
-		os << '(' << pt.m_x << ',' << pt.m_y << ',' << pt.m_z << ')';
-		return os;
-	}
-};
-
-/// \brief By default we use points with float value_type
-using point = point_type<float>;
 
 // --------------------------------------------------------------------
 // several standard 3d operations
@@ -715,43 +90,16 @@ using point = point_type<float>;
 template <typename F1, typename F2>
 constexpr auto distance_squared(const point_type<F1> &a, const point_type<F2> &b)
 {
-	return (a.m_x - b.m_x) * (a.m_x - b.m_x) +
+	return (a.x - b.x) * (a.x - b.x) +
-	       (a.m_y - b.m_y) * (a.m_y - b.m_y) +
+	       (a.y - b.y) * (a.y - b.y) +
-	       (a.m_z - b.m_z) * (a.m_z - b.m_z);
+	       (a.z - b.z) * (a.z - b.z);
-}
-
-/// \brief return the distance between points @a a and @a b
-template <typename F1, typename F2>
-constexpr auto distance(const point_type<F1> &a, const point_type<F2> &b)
-{
-	return std::sqrt(
-		(a.m_x - b.m_x) * (a.m_x - b.m_x) +
-		(a.m_y - b.m_y) * (a.m_y - b.m_y) +
-		(a.m_z - b.m_z) * (a.m_z - b.m_z));
-}
-
-/// \brief return the dot product between the vectors @a a and @a b
-template <typename F1, typename F2>
-inline constexpr auto dot_product(const point_type<F1> &a, const point_type<F2> &b)
-{
-	return a.m_x * b.m_x + a.m_y * b.m_y + a.m_z * b.m_z;
-}
-
-/// \brief return the cross product between the vectors @a a and @a b
-template <typename F1, typename F2>
-inline constexpr auto cross_product(const point_type<F1> &a, const point_type<F2> &b)
-{
-	return point_type<std::common_type_t<F1, F2>>(
-		a.m_y * b.m_z - b.m_y * a.m_z,
-		a.m_z * b.m_x - b.m_z * a.m_x,
-		a.m_x * b.m_y - b.m_x * a.m_y);
 }
 
 /// \brief return the squared norm of point @a p
 template <typename F>
 constexpr F norm_squared(const point_type<F> &p)
 {
-	return p.m_x * p.m_x + p.m_y * p.m_y + p.m_z * p.m_z;
+	return p.x * p.x + p.y * p.y + p.z * p.z;
 }
 
 /// \brief return the norm of point @a p
@@ -763,23 +111,23 @@ constexpr point_type<F> norm(const point_type<F> &p)
 
 /// \brief return the point where two lines intersect, or an empty value if they don't intersect at all
 template <typename F>
-std::optional<cif::point> line_line_intersection(const point_type<F> &p1,
+std::optional<point> line_line_intersection(const point_type<F> &p1,
 	const point_type<F> &p2, const point_type<F> &p3, const point_type<F> &p4)
 {
 	auto p13 = p1 - p3;
 	auto p43 = p4 - p3;
-	if (std::abs(p43.m_x) < std::numeric_limits<F>::epsilon() and std::abs(p43.m_y) < std::numeric_limits<F>::epsilon() and std::abs(p43.m_z) < std::numeric_limits<F>::epsilon())
+	if (std::abs(p43.x) < std::numeric_limits<F>::epsilon() and std::abs(p43.y) < std::numeric_limits<F>::epsilon() and std::abs(p43.z) < std::numeric_limits<F>::epsilon())
 		return std::nullopt;
 
 	auto p21 = p2 - p1;
-	if (std::abs(p21.m_x) < std::numeric_limits<F>::epsilon() and std::abs(p21.m_y) < std::numeric_limits<F>::epsilon() and std::abs(p21.m_z) < std::numeric_limits<F>::epsilon())
+	if (std::abs(p21.x) < std::numeric_limits<F>::epsilon() and std::abs(p21.y) < std::numeric_limits<F>::epsilon() and std::abs(p21.z) < std::numeric_limits<F>::epsilon())
 		return std::nullopt;
 
-	auto d1343 = cif::dot_product(p43, p13);
+	auto d1343 = dot(p43, p13);
-	auto d4321 = cif::dot_product(p43, p21);
+	auto d4321 = dot(p43, p21);
-	auto d1321 = cif::dot_product(p13, p21);
+	auto d1321 = dot(p13, p21);
-	auto d4343 = cif::dot_product(p43, p43);
+	auto d4343 = dot(p43, p43);
-	auto d2121 = cif::dot_product(p21, p21);
+	auto d2121 = dot(p21, p21);
 
 	auto denom = d2121 * d4343 - d4321 * d4321;
 	if (std::abs(denom) < std::numeric_limits<F>::epsilon())
@@ -793,7 +141,7 @@ std::optional<cif::point> line_line_intersection(const point_type<F> &p1,
 	auto pa = p1 + mua * p21;
 	auto pb = p3 + mub * p43;
 
-	return { (pa + pb) / 2 };
+	return { (pa + pb) / 2.0f };
 }
 
 /// \brief return the angle in degrees between the vectors from point @a p2 to @a p1 and @a p2 to @a p3
@@ -803,7 +151,7 @@ constexpr auto angle(const point_type<F> &p1, const point_type<F> &p2, const poi
 	point_type<F> v1 = p1 - p2;
 	point_type<F> v2 = p3 - p2;
 
-	return std::acos(dot_product(v1, v2) / (v1.length() * v2.length())) * 180 / std::numbers::pi_v<F>;
+	return std::acos(glm::dot(v1, v2) / (glm::length(v1) * glm::length(v2))) * 180 / std::numbers::pi_v<F>;
 }
 
 /// \brief return the dihedral angle in degrees for the four points @a p1, @a p2, @a p3 and @a p4
@@ -817,18 +165,18 @@ constexpr auto dihedral_angle(const point_type<F> &p1, const point_type<F> &p2,
 
 	point_type<F> z = p2 - p3; // vector from p3 to p2
 
-	point_type<F> p = cross_product(z, v12);
+	point_type<F> p = glm::cross(z, v12);
-	point_type<F> x = cross_product(z, v43);
+	point_type<F> x = glm::cross(z, v43);
-	point_type<F> y = cross_product(z, x);
+	point_type<F> y = glm::cross(z, x);
 
-	auto u = dot_product(x, x);
+	auto u = glm::dot(x, x);
-	auto v = dot_product(y, y);
+	auto v = glm::dot(y, y);
 
 	F result = 360;
 	if (u > 0 and v > 0)
 	{
-		u = dot_product(p, x) / std::sqrt(u);
+		u = glm::dot(p, x) / std::sqrt(u);
-		v = dot_product(p, y) / std::sqrt(v);
+		v = glm::dot(p, y) / std::sqrt(v);
 		if (u != 0 or v != 0)
 			result = std::atan2(v, u) * static_cast<F>(180 / std::numbers::pi_v<F>);
 	}
@@ -843,9 +191,9 @@ constexpr auto cosinus_angle(const point_type<F> &p1, const point_type<F> &p2, c
 	point_type<F> v12 = p1 - p2;
 	point_type<F> v34 = p3 - p4;
 
-	auto x = dot_product(v12, v12) * dot_product(v34, v34);
+	auto x = glm::dot(v12, v12) * glm::dot(v34, v34);
 
-	return x > 0 ? dot_product(v12, v34) / std::sqrt(x) : 0;
+	return x > 0 ? glm::dot(v12, v34) / std::sqrt(x) : 0;
 }
 
 /// \brief return the distance from point @a p to the line from @a l1 to @a l2
@@ -855,28 +203,26 @@ constexpr auto distance_point_to_line(const point_type<F> &l1, const point_type<
 	auto line = l2 - l1;
 	auto p_to_l1 = p - l1;
 	auto p_to_l2 = p - l2;
-	auto cross = cross_product(p_to_l1, p_to_l2);
+	auto cross = glm::cross(p_to_l1, p_to_l2);
-	return cross.length() / line.length();
+	return glm::length(cross) / glm::length(line);
 }
 
 /// \brief return the smallest sphere around the points in @a pts
 std::tuple<point, float> smallest_sphere_around_points(std::vector<point> pts);
 
-// --------------------------------------------------------------------
-/**
- * @brief For e.g. simulated annealing, returns a new point that is moved in
- * a random direction with a distance randomly chosen from a normal
- * distribution with a stddev of offset.
- */
-point nudge(point p, float offset);
-
 // --------------------------------------------------------------------
 
 /// \brief Return a quaternion created from angle @a angle and axis @a axis
-quaternion construct_from_angle_axis(float angle, point axis);
+constexpr quaternion construct_from_angle_axis(float angle, point axis)
+{
+	return glm::angleAxis(glm::radians(angle), glm::normalize(axis));
+}
 
 /// \brief Return a tuple of an angle and an axis for quaternion @a q
-std::tuple<float, point> quaternion_to_angle_axis(quaternion q);
+constexpr std::tuple<float, point> quaternion_to_angle_axis(quaternion q)
+{
+	return { glm::degrees(glm::angle(q)), glm::axis(q) };
+}
 
 /// @brief Given four points and an angle, return the quaternion required to rotate
 /// point p4 along the p2-p3 axis and around point p3 to obtain the required within
@@ -902,74 +248,19 @@ quaternion align_points(const std::vector<point> &a, const std::vector<point> &b
 /// \brief The RMSd for the points in \a a and \a b
 double RMSd(const std::vector<point> &a, const std::vector<point> &b);
 
-// --------------------------------------------------------------------
-/**
- * @brief Helper class to generate evenly divided points on a sphere
- *
- * We use a fibonacci sphere to calculate even distribution of the dots
- *
- * @tparam N The number of points on the sphere is 2 * N + 1
- */
-template <int N>
-class spherical_dots
-{
-  public:
-	/// \brief the number of points
-	constexpr static int P = 2 * N * 1;
-
-	/// \brief the *weight* of the fibonacci sphere
-	constexpr static double W = (4 * std::numbers::pi) / P;
-
-	/// \brief the internal storage type
-	using array_type = typename std::array<point, P>;
-
-	/// \brief iterator type
-	using iterator = typename array_type::const_iterator;
-
-	/// \brief singleton instance
-	static spherical_dots &instance()
-	{
-		static spherical_dots sInstance;
-		return sInstance;
-	}
-
-	/// \brief The number of points
-	[[nodiscard]] std::size_t size() const { return P; }
-
-	/// \brief Access a point by index
-	const point operator[](uint32_t inIx) const { return m_points[inIx]; }
-
-	/// \brief iterator pointing to the first point
-	[[nodiscard]] iterator begin() const { return m_points.begin(); }
-
-	/// \brief iterator pointing past the last point
-	[[nodiscard]] iterator end() const { return m_points.end(); }
-
-	/// \brief return the *weight*,
-	[[nodiscard]] double weight() const { return W; }
-
-	spherical_dots()
-	{
-		const double
-			kGoldenRatio = std::numbers::phi;
-
-		auto p = m_points.begin();
-
-		for (int32_t i = -N; i <= N; ++i)
-		{
-			double lat = std::asin((2.0 * i) / P);
-			double lon = std::fmod(i, kGoldenRatio) * 2 * std::numbers::pi / kGoldenRatio;
-
-			p->m_x = std::sin(lon) * std::cos(lat);
-			p->m_y = std::cos(lon) * std::cos(lat);
-			p->m_z = std::sin(lat);
-
-			++p;
-		}
-	}
-
-  private:
-	array_type m_points;
-};
 
 } // namespace cif
+
+template <>
+struct std::formatter<glm::vec3> : std::formatter<std::string_view> // NOLINT
+{
+	template <class FmtContext>
+	auto format(glm::vec3 pt, FmtContext &ctx) const
+	{
+		std::string temp;
+		std::format_to(std::back_inserter(temp), "( {}, {}, {} )", pt.x, pt.y, pt.z);
+
+		return std::formatter<std::string_view>::format(temp, ctx);
+	}
+};
+

include/cif++/symmetry.hpp

@@ -27,11 +27,11 @@
 #pragma once
 
 #include "cif++/exports.hpp"
-#include "cif++/matrix.hpp"
 #include "cif++/point.hpp"
 
 #include <array>
 #include <cstdint>
+#include <glm/ext/matrix_float3x3.hpp>
 #include <string>
 
 #if defined(__cpp_impl_three_way_comparison)
@@ -49,18 +49,6 @@ namespace cif
 
 // --------------------------------------------------------------------
 
-/// \brief Apply matrix transformation @a m on point @a pt and return the result
-inline point operator*(const matrix3x3<float> &m, const point &pt)
-{
-	return {
-		m(0, 0) * pt.m_x + m(0, 1) * pt.m_y + m(0, 2) * pt.m_z,
-		m(1, 0) * pt.m_x + m(1, 1) * pt.m_y + m(1, 2) * pt.m_z,
-		m(2, 0) * pt.m_x + m(2, 1) * pt.m_y + m(2, 2) * pt.m_z
-	};
-}
-
-// --------------------------------------------------------------------
-
 /// \brief the space groups we know
 enum class space_group_name
 {
@@ -326,7 +314,7 @@ class transformation
 	transformation(const symop_data &data);
 
 	/// \brief constructor taking a rotation matrix @a r and a translation vector @a t
-	transformation(const matrix3x3<float> &r, const cif::point &t);
+	transformation(const glm::mat3 &r, const cif::point &t);
 
 	/** @cond */
 	transformation(const transformation &) = default;
@@ -338,8 +326,8 @@ class transformation
 	/// \brief operator() to perform the transformation on point @a pt and return the result
 	point operator()(point pt) const
 	{
-		if (m_q)
+		if (m_q != quaternion{})
-			pt.rotate(m_q);
+			pt = m_q * pt;
 		else
 			pt = m_rotation * pt;
 
@@ -367,9 +355,9 @@ class transformation
 
 	void try_create_quaternion();
 
-	matrix3x3<float> m_rotation;
+	glm::mat3 m_rotation{};
-	quaternion m_q;
+	quaternion m_q{};
-	point m_translation;
+	point m_translation{};
 };
 
 // --------------------------------------------------------------------
@@ -399,14 +387,14 @@ class cell
 
 	[[nodiscard]] float get_volume() const; ///< return the calculated volume for this cell
 
-	[[nodiscard]] matrix3x3<float> get_orthogonal_matrix() const { return m_orthogonal; } ///< return the matrix to use to transform coordinates from fractional to orthogonal
+	[[nodiscard]] glm::mat3 get_orthogonal_matrix() const { return m_orthogonal; } ///< return the matrix to use to transform coordinates from fractional to orthogonal
-	[[nodiscard]] matrix3x3<float> get_fractional_matrix() const { return m_fractional; } ///< return the matrix to use to transform coordinates from orthogonal to fractional
+	[[nodiscard]] glm::mat3 get_fractional_matrix() const { return m_fractional; } ///< return the matrix to use to transform coordinates from orthogonal to fractional
 
   private:
 	void init();
 
 	float m_a, m_b, m_c, m_alpha, m_beta, m_gamma;
-	matrix3x3<float> m_orthogonal, m_fractional;
+	glm::mat3 m_orthogonal, m_fractional;
 };
 
 // --------------------------------------------------------------------

src/model.cpp

@@ -67,9 +67,9 @@ void atom::atom_impl::moveTo(const point &p)
 
 	auto r = row();
 
-	r.assign("Cartn_x", { p.m_x, 3 }, false, false);
+	r.assign("Cartn_x", { p.x, 3 }, false, false);
-	r.assign("Cartn_y", { p.m_y, 3 }, false, false);
+	r.assign("Cartn_y", { p.y, 3 }, false, false);
-	r.assign("Cartn_z", { p.m_z, 3 }, false, false);
+	r.assign("Cartn_z", { p.z, 3 }, false, false);
 
 	m_location = p;
 }
@@ -689,8 +689,8 @@ float monomer::chiral_volume() const
 		auto atom2 = get_atom_by_atom_id("CD1");
 		auto atom3 = get_atom_by_atom_id("CD2");
 
-		result = dot_product(atom1.get_location() - centre.get_location(),
+		result = dot(atom1.get_location() - centre.get_location(),
-			cross_product(atom2.get_location() - centre.get_location(), atom3.get_location() - centre.get_location()));
+			cross(atom2.get_location() - centre.get_location(), atom3.get_location() - centre.get_location()));
 	}
 	else if (m_compound_id == "VAL")
 	{
@@ -699,8 +699,8 @@ float monomer::chiral_volume() const
 		auto atom2 = get_atom_by_atom_id("CG1");
 		auto atom3 = get_atom_by_atom_id("CG2");
 
-		result = dot_product(atom1.get_location() - centre.get_location(),
+		result = dot(atom1.get_location() - centre.get_location(),
-			cross_product(atom2.get_location() - centre.get_location(), atom3.get_location() - centre.get_location()));
+			cross(atom2.get_location() - centre.get_location(), atom3.get_location() - centre.get_location()));
 	}
 
 	return result;
@@ -2288,17 +2288,15 @@ std::string structure::create_non_poly(const std::string &compound_id, bool skip
 		if (skip_hydrogen and cif::atom_type_traits(a.type_symbol).symbol() == "H")
 			continue;
 
-		auto ax = a.get_location().get_x();
+		auto aloc = a.get_location();
-		auto ay = a.get_location().get_y();
-		auto az = a.get_location().get_z();
 
 		atoms.emplace_back(cif::row_initializer{
 			{ "type_symbol", cif::atom_type_traits(a.type_symbol).symbol() },
 			{ "label_atom_id", a.id },
 			{ "auth_atom_id", a.id },
-			{ "Cartn_x", ax },
+			{ "Cartn_x", aloc.x },
-			{ "Cartn_y", ay },
+			{ "Cartn_y", aloc.y },
-			{ "Cartn_z", az },
+			{ "Cartn_z", aloc.z },
 			{ "B_iso_or_equiv", 30.00 } });
 	}
 

src/point.cpp

@@ -34,8 +34,9 @@
 #include <cstdint>
 #include <cstdlib>
 #include <initializer_list>
-#include <numbers>
 #include <optional>
+#include <glm/ext/matrix_double3x3.hpp>
+#include <glm/ext/matrix_double4x4.hpp>
 #include <random>
 #include <stdexcept>
 #include <tuple>
@@ -47,84 +48,26 @@ namespace cif
 
 // --------------------------------------------------------------------
 
-template <typename T>
-quaternion_type<T> normalize(quaternion_type<T> q)
-{
-	std::valarray<double> t(4);
-
-	t[0] = q.get_a();
-	t[1] = q.get_b();
-	t[2] = q.get_c();
-	t[3] = q.get_d();
-
-	t *= t;
-
-	double length = std::sqrt(t.sum());
-
-	if (length > 0.001)
-		q /= static_cast<quaternion::value_type>(length);
-	else
-		q = quaternion(1, 0, 0, 0);
-
-	return q;
-}
-
-// --------------------------------------------------------------------
-
-quaternion construct_from_angle_axis(float angle, point axis)
-{
-	angle = (angle * std::numbers::pi_v<float> / 180) / 2;
-	auto s = std::sin(angle);
-	auto c = std::cos(angle);
-
-	axis.normalize();
-
-	return normalize(quaternion{
-		static_cast<float>(c),
-		static_cast<float>(s * axis.m_x),
-		static_cast<float>(s * axis.m_y),
-		static_cast<float>(s * axis.m_z) });
-}
-
-std::tuple<float, point> quaternion_to_angle_axis(quaternion q)
-{
-	if (q.get_a() > 1)
-		q = normalize(q);
-
-	// angle:
-	float angle = 2 * std::acos(q.get_a());
-	angle = angle * 180 / std::numbers::pi_v<float>;
-
-	// axis:
-	float s = std::sqrt(1 - q.get_a() * q.get_a());
-	if (s < 0.001)
-		s = 1;
-
-	point axis(q.get_b() / s, q.get_c() / s, q.get_d() / s);
-
-	return { angle, axis };
-}
-
 point center_points(std::vector<point> &Points)
 {
-	point t;
+	point t{};
 
 	for (point &pt : Points)
 	{
-		t.m_x += pt.m_x;
+		t.x += pt.x;
-		t.m_y += pt.m_y;
+		t.y += pt.y;
-		t.m_z += pt.m_z;
+		t.z += pt.z;
 	}
 
-	t.m_x /= static_cast<float>(Points.size());
+	t.x /= static_cast<float>(Points.size());
-	t.m_y /= static_cast<float>(Points.size());
+	t.y /= static_cast<float>(Points.size());
-	t.m_z /= static_cast<float>(Points.size());
+	t.z /= static_cast<float>(Points.size());
 
 	for (point &pt : Points)
 	{
-		pt.m_x -= t.m_x;
+		pt.x -= t.x;
-		pt.m_y -= t.m_y;
+		pt.y -= t.y;
-		pt.m_z -= t.m_z;
+		pt.z -= t.z;
 	}
 
 	return t;
@@ -163,9 +106,9 @@ double RMSd(const std::vector<point> &a, const std::vector<point> &b)
 	{
 		std::valarray<double> d(3);
 
-		d[0] = b[i].m_x - a[i].m_x;
+		d[0] = b[i].x - a[i].x;
-		d[1] = b[i].m_y - a[i].m_y;
+		d[1] = b[i].y - a[i].y;
-		d[2] = b[i].m_z - a[i].m_z;
+		d[2] = b[i].z - a[i].z;
 
 		d *= d;
 
@@ -214,43 +157,82 @@ double LargestDepressedQuarticSolution(double a, double b, double c)
 	return t.max();
 }
 
+/**
+ * @brief Implementation of a cofactor calculation as a matrix expression
+ *
+ * @tparam M Type of matrix
+ */
+
+auto cofactors(glm::dmat4 &m)
+{
+	glm::dmat4 result{};
+
+	const std::size_t ixs[4][3] = {
+		{ 1, 2, 3 },
+		{ 0, 2, 3 },
+		{ 0, 1, 3 },
+		{ 0, 1, 2 }
+	};
+
+	for (size_t i = 0; i < 4; ++i)
+	{
+		for (size_t j = 0; j < 4; ++j)
+		{
+			const std::size_t *ix = ixs[i];
+			const std::size_t *iy = ixs[j];
+
+			auto e =
+				m[ix[0]][iy[0]] * m[ix[1]][iy[1]] * m[ix[2]][iy[2]] +
+				m[ix[0]][iy[1]] * m[ix[1]][iy[2]] * m[ix[2]][iy[0]] +
+				m[ix[0]][iy[2]] * m[ix[1]][iy[0]] * m[ix[2]][iy[1]] -
+				m[ix[0]][iy[2]] * m[ix[1]][iy[1]] * m[ix[2]][iy[0]] -
+				m[ix[0]][iy[1]] * m[ix[1]][iy[0]] * m[ix[2]][iy[2]] -
+				m[ix[0]][iy[0]] * m[ix[1]][iy[2]] * m[ix[2]][iy[1]];
+
+			result[i][j] = (i + j) % 2 == 1 ? -e : e;
+		}
+	}
+
+	return result;
+}
+
 quaternion align_points(const std::vector<point> &pa, const std::vector<point> &pb)
 {
 	// First calculate M, a 3x3 matrix containing the sums of products of the coordinates of A and B
-	matrix3x3<double> M;
+	glm::dmat3x3 M{};
 
 	for (uint32_t i = 0; i < pa.size(); ++i)
 	{
 		const point &a = pa[i];
 		const point &b = pb[i];
 
-		M(0, 0) += a.m_x * b.m_x;
+		M[0][0] += a.x * b.x;
-		M(0, 1) += a.m_x * b.m_y;
+		M[0][1] += a.x * b.y;
-		M(0, 2) += a.m_x * b.m_z;
+		M[0][2] += a.x * b.z;
-		M(1, 0) += a.m_y * b.m_x;
+		M[1][0] += a.y * b.x;
-		M(1, 1) += a.m_y * b.m_y;
+		M[1][1] += a.y * b.y;
-		M(1, 2) += a.m_y * b.m_z;
+		M[1][2] += a.y * b.z;
-		M(2, 0) += a.m_z * b.m_x;
+		M[2][0] += a.z * b.x;
-		M(2, 1) += a.m_z * b.m_y;
+		M[2][1] += a.z * b.y;
-		M(2, 2) += a.m_z * b.m_z;
+		M[2][2] += a.z * b.z;
 	}
 
 	// Now calculate N, a symmetric 4x4 matrix
-	symmetric_matrix4x4<double> N(4);
+	glm::dmat4x4 N;
 
-	N(0, 0) = M(0, 0) + M(1, 1) + M(2, 2);
+	N[0][0] = M[0][0] + M[1][1] + M[2][2];
-	N(0, 1) = M(1, 2) - M(2, 1);
+	N[0][1] = N[1][0] = M[1][2] - M[2][1];
-	N(0, 2) = M(2, 0) - M(0, 2);
+	N[0][2] = N[2][0] = M[2][0] - M[0][2];
-	N(0, 3) = M(0, 1) - M(1, 0);
+	N[0][3] = N[3][0] = M[0][1] - M[1][0];
 
-	N(1, 1) = M(0, 0) - M(1, 1) - M(2, 2);
+	N[1][1] = M[0][0] - M[1][1] - M[2][2];
-	N(1, 2) = M(0, 1) + M(1, 0);
+	N[1][2] = N[2][1] = M[0][1] + M[1][0];
-	N(1, 3) = M(0, 2) + M(2, 0);
+	N[1][3] = N[3][1] = M[0][2] + M[2][0];
 
-	N(2, 2) = -M(0, 0) + M(1, 1) - M(2, 2);
+	N[2][2] = -M[0][0] + M[1][1] - M[2][2];
-	N(2, 3) = M(1, 2) + M(2, 1);
+	N[2][3] = N[3][2] = M[1][2] + M[2][1];
 
-	N(3, 3) = -M(0, 0) - M(1, 1) + M(2, 2);
+	N[3][3] = -M[0][0] - M[1][1] + M[2][2];
 
 	// det(N - λI) = 0
 	// find the largest λ (λm)
@@ -261,41 +243,41 @@ quaternion align_points(const std::vector<point> &pa, const std::vector<point> &
 	// and so this is a so-called depressed quartic
 	// solve it using Ferrari's algorithm
 
-	double C = -2 * (M(0, 0) * M(0, 0) + M(0, 1) * M(0, 1) + M(0, 2) * M(0, 2) +
+	double C = -2 * (M[0][0] * M[0][0] + M[0][1] * M[0][1] + M[0][2] * M[0][2] +
-						M(1, 0) * M(1, 0) + M(1, 1) * M(1, 1) + M(1, 2) * M(1, 2) +
+						M[1][0] * M[1][0] + M[1][1] * M[1][1] + M[1][2] * M[1][2] +
-						M(2, 0) * M(2, 0) + M(2, 1) * M(2, 1) + M(2, 2) * M(2, 2));
+						M[2][0] * M[2][0] + M[2][1] * M[2][1] + M[2][2] * M[2][2]);
 
-	double D = 8 * (M(0, 0) * M(1, 2) * M(2, 1) +
+	double D = 8 * (M[0][0] * M[1][2] * M[2][1] +
-					   M(1, 1) * M(2, 0) * M(0, 2) +
+					   M[1][1] * M[2][0] * M[0][2] +
-					   M(2, 2) * M(0, 1) * M(1, 0)) -
+					   M[2][2] * M[0][1] * M[1][0]) -
-	           8 * (M(0, 0) * M(1, 1) * M(2, 2) +
+	           8 * (M[0][0] * M[1][1] * M[2][2] +
-					   M(1, 2) * M(2, 0) * M(0, 1) +
+					   M[1][2] * M[2][0] * M[0][1] +
-					   M(2, 1) * M(1, 0) * M(0, 2));
+					   M[2][1] * M[1][0] * M[0][2]);
 
 	// E is the determinant of N:
 	double E =
-		(N(0, 0) * N(1, 1) - N(0, 1) * N(0, 1)) * (N(2, 2) * N(3, 3) - N(2, 3) * N(2, 3)) +
+		(N[0][0] * N[1][1] - N[0][1] * N[0][1]) * (N[2][2] * N[3][3] - N[2][3] * N[2][3]) +
-		(N(0, 1) * N(0, 2) - N(0, 0) * N(2, 1)) * (N(2, 1) * N(3, 3) - N(2, 3) * N(1, 3)) +
+		(N[0][1] * N[0][2] - N[0][0] * N[2][1]) * (N[2][1] * N[3][3] - N[2][3] * N[1][3]) +
-		(N(0, 0) * N(1, 3) - N(0, 1) * N(0, 3)) * (N(2, 1) * N(2, 3) - N(2, 2) * N(1, 3)) +
+		(N[0][0] * N[1][3] - N[0][1] * N[0][3]) * (N[2][1] * N[2][3] - N[2][2] * N[1][3]) +
-		(N(0, 1) * N(2, 1) - N(1, 1) * N(0, 2)) * (N(0, 2) * N(3, 3) - N(2, 3) * N(0, 3)) +
+		(N[0][1] * N[2][1] - N[1][1] * N[0][2]) * (N[0][2] * N[3][3] - N[2][3] * N[0][3]) +
-		(N(1, 1) * N(0, 3) - N(0, 1) * N(1, 3)) * (N(0, 2) * N(2, 3) - N(2, 2) * N(0, 3)) +
+		(N[1][1] * N[0][3] - N[0][1] * N[1][3]) * (N[0][2] * N[2][3] - N[2][2] * N[0][3]) +
-		(N(0, 2) * N(1, 3) - N(2, 1) * N(0, 3)) * (N(0, 2) * N(1, 3) - N(2, 1) * N(0, 3));
+		(N[0][2] * N[1][3] - N[2][1] * N[0][3]) * (N[0][2] * N[1][3] - N[2][1] * N[0][3]);
 
 	// solve quartic
 	double lambda = LargestDepressedQuarticSolution(C, D, E);
 
 	// calculate t = (N - λI)
-	matrix<double> t(N - identity_matrix(4) * lambda);
+	auto t = N - glm::dmat4x4(1.0) * lambda;
 
 	// calculate a matrix of cofactors for t
-	auto cf = matrix_cofactors(t);
+	auto cf = cofactors(t);
 
 	int maxR = 0;
-	double maxCF = std::abs(cf(0, 0));
+	double maxCF = std::abs(cf[0][0]);
 
 	for (int r = 1; r < 4; ++r)
 	{
-		auto cfr = std::abs(cf(r, 0));
+		auto cfr = std::abs(cf[r][0]);
 		if (maxCF < cfr)
 		{
 			maxCF = cfr;
@@ -304,10 +286,10 @@ quaternion align_points(const std::vector<point> &pa, const std::vector<point> &
 	}
 
 	quaternion q(
-		static_cast<float>(cf(maxR, 0)),
+		static_cast<float>(cf[maxR][0]),
-		static_cast<float>(cf(maxR, 1)),
+		static_cast<float>(cf[maxR][1]),
-		static_cast<float>(cf(maxR, 2)),
+		static_cast<float>(cf[maxR][2]),
-		static_cast<float>(cf(maxR, 3)));
+		static_cast<float>(cf[maxR][3]));
 	q = normalize(q);
 
 	return q;
@@ -315,52 +297,29 @@ quaternion align_points(const std::vector<point> &pa, const std::vector<point> &
 
 // --------------------------------------------------------------------
 
-point nudge(point p, float offset)
+std::tuple<point, float> smallest_sphere_around_2_points(std::array<point, 2> pts)
 {
-	static std::random_device rd;
+	return { (pts[0] + pts[1]) / 2.0f, distance(pts[0], pts[1]) / 2.0f };
-	static std::mt19937_64 rng(rd());
-
-	std::uniform_real_distribution<float> randomAngle(0, 2 * std::numbers::pi);
-	std::normal_distribution<float> randomOffset(0, offset);
-
-	float theta = randomAngle(rng);
-	float phi1 = randomAngle(rng) - static_cast<float>(std::numbers::pi);
-	float phi2 = randomAngle(rng) - static_cast<float>(std::numbers::pi);
-
-	quaternion q = spherical(1.0f, theta, phi1, phi2);
-
-	point r{ 0, 0, 1 };
-	r.rotate(q);
-	r *= randomOffset(rng);
-
-	return p + r;
 }
 
-// --------------------------------------------------------------------
+std::tuple<point, float> smallest_sphere_around_3_points(std::array<point, 3> pts)
-
-std::tuple<point, float> smallest_sphere_around_2_points(std::array<cif::point, 2> pts)
-{
-	return { (pts[0] + pts[1]) / 2, distance(pts[0], pts[1]) / 2 };
-}
-
-std::tuple<point, float> smallest_sphere_around_3_points(std::array<cif::point, 3> pts)
 {
 	// Find two bisectors
-	auto vz = cross_product(pts[1] - pts[0], pts[2] - pts[0]);
+	auto vz = cross(pts[1] - pts[0], pts[2] - pts[0]);
 
-	auto bs1 = cross_product(vz, pts[1] - pts[0]);
+	auto bs1 = cross(vz, pts[1] - pts[0]);
-	bs1.normalize();
+	bs1 = glm::normalize(bs1);
 
 	auto v1 = (pts[1] - pts[0]);
-	v1.normalize();
+	v1 = glm::normalize(v1);
 
 	auto s1 = pts[0] + (distance(pts[1], pts[0]) / 2) * v1;
 
-	auto bs2 = cross_product(vz, pts[2] - pts[0]);
+	auto bs2 = cross(vz, pts[2] - pts[0]);
-	bs2.normalize();
+	bs2 = glm::normalize(bs2);
 
 	auto v2 = (pts[2] - pts[0]);
-	v2.normalize();
+	v2 = glm::normalize(v2);
 
 	auto s2 = pts[0] + (distance(pts[2], pts[0]) / 2) * v2;
 
@@ -381,7 +340,7 @@ std::tuple<point, float> smallest_sphere_around_3_points(std::array<cif::point,
 		return smallest_sphere_around_2_points({ pts[1], pts[2] });
 }
 
-std::tuple<point, float> smallest_sphere_around_4_points(std::array<cif::point, 4> pts)
+std::tuple<point, float> smallest_sphere_around_4_points(std::array<point, 4> pts)
 {
 	auto t0 = -norm_squared(pts[0]);
 	auto t1 = -norm_squared(pts[1]);
@@ -389,46 +348,46 @@ std::tuple<point, float> smallest_sphere_around_4_points(std::array<cif::point,
 	auto t3 = -norm_squared(pts[3]);
 
 	// clang-format off
-	matrix4x4<float> Tm({
+	glm::mat4x4 Tm({
-		pts[0].m_x, pts[0].m_y, pts[0].m_z, 1,
+		pts[0].x, pts[0].y, pts[0].z, 1,
-		pts[1].m_x, pts[1].m_y, pts[1].m_z, 1,
+		pts[1].x, pts[1].y, pts[1].z, 1,
-		pts[2].m_x, pts[2].m_y, pts[2].m_z, 1,
+		pts[2].x, pts[2].y, pts[2].z, 1,
-		pts[3].m_x, pts[3].m_y, pts[3].m_z, 1
+		pts[3].x, pts[3].y, pts[3].z, 1
 	});
 	auto T = determinant(Tm);
 
 	if (T != 0)
 	{
-		matrix4x4<float> Dm({
+		glm::mat4x4 Dm({
-			t0, pts[0].m_y, pts[0].m_z, 1,
+			t0, pts[0].y, pts[0].z, 1,
-			t1, pts[1].m_y, pts[1].m_z, 1,
+			t1, pts[1].y, pts[1].z, 1,
-			t2, pts[2].m_y, pts[2].m_z, 1,
+			t2, pts[2].y, pts[2].z, 1,
-			t3, pts[3].m_y, pts[3].m_z, 1
+			t3, pts[3].y, pts[3].z, 1
 		});
 		auto D = determinant(Dm) / T;
 		
-		matrix4x4<float> Em({
+		glm::mat4x4 Em({
-			pts[0].m_x, t0, pts[0].m_z, 1,
+			pts[0].x, t0, pts[0].z, 1,
-			pts[1].m_x, t1, pts[1].m_z, 1,
+			pts[1].x, t1, pts[1].z, 1,
-			pts[2].m_x, t2, pts[2].m_z, 1,
+			pts[2].x, t2, pts[2].z, 1,
-			pts[3].m_x, t3, pts[3].m_z, 1
+			pts[3].x, t3, pts[3].z, 1
 		});
 		auto E = determinant(Em) / T;
 		
-		matrix4x4<float> Fm({
+		glm::mat4x4 Fm({
-			pts[0].m_x, pts[0].m_y, t0, 1,
+			pts[0].x, pts[0].y, t0, 1,
-			pts[1].m_x, pts[1].m_y, t1, 1,
+			pts[1].x, pts[1].y, t1, 1,
-			pts[2].m_x, pts[2].m_y, t2, 1,
+			pts[2].x, pts[2].y, t2, 1,
-			pts[3].m_x, pts[3].m_y, t3, 1
+			pts[3].x, pts[3].y, t3, 1
 		});
 		
 		auto F = determinant(Fm) / T;
 		
-		matrix4x4<float> Gm({
+		glm::mat4x4 Gm({
-			pts[0].m_x, pts[0].m_y, pts[0].m_z, t0,
+			pts[0].x, pts[0].y, pts[0].z, t0,
-			pts[1].m_x, pts[1].m_y, pts[1].m_z, t1,
+			pts[1].x, pts[1].y, pts[1].z, t1,
-			pts[2].m_x, pts[2].m_y, pts[2].m_z, t2,
+			pts[2].x, pts[2].y, pts[2].z, t2,
-			pts[3].m_x, pts[3].m_y, pts[3].m_z, t3
+			pts[3].x, pts[3].y, pts[3].z, t3
 		});
 		auto G = determinant(Gm) / T;
 		
@@ -508,19 +467,19 @@ bool point_in_circle(point p, std::vector<point> c)
 		case 2:
 		{
 			auto [center, radius] = smallest_sphere_around_2_points({ c[0], c[1] });
-			return cif::distance_squared(p, center) <= radius * radius;
+			return distance_squared(p, center) <= radius * radius;
 		}
 
 		case 3:
 		{
 			auto [center, radius] = smallest_sphere_around_3_points({ c[0], c[1], c[2] });
-			return cif::distance_squared(p, center) <= radius * radius;
+			return distance_squared(p, center) <= radius * radius;
 		}
 
 		case 4:
 		{
 			auto [center, radius] = smallest_sphere_around_4_points({ c[0], c[1], c[2], c[3] });
-			return cif::distance_squared(p, center) <= radius * radius;
+			return distance_squared(p, center) <= radius * radius;
 		}
 
 		default:

src/symmetry.cpp

@@ -81,14 +81,15 @@ void cell::init()
 
 	auto alpha_star = std::acos((std::cos(gamma) * std::cos(beta) - std::cos(alpha)) / (std::sin(beta) * std::sin(gamma)));
 
-	m_orthogonal = identity_matrix(3);
+	m_orthogonal = glm::mat3(1.0f);
 
-	m_orthogonal(0, 0) = m_a;
+	// WARNING: glm matrices are column major, by default
-	m_orthogonal(0, 1) = m_b * std::cos(gamma);
+	m_orthogonal[0][0] = m_a;
-	m_orthogonal(0, 2) = m_c * std::cos(beta);
+	m_orthogonal[1][0] = m_b * std::cos(gamma);
-	m_orthogonal(1, 1) = m_b * std::sin(gamma);
+	m_orthogonal[2][0] = m_c * std::cos(beta);
-	m_orthogonal(1, 2) = m_c * std::sin(beta) * std::cos(alpha_star);
+	m_orthogonal[1][1] = m_b * std::sin(gamma);
-	m_orthogonal(2, 2) = m_c * std::sin(beta) * std::sin(alpha_star);
+	m_orthogonal[2][1] = m_c * std::sin(beta) * std::cos(alpha_star);
+	m_orthogonal[2][2] = m_c * std::sin(beta) * std::sin(alpha_star);
 
 	m_fractional = inverse(m_orthogonal);
 }
@@ -150,24 +151,25 @@ transformation::transformation(const symop_data &data)
 {
 	const auto &d = data.data();
 
-	m_rotation(0, 0) = static_cast<float>(d[0]);
+	// WARNING: glm::mat is column major
-	m_rotation(0, 1) = static_cast<float>(d[1]);
+	m_rotation[0][0] = static_cast<float>(d[0]);
-	m_rotation(0, 2) = static_cast<float>(d[2]);
+	m_rotation[1][0] = static_cast<float>(d[1]);
-	m_rotation(1, 0) = static_cast<float>(d[3]);
+	m_rotation[2][0] = static_cast<float>(d[2]);
-	m_rotation(1, 1) = static_cast<float>(d[4]);
+	m_rotation[0][1] = static_cast<float>(d[3]);
-	m_rotation(1, 2) = static_cast<float>(d[5]);
+	m_rotation[1][1] = static_cast<float>(d[4]);
-	m_rotation(2, 0) = static_cast<float>(d[6]);
+	m_rotation[2][1] = static_cast<float>(d[5]);
-	m_rotation(2, 1) = static_cast<float>(d[7]);
+	m_rotation[0][2] = static_cast<float>(d[6]);
-	m_rotation(2, 2) = static_cast<float>(d[8]);
+	m_rotation[1][2] = static_cast<float>(d[7]);
+	m_rotation[2][2] = static_cast<float>(d[8]);
 
 	try_create_quaternion();
 
-	m_translation.m_x = static_cast<float>(d[9] == 0 ? 0 : 1.0 * d[9] / d[10]);
+	m_translation.x = static_cast<float>(d[9] == 0 ? 0 : 1.0 * d[9] / d[10]);
-	m_translation.m_y = static_cast<float>(d[11] == 0 ? 0 : 1.0 * d[11] / d[12]);
+	m_translation.y = static_cast<float>(d[11] == 0 ? 0 : 1.0 * d[11] / d[12]);
-	m_translation.m_z = static_cast<float>(d[13] == 0 ? 0 : 1.0 * d[13] / d[14]);
+	m_translation.z = static_cast<float>(d[13] == 0 ? 0 : 1.0 * d[13] / d[14]);
 }
 
-transformation::transformation(const matrix3x3<float> &r, const cif::point &t)
+transformation::transformation(const glm::mat3 &r, const cif::point &t)
 	: m_rotation(r)
 	, m_translation(t)
 {
@@ -178,9 +180,9 @@ void transformation::try_create_quaternion()
 {
 	Eigen::Matrix3f rot;
 
-	rot << m_rotation(0, 0), m_rotation(0, 1), m_rotation(0, 2),
+	rot << m_rotation[0][0], m_rotation[0][1], m_rotation[0][2],
-		m_rotation(1, 0), m_rotation(1, 1), m_rotation(1, 2),
+		m_rotation[1][0], m_rotation[1][1], m_rotation[1][2],
-		m_rotation(2, 0), m_rotation(2, 1), m_rotation(2, 2);
+		m_rotation[2][0], m_rotation[2][1], m_rotation[2][2];
 
 	if (rot * rot.transpose() == Eigen::Matrix3f::Identity() and rot.determinant() == 1)
 	{
@@ -241,20 +243,20 @@ point offsetToOrigin(const cell &c, const point &p)
 {
 	point d{};
 
-	while (p.m_x + d.m_x < -(c.get_a()))
+	while (p.x + d.x < -(c.get_a()))
-		d.m_x += c.get_a();
+		d.x += c.get_a();
-	while (p.m_x + d.m_x > (c.get_a()))
+	while (p.x + d.x > (c.get_a()))
-		d.m_x -= c.get_a();
+		d.x -= c.get_a();
 
-	while (p.m_y + d.m_y < -(c.get_b()))
+	while (p.y + d.y < -(c.get_b()))
-		d.m_y += c.get_b();
+		d.y += c.get_b();
-	while (p.m_y + d.m_y > (c.get_b()))
+	while (p.y + d.y > (c.get_b()))
-		d.m_y -= c.get_b();
+		d.y -= c.get_b();
 
-	while (p.m_z + d.m_z < -(c.get_c()))
+	while (p.z + d.z < -(c.get_c()))
-		d.m_z += c.get_c();
+		d.z += c.get_c();
-	while (p.m_z + d.m_z > (c.get_c()))
+	while (p.z + d.z > (c.get_c()))
-		d.m_z -= c.get_c();
+		d.z -= c.get_c();
 
 	return d;
 };
@@ -263,20 +265,20 @@ point offsetToOriginFractional(const point &p)
 {
 	point d{};
 
-	while (p.m_x + d.m_x < -0.5f)
+	while (p.x + d.x < -0.5f)
-		d.m_x += 1;
+		d.x += 1;
-	while (p.m_x + d.m_x > 0.5f)
+	while (p.x + d.x > 0.5f)
-		d.m_x -= 1;
+		d.x -= 1;
 
-	while (p.m_y + d.m_y < -0.5f)
+	while (p.y + d.y < -0.5f)
-		d.m_y += 1;
+		d.y += 1;
-	while (p.m_y + d.m_y > 0.5f)
+	while (p.y + d.y > 0.5f)
-		d.m_y -= 1;
+		d.y -= 1;
 
-	while (p.m_z + d.m_z < -0.5f)
+	while (p.z + d.z < -0.5f)
-		d.m_z += 1;
+		d.z += 1;
-	while (p.m_z + d.m_z > 0.5f)
+	while (p.z + d.z > 0.5f)
-		d.m_z -= 1;
+		d.z -= 1;
 
 	return d;
 };
@@ -288,16 +290,11 @@ point spacegroup::operator()(const point &pt, const cell &c, sym_op symop) const
 
 	transformation t = at(symop.m_nr - 1);
 
-	t.m_translation.m_x += symop.m_ta - 5;
+	t.m_translation.x += symop.m_ta - 5;
-	t.m_translation.m_y += symop.m_tb - 5;
+	t.m_translation.y += symop.m_tb - 5;
-	t.m_translation.m_z += symop.m_tc - 5;
+	t.m_translation.z += symop.m_tc - 5;
 
-	auto fpt = fractional(pt, c);
+	return orthogonal(t(fractional(pt, c)), c);
-	auto o = offsetToOriginFractional(fpt);
-
-	auto spt = t(fpt + o) - o;
-
-	return orthogonal(spt, c);
 }
 
 point spacegroup::inverse(const point &pt, const cell &c, sym_op symop) const
@@ -307,17 +304,12 @@ point spacegroup::inverse(const point &pt, const cell &c, sym_op symop) const
 
 	transformation t = at(symop.m_nr - 1);
 
-	t.m_translation.m_x += symop.m_ta - 5;
+	t.m_translation.x += symop.m_ta - 5;
-	t.m_translation.m_y += symop.m_tb - 5;
+	t.m_translation.y += symop.m_tb - 5;
-	t.m_translation.m_z += symop.m_tc - 5;
+	t.m_translation.z += symop.m_tc - 5;
-
-	auto fpt = fractional(pt, c);
-	auto o = offsetToOriginFractional(fpt);
 
 	auto it = cif::inverse(t);
-	auto spt = it(fpt + o) - o;
+	return orthogonal(it(fractional(pt, c)), c);
-
-	return orthogonal(spt, c);
 }
 
 // --------------------------------------------------------------------
@@ -448,7 +440,7 @@ std::tuple<float, point, sym_op> crystal::closest_symmetry_copy(point a, point b
 	if (m_cell.get_a() == 0 or m_cell.get_b() == 0 or m_cell.get_c() == 0)
 		throw std::runtime_error("Invalid cell, contains a dimension that is zero");
 
-	point result_fsb;
+	point result_fsb{};
 	float result_d = std::numeric_limits<float>::max();
 	sym_op result_s;
 
@@ -469,39 +461,39 @@ std::tuple<float, point, sym_op> crystal::closest_symmetry_copy(point a, point b
 
 		auto fsb = t(fb);
 
-		while (fsb.m_x - 0.5f > fa.m_x)
+		while (fsb.x - 0.5f > fa.x)
 		{
-			fsb.m_x -= 1;
+			fsb.x -= 1;
 			s.m_ta -= 1;
 		}
 
-		while (fsb.m_x + 0.5f < fa.m_x)
+		while (fsb.x + 0.5f < fa.x)
 		{
-			fsb.m_x += 1;
+			fsb.x += 1;
 			s.m_ta += 1;
 		}
 
-		while (fsb.m_y - 0.5f > fa.m_y)
+		while (fsb.y - 0.5f > fa.y)
 		{
-			fsb.m_y -= 1;
+			fsb.y -= 1;
 			s.m_tb -= 1;
 		}
 
-		while (fsb.m_y + 0.5f < fa.m_y)
+		while (fsb.y + 0.5f < fa.y)
 		{
-			fsb.m_y += 1;
+			fsb.y += 1;
 			s.m_tb += 1;
 		}
 
-		while (fsb.m_z - 0.5f > fa.m_z)
+		while (fsb.z - 0.5f > fa.z)
 		{
-			fsb.m_z -= 1;
+			fsb.z -= 1;
 			s.m_tc -= 1;
 		}
 
-		while (fsb.m_z + 0.5f < fa.m_z)
+		while (fsb.z + 0.5f < fa.z)
 		{
-			fsb.m_z += 1;
+			fsb.z += 1;
 			s.m_tc += 1;
 		}
 

test/CMakeLists.txt

@@ -50,7 +50,6 @@ list(
 	reconstruction
 	validate-pdbx
 	cql
-	matrix
 )
 
 add_library(test-main OBJECT "${CMAKE_CURRENT_SOURCE_DIR}/test-main.cpp")

test/matrix-test.cpp

@@ -1,115 +0,0 @@
-/*-
- * SPDX-License-Identifier: BSD-2-Clause
- *
- * Copyright (c) 2025 NKI/AVL, Netherlands Cancer Institute
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions are met:
- *
- * 1. Redistributions of source code must retain the above copyright notice, this
- *    list of conditions and the following disclaimer
- * 2. Redistributions in binary form must reproduce the above copyright notice,
- *    this list of conditions and the following disclaimer in the documentation
- *    and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
- * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
- * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
- * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
- * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
- * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
- * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
- * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
-
-#include "cif++/matrix.hpp"
-#include "test-main.hpp"
-
-#include <catch2/catch_test_macros.hpp>
-#include <cif++/cif++.hpp>
-
-TEST_CASE("m1")
-{
-	cif::matrix3x3<int> m = cif::identity_matrix<int>(3);
-
-	CHECK(cif::determinant(m) == 1);
-}
-
-TEST_CASE("m2")
-{
-	cif::matrix4x4<int> m = cif::identity_matrix<int>(4);
-
-	cif::sub_matrix<cif::matrix4x4<int>> ms(m, 1, 1);
-	CHECK(ms == cif::identity_matrix<int>(3));
-}
-
-TEST_CASE("m3")
-{
-	cif::matrix4x4<int> m{
-		{ 1, 2, 3, 4,      //
-			5, 6, 7, 8,    //
-			9, 10, 11, 12, //
-			13, 14, 15, 16 }
-	};
-	cif::sub_matrix<cif::matrix4x4<int>> ms(m, 1, 1);
-
-	cif::matrix3x3<int> t{
-		{ 1, 3, 4, 9, 11, 12, 13, 15, 16 }
-	};
-
-	CHECK(ms == t);
-}
-
-TEST_CASE("m4")
-{
-	cif::matrix4x4<int> m{
-		{
-			-2,
-			3,
-			1,
-			0,
-			4,
-			1,
-			-3,
-			2,
-			0,
-			-1,
-			2,
-			5,
-			3,
-			2,
-			0,
-			-4,
-		}
-	};
-
-	std::cout << m << "\n\n";
-
-	// std::cout << cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 0)) << "\n\n";
-	// std::cout << cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 1)) << "\n\n";
-	// std::cout << cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 2)) << "\n\n";
-	// std::cout << cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 3)) << "\n\n";
-
-	// std::cout << cif::determinant(cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 0))) << "\n\n";
-	// std::cout << cif::determinant(cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 1))) << "\n\n";
-	// std::cout << cif::determinant(cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 2))) << "\n\n";
-	// std::cout << cif::determinant(cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 3))) << "\n\n";
-
-	CHECK(cif::determinant(m) == 332);
-}
-
-// --------------------------------------------------------------------
-
-TEST_CASE("m5")
-{
-    cif::matrix4x4<float> m = cif::identity_matrix<float>(4);
-    cif::matrix_fixed<float, 1, 4> v({ 0, 0.5f, 0, 1.0f });
-
-    auto mv = v * m;
-
-    CHECK(mv == v);
-
-}
-

test/model-test.cpp

@@ -443,7 +443,7 @@ TEST_CASE("test_alternates_1")
 	const std::filesystem::path example(gTestDir / ".." / "examples" / "1cbs.cif.gz");
 	cif::file file(example.string());
 
-	auto &db = file.front();
+	// auto &db = file.front();
 
 	cif::mm::structure s(file);
 

test/reconstruction-test.cpp

@@ -41,10 +41,6 @@ TEST_CASE("reconstruct")
 
 	for (std::filesystem::directory_iterator i(gTestDir / "reconstruct"); i != std::filesystem::directory_iterator{}; ++i)
 	{
-if (i->path().filename().string() != "9d3j_fixDMC.cif")
-		continue;
-
-
 		std::cout << i->path() << '\n';
 
 		if (i->path().extension() == ".pdb")

test/unit-3d-test.cpp

@@ -87,7 +87,7 @@ TEST_CASE("t1")
 	cif::center_points(p1);
 
 	for (auto &p : p2)
-		p.rotate(q);
+		p = cif::rotate(p, q);
 
 	cif::center_points(p2);
 
@@ -98,7 +98,7 @@ TEST_CASE("t1")
 	CHECK_THAT(std::fmod(360 + angle, 360), Catch::Matchers::WithinRel(std::fmod(360 - angle0, 360), 0.01));
 
 	for (auto &p : p1)
-		p.rotate(q2);
+		p = q2 * p;
 
 	auto rmsd = cif::RMSd(p1, p2);
 
@@ -115,7 +115,7 @@ TEST_CASE("t2")
 		{ 1, 2, 0 }
 	};
 
-	cif::point xp = cif::cross_product(p[1] - p[0], p[2] - p[0]);
+	auto xp = glm::cross(p[1] - p[0], p[2] - p[0]);
 
 	auto q = cif::construct_from_angle_axis(45, xp);
 
@@ -132,16 +132,16 @@ TEST_CASE("t3")
 		{ 1, 2, 0 }
 	};
 
-	cif::point xp = cif::cross_product(p[1] - p[0], p[2] - p[0]);
+	cif::point xp = glm::cross(p[1] - p[0], p[2] - p[0]);
 
 	auto q = cif::construct_from_angle_axis(45, xp);
 
 	auto v = p[1];
 	v -= p[0];
-	v.rotate(q);
+	v = q * v;
 	v += p[0];
 
-	std::cout << v << '\n';
+	std::println(std::cout, "{}", v);
 
 	double a = cif::angle(v, p[0], p[1]);
 
@@ -165,22 +165,22 @@ TEST_CASE("dh_q_0")
 
 	auto q = cif::construct_from_angle_axis(90, axis);
 
-	p.rotate(q);
+	p = q * p;
 
-	REQUIRE(std::abs(p.m_x - 1.f) < 0.01f);
+	REQUIRE(std::abs(p.x - 1.f) < 0.01f);
-	REQUIRE(std::abs(p.m_y - 0.f) < 0.01f);
+	REQUIRE(std::abs(p.y - 0.f) < 0.01f);
-	REQUIRE(std::abs(p.m_z - 1.f) < 0.01f);
+	REQUIRE(std::abs(p.z - 1.f) < 0.01f);
 
 	a = cif::dihedral_angle(t[0], t[1], t[2], p);
 	REQUIRE(std::abs(a - 90.f) < 0.01f);
 
 	q = cif::construct_from_angle_axis(-90, axis);
 
-	p.rotate(q);
+	p = q * p;
 
-	REQUIRE(std::abs(p.m_x - 1.f) < 0.01f);
+	REQUIRE(std::abs(p.x - 1.f) < 0.01f);
-	REQUIRE(std::abs(p.m_y - 1.f) < 0.01f);
+	REQUIRE(std::abs(p.y - 1.f) < 0.01f);
-	REQUIRE(std::abs(p.m_z - 0.f) < 0.01f);
+	REQUIRE(std::abs(p.z - 0.f) < 0.01f);
 
 	a = cif::dihedral_angle(t[0], t[1], t[2], p);
 	REQUIRE(std::abs(a - 0.f) < 0.01f);
@@ -222,7 +222,7 @@ TEST_CASE("dh_q_1")
 	{
 		auto q = cif::construct_for_dihedral_angle(pts[0], pts[1], pts[2], pts[3], angle, 1);
 
-		pts[3].rotate(q, pts[2]);
+		pts[3] = cif::rotate(pts[3], q, pts[2]);
 
 		auto dh = cif::dihedral_angle(pts[0], pts[1], pts[2], pts[3]);
 		CHECK_THAT(dh, Catch::Matchers::WithinRel(angle, 0.1f));
@@ -248,9 +248,9 @@ TEST_CASE("dh_q_1")
 
         cif::point p1{ 1, 1, 1 };
         cif::point p2 = p1;
-        p2.rotate(q);
+        p2 = q * p2;
 
-        cif::matrix3x3<float> rot_c({ static_cast<float>(d[0]),
+        cif::glm::mat3 rot_c({ static_cast<float>(d[0]),
             static_cast<float>(d[1]),
             static_cast<float>(d[2]),
             static_cast<float>(d[3]),
@@ -262,9 +262,9 @@ TEST_CASE("dh_q_1")
 
         cif::point p3 = rot_c * p1;
 
-        CHECK_THAT(p2.m_x, Catch::Matchers::WithinRel(p3.m_x, 0.01f));
+        CHECK_THAT(p2.x, Catch::Matchers::WithinRel(p3.x, 0.01f));
-        CHECK_THAT(p2.m_y, Catch::Matchers::WithinRel(p3.m_y, 0.01f));
+        CHECK_THAT(p2.y, Catch::Matchers::WithinRel(p3.y, 0.01f));
-        CHECK_THAT(p2.m_z, Catch::Matchers::WithinRel(p3.m_z, 0.01f));
+        CHECK_THAT(p2.z, Catch::Matchers::WithinRel(p3.z, 0.01f));
     }
 }
  */
@@ -275,7 +275,7 @@ TEST_CASE("dh_q_1")
 // 	{
 // 		auto d = cif::kSymopNrTable[i].symop().data();
 
-// 		cif::matrix3x3<float> rot;
+// 		cif::glm::mat3 rot;
 // 		float Qxx = rot(0, 0) = d[0];
 // 		float Qxy = rot(0, 1) = d[1];
 // 		float Qxz = rot(0, 2) = d[2];
@@ -318,13 +318,13 @@ TEST_CASE("dh_q_1")
 
 // 		cif::point p1{ 1, 1, 1 };
 // 		cif::point p2 = p1;
-// 		p2.rotate(q);
+// 		p2 = q * p2;
 
 // 		cif::point p3 = rot * p1;
 
-// 		REQUIRE(p2.m_x == p3.m_x);
+// 		REQUIRE(p2.x == p3.x);
-// 		REQUIRE(p2.m_y == p3.m_y);
+// 		REQUIRE(p2.y == p3.y);
-// 		REQUIRE(p2.m_z == p3.m_z);
+// 		REQUIRE(p2.z == p3.z);
 // 	}
 // }
 
@@ -338,15 +338,15 @@ TEST_CASE("symm_1")
 
 	cif::point f = fractional(p, c);
 
-	CHECK_THAT(f.m_x, Catch::Matchers::WithinRel(0.1f, 0.01f));
+	CHECK_THAT(f.x, Catch::Matchers::WithinRel(0.1f, 0.01f));
-	CHECK_THAT(f.m_y, Catch::Matchers::WithinRel(0.1f, 0.01f));
+	CHECK_THAT(f.y, Catch::Matchers::WithinRel(0.1f, 0.01f));
-	CHECK_THAT(f.m_z, Catch::Matchers::WithinRel(0.1f, 0.01f));
+	CHECK_THAT(f.z, Catch::Matchers::WithinRel(0.1f, 0.01f));
 
 	cif::point o = orthogonal(f, c);
 
-	CHECK_THAT(o.m_x, Catch::Matchers::WithinRel(1.f, 0.01f));
+	CHECK_THAT(o.x, Catch::Matchers::WithinRel(1.f, 0.01f));
-	CHECK_THAT(o.m_y, Catch::Matchers::WithinRel(1.f, 0.01f));
+	CHECK_THAT(o.y, Catch::Matchers::WithinRel(1.f, 0.01f));
-	CHECK_THAT(o.m_z, Catch::Matchers::WithinRel(1.f, 0.01f));
+	CHECK_THAT(o.z, Catch::Matchers::WithinRel(1.f, 0.01f));
 }
 
 TEST_CASE("symm_2")
@@ -385,44 +385,15 @@ TEST_CASE("symm_4")
 	CHECK_THAT(distance(a, sg(a, c, "1_554"_symop)), Catch::Matchers::WithinRel(static_cast<float>(c.get_c()), 0.01f));
 
 	auto sb2 = sg(b, c, "4_565"_symop);
-	CHECK_THAT(sb.m_x, Catch::Matchers::WithinRel(sb2.m_x, 0.01f));
+	CHECK_THAT(sb.x, Catch::Matchers::WithinRel(sb2.x, 0.01f));
-	CHECK_THAT(sb.m_y, Catch::Matchers::WithinRel(sb2.m_y, 0.01f));
+	CHECK_THAT(sb.y, Catch::Matchers::WithinRel(sb2.y, 0.01f));
-	CHECK_THAT(sb.m_z, Catch::Matchers::WithinRel(sb2.m_z, 0.01f));
+	CHECK_THAT(sb.z, Catch::Matchers::WithinRel(sb2.z, 0.01f));
 
 	CHECK_THAT(distance(a, sb2), Catch::Matchers::WithinRel(7.42f, 0.01f));
 }
 
 // --------------------------------------------------------------------
 
-TEST_CASE("symm_4wvp_1")
-{
-	using namespace cif::literals;
-
-	cif::file f(gTestDir / "4wvp.cif.gz");
-	f.front().set_validator(cif::validator_factory::instance().get("mmcif_pdbx.dic"));
-
-	auto &db = f.front();
-	cif::mm::structure s(db);
-
-	cif::crystal c(db);
-
-	cif::point p{ -78.722f, 98.528f, 11.994f };
-	auto a = s.get_residue("A", 10, "").get_atom_by_atom_id("O");
-
-	auto sp1 = c.symmetry_copy(a.get_location(), "2_565"_symop);
-	CHECK_THAT(sp1.m_x, Catch::Matchers::WithinAbs(p.m_x, 0.5f));
-	CHECK_THAT(sp1.m_y, Catch::Matchers::WithinAbs(p.m_y, 0.5f));
-	CHECK_THAT(sp1.m_z, Catch::Matchers::WithinAbs(p.m_z, 0.5f));
-
-	const auto &[d, sp2, so] = c.closest_symmetry_copy(p, a.get_location());
-
-	REQUIRE(d < 1);
-
-	CHECK_THAT(sp2.m_x, Catch::Matchers::WithinAbs(p.m_x, 0.5f));
-	CHECK_THAT(sp2.m_y, Catch::Matchers::WithinAbs(p.m_y, 0.5f));
-	CHECK_THAT(sp2.m_z, Catch::Matchers::WithinAbs(p.m_z, 0.5f));
-}
-
 TEST_CASE("symm_2bi3_1")
 {
 	cif::file f(gTestDir / "2bi3.cif.gz");
@@ -453,17 +424,17 @@ TEST_CASE("symm_2bi3_1")
 		auto sa1 = c.symmetry_copy(a1.get_location(), cif::sym_op(symm1));
 		auto sa2 = c.symmetry_copy(a2.get_location(), cif::sym_op(symm2));
 
-		CHECK_THAT(cif::distance(sa1, sa2), Catch::Matchers::WithinAbs(dist, 0.5f));
+		CHECK_THAT(glm::distance(sa1, sa2), Catch::Matchers::WithinAbs(dist, 0.01f));
 
 		auto pa1 = a1.get_location();
 
 		const auto &[d, p, so] = c.closest_symmetry_copy(pa1, a2.get_location());
 
-		CHECK_THAT(p.m_x, Catch::Matchers::WithinAbs(sa2.m_x, 0.5f));
+		CHECK_THAT(p.x, Catch::Matchers::WithinAbs(sa2.x, 0.01f));
-		CHECK_THAT(p.m_y, Catch::Matchers::WithinAbs(sa2.m_y, 0.5f));
+		CHECK_THAT(p.y, Catch::Matchers::WithinAbs(sa2.y, 0.01f));
-		CHECK_THAT(p.m_z, Catch::Matchers::WithinAbs(sa2.m_z, 0.5f));
+		CHECK_THAT(p.z, Catch::Matchers::WithinAbs(sa2.z, 0.01f));
 
-		CHECK_THAT(d, Catch::Matchers::WithinAbs(dist, 0.5f));
+		CHECK_THAT(d, Catch::Matchers::WithinAbs(dist, 0.01f));
 		REQUIRE(so.string() == symm2);
 	}
 }
@@ -491,25 +462,28 @@ TEST_CASE("symm_2bi3_1a")
 			 "ptnr2_label_asym_id", "ptnr2_label_seq_id", "ptnr2_auth_seq_id", "ptnr2_label_atom_id", "ptnr2_symmetry",
 			 "pdbx_dist_value"))
 	{
-		cif::point p1 = atom_site.find1<float, float, float>(
+		auto t1 = atom_site.find1<float, float, float>(
 			"label_asym_id"_key == asym1 and "label_seq_id"_key == seqid1 and "auth_seq_id"_key == authseqid1 and "label_atom_id"_key == atomid1,
 			"cartn_x", "cartn_y", "cartn_z");
-		cif::point p2 = atom_site.find1<float, float, float>(
+		auto t2 = atom_site.find1<float, float, float>(
 			"label_asym_id"_key == asym2 and "label_seq_id"_key == seqid2 and "auth_seq_id"_key == authseqid2 and "label_atom_id"_key == atomid2,
 			"cartn_x", "cartn_y", "cartn_z");
 
+		cif::point p1{ std::get<0>(t1), std::get<1>(t1), std::get<2>(t1) };
+		cif::point p2{ std::get<0>(t2), std::get<1>(t2), std::get<2>(t2) };
+
 		auto sa1 = c.symmetry_copy(p1, cif::sym_op(symm1));
 		auto sa2 = c.symmetry_copy(p2, cif::sym_op(symm2));
 
-		CHECK_THAT(cif::distance(sa1, sa2), Catch::Matchers::WithinAbs(dist, 0.5f));
+		CHECK_THAT(glm::distance(sa1, sa2), Catch::Matchers::WithinAbs(dist, 0.01f));
 
 		const auto &[d, p, so] = c.closest_symmetry_copy(p1, p2);
 
-		CHECK_THAT(p.m_x, Catch::Matchers::WithinAbs(sa2.m_x, 0.5f));
+		CHECK_THAT(p.x, Catch::Matchers::WithinAbs(sa2.x, 0.01f));
-		CHECK_THAT(p.m_y, Catch::Matchers::WithinAbs(sa2.m_y, 0.5f));
+		CHECK_THAT(p.y, Catch::Matchers::WithinAbs(sa2.y, 0.01f));
-		CHECK_THAT(p.m_z, Catch::Matchers::WithinAbs(sa2.m_z, 0.5f));
+		CHECK_THAT(p.z, Catch::Matchers::WithinAbs(sa2.z, 0.01f));
 
-		CHECK_THAT(d, Catch::Matchers::WithinAbs(dist, 0.5f));
+		CHECK_THAT(d, Catch::Matchers::WithinAbs(dist, 0.01f));
 		REQUIRE(so.string() == symm2);
 	}
 }
@@ -533,7 +507,7 @@ TEST_CASE("symm_3bwh_1")
 
 			const auto &[d, p, so] = c.closest_symmetry_copy(a1.get_location(), a2.get_location());
 
-			CHECK_THAT(d, Catch::Matchers::WithinAbs(distance(a1.get_location(), p), 0.5f));
+			CHECK_THAT(d, Catch::Matchers::WithinAbs(distance(a1.get_location(), p), 0.01f));
 		}
 	}
 }
@@ -582,7 +556,7 @@ TEST_CASE("smallest_sphere-1")
 	for (int i = 0; i < 1000; ++i)
 	{
 		auto [c, r] = cif::smallest_sphere_around_points(pts);
-		CHECK_THAT(cif::distance(c, cif::point{ 0, 0.743099928, 51.1741028 }), Catch::Matchers::WithinAbs(0.f, 0.01f));
+		CHECK_THAT(glm::distance(c, cif::point{ 0, 0.743099928, 51.1741028 }), Catch::Matchers::WithinAbs(0.f, 0.01f));
 		CHECK_THAT(r, Catch::Matchers::WithinAbs(7.31248331f, 0.01f));
 	}
 }