// Copyright (c) 2024-2026 The Regents of the University of Michigan.
// Part of hoomd-rs, released under the BSD 3-Clause License.

//! Implement `AppendMicrostate` for built-in site and boundary types.

use hoomd_geometry::shape::Hypercuboid;
use hoomd_gsd::hoomd::{AppendError, Dimensions, Frame, HoomdGsdFile};
use hoomd_manifold::{Hyperbolic, Spherical};
use hoomd_vector::{Angle, Cartesian, Versor};

use crate::{
    AppendMicrostate, Microstate,
    boundary::{Closed, Periodic},
    property::{OrientedHyperbolicPoint, OrientedPoint, Point},
};

impl<B, X> AppendMicrostate<B, Point<Cartesian<2>>, X, Closed<Hypercuboid<2>>> for HoomdGsdFile {
    #[inline]
    fn append_microstate(
        &mut self,
        microstate: &Microstate<B, Point<Cartesian<2>>, X, Closed<Hypercuboid<2>>>,
    ) -> Result<Frame<'_>, AppendError> {
        self.append_frame(microstate.step())?
            .configuration_box(microstate.boundary().0.to_gsd_box())?
            .configuration_dimensions(Dimensions::Two)?
            .particles_position(
                microstate
                    .iter_sites_tag_order()
                    .map(|s| s.properties.position)
                    .map(|p| [p[0], p[1], 0.0].into()),
            )
    }
}

impl<B, X> AppendMicrostate<B, Point<Cartesian<2>>, X, Periodic<Hypercuboid<2>>> for HoomdGsdFile {
    #[inline]
    fn append_microstate(
        &mut self,
        microstate: &Microstate<B, Point<Cartesian<2>>, X, Periodic<Hypercuboid<2>>>,
    ) -> Result<Frame<'_>, AppendError> {
        self.append_frame(microstate.step())?
            .configuration_box(microstate.boundary().shape().to_gsd_box())?
            .configuration_dimensions(Dimensions::Two)?
            .particles_position(
                microstate
                    .iter_sites_tag_order()
                    .map(|s| s.properties.position)
                    .map(|p| [p[0], p[1], 0.0].into()),
            )
    }
}

impl<B, X> AppendMicrostate<B, OrientedPoint<Cartesian<2>, Angle>, X, Closed<Hypercuboid<2>>>
    for HoomdGsdFile
{
    #[inline]
    fn append_microstate(
        &mut self,
        microstate: &Microstate<B, OrientedPoint<Cartesian<2>, Angle>, X, Closed<Hypercuboid<2>>>,
    ) -> Result<Frame<'_>, AppendError> {
        self.append_frame(microstate.step())?
            .configuration_box(microstate.boundary().0.to_gsd_box())?
            .configuration_dimensions(Dimensions::Two)?
            .particles_position(
                microstate
                    .iter_sites_tag_order()
                    .map(|s| s.properties.position)
                    .map(|p| [p[0], p[1], 0.0].into()),
            )?
            .particles_orientation(
                microstate
                    .iter_sites_tag_order()
                    .map(|s| s.properties.orientation.theta)
                    .map(|a| {
                        Versor::from_axis_angle(
                            [0.0, 0.0, 1.0]
                                .try_into()
                                .expect("hard-coded vector can be normalized"),
                            a,
                        )
                    }),
            )
    }
}

impl<B, X> AppendMicrostate<B, OrientedPoint<Cartesian<2>, Angle>, X, Periodic<Hypercuboid<2>>>
    for HoomdGsdFile
{
    #[inline]
    fn append_microstate(
        &mut self,
        microstate: &Microstate<B, OrientedPoint<Cartesian<2>, Angle>, X, Periodic<Hypercuboid<2>>>,
    ) -> Result<Frame<'_>, AppendError> {
        self.append_frame(microstate.step())?
            .configuration_box(microstate.boundary().shape().to_gsd_box())?
            .configuration_dimensions(Dimensions::Two)?
            .particles_position(
                microstate
                    .iter_sites_tag_order()
                    .map(|s| s.properties.position)
                    .map(|p| [p[0], p[1], 0.0].into()),
            )?
            .particles_orientation(
                microstate
                    .iter_sites_tag_order()
                    .map(|s| s.properties.orientation.theta)
                    .map(|a| {
                        Versor::from_axis_angle(
                            [0.0, 0.0, 1.0]
                                .try_into()
                                .expect("hard-coded vector can be normalized"),
                            a,
                        )
                    }),
            )
    }
}

impl<B, X> AppendMicrostate<B, Point<Cartesian<3>>, X, Closed<Hypercuboid<3>>> for HoomdGsdFile {
    #[inline]
    fn append_microstate(
        &mut self,
        microstate: &Microstate<B, Point<Cartesian<3>>, X, Closed<Hypercuboid<3>>>,
    ) -> Result<Frame<'_>, AppendError> {
        self.append_frame(microstate.step())?
            .configuration_box(microstate.boundary().0.to_gsd_box())?
            .configuration_dimensions(Dimensions::Three)?
            .particles_position(
                microstate
                    .iter_sites_tag_order()
                    .map(|s| s.properties.position),
            )
    }
}

