diff --git a/README.md b/README.md
deleted file mode 100644
index aef82d8..0000000
--- a/README.md
+++ /dev/null
@@ -1,19 +0,0 @@
-# Stage JSB : Simulation spatiale
-
-## Aide
-
-[Documentation de Bevy](https://docs.rs/bevy/latest/bevy/)
-
-## Commandes
-
-Vérifier la validité du code :
-
-    cargo check
-
-Compiler et lancer le programme :
-
-    cargo run
-
-Compiler et lancer le programme optimisé (plus lent à compiler, mais plus léger et rapide à exécuter) :
-
-    cargo run --release
diff --git a/src/gen.rs b/src/gen.rs
index ee57f63..32243d7 100644
--- a/src/gen.rs
+++ b/src/gen.rs
@@ -2,14 +2,28 @@ use bevy::{prelude::*, render::render_resource::PrimitiveTopology};
 use opensimplex_noise_rs::OpenSimplexNoise;
 use rand::Rng;
 
-pub fn _planet() -> Mesh {
+pub fn planet() -> Mesh {
 	let mut rng = rand::thread_rng();
-	let _simplex = OpenSimplexNoise::new(Some(rng.gen()));
+	let simplex = OpenSimplexNoise::new(Some(rng.gen()));
 	let mut mesh = Mesh::new(PrimitiveTopology::TriangleList);
+	let perimeter: u32 = 1000;
 	mesh.insert_attribute(
 		Mesh::ATTRIBUTE_POSITION,
-		vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]],
+		(0..perimeter)
+			.map(|i| {
+				let mut r = simplex.eval_2d(i as f64 * 0.02, 0.) as f32 * 20.0 + 100.0;
+				r += simplex.eval_2d(i as f64 * 0.05, 10.) as f32 * 10.0;
+				r += simplex.eval_2d(i as f64 * 0.2, 10.) as f32 * 4.0;
+				let a = std::f32::consts::TAU * i as f32 / perimeter as f32;
+				[r * a.cos(), r * a.sin(), 0.]
+			})
+			.chain([[0., 0., 0.]].into_iter())
+			.collect::<Vec<_>>(),
 	);
-	mesh.set_indices(Some(bevy::render::mesh::Indices::U32(vec![0, 1, 2])));
+	let mut triangles = Vec::new();
+	for i in 0..perimeter {
+		triangles.extend_from_slice(&[i, perimeter, (i + 1) % perimeter]);
+	}
+	mesh.set_indices(Some(bevy::render::mesh::Indices::U32(triangles)));
 	mesh
 }
diff --git a/src/main.rs b/src/main.rs
index 8ec9ad3..097c611 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -1,11 +1,14 @@
 mod gen;
+mod quadtree;
 
-use bevy::{prelude::*, sprite::MaterialMesh2dBundle};
+use bevy::{ecs::query::BatchingStrategy, prelude::*, sprite::MaterialMesh2dBundle};
 
 #[global_allocator]
 static ALLOCATOR: cap::Cap<std::alloc::System> =
 	cap::Cap::new(std::alloc::System, 1024 * 1024 * 1024);
 
