//! # Nbody barnes hut
//!
//! `nbody_barnes_hut` is designed to facilitate the simulation of N-body systems in `O(nlogn)` time.
//! This is useful for many applications: common ones are gravitational simulations and electrostatic simulations.
//! Simulations in 2D and 3D are both supported.
//!
//! This crate is not multithreaded. Rather, call `calc_forces_on_particle` in a multithreaded loop (for example, with `rayon`).
//! The time to create the tree is negligible in comparison to the time used calculating forces.
//!
//! # Example
//!
//! Here is a basic 3D gravitational simulator:
//! ```
//!
//! use rand::Rng;
//! use nbody_barnes_hut::particle_3d::Particle3D;
//! use nbody_barnes_hut::vector_3d::Vector3D;
//! use nbody_barnes_hut::barnes_hut_3d::OctTree;
//!
//! const G: f64 = 6.67E-11; // Newton's Gravitational Constant
//!
//! // Create 10 000 random points
//! let mut rng = rand::thread_rng();
//! let points: Vec<Particle3D> = (0..10_000)
//!     .map(|_| {
//!         let pos = Vector3D::new(
//!             rng.gen_range(-1000.0, 1000.0),
//!             rng.gen_range(-1000.0, 1000.0),
//!             rng.gen_range(-1000.0, 1000.0),
//!         );
//!         Particle3D::new(pos, 30.0)
//!     })
//!     .collect();
//!
//! // This is pretty hacky
//! let points_ref = &points.iter().collect::<Vec<&Particle3D>>()[..];
//!
//! let tree = OctTree::new(points_ref, 0.5);
//!
//! for p in &points {
//!     // Do something with this value
//!     let acceleration_on_particle = tree.calc_forces_on_particle(
//!         p.position,
//!         (),
//!         |d_squared, mass, dis_vec, _| {
//!             // dis_vec is not normalized, so we have to normalize it here
//!             G * mass * dis_vec / (d_squared * d_squared.sqrt())
//!         },
//!     );
//! }
//! ```
pub mod barnes_hut_2d;
pub mod barnes_hut_3d;
pub mod particle_2d;
pub mod particle_3d;
mod vector;
pub mod vector_2d;
pub mod vector_3d;

#[cfg(test)]
mod tests {
	use crate::barnes_hut_2d::QuadTree;
	use crate::barnes_hut_3d::OctTree;
	use crate::particle_2d::Particle2D;
	use crate::particle_3d::Particle3D;
	use crate::vector::Dist;
	use crate::vector_2d::Vector2D;
	use crate::vector_3d::Vector3D;
	use std::ops::Sub;

	#[test]
	fn gravity_tests_3d() {
		let particle_1 = Particle3D::new(Vector3D::ZERO, 10.0);
		let p2_position = Vector3D::new(5.0, 10.0, 15.0);
		let particle_2 = Particle3D::new(p2_position, 100.0);

		let dist_squared = 5.0 * 5.0 + 10.0 * 10.0 + 15.0 * 15.0;

		let force_expected = -10.0 * p2_position.normalize() / dist_squared;

		let tree = OctTree::new(&[&particle_1, &particle_2][..], 0.5);
		let force_actual = tree.calc_forces_on_particle(
			particle_2.position,
			(),
			|dist_squared, mass, dis_vec, _| mass * dis_vec / (dist_squared * dist_squared.sqrt()),
		);

		assert_eq!(force_expected, force_actual);

		let tree = OctTree::new(&[&particle_1][..], 0.5);
		assert_vector_eq(
			Vector3D::ZERO,
			tree.calc_forces_on_particle(
				particle_1.position,
				(),
				|dist_squared, mass, dis_vec, _| {
					mass * dis_vec / (dist_squared * dist_squared.sqrt())
				},
			),
		)
	}

	#[test]
	fn electrostatics_3d() {
		let charge_one = 100.0;
		// NOTE we are using the charge as mass
		let particle_1 = Particle3D::new(Vector3D::ZERO, charge_one);

		let charge_two = -50.0;
		let particle_2_pos = Vector3D::new(15.0, 20.0, -1.1);
		let particle_2 = Particle3D::new(particle_2_pos, charge_two);

		let dis_squared = 15.0 * 15.0 + 20.0 * 20.0 + 1.1 * 1.1;

		// F = Kq1q2/r^2
		let accel_expected = -charge_one * charge_two * particle_2_pos.normalize() / dis_squared;

		let tree = OctTree::new(&[&particle_1, &particle_2][..], 0.5);
		let accel_actual = tree.calc_forces_on_particle(
			particle_2_pos,
			charge_two,
			|dist_squared, charge_other, dis_vec, charge_self| {
				charge_other * charge_self * dis_vec / (dist_squared * dist_squared.sqrt())
			},
		);

		assert_vector_eq(accel_expected, accel_actual);
	}

	#[test]
	fn gravity_tests_2d() {
		let particle_1 = Particle2D::new(Vector2D::ZERO, 10.0);
		let p2_position = Vector2D::new(5.0, 10.0);
		let particle_2 = Particle2D::new(p2_position, 100.0);

		let dist_squared = 5.0 * 5.0 + 10.0 * 10.0;

		let force_expected = -10.0 * p2_position.normalize() / dist_squared;

		let tree = QuadTree::new(&[&particle_1, &particle_2][..], 0.5);
		let force_actual = tree.calc_forces_on_particle(
			particle_2.position,
			(),
			|dist_squared, mass, dis_vec, _| mass * dis_vec / (dist_squared * dist_squared.sqrt()),
		);

		assert_vector_eq(force_expected, force_actual);

		let tree = QuadTree::new(&[&particle_1][..], 0.5);
		assert_vector_eq(
			Vector2D::ZERO,
			tree.calc_forces_on_particle(
				particle_1.position,
				(),
				|dist_squared, mass, dis_vec, _| {
					mass * dis_vec / (dist_squared * dist_squared.sqrt())
				},
			),
		)
	}

	#[test]
	fn electrostatics_2d() {
		let charge_one = 100.0;
		// NOTE we are using the charge as mass
		let particle_1 = Particle2D::new(Vector2D::ZERO, charge_one);

		let charge_two = -50.0;
		let particle_2_pos = Vector2D::new(15.0, -20.0);
		let particle_2 = Particle2D::new(particle_2_pos, charge_two);

		let dis_squared = 15.0 * 15.0 + 20.0 * 20.0;

		// F = Kq1q2/r^2
		let accel_expected = -charge_one * charge_two * particle_2_pos.normalize() / dis_squared;

		let tree = QuadTree::new(&[&particle_1, &particle_2][..], 0.5);
		let accel_actual = tree.calc_forces_on_particle(
			particle_2_pos,
			charge_two,
			|dist_squared, charge_other, dis_vec, charge_self| {
				charge_other * charge_self * dis_vec / (dist_squared * dist_squared.sqrt())
			},
		);

		assert_vector_eq(accel_expected, accel_actual);
	}

	fn assert_vector_eq<T: Sub<Output = T> + Dist>(a: T, b: T) {
		let d = (a - b).dist();
		assert!(d < 0.00001);
	}
}