impl<B, X> AppendMicrostate<B, Point<Cartesian<3>>, X, Periodic<Hypercuboid<3>>> for HoomdGsdFile {
    #[inline]
    fn append_microstate(
        &mut self,
        microstate: &Microstate<B, Point<Cartesian<3>>, X, Periodic<Hypercuboid<3>>>,
    ) -> Result<Frame<'_>, AppendError> {
        self.append_frame(microstate.step())?
            .configuration_box(microstate.boundary().shape().to_gsd_box())?
            .configuration_dimensions(Dimensions::Three)?
            .particles_position(
                microstate
                    .iter_sites_tag_order()
                    .map(|s| s.properties.position),
            )
    }
}

impl<B, X> AppendMicrostate<B, OrientedPoint<Cartesian<3>, Versor>, X, Closed<Hypercuboid<3>>>
    for HoomdGsdFile
{
    #[inline]
    fn append_microstate(
        &mut self,
        microstate: &Microstate<B, OrientedPoint<Cartesian<3>, Versor>, X, Closed<Hypercuboid<3>>>,
    ) -> Result<Frame<'_>, AppendError> {
        self.append_frame(microstate.step())?
            .configuration_box(microstate.boundary().0.to_gsd_box())?
            .configuration_dimensions(Dimensions::Three)?
            .particles_position(
                microstate
                    .iter_sites_tag_order()
                    .map(|s| s.properties.position),
            )?
            .particles_orientation(
                microstate
                    .iter_sites_tag_order()
                    .map(|s| s.properties.orientation),
            )
    }
}

impl<B, X> AppendMicrostate<B, OrientedPoint<Cartesian<3>, Versor>, X, Periodic<Hypercuboid<3>>>
    for HoomdGsdFile
{
    #[inline]
    fn append_microstate(
        &mut self,
        microstate: &Microstate<
            B,
            OrientedPoint<Cartesian<3>, Versor>,
            X,
            Periodic<Hypercuboid<3>>,
        >,
    ) -> Result<Frame<'_>, AppendError> {
        self.append_frame(microstate.step())?
            .configuration_box(microstate.boundary().shape().to_gsd_box())?
            .configuration_dimensions(Dimensions::Three)?
            .particles_position(
                microstate
                    .iter_sites_tag_order()
                    .map(|s| s.properties.position),
            )?
            .particles_orientation(
                microstate
                    .iter_sites_tag_order()
                    .map(|s| s.properties.orientation),
            )
    }
}

impl<B, X, C> AppendMicrostate<B, Point<Spherical<3>>, X, C> for HoomdGsdFile {
    #[inline]
    fn append_microstate(
        &mut self,
        microstate: &Microstate<B, Point<Spherical<3>>, X, C>,
    ) -> Result<Frame<'_>, AppendError> {
        self.append_frame(microstate.step())?
            .configuration_box([2.0, 2.0, 2.0, 0.0, 0.0, 0.0])?
            .configuration_dimensions(Dimensions::Three)?
            .particles_position(
                microstate
                    .iter_sites_tag_order()
                    .map(|s| *s.properties.position.point()),
            )
    }
}

impl<B, X, C> AppendMicrostate<B, Point<Spherical<4>>, X, C> for HoomdGsdFile {
    #[inline]
    fn append_microstate(
        &mut self,
        microstate: &Microstate<B, Point<Spherical<4>>, X, C>,
    ) -> Result<Frame<'_>, AppendError> {
        self.append_frame(microstate.step())?
            .configuration_box([2.0, 2.0, 2.0, 0.0, 0.0, 0.0])?
            .configuration_dimensions(Dimensions::Three)?
            .particles_position(microstate.iter_sites_tag_order().map(|s| {
                let proj: Vec<f64> = s.properties.position.stereographic_projection();
                [proj[0], proj[1], proj[2]].into()
            }))
    }
}

impl<B, X, C> AppendMicrostate<B, OrientedHyperbolicPoint<3, Angle>, X, C> for HoomdGsdFile {
    #[inline]
    fn append_microstate(
        &mut self,
        microstate: &Microstate<B, OrientedHyperbolicPoint<3, Angle>, X, C>,
    ) -> Result<Frame<'_>, AppendError> {
        self.append_frame(microstate.step())?
            .configuration_box([2.0, 2.0, 0.0, 0.0, 0.0, 0.0])?
            .configuration_dimensions(Dimensions::Two)?
            .particles_position(
                microstate
                    .iter_sites_tag_order()
                    .map(|s| s.properties.position.to_poincare())
                    .map(|p| [p[0], p[1], 0.0].into()),
            )?
            .particles_orientation(
                microstate
                    .iter_sites_tag_order()
                    .map(|s| s.properties.orientation.theta)
                    .map(|a| {
                        Versor::from_axis_angle(
                            [0.0, 0.0, 1.0]
                                .try_into()
                                .expect("hard-coded vector can be normalized"),
                            a,
                        )
                    }),
            )
    }
}

