Particular is a crate providing a simple way to compute N-body interaction of particles in Rust.

Please note that this branch is for development purposes and may not represent the latest stable release of the library. For the most recent stable version, refer to the latest branch.

The main goal of this crate is to provide users with a simple API to set up N-body simulations that can easily be integrated into existing game and physics engines. Thus it does not include anything related to numerical integration or other similar tools and instead only focuses on the calculations of the interactions between particles.

In order to do this, particular provides multiple algorithms to compute the interactions between particles, allowing users to choose the one that best fits their needs in terms of performance and accuracy.

There are currently 2 algorithms: Brute-force and Barnes-Hut.

Generally speaking, the Brute-force algorithm is more accurate, but slower. The Barnes-Hut algorithm allows trading accuracy for speed by increasing the theta parameter.
You can see more about their relative performance here.

Particular uses rayon for parallelization and wgpu for GPU computation.
Enable the respective parallel and gpu features to access the relevant algorithms.

At its core, Particular is a simple set of traits designed to simplify the implementation of different algorithms to compute interactions between particles. These interactions can be of any nature, but a common example is the gravitational interaction between particles, which is provided by the gravity module and the main focus of the library.
This module provides implementations using popular vector math libraries such as glam and nalgebra. You you will need to enable the corresponding features to use them.

The Position and Mass traits provide methods to retrieve the position, mass and gravitational parameter of a particle. Implementing these traits for a type will result in multiple blanket implementations that will allow it to be used with all available algorithms for computing different interactions.

When the type has fields named position and mass or mu, you can derive these traits. An optional attribute #[G = ...] defining the gravitational constant can be added on the field. If the attribute is missing, the value of the gravitational constant is 1.0.

use particular::prelude::*;
use glam::DVec3;

#[derive(Position, Mass)]
struct Planet {
    velocity: DVec3,
    position: DVec3,
    #[G = 6.67430e-11]
    mass: f64,
}

If you can't implement these traits, you can use the fact that it is implemented for tuples of a position and a mass or gravitational parameter instead of creating an intermediate type.

use particular::prelude::*;
use glam::DVec3;

let particle = (DVec3::new(-1.0, 0.0, 1.0), 5.0);

assert_eq!(particle.position(), DVec3::new(-1.0, 0.0, 1.0));
assert_eq!(particle.mass(), 5.0);
assert_eq!(particle.mu(), 5.0);

In order to compute the gravitational interaction between particles, you can use the various methods implemented on storages of particles.
Note that most algorithms return iterators that borrow the particles, and as such, rust's borrowing rules will prevent you from mutating these particles as you iterate the computed interactions, even if the interaction does not use the fields you are mutating. You can either immediately collect the interactions or use separate collections for the properties you need to mutate.

use particular::prelude::*;
use particular::gravity::newtonian::Acceleration;

let accelerations: Vec<_> = planets.brute_force(Acceleration::checked()).collect::<Vec<_>>();

for (planet, acceleration) in planets.iter_mut().zip(accelerations) {
    planet.velocity += acceleration * DT;
    planet.position += planet.velocity * DT;
}

You can alternatively use indices, but note that this requires making sure the positions don't change as the interaction is computed or you will get unexpected behaviour.

use particular::prelude::*;
use particular::gravity::newtonian::AccelerationSoftened;

for i in 0..planets.len() {
    let between = Between(&planets[i], planets.as_slice());
    let acceleration = between.brute_force(AccelerationSoftened::checked(1e-3));
    planets[i].velocity += acceleration * DT;
}

for planet in planets.iter_mut() {
    planet.position += planet.velocity * DT;
}

use particular::prelude::*;
use particular::gravity::newtonian::{Acceleration, TreeData};

use glam::Vec3;

// Particles here are simple gravitational fields.
let particles: Vec<GravitationalField<Vec3, f32>> = vec![
    // ...
];

// Function to determine if a field will affect particles.
fn is_massive(p: &GravitationalField<Vec3, f32>) -> bool {
    p.m != 0.0
}

// Create an `Ordered` storage to split massive and massless particles.
let ordered = Ordered::new(&particles, is_massive);

// We can build a `RootedOrthtree` from just the massive particles. Note that this is done
// automatically when computing over an [`Ordered`] or [`Reordered`] storage for the Barnes-Hut
// algorithm.
let tree = RootedOrthtree::new(
    ordered.affecting(),
    |p| p.position().to_array(),
    |slice| GravitationalField::centre_of_mass(slice.iter().cloned()),
);

// Do something with the tree.
for (node, data) in std::iter::zip(&tree.get().nodes, &tree.get().data) {
    // ...
}

// Compute the acceleration the massive particles in a tree exert on the massless particles.
let between = Between(ordered.non_affecting(), &tree);
let accelerations = between.barnes_hut(0.5, Acceleration::checked()).collect::<Vec<_>>();

use particular::prelude::*;

use glam::DVec3;

#[derive(Clone, Copy)]
struct Body {
    position: DVec3,
    mass: f64,
}

#[derive(Clone)]
struct GravitationalForce(pub f64);

impl Interaction<Between<&Body, &Body>> for GravitationalForce {
   type Output = DVec3;

    fn compute(&mut self, Between(affected, affecting): Between<&Body, &Body>) -> DVec3 {
        if affected.position == affecting.position {
            return DVec3::ZERO;
        }

        let r = affecting.position - affected.position;
        let l = r.length_squared();
        let f = self.0 * affected.mass * affecting.mass / (l * l.sqrt());

        r * f
    }
}

// Distance: AU, Mass: M☉ (solar mass)
const G: f64 = 4.0 * std::f64::consts::PI * std::f64::consts::PI;
let sun = Body::new(DVec3::ZERO, 1.0);
let earth = Body::new(DVec3::new(1.0, 0.0, 0.0), 3.0027e-6);
let jupiter = Body::new(DVec3::new(5.2, 0.0, 0.0), 0.000954588);

let solar_system = [sun, earth, jupiter];
let mut forces = solar_system.brute_force(GravitationalForce(G));

let sun_earth = GravitationalForce(G).compute(Between(&sun, &earth));
let sun_jupiter = GravitationalForce(G).compute(Between(&sun, &jupiter));
let earth_jupiter = GravitationalForce(G).compute(Between(&earth, &jupiter));
assert_eq!(forces.next(), Some(sun_earth + sun_jupiter));
assert_eq!(forces.next(), Some(-sun_earth + earth_jupiter));
assert_eq!(forces.next(), Some(-sun_jupiter - earth_jupiter));

This project is licensed under either of Apache License, Version 2.0 or MIT license, at your option.

Unless you explicitly state otherwise, any contribution intentionally submitted for inclusion in this project by you, as defined in the Apache 2.0 license, shall be dual licensed as above, without any additional terms or conditions.