+static UNIVERSE_POS: (Vec2, Vec2) = (Vec2::new(-1e6, -1e6), Vec2::new(1e6, 1e6));
+
 fn main() {
 	App::new()
 		.add_plugins((
@@ -13,9 +16,16 @@ fn main() {
 			bevy_fps_counter::FpsCounterPlugin,
 			bevy_pancam::PanCamPlugin,
 		))
+		.insert_resource(Constants { g: 6.674e-11 })
 		.add_systems(Startup, setup)
-		.configure_sets(Update, (Set::Demo).chain())
-		.add_systems(Update, (demo_system.in_set(Set::Demo),))
+		.configure_sets(Update, (Set::Force, Set::Apply).chain())
+		.add_systems(
+			Update,
+			(
+				weight_system.in_set(Set::Force),
+				apply_system.in_set(Set::Apply),
+			),
+		)
 		.run();
 }
 
@@ -28,34 +38,150 @@ fn setup(
 		.spawn(Camera2dBundle::default())
 		.insert(bevy_pancam::PanCam::default());
 
-	commands
-		.spawn(Ball {
-			pos: TransformBundle::from(Transform::from_translation(Vec3::new(0., 0., 0.))),
-			visibility: InheritedVisibility::VISIBLE,
-		})
-		.with_children(|parent| {
-			parent.spawn(MaterialMesh2dBundle {
-				mesh: meshes.add(shape::Circle::new(10.).into()).into(),
-				material: materials.add(ColorMaterial::from(Color::ORANGE)),
+	commands.spawn(Planet {
+		mass: Mass(1.988e18),
+		speed: Speed(Vec2::new(0.0, 0.0)),
+		mesh: MaterialMesh2dBundle {
+			mesh: meshes.add(gen::planet()).into(),
+			material: materials.add(ColorMaterial::from(Color::YELLOW)),
+			transform: Transform::from_translation(Vec3::new(0., 0., 0.)),
+			..default()
+		},
+	});
+	commands.spawn(Planet {
+		mass: Mass(5.9736e14),
+		speed: Speed(Vec2::new(0.0, 500.0)),
+		mesh: MaterialMesh2dBundle {
+			mesh: meshes.add(shape::Circle::new(10.).into()).into(),
+			material: materials.add(ColorMaterial::from(Color::BLUE)),
+			transform: Transform::from_translation(Vec3::new(400., 0., 0.)),
+			..default()
+		},
+	});
+	commands.spawn(Planet {
+		mass: Mass(5.9736e14),
+		speed: Speed(Vec2::new(0.0, -500.0)),
+		mesh: MaterialMesh2dBundle {
+			mesh: meshes.add(shape::Circle::new(10.).into()).into(),
+			material: materials.add(ColorMaterial::from(Color::BLUE)),
+			transform: Transform::from_translation(Vec3::new(-400., 0., 0.)),
+			..default()
+		},
+	});
+	for i in 0..4000u32 {
+		commands.spawn(Planet {
+			mass: Mass(1.),
+			speed: Speed(Vec2::new(0.0, -500.0)),
+			mesh: MaterialMesh2dBundle {
+				mesh: meshes.add(shape::Circle::new(5.).into()).into(),
+				material: materials.add(ColorMaterial::from(Color::RED)),
+				transform: Transform::from_translation(Vec3::new(-450. - i as f32 / 4., 0., 0.)),
 				..default()
-			});
+			},
 		});
+	}
 }
 
 #[derive(Clone, Debug, Eq, Hash, PartialEq, SystemSet)]
 enum Set {
-	Demo,
+	Force,
+	Apply,
 }
 
+#[derive(Component)]
+struct Speed(Vec2);
+
+#[derive(Component)]
+struct Mass(f32);
+
 #[derive(Bundle)]
-struct Ball {
-	pos: TransformBundle,
-	visibility: InheritedVisibility,
+struct Planet {
+	mass: Mass,
+	speed: Speed,
+	mesh: MaterialMesh2dBundle<ColorMaterial>,
 }
 
-fn demo_system(query: Query<&Transform>, time: Res<Time>) {
-	println!("Temps écoulé : {}s", time.delta_seconds());
-	for pos in query.iter() {
-		println!("Il y a un objet aux coordonnées {:?}", pos.translation);
+#[derive(Resource)]
+struct Constants {
+	g: f32,
+}
+
+struct Body {
+	mass: f32,
+	pos: Vec2,
+}
+
+impl quadtree::Body for Body {
+	fn mass(&self) -> f32 {
+		self.mass
+	}
+	fn pos(&self) -> Vec2 {
+		self.pos
+	}
+	fn add_mass(&mut self, mass: f32) {
+		self.mass += mass;
 	}
 }
+
+/*fn weight_system(
+	constants: Res<Constants>,
+	query1: Query<(&Transform, &Mass)>,
+	mut query2: Query<(&Transform, &mut Speed)>,
+	time: Res<Time>,
+) {
+	let gdt = constants.g * time.delta_seconds();
+	for (n1_pos, n1_mass) in query1.iter() {
+		for (n2_pos, mut n2_speed) in query2.iter_mut() {
+			let d2 = (n1_pos.translation.x - n2_pos.translation.x).powi(2)
+				+ (n1_pos.translation.y - n2_pos.translation.y).powi(2);
+			if d2 == 0.0 {
+				continue;
+			}
+			let f = n1_mass.0 * gdt / (d2 * d2.sqrt());
+			n2_speed.x += (n1_pos.translation.x - n2_pos.translation.x) * f;
+			n2_speed.y += (n1_pos.translation.y - n2_pos.translation.y) * f;
+		}
+	}
+}*/
+
+fn weight_system(
+	constants: Res<Constants>,
+	mut query: Query<(&Transform, &Mass, &mut Speed)>,
+	time: Res<Time>,
+) {
+	let mut tree = quadtree::Node::new(UNIVERSE_POS);
+	let gdt = constants.g * time.delta_seconds();
+	for (pos, mass, _speed) in query.iter() {
+		tree.add_body(Body {
+			mass: mass.0,
+			pos: pos.translation.xy(),
+		});
+	}
+	query.par_iter_mut().for_each(|(pos, _mass, mut speed)| {
+		speed.0 += gdt * tree.apply(pos.translation.xy(), 0.5);
+	});
+	// let gdt = constants.g * time.delta_seconds();
+	// let mut iter = query.iter_combinations_mut();
+	// while let Some([(n1_pos, n1_mass, mut n1_speed), (n2_pos, n2_mass, mut n2_speed)]) =
+	// 	iter.fetch_next()
+	// {
+	// 	let d2 = (n1_pos.translation.x - n2_pos.translation.x).powi(2)
+	// 		+ (n1_pos.translation.y - n2_pos.translation.y).powi(2);
+	// 	let f = gdt / (d2 * d2.sqrt());
+	// 	n1_speed.x -= (n1_pos.translation.x - n2_pos.translation.x) * f * n2_mass.0;
+	// 	n1_speed.y -= (n1_pos.translation.y - n2_pos.translation.y) * f * n2_mass.0;
+	// 	n2_speed.x += (n1_pos.translation.x - n2_pos.translation.x) * f * n1_mass.0;
+	// 	n2_speed.y += (n1_pos.translation.y - n2_pos.translation.y) * f * n1_mass.0;
+	// }
+}
+
+fn apply_system(mut query: Query<(&mut Transform, &Speed)>, time: Res<Time>) {
+	let dt = time.delta_seconds();
+	query
+		.par_iter_mut()
+		.batching_strategy(BatchingStrategy::fixed(128))
+		.for_each(|(mut pos, speed)| {
+			pos.translation.x += speed.0.x * dt;
+			pos.translation.y += speed.0.y * dt;
+		});
+}
diff --git a/src/quadtree.rs b/src/quadtree.rs
new file mode 100644
index 0000000..c624d15
--- /dev/null
+++ b/src/quadtree.rs
@@ -0,0 +1,160 @@
+use bevy::prelude::*;
+
+pub trait Body {
+	fn mass(&self) -> f32;
+	fn pos(&self) -> Vec2;
+	fn add_mass(&mut self, mass: f32);
+}
+
+pub enum Node<L> {
+	Branch {
+		nodes: Box<[Node<L>; 4]>,
+		center: Vec2,
+		mass: f32,
+		center_of_mass: Vec2,
+		width: f32,
+	},
+	Leaf {
+		body: Option<L>,
+		pos: (Vec2, Vec2),
+	},
+}
+
+impl<L: Body> Node<L> {
+	pub fn new(pos: (Vec2, Vec2)) -> Self {
+		Node::Leaf { body: None, pos }
+		// let center = (pos.1 - pos.0) / 2.0;
+		// Node::Branch {
+		// 	nodes: [
+		// 		Box::new(Node::Leaf {
+		// 			body: None,
+		// 			pos: (pos.0, center),
+		// 		}),
+		// 		Box::new(Node::Leaf {
+		// 			body: None,
+		// 			pos: (Vec2::new(center.x, pos.0.y), Vec2::new(pos.1.x, center.y)),
+		// 		}),
+		// 		Box::new(Node::Leaf {
+		// 			body: None,
+		// 			pos: (Vec2::new(pos.0.x, center.y), Vec2::new(center.x, pos.1.y)),
+		// 		}),
+		// 		Box::new(Node::Leaf {
+		// 			body: None,
+		// 			pos: (center, pos.1),
+		// 		}),
+		// 	],
+		// 	center,
+		// 	mass: 0.0,
+		// 	center_of_mass: center,
+		// 	width: pos.1.x - pos.0.x,
+		// }
+	}
+	pub fn add_body(&mut self, new_body: L) {
+		match self {
+			Node::Branch {
+				nodes,
+				center,
+				mass,
+				center_of_mass,
+				..
+			} => {
+				let new_body_pos = new_body.pos();
+				let new_body_mass = new_body.mass();
+				*center_of_mass = (*center_of_mass * *mass + new_body_mass * new_body_pos)
+					/ (*mass + new_body_mass);
+				*mass += new_body_mass;
+				nodes[if new_body_pos.x < center.x {
+					if new_body_pos.y < center.y {
+						0
+					} else {
+						2
+					}
+				} else {
+					if new_body_pos.y < center.y {
+						1
+					} else {
+						3
+					}
+				}]
+				.add_body(new_body)
+			}
+			Node::Leaf { body, pos } => {
+				if let Some(mut body) = body.take() {
+					if body.pos().distance_squared(new_body.pos()) < 1.0 {
+						body.add_mass(new_body.mass());
+						*self = Node::Leaf {
+							body: Some(body),
+							pos: *pos,
+						};
+						return;
+					}
+					let center = (pos.0 + pos.1) / 2.0;
+					*self = Node::Branch {
+						nodes: Box::new([
+							Node::Leaf {
+								body: None,
+								pos: (pos.0, center),
+							},
+							Node::Leaf {
+								body: None,
+								pos: (Vec2::new(center.x, pos.0.y), Vec2::new(pos.1.x, center.y)),
+							},
+							Node::Leaf {
+								body: None,
+								pos: (Vec2::new(pos.0.x, center.y), Vec2::new(center.x, pos.1.y)),
+							},
+							Node::Leaf {
+								body: None,
+								pos: (center, pos.1),
+							},
+						]),
+						center,
+						mass: 0.0,
+						center_of_mass: center,
+						width: pos.1.x - pos.0.x,
+					};
+					self.add_body(body);
+					self.add_body(new_body)
+				} else {
+					*body = Some(new_body);
+				}
+			}
+		}
+	}
+
+	pub fn apply(&self, on: Vec2, theta: f32) -> Vec2 {
+		match self {
+			Node::Branch {
+				nodes,
+				mass,
+				center_of_mass,
+				width,
+				..
+			} => {
+				if on == *center_of_mass {
+					return Vec2::ZERO;
+				}
+				let dist = on.distance(*center_of_mass);
+				if width / dist < theta {
+					*mass * (*center_of_mass - on) / (dist * dist * dist)
+				} else {
+					nodes[0].apply(on, theta)
+						+ nodes[1].apply(on, theta)
+						+ nodes[2].apply(on, theta)
+						+ nodes[3].apply(on, theta)
+				}
+			}
+			Node::Leaf { body, .. } => {
+				if let Some(body) = body {
+					if on == body.pos() {
+						return Vec2::ZERO;
+					}
+					let dist = on.distance(body.pos());
+					body.mass() * (body.pos() - on) / (dist * dist * dist)
+				} else {
+					Vec2::ZERO
+				}
+			}
+		}
+	}
+}