impl<B, X, C> AppendMicrostate<B, Point<Hyperbolic<3>>, X, C> for HoomdGsdFile {
    #[inline]
    fn append_microstate(
        &mut self,
        microstate: &Microstate<B, Point<Hyperbolic<3>>, X, C>,
    ) -> Result<Frame<'_>, AppendError> {
        self.append_frame(microstate.step())?
            .configuration_box([2.0, 2.0, 0.0, 0.0, 0.0, 0.0])?
            .configuration_dimensions(Dimensions::Two)?
            .particles_position(
                microstate
                    .iter_sites_tag_order()
                    .map(|s| s.properties.position.to_poincare())
                    .map(|p| [p[0], p[1], 0.0].into()),
            )
    }
}

#[cfg(test)]
mod test {
    use approxim::assert_relative_eq;
    use std::f64::consts::PI;
    use tempfile::tempdir;

    use super::*;
    use crate::{Body, boundary::Open};
    use hoomd_geometry::shape::{EightEight, Rectangle};
    use hoomd_gsd::file_layer::{GsdFile, Mode};

    #[test]
    fn point_closed_rectangle_2d() -> anyhow::Result<()> {
        let boundary = Closed(Rectangle {
            edge_lengths: [12.0.try_into()?, 18.0.try_into()?],
        });

        let microstate = Microstate::builder()
            .boundary(boundary)
            .bodies([
                Body::point(Cartesian::from([1.0, 0.0])),
                Body::point(Cartesian::from([-1.0, 2.0])),
            ])
            .step(1234)
            .try_build()?;

        let tmp_dir = tempdir()?;
        let path = tmp_dir.path().join("test.gsd");
        let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?;
        hoomd_gsd_file.append_microstate(&microstate)?;

        drop(hoomd_gsd_file);

        let gsd_file = GsdFile::open(path, Mode::Read)?;

        assert!(gsd_file.n_frames() == 1);

        let step = gsd_file.iter_scalars::<u64>(0, "configuration/step")?;
        itertools::assert_equal(step, [1234]);

        let positions = gsd_file.iter_arrays::<f32, 3>(0, "particles/position")?;
        itertools::assert_equal(positions, [[1.0, 0.0, 0.0], [-1.0, 2.0, 0.0]]);

        let dimensions = gsd_file.iter_scalars::<u8>(0, "configuration/dimensions")?;
        itertools::assert_equal(dimensions, [2]);

        let box_ = gsd_file.iter_scalars::<f32>(0, "configuration/box")?;
        itertools::assert_equal(box_, [12.0_f32, 18.0, 0.0, 0.0, 0.0, 0.0]);

        Ok(())
    }

    #[test]
    fn point_periodic_rectangle_2d() -> anyhow::Result<()> {
        let boundary = Periodic::new(
            0.0,
            Rectangle {
                edge_lengths: [12.0.try_into()?, 18.0.try_into()?],
            },
        )?;

        let microstate = Microstate::builder()
            .boundary(boundary)
            .bodies([
                Body::point(Cartesian::from([1.0, 0.0])),
                Body::point(Cartesian::from([-1.0, 2.0])),
            ])
            .step(1234)
            .try_build()?;

        let tmp_dir = tempdir()?;
        let path = tmp_dir.path().join("test.gsd");
        let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?;
        hoomd_gsd_file.append_microstate(&microstate)?;

        drop(hoomd_gsd_file);

        let gsd_file = GsdFile::open(path, Mode::Read)?;

        assert!(gsd_file.n_frames() == 1);

        let step = gsd_file.iter_scalars::<u64>(0, "configuration/step")?;
        itertools::assert_equal(step, [1234]);

        let positions = gsd_file.iter_arrays::<f32, 3>(0, "particles/position")?;
        itertools::assert_equal(positions, [[1.0, 0.0, 0.0], [-1.0, 2.0, 0.0]]);

        let dimensions = gsd_file.iter_scalars::<u8>(0, "configuration/dimensions")?;
        itertools::assert_equal(dimensions, [2]);

        let box_ = gsd_file.iter_scalars::<f32>(0, "configuration/box")?;
        itertools::assert_equal(box_, [12.0_f32, 18.0, 0.0, 0.0, 0.0, 0.0]);

        Ok(())
    }

    #[test]
    fn point_spherical_2d() -> anyhow::Result<()> {
        let microstate = Microstate::builder()
            .boundary(Open)
            .bodies([
                Body::point(Spherical::from_cartesian_coordinates(Cartesian::from([
                    -1.0, 0.0, 0.0,
                ]))),
                Body::point(Spherical::from_cartesian_coordinates(Cartesian::from([
                    0.0,
                    f64::sqrt(0.5),
                    f64::sqrt(0.5),
                ]))),
                Body::point(Spherical::from_cartesian_coordinates(Cartesian::from([
                    -f64::sqrt(0.25),
                    0.0,
                    f64::sqrt(0.75),
                ]))),
            ])
            .step(1234)
            .try_build()?;

        let tmp_dir = tempdir()?;
        let path = tmp_dir.path().join("test.gsd");
        let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?;
        hoomd_gsd_file.append_microstate(&microstate)?;

        drop(hoomd_gsd_file);

        let gsd_file = GsdFile::open(path, Mode::Read)?;

        assert!(gsd_file.n_frames() == 1);

        let step = gsd_file.iter_scalars::<u64>(0, "configuration/step")?;
        itertools::assert_equal(step, [1234]);

        let positions: Vec<[f32; 3]> = gsd_file
            .iter_arrays::<f32, 3>(0, "particles/position")?
            .collect();
        assert_relative_eq!(positions[0][0], -1.0);
        assert_relative_eq!(positions[0][1], 0.0);
        assert_relative_eq!(positions[0][2], 0.0);
        assert_relative_eq!(positions[1][0], 0.0);
        assert_relative_eq!(positions[1][1], f32::sqrt(0.5));
        assert_relative_eq!(positions[1][2], f32::sqrt(0.5));
        assert_relative_eq!(positions[2][0], -f32::sqrt(0.25));
        assert_relative_eq!(positions[2][1], 0.0);
        assert_relative_eq!(positions[2][2], f32::sqrt(0.75));
        Ok(())
    }

    #[test]
    fn point_spherical_3d() -> anyhow::Result<()> {
        let microstate = Microstate::builder()
            .boundary(Open)
            .bodies([
                Body::point(Spherical::from_cartesian_coordinates(Cartesian::from([
                    1.0, 0.0, 0.0, 0.0,
                ]))),
                Body::point(Spherical::from_cartesian_coordinates(Cartesian::from([
                    f64::sqrt(1.0 / 3.0),
                    -f64::sqrt(1.0 / 3.0),
                    f64::sqrt(1.0 / 3.0),
                    0.0,
                ]))),
                Body::point(Spherical::from_cartesian_coordinates(Cartesian::from([
                    0.0,
                    0.0,
                    f64::sqrt(0.5),
                    f64::sqrt(0.5),
                ]))),
            ])
            .step(1234)
            .try_build()?;

        let tmp_dir = tempdir()?;
        let path = tmp_dir.path().join("test.gsd");
        let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?;
        hoomd_gsd_file.append_microstate(&microstate)?;

        drop(hoomd_gsd_file);

        let gsd_file = GsdFile::open(path, Mode::Read)?;

        assert!(gsd_file.n_frames() == 1);

        let step = gsd_file.iter_scalars::<u64>(0, "configuration/step")?;
        itertools::assert_equal(step, [1234]);

        let positions: Vec<[f32; 3]> = gsd_file
            .iter_arrays::<f32, 3>(0, "particles/position")?
            .collect();
        assert_relative_eq!(positions[0][0], 1.0);
        assert_relative_eq!(positions[0][1], 0.0);
        assert_relative_eq!(positions[0][2], 0.0);
        assert_relative_eq!(positions[1][0], f32::sqrt(1.0 / 3.0));
        assert_relative_eq!(positions[1][1], -f32::sqrt(1.0 / 3.0));
        assert_relative_eq!(positions[1][2], f32::sqrt(1.0 / 3.0));
        assert_relative_eq!(positions[2][0], 0.0);
        assert_relative_eq!(positions[2][1], 0.0);
        assert_relative_eq!(positions[2][2], f32::sqrt(0.5) / (1.0 - f32::sqrt(0.5)));
        Ok(())
    }

    #[test]
    fn point_hyperbolic_2d_open() -> anyhow::Result<()> {
        let microstate = Microstate::builder()
            .boundary(Open)
            .bodies([
                Body::point(Hyperbolic::<3>::from_polar_coordinates(1.2, 0.0)),
                Body::point(Hyperbolic::<3>::from_polar_coordinates(0.6, 1.5)),
            ])
            .step(1234)
            .try_build()?;

        let tmp_dir = tempdir()?;
        let path = tmp_dir.path().join("test.gsd");
        let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?;
        hoomd_gsd_file.append_microstate(&microstate)?;

        drop(hoomd_gsd_file);

        let gsd_file = GsdFile::open(path, Mode::Read)?;

        assert!(gsd_file.n_frames() == 1);

        let step = gsd_file.iter_scalars::<u64>(0, "configuration/step")?;
        itertools::assert_equal(step, [1234]);

        let positions: Vec<[f32; 3]> = gsd_file
            .iter_arrays::<f32, 3>(0, "particles/position")?
            .collect();
        assert_relative_eq!(positions[0][0], (f32::sinh(1.2)) / (1.0 + f32::cosh(1.2)));
        assert_relative_eq!(positions[0][1], 0.0);
        assert_relative_eq!(
            positions[1][0],
            (f32::sinh(0.6) * f32::cos(1.5)) / (1.0 + f32::cosh(0.6))
        );
        assert_relative_eq!(
            positions[1][1],
            (f32::sinh(0.6) * f32::sin(1.5)) / (1.0 + f32::cosh(0.6))
        );
        Ok(())
    }

    #[test]
    fn point_hyperbolic_2d_eighteight_periodic() -> anyhow::Result<()> {
        let boundary = Periodic::new(0.6, EightEight {})?;
        let microstate = Microstate::builder()
            .boundary(boundary)
            .bodies([
                Body::point(Hyperbolic::<3>::from_polar_coordinates(0.2, 0.3)),
                Body::point(Hyperbolic::<3>::from_polar_coordinates(0.6, 0.7)),
            ])
            .step(1234)
            .try_build()?;

        let tmp_dir = tempdir()?;
        let path = tmp_dir.path().join("test.gsd");
        let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?;
        hoomd_gsd_file.append_microstate(&microstate)?;

        drop(hoomd_gsd_file);

        let gsd_file = GsdFile::open(path, Mode::Read)?;

        assert!(gsd_file.n_frames() == 1);

        let step = gsd_file.iter_scalars::<u64>(0, "configuration/step")?;
        itertools::assert_equal(step, [1234]);

        let positions: Vec<[f32; 3]> = gsd_file
            .iter_arrays::<f32, 3>(0, "particles/position")?
            .collect();
        assert_relative_eq!(
            positions[0][0],
            (f32::sinh(0.2) * f32::cos(0.3)) / (1.0 + f32::cosh(0.2))
        );
        assert_relative_eq!(
            positions[0][1],
            (f32::sinh(0.2) * f32::sin(0.3)) / (1.0 + f32::cosh(0.2))
        );
        assert_relative_eq!(
            positions[1][0],
            (f32::sinh(0.6) * f32::cos(0.7)) / (1.0 + f32::cosh(0.6))
        );
        assert_relative_eq!(
            positions[1][1],
            (f32::sinh(0.6) * f32::sin(0.7)) / (1.0 + f32::cosh(0.6))
        );
        Ok(())
    }

    #[test]
    fn oriented_point_hyperbolic_2d_open() -> anyhow::Result<()> {
        let microstate = Microstate::builder()
            .boundary(Open)
            .bodies([
                Body {
                    properties: OrientedHyperbolicPoint {
                        position: Hyperbolic::<3>::from_polar_coordinates(0.5, 0.6),
                        orientation: Angle::from(0.3),
                    },
                    sites: vec![OrientedHyperbolicPoint {
                        position: Hyperbolic::<3>::default(),
                        orientation: Angle::default(),
                    }],
                },
                Body {
                    properties: OrientedHyperbolicPoint {
                        position: Hyperbolic::<3>::from_polar_coordinates(0.9, 0.3),
                        orientation: Angle::from(1.2),
                    },
                    sites: vec![OrientedHyperbolicPoint {
                        position: Hyperbolic::<3>::default(),
                        orientation: Angle::default(),
                    }],
                },
            ])
            .step(1234)
            .try_build()?;

        let tmp_dir = tempdir()?;
        let path = tmp_dir.path().join("test.gsd");
        let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?;
        hoomd_gsd_file.append_microstate(&microstate)?;

        drop(hoomd_gsd_file);

        let gsd_file = GsdFile::open(path, Mode::Read)?;

        assert!(gsd_file.n_frames() == 1);

        let step = gsd_file.iter_scalars::<u64>(0, "configuration/step")?;
        itertools::assert_equal(step, [1234]);

        let positions: Vec<[f32; 3]> = gsd_file
            .iter_arrays::<f32, 3>(0, "particles/position")?
            .collect();
        assert_relative_eq!(
            positions[0][0],
            (f32::sinh(0.5) * f32::cos(0.6)) / (1.0 + f32::cosh(0.5))
        );
        assert_relative_eq!(
            positions[0][1],
            (f32::sinh(0.5) * f32::sin(0.6)) / (1.0 + f32::cosh(0.5))
        );
        assert_relative_eq!(
            positions[1][0],
            (f32::sinh(0.9) * f32::cos(0.3)) / (1.0 + f32::cosh(0.9))
        );
        assert_relative_eq!(
            positions[1][1],
            (f32::sinh(0.9) * f32::sin(0.3)) / (1.0 + f32::cosh(0.9))
        );
        Ok(())
    }

    #[test]
    fn oriented_point_hyperbolic_2d_periodic_eighteight() -> anyhow::Result<()> {
        let boundary = Periodic::new(0.6, EightEight {})?;
        let microstate = Microstate::builder()
            .boundary(boundary)
            .bodies([
                Body {
                    properties: OrientedHyperbolicPoint {
                        position: Hyperbolic::<3>::from_polar_coordinates(0.5, 0.6),
                        orientation: Angle::from(0.3),
                    },
                    sites: vec![OrientedHyperbolicPoint {
                        position: Hyperbolic::<3>::default(),
                        orientation: Angle::default(),
                    }],
                },
                Body {
                    properties: OrientedHyperbolicPoint {
                        position: Hyperbolic::<3>::from_polar_coordinates(0.9, 0.3),
                        orientation: Angle::from(1.2),
                    },
                    sites: vec![OrientedHyperbolicPoint {
                        position: Hyperbolic::<3>::default(),
                        orientation: Angle::default(),
                    }],
                },
            ])
            .step(1234)
            .try_build()?;

        let tmp_dir = tempdir()?;
        let path = tmp_dir.path().join("test.gsd");
        let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?;
        hoomd_gsd_file.append_microstate(&microstate)?;

        drop(hoomd_gsd_file);

        let gsd_file = GsdFile::open(path, Mode::Read)?;

        assert!(gsd_file.n_frames() == 1);

        let step = gsd_file.iter_scalars::<u64>(0, "configuration/step")?;
        itertools::assert_equal(step, [1234]);

        let positions: Vec<[f32; 3]> = gsd_file
            .iter_arrays::<f32, 3>(0, "particles/position")?
            .collect();
        assert_relative_eq!(
            positions[0][0],
            (f32::sinh(0.5) * f32::cos(0.6)) / (1.0 + f32::cosh(0.5))
        );
        assert_relative_eq!(
            positions[0][1],
            (f32::sinh(0.5) * f32::sin(0.6)) / (1.0 + f32::cosh(0.5))
        );
        assert_relative_eq!(
            positions[1][0],
            (f32::sinh(0.9) * f32::cos(0.3)) / (1.0 + f32::cosh(0.9))
        );
        assert_relative_eq!(
            positions[1][1],
            (f32::sinh(0.9) * f32::sin(0.3)) / (1.0 + f32::cosh(0.9))
        );
        Ok(())
    }

    #[test]
    fn oriented_point_closed_rectangle_2d() -> anyhow::Result<()> {
        let boundary = Closed(Rectangle {
            edge_lengths: [12.0.try_into()?, 18.0.try_into()?],
        });

        let site = OrientedPoint {
            position: Cartesian::from([0.0, 0.0]),
            orientation: Angle::default(),
        };
        let a = OrientedPoint {
            position: Cartesian::from([1.0, 0.0]),
            orientation: Angle::from(PI / 2.0),
        };
        let b = OrientedPoint {
            position: Cartesian::from([-1.0, 2.0]),
            orientation: Angle::from(PI),
        };
        let body_a = Body {
            properties: a,
            sites: [site].into(),
        };
        let body_b = Body {
            properties: b,
            sites: [site].into(),
        };

        let microstate = Microstate::builder()
            .boundary(boundary)
            .bodies([body_a, body_b])
            .step(1234)
            .try_build()?;

        let tmp_dir = tempdir()?;
        let path = tmp_dir.path().join("test.gsd");
        let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?;
        hoomd_gsd_file.append_microstate(&microstate)?;

        drop(hoomd_gsd_file);

        let gsd_file = GsdFile::open(path, Mode::Read)?;

        assert!(gsd_file.n_frames() == 1);

        let step = gsd_file.iter_scalars::<u64>(0, "configuration/step")?;
        itertools::assert_equal(step, [1234]);

        let positions = gsd_file.iter_arrays::<f32, 3>(0, "particles/position")?;
        itertools::assert_equal(positions, [[1.0, 0.0, 0.0], [-1.0, 2.0, 0.0]]);

        assert!(
            gsd_file
                .iter_arrays::<f32, 4>(0, "particles/orientation")?
                .count()
                == 2
        );

        let dimensions = gsd_file.iter_scalars::<u8>(0, "configuration/dimensions")?;
        itertools::assert_equal(dimensions, [2]);

        let box_ = gsd_file.iter_scalars::<f32>(0, "configuration/box")?;
        itertools::assert_equal(box_, [12.0_f32, 18.0, 0.0, 0.0, 0.0, 0.0]);

        Ok(())
    }

    #[test]
    fn oriented_point_periodic_rectangle_2d() -> anyhow::Result<()> {
        let boundary = Periodic::new(
            0.0,
            Rectangle {
                edge_lengths: [12.0.try_into()?, 18.0.try_into()?],
            },
        )?;

        let site = OrientedPoint {
            position: Cartesian::from([0.0, 0.0]),
            orientation: Angle::default(),
        };
        let a = OrientedPoint {
            position: Cartesian::from([1.0, 0.0]),
            orientation: Angle::from(PI / 2.0),
        };
        let b = OrientedPoint {
            position: Cartesian::from([-1.0, 2.0]),
            orientation: Angle::from(PI),
        };
        let body_a = Body {
            properties: a,
            sites: [site].into(),
        };
        let body_b = Body {
            properties: b,
            sites: [site].into(),
        };

        let microstate = Microstate::builder()
            .boundary(boundary)
            .bodies([body_a, body_b])
            .step(1234)
            .try_build()?;

        let tmp_dir = tempdir()?;
        let path = tmp_dir.path().join("test.gsd");
        let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?;
        hoomd_gsd_file.append_microstate(&microstate)?;

        drop(hoomd_gsd_file);

        let gsd_file = GsdFile::open(path, Mode::Read)?;

        assert!(gsd_file.n_frames() == 1);

        let step = gsd_file.iter_scalars::<u64>(0, "configuration/step")?;
        itertools::assert_equal(step, [1234]);

        let positions = gsd_file.iter_arrays::<f32, 3>(0, "particles/position")?;
        itertools::assert_equal(positions, [[1.0, 0.0, 0.0], [-1.0, 2.0, 0.0]]);

        assert!(
            gsd_file
                .iter_arrays::<f32, 4>(0, "particles/orientation")?
                .count()
                == 2
        );

        let dimensions = gsd_file.iter_scalars::<u8>(0, "configuration/dimensions")?;
        itertools::assert_equal(dimensions, [2]);

        let box_ = gsd_file.iter_scalars::<f32>(0, "configuration/box")?;
        itertools::assert_equal(box_, [12.0_f32, 18.0, 0.0, 0.0, 0.0, 0.0]);

        Ok(())
    }

    #[test]
    fn point_closed_cuboid_3d() -> anyhow::Result<()> {
        let boundary = Closed(Hypercuboid {
            edge_lengths: [12.0.try_into()?, 18.0.try_into()?, 24.0.try_into()?],
        });

        let microstate = Microstate::builder()
            .boundary(boundary)
            .bodies([
                Body::point(Cartesian::from([1.0, 0.0, 4.0])),
                Body::point(Cartesian::from([-1.0, 2.0, -2.0])),
            ])
            .step(1234)
            .try_build()?;

        let tmp_dir = tempdir()?;
        let path = tmp_dir.path().join("test.gsd");
        let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?;
        hoomd_gsd_file.append_microstate(&microstate)?;

        drop(hoomd_gsd_file);

        let gsd_file = GsdFile::open(path, Mode::Read)?;

        assert!(gsd_file.n_frames() == 1);

        let step = gsd_file.iter_scalars::<u64>(0, "configuration/step")?;
        itertools::assert_equal(step, [1234]);

        let positions = gsd_file.iter_arrays::<f32, 3>(0, "particles/position")?;
        itertools::assert_equal(positions, [[1.0, 0.0, 4.0], [-1.0, 2.0, -2.0]]);

        let dimensions = gsd_file.iter_scalars::<u8>(0, "configuration/dimensions")?;
        itertools::assert_equal(dimensions, [3]);

        let box_ = gsd_file.iter_scalars::<f32>(0, "configuration/box")?;
        itertools::assert_equal(box_, [12.0_f32, 18.0, 24.0, 0.0, 0.0, 0.0]);

        Ok(())
    }

    #[test]
    fn point_periodic_cuboid_3d() -> anyhow::Result<()> {
        let boundary = Periodic::new(
            0.0,
            Hypercuboid {
                edge_lengths: [12.0.try_into()?, 18.0.try_into()?, 24.0.try_into()?],
            },
        )?;

        let microstate = Microstate::builder()
            .boundary(boundary)
            .bodies([
                Body::point(Cartesian::from([1.0, 0.0, 4.0])),
                Body::point(Cartesian::from([-1.0, 2.0, -2.0])),
            ])
            .step(1234)
            .try_build()?;

        let tmp_dir = tempdir()?;
        let path = tmp_dir.path().join("test.gsd");
        let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?;
        hoomd_gsd_file.append_microstate(&microstate)?;

        drop(hoomd_gsd_file);

        let gsd_file = GsdFile::open(path, Mode::Read)?;

        assert!(gsd_file.n_frames() == 1);

        let step = gsd_file.iter_scalars::<u64>(0, "configuration/step")?;
        itertools::assert_equal(step, [1234]);

        let positions = gsd_file.iter_arrays::<f32, 3>(0, "particles/position")?;
        itertools::assert_equal(positions, [[1.0, 0.0, 4.0], [-1.0, 2.0, -2.0]]);

        let dimensions = gsd_file.iter_scalars::<u8>(0, "configuration/dimensions")?;
        itertools::assert_equal(dimensions, [3]);

        let box_ = gsd_file.iter_scalars::<f32>(0, "configuration/box")?;
        itertools::assert_equal(box_, [12.0_f32, 18.0, 24.0, 0.0, 0.0, 0.0]);

        Ok(())
    }

    #[test]
    fn oriented_point_closed_cuboid_3d() -> anyhow::Result<()> {
        let boundary = Closed(Hypercuboid {
            edge_lengths: [12.0.try_into()?, 18.0.try_into()?, 24.0.try_into()?],
        });

        let site = OrientedPoint {
            position: Cartesian::default(),
            orientation: Versor::default(),
        };
        let a = OrientedPoint {
            position: Cartesian::from([1.0, 0.0, 4.0]),
            orientation: Versor::from_axis_angle([1.0, 0.0, 0.0].try_into()?, PI / 2.0),
        };
        let b = OrientedPoint {
            position: Cartesian::from([-1.0, 2.0, -2.0]),
            orientation: Versor::from_axis_angle([0.0, 1.0, 0.0].try_into()?, PI / 2.0),
        };
        let body_a = Body {
            properties: a,
            sites: [site].into(),
        };
        let body_b = Body {
            properties: b,
            sites: [site].into(),
        };

        let microstate = Microstate::builder()
            .boundary(boundary)
            .bodies([body_a, body_b])
            .step(1234)
            .try_build()?;

        let tmp_dir = tempdir()?;
        let path = tmp_dir.path().join("test.gsd");
        let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?;
        hoomd_gsd_file.append_microstate(&microstate)?;

        drop(hoomd_gsd_file);

        let gsd_file = GsdFile::open(path, Mode::Read)?;

        assert!(gsd_file.n_frames() == 1);

        let step = gsd_file.iter_scalars::<u64>(0, "configuration/step")?;
        itertools::assert_equal(step, [1234]);

        let positions = gsd_file.iter_arrays::<f32, 3>(0, "particles/position")?;
        itertools::assert_equal(positions, [[1.0, 0.0, 4.0], [-1.0, 2.0, -2.0]]);

        assert!(
            gsd_file
                .iter_arrays::<f32, 4>(0, "particles/orientation")?
                .count()
                == 2
        );

        let dimensions = gsd_file.iter_scalars::<u8>(0, "configuration/dimensions")?;
        itertools::assert_equal(dimensions, [3]);

        let box_ = gsd_file.iter_scalars::<f32>(0, "configuration/box")?;
        itertools::assert_equal(box_, [12.0_f32, 18.0, 24.0, 0.0, 0.0, 0.0]);

        Ok(())
    }

    #[test]
    fn oriented_point_periodic_cuboid_3d() -> anyhow::Result<()> {
        let boundary = Periodic::new(
            0.0,
            Hypercuboid {
                edge_lengths: [12.0.try_into()?, 18.0.try_into()?, 24.0.try_into()?],
            },
        )?;

        let site = OrientedPoint {
            position: Cartesian::from([0.0, 0.0, 0.0]),
            orientation: Versor::default(),
        };
        let a = OrientedPoint {
            position: Cartesian::from([1.0, 0.0, 4.0]),
            orientation: Versor::from_axis_angle([1.0, 0.0, 0.0].try_into()?, PI / 2.0),
        };
        let b = OrientedPoint {
            position: Cartesian::from([-1.0, 2.0, -2.0]),
            orientation: Versor::from_axis_angle([0.0, 1.0, 0.0].try_into()?, PI / 2.0),
        };
        let body_a = Body {
            properties: a,
            sites: [site].into(),
        };
        let body_b = Body {
            properties: b,
            sites: [site].into(),
        };

        let microstate = Microstate::builder()
            .boundary(boundary)
            .bodies([body_a, body_b])
            .step(1234)
            .try_build()?;

        let tmp_dir = tempdir()?;
        let path = tmp_dir.path().join("test.gsd");
        let mut hoomd_gsd_file = HoomdGsdFile::create(path.clone())?;
        hoomd_gsd_file.append_microstate(&microstate)?;

        drop(hoomd_gsd_file);

        let gsd_file = GsdFile::open(path, Mode::Read)?;

        assert!(gsd_file.n_frames() == 1);

        let step = gsd_file.iter_scalars::<u64>(0, "configuration/step")?;
        itertools::assert_equal(step, [1234]);

        let positions = gsd_file.iter_arrays::<f32, 3>(0, "particles/position")?;
        itertools::assert_equal(positions, [[1.0, 0.0, 4.0], [-1.0, 2.0, -2.0]]);

        assert!(
            gsd_file
                .iter_arrays::<f32, 4>(0, "particles/orientation")?
                .count()
                == 2
        );

        let dimensions = gsd_file.iter_scalars::<u8>(0, "configuration/dimensions")?;
        itertools::assert_equal(dimensions, [3]);

        let box_ = gsd_file.iter_scalars::<f32>(0, "configuration/box")?;
        itertools::assert_equal(box_, [12.0_f32, 18.0, 24.0, 0.0, 0.0, 0.0]);

        Ok(())
    }
}
