diff --git a/.gitignore b/.gitignore
index 99be97d9..c2f90191 100644
--- a/.gitignore
+++ b/.gitignore
@@ -16,3 +16,7 @@
 *.blend
 *.blend1
 .idea/
+/images
+/data_analysis
+*.png
+*.gif
diff --git a/Cargo.lock b/Cargo.lock
index 378f59d4..45ba342d 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -3,22 +3,68 @@
 version = 4
 
 [[package]]
-name = "ahash"
-version = "0.7.8"
+name = "accesskit"
+version = "0.18.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "891477e0c6a8957309ee5c45a6368af3ae14bb510732d2684ffa19af310920f9"
+checksum = "becf0eb5215b6ecb0a739c31c21bd83c4f326524c9b46b7e882d77559b60a529"
+
+[[package]]
+name = "aho-corasick"
+version = "1.1.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8e60d3430d3a69478ad0993f19238d2df97c507009a52b3c10addcd7f6bcb916"
 dependencies = [
- "getrandom 0.2.6",
- "once_cell",
- "version_check",
+ "memchr",
+]
+
+[[package]]
+name = "allocator-api2"
+version = "0.2.21"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "683d7910e743518b0e34f1186f92494becacb047c7b6bf616c96772180fef923"
+
+[[package]]
+name = "android-activity"
+version = "0.6.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ef6978589202a00cd7e118380c448a08b6ed394c3a8df3a430d0898e3a42d046"
+dependencies = [
+ "android-properties",
+ "bitflags 2.9.1",
+ "cc",
+ "cesu8",
+ "jni",
+ "jni-sys",
+ "libc",
+ "log",
+ "ndk",
+ "ndk-context",
+ "ndk-sys 0.6.0+11769913",
+ "num_enum",
+ "thiserror 1.0.69",
 ]
 
+[[package]]
+name = "android-properties"
+version = "0.2.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "fc7eb209b1518d6bb87b283c20095f5228ecda460da70b44f0802523dea6da04"
+
 [[package]]
 name = "android_log-sys"
 version = "0.3.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "84521a3cf562bc62942e294181d9eef17eb38ceb8c68677bc49f144e4c3d4f8d"
 
+[[package]]
+name = "android_system_properties"
+version = "0.1.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "819e7219dbd41043ac279b19830f2efc897156490d7fd6ea916720117ee66311"
+dependencies = [
+ "libc",
+]
+
 [[package]]
 name = "anes"
 version = "0.1.6"
@@ -49,12 +95,27 @@ dependencies = [
  "num-traits",
 ]
 
+[[package]]
+name = "arrayref"
+version = "0.3.9"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "76a2e8124351fda1ef8aaaa3bbd7ebbcb486bbcd4225aca0aa0d84bb2db8fecb"
+
 [[package]]
 name = "arrayvec"
 version = "0.7.6"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "7c02d123df017efcdfbd739ef81735b36c5ba83ec3c59c80a9d7ecc718f92e50"
 
+[[package]]
+name = "ash"
+version = "0.38.0+1.3.281"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0bb44936d800fea8f016d7f2311c6a4f97aebd5dc86f09906139ec848cf3a46f"
+dependencies = [
+ "libloading",
+]
+
 [[package]]
 name = "assert_approx_eq"
 version = "1.1.0"
@@ -72,6 +133,30 @@ dependencies = [
  "syn 2.0.104",
 ]
 
+[[package]]
+name = "async-broadcast"
+version = "0.7.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "435a87a52755b8f27fcf321ac4f04b2802e337c8c4872923137471ec39c37532"
+dependencies = [
+ "event-listener",
+ "event-listener-strategy",
+ "futures-core",
+ "pin-project-lite",
+]
+
+[[package]]
+name = "async-channel"
+version = "2.5.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "924ed96dd52d1b75e9c1a3e6275715fd320f5f9439fb5a4a11fa51f4221158d2"
+dependencies = [
+ "concurrent-queue",
+ "event-listener-strategy",
+ "futures-core",
+ "pin-project-lite",
+]
+
 [[package]]
 name = "async-executor"
 version = "1.13.2"
@@ -86,6 +171,28 @@ dependencies = [
  "slab",
 ]
 
+[[package]]
+name = "async-fs"
+version = "2.1.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "09f7e37c0ed80b2a977691c47dae8625cfb21e205827106c64f7c588766b2e50"
+dependencies = [
+ "async-lock",
+ "blocking",
+ "futures-lite",
+]
+
+[[package]]
+name = "async-lock"
+version = "3.4.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5fd03604047cee9b6ce9de9f70c6cd540a0520c813cbd49bae61f33ab80ed1dc"
+dependencies = [
+ "event-listener",
+ "event-listener-strategy",
+ "pin-project-lite",
+]
+
 [[package]]
 name = "async-task"
 version = "4.7.1"
@@ -99,18 +206,21 @@ dependencies = [
 name = "atomecs"
 version = "0.8.1"
 dependencies = [
+ "approx",
  "assert_approx_eq",
  "atomecs-derive",
  "bevy",
+ "bytemuck",
  "byteorder",
  "criterion",
  "csv",
- "hashbrown 0.12.3",
+ "hashbrown",
  "multimap",
  "nalgebra",
- "rand 0.9.1",
+ "rand 0.9.2",
  "rand_distr 0.5.1",
  "rayon",
+ "regex",
  "serde",
  "serde_arrays",
  "serde_json",
@@ -136,12 +246,28 @@ dependencies = [
  "portable-atomic",
 ]
 
+[[package]]
+name = "atomicow"
+version = "1.1.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f52e8890bb9844440d0c412fa74b67fd2f14e85248b6e00708059b6da9e5f8bf"
+dependencies = [
+ "portable-atomic",
+ "portable-atomic-util",
+]
+
 [[package]]
 name = "autocfg"
 version = "1.1.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa"
 
+[[package]]
+name = "base64"
+version = "0.21.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9d297deb1925b89f2ccc13d7635fa0714f12c87adce1c75356b39ca9b7178567"
+
 [[package]]
 name = "bevy"
 version = "0.16.1"
@@ -151,6 +277,19 @@ dependencies = [
  "bevy_internal",
 ]
 
+[[package]]
+name = "bevy_a11y"
+version = "0.16.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ed3561712cf49074d89e9989bfc2e6c6add5d33288f689db9a0c333300d2d004"
+dependencies = [
+ "accesskit",
+ "bevy_app",
+ "bevy_derive",
+ "bevy_ecs",
+ "bevy_reflect",
+]
+
 [[package]]
 name = "bevy_app"
 version = "0.16.1"
@@ -168,12 +307,80 @@ dependencies = [
  "ctrlc",
  "downcast-rs",
  "log",
- "thiserror",
+ "thiserror 2.0.12",
  "variadics_please",
  "wasm-bindgen",
  "web-sys",
 ]
 
+[[package]]
+name = "bevy_asset"
+version = "0.16.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f56111d9b88d8649f331a667d9d72163fb26bd09518ca16476d238653823db1e"
+dependencies = [
+ "async-broadcast",
+ "async-fs",
+ "async-lock",
+ "atomicow",
+ "bevy_app",
+ "bevy_asset_macros",
+ "bevy_ecs",
+ "bevy_platform",
+ "bevy_reflect",
+ "bevy_tasks",
+ "bevy_utils",
+ "bevy_window",
+ "bitflags 2.9.1",
+ "blake3",
+ "crossbeam-channel",
+ "derive_more",
+ "disqualified",
+ "downcast-rs",
+ "either",
+ "futures-io",
+ "futures-lite",
+ "js-sys",
+ "parking_lot",
+ "ron",
+ "serde",
+ "stackfuture",
+ "thiserror 2.0.12",
+ "tracing",
+ "uuid",
+ "wasm-bindgen",
+ "wasm-bindgen-futures",
+ "web-sys",
+]
+
+[[package]]
+name = "bevy_asset_macros"
+version = "0.16.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a4cca3e67c0ec760d8889d42293d987ce5da92eaf9c592bf5d503728a63b276d"
+dependencies = [
+ "bevy_macro_utils",
+ "proc-macro2",
+ "quote",
+ "syn 2.0.104",
+]
+
+[[package]]
+name = "bevy_color"
+version = "0.16.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5c101cbe1e26b8d701eb77263b14346e2e0cbbd2a6e254b9b1aead814e5ca8d3"
+dependencies = [
+ "bevy_math",
+ "bevy_reflect",
+ "bytemuck",
+ "derive_more",
+ "encase",
+ "serde",
+ "thiserror 2.0.12",
+ "wgpu-types",
+]
+
 [[package]]
 name = "bevy_derive"
 version = "0.16.1"
@@ -199,6 +406,7 @@ dependencies = [
  "bevy_utils",
  "const-fnv1a-hash",
  "log",
+ "serde",
 ]
 
 [[package]]
@@ -225,7 +433,7 @@ dependencies = [
  "nonmax",
  "serde",
  "smallvec",
- "thiserror",
+ "thiserror 2.0.12",
  "variadics_please",
 ]
 
@@ -241,6 +449,42 @@ dependencies = [
  "syn 2.0.104",
 ]
 
+[[package]]
+name = "bevy_encase_derive"
+version = "0.16.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8148f4edee470a2ea5cad010184c492a4c94c36d7a7158ea28e134ea87f274ab"
+dependencies = [
+ "bevy_macro_utils",
+ "encase_derive_impl",
+]
+
+[[package]]
+name = "bevy_image"
+version = "0.16.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "65e6e900cfecadbc3149953169e36b9e26f922ed8b002d62339d8a9dc6129328"
+dependencies = [
+ "bevy_app",
+ "bevy_asset",
+ "bevy_color",
+ "bevy_math",
+ "bevy_platform",
+ "bevy_reflect",
+ "bevy_utils",
+ "bitflags 2.9.1",
+ "bytemuck",
+ "futures-lite",
+ "guillotiere",
+ "half",
+ "image",
+ "rectangle-pack",
+ "serde",
+ "thiserror 2.0.12",
+ "tracing",
+ "wgpu-types",
+]
+
 [[package]]
 name = "bevy_input"
 version = "0.16.1"
@@ -255,7 +499,23 @@ dependencies = [
  "bevy_utils",
  "derive_more",
  "log",
- "thiserror",
+ "thiserror 2.0.12",
+]
+
+[[package]]
+name = "bevy_input_focus"
+version = "0.16.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2e2d079fda74d1416e0a57dac29ea2b79ff77f420cd6b87f833d3aa29a46bc4d"
+dependencies = [
+ "bevy_app",
+ "bevy_ecs",
+ "bevy_input",
+ "bevy_math",
+ "bevy_reflect",
+ "bevy_window",
+ "log",
+ "thiserror 2.0.12",
 ]
 
 [[package]]
@@ -264,20 +524,27 @@ version = "0.16.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "857da8785678fde537d02944cd20dec9cafb7d4c447efe15f898dc60e733cacd"
 dependencies = [
+ "bevy_a11y",
  "bevy_app",
+ "bevy_asset",
+ "bevy_color",
  "bevy_derive",
  "bevy_diagnostic",
  "bevy_ecs",
  "bevy_input",
+ "bevy_input_focus",
  "bevy_log",
  "bevy_math",
  "bevy_platform",
  "bevy_ptr",
  "bevy_reflect",
+ "bevy_render",
+ "bevy_state",
  "bevy_tasks",
  "bevy_time",
  "bevy_transform",
  "bevy_utils",
+ "bevy_window",
 ]
 
 [[package]]
@@ -316,19 +583,54 @@ version = "0.16.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "68553e0090fe9c3ba066c65629f636bd58e4ebd9444fdba097b91af6cd3e243f"
 dependencies = [
+ "approx",
  "bevy_reflect",
  "derive_more",
- "glam",
+ "glam 0.29.3",
  "itertools 0.14.0",
  "libm",
  "rand 0.8.5",
  "rand_distr 0.4.3",
  "serde",
  "smallvec",
- "thiserror",
+ "thiserror 2.0.12",
  "variadics_please",
 ]
 
+[[package]]
+name = "bevy_mesh"
+version = "0.16.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b10399c7027001edbc0406d7d0198596b1f07206c1aae715274106ba5bdcac40"
+dependencies = [
+ "bevy_asset",
+ "bevy_derive",
+ "bevy_ecs",
+ "bevy_image",
+ "bevy_math",
+ "bevy_mikktspace",
+ "bevy_platform",
+ "bevy_reflect",
+ "bevy_transform",
+ "bevy_utils",
+ "bitflags 2.9.1",
+ "bytemuck",
+ "hexasphere",
+ "serde",
+ "thiserror 2.0.12",
+ "tracing",
+ "wgpu-types",
+]
+
+[[package]]
+name = "bevy_mikktspace"
+version = "0.16.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8bb60c753b968a2de0fd279b76a3d19517695e771edb4c23575c7f92156315de"
+dependencies = [
+ "glam 0.29.3",
+]
+
 [[package]]
 name = "bevy_platform"
 version = "0.16.1"
@@ -339,7 +641,7 @@ dependencies = [
  "critical-section",
  "foldhash",
  "getrandom 0.2.6",
- "hashbrown 0.15.4",
+ "hashbrown",
  "portable-atomic",
  "portable-atomic-util",
  "serde",
@@ -369,11 +671,11 @@ dependencies = [
  "downcast-rs",
  "erased-serde",
  "foldhash",
- "glam",
+ "glam 0.29.3",
  "serde",
  "smallvec",
  "smol_str",
- "thiserror",
+ "thiserror 2.0.12",
  "uuid",
  "variadics_please",
  "wgpu-types",
@@ -393,52 +695,149 @@ dependencies = [
 ]
 
 [[package]]
-name = "bevy_tasks"
+name = "bevy_render"
 version = "0.16.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "5b674242641cab680688fc3b850243b351c1af49d4f3417a576debd6cca8dcf5"
+checksum = "ef91fed1f09405769214b99ebe4390d69c1af5cdd27967deae9135c550eb1667"
 dependencies = [
- "async-executor",
- "async-task",
- "atomic-waker",
+ "async-channel",
+ "bevy_app",
+ "bevy_asset",
+ "bevy_color",
+ "bevy_derive",
+ "bevy_diagnostic",
+ "bevy_ecs",
+ "bevy_encase_derive",
+ "bevy_image",
+ "bevy_math",
+ "bevy_mesh",
  "bevy_platform",
- "cfg-if",
- "crossbeam-queue",
+ "bevy_reflect",
+ "bevy_render_macros",
+ "bevy_tasks",
+ "bevy_time",
+ "bevy_transform",
+ "bevy_utils",
+ "bevy_window",
+ "bitflags 2.9.1",
+ "bytemuck",
+ "codespan-reporting",
  "derive_more",
- "futures-channel",
+ "downcast-rs",
+ "encase",
+ "fixedbitset",
  "futures-lite",
- "heapless",
- "pin-project",
- "wasm-bindgen-futures",
+ "image",
+ "indexmap",
+ "js-sys",
+ "naga",
+ "naga_oil",
+ "nonmax",
+ "offset-allocator",
+ "send_wrapper",
+ "serde",
+ "smallvec",
+ "thiserror 2.0.12",
+ "tracing",
+ "variadics_please",
+ "wasm-bindgen",
+ "web-sys",
+ "wgpu",
 ]
 
 [[package]]
-name = "bevy_time"
+name = "bevy_render_macros"
 version = "0.16.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "bc98eb356c75be04fbbc77bb3d8ffa24c8bacd99f76111cee23d444be6ac8c9c"
+checksum = "abd42cf6c875bcf38da859f8e731e119a6aff190d41dd0a1b6000ad57cf2ed3d"
+dependencies = [
+ "bevy_macro_utils",
+ "proc-macro2",
+ "quote",
+ "syn 2.0.104",
+]
+
+[[package]]
+name = "bevy_state"
+version = "0.16.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "155d3cd97b900539008cdcaa702f88b724d94b08977b8e591a32536ce66faa8c"
 dependencies = [
  "bevy_app",
  "bevy_ecs",
  "bevy_platform",
  "bevy_reflect",
+ "bevy_state_macros",
+ "bevy_utils",
  "log",
+ "variadics_please",
 ]
 
 [[package]]
-name = "bevy_transform"
+name = "bevy_state_macros"
 version = "0.16.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "df218e440bb9a19058e1b80a68a031c887bcf7bd3a145b55f361359a2fa3100d"
+checksum = "2481c1304fd2a1851a0d4cb63a1ce6421ae40f3f0117cbc9882963ee4c9bb609"
 dependencies = [
- "bevy_app",
+ "bevy_macro_utils",
+ "proc-macro2",
+ "quote",
+ "syn 2.0.104",
+]
+
+[[package]]
+name = "bevy_tasks"
+version = "0.16.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5b674242641cab680688fc3b850243b351c1af49d4f3417a576debd6cca8dcf5"
+dependencies = [
+ "async-channel",
+ "async-executor",
+ "async-task",
+ "atomic-waker",
+ "bevy_platform",
+ "cfg-if",
+ "concurrent-queue",
+ "crossbeam-queue",
+ "derive_more",
+ "futures-channel",
+ "futures-lite",
+ "heapless",
+ "pin-project",
+ "wasm-bindgen-futures",
+]
+
+[[package]]
+name = "bevy_time"
+version = "0.16.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "bc98eb356c75be04fbbc77bb3d8ffa24c8bacd99f76111cee23d444be6ac8c9c"
+dependencies = [
+ "bevy_app",
  "bevy_ecs",
+ "bevy_platform",
+ "bevy_reflect",
+ "crossbeam-channel",
+ "log",
+ "serde",
+]
+
+[[package]]
+name = "bevy_transform"
+version = "0.16.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "df218e440bb9a19058e1b80a68a031c887bcf7bd3a145b55f361359a2fa3100d"
+dependencies = [
+ "bevy_app",
+ "bevy_ecs",
+ "bevy_log",
  "bevy_math",
  "bevy_reflect",
  "bevy_tasks",
+ "bevy_utils",
  "derive_more",
  "serde",
- "thiserror",
+ "thiserror 2.0.12",
 ]
 
 [[package]]
@@ -451,6 +850,26 @@ dependencies = [
  "thread_local",
 ]
 
+[[package]]
+name = "bevy_window"
+version = "0.16.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "df7e8ad0c17c3cc23ff5566ae2905c255e6986037fb041f74c446216f5c38431"
+dependencies = [
+ "android-activity",
+ "bevy_app",
+ "bevy_ecs",
+ "bevy_input",
+ "bevy_math",
+ "bevy_platform",
+ "bevy_reflect",
+ "bevy_utils",
+ "log",
+ "raw-window-handle",
+ "serde",
+ "smol_str",
+]
+
 [[package]]
 name = "bindgen"
 version = "0.70.1"
@@ -471,6 +890,36 @@ dependencies = [
  "syn 2.0.104",
 ]
 
+[[package]]
+name = "bit-set"
+version = "0.5.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0700ddab506f33b20a03b13996eccd309a48e5ff77d0d95926aa0210fb4e95f1"
+dependencies = [
+ "bit-vec 0.6.3",
+]
+
+[[package]]
+name = "bit-set"
+version = "0.8.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "08807e080ed7f9d5433fa9b275196cfc35414f66a0c79d864dc51a0d825231a3"
+dependencies = [
+ "bit-vec 0.8.0",
+]
+
+[[package]]
+name = "bit-vec"
+version = "0.6.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "349f9b6a179ed607305526ca489b34ad0a41aed5f7980fa90eb03160b69598fb"
+
+[[package]]
+name = "bit-vec"
+version = "0.8.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5e764a1d40d510daf35e07be9eb06e75770908c27d411ee6c92109c9840eaaf7"
+
 [[package]]
 name = "bitflags"
 version = "1.3.2"
@@ -486,6 +935,38 @@ dependencies = [
  "serde",
 ]
 
+[[package]]
+name = "blake3"
+version = "1.8.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3888aaa89e4b2a40fca9848e400f6a658a5a3978de7be858e209cafa8be9a4a0"
+dependencies = [
+ "arrayref",
+ "arrayvec",
+ "cc",
+ "cfg-if",
+ "constant_time_eq",
+]
+
+[[package]]
+name = "block"
+version = "0.1.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0d8c1fef690941d3e7788d328517591fecc684c084084702d6ff1641e993699a"
+
+[[package]]
+name = "blocking"
+version = "1.6.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e83f8d02be6967315521be875afa792a316e28d57b5a2d401897e2a7921b7f21"
+dependencies = [
+ "async-channel",
+ "async-task",
+ "futures-io",
+ "futures-lite",
+ "piper",
+]
+
 [[package]]
 name = "bumpalo"
 version = "3.9.1"
@@ -494,9 +975,23 @@ checksum = "a4a45a46ab1f2412e53d3a0ade76ffad2025804294569aae387231a0cd6e0899"
 
 [[package]]
 name = "bytemuck"
-version = "1.9.1"
+version = "1.23.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "cdead85bdec19c194affaeeb670c0e41fe23de31459efd1c174d049269cf02cc"
+checksum = "3995eaeebcdf32f91f980d360f78732ddc061097ab4e39991ae7a6ace9194677"
+dependencies = [
+ "bytemuck_derive",
+]
+
+[[package]]
+name = "bytemuck_derive"
+version = "1.10.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4f154e572231cb6ba2bd1176980827e3d5dc04cc183a75dea38109fbdd672d29"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.104",
+]
 
 [[package]]
 name = "byteorder"
@@ -504,6 +999,18 @@ version = "1.5.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b"
 
+[[package]]
+name = "byteorder-lite"
+version = "0.1.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8f1fe948ff07f4bd06c30984e69f5b4899c516a3ef74f34df92a2df2ab535495"
+
+[[package]]
+name = "bytes"
+version = "1.10.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d71b6127be86fdcfddb610f7182ac57211d4b18a3e9c82eb2d17662f2227ad6a"
+
 [[package]]
 name = "cast"
 version = "0.3.0"
@@ -512,9 +1019,20 @@ checksum = "37b2a672a2cb129a2e41c10b1224bb368f9f37a2b16b612598138befd7b37eb5"
 
 [[package]]
 name = "cc"
-version = "1.0.73"
+version = "1.2.32"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "2fff2a6927b3bb87f9595d67196a70493f627687a71d87a0d692242c33f58c11"
+checksum = "2352e5597e9c544d5e6d9c95190d5d27738ade584fa8db0a16e130e5c2b5296e"
+dependencies = [
+ "jobserver",
+ "libc",
+ "shlex",
+]
+
+[[package]]
+name = "cesu8"
+version = "1.1.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6d43a04d8753f35258c91f8ec639f792891f748a1edbd759cf1dcea3382ad83c"
 
 [[package]]
 name = "cexpr"
@@ -600,6 +1118,26 @@ version = "0.7.5"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "b94f61472cee1439c0b966b47e3aca9ae07e45d070759512cd390ea2bebc6675"
 
+[[package]]
+name = "codespan-reporting"
+version = "0.11.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3538270d33cc669650c4b093848450d380def10c331d38c768e34cac80576e6e"
+dependencies = [
+ "termcolor",
+ "unicode-width",
+]
+
+[[package]]
+name = "combine"
+version = "4.6.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ba5a308b75df32fe02788e748662718f03fde005016435c444eea572398219fd"
+dependencies = [
+ "bytes",
+ "memchr",
+]
+
 [[package]]
 name = "concurrent-queue"
 version = "2.5.0"
@@ -626,11 +1164,65 @@ version = "1.1.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "32b13ea120a812beba79e34316b3942a857c86ec1593cb34f27bb28272ce2cca"
 
+[[package]]
+name = "const_panic"
+version = "0.2.13"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b98d1483e98c9d67f341ab4b3915cfdc54740bd6f5cccc9226ee0535d86aa8fb"
+
+[[package]]
+name = "const_soft_float"
+version = "0.1.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "87ca1caa64ef4ed453e68bb3db612e51cf1b2f5b871337f0fcab1c8f87cc3dff"
+
+[[package]]
+name = "constant_time_eq"
+version = "0.3.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7c74b8349d32d297c9134b8c88677813a227df8f779daa29bfc29c183fe3dca6"
+
+[[package]]
+name = "constgebra"
+version = "0.1.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e1aaf9b65849a68662ac6c0810c8893a765c960b907dd7cfab9c4a50bf764fbc"
+dependencies = [
+ "const_soft_float",
+]
+
+[[package]]
+name = "core-foundation"
+version = "0.9.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "91e195e091a93c46f7102ec7818a2aa394e1e1771c3ab4825963fa03e45afb8f"
+dependencies = [
+ "core-foundation-sys",
+ "libc",
+]
+
+[[package]]
+name = "core-foundation-sys"
+version = "0.8.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "773648b94d0e5d620f64f280777445740e61fe701025087ec8b57f45c791888b"
+
+[[package]]
+name = "core-graphics-types"
+version = "0.1.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "45390e6114f68f718cc7a830514a96f903cccd70d02a8f6d9f643ac4ba45afaf"
+dependencies = [
+ "bitflags 1.3.2",
+ "core-foundation",
+ "libc",
+]
+
 [[package]]
 name = "criterion"
-version = "0.6.0"
+version = "0.7.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "3bf7af66b0989381bd0be551bd7cc91912a655a58c6918420c9527b1fd8b4679"
+checksum = "e1c047a62b0cc3e145fa84415a3191f628e980b194c2755aa12300a4e6cbd928"
 dependencies = [
  "anes",
  "cast",
@@ -651,12 +1243,12 @@ dependencies = [
 
 [[package]]
 name = "criterion-plot"
-version = "0.5.0"
+version = "0.6.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "6b50826342786a51a89e2da3a28f1c32b06e387201bc2d19791f622c673706b1"
+checksum = "9b1bcc0dc7dfae599d84ad0b1a55f80cde8af3725da8313b528da95ef783e338"
 dependencies = [
  "cast",
- "itertools 0.10.3",
+ "itertools 0.13.0",
 ]
 
 [[package]]
@@ -756,6 +1348,12 @@ dependencies = [
  "windows-sys 0.59.0",
 ]
 
+[[package]]
+name = "data-encoding"
+version = "2.9.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2a2330da5de22e8a3cb63252ce2abb30116bf5265e89c0e01bc17015ce30a476"
+
 [[package]]
 name = "derive_more"
 version = "1.0.0"
@@ -783,6 +1381,15 @@ version = "1.0.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "c9c272297e804878a2a4b707cfcfc6d2328b5bb936944613b4fdf2b9269afdfd"
 
+[[package]]
+name = "document-features"
+version = "0.2.11"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "95249b50c6c185bee49034bcb378a49dc2b5dff0be90ff6616d31d64febab05d"
+dependencies = [
+ "litrs",
+]
+
 [[package]]
 name = "downcast-rs"
 version = "2.0.1"
@@ -791,9 +1398,41 @@ checksum = "ea8a8b81cacc08888170eef4d13b775126db426d0b348bee9d18c2c1eaf123cf"
 
 [[package]]
 name = "either"
-version = "1.6.1"
+version = "1.15.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "e78d4f1cc4ae33bbfc157ed5d5a5ef3bc29227303d595861deb238fcec4e9457"
+checksum = "48c757948c5ede0e46177b7add2e67155f70e33c07fea8284df6576da70b3719"
+
+[[package]]
+name = "encase"
+version = "0.10.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b0a05902cf601ed11d564128448097b98ebe3c6574bd7b6a653a3d56d54aa020"
+dependencies = [
+ "const_panic",
+ "encase_derive",
+ "glam 0.29.3",
+ "thiserror 1.0.69",
+]
+
+[[package]]
+name = "encase_derive"
+version = "0.10.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "181d475b694e2dd56ae919ce7699d344d1fd259292d590c723a50d1189a2ea85"
+dependencies = [
+ "encase_derive_impl",
+]
+
+[[package]]
+name = "encase_derive_impl"
+version = "0.10.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f97b51c5cc57ef7c5f7a0c57c250251c49ee4c28f819f87ac32f4aceabc36792"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.104",
+]
 
 [[package]]
 name = "equivalent"
@@ -811,6 +1450,36 @@ dependencies = [
  "typeid",
 ]
 
+[[package]]
+name = "euclid"
+version = "0.22.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b52c2ef4a78da0ba68fbe1fd920627411096d2ac478f7f4c9f3a54ba6705bade"
+dependencies = [
+ "num-traits",
+]
+
+[[package]]
+name = "event-listener"
+version = "5.4.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e13b66accf52311f30a0db42147dadea9850cb48cd070028831ae5f5d4b856ab"
+dependencies = [
+ "concurrent-queue",
+ "parking",
+ "pin-project-lite",
+]
+
+[[package]]
+name = "event-listener-strategy"
+version = "0.5.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8be9f3dfaaffdae2972880079a491a1a8bb7cbed0b8dd7a347f668b4150a3b93"
+dependencies = [
+ "event-listener",
+ "pin-project-lite",
+]
+
 [[package]]
 name = "fastrand"
 version = "2.3.0"
@@ -818,92 +1487,303 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "37909eebbb50d72f9059c3b6d82c0463f2ff062c9e95845c43a6c9c0355411be"
 
 [[package]]
-name = "fixedbitset"
-version = "0.5.7"
+name = "fixedbitset"
+version = "0.5.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1d674e81391d1e1ab681a28d99df07927c6d4aa5b027d7da16ba32d1d21ecd99"
+
+[[package]]
+name = "foldhash"
+version = "0.1.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d9c4f5dac5e15c24eb999c26181a6ca40b39fe946cbe4c263c7209467bc83af2"
+
+[[package]]
+name = "foreign-types"
+version = "0.5.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d737d9aa519fb7b749cbc3b962edcf310a8dd1f4b67c91c4f83975dbdd17d965"
+dependencies = [
+ "foreign-types-macros",
+ "foreign-types-shared",
+]
+
+[[package]]
+name = "foreign-types-macros"
+version = "0.2.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1a5c6c585bc94aaf2c7b51dd4c2ba22680844aba4c687be581871a6f518c5742"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.104",
+]
+
+[[package]]
+name = "foreign-types-shared"
+version = "0.3.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "aa9a19cbb55df58761df49b23516a86d432839add4af60fc256da840f66ed35b"
+
+[[package]]
+name = "futures-channel"
+version = "0.3.21"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c3083ce4b914124575708913bca19bfe887522d6e2e6d0952943f5eac4a74010"
+dependencies = [
+ "futures-core",
+]
+
+[[package]]
+name = "futures-core"
+version = "0.3.21"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0c09fd04b7e4073ac7156a9539b57a484a8ea920f79c7c675d05d289ab6110d3"
+
+[[package]]
+name = "futures-io"
+version = "0.3.31"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9e5c1b78ca4aae1ac06c48a526a655760685149f0d465d21f37abfe57ce075c6"
+
+[[package]]
+name = "futures-lite"
+version = "2.6.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f5edaec856126859abb19ed65f39e90fea3a9574b9707f13539acf4abf7eb532"
+dependencies = [
+ "fastrand",
+ "futures-core",
+ "futures-io",
+ "parking",
+ "pin-project-lite",
+]
+
+[[package]]
+name = "getrandom"
+version = "0.2.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9be70c98951c83b8d2f8f60d7065fa6d5146873094452a1008da8c2f1e4205ad"
+dependencies = [
+ "cfg-if",
+ "js-sys",
+ "libc",
+ "wasi 0.10.2+wasi-snapshot-preview1",
+ "wasm-bindgen",
+]
+
+[[package]]
+name = "getrandom"
+version = "0.3.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "26145e563e54f2cadc477553f1ec5ee650b00862f0a58bcd12cbdc5f0ea2d2f4"
+dependencies = [
+ "cfg-if",
+ "libc",
+ "r-efi",
+ "wasi 0.14.2+wasi-0.2.4",
+]
+
+[[package]]
+name = "gl_generator"
+version = "0.14.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1a95dfc23a2b4a9a2f5ab41d194f8bfda3cabec42af4e39f08c339eb2a0c124d"
+dependencies = [
+ "khronos_api",
+ "log",
+ "xml-rs",
+]
+
+[[package]]
+name = "glam"
+version = "0.14.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "333928d5eb103c5d4050533cec0384302db6be8ef7d3cebd30ec6a35350353da"
+
+[[package]]
+name = "glam"
+version = "0.15.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3abb554f8ee44336b72d522e0a7fe86a29e09f839a36022fa869a7dfe941a54b"
+
+[[package]]
+name = "glam"
+version = "0.16.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4126c0479ccf7e8664c36a2d719f5f2c140fbb4f9090008098d2c291fa5b3f16"
+
+[[package]]
+name = "glam"
+version = "0.17.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e01732b97afd8508eee3333a541b9f7610f454bb818669e66e90f5f57c93a776"
+
+[[package]]
+name = "glam"
+version = "0.18.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "525a3e490ba77b8e326fb67d4b44b4bd2f920f44d4cc73ccec50adc68e3bee34"
+
+[[package]]
+name = "glam"
+version = "0.19.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2b8509e6791516e81c1a630d0bd7fbac36d2fa8712a9da8662e716b52d5051ca"
+
+[[package]]
+name = "glam"
+version = "0.20.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f43e957e744be03f5801a55472f593d43fabdebf25a4585db250f04d86b1675f"
+
+[[package]]
+name = "glam"
+version = "0.21.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "518faa5064866338b013ff9b2350dc318e14cc4fcd6cb8206d7e7c9886c98815"
+
+[[package]]
+name = "glam"
+version = "0.22.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "12f597d56c1bd55a811a1be189459e8fad2bbc272616375602443bdfb37fa774"
+
+[[package]]
+name = "glam"
+version = "0.23.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8e4afd9ad95555081e109fe1d21f2a30c691b5f0919c67dfa690a2e1eb6bd51c"
+
+[[package]]
+name = "glam"
+version = "0.24.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b5418c17512bdf42730f9032c74e1ae39afc408745ebb2acf72fbc4691c17945"
+
+[[package]]
+name = "glam"
+version = "0.25.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "151665d9be52f9bb40fc7966565d39666f2d1e69233571b71b87791c7e0528b3"
+
+[[package]]
+name = "glam"
+version = "0.27.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9e05e7e6723e3455f4818c7b26e855439f7546cf617ef669d1adedb8669e5cb9"
+
+[[package]]
+name = "glam"
+version = "0.28.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "779ae4bf7e8421cf91c0b3b64e7e8b40b862fba4d393f59150042de7c4965a94"
+
+[[package]]
+name = "glam"
+version = "0.29.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8babf46d4c1c9d92deac9f7be466f76dfc4482b6452fc5024b5e8daf6ffeb3ee"
+dependencies = [
+ "bytemuck",
+ "libm",
+ "rand 0.8.5",
+ "serde",
+]
+
+[[package]]
+name = "glam"
+version = "0.30.5"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "1d674e81391d1e1ab681a28d99df07927c6d4aa5b027d7da16ba32d1d21ecd99"
+checksum = "f2d1aab06663bdce00d6ca5e5ed586ec8d18033a771906c993a1e3755b368d85"
 
 [[package]]
-name = "foldhash"
-version = "0.1.5"
+name = "glob"
+version = "0.3.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "d9c4f5dac5e15c24eb999c26181a6ca40b39fe946cbe4c263c7209467bc83af2"
+checksum = "a8d1add55171497b4705a648c6b583acafb01d58050a51727785f0b2c8e0a2b2"
 
 [[package]]
-name = "futures-channel"
-version = "0.3.21"
+name = "glow"
+version = "0.16.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "c3083ce4b914124575708913bca19bfe887522d6e2e6d0952943f5eac4a74010"
+checksum = "c5e5ea60d70410161c8bf5da3fdfeaa1c72ed2c15f8bbb9d19fe3a4fad085f08"
 dependencies = [
- "futures-core",
+ "js-sys",
+ "slotmap",
+ "wasm-bindgen",
+ "web-sys",
 ]
 
 [[package]]
-name = "futures-core"
-version = "0.3.21"
+name = "glutin_wgl_sys"
+version = "0.6.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "0c09fd04b7e4073ac7156a9539b57a484a8ea920f79c7c675d05d289ab6110d3"
+checksum = "2c4ee00b289aba7a9e5306d57c2d05499b2e5dc427f84ac708bd2c090212cf3e"
+dependencies = [
+ "gl_generator",
+]
 
 [[package]]
-name = "futures-io"
-version = "0.3.21"
+name = "gpu-alloc"
+version = "0.6.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "fc4045962a5a5e935ee2fdedaa4e08284547402885ab326734432bed5d12966b"
+checksum = "fbcd2dba93594b227a1f57ee09b8b9da8892c34d55aa332e034a228d0fe6a171"
+dependencies = [
+ "bitflags 2.9.1",
+ "gpu-alloc-types",
+]
 
 [[package]]
-name = "futures-lite"
-version = "2.6.0"
+name = "gpu-alloc-types"
+version = "0.3.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "f5edaec856126859abb19ed65f39e90fea3a9574b9707f13539acf4abf7eb532"
+checksum = "98ff03b468aa837d70984d55f5d3f846f6ec31fe34bbb97c4f85219caeee1ca4"
 dependencies = [
- "fastrand",
- "futures-core",
- "futures-io",
- "parking",
- "pin-project-lite",
+ "bitflags 2.9.1",
 ]
 
 [[package]]
-name = "getrandom"
-version = "0.2.6"
+name = "gpu-allocator"
+version = "0.27.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "9be70c98951c83b8d2f8f60d7065fa6d5146873094452a1008da8c2f1e4205ad"
+checksum = "c151a2a5ef800297b4e79efa4f4bec035c5f51d5ae587287c9b952bdf734cacd"
 dependencies = [
- "cfg-if",
- "js-sys",
- "libc",
- "wasi 0.10.2+wasi-snapshot-preview1",
- "wasm-bindgen",
+ "log",
+ "presser",
+ "thiserror 1.0.69",
+ "windows",
 ]
 
 [[package]]
-name = "getrandom"
-version = "0.3.3"
+name = "gpu-descriptor"
+version = "0.3.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "26145e563e54f2cadc477553f1ec5ee650b00862f0a58bcd12cbdc5f0ea2d2f4"
+checksum = "b89c83349105e3732062a895becfc71a8f921bb71ecbbdd8ff99263e3b53a0ca"
 dependencies = [
- "cfg-if",
- "libc",
- "r-efi",
- "wasi 0.14.2+wasi-0.2.4",
+ "bitflags 2.9.1",
+ "gpu-descriptor-types",
+ "hashbrown",
 ]
 
 [[package]]
-name = "glam"
-version = "0.29.3"
+name = "gpu-descriptor-types"
+version = "0.2.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "8babf46d4c1c9d92deac9f7be466f76dfc4482b6452fc5024b5e8daf6ffeb3ee"
+checksum = "fdf242682df893b86f33a73828fb09ca4b2d3bb6cc95249707fc684d27484b91"
 dependencies = [
- "bytemuck",
- "libm",
- "serde",
+ "bitflags 2.9.1",
 ]
 
 [[package]]
-name = "glob"
-version = "0.3.2"
+name = "guillotiere"
+version = "0.6.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "a8d1add55171497b4705a648c6b583acafb01d58050a51727785f0b2c8e0a2b2"
+checksum = "b62d5865c036cb1393e23c50693df631d3f5d7bcca4c04fe4cc0fd592e74a782"
+dependencies = [
+ "euclid",
+ "svg_fmt",
+]
 
 [[package]]
 name = "half"
@@ -926,21 +1806,14 @@ dependencies = [
 
 [[package]]
 name = "hashbrown"
-version = "0.12.3"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "8a9ee70c43aaf417c914396645a0fa852624801b24ebb7ae78fe8272889ac888"
-dependencies = [
- "ahash",
- "rayon",
-]
-
-[[package]]
-name = "hashbrown"
-version = "0.15.4"
+version = "0.15.5"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "5971ac85611da7067dbfcabef3c70ebb5606018acd9e2a3903a0da507521e0d5"
+checksum = "9229cfe53dfd69f0609a49f65461bd93001ea1ef889cd5529dd176593f5338a1"
 dependencies = [
+ "allocator-api2",
  "equivalent",
+ "foldhash",
+ "rayon",
  "serde",
 ]
 
@@ -956,31 +1829,47 @@ dependencies = [
 ]
 
 [[package]]
-name = "hermit-abi"
-version = "0.1.19"
+name = "heck"
+version = "0.5.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea"
+
+[[package]]
+name = "hexasphere"
+version = "15.1.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "62b467343b94ba476dcb2500d242dadbb39557df889310ac77c5d99100aaac33"
+checksum = "d9c9e718d32b6e6b2b32354e1b0367025efdd0b11d6a740b905ddf5db1074679"
 dependencies = [
- "libc",
+ "constgebra",
+ "glam 0.29.3",
+ "tinyvec",
 ]
 
 [[package]]
-name = "indexmap"
-version = "2.10.0"
+name = "hexf-parse"
+version = "0.2.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "fe4cd85333e22411419a0bcae1297d25e58c9443848b11dc6a86fefe8c78a661"
+checksum = "dfa686283ad6dd069f105e5ab091b04c62850d3e4cf5d67debad1933f55023df"
+
+[[package]]
+name = "image"
+version = "0.25.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "db35664ce6b9810857a38a906215e75a9c879f0696556a39f59c62829710251a"
 dependencies = [
- "equivalent",
- "hashbrown 0.15.4",
+ "bytemuck",
+ "byteorder-lite",
+ "num-traits",
 ]
 
 [[package]]
-name = "itertools"
-version = "0.10.3"
+name = "indexmap"
+version = "2.10.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "a9a9d19fa1e79b6215ff29b9d6880b706147f16e9b1dbb1e4e5947b5b02bc5e3"
+checksum = "fe4cd85333e22411419a0bcae1297d25e58c9443848b11dc6a86fefe8c78a661"
 dependencies = [
- "either",
+ "equivalent",
+ "hashbrown",
 ]
 
 [[package]]
@@ -1007,6 +1896,38 @@ version = "1.0.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "1aab8fc367588b89dcee83ab0fd66b72b50b72fa1904d7095045ace2b0c81c35"
 
+[[package]]
+name = "jni"
+version = "0.21.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1a87aa2bb7d2af34197c04845522473242e1aa17c12f4935d5856491a7fb8c97"
+dependencies = [
+ "cesu8",
+ "cfg-if",
+ "combine",
+ "jni-sys",
+ "log",
+ "thiserror 1.0.69",
+ "walkdir",
+ "windows-sys 0.45.0",
+]
+
+[[package]]
+name = "jni-sys"
+version = "0.3.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8eaf4bc02d17cbdd7ff4c7438cafcdf7fb9a4613313ad11b4f8fefe7d3fa0130"
+
+[[package]]
+name = "jobserver"
+version = "0.1.33"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "38f262f097c174adebe41eb73d66ae9c06b2844fb0da69969647bbddd9b0538a"
+dependencies = [
+ "getrandom 0.3.3",
+ "libc",
+]
+
 [[package]]
 name = "js-sys"
 version = "0.3.77"
@@ -1017,6 +1938,23 @@ dependencies = [
  "wasm-bindgen",
 ]
 
+[[package]]
+name = "khronos-egl"
+version = "6.0.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6aae1df220ece3c0ada96b8153459b67eebe9ae9212258bb0134ae60416fdf76"
+dependencies = [
+ "libc",
+ "libloading",
+ "pkg-config",
+]
+
+[[package]]
+name = "khronos_api"
+version = "3.1.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e2db585e1d738fc771bf08a151420d3ed193d9d895a36df7f6f8a9456b911ddc"
+
 [[package]]
 name = "lazy_static"
 version = "1.4.0"
@@ -1045,6 +1983,12 @@ version = "0.2.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "33a33a362ce288760ec6a508b94caaec573ae7d3bbbd91b87aa0bad4456839db"
 
+[[package]]
+name = "litrs"
+version = "0.4.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f5e54036fe321fd421e10d732f155734c4e4afd610dd556d9a82833ab3ee0bed"
+
 [[package]]
 name = "lock_api"
 version = "0.4.7"
@@ -1061,13 +2005,22 @@ version = "0.4.27"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "13dc2df351e3202783a1fe0d44375f7295ffb4049267b0f3018346dc122a1d94"
 
+[[package]]
+name = "malloc_buf"
+version = "0.0.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "62bb907fe88d54d8d9ce32a3cceab4218ed2f6b7d35617cafe9adf84e43919cb"
+dependencies = [
+ "libc",
+]
+
 [[package]]
 name = "matchers"
 version = "0.1.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "8263075bb86c5a1b1427b5ae862e8889656f126e9f77c484496e8b47cf5c5558"
 dependencies = [
- "regex-automata",
+ "regex-automata 0.1.10",
 ]
 
 [[package]]
@@ -1094,6 +2047,21 @@ dependencies = [
  "autocfg",
 ]
 
+[[package]]
+name = "metal"
+version = "0.31.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f569fb946490b5743ad69813cb19629130ce9374034abe31614a36402d18f99e"
+dependencies = [
+ "bitflags 2.9.1",
+ "block",
+ "core-graphics-types",
+ "foreign-types",
+ "log",
+ "objc",
+ "paste",
+]
+
 [[package]]
 name = "minimal-lexical"
 version = "0.2.1"
@@ -1109,13 +2077,72 @@ dependencies = [
  "serde",
 ]
 
+[[package]]
+name = "naga"
+version = "24.0.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e380993072e52eef724eddfcde0ed013b0c023c3f0417336ed041aa9f076994e"
+dependencies = [
+ "arrayvec",
+ "bit-set 0.8.0",
+ "bitflags 2.9.1",
+ "cfg_aliases",
+ "codespan-reporting",
+ "hexf-parse",
+ "indexmap",
+ "log",
+ "pp-rs",
+ "rustc-hash",
+ "spirv",
+ "strum",
+ "termcolor",
+ "thiserror 2.0.12",
+ "unicode-xid",
+]
+
+[[package]]
+name = "naga_oil"
+version = "0.17.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f2464f7395decfd16bb4c33fb0cb3b2c645cc60d051bc7fb652d3720bfb20f18"
+dependencies = [
+ "bit-set 0.5.3",
+ "codespan-reporting",
+ "data-encoding",
+ "indexmap",
+ "naga",
+ "once_cell",
+ "regex",
+ "regex-syntax 0.8.5",
+ "rustc-hash",
+ "thiserror 1.0.69",
+ "tracing",
+ "unicode-ident",
+]
+
 [[package]]
 name = "nalgebra"
-version = "0.33.2"
+version = "0.34.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "26aecdf64b707efd1310e3544d709c5c0ac61c13756046aaaba41be5c4f66a3b"
+checksum = "9cd59afb6639828b33677758314a4a1a745c15c02bc597095b851c8fd915cf49"
 dependencies = [
  "approx",
+ "glam 0.14.0",
+ "glam 0.15.2",
+ "glam 0.16.0",
+ "glam 0.17.3",
+ "glam 0.18.0",
+ "glam 0.19.0",
+ "glam 0.20.5",
+ "glam 0.21.3",
+ "glam 0.22.0",
+ "glam 0.23.0",
+ "glam 0.24.2",
+ "glam 0.25.0",
+ "glam 0.27.0",
+ "glam 0.28.0",
+ "glam 0.29.3",
+ "glam 0.30.5",
  "matrixmultiply",
  "nalgebra-macros",
  "num-complex",
@@ -1128,15 +2155,53 @@ dependencies = [
 
 [[package]]
 name = "nalgebra-macros"
-version = "0.2.2"
+version = "0.3.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "254a5372af8fc138e36684761d3c0cdb758a4410e938babcff1c860ce14ddbfc"
+checksum = "973e7178a678cfd059ccec50887658d482ce16b0aa9da3888ddeab5cd5eb4889"
 dependencies = [
  "proc-macro2",
  "quote",
  "syn 2.0.104",
 ]
 
+[[package]]
+name = "ndk"
+version = "0.9.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c3f42e7bbe13d351b6bead8286a43aac9534b82bd3cc43e47037f012ebfd62d4"
+dependencies = [
+ "bitflags 2.9.1",
+ "jni-sys",
+ "log",
+ "ndk-sys 0.6.0+11769913",
+ "num_enum",
+ "thiserror 1.0.69",
+]
+
+[[package]]
+name = "ndk-context"
+version = "0.1.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "27b02d87554356db9e9a873add8782d4ea6e3e58ea071a9adb9a2e8ddb884a8b"
+
+[[package]]
+name = "ndk-sys"
+version = "0.5.0+25.2.9519653"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8c196769dd60fd4f363e11d948139556a344e79d451aeb2fa2fd040738ef7691"
+dependencies = [
+ "jni-sys",
+]
+
+[[package]]
+name = "ndk-sys"
+version = "0.6.0+11769913"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ee6cda3051665f1fb8d9e08fc35c96d5a244fb1be711a03b71118828afc9a873"
+dependencies = [
+ "jni-sys",
+]
+
 [[package]]
 name = "nix"
 version = "0.30.1"
@@ -1207,13 +2272,44 @@ dependencies = [
 ]
 
 [[package]]
-name = "num_cpus"
-version = "1.13.1"
+name = "num_enum"
+version = "0.7.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "19e64526ebdee182341572e50e9ad03965aa510cd94427a4549448f285e957a1"
+checksum = "a973b4e44ce6cad84ce69d797acf9a044532e4184c4f267913d1b546a0727b7a"
 dependencies = [
- "hermit-abi",
- "libc",
+ "num_enum_derive",
+ "rustversion",
+]
+
+[[package]]
+name = "num_enum_derive"
+version = "0.7.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "77e878c846a8abae00dd069496dbe8751b16ac1c3d6bd2a7283a938e8228f90d"
+dependencies = [
+ "proc-macro-crate",
+ "proc-macro2",
+ "quote",
+ "syn 2.0.104",
+]
+
+[[package]]
+name = "objc"
+version = "0.2.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "915b1b472bc21c53464d6c8461c9d3af805ba1ef837e1cac254428f4a77177b1"
+dependencies = [
+ "malloc_buf",
+]
+
+[[package]]
+name = "offset-allocator"
+version = "0.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e234d535da3521eb95106f40f0b73483d80bfb3aacf27c40d7e2b72f1a3e00a2"
+dependencies = [
+ "log",
+ "nonmax",
 ]
 
 [[package]]
@@ -1228,6 +2324,15 @@ version = "11.1.3"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "0ab1bc2a289d34bd04a330323ac98a1b4bc82c9d9fcb1e66b63caa84da26b575"
 
+[[package]]
+name = "ordered-float"
+version = "4.6.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7bb71e1b3fa6ca1c61f383464aaf2bb0e2f8e772a1f01d486832464de363b951"
+dependencies = [
+ "num-traits",
+]
+
 [[package]]
 name = "parking"
 version = "2.2.1"
@@ -1285,9 +2390,26 @@ dependencies = [
 
 [[package]]
 name = "pin-project-lite"
-version = "0.2.9"
+version = "0.2.16"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3b3cff922bd51709b605d9ead9aa71031d81447142d828eb4a6eba76fe619f9b"
+
+[[package]]
+name = "piper"
+version = "0.2.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "96c8c490f422ef9a4efd2cb5b42b76c8613d7e7dfc1caf667b8a3350a5acc066"
+dependencies = [
+ "atomic-waker",
+ "fastrand",
+ "futures-io",
+]
+
+[[package]]
+name = "pkg-config"
+version = "0.3.32"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "e0a7ae3ac2f1173085d398531c705756c94a4c56843785df85a60c1a0afac116"
+checksum = "7edddbd0b52d732b21ad9a5fab5c704c14cd949e5e9a1ec5929a24fded1b904c"
 
 [[package]]
 name = "plotters"
@@ -1332,12 +2454,27 @@ dependencies = [
  "portable-atomic",
 ]
 
+[[package]]
+name = "pp-rs"
+version = "0.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "bb458bb7f6e250e6eb79d5026badc10a3ebb8f9a15d1fff0f13d17c71f4d6dee"
+dependencies = [
+ "unicode-xid",
+]
+
 [[package]]
 name = "ppv-lite86"
 version = "0.2.16"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "eb9f9e6e233e5c4a35559a617bf40a4ec447db2e84c20b55a6f83167b7e57872"
 
+[[package]]
+name = "presser"
+version = "0.3.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e8cf8e6a8aa66ce33f63993ffc4ea4271eb5b0530a9002db8455ea6050c77bfa"
+
 [[package]]
 name = "prettyplease"
 version = "0.2.35"
@@ -1348,6 +2485,15 @@ dependencies = [
  "syn 2.0.104",
 ]
 
+[[package]]
+name = "proc-macro-crate"
+version = "3.3.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "edce586971a4dfaa28950c6f18ed55e0406c1ab88bbce2c6f6293a7aaba73d35"
+dependencies = [
+ "toml_edit",
+]
+
 [[package]]
 name = "proc-macro2"
 version = "1.0.95"
@@ -1357,6 +2503,12 @@ dependencies = [
  "unicode-ident",
 ]
 
+[[package]]
+name = "profiling"
+version = "1.0.17"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3eb8486b569e12e2c32ad3e204dbaba5e4b5b216e9367044f25f1dba42341773"
+
 [[package]]
 name = "quote"
 version = "1.0.40"
@@ -1385,9 +2537,9 @@ dependencies = [
 
 [[package]]
 name = "rand"
-version = "0.9.1"
+version = "0.9.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "9fbfd9d094a40bf3ae768db9361049ace4c0e04a4fd6b359518bd7b73a73dd97"
+checksum = "6db2770f06117d490610c7488547d543617b21bfa07796d7a12f6f1bd53850d1"
 dependencies = [
  "rand_chacha 0.9.0",
  "rand_core 0.9.3",
@@ -1448,9 +2600,21 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "6a8615d50dcf34fa31f7ab52692afec947c4dd0ab803cc87cb3b0b4570ff7463"
 dependencies = [
  "num-traits",
- "rand 0.9.1",
+ "rand 0.9.2",
 ]
 
+[[package]]
+name = "range-alloc"
+version = "0.1.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c3d6831663a5098ea164f89cff59c6284e95f4e3c76ce9848d4529f5ccca9bde"
+
+[[package]]
+name = "raw-window-handle"
+version = "0.6.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "20675572f6f24e9e76ef639bc5552774ed45f1c30e2951e1e99c59888861c539"
+
 [[package]]
 name = "rawpointer"
 version = "0.2.1"
@@ -1459,28 +2623,30 @@ checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3"
 
 [[package]]
 name = "rayon"
-version = "1.5.2"
+version = "1.10.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "fd249e82c21598a9a426a4e00dd7adc1d640b22445ec8545feef801d1a74c221"
+checksum = "b418a60154510ca1a002a752ca9714984e21e4241e804d32555251faf8b78ffa"
 dependencies = [
- "autocfg",
- "crossbeam-deque",
  "either",
  "rayon-core",
 ]
 
 [[package]]
 name = "rayon-core"
-version = "1.9.2"
+version = "1.12.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "9f51245e1e62e1f1629cbfec37b5793bbabcaeb90f30e94d2ba03564687353e4"
+checksum = "1465873a3dfdaa8ae7cb14b4383657caab0b3e8a0aa9ae8e04b044854c8dfce2"
 dependencies = [
- "crossbeam-channel",
  "crossbeam-deque",
  "crossbeam-utils",
- "num_cpus",
 ]
 
+[[package]]
+name = "rectangle-pack"
+version = "0.4.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a0d463f2884048e7153449a55166f91028d5b0ea53c79377099ce4e8cf0cf9bb"
+
 [[package]]
 name = "redox_syscall"
 version = "0.2.13"
@@ -1492,11 +2658,14 @@ dependencies = [
 
 [[package]]
 name = "regex"
-version = "1.5.5"
+version = "1.11.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "1a11647b6b25ff05a515cb92c365cec08801e83423a235b51e231e1808747286"
+checksum = "b544ef1b4eac5dc2db33ea63606ae9ffcfac26c1416a2806ae0bf5f56b201191"
 dependencies = [
- "regex-syntax",
+ "aho-corasick",
+ "memchr",
+ "regex-automata 0.4.9",
+ "regex-syntax 0.8.5",
 ]
 
 [[package]]
@@ -1505,7 +2674,18 @@ version = "0.1.10"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "6c230d73fb8d8c1b9c0b3135c5142a8acee3a0558fb8db5cf1cb65f8d7862132"
 dependencies = [
- "regex-syntax",
+ "regex-syntax 0.6.25",
+]
+
+[[package]]
+name = "regex-automata"
+version = "0.4.9"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "809e8dc61f6de73b46c85f4c96486310fe304c434cfa43669d7b40f711150908"
+dependencies = [
+ "aho-corasick",
+ "memchr",
+ "regex-syntax 0.8.5",
 ]
 
 [[package]]
@@ -1514,6 +2694,30 @@ version = "0.6.25"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "f497285884f3fcff424ffc933e56d7cbca511def0c9831a7f9b5f6153e3cc89b"
 
+[[package]]
+name = "regex-syntax"
+version = "0.8.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2b15c43186be67a4fd63bee50d0303afffcef381492ebe2c5d87f324e1b8815c"
+
+[[package]]
+name = "renderdoc-sys"
+version = "1.1.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "19b30a45b0cd0bcca8037f3d0dc3421eaf95327a17cad11964fb8179b4fc4832"
+
+[[package]]
+name = "ron"
+version = "0.8.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b91f7eff05f748767f183df4320a63d6936e9c6107d97c9e6bdd9784f4289c94"
+dependencies = [
+ "base64",
+ "bitflags 2.9.1",
+ "serde",
+ "serde_derive",
+]
+
 [[package]]
 name = "rustc-hash"
 version = "1.1.0"
@@ -1556,6 +2760,12 @@ version = "1.1.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "d29ab0c6d3fc0ee92fe66e2d99f700eab17a8d57d1c1d3b748380fb20baa78cd"
 
+[[package]]
+name = "send_wrapper"
+version = "0.6.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "cd0b0ec5f1c1ca621c432a25813d8d60c88abe6d3e08a3eb9cf37d97a0fe3d73"
+
 [[package]]
 name = "serde"
 version = "1.0.219"
@@ -1587,9 +2797,9 @@ dependencies = [
 
 [[package]]
 name = "serde_json"
-version = "1.0.140"
+version = "1.0.142"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "20068b6e96dc6c9bd23e01df8827e6c7e1f2fddd43c21810382803c136b99373"
+checksum = "030fedb782600dcbd6f02d479bf0d817ac3bb40d644745b769d6a96bc3afc5a7"
 dependencies = [
  "itoa",
  "memchr",
@@ -1644,6 +2854,15 @@ version = "0.4.10"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "04dc19736151f35336d325007ac991178d504a119863a2fcb3758cdb5e52c50d"
 
+[[package]]
+name = "slotmap"
+version = "1.0.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "dbff4acf519f630b3a3ddcfaea6c06b42174d9a44bc70c620e9ed1649d58b82a"
+dependencies = [
+ "version_check",
+]
+
 [[package]]
 name = "smallvec"
 version = "1.15.1"
@@ -1668,12 +2887,61 @@ dependencies = [
  "portable-atomic",
 ]
 
+[[package]]
+name = "spirv"
+version = "0.3.0+sdk-1.3.268.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "eda41003dc44290527a59b13432d4a0379379fa074b70174882adfbdfd917844"
+dependencies = [
+ "bitflags 2.9.1",
+]
+
 [[package]]
 name = "stable_deref_trait"
 version = "1.2.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "a8f112729512f8e442d81f95a8a7ddf2b7c6b8a1a6f509a95864142b30cab2d3"
 
+[[package]]
+name = "stackfuture"
+version = "0.3.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6eae92052b72ef70dafa16eddbabffc77e5ca3574be2f7bc1127b36f0a7ad7f2"
+
+[[package]]
+name = "static_assertions"
+version = "1.1.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a2eb9349b6444b326872e140eb1cf5e7c522154d69e7a0ffb0fb81c06b37543f"
+
+[[package]]
+name = "strum"
+version = "0.26.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8fec0f0aef304996cf250b31b5a10dee7980c85da9d759361292b8bca5a18f06"
+dependencies = [
+ "strum_macros",
+]
+
+[[package]]
+name = "strum_macros"
+version = "0.26.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4c6bee85a5a24955dc440386795aa378cd9cf82acd5f764469152d2270e581be"
+dependencies = [
+ "heck",
+ "proc-macro2",
+ "quote",
+ "rustversion",
+ "syn 2.0.104",
+]
+
+[[package]]
+name = "svg_fmt"
+version = "0.4.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0193cc4331cfd2f3d2011ef287590868599a2f33c3e69bc22c1a3d3acf9e02fb"
+
 [[package]]
 name = "syn"
 version = "1.0.91"
@@ -1696,13 +2964,42 @@ dependencies = [
  "unicode-ident",
 ]
 
+[[package]]
+name = "termcolor"
+version = "1.4.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "06794f8f6c5c898b3275aebefa6b8a1cb24cd2c6c79397ab15774837a0bc5755"
+dependencies = [
+ "winapi-util",
+]
+
+[[package]]
+name = "thiserror"
+version = "1.0.69"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b6aaf5339b578ea85b50e080feb250a3e8ae8cfcdff9a461c9ec2904bc923f52"
+dependencies = [
+ "thiserror-impl 1.0.69",
+]
+
 [[package]]
 name = "thiserror"
 version = "2.0.12"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "567b8a2dae586314f7be2a752ec7474332959c6460e02bde30d702a66d488708"
 dependencies = [
- "thiserror-impl",
+ "thiserror-impl 2.0.12",
+]
+
+[[package]]
+name = "thiserror-impl"
+version = "1.0.69"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4fee6c4efc90059e10f81e6d42c60a18f76588c3d74cb83a0b242a2b6c7504c1"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.104",
 ]
 
 [[package]]
@@ -1735,6 +3032,21 @@ dependencies = [
  "serde_json",
 ]
 
+[[package]]
+name = "tinyvec"
+version = "1.9.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "09b3661f17e86524eccd4371ab0429194e0d7c008abb45f7a7495b1719463c71"
+dependencies = [
+ "tinyvec_macros",
+]
+
+[[package]]
+name = "tinyvec_macros"
+version = "0.1.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1f3ccbac311fea05f86f61904b462b55fb3df8837a366dfc601a0161d0532f20"
+
 [[package]]
 name = "toml_datetime"
 version = "0.6.11"
@@ -1869,11 +3181,17 @@ version = "1.0.18"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "5a5f39404a5da50712a4c1eecf25e90dd62b613502b7e925fd4e4d19b5c96512"
 
+[[package]]
+name = "unicode-width"
+version = "0.1.14"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7dd6e30e90baa6f72411720665d41d89b9a3d039dc45b8faea1ddd07f617f6af"
+
 [[package]]
 name = "unicode-xid"
-version = "0.2.3"
+version = "0.2.6"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "957e51f3646910546462e67d5f7599b9e4fb8acdd304b087a6494730f9eebf04"
+checksum = "ebc1c04c71510c7f702b52b7c350734c9ff1295c464a03335b00bb84fc54f853"
 
 [[package]]
 name = "unsafe-libyaml"
@@ -1970,12 +3288,13 @@ dependencies = [
 
 [[package]]
 name = "wasm-bindgen-futures"
-version = "0.4.34"
+version = "0.4.50"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "f219e0d211ba40266969f6dbdd90636da12f75bee4fc9d6c23d1260dadb51454"
+checksum = "555d470ec0bc3bb57890405e5d4322cc9ea83cebb085523ced7be4144dac1e61"
 dependencies = [
  "cfg-if",
  "js-sys",
+ "once_cell",
  "wasm-bindgen",
  "web-sys",
 ]
@@ -2032,6 +3351,103 @@ dependencies = [
  "wasm-bindgen",
 ]
 
+[[package]]
+name = "wgpu"
+version = "24.0.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6b0b3436f0729f6cdf2e6e9201f3d39dc95813fad61d826c1ed07918b4539353"
+dependencies = [
+ "arrayvec",
+ "bitflags 2.9.1",
+ "cfg_aliases",
+ "document-features",
+ "js-sys",
+ "log",
+ "naga",
+ "parking_lot",
+ "profiling",
+ "raw-window-handle",
+ "smallvec",
+ "static_assertions",
+ "wasm-bindgen",
+ "wasm-bindgen-futures",
+ "web-sys",
+ "wgpu-core",
+ "wgpu-hal",
+ "wgpu-types",
+]
+
+[[package]]
+name = "wgpu-core"
+version = "24.0.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7f0aa306497a238d169b9dc70659105b4a096859a34894544ca81719242e1499"
+dependencies = [
+ "arrayvec",
+ "bit-vec 0.8.0",
+ "bitflags 2.9.1",
+ "cfg_aliases",
+ "document-features",
+ "indexmap",
+ "log",
+ "naga",
+ "once_cell",
+ "parking_lot",
+ "profiling",
+ "raw-window-handle",
+ "rustc-hash",
+ "smallvec",
+ "thiserror 2.0.12",
+ "wgpu-hal",
+ "wgpu-types",
+]
+
+[[package]]
+name = "wgpu-hal"
+version = "24.0.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f112f464674ca69f3533248508ee30cb84c67cf06c25ff6800685f5e0294e259"
+dependencies = [
+ "android_system_properties",
+ "arrayvec",
+ "ash",
+ "bit-set 0.8.0",
+ "bitflags 2.9.1",
+ "block",
+ "bytemuck",
+ "cfg_aliases",
+ "core-graphics-types",
+ "glow",
+ "glutin_wgl_sys",
+ "gpu-alloc",
+ "gpu-allocator",
+ "gpu-descriptor",
+ "js-sys",
+ "khronos-egl",
+ "libc",
+ "libloading",
+ "log",
+ "metal",
+ "naga",
+ "ndk-sys 0.5.0+25.2.9519653",
+ "objc",
+ "once_cell",
+ "ordered-float",
+ "parking_lot",
+ "profiling",
+ "range-alloc",
+ "raw-window-handle",
+ "renderdoc-sys",
+ "rustc-hash",
+ "smallvec",
+ "thiserror 2.0.12",
+ "wasm-bindgen",
+ "web-sys",
+ "wgpu-types",
+ "windows",
+ "windows-core",
+]
+
 [[package]]
 name = "wgpu-types"
 version = "24.0.0"
@@ -2086,6 +3502,70 @@ version = "0.4.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f"
 
+[[package]]
+name = "windows"
+version = "0.58.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "dd04d41d93c4992d421894c18c8b43496aa748dd4c081bac0dc93eb0489272b6"
+dependencies = [
+ "windows-core",
+ "windows-targets 0.52.6",
+]
+
+[[package]]
+name = "windows-core"
+version = "0.58.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6ba6d44ec8c2591c134257ce647b7ea6b20335bf6379a27dac5f1641fcf59f99"
+dependencies = [
+ "windows-implement",
+ "windows-interface",
+ "windows-result",
+ "windows-strings",
+ "windows-targets 0.52.6",
+]
+
+[[package]]
+name = "windows-implement"
+version = "0.58.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2bbd5b46c938e506ecbce286b6628a02171d56153ba733b6c741fc627ec9579b"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.104",
+]
+
+[[package]]
+name = "windows-interface"
+version = "0.58.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "053c4c462dc91d3b1504c6fe5a726dd15e216ba718e84a0e46a88fbe5ded3515"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.104",
+]
+
+[[package]]
+name = "windows-result"
+version = "0.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1d1043d8214f791817bab27572aaa8af63732e11bf84aa21a45a78d6c317ae0e"
+dependencies = [
+ "windows-targets 0.52.6",
+]
+
+[[package]]
+name = "windows-strings"
+version = "0.1.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4cd9b125c486025df0eabcb585e62173c6c9eddcec5d117d3b6e8c30e2ee4d10"
+dependencies = [
+ "windows-result",
+ "windows-targets 0.52.6",
+]
+
 [[package]]
 name = "windows-sys"
 version = "0.36.1"
@@ -2099,6 +3579,15 @@ dependencies = [
  "windows_x86_64_msvc 0.36.1",
 ]
 
+[[package]]
+name = "windows-sys"
+version = "0.45.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "75283be5efb2831d37ea142365f009c02ec203cd29a3ebecbc093d52315b66d0"
+dependencies = [
+ "windows-targets 0.42.2",
+]
+
 [[package]]
 name = "windows-sys"
 version = "0.59.0"
@@ -2108,6 +3597,21 @@ dependencies = [
  "windows-targets 0.52.6",
 ]
 
+[[package]]
+name = "windows-targets"
+version = "0.42.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8e5180c00cd44c9b1c88adb3693291f1cd93605ded80c250a75d472756b4d071"
+dependencies = [
+ "windows_aarch64_gnullvm 0.42.2",
+ "windows_aarch64_msvc 0.42.2",
+ "windows_i686_gnu 0.42.2",
+ "windows_i686_msvc 0.42.2",
+ "windows_x86_64_gnu 0.42.2",
+ "windows_x86_64_gnullvm 0.42.2",
+ "windows_x86_64_msvc 0.42.2",
+]
+
 [[package]]
 name = "windows-targets"
 version = "0.52.6"
@@ -2140,6 +3644,12 @@ dependencies = [
  "windows_x86_64_msvc 0.53.0",
 ]
 
+[[package]]
+name = "windows_aarch64_gnullvm"
+version = "0.42.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "597a5118570b68bc08d8d59125332c54f1ba9d9adeedeef5b99b02ba2b0698f8"
+
 [[package]]
 name = "windows_aarch64_gnullvm"
 version = "0.52.6"
@@ -2158,6 +3668,12 @@ version = "0.36.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "9bb8c3fd39ade2d67e9874ac4f3db21f0d710bee00fe7cab16949ec184eeaa47"
 
+[[package]]
+name = "windows_aarch64_msvc"
+version = "0.42.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e08e8864a60f06ef0d0ff4ba04124db8b0fb3be5776a5cd47641e942e58c4d43"
+
 [[package]]
 name = "windows_aarch64_msvc"
 version = "0.52.6"
@@ -2176,6 +3692,12 @@ version = "0.36.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "180e6ccf01daf4c426b846dfc66db1fc518f074baa793aa7d9b9aaeffad6a3b6"
 
+[[package]]
+name = "windows_i686_gnu"
+version = "0.42.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c61d927d8da41da96a81f029489353e68739737d3beca43145c8afec9a31a84f"
+
 [[package]]
 name = "windows_i686_gnu"
 version = "0.52.6"
@@ -2206,6 +3728,12 @@ version = "0.36.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "e2e7917148b2812d1eeafaeb22a97e4813dfa60a3f8f78ebe204bcc88f12f024"
 
+[[package]]
+name = "windows_i686_msvc"
+version = "0.42.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "44d840b6ec649f480a41c8d80f9c65108b92d89345dd94027bfe06ac444d1060"
+
 [[package]]
 name = "windows_i686_msvc"
 version = "0.52.6"
@@ -2224,6 +3752,12 @@ version = "0.36.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "4dcd171b8776c41b97521e5da127a2d86ad280114807d0b2ab1e462bc764d9e1"
 
+[[package]]
+name = "windows_x86_64_gnu"
+version = "0.42.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8de912b8b8feb55c064867cf047dda097f92d51efad5b491dfb98f6bbb70cb36"
+
 [[package]]
 name = "windows_x86_64_gnu"
 version = "0.52.6"
@@ -2236,6 +3770,12 @@ version = "0.53.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "2e55b5ac9ea33f2fc1716d1742db15574fd6fc8dadc51caab1c16a3d3b4190ba"
 
+[[package]]
+name = "windows_x86_64_gnullvm"
+version = "0.42.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "26d41b46a36d453748aedef1486d5c7a85db22e56aff34643984ea85514e94a3"
+
 [[package]]
 name = "windows_x86_64_gnullvm"
 version = "0.52.6"
@@ -2254,6 +3794,12 @@ version = "0.36.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "c811ca4a8c853ef420abd8592ba53ddbbac90410fab6903b3e79972a631f7680"
 
+[[package]]
+name = "windows_x86_64_msvc"
+version = "0.42.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9aec5da331524158c6d1a4ac0ab1541149c0b9505fde06423b02f5ef0106b9f0"
+
 [[package]]
 name = "windows_x86_64_msvc"
 version = "0.52.6"
@@ -2283,3 +3829,9 @@ checksum = "6f42320e61fe2cfd34354ecb597f86f413484a798ba44a8ca1165c58d42da6c1"
 dependencies = [
  "bitflags 2.9.1",
 ]
+
+[[package]]
+name = "xml-rs"
+version = "0.8.27"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6fd8403733700263c6eb89f192880191f1b83e332f7a20371ddcf421c4a337c7"
diff --git a/Cargo.toml b/Cargo.toml
index f513324c..8f53439d 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -16,23 +16,34 @@ resolver = "2"
 
 [dependencies]
 atomecs-derive = "0.1.0"
-bevy = { version = "0.16.1", default-features = false, features = ["bevy_log"] }
-nalgebra = { version = "0.33.2", features = ["serde-serialize"] }
-rand = "0.9.1"
+bevy = { version = "0.16.1", default-features = false, features = [
+    "bevy_log",
+    "multi_threaded",
+] }
+bytemuck = "1.9.1" # I added this with the feature multi-threaded, I have no idea what it does.
+nalgebra = { version = "0.34.0", features = ["serde-serialize"] }
+rand = "0.9.2"
 rand_distr = "0.5.1"
 serde = { version = "1.0", features = ["derive"] }
-serde_json = "1.0.140"
+serde_json = "1.0.142"
 serde_yaml = "0.9.34-deprecated"
 assert_approx_eq = "1.1.0"
+approx = "0.5.1"
 csv = "1.3.1"
 byteorder = "1.5.0"
 multimap = "0.10.1"
 serde_arrays = "0.2.0"
-hashbrown = { version = "^0.12.1", features = ["rayon"] }
-rayon = "1.5.0"
+hashbrown = { version = "^0.15.5", features = ["rayon"] }
+rayon = "1.10.0"
+regex = { version = "1.10.4", default-features = false, features = [
+    "std",
+    "unicode-perl",
+    "unicode-case",
+] } # This was added to fix panics with tracing subscriber, 
+# after adding multi-threaded feature to bevy. Is unnecesary if you dont attempt to add TaskPoolPlugin, which seems to be just as fast.
 
 [dev-dependencies]
-criterion = "0.6"
+criterion = "0.7.0"
 
 [profile.release]
 opt-level = 3
diff --git a/examples/benchmark.rs b/examples/benchmark.rs
index 29a65c88..71abf733 100644
--- a/examples/benchmark.rs
+++ b/examples/benchmark.rs
@@ -20,6 +20,8 @@ use std::fs::read_to_string;
 use std::fs::File;
 use std::time::Instant;
 
+use bevy::app::TaskPoolThreadAssignmentPolicy;
+
 extern crate serde;
 use serde::{Deserialize, Serialize};
 
@@ -32,7 +34,7 @@ pub struct BenchmarkConfiguration {
 impl Default for BenchmarkConfiguration {
     fn default() -> Self {
         BenchmarkConfiguration {
-            n_atoms: 10000,
+            n_atoms: 1_000_000,
             n_threads: 12,
             n_steps: 5000,
         }
@@ -56,14 +58,6 @@ fn main() {
     let mut app = App::new();
     app.add_plugins(LaserPlugin::<{ BEAM_NUMBER }>);
     app.add_plugins(LaserCoolingPlugin::<Rubidium87_780D2, { BEAM_NUMBER }>::default());
-    app.add_plugins(
-        DefaultPlugins.set(TaskPoolPlugin {
-            task_pool_options: TaskPoolOptions::with_num_threads(10),
-        }), // .set(LogPlugin {
-            //     level: Level::DEBUG,
-            //     filter: "bevy_core=trace".to_string(),
-            // }),
-    );
     app.add_plugins(atomecs::integrator::IntegrationPlugin);
     app.add_plugins(atomecs::initiate::InitiatePlugin);
     app.add_plugins(atomecs::magnetic::MagneticsPlugin);
@@ -72,6 +66,60 @@ fn main() {
     //app.add_startup_system(setup_world);
 
     // TODO: Configure bevy compute pool size
+    
+    // Enabling the task pool plugin with custom thread assignment policy settings,
+    // requires the regex features to be enabled in Cargo.toml. This seems to perform the same as the default settings.
+    let task_pool_options = TaskPoolOptions {
+        // Use 25% of cores for IO, at least 1, no more than 4
+        io: TaskPoolThreadAssignmentPolicy {
+            min_threads: 0,
+            max_threads: 0,
+            percent: 0.0,
+            on_thread_spawn: None,
+            on_thread_destroy: None,
+        },
+
+        // Use 25% of cores for async compute, at least 1, no more than 4
+        async_compute: TaskPoolThreadAssignmentPolicy {
+            min_threads: 0,
+            max_threads: 0,
+            percent: 0.0,
+            on_thread_spawn: None,
+            on_thread_destroy: None,
+        },
+        min_total_threads: 1,
+        max_total_threads: usize::MAX,
+        compute: TaskPoolThreadAssignmentPolicy {
+            min_threads: 1,
+            max_threads: usize::MAX,
+            percent: 100.0,
+            on_thread_spawn: None,
+            on_thread_destroy: None,
+        },
+    };
+    
+    app.add_plugins(
+        DefaultPlugins
+        .set(TaskPoolPlugin {
+            task_pool_options
+        }),
+        // .set(LogPlugin {
+        //     level: Level::DEBUG,
+        //     filter: "bevy_core=trace".to_string(),
+        // }),
+    );
+
+    // This (original) however seems to be worse
+    // app.add_plugins(
+    //     DefaultPlugins
+    //     .set(TaskPoolPlugin {
+    //         task_pool_options: TaskPoolOptions::with_num_threads(10),
+    //     }),
+    //     // .set(LogPlugin {
+    //     //     level: Level::DEBUG,
+    //     //     filter: "bevy_core=trace".to_string(),
+    //     // }),
+    // );
 
     // Create magnetic field.
     app.world_mut()
diff --git a/examples/jtheta_test.rs b/examples/jtheta_test.rs
new file mode 100644
index 00000000..407d5345
--- /dev/null
+++ b/examples/jtheta_test.rs
@@ -0,0 +1,91 @@
+extern crate atomecs as lib;
+extern crate nalgebra;
+use nalgebra::Vector3;
+use bevy::prelude::*;
+use std::time::Instant;
+use std::marker::PhantomData;
+use lib::constant::PI;
+use lib::atom::{Atom, Position, Velocity};
+use lib::integrator::Timestep;
+use lib::output::file::FileOutputPlugin;
+use lib::output::file::Text;
+use lib::simulation::SimulationBuilder;
+use lib::sim_region::{SimulationVolume, VolumeType};
+use lib::shapes::{Cylinder as MyCylinder, CylindricalPipe};
+use lib::species::Strontium88;
+use lib::atom_sources::{AtomSourcePlugin, WeightedProbabilityDistribution, VelocityCap};
+use lib::atom_sources::emit::{EmitFixedRate, AtomNumberToEmit};
+use lib::atom_sources::mass::{MassRatio, MassDistribution};
+use lib::atom_sources::oven::{OvenAperture, Oven};
+use lib::collisions::{CollisionPlugin, ApplyAtomCollisions, ApplyWallCollisions};
+use lib::collisions::wall_collisions::{WallData, WallType};
+
+fn main() {
+    let now = Instant::now();
+    let mut sim_builder = SimulationBuilder::default();
+    sim_builder.add_plugins(AtomSourcePlugin::<Strontium88>::default());
+    sim_builder.add_plugins(FileOutputPlugin::<Position, Text, Atom>::new("pos.txt".to_string(),10));
+    sim_builder.add_plugins(FileOutputPlugin::<Velocity, Text, Atom>::new("vel.txt".to_string(),10));
+    // sim_builder.add_plugins(FileOutputPlugin::<NumberOfWallCollisions, Text, Atom>::new("wall_collisions.txt".to_string(),1));
+    sim_builder.add_plugins(CollisionPlugin);
+
+    let mut sim = sim_builder.build();
+
+    sim.world_mut().insert_resource(ApplyAtomCollisions(false));
+    sim.world_mut().insert_resource(ApplyWallCollisions(true));
+
+    sim.world_mut()
+        .spawn(WallData{
+            wall_type: WallType::Rough})
+        .insert(CylindricalPipe::new(1e-5, 400e-6, Vector3::new(1.0, 0.0, 0.0)))
+        .insert(Position{pos: Vector3::new(-200e-6, 0.0, 0.0)});
+
+    sim.world_mut()
+        .spawn(SimulationVolume{volume_type: VolumeType::Inclusive})
+        .insert(MyCylinder::new(110e-5, 2000e-6, Vector3::new(1.0, 0.0, 0.0)))
+        .insert(Position{pos: Vector3::new(600e-6, 0.0, 0.0)});
+
+    let number_to_emit = 1e10;
+
+    let mut thetas = Vec::<f64>::new();
+    let mut weights = Vec::<f64>::new();
+
+    let n = 10000;
+    for i in 0..n {
+        let theta = (i as f64) / (n as f64 + 1.0) * PI / 2.0;
+        let weight = theta.sin() * theta.cos();
+        thetas.push(theta);
+        weights.push(weight);
+    }
+
+    let uniform_distribution = WeightedProbabilityDistribution::new(thetas, weights);
+
+    sim.world_mut()
+        .spawn(Oven::<Strontium88> {
+            temperature: 700.0,
+            aperture: OvenAperture::Circular{radius: 1e-5, thickness: 1e-9},
+            direction: Vector3::new(1.0, 0.0, 0.0),
+            theta_distribution: uniform_distribution,
+            max_theta: PI/2.0,
+            phantom: PhantomData,
+        })
+        .insert(Position {
+            pos: Vector3::new(-39999e-8, 0.0, 0.0),
+        })
+        .insert(MassDistribution::new(vec![MassRatio {
+            mass: 88.0,
+            ratio: 1.0,
+        }]))
+        .insert(AtomNumberToEmit{number: 0})
+        .insert(EmitFixedRate{rate: number_to_emit});
+
+    // Define timestep
+    sim.world_mut().insert_resource(Timestep { delta: 1e-7 });
+    sim.world_mut().insert_resource(VelocityCap {value: f64::MAX});
+
+    // Run the simulation for a number of steps.
+    for _i in 0..1_000 {
+        sim.update();
+    }
+    println!("Simulation completed in {} ms.", now.elapsed().as_millis());
+}
diff --git a/examples/thermalization.rs b/examples/thermalization.rs
new file mode 100644
index 00000000..9262896d
--- /dev/null
+++ b/examples/thermalization.rs
@@ -0,0 +1,133 @@
+extern crate atomecs as lib;
+extern crate nalgebra;
+use nalgebra::Vector3;
+use rand_distr::{Distribution, Uniform};
+use bevy::prelude::*;
+use lib::atom::{Atom, Force, Mass, Position, Velocity};
+use lib::initiate::NewlyCreated;
+use lib::integrator::Timestep;
+use lib::output::file::FileOutputPlugin;
+use lib::output::file::Text;
+use lib::simulation::SimulationBuilder;
+use std::fs::File;
+use std::io::{Error, Write};
+use std::time::Instant;
+use lib::shapes::Sphere as MySphere;
+use lib::collisions::wall_collisions::{WallData, WallType};
+use lib::collisions::CollisionPlugin;
+use lib::collisions::atom_collisions::{CollisionParameters, CollisionsTracker};
+use lib::sim_region::{SimulationVolume, VolumeType};
+
+fn main() {
+    let now = Instant::now();
+    let mut sim_builder = SimulationBuilder::default();
+    sim_builder.add_plugins(FileOutputPlugin::<Position, Text, Atom>::new("pos.txt".to_string(),150));
+    sim_builder.add_plugins(FileOutputPlugin::<Velocity, Text, Atom>::new("vel.txt".to_string(),150));
+    sim_builder.add_plugins(CollisionPlugin);
+
+    let mut sim = sim_builder.build();
+
+    let p_dist = Uniform::new(-25e-3, 25e-3).unwrap();
+    let v_dist = Uniform::new(-5e-1, 5e-1).unwrap();
+    for _i in 0..25000 {
+        sim.world_mut()
+            .spawn(Position {
+                pos: Vector3::new(
+                    p_dist.sample(&mut rand::rng()),
+                    p_dist.sample(&mut rand::rng()),
+                    p_dist.sample(&mut rand::rng()),
+                ),
+            })
+            .insert(Atom)
+            .insert(Force::default())
+            .insert(Velocity {
+                vel: Vector3::new(
+                    v_dist.sample(&mut rand::rng()),
+                    v_dist.sample(&mut rand::rng()),
+                    v_dist.sample(&mut rand::rng()),
+                ).normalize(),
+            })
+            .insert(NewlyCreated)
+            .insert(Mass { value: 87.0 });
+    }
+
+    sim.world_mut().insert_resource(CollisionParameters {
+        macroparticle: 5e8,
+        box_number: 100, //Any number large enough to cover entire cloud with collision boxes. Overestimating box number will not affect performance.
+        box_width: 1e-2, //Too few particles per box will both underestimate collision rate and cause large statistical fluctuations.
+        //Boxes must also be smaller than typical length scale of density variations within the cloud, since the collisions model treats gas within a box as homogeneous.
+        sigma: 3.5e-16, //Approximate collisional cross section of Rb87
+        collision_limit: 10_000_000.0, //Maximum number of collisions that can be calculated in one frame.
+                                       //This avoids absurdly high collision numbers if many atoms are initialised with the same position, for example.
+    });
+
+    sim.world_mut().insert_resource(CollisionsTracker {
+        num_collisions: Vec::new(),
+        num_atoms: Vec::new(),
+        num_particles: Vec::new(),
+    });
+
+    sim.world_mut()
+        .spawn(WallData{
+            wall_type: WallType::Rough})
+        .insert(MySphere{radius: 5e-2})
+        .insert(Position{pos: Vector3::new(0.0, 0.0, 0.0)});
+
+    sim.world_mut()
+        .spawn(Position {
+            pos: Vector3::new(0.0, 0.0, 0.0),
+        })
+        .insert(MySphere {
+            radius: 502e-4, 
+        })
+        .insert(SimulationVolume {
+            volume_type: VolumeType::Inclusive,
+        });
+
+    // Define timestep
+    sim.world_mut().insert_resource(Timestep { delta: 5e-5 });
+
+    let mut filename = File::create("collisions.txt").expect("Cannot create file.");
+
+    // Run the simulation for a number of steps.
+    for _i in 0..25_000 {
+        sim.update();
+
+        if (_i > 0) && (_i % 50_i32 == 0) {
+            let tracker = sim.world().get_resource::<CollisionsTracker>().unwrap();
+            let _result = write_collisions_tracker(
+                &mut filename,
+                &_i,
+                &tracker.num_collisions,
+                &tracker.num_atoms,
+                &tracker.num_particles,
+            )
+            .expect("Could not write collision stats file.");
+        }
+    }
+    println!("Simulation completed in {} ms.", now.elapsed().as_millis());
+}
+
+// Write collision stats to file
+
+fn write_collisions_tracker(
+    filename: &mut File,
+    step: &i32,
+    num_collisions: &Vec<i32>,
+    num_atoms: &Vec<f64>,
+    num_particles: &Vec<i32>,
+) -> Result<(), Error> {
+    let str_collisions: Vec<String> = num_collisions.iter().map(|n| n.to_string()).collect();
+    let str_atoms: Vec<String> = num_atoms.iter().map(|n| format!("{:.2}", n)).collect();
+    let str_particles: Vec<String> = num_particles.iter().map(|n| n.to_string()).collect();
+    write!(
+        filename,
+        "{:?}\r\n{:}\r\n{:}\r\n{:}\r\n",
+        step,
+        str_collisions.join(" "),
+        str_atoms.join(" "),
+        str_particles.join(" ")
+    )?;
+    Ok(())
+}
+
diff --git a/examples/top_trap_with_collisions.rs b/examples/top_trap_with_collisions.rs
index a6b89235..c3fa25b1 100644
--- a/examples/top_trap_with_collisions.rs
+++ b/examples/top_trap_with_collisions.rs
@@ -4,12 +4,14 @@ extern crate atomecs as lib;
 extern crate nalgebra;
 use lib::atom::{Atom, Force, Mass, Position, Velocity};
 use lib::collisions::CollisionPlugin;
-use lib::collisions::{ApplyCollisionsOption, CollisionParameters, CollisionsTracker};
+use lib::collisions::atom_collisions::{CollisionParameters, CollisionsTracker};
+use lib::collisions::{ApplyAtomCollisions, ApplyWallCollisions};
 use lib::initiate::NewlyCreated;
 use lib::integrator::Timestep;
 use lib::magnetic::force::{MagneticDipole};
 use lib::magnetic::quadrupole::QuadrupoleField3D;
 use lib::magnetic::top::UniformFieldRotator;
+use lib::magnetic::uniform::UniformMagneticField;
 use lib::output::file::FileOutputPlugin;
 use lib::output::file::Text;
 use lib::simulation::SimulationBuilder;
@@ -38,7 +40,9 @@ fn main() {
         .spawn(UniformFieldRotator{
             amplitude: 2e-3,
             frequency: 3000.0,
-        }); // Time averaged TOP theory assumes rotation frequency much greater than velocity of atoms
+        })
+        .insert(UniformMagneticField {field: Vector3::new(0.0, 0.0, 0.0)})
+        .insert(Position::default()); // Time averaged TOP theory assumes rotation frequency much greater than velocity of atoms
 
     let p_dist = Normal::new(0.0, 50e-6).unwrap();
     let v_dist = Normal::new(0.0, 0.004).unwrap(); // ~100nK
@@ -66,7 +70,8 @@ fn main() {
             .insert(Mass { value: 87.0 });
     }
 
-    sim.world_mut().insert_resource(ApplyCollisionsOption{apply_collision: true });
+    sim.world_mut().insert_resource(ApplyAtomCollisions(true));
+    sim.world_mut().insert_resource(ApplyWallCollisions(false));
     sim.world_mut().insert_resource(CollisionParameters {
         macroparticle: 4e2,
         box_number: 200, //Any number large enough to cover entire cloud with collision boxes. Overestimating box number will not affect performance.
diff --git a/src/atom.rs b/src/atom.rs
index 41c4a832..3735158e 100644
--- a/src/atom.rs
+++ b/src/atom.rs
@@ -22,7 +22,6 @@ pub struct Position {
 impl Default for Position {
     fn default() -> Self {
         Position {
-            /// position in 3D in units of m
             pos: Vector3::new(0.0, 0.0, 0.0),
         }
     }
diff --git a/src/atom_sources/emit.rs b/src/atom_sources/emit.rs
index a73498ba..832664f3 100644
--- a/src/atom_sources/emit.rs
+++ b/src/atom_sources/emit.rs
@@ -35,7 +35,6 @@ pub struct AtomNumberToEmit {
 
 
 /// Calculates the number of atoms to emit per frame for fixed atoms-per-timestep ovens
-/// Not doing parallel iteration here, expect few ovens.
 pub fn emit_number_per_frame_system(
     mut query: Query<(&EmitNumberPerFrame, &mut AtomNumberToEmit)>,
 ) {
@@ -77,17 +76,8 @@ pub fn emit_once_system(
 }
 
 #[cfg(test)]
-pub mod tests {
-    // These imports are actually needed! The compiler is getting confused and warning they are not.
-    #[allow(unused_imports)]
+mod tests {
     use super::*;
-    // extern crate specs;
-    #[allow(unused_imports)]
-    // use specs::{Builder, Entity, RunNow, World};
-    use bevy::prelude::*;
-    extern crate nalgebra;
-    #[allow(unused_imports)]
-    use nalgebra::Vector3;
 
     #[test]
     fn test_fixed_rate_emitter() {
diff --git a/src/atom_sources/mass.rs b/src/atom_sources/mass.rs
index 09df50f1..834ce560 100644
--- a/src/atom_sources/mass.rs
+++ b/src/atom_sources/mass.rs
@@ -48,7 +48,7 @@ impl MassDistribution {
             total += mr.ratio;
         }
 
-        for mut mr in &mut self.distribution {
+        for mr in &mut self.distribution {
             mr.ratio /= total;
         }
         self.normalised = true
@@ -74,10 +74,9 @@ impl MassDistribution {
     }
 }
 
-pub mod tests {
-    #[allow(unused_imports)]
+#[cfg(test)]
+mod tests {
     use super::*;
-    #[allow(unused_imports)]
     use assert_approx_eq::assert_approx_eq;
 
     #[test]
diff --git a/src/atom_sources/oven.rs b/src/atom_sources/oven.rs
index 52f493f5..822c9935 100644
--- a/src/atom_sources/oven.rs
+++ b/src/atom_sources/oven.rs
@@ -7,7 +7,6 @@ use std::marker::PhantomData;
 use super::emit::AtomNumberToEmit;
 use super::precalc::{MaxwellBoltzmannSource, PrecalculatedSpeciesInformation};
 use super::species::{AtomCreator};
-use crate::constant;
 use crate::constant::PI;
 use crate::initiate::*;
 
@@ -132,12 +131,12 @@ pub struct Oven<T> where T : AtomCreator {
     pub direction: Vector3<f64>,
 
     /// Angular distribution for atoms emitted by the oven.
-    theta_distribution: WeightedProbabilityDistribution,
+    pub theta_distribution: WeightedProbabilityDistribution,
 
     /// The maximum angle theta at which atoms can be emitted from the oven. This can be constricted eg by a heat shield, or 'hot lip'.
     pub max_theta: f64,
 
-    phantom : PhantomData<T>
+    pub phantom : PhantomData<T>
 }
 impl<T> MaxwellBoltzmannSource for Oven<T> where T : AtomCreator {
     fn get_temperature(&self) -> f64 {
@@ -161,12 +160,17 @@ impl<T> Oven<T> where T : AtomCreator {
             }
             OvenAperture::Circular { radius, thickness } => {
                 let dir = self.direction.normalize();
-                let dir_1 = dir.cross(&Vector3::new(2.0, 1.0, 0.5)).normalize();
+                let dir_1 = dir.cross(&Vector3::new(2.3, 1.2, 0.3)).normalize();
                 let dir_2 = dir.cross(&dir_1).normalize();
-                let theta = rng.random_range(0.0..2. * constant::PI);
-                let r = rng.random_range(0.0..radius);
+                let (x, y) = loop {
+                    let x = rng.random_range(-radius..radius);
+                    let y = rng.random_range(-radius..radius);
+                    if x.powi(2) + y.powi(2) < radius.powi(2){
+                        break (x, y);
+                    }
+                };
                 let h = rng.random_range(-0.5 * thickness..0.5 * thickness);
-                dir * h + r * dir_1 * theta.sin() + r * dir_2 * theta.cos()
+                dir * h + x * dir_1 + y * dir_2
             }
         }
     }
@@ -294,5 +298,3 @@ fn create_jtheta_distribution(
 
     WeightedProbabilityDistribution::new(thetas, weights)
 }
-
-
diff --git a/src/atom_sources/precalc.rs b/src/atom_sources/precalc.rs
index 51b709ff..5b99aa4a 100644
--- a/src/atom_sources/precalc.rs
+++ b/src/atom_sources/precalc.rs
@@ -8,7 +8,6 @@ use rand;
 use rand::distr::Distribution;
 use rand::distr::weighted::WeightedIndex;
 use rand::Rng;
-use std::marker::PhantomData;
 
 use bevy::prelude::*;
 
diff --git a/src/collisions.rs b/src/collisions/atom_collisions.rs
similarity index 56%
rename from src/collisions.rs
rename to src/collisions/atom_collisions.rs
index bc8696e9..62b615bd 100644
--- a/src/collisions.rs
+++ b/src/collisions/atom_collisions.rs
@@ -20,21 +20,7 @@ use hashbrown::HashMap;
 use nalgebra::Vector3;
 use rand::Rng;
 use bevy::prelude::*;
-use crate::integrator::IntegrationSet;
-use crate::integrator::AtomECSBatchStrategy;
-
-/// A resource that indicates that the simulation should apply scattering
-#[derive(Resource)]
-pub struct ApplyCollisionsOption{
-    pub apply_collision: bool,
-}
-
-/// Component that marks which box an atom is in for spatial partitioning
-#[derive(Component)]
-pub struct BoxID {
-    /// ID of the box
-    pub id: i64,
-}
+use crate::collisions::spatial_grid::BoxID;
 
 /// A patition of space within which collisions can occur
 pub struct CollisionBox<> {
@@ -149,86 +135,54 @@ pub struct CollisionsTracker {
     pub num_atoms: Vec<f64>,
 }
 
-/// Initializes boxid component
-pub fn init_boxid_system(
-    mut query: Query<(Entity, &Position), (Without<BoxID>, With<Atom>)>,
-    collisions_option: Res<ApplyCollisionsOption>,
-    mut commands: Commands,
-) {
-    match collisions_option.apply_collision {
-        false => return, // No need to initialize box IDs if collisions are not applied
-        true => {
-            query
-                .iter_mut()
-                .for_each(|(entity, _)| {
-                    commands.entity(entity).insert(BoxID { id: 0 });
-                });},
-    }
-}
-
 /// Performs collisions within the atom cloud using a spatially partitioned Monte-Carlo approach.
 pub fn apply_collisions_system(
     mut query: Query<(Entity, &Position, &mut Velocity, &mut BoxID), With<Atom>>,
-    collisions_option: Res<ApplyCollisionsOption>,
     timestep: Res<Timestep>,
     params: Res<CollisionParameters>,
     mut tracker: ResMut<CollisionsTracker>,
-    batch_strategy: Res<AtomECSBatchStrategy>,
 ) {
     use rayon::prelude::*;
-    match collisions_option.apply_collision {
-        false => (),
-        true => {
-            // Assign box IDs to atoms
-            query
-                .par_iter_mut()
-                .batching_strategy(batch_strategy.0.clone())
-                .for_each(|(_, position, _, mut boxid)| {
-                    boxid.id = pos_to_id(position.pos, params.box_number, params.box_width);
-                });
-
-            // insert atom velocity into hash
-            let mut map: HashMap<i64, CollisionBox> = HashMap::new();
-            for (entity, _, velocity, boxid) in query.iter() {
-                if boxid.id == i64::MAX {
-                    continue;
-                } else {
-                    map.entry(boxid.id).or_default()
-                        .entity_velocities.push((entity, velocity.vel));
-                }
-            }
+    // insert atom velocity into hash
+    let mut map: HashMap<i64, CollisionBox> = HashMap::new();
+    for (entity, _, velocity, boxid) in query.iter() {
+        if boxid.id == i64::MAX {
+            continue;
+        } else {
+            map.entry(boxid.id).or_default()
+                .entity_velocities.push((entity, velocity.vel));
+        }
+    }
 
-            // get immutable list of boxes and iterate in parallel
-            // (Note that using hashmap parallel values mut does not work in parallel, tested.)
-            let boxes: Vec<&mut CollisionBox> = map.values_mut().collect();
-            boxes.into_par_iter().for_each(|collision_box| {
-                collision_box.do_collisions(*params, timestep.delta);
-            });
-
-            // Write back values.
-            for collision_box in map.values() {
-                for &(entity, new_velocity) in &collision_box.entity_velocities {
-                    if let Ok((_, _, mut velocity, _)) = query.get_mut(entity) {
-                        velocity.vel = new_velocity;
-                    }
-                }
+    // get immutable list of boxes and iterate in parallel
+    // (Note that using hashmap parallel values mut does not work in parallel, tested.)
+    let boxes: Vec<&mut CollisionBox> = map.values_mut().collect();
+    boxes.into_par_iter().for_each(|collision_box| {
+        collision_box.do_collisions(*params, timestep.delta);
+    });
+
+    // Write back values.
+    for collision_box in map.values() {
+        for &(entity, new_velocity) in &collision_box.entity_velocities {
+            if let Ok((_, _, mut velocity, _)) = query.get_mut(entity) {
+                velocity.vel = new_velocity;
             }
-
-            // Update tracker
-            tracker.num_atoms = map
-                .values()
-                .map(|collision_box| collision_box.atom_number)
-                .collect();
-            tracker.num_collisions = map
-                .values()
-                .map(|collision_box| collision_box.collision_number)
-                .collect();
-            tracker.num_particles = map
-                .values()
-                .map(|collision_box| collision_box.particle_number)
-                .collect();
         }
     }
+
+    // Update tracker
+    tracker.num_atoms = map
+        .values()
+        .map(|collision_box| collision_box.atom_number)
+        .collect();
+    tracker.num_collisions = map
+        .values()
+        .map(|collision_box| collision_box.collision_number)
+        .collect();
+    tracker.num_particles = map
+        .values()
+        .map(|collision_box| collision_box.particle_number)
+        .collect();
 }
 
 
@@ -254,102 +208,25 @@ fn do_collision(mut v1: Vector3<f64>, mut v2: Vector3<f64>) -> (Vector3<f64>, Ve
     (v1, v2)
 }
 
-fn pos_to_id(pos: Vector3<f64>, n: i64, width: f64) -> i64 {
-    //Assume that atoms that leave the grid are too sparse to collide, so disregard them
-    //We'll assign them the max value of i64, and then check for this value when we do a collision and ignore them
-    let bound = (n as f64) / 2.0 * width;
-
-    let id: i64;
-    if pos[0].abs() > bound || pos[1].abs() > bound || pos[2].abs() > bound {
-        id = i64::MAX;
-    } else {
-        let xp: i64;
-        let yp: i64;
-        let zp: i64;
-
-        // even number of boxes, vertex of a box is on origin
-        // odd number of boxes, centre of a box is on the origin
-        // grid cells run from [0, width), i.e include lower bound but exclude upper
-
-        xp = (pos[0] / width + 0.5 * (n as f64)).floor() as i64;
-        yp = (pos[1] / width + 0.5 * (n as f64)).floor() as i64;
-        zp = (pos[2] / width + 0.5 * (n as f64)).floor() as i64;
-        //convert position to box id
-        id = xp + n * yp + n.pow(2) * zp;
-    }
-
-    id
-}
-
-pub struct CollisionPlugin;
-impl Plugin for CollisionPlugin {
-    fn build(&self, app: &mut App) {
-        // Note that the collisions system must be applied after the velocity integrator or it will violate conservation of energy and cause heating
-        app.add_systems(PostUpdate, init_boxid_system.after(IntegrationSet::EndIntegration));
-        app.add_systems(PostUpdate, apply_collisions_system.after(IntegrationSet::EndIntegration).after(init_boxid_system));
-    }
-}
-
-pub mod tests {
-    #[allow(unused_imports)]
+#[cfg(test)]
+mod tests {
     use super::*;
-    #[allow(unused_imports)]
     use crate::atom::{Atom, Force, Mass, Position, Velocity};
-    #[allow(unused_imports)]
     use crate::initiate::NewlyCreated;
-    #[allow(unused_imports)]
-    use crate::integrator::{
-        Step, Timestep,
-    };
-
-    #[allow(unused_imports)]
+    use crate::integrator::Timestep;
+    use crate::collisions::CollisionPlugin;
     use nalgebra::Vector3;
-    #[allow(unused_imports)]
-    use bevy::prelude::*;
-    #[allow(unused_imports)]
     use crate::simulation::SimulationBuilder;
-    extern crate bevy;
-
-    #[test]
-    fn test_pos_to_id() {
-        let n: i64 = 10;
-        let width: f64 = 2.0;
-
-        let pos1 = Vector3::new(0.0, 0.0, 0.0);
-        let pos2 = Vector3::new(1.0, 0.0, 0.0);
-        let pos3 = Vector3::new(2.0, 0.0, 0.0);
-        let pos4 = Vector3::new(9.9, 0.0, 0.0);
-        let pos5 = Vector3::new(-9.9, 0.0, 0.0);
-        let pos6 = Vector3::new(10.1, 0.0, 0.0);
-        let pos7 = Vector3::new(-9.9, -9.9, -9.9);
-
-        let id1 = pos_to_id(pos1, n, width);
-        let id2 = pos_to_id(pos2, n, width);
-        let id3 = pos_to_id(pos3, n, width);
-        let id4 = pos_to_id(pos4, n, width);
-        let id5 = pos_to_id(pos5, n, width);
-        let id6 = pos_to_id(pos6, n, width);
-        let id7 = pos_to_id(pos7, n, width);
-
-        assert_eq!(id1, 555);
-        assert_eq!(id2, 555);
-        assert_eq!(id3, 556);
-        assert_eq!(id4, 559);
-        assert_eq!(id5, 550);
-        assert_eq!(id6, i64::MAX);
-        assert_eq!(id7, 0);
-    }
 
     #[test]
     fn test_do_collision() {
         // do this test muliple times since there is a random element involved in do_collision
+        let v1 = Vector3::new(0.5, 1.0, 0.75);
+        let v2 = Vector3::new(0.2, 0.0, 1.25);
+        //calculate energy and momentum before
+        let ptoti = v1 + v2;
+        let energyi = 0.5 * (v1.norm_squared() + v2.norm_squared());
         for _i in 0..50 {
-            let v1 = Vector3::new(0.5, 1.0, 0.75);
-            let v2 = Vector3::new(0.2, 0.0, 1.25);
-            //calculate energy and momentum before
-            let ptoti = v1 + v2;
-            let energyi = 0.5 * (v1.norm_squared() + v2.norm_squared());
-
             let (v1new, v2new) = do_collision(v1, v2);
 
             //energy and momentum after
@@ -366,44 +243,44 @@ pub mod tests {
     /// Test that the expected number of collisions in a CollisionBox is correct.
     #[test]
     fn collision_rate() {
-    use assert_approx_eq::assert_approx_eq;
-    let vel = Vector3::new(1.0, 0.0, 0.0);
-    const MACRO_ATOM_NUMBER: usize = 100;
-    
-    // Create dummy entities for testing
-    let mut entity_velocities: Vec<(Entity, Vector3<f64>)> = Vec::new();
-    for i in 0..MACRO_ATOM_NUMBER {
-        let entity = Entity::PLACEHOLDER; // Or just use the same placeholder for all
-        entity_velocities.push((entity, vel));
+        use assert_approx_eq::assert_approx_eq;
+        let vel = Vector3::new(1.0, 0.0, 0.0);
+        const MACRO_ATOM_NUMBER: usize = 100;
+        
+        // Create dummy entities for testing
+        let mut entity_velocities: Vec<(Entity, Vector3<f64>)> = Vec::new();
+        for _i in 0..MACRO_ATOM_NUMBER {
+            let entity = Entity::PLACEHOLDER; // Or just use the same placeholder for all
+            entity_velocities.push((entity, vel));
+        }
+        
+        let mut collision_box = CollisionBox {
+            entity_velocities,
+            ..Default::default()
+        };
+        
+        let params = CollisionParameters {
+            macroparticle: 10.0,
+            box_number: 1,
+            box_width: 1e-3,
+            sigma: 1e-8,
+            collision_limit: 10_000.0,
+        };
+        let dt = 1e-3;
+        collision_box.do_collisions(params, dt);
+        
+        assert_eq!(collision_box.particle_number, MACRO_ATOM_NUMBER as i32);
+        let atom_number = params.macroparticle * MACRO_ATOM_NUMBER as f64;
+        assert_eq!(collision_box.atom_number, atom_number);
+        let density = atom_number / params.box_width.powi(3);
+        let expected_number =
+            (1.0 / SQRT2) * MACRO_ATOM_NUMBER as f64 * density * params.sigma * vel.norm() * dt;
+        assert_approx_eq!(
+            collision_box.expected_collision_number,
+            expected_number,
+            0.01
+        );
     }
-    
-    let mut collision_box = CollisionBox {
-        entity_velocities,
-        ..Default::default()
-    };
-    
-    let params = CollisionParameters {
-        macroparticle: 10.0,
-        box_number: 1,
-        box_width: 1e-3,
-        sigma: 1e-8,
-        collision_limit: 10_000.0,
-    };
-    let dt = 1e-3;
-    collision_box.do_collisions(params, dt);
-    
-    assert_eq!(collision_box.particle_number, MACRO_ATOM_NUMBER as i32);
-    let atom_number = params.macroparticle * MACRO_ATOM_NUMBER as f64;
-    assert_eq!(collision_box.atom_number, atom_number);
-    let density = atom_number / params.box_width.powi(3);
-    let expected_number =
-        (1.0 / SQRT2) * MACRO_ATOM_NUMBER as f64 * density * params.sigma * vel.norm() * dt;
-    assert_approx_eq!(
-        collision_box.expected_collision_number,
-        expected_number,
-        0.01
-    );
-}
 
 
     /// Test that the system runs and causes nearby atoms to collide. More of an integration test than a unit test.
@@ -444,7 +321,6 @@ pub mod tests {
 
         let dt = 1.0;
         sim.world_mut().insert_resource(Timestep { delta: dt });
-        sim.world_mut().insert_resource(ApplyCollisionsOption{apply_collision: true});
         sim.world_mut().insert_resource(CollisionsTracker {
             num_collisions: Vec::new(),
             num_atoms: Vec::new(),
diff --git a/src/collisions/mod.rs b/src/collisions/mod.rs
new file mode 100644
index 00000000..5b9c7cd4
--- /dev/null
+++ b/src/collisions/mod.rs
@@ -0,0 +1,144 @@
+pub mod atom_collisions;
+pub mod wall_collisions;
+pub mod spatial_grid;
+pub mod wall_collision_info;
+
+use bevy::prelude::*;
+use std::fmt;
+use crate::collisions::{atom_collisions::*, wall_collisions::*, spatial_grid::*};
+use crate::atom::Atom;
+use crate::integrator::IntegrationSet;
+use crate::shapes::{
+    Cylinder as MyCylinder,
+    Cuboid as MyCuboid,
+    Sphere as MySphere,
+    CylindricalPipe,
+};
+use rand;
+use rand::distr::Distribution;
+use rand::distr::weighted::WeightedIndex;
+use rand::Rng;
+
+/// A resource that indicates that the simulation should apply atom collisions
+#[derive(Resource)]
+pub struct ApplyAtomCollisions(pub bool);
+
+/// A resource that indicates that the simulation should apply wall collisions
+#[derive(Resource)]
+pub struct ApplyWallCollisions(pub bool);
+
+#[derive(Debug, Hash, PartialEq, Eq, Clone, SystemSet)]
+pub enum CollisionsSet {
+    Set,
+    WallCollisionSystems,
+    AtomCollisionSystems,
+}
+
+fn apply_atom_collisions(apply_atom_collision: Res<ApplyAtomCollisions>) -> bool {
+    apply_atom_collision.0
+}
+
+fn apply_wall_collisions(apply_wall_collision: Res<ApplyWallCollisions>) -> bool {
+    apply_wall_collision.0
+}
+
+#[derive(Component, Clone)]
+pub struct NumberOfWallCollisions {
+    pub value: i32
+}
+
+impl fmt::Display for NumberOfWallCollisions {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(f, "{:?}", self.value)
+    }
+}
+
+/// Initializes number of wall collisions component
+pub fn init_number_of_collisions_system (
+    mut query: Query<Entity, (With<Atom>, Without<NumberOfWallCollisions>)>,
+    mut commands: Commands,
+) {
+    for atom_entity in query.iter_mut() {
+        commands.entity(atom_entity).insert(NumberOfWallCollisions {
+            value: 0
+        });
+    }
+}
+
+pub struct CollisionPlugin;
+impl Plugin for CollisionPlugin {
+    fn build(&self, app: &mut App) {
+        app.world_mut().insert_resource(ApplyAtomCollisions(true));
+        app.world_mut().insert_resource(ApplyWallCollisions(true));
+        app.world_mut().insert_resource(MaxSteps(1000));
+        app.world_mut().insert_resource(SurfaceThreshold(1e-9));
+
+        // Atom Collision Systems
+        // Note that the atom collisions system must be applied after the velocity integrator or it will violate conservation of energy and cause heating
+        app.add_systems(PostUpdate, init_boxid_system
+            .in_set(CollisionsSet::AtomCollisionSystems)
+            .after(IntegrationSet::EndIntegration)
+            .run_if(apply_atom_collisions));
+        app.add_systems(PostUpdate, assign_boxid_system
+            .in_set(CollisionsSet::AtomCollisionSystems)
+            .after(init_boxid_system)
+            .run_if(apply_atom_collisions));
+        app.add_systems(PostUpdate, apply_collisions_system
+            .in_set(CollisionsSet::AtomCollisionSystems)
+            .after(assign_boxid_system)
+            .run_if(apply_atom_collisions));
+        
+        // Wall Collisions Systems
+        app.add_systems(Startup, create_cosine_distribution
+            .in_set(CollisionsSet::Set)
+            .run_if(apply_wall_collisions));
+        app.add_systems(PreUpdate, assign_location_status_system
+            .in_set(CollisionsSet::WallCollisionSystems)
+            .run_if(apply_wall_collisions));
+        app.add_systems(PreUpdate, init_number_of_collisions_system
+            .in_set(CollisionsSet::WallCollisionSystems)
+            .run_if(apply_wall_collisions));
+        app.add_systems(PreUpdate, update_distance_to_travel_system
+            .in_set(CollisionsSet::WallCollisionSystems)
+            .run_if(apply_wall_collisions));
+        app.add_systems(PreUpdate, wall_collision_system::<MySphere>
+            .in_set(CollisionsSet::WallCollisionSystems)
+            .after(IntegrationSet::BeginIntegration)
+            .run_if(apply_wall_collisions));
+        app.add_systems(PreUpdate, wall_collision_system::<MyCuboid>
+            .in_set(CollisionsSet::WallCollisionSystems)
+            .after(IntegrationSet::BeginIntegration)
+            .run_if(apply_wall_collisions));
+        app.add_systems(PreUpdate, wall_collision_system::<MyCylinder>
+            .in_set(CollisionsSet::WallCollisionSystems)
+            .after(IntegrationSet::BeginIntegration)
+            .run_if(apply_wall_collisions));
+        app.add_systems(PreUpdate, wall_collision_system::<CylindricalPipe>
+            .in_set(CollisionsSet::WallCollisionSystems)
+            .after(IntegrationSet::BeginIntegration)
+            .run_if(apply_wall_collisions));
+    }
+}
+
+/// A copy of weighted probability distribution, used for the Lambertian cosine distribution
+#[derive(Resource)]
+pub struct LambertianProbabilityDistribution {
+    values: Vec<f64>,
+    weighted_index: WeightedIndex<f64>,
+}
+
+impl LambertianProbabilityDistribution {
+    pub fn new(values: Vec<f64>, weights: Vec<f64>) -> Self {
+        LambertianProbabilityDistribution {
+            values,
+            weighted_index: WeightedIndex::new(&weights).unwrap(),
+        }
+    }
+}
+
+impl Distribution<f64> for LambertianProbabilityDistribution {
+    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
+        let index = self.weighted_index.sample(rng);
+        self.values[index]
+    }
+}
\ No newline at end of file
diff --git a/src/collisions/spatial_grid.rs b/src/collisions/spatial_grid.rs
new file mode 100644
index 00000000..f95bcc9e
--- /dev/null
+++ b/src/collisions/spatial_grid.rs
@@ -0,0 +1,102 @@
+use bevy::prelude::*;
+use crate::integrator::AtomECSBatchStrategy;
+use crate::collisions::atom_collisions::CollisionParameters;
+use crate::atom::{Position, Atom};
+use nalgebra::Vector3;
+
+/// Component that marks which box an atom is in for spatial partitioning
+#[derive(Component)]
+pub struct BoxID {
+    /// ID of the box
+    pub id: i64,
+}
+
+/// Initializes boxid component
+pub fn init_boxid_system(
+    mut query: Query<Entity, (Without<BoxID>, With<Atom>)>,
+    mut commands: Commands,
+) {
+    query
+        .iter_mut()
+        .for_each(|entity| {
+            commands.entity(entity).insert(BoxID { id: 0 });
+        });
+}
+
+
+/// Assigns box IDs to atoms
+pub fn assign_boxid_system(
+    mut query: Query<(&Position, &mut BoxID), With<Atom>>,
+    params: Res<CollisionParameters>,
+    batch_strategy: Res<AtomECSBatchStrategy>,
+) {
+    query
+        .par_iter_mut()
+        .batching_strategy(batch_strategy.0.clone())
+        .for_each(|(position, mut boxid)| {
+            boxid.id = pos_to_id(position.pos, params.box_number, params.box_width);
+        });
+}
+
+
+fn pos_to_id(pos: Vector3<f64>, n: i64, width: f64) -> i64 {
+    //Assume that atoms that leave the grid are too sparse to collide, so disregard them
+    //We'll assign them the max value of i64, and then check for this value when we do a collision and ignore them
+    let bound = (n as f64) / 2.0 * width;
+
+    let id: i64;
+    if pos[0].abs() > bound || pos[1].abs() > bound || pos[2].abs() > bound {
+        id = i64::MAX;
+    } else {
+        let xp: i64;
+        let yp: i64;
+        let zp: i64;
+
+        // even number of boxes, vertex of a box is on origin
+        // odd number of boxes, centre of a box is on the origin
+        // grid cells run from [0, width), i.e include lower bound but exclude upper
+
+        xp = (pos[0] / width + 0.5 * (n as f64)).floor() as i64;
+        yp = (pos[1] / width + 0.5 * (n as f64)).floor() as i64;
+        zp = (pos[2] / width + 0.5 * (n as f64)).floor() as i64;
+        //convert position to box id
+        id = xp * 9803 + n * yp * 5213 + n.pow(2) * zp * 7789;
+    }
+
+    id
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_pos_to_id() {
+        let n: i64 = 10;
+        let width: f64 = 2.0;
+
+        let pos1 = Vector3::new(0.0, 0.0, 0.0);
+        let pos2 = Vector3::new(1.0, 0.0, 0.0);
+        let pos3 = Vector3::new(2.0, 0.0, 0.0);
+        let pos4 = Vector3::new(9.9, 0.0, 0.0);
+        let pos5 = Vector3::new(-9.9, 0.0, 0.0);
+        let pos6 = Vector3::new(10.1, 0.0, 0.0);
+        let pos7 = Vector3::new(-9.9, -9.9, -9.9);
+
+        let id1 = pos_to_id(pos1, n, width);
+        let id2 = pos_to_id(pos2, n, width);
+        let id3 = pos_to_id(pos3, n, width);
+        let id4 = pos_to_id(pos4, n, width);
+        let id5 = pos_to_id(pos5, n, width);
+        let id6 = pos_to_id(pos6, n, width);
+        let id7 = pos_to_id(pos7, n, width);
+
+        assert_eq!(id1, 4204165);
+        assert_eq!(id2, 4204165);
+        assert_eq!(id3, 4213968);
+        assert_eq!(id4, 4243377);
+        assert_eq!(id5, 4155150);
+        assert_eq!(id6, i64::MAX);
+        assert_eq!(id7, 0);
+    }
+}
\ No newline at end of file
diff --git a/src/collisions/wall_collision_info.rs b/src/collisions/wall_collision_info.rs
new file mode 100644
index 00000000..9b12f34b
--- /dev/null
+++ b/src/collisions/wall_collision_info.rs
@@ -0,0 +1,600 @@
+// An attempt at making wall_collision not so horrendously big
+
+use nalgebra::Vector3;
+use crate::shapes::{
+    Cylinder,
+    Cuboid,
+    Sphere, 
+    CylindricalPipe,
+};
+
+pub trait Intersect{
+    /// Calculate intersection point between atom trajectory and shape and return point of intersection
+    fn calculate_intersect(&self,
+        atom_pos: &Vector3<f64>,
+        atom_vel: &Vector3<f64>,
+        wall_pos: &Vector3<f64>,
+        time: f64,
+        tolerance: f64,
+        max_steps: i32,)
+        -> Option<Vector3<f64>>;
+}
+
+pub trait Normal{
+    /// Calculate normal at point of intersection (collision point)
+    fn calculate_normal(&self, point: &Vector3<f64>, wall_pos: &Vector3<f64>, tolerance: f64) -> Option<Vector3<f64>>;
+}
+
+trait SDF {
+    /// Calculate signed distance
+    fn signed_distance(&self, point: &Vector3<f64>) -> f64;
+}
+
+/// Credit to Inigo Quilez's for the SDFs. Found at https://www.shadertoy.com/view/Xds3zN
+
+/// Signed distance to a sphere.
+impl SDF for Sphere{
+    fn signed_distance(&self, point: &Vector3<f64>) -> f64 {
+        point.norm() - &self.radius
+    }
+}
+
+/// Signed distance to an axis-aligned box centered at the origin.
+impl SDF for Cuboid{
+    fn signed_distance(&self, point: &Vector3<f64>) -> f64 {
+        let delta = Vector3::new(
+            point.x.abs() - &self.half_width.x,
+            point.y.abs() - &self.half_width.y,
+            point.z.abs() - &self.half_width.z,
+        );
+
+        let outside = Vector3::new(
+            delta.x.max(0.0),
+            delta.y.max(0.0),
+            delta.z.max(0.0),
+        );
+
+        let outside_dist = outside.norm();
+        let inside_dist = f64::min(f64::max(delta.x, f64::max(delta.y, delta.z)), 0.0);
+
+        outside_dist + inside_dist
+    }
+}
+
+/// Signed distance to a capped cylinder.
+/// `cyl_start` and `cyl_end` define the axis of the cylinder (centre of the caps).
+impl SDF for Cylinder{
+    fn signed_distance(&self, point: &Vector3<f64>) -> f64 {
+        let cyl_start = -self.direction*self.length*0.5;
+        let local_point = point - cyl_start;
+
+        let point_dot_axis = local_point.dot(&self.direction);
+
+        let radial_dist = (local_point - self.direction * point_dot_axis).norm() - self.radius;
+
+        let axial_dist = (point_dot_axis - self.length * 0.5).abs() - self.length * 0.5;
+
+        let radial_sq = radial_dist * radial_dist;
+        let axial_sq = axial_dist * axial_dist;
+
+        let dist = if radial_dist.max(axial_dist) < 0.0 {
+            -radial_sq.min(axial_sq)
+        } else {
+            (if radial_dist > 0.0 { radial_sq } else { 0.0 }) + (if axial_dist > 0.0 { axial_sq } else { 0.0 })
+        };
+
+        dist.signum() * dist.abs().sqrt()
+    }
+}
+
+// Uses -velocity as direction, ignores effect of acceleration.
+impl<T : SDF> Intersect for T {
+    fn calculate_intersect(
+        &self,
+        atom_pos: &Vector3<f64>,
+        atom_vel: &Vector3<f64>,
+        wall_pos: &Vector3<f64>,
+        max_time: f64,
+        tolerance: f64,
+        max_steps: i32,
+    ) -> Option<Vector3<f64>> {
+        let speed = atom_vel.norm();
+        if speed == 0.0 {
+            // eprintln!("Speed 0 atom");
+            return None;
+        }
+
+        let direction = -atom_vel / speed;
+        let max_distance = speed * max_time;
+
+        let mut distance_traveled = 0.0;
+
+        for _ in 0..max_steps {
+            let curr_pos = atom_pos + direction * distance_traveled;
+            let local_pos = curr_pos - wall_pos;
+
+            let distance_to_surface = self.signed_distance(&local_pos);
+
+            if distance_to_surface < tolerance {
+                if distance_traveled + distance_to_surface < max_distance + tolerance {
+                    return Some(curr_pos);
+                }
+                else {
+                    // eprintln!("Took too long to collide");
+                    return None;
+                }
+            }
+
+            distance_traveled += distance_to_surface;
+        }
+        // eprintln!("Too few steps for convergence. collision");
+        None
+    }
+}
+
+impl Intersect for CylindricalPipe {
+    fn calculate_intersect(
+        &self,
+        atom_pos: &Vector3<f64>,
+        atom_vel: &Vector3<f64>,
+        wall_pos: &Vector3<f64>,
+        max_time: f64,
+        _tolerance: f64,
+        _max_steps: i32,
+    ) -> Option<Vector3<f64>> {
+        let prev_pos = atom_pos - atom_vel * max_time;
+        let delta_prev = prev_pos - wall_pos;
+        let delta_curr = atom_pos - wall_pos;
+
+        let axial_prev = delta_prev.dot(&self.direction);
+        let radial_prev = delta_prev - axial_prev * self.direction;
+        let axial_curr = delta_curr.dot(&self.direction);
+        let radial_curr = delta_curr - axial_curr * self.direction;
+
+        let a = (radial_curr - radial_prev).norm_squared();
+        let b = 2.0 * radial_prev.dot(&(&radial_curr - &radial_prev));
+        let c = radial_prev.norm_squared() - self.radius.powi(2);
+
+        // Solve quadratic equation for t
+        let discriminant = b * b - 4.0 * a * c;
+
+        if discriminant < 0.0 {
+            // No intersection, should never happen
+            return None;
+        }
+        let sqrt_discriminant = discriminant.sqrt();
+        let t1 = (-b + sqrt_discriminant) / (2.0 * a);
+        let t2 = (-b - sqrt_discriminant) / (2.0 * a);
+        if (t1 <= 0.0 || t1 >= 1.0 ) && (t2 <= 0.0 || t2 >= 1.0 ) {
+            return None;
+        }
+
+        let t = if t1 <= 0.0 || t1 >= 1.0  {
+            t2
+        } else if t2 <= 0.0 || t2 >= 1.0 {
+            t1
+        } else {
+            f64::min(t1, t2)
+        };
+
+        // Calculate the collision point
+        let collision_point = prev_pos * (1.0 - t) + atom_pos * t;
+        if (collision_point - wall_pos).dot(&self.direction).abs() <= self.length * 0.5 {
+            Some(collision_point)
+        } else {
+            None
+        }
+    }
+}
+
+impl Normal for Sphere {
+    fn calculate_normal(&self, point: &Vector3<f64>, wall_pos: &Vector3<f64>, tolerance: f64) -> Option<Vector3<f64>> {
+        let normal_vector = point - wall_pos;
+
+        // Gives the normal assuming collision from the inside
+        if (normal_vector.norm() - self.radius).abs() < tolerance  {
+            Some(-normal_vector.normalize())
+        } else { None } // Should never happen
+    }
+}
+
+impl Normal for Cuboid {
+    fn calculate_normal(&self, point: &Vector3<f64>, wall_pos: &Vector3<f64>, tolerance: f64) -> Option<Vector3<f64>> {
+        let local = point - wall_pos;
+
+        // Gives the normal assuming collision from the inside
+        if (local.x.abs() - self.half_width.x).abs() < tolerance {
+            Some(-Vector3::new(local.x.signum(), 0.0, 0.0))
+        }
+        else if (local.y.abs() - self.half_width.y).abs() < tolerance{
+            Some(-Vector3::new(0.0, local.y.signum(), 0.0))
+        }
+        else if (local.z.abs() - self.half_width.z).abs() < tolerance{
+            Some(-Vector3::new(0.0, 0.0, local.z.signum()))
+        }
+        else {
+            None
+        }
+    }
+}
+
+impl Normal for Cylinder {
+    fn calculate_normal(&self, point: &Vector3<f64>, wall_pos: &Vector3<f64>, tolerance: f64) -> Option<Vector3<f64>> {
+        // Convert to local space
+        let rel = point - wall_pos;
+        let local = Vector3::new(
+            rel.dot(&self.perp_x),
+            rel.dot(&self.perp_y),
+            rel.dot(&self.direction),
+        );
+
+        let axial_dist_to_surface = local.z.abs() - self.length * 0.5;
+        let radial_dist_to_surface = (local.x.powi(2) + local.y.powi(2)).sqrt() - self.radius;
+        // Gives the normal assuming collision from the inside
+        if axial_dist_to_surface.abs() < tolerance && radial_dist_to_surface < 0.0 {
+            let normal = self.direction * local.z.signum();
+            Some(-normal) 
+        } else if (radial_dist_to_surface).abs() < tolerance {
+            let normal = self.perp_x * local.x + self.perp_y * local.y;
+            Some(-normal.normalize()) 
+        }
+        else{
+            None
+        }
+    }
+}
+
+impl Normal for CylindricalPipe {
+    fn calculate_normal(&self, point: &Vector3<f64>, wall_pos: &Vector3<f64>, tolerance: f64) -> Option<Vector3<f64>> {
+        // Convert to local space
+        let rel = point - wall_pos;
+        let local = Vector3::new(
+            rel.dot(&self.perp_x),
+            rel.dot(&self.perp_y),
+            rel.dot(&self.direction),
+        );
+
+        let radial_dist_to_surface = (local.x.powi(2) + local.y.powi(2)).sqrt() - self.radius;
+        if (radial_dist_to_surface).abs() < tolerance {
+            let normal = self.perp_x * local.x + self.perp_y * local.y;
+            Some(-normal.normalize()) 
+        }
+        else{
+            None
+        }
+    }
+}
+
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::shapes::{Cylinder, Cuboid, Sphere, Volume};
+    use crate::constant::PI;
+    use rand::Rng;
+    use assert_approx_eq::assert_approx_eq;
+
+
+    const TEST_RUNS: usize = 10_000;
+    const TOLERANCE: f64 = 1e-9;
+    const MAX_DT: f64 = 1e3;
+    const MAX_STEPS: i32 = 1000;
+
+    fn rand_vec(min: f64, max: f64) -> Vector3<f64> {
+        let mut rng = rand::rng();
+        Vector3::new(
+            rng.random_range(min..=max),
+            rng.random_range(min..=max),
+            rng.random_range(min..=max),
+        )
+    }
+
+    #[test]
+    fn test_intersect_cuboid() {
+        let mut success = 0;
+        let mut failed = 0;
+        for _ in 0..TEST_RUNS {
+            let mut rng = rand::rng();
+            let length = 1.0;
+            let breadth = 1.5;
+            let height = 2.0;
+            let half_width = Vector3::new(length, breadth, height);
+
+            let cuboid = Cuboid { half_width };
+            let wall_pos = rand_vec(-5.0, 5.0);
+
+            let inside_local = Vector3::new(
+                rng.random_range(-half_width.x ..=half_width.x),
+                rng.random_range(-half_width.y ..=half_width.y),
+                rng.random_range(-half_width.z..=half_width.z),
+            );
+
+            let velocity = rand_vec(-1.0, 1.0).normalize();
+            let outside_local = inside_local + velocity*10.0;
+
+            let inside = wall_pos + inside_local;
+            let outside = wall_pos + outside_local;
+
+            if cuboid.contains(&wall_pos, &outside) || !cuboid.contains(&wall_pos, &inside) {
+                panic!("Wrong initialization in test");
+            }
+
+            let result = cuboid.calculate_intersect(&outside, &velocity, &wall_pos, MAX_DT, TOLERANCE, MAX_STEPS);
+
+            match result {
+                Some(hit) => {
+                    let sdf_val = cuboid.signed_distance(&(hit - wall_pos));
+                    if sdf_val.abs() < TOLERANCE {
+                        success += 1;
+                    } else {
+                        panic!("[Cuboid] Intersection inaccurate: SDF = {}", sdf_val)
+                    }
+                }
+                None => {
+                    failed += 1;
+                }
+            }
+        }
+
+        let failure_rate = (failed as f64 / TEST_RUNS as f64) * 100.0;
+
+        println!(
+            "[Cuboid Test] Success: {}, Failed: {}, Failure Rate: {:.2}%",
+            success,
+            failed,
+            failure_rate,
+        );
+
+        if failure_rate > 1.0 {
+            panic!("Failure rate too high! More than 1%")
+        }
+    }
+
+    #[test]
+    fn test_intersect_sphere() {
+        let mut success = 0;
+        let mut failed = 0;
+        for _ in 0..TEST_RUNS {
+            let mut rng = rand::rng();
+            let radius = rng.random_range(0.0..10.0);
+            let sphere = Sphere { radius };
+            let wall_pos = rand_vec(-5.0, 5.0);
+
+            let radius_in = rng.random_range(0.0..radius);
+            let inside_local = rand_vec(-1.0,1.0).normalize() * radius_in;
+
+            let velocity = rand_vec(-1.0, 1.0).normalize();
+            let outside_local = inside_local + velocity*radius*2.0;
+
+            let inside = wall_pos + inside_local;
+            let outside = wall_pos + outside_local;
+
+            if sphere.contains(&wall_pos, &outside) || !sphere.contains(&wall_pos, &inside) {
+                panic!("Wrong initialization in test");
+            }
+
+            let result = sphere.calculate_intersect(&outside, &velocity, &wall_pos, MAX_DT, TOLERANCE, MAX_STEPS);
+
+            match result {
+                Some(hit) => {
+                    let sdf_val = sphere.signed_distance(&(hit - wall_pos));
+                    if sdf_val.abs() < TOLERANCE {
+                        success += 1;
+                    } else {
+                        panic!("[Sphere] Intersection inaccurate: SDF = {}", sdf_val)
+                    }
+                }
+                None => {
+                    failed += 1;
+                }
+            }
+        }
+
+        let failure_rate = (failed as f64 / TEST_RUNS as f64) * 100.0;
+
+        println!(
+            "[Sphere Test] Success: {}, Failed: {}, Failure Rate: {:.2}%",
+            success,
+            failed,
+            failure_rate,
+        );
+
+        if failure_rate > 1.0 {
+            panic!("Failure rate too high! More than 1%")
+        }
+
+    }
+
+    #[test]
+    fn test_intersect_cylinder() {
+        let mut success = 0;
+        let mut failed = 0;
+        for _ in 0..TEST_RUNS {
+            let mut rng = rand::rng();
+            let radius = rng.random_range(0.1..10.0);
+            let length = rng.random_range(1.0..10.0);
+
+            // Random direction vector
+            let random_dir = rand_vec(-1.0, 1.0);
+
+            let wall_pos = rand_vec(-5.0, 5.0);
+            let cylinder = Cylinder::new(radius, length, random_dir);
+
+            let inside_local = {
+                let radial_in = rng.random_range(0.0..radius);
+                let angle_in = rng.random_range(0.0..std::f64::consts::TAU);
+                let x_in = radial_in * angle_in.cos();
+                let y_in = radial_in * angle_in.sin();
+                let z_in = rng.random_range(-length * 0.5..length * 0.5);
+                let inside = cylinder.perp_x * x_in + cylinder.perp_y * y_in + cylinder.direction * z_in;
+                inside
+            };
+
+            let velocity = rand_vec(-1.0, 1.0).normalize();
+            let outside_local = inside_local + velocity*22.5;
+
+            let inside = wall_pos + inside_local;
+            let outside = wall_pos + outside_local;
+
+            if cylinder.contains(&wall_pos, &outside) || !cylinder.contains(&wall_pos, &inside) {
+                panic!("Wrong initialization in test");
+            }
+
+            let result = cylinder.calculate_intersect(&outside, &velocity, &wall_pos, MAX_DT, TOLERANCE, MAX_STEPS);
+
+            match result {
+                Some(hit) => {
+                    let sdf_val = cylinder.signed_distance(&(hit - wall_pos));
+                    if sdf_val.abs() < TOLERANCE {
+                        success += 1;
+                    } else {
+                        panic!("[Cylinder] Intersection inaccurate: SDF = {}", sdf_val)
+                    }
+                }
+                None => {
+                    failed += 1;
+                }
+            }
+        }
+
+        let failure_rate = (failed as f64 / TEST_RUNS as f64) * 100.0;
+
+        println!(
+            "[Cylinder Test] Success: {}, Failed: {}, Failure Rate: {:.2}%",
+            success,
+            failed,
+            failure_rate,
+        );
+
+        if failure_rate > 1.0 {
+            panic!("Failure rate too high! More than 1%")
+        }
+    }
+
+    #[test]
+    fn test_normal_sphere() {
+        for _ in 0..TEST_RUNS {
+            let mut rng = rand::rng();
+            let radius = rng.random_range(0.1..10.0);
+            let sphere = Sphere { radius };
+            let wall_pos = rand_vec(-5.0, 5.0);
+            
+            let inside_local = rand_vec(-1.0, 1.0).normalize() * radius;
+            let inside = wall_pos + inside_local;
+
+            let result = sphere.calculate_normal(&inside, &wall_pos, TOLERANCE);
+            
+            match result {
+                Some(normal) => {
+                    assert_approx_eq!(-inside_local.normalize()[0], normal[0], 1e-12);
+                    assert_approx_eq!(-inside_local.normalize()[1], normal[1], 1e-12);
+                    assert_approx_eq!(-inside_local.normalize()[2], normal[2], 1e-12);
+                }
+                None => panic!("Normal calculation failed for sphere"),
+            }
+        }
+    }
+
+    #[test]
+    fn test_normal_cuboid() {
+        for _ in 0..TEST_RUNS {
+            let mut rng = rand::rng();
+            let length = rng.random_range(0.0..10.0);
+            let breadth = rng.random_range(0.0..10.0);
+            let height = rng.random_range(0.0..10.0);
+            let half_width = Vector3::new(length, breadth, height);
+
+            let cuboid = Cuboid { half_width };
+            let wall_pos = rand_vec(-5.0, 5.0);
+        
+            let face = rng.random_range(0..6); 
+        
+            let (inside_local, expected_normal) = match face {
+                0 => (Vector3::new(length, rng.random_range(-breadth.. breadth), rng.random_range(-height..height)), Vector3::new(-1.0,0.0,0.0)),// Positive X face
+                1 => (Vector3::new(-length, rng.random_range(-breadth.. breadth), rng.random_range(-height..height)),Vector3::new(1.0,0.0,0.0)), // Negative X face
+                2 => (Vector3::new(rng.random_range(-length..length), breadth, rng.random_range(-height..height)), Vector3::new(0.0,-1.0,0.0)),// Positive Y face
+                3 => (Vector3::new(rng.random_range(-length..length), -breadth, rng.random_range(-height..height)),Vector3::new(0.0,1.0,0.0)), // Negative Y face
+                4 => (Vector3::new(rng.random_range(-length..length), rng.random_range(-breadth.. breadth), height), Vector3::new(0.0,0.0,-1.0)),// Positive Z face
+                5 => (Vector3::new(rng.random_range(-length..length), rng.random_range(-breadth.. breadth), -height),Vector3::new(0.0,0.0,1.0)), // Negative Z face
+                _ => panic!("Invalid face selection"), // This should never happen
+            };
+        
+            let inside = wall_pos + inside_local;
+        
+            let result = cuboid.calculate_normal(&inside, &wall_pos, TOLERANCE);
+        
+            match result {
+                Some(normal) => {
+                    assert_approx_eq!(normal[0], expected_normal[0], 1e-12);
+                    assert_approx_eq!(normal[1], expected_normal[1], 1e-12);
+                    assert_approx_eq!(normal[2], expected_normal[2], 1e-12);
+                }
+                None => panic!("Normal calculation failed for cuboid"),
+            }
+        }
+    }
+    
+    // This is more complicated than the thing its testing
+    #[test]
+    fn test_normal_cylinder() {
+        for _ in 0..TEST_RUNS {
+            let mut rng = rand::rng();
+            let radius = 1.0;
+            let length = 2.0;
+            let random_dir = Vector3::new(0.0,0.0,1.0); // Random direction
+            let cylinder = Cylinder::new(radius, length, random_dir);
+            let wall_pos = rand_vec(-5.0, 5.0);
+    
+            let face = rng.random_range(0..3);
+    
+            let (surface_local, expected_normal) = match face {
+                0 => {  // Curved face
+                    let angle = rng.random_range(0.0..2.0 * PI);
+                    let x = radius * angle.cos();
+                    let y = radius * angle.sin();
+                    let z = rng.random_range(-length*0.5..length*0.5);
+                    let surface_local = x * cylinder.perp_x + y * cylinder.perp_y + z * cylinder.direction;
+                    let normal = -(x * cylinder.perp_x + y * cylinder.perp_y).normalize();
+                    (surface_local, normal) 
+                },
+                1 => {  // Top cap
+                    let radial = rng.random_range(0.0..radius);
+                    let angle = rng.random_range(0.0..2.0 * PI);
+                    let x = radial * angle.cos();
+                    let y = radial * angle.sin();
+                    let z = length * 0.5;
+                    let surface_local = x*cylinder.perp_x + y*cylinder.perp_y + z*cylinder.direction;
+                    (surface_local, -random_dir)
+                },
+                2 => {  // Bottom cap 
+                    let radial = rng.random_range(0.0..radius);
+                    let angle = rng.random_range(0.0..2.0 * PI);
+                    let x = radial * angle.cos();
+                    let y = radial * angle.sin();
+                    let z = -length * 0.5; 
+                    let surface_local = x*cylinder.perp_x + y*cylinder.perp_y + z*cylinder.direction;
+                    (surface_local, random_dir) 
+                },
+                _ => panic!("Invalid face selection"),
+            };
+            
+            if cylinder.signed_distance(&surface_local).abs() > TOLERANCE {
+                println!("{}", cylinder.signed_distance(&surface_local));
+                panic!("Surface point not initialized correctly fix test")
+            }
+
+            let surface = wall_pos + surface_local;
+    
+            let result = cylinder.calculate_normal(&surface, &wall_pos, TOLERANCE);
+    
+            match result {
+                Some(normal) => {
+                    assert_approx_eq!(normal[0], expected_normal[0], 1e-12);
+                    assert_approx_eq!(normal[1], expected_normal[1], 1e-12);
+                    assert_approx_eq!(normal[2], expected_normal[2], 1e-12);
+                }
+                None => panic!("Normal calculation failed for cylinder"),
+            }
+        }
+    }
+}
\ No newline at end of file
diff --git a/src/collisions/wall_collisions.rs b/src/collisions/wall_collisions.rs
new file mode 100644
index 00000000..cd885e91
--- /dev/null
+++ b/src/collisions/wall_collisions.rs
@@ -0,0 +1,621 @@
+/// Wall collision logic. Does not work with forces/accelertions, or multiple entities.
+
+use bevy::prelude::*;
+use nalgebra::Vector3;
+use rand::Rng;
+use rand::distr::Distribution;
+use crate::shapes::{
+    Cylinder as MyCylinder,
+    Cuboid as MyCuboid,
+    Sphere as MySphere,
+    CylindricalPipe,
+    Volume
+}; // Aliasing issues with bevy.
+use crate::constant::{PI, EXP};
+use crate::collisions::NumberOfWallCollisions;
+use crate::atom::{Atom, Position, Velocity};
+use crate::integrator::{Timestep, AtomECSBatchStrategy,};
+use crate::initiate::NewlyCreated;
+use crate::collisions::wall_collision_info::*;
+use super::LambertianProbabilityDistribution;
+
+#[derive(Resource)]
+pub struct MaxSteps(pub i32); // Max steps for sdf convergence
+
+#[derive(Resource)]
+pub struct SurfaceThreshold(pub f64); // Tolerance for convergence
+
+/// Struct to signify shape is a wall
+#[derive(Component)]
+#[component(storage = "SparseSet")]
+pub struct WallData{
+    pub wall_type: WallType,
+}
+
+/// Enum for designating wall type for collisions
+pub enum WallType {
+    // Specular
+    Smooth,
+    // Diffuse
+    Rough,
+    // choose randomly between specular and diffuse. Probability of specular reflection.
+    Random {spec_prob: f64},
+    // choose randomly depending on a coefficient that is exponential in angle of incidence.
+    Exponential{value: f64},
+}
+
+#[derive(Component)]
+pub struct DistanceToTravel {
+    pub distance_to_travel: f64
+}
+
+#[derive(Component)]
+pub enum VolumeStatus{
+    Inside,
+    Outside,
+    Open,
+}
+
+pub struct CollisionInfo{
+    pub atom: Entity,
+    pub wall: Entity,
+    pub collision_point: Vector3<f64>,
+    pub collision_normal: Vector3<f64>
+}
+
+pub trait Wall {
+    /// Collision check for open walls
+    fn preliminary_collision_check(&self, 
+        atom_pos: &Vector3<f64>, 
+        atom_vel: &Vector3<f64>, 
+        atom_location: &VolumeStatus,
+        wall_pos: &Vector3<f64>, 
+        dt: f64) -> bool;
+}
+
+impl<T : Volume> Wall for T {
+    fn preliminary_collision_check(&self, 
+        atom_pos: &Vector3<f64>, 
+        _atom_vel: &Vector3<f64>, 
+        atom_location: &VolumeStatus,
+        wall_pos: &Vector3<f64>, 
+        _dt: f64) -> bool {
+        match atom_location {
+            VolumeStatus::Inside => {
+                return !self.contains(wall_pos, atom_pos);
+            }
+            VolumeStatus::Outside => {
+                return self.contains(wall_pos, atom_pos);
+            }
+            VolumeStatus::Open => {
+                return false;
+            }
+        }
+    }
+}
+
+impl Wall for CylindricalPipe {
+    fn preliminary_collision_check(&self, 
+        atom_pos: &Vector3<f64>, 
+        atom_vel: &Vector3<f64>, 
+        _atom_location: &VolumeStatus,
+        wall_pos: &Vector3<f64>, 
+        dt: f64) -> bool {
+        let delta_prev = (atom_pos - atom_vel * dt) - wall_pos;
+        let axial_prev = delta_prev.dot(&self.direction);
+        let radial_prev = delta_prev - axial_prev * self.direction;
+
+        let delta_curr = atom_pos - wall_pos;
+        let axial_curr = delta_curr.dot(&self.direction);
+        let radial_curr = delta_curr - axial_curr * self.direction;
+
+        let is_within_radius_prev = radial_prev.norm() <= self.radius;
+        let is_within_length_prev = f64::abs(axial_prev) <= self.length / 2.0;
+        let is_within_radius_curr = radial_curr.norm() <= self.radius;
+        // let is_within_length_curr = f64::abs(axial_curr) <= self.length / 2.0;
+
+        if is_within_radius_prev ^ is_within_radius_curr {
+            if !is_within_length_prev {return false;}
+            else {return true;}
+        } else {
+            return false;
+        }
+    }
+}
+
+/// Checks which atoms have collided.
+fn collision_check<T: Wall + Intersect + Normal>(
+    atom: (Entity, &Position, &Velocity, &VolumeStatus),
+    wall: (Entity, &T, &Position),       
+    dt: f64,
+    tolerance: f64, 
+    max_steps: i32,
+) -> Option<CollisionInfo> {
+    let (atom_entity, atom_pos, atom_vel, atom_location) = atom;
+    let (wall_entity, shape, wall_pos) = wall;
+
+    if !shape.preliminary_collision_check(&atom_pos.pos, &atom_vel.vel, atom_location, &wall_pos.pos, dt) {
+        return None;
+    }
+
+    // Do collision check
+    if let Some(mut collision_point) = shape.calculate_intersect(&atom_pos.pos, &atom_vel.vel, &wall_pos.pos, dt, tolerance, max_steps) {
+        if let Some(collision_normal) = shape.calculate_normal(&collision_point, &wall_pos.pos, tolerance) {
+            collision_point += 1e-10 * collision_normal; // Offset collision point along normal to avoid numerical issues
+            return Some(CollisionInfo {
+                atom: atom_entity,
+                wall: wall_entity,
+                collision_point,
+                collision_normal,
+            });
+        } else {
+            eprintln!(
+                "Warning: Could not calculate normal. Atom: {:?}, Wall Pos: {:?}",
+                atom_pos.pos, wall_pos.pos
+            );
+            None
+        }
+    } else {
+        None
+    }
+}
+
+/// Specular reflection
+fn specular(
+    collision_normal: &Vector3<f64>,
+    velocity: &Vector3<f64>,
+) -> Vector3<f64> {
+    let reflected = velocity - 2.0 * velocity.dot(collision_normal) * collision_normal;
+    reflected
+}
+
+/// Diffuse collision (Lambertian)
+fn diffuse( 
+    collision_normal: &Vector3<f64>, 
+    velocity: &Vector3<f64>,
+    distribution: &LambertianProbabilityDistribution) -> Vector3<f64> {
+    // Get random weighted angles
+    let mut rng = rand::rng();
+    let phi = rng.random_range(0.0..2.0 * PI);
+    let theta = distribution.sample(&mut rng);
+
+    // Get basis vectors
+    let dir = Vector3::new(65.514, 5.54, 41.5).normalize();
+    let perp_x = collision_normal.cross(&dir).normalize();
+    let perp_y = collision_normal.cross(&perp_x).normalize();
+    
+    // Get scattered velocity
+    let x = phi.cos() * theta.sin();
+    let y = phi.sin() * theta.sin();
+    let z = theta.cos();
+
+    (x * perp_x + y * perp_y + z * collision_normal) * velocity.norm()
+}
+
+// To be run once at startup
+pub fn create_cosine_distribution(mut commands: Commands) {
+    // tuple list of (theta, weight)
+    let mut thetas = Vec::<f64>::new();
+    let mut weights = Vec::<f64>::new();
+
+    // precalculate the discretized cosine distribution.
+    let n = 1000; // resolution over which to discretize `theta`.
+    for i in 0..n {
+        let theta = (i as f64) / (n as f64) * PI / 2.0;
+        let weight = theta.cos() * theta.sin();
+        thetas.push(theta);
+        weights.push(weight);
+        // Note: we can exclude d_theta because it is constant and the distribution will be normalized.
+    }
+    let cosine_distribution = LambertianProbabilityDistribution::new(thetas, weights);
+    commands.insert_resource::<LambertianProbabilityDistribution>(cosine_distribution);
+    println!("Cosine distribution created!")
+}
+
+/// Do wall collision
+fn do_wall_collision(
+    atom: (&mut Position, &mut Velocity, &mut DistanceToTravel, &mut NumberOfWallCollisions),
+    wall: &WallData,
+    collision: &CollisionInfo,
+    distribution: &LambertianProbabilityDistribution,
+    dt: f64,
+    tolerance: f64,
+) -> f64 {
+    let (atom_pos, atom_vel, distance, num_of_collisions) = atom;
+    num_of_collisions.value += 1;
+    // subtract distance traveled to the collision point from prev position
+    let traveled = ((atom_pos.pos - atom_vel.vel*dt) - collision.collision_point).norm();
+    distance.distance_to_travel -= traveled;
+    if - distance.distance_to_travel > tolerance {
+        distance.distance_to_travel = 0.0;
+        eprintln!("Distance to travel set to a negative value somehow");
+    }
+    // println!("original pos {}", (atom_pos.pos - atom_vel.vel*dt));
+
+    // do collision
+    match wall.wall_type {
+        WallType::Smooth =>
+            {atom_vel.vel = specular(&collision.collision_normal, &atom_vel.vel)}
+        WallType::Rough =>
+            {atom_vel.vel = diffuse(&collision.collision_normal, &atom_vel.vel, &distribution)}
+        WallType::Random{spec_prob} => {
+            if rand::rng().random_bool(spec_prob) {
+                atom_vel.vel = specular(&collision.collision_normal, &atom_vel.vel)
+            } else {
+                atom_vel.vel = diffuse(&collision.collision_normal, &atom_vel.vel, &distribution)
+            }
+        }
+        WallType::Exponential{value} => {
+            let angle_of_incidence = (atom_vel.vel.dot(&collision.collision_normal)/atom_vel.vel.norm()).acos();
+            let prob_spec = (EXP.powf(value * angle_of_incidence) - 1.0)/(EXP.powf(value * PI/2.0) - 1.0);
+            if rand::rng().random_bool(prob_spec) {
+                atom_vel.vel = specular(&collision.collision_normal, &atom_vel.vel)
+            } else {
+                atom_vel.vel = diffuse(&collision.collision_normal, &atom_vel.vel, &distribution)
+            }
+        } 
+    }
+
+    if distance.distance_to_travel > 0.0 {
+        // propagate along chosen direction
+        atom_pos.pos = collision.collision_point + atom_vel.vel.normalize() * distance.distance_to_travel;
+    } else {
+        atom_pos.pos = collision.collision_point;
+    }
+    // Return time used
+    traveled / atom_vel.vel.norm()
+}
+
+/// Do the collisions
+/// To be run after verlet integrate position
+pub fn wall_collision_system<T: Wall + Component + Intersect + Normal>(
+    wall_query: Query<(Entity, &T, &WallData, &Position), Without<Atom>>,
+    mut atom_query: Query<(Entity, 
+        &mut Position, 
+        &mut Velocity, 
+        &mut DistanceToTravel, 
+        &mut NumberOfWallCollisions, 
+        &VolumeStatus), 
+        With<Atom>>,
+    batch_strategy: Res<AtomECSBatchStrategy>,
+    timestep: Res<Timestep>,
+    threshold: Res<SurfaceThreshold>, 
+    max_steps: Res<MaxSteps>,
+    distribution: Res<LambertianProbabilityDistribution>,
+) {
+    let tolerance = threshold.0;
+    let max_steps = max_steps.0;
+    let dt = timestep.delta;
+    let walls: Vec<_> = wall_query.iter().collect();
+    atom_query
+        .par_iter_mut()
+        .batching_strategy(batch_strategy.0.clone())
+        .for_each(|(atom_entity, mut pos, mut vel, mut distance, mut collisions, location)| {
+
+            let mut dt_atom = dt;
+            let mut num_of_collisions = 0;
+
+            while distance.distance_to_travel > 0.0 && num_of_collisions < 10000 { 
+                let mut collided = false;
+                for (wall_entity, shape, wall, wall_pos) in &walls {
+                    if let Some(collision_info) = collision_check(
+                        (atom_entity, &pos, &vel, &location),
+                        (*wall_entity, *shape, *wall_pos),
+                        dt_atom,
+                        tolerance,
+                        max_steps,
+                    ) {
+                        // Handle the collision
+                        let time_used = do_wall_collision((&mut pos, &mut vel, &mut distance, &mut collisions), 
+                            *wall,
+                            &collision_info, 
+                            &distribution, 
+                            dt_atom, 
+                            tolerance);
+                        collided = true;
+                        dt_atom -= time_used;
+                        if dt_atom < 0.0 {
+                            dt_atom = 0.0;
+                            // eprintln!("Time set to a negative value somehow, wall collision system");
+                        }
+                    }
+                }
+                
+                if !collided {
+                    break;
+                }
+                num_of_collisions += 1;
+            }
+        });
+}
+
+/// Updates distance to travel component for atoms
+/// If an atom doesn't have DistanceToTravel, it will be initialized
+/// Otherwise, it will be reset for the next frame
+pub fn update_distance_to_travel_system(
+    mut query_new: Query<(Entity, &Velocity), (With<Atom>, Without<DistanceToTravel>)>,
+    mut query_existing: Query<(&Velocity, &mut DistanceToTravel), With<Atom>>,
+    mut commands: Commands,
+    timestep: Res<Timestep>,
+) {
+    // Initialize for new atoms
+    for (atom_entity, vel) in query_new.iter_mut() {
+        commands.entity(atom_entity).insert(DistanceToTravel {
+            distance_to_travel: vel.vel.norm() * timestep.delta
+        });
+    }
+
+    // Reset for existing atoms
+    for (vel, mut distance) in query_existing.iter_mut() {
+        distance.distance_to_travel = vel.vel.norm() * timestep.delta;
+    }
+}
+
+/// Assign location status to atoms
+pub fn assign_location_status_system(
+    mut query: Query<(Entity, &Position), (With<Atom>, With<NewlyCreated>)>,
+    spheres: Query<(&MySphere, &Position), With<WallData>>,
+    cylinders: Query<(&MyCylinder, &Position), With<WallData>>,
+    cuboids: Query<(&MyCuboid, &Position), With<WallData>>,
+    mut commands: Commands,
+) {
+    for (atom_entity, pos) in query.iter_mut() {
+        let has_walls = spheres.iter().count() > 0 || 
+                        cylinders.iter().count() > 0 || 
+                        cuboids.iter().count() > 0;
+
+        let inside = spheres.iter().any(|(wall, wall_pos)| wall.contains(&wall_pos.pos, &pos.pos)) ||
+                     cylinders.iter().any(|(wall, wall_pos)| wall.contains(&wall_pos.pos, &pos.pos)) ||
+                     cuboids.iter().any(|(wall, wall_pos)| wall.contains(&wall_pos.pos, &pos.pos));
+
+        let status = if has_walls {
+            if inside {
+                VolumeStatus::Inside
+            } else {
+                VolumeStatus::Outside
+            }
+        } else {
+            // No closed walls, 
+            VolumeStatus::Open
+        };
+
+        commands.entity(atom_entity).insert(status);
+    }
+}
+
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::atom::{Atom, Position, Velocity};
+    use crate::shapes::{
+        Cylinder as MyCylinder,
+        Cuboid as MyCuboid,
+        Sphere as MySphere,
+    };
+    use assert_approx_eq::assert_approx_eq;
+
+    const TOLERANCE: f64 = 1e-9;
+    const MAX_STEPS: i32 = 1000;
+
+    #[test]
+    fn test_collision_check() {
+        let mut app = App::new();
+        
+        // Helper function to test collisions
+        fn test_collision<T>(
+            app: &mut App,
+            atom_pos: Vector3<f64>,
+            wall_pos: Vector3<f64>,
+            wall: T,
+            expected_point: Vector3<f64>,
+        ) where T: Component + Wall + Intersect + Normal {
+            let velocity = Velocity { vel: Vector3::new(1.0, 0.0, 0.0) };
+            let position = Position { pos: atom_pos };
+            let dt = 1.0 + 1e-15;
+
+            let atom = app.world_mut().spawn(position.clone())
+                .insert(velocity.clone())
+                .insert(Atom)
+                .insert(VolumeStatus::Inside)
+                .id();
+
+            let wall_position = Position { pos: wall_pos };
+            let wall_entity = app.world_mut().spawn(wall_position.clone())
+                .insert(WallData { wall_type: WallType::Smooth })
+                .insert(wall)
+                .id();
+
+            let wall_component = app.world().entity(wall_entity).get::<T>().unwrap();
+
+            let result = collision_check(
+                (atom, &position, &velocity, &VolumeStatus::Inside),
+                (wall_entity, wall_component, &wall_position),
+                dt, TOLERANCE, MAX_STEPS
+            );
+
+            match result {
+                Some(collision) => {
+                    assert_eq!(atom, collision.atom);
+                    assert_eq!(wall_entity, collision.wall);
+                    for i in 0..3 {
+                        assert_approx_eq!(expected_point[i], collision.collision_point[i], TOLERANCE);
+                        assert_approx_eq!([-1.0, 0.0, 0.0][i], collision.collision_normal[i], 1e-14);
+                    }
+                }
+                None => panic!("Collision not detected!"),
+            }
+
+            app.world_mut().despawn(atom);
+            app.world_mut().despawn(wall_entity);
+        }
+
+
+        // Test cases for different wall types
+        test_collision(
+            &mut app,
+            Vector3::new(5.0, 0.0, 0.0),
+            Vector3::new(0.0, 0.0, 0.0),
+            MySphere { radius: 4.0 },
+            Vector3::new(4.0, 0.0, 0.0),
+        );
+
+        test_collision(
+            &mut app,
+            Vector3::new(5.0, 0.0, 0.0),
+            Vector3::new(0.0, 0.0, 0.0),
+            MyCuboid { half_width: Vector3::new(4.0, 1.0, 1.0) },
+            Vector3::new(4.0, 0.0, 0.0),
+        );
+
+        test_collision(
+            &mut app,
+            Vector3::new(5.0, 0.0, 0.0),
+            Vector3::new(0.0, 0.0, 0.0),
+            MyCylinder::new(1.0, 8.0, Vector3::new(1.0, 0.0, 0.0)),
+            Vector3::new(4.0, 0.0, 0.0),
+        );
+
+        test_collision(
+            &mut app,
+            Vector3::new(2.5, 1.0, 1.0),
+            Vector3::new(1.0, 1.0, 1.0),
+            CylindricalPipe::new(1.0, 8.0, Vector3::new(0.0, 1.0, 0.0)),
+            Vector3::new(2.0, 1.0, 1.0),
+        );
+    }
+
+    // this is ugly as sin
+    #[test]
+    fn test_do_wall_collision() {
+        let mut app = App::new();
+        let dt = 1.0 - 1e-10;
+        let position = Position {pos: Vector3::new(1.0, 1.0, 0.0)};
+        let velocity = Velocity {vel: Vector3::new(2.0,2.0,0.0)};
+        let length = velocity.vel.norm() * dt;
+        let distance = DistanceToTravel{distance_to_travel: length};
+        let atom = app.world_mut()
+            .spawn(Atom)
+            .insert(position.clone())
+            .insert(velocity.clone())
+            .insert(DistanceToTravel{distance_to_travel: length})
+            .insert(NumberOfWallCollisions{value: 0})
+            .insert(VolumeStatus::Inside)
+            .id();
+        app.world_mut().spawn(WallData {wall_type: WallType::Smooth});
+        let collision_point = Vector3::new(0.0,0.0,0.0);
+
+        fn test_system(
+            mut query: Query<(Entity, &mut Position, &mut Velocity, &mut DistanceToTravel, &mut NumberOfWallCollisions)>,
+            wall: Query<(Entity, &WallData)>,
+            distribution: Res<LambertianProbabilityDistribution>,
+            tolerance: Res<SurfaceThreshold>) {
+            query.iter_mut().for_each(|(atom, mut atom_pos, mut atom_vel, mut distance, mut collisions)| {
+                for (wall_entity, wall) in wall.iter() {
+                    let dt = 1.0 - 1e-10; // If you change this make sure to change the dt above as well
+                    let collision_point = Vector3::new(0.0,0.0,0.0);
+                    let collision_normal = Vector3::new(-1.0,0.0,0.0);
+                    let collision_info = CollisionInfo {
+                        atom,
+                        wall: wall_entity,
+                        collision_point,
+                        collision_normal,
+                    };
+                    let time_used = do_wall_collision((&mut atom_pos, &mut atom_vel, &mut distance, &mut collisions), 
+                        wall,  
+                        &collision_info, 
+                        &distribution, 
+                        dt, 
+                        tolerance.0);
+                    assert_approx_eq!(time_used, ((atom_pos.pos - atom_vel.vel*dt) - collision_point).norm()/atom_vel.vel.norm())
+                }
+            });
+        }
+        app.add_systems(Startup, create_cosine_distribution);
+        app.add_systems(Update, test_system);
+        app.world_mut().insert_resource(SurfaceThreshold(1e-9));
+        app.update();
+
+        let new_velocity = app.world()
+            .entity(atom)
+            .get::<Velocity>()
+            .expect("Entity not found!")
+            .vel;
+        let new_position = app.world()
+            .entity(atom)
+            .get::<Position>()
+            .expect("Entity not found")
+            .pos;
+        let new_distance = app.world()
+            .entity(atom)
+            .get::<DistanceToTravel>()
+            .expect("Entity not found")
+            .distance_to_travel;
+        let new_number_of_wall_collisions = app.world()
+            .entity(atom)
+            .get::<NumberOfWallCollisions>()
+            .expect("Entity not found")
+            .value;
+            
+        let time_used = ((position.pos - velocity.vel*dt) - collision_point).norm()/velocity.vel.norm();
+        let distance_travelled = ((position.pos - velocity.vel*dt) - collision_point).norm();
+
+        assert_ne!(new_velocity, velocity.vel);
+        assert_eq!(new_number_of_wall_collisions, 1);
+        assert_approx_eq!(new_position[0], collision_point[0] + new_velocity[0]*(dt - time_used), 1e-15);
+        assert_approx_eq!(new_position[1], collision_point[1] + new_velocity[1]*(dt - time_used), 1e-15);
+        assert_approx_eq!(new_position[2], collision_point[2] + new_velocity[2]*(dt - time_used), 1e-15);
+        assert_approx_eq!(distance.distance_to_travel - distance_travelled, new_distance, 1e-15);
+    }
+
+    // Test multiple collisions
+    // Basically an integration test
+    #[test]
+    fn test_systems() {
+        use crate::integrator::{IntegrationPlugin, Timestep, OldForce};
+        use crate::initiate::NewlyCreated;
+        use crate::atom::{Force, Mass};
+        use crate::collisions::ApplyAtomCollisions;
+        use crate::collisions::CollisionPlugin;
+
+        let mut app = App::new();
+        app.world_mut()
+            .spawn(WallData{wall_type: WallType::Smooth})
+            .insert(CylindricalPipe::new(1.0, 8.0, Vector3::new(0.0, 0.0, 1.0)))
+            .insert(Position{pos: Vector3::new(0.0, 0.0, 0.0)});
+        let atom = app.world_mut()
+            .spawn(Atom)
+            .insert(Position{pos: Vector3::new(-1.0 + 1e-8, 0.0, 0.0)})
+            .insert(Velocity{vel: Vector3::new(0.0, -2.5 * PI, 0.0)})
+            .insert(NewlyCreated)
+            .insert(Force::default())
+            .insert(Mass { value: 1.0 })
+            .insert(OldForce(Force::default()))
+            .id();
+        app.add_plugins(IntegrationPlugin);
+        app.add_plugins(CollisionPlugin);
+        app.world_mut().insert_resource(ApplyAtomCollisions(false));
+        app.world_mut().insert_resource(Timestep {delta: 1.0});
+        app.update();
+
+        let new_position = app.world()
+            .entity(atom)
+            .get::<Position>()
+            .expect("Entity not found")
+            .pos;
+
+        let wall_collisions = app.world()
+            .entity(atom)
+            .get::<NumberOfWallCollisions>()
+            .expect("Entity not found")
+            .value;
+
+        println!("{}", new_position);
+        println!("wall collisions {}", wall_collisions); 
+        assert_approx_eq!(new_position[0], 0.0, 1e-5);
+        assert_approx_eq!(new_position[1], -1.0, 1e-5);
+        assert_approx_eq!(new_position[2], 0.0, 1e-5);
+    }
+}
+
+
diff --git a/src/destructor.rs b/src/destructor.rs
index 1e28fa70..dc9fe001 100644
--- a/src/destructor.rs
+++ b/src/destructor.rs
@@ -34,12 +34,9 @@ impl Plugin for DestroyAtomsPlugin {
 #[derive(Default, Component)]
 pub struct ToBeDestroyed;
 
-pub mod tests {
-    // These imports are actually needed! The compiler is getting confused and warning they are not.
-    #[allow(unused_imports)]
+#[cfg(test)]
+mod tests {
     use super::*;
-
-    #[allow(unused_imports)]
     use crate::atom::Position;
 
     #[test]
diff --git a/src/dipole/force.rs b/src/dipole/force.rs
index f78aacce..94ddacdd 100644
--- a/src/dipole/force.rs
+++ b/src/dipole/force.rs
@@ -1,6 +1,5 @@
 use crate::laser::intensity_gradient::LaserIntensityGradientSamplers;
 use bevy::prelude::*;
-extern crate nalgebra;
 use crate::atom::Force;
 use crate::dipole::DipoleLight;
 use crate::dipole::Polarizability;
@@ -30,14 +29,12 @@ pub fn apply_dipole_force_system<const N: usize>(
 
 
 #[cfg(test)]
-pub mod tests {
+mod tests {
     use super::*;
-
-    extern crate bevy;
     use assert_approx_eq::assert_approx_eq;
+    use nalgebra::Vector3;
+    // use bevy::prelude::*;
     use crate::integrator::AtomECSBatchStrategy;
-    use bevy::prelude::*;
-    extern crate nalgebra;
     use crate::constant;
     use crate::laser;
     use crate::laser::*;
@@ -45,9 +42,8 @@ pub mod tests {
     use crate::laser::gaussian::GaussianBeam;
     use crate::laser::DEFAULT_BEAM_LIMIT;
     use crate::laser::intensity_gradient::LaserIntensityGradientSampler;
-    use nalgebra::Vector3;
     use crate::atom::Force;
-    use crate::atom::Position;
+    // use crate::atom::Position;
 
     #[test]
     fn test_apply_dipole_force_system() {
diff --git a/src/dipole/mod.rs b/src/dipole/mod.rs
index f4e62ff3..5981fccf 100644
--- a/src/dipole/mod.rs
+++ b/src/dipole/mod.rs
@@ -4,7 +4,6 @@
 
 use crate::{constant};
 use crate::laser::index::LaserIndex;
-use crate::integrator::AtomECSBatchStrategy;
 use crate::dipole::force::apply_dipole_force_system;
 use crate::laser::LaserSystemsSet;
 
diff --git a/src/initiate.rs b/src/initiate.rs
index b76c3987..d23839e1 100644
--- a/src/initiate.rs
+++ b/src/initiate.rs
@@ -36,12 +36,12 @@ fn deflag_new_atoms(mut commands: Commands, query: Query<Entity, With<NewlyCreat
 pub struct InitiatePlugin;
 impl Plugin for InitiatePlugin {
     fn build(&self, app: &mut App) {
-        app.add_systems(Update, deflag_new_atoms.before(crate::atom_sources::emit::emit_once_system));
+        app.add_systems(Update, deflag_new_atoms.before(crate::atom_sources::emit::emit_number_per_frame_system));
     }
 }
 
-pub mod tests {
-    #[allow(unused_imports)]
+#[cfg(test)]
+mod tests {
     use super::*;
 
     /// Test the [NewlyCreated] component is properly removed from atoms after an update.
diff --git a/src/integration_tests/rate_equation.rs b/src/integration_tests/rate_equation.rs
index fc6707a6..45ced8da 100644
--- a/src/integration_tests/rate_equation.rs
+++ b/src/integration_tests/rate_equation.rs
@@ -3,7 +3,7 @@
 //! This module tests the rate equation implementation in atomecs by comparison to the exact analytic results for a single beam.
 
 #[cfg(test)]
-pub mod tests {
+mod tests {
     use crate::atom::{Atom, Force, Mass, Position, Velocity};
     use crate::initiate::NewlyCreated;
     use crate::integrator::Timestep;
@@ -17,9 +17,7 @@ pub mod tests {
     use crate::laser_cooling::{CoolingLight, LaserCoolingPlugin};
     use crate::simulation;
     use crate::species::Rubidium87_780D2;
-    extern crate nalgebra;
     use assert_approx_eq::assert_approx_eq;
-    use bevy::prelude::*;
     use nalgebra::Vector3;
 
     #[test]
diff --git a/src/integrator.rs b/src/integrator.rs
index a5cfa37b..c4fcef1b 100644
--- a/src/integrator.rs
+++ b/src/integrator.rs
@@ -101,7 +101,7 @@ fn clear_force(mut query: Query<&mut Force>, batch_strategy: Res<AtomECSBatchStr
 
 /// Stores the value of the force calculation from the previous frame.
 #[derive(Component, Default)]
-pub struct OldForce(Force);
+pub struct OldForce(pub Force);
 
 #[derive(Debug, Hash, PartialEq, Eq, Clone, SystemSet)]
 pub enum IntegrationSet {
@@ -135,8 +135,8 @@ impl Plugin for IntegrationPlugin {
     }
 }
 
-pub mod tests {
-    #[allow(unused_imports)]
+#[cfg(test)]
+mod tests {
     use super::*;
 
     #[test]
@@ -192,8 +192,7 @@ pub mod tests {
         app.world_mut().insert_resource(Timestep { delta: dt });
 
         // run simulation loop 1_000 times.
-        
-        let n_steps = 50000;
+        let n_steps = 1000;
         for _i in 0..n_steps {
             app.update()
         }
diff --git a/src/laser/intensity.rs b/src/laser/intensity.rs
index 26234d3a..ef7cd2f5 100644
--- a/src/laser/intensity.rs
+++ b/src/laser/intensity.rs
@@ -20,7 +20,7 @@ pub struct LaserIntensitySampler {
 impl Default for LaserIntensitySampler {
     fn default() -> Self {
         LaserIntensitySampler {
-            /// Intensity in SI units of W/m^2
+            // Intensity in SI units of W/m^2
             intensity: f64::NAN,
         }
     }
diff --git a/src/laser/intensity_gradient.rs b/src/laser/intensity_gradient.rs
index a021221b..9ac22cdb 100644
--- a/src/laser/intensity_gradient.rs
+++ b/src/laser/intensity_gradient.rs
@@ -21,7 +21,7 @@ pub struct LaserIntensityGradientSampler {
 impl Default for LaserIntensityGradientSampler {
     fn default() -> Self {
         LaserIntensityGradientSampler {
-            /// Intensity in SI units of W/m^2
+            // Intensity in SI units of W/m^2
             gradient: Vector3::new(f64::NAN, f64::NAN, f64::NAN),
         }
     }
diff --git a/src/laser/mod.rs b/src/laser/mod.rs
index de24882a..e0e144c6 100644
--- a/src/laser/mod.rs
+++ b/src/laser/mod.rs
@@ -75,8 +75,8 @@ pub enum LaserSystemsSet {
     IndexLasers,
 }
 
-pub mod tests {
-    #[allow(unused_imports)]
+#[cfg(test)]
+mod tests {
     use super::*;
 
     /// Test samplers are added to [NewlyCreated] entities.
diff --git a/src/laser_cooling/doppler.rs b/src/laser_cooling/doppler.rs
index c26a4674..c53f9c10 100644
--- a/src/laser_cooling/doppler.rs
+++ b/src/laser_cooling/doppler.rs
@@ -19,7 +19,7 @@ pub struct DopplerShiftSampler {
 impl Default for DopplerShiftSampler {
     fn default() -> Self {
         DopplerShiftSampler {
-            /// Doppler shift with respect to laser beam, in SI units of rad/s.
+            // Doppler shift with respect to laser beam, in SI units of rad/s.
             doppler_shift: f64::NAN,
         }
     }
diff --git a/src/laser_cooling/photons_scattered.rs b/src/laser_cooling/photons_scattered.rs
index f26766b1..572df05c 100644
--- a/src/laser_cooling/photons_scattered.rs
+++ b/src/laser_cooling/photons_scattered.rs
@@ -35,7 +35,7 @@ where
 {
     fn default() -> Self {
         TotalPhotonsScattered {
-            /// Number of photons scattered from all beams
+            // Number of photons scattered from all beams
             total: f64::NAN,
             phantom: PhantomData,
         }
@@ -74,7 +74,7 @@ where
 {
     fn default() -> Self {
         ExpectedPhotonsScattered {
-            ///photons scattered by the atom from a specific beam
+            // photons scattered by the atom from a specific beam
             scattered: f64::NAN,
             phantom: PhantomData,
         }
@@ -163,7 +163,7 @@ where
 {
     fn default() -> Self {
         ActualPhotonsScattered {
-            ///  number of photons actually scattered by the atom from a specific beam
+            //  number of photons actually scattered by the atom from a specific beam
             scattered: 0.0,
             phantom: PhantomData,
         }
diff --git a/src/laser_cooling/rate.rs b/src/laser_cooling/rate.rs
index acc3d287..b584d63f 100644
--- a/src/laser_cooling/rate.rs
+++ b/src/laser_cooling/rate.rs
@@ -21,7 +21,7 @@ pub struct RateCoefficient<T>
 where
     T: TransitionComponent,
 {
-    /// rate coefficient in Hz
+    // rate coefficient in Hz
     pub rate: f64,
     phantom: PhantomData<T>,
 }
@@ -31,7 +31,7 @@ where
 {
     fn default() -> Self {
         RateCoefficient {
-            /// rate coefficient in Hz
+            // rate coefficient in Hz
             rate: f64::NAN,
             phantom: PhantomData,
         }
diff --git a/src/laser_cooling/twolevel.rs b/src/laser_cooling/twolevel.rs
index c4be6d10..67b00a04 100644
--- a/src/laser_cooling/twolevel.rs
+++ b/src/laser_cooling/twolevel.rs
@@ -15,11 +15,11 @@ pub struct TwoLevelPopulation<T>
 where
     T: TransitionComponent,
 {
-    /// steady-state population density of the ground state, a LASER_COUNT
-    ///  in [0,1]
+    // steady-state population density of the ground state, a LASER_COUNT
+    //  in [0,1]
     pub ground: f64,
-    /// steady-state population density of the excited state, a LASER_COUNT
-    ///  in [0,1]
+    // steady-state population density of the excited state, a LASER_COUNT
+    //  in [0,1]
     pub excited: f64,
     marker: PhantomData<T>,
 }
@@ -37,11 +37,7 @@ where
 {
     fn default() -> Self {
         TwoLevelPopulation {
-            /// steady-state population density of the ground state, a LASER_COUNT
-            ///  in [0,1]
             ground: f64::NAN,
-            /// steady-state population density of the excited state, a LASER_COUNT
-            ///  in [0,1]
             excited: f64::NAN,
             marker: PhantomData,
         }
diff --git a/src/laser_cooling/zeeman.rs b/src/laser_cooling/zeeman.rs
index e7ca2a60..027b453a 100644
--- a/src/laser_cooling/zeeman.rs
+++ b/src/laser_cooling/zeeman.rs
@@ -30,11 +30,8 @@ where
 {
     fn default() -> Self {
         ZeemanShiftSampler::<T> {
-            /// Zeemanshift for sigma plus transition in rad/s
             sigma_plus: f64::NAN,
-            /// Zeemanshift for sigma minus transition in rad/s
             sigma_minus: f64::NAN,
-            /// Zeemanshift for pi transition in rad/s
             sigma_pi: f64::NAN,
             phantom: PhantomData,
         }
diff --git a/src/lib.rs b/src/lib.rs
index 986f52e9..2e9e9ce1 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -9,11 +9,11 @@ extern crate atomecs_derive;
 
 pub mod atom;
 pub mod atom_sources;
+pub mod bevy_bridge;
 pub mod collisions;
 pub mod constant;
 pub mod destructor;
 pub mod dipole;
-pub mod bevy_bridge;
 pub mod gravity;
 pub mod initiate;
 pub mod integration_tests;
@@ -21,6 +21,7 @@ pub mod integrator;
 pub mod laser;
 pub mod laser_cooling;
 pub mod magnetic;
+pub mod marker;
 pub mod maths;
 pub mod output;
 pub mod ramp;
diff --git a/src/marker.rs b/src/marker.rs
new file mode 100644
index 00000000..79a49c8a
--- /dev/null
+++ b/src/marker.rs
@@ -0,0 +1,30 @@
+use bevy::prelude::*;
+use crate::atom::{Atom, Position, Velocity};
+use crate::collisions::CollisionsSet;
+
+#[derive(Component)]
+pub struct Marker;
+
+pub fn add_marker_system(
+    mut commands: Commands,
+    query: Query<(Entity, &Position, &Velocity), With<Atom>>,
+) {
+    for (entity, pos, vel) in query.iter() {
+        // if pos.pos.x < 210e-6 && pos.pos.x > 200e-6 {
+        //     commands.entity(entity).insert(Marker);
+        // }
+        // else if pos.pos.x > 210e-6 {
+        //     commands.entity(entity).remove::<Marker>();
+        // }
+        if pos.pos.x > 0.0 {
+            commands.entity(entity).insert(Marker);
+        }
+    }
+}
+
+pub struct MarkerPlugin;
+impl Plugin for MarkerPlugin {
+    fn build(&self, app: &mut App) {
+        app.add_systems(PreUpdate, add_marker_system.after(CollisionsSet::WallCollisionSystems));
+    }
+}
\ No newline at end of file
diff --git a/src/shapes.rs b/src/shapes.rs
index ffcbad2f..abcf76b8 100644
--- a/src/shapes.rs
+++ b/src/shapes.rs
@@ -36,8 +36,8 @@ pub struct Cylinder {
 impl Cylinder {
     pub fn new(radius: f64, length: f64, direction: Vector3<f64>) -> Cylinder {
         let dir = Vector3::new(0.23, 1.2, 0.4563).normalize();
-        let perp_x = direction.normalize().cross(&dir);
-        let perp_y = direction.normalize().cross(&perp_x);
+        let perp_x = direction.normalize().cross(&dir).normalize();
+        let perp_y = direction.normalize().cross(&perp_x).normalize();
         Cylinder {
             radius,
             length,
@@ -183,4 +183,49 @@ impl Surface for Cuboid {
         let position = surface_position + point;
         (position, normal)
     }
-}
\ No newline at end of file
+}
+
+/// A cylindrical Pipe, open at both ends.
+#[derive(Component)]
+#[component(storage = "SparseSet")]
+pub struct CylindricalPipe {
+    /// Radius of the cylindrical volume.
+    pub radius: f64,
+    /// Length of the cylindrical volume.
+    pub length: f64,
+    /// A normalised vector, aligned to the direction of the cylinder.
+    pub direction: Vector3<f64>,
+    /// A normalised vector, aligned perpendicular to the cylinder.
+    pub perp_x: Vector3<f64>,
+    /// A normalised vector, aligned perpendicular to the cylinder.
+    pub perp_y: Vector3<f64>,
+}
+
+impl CylindricalPipe {
+    pub fn new(radius: f64, length: f64, direction: Vector3<f64>) -> CylindricalPipe {
+        let dir = Vector3::new(0.23, 1.2, 0.4563).normalize();
+        let perp_x = direction.normalize().cross(&dir).normalize();
+        let perp_y = direction.normalize().cross(&perp_x).normalize();
+        CylindricalPipe {
+            radius,
+            length,
+            direction: direction.normalize(),
+            perp_x,
+            perp_y,
+        }
+    }
+}
+
+impl Surface for CylindricalPipe {
+    fn get_random_point_on_surface(
+        &self,
+        surface_position: &Vector3<f64>,
+    ) -> (Vector3<f64>, Vector3<f64>) {
+        let mut rng = rand::rng();
+        let angle = rng.random_range(0.0..2.0 * std::f64::consts::PI);
+        let axial = rng.random_range(-self.length..self.length) / 2.0;
+        let normal = self.perp_x * angle.cos() + self.perp_y * angle.sin();
+        let point = surface_position + normal * self.radius + self.direction * axial;
+        (point, normal)
+    }
+}
diff --git a/src/simulation.rs b/src/simulation.rs
index 9177e597..67650047 100644
--- a/src/simulation.rs
+++ b/src/simulation.rs
@@ -5,7 +5,7 @@ use bevy::{app::TaskPoolThreadAssignmentPolicy, log::LogPlugin, prelude::*};
 use crate::{
     destructor::DestroyAtomsPlugin, gravity::GravityPlugin, initiate::InitiatePlugin,
     integrator::IntegrationPlugin, magnetic::MagneticsPlugin,
-    output::console_output::console_output, sim_region::SimulationRegionPlugin,
+    output::console_output::console_output, sim_region::SimulationRegionPlugin, marker::MarkerPlugin,
 };
 
 /// Used to construct a simulation in AtomECS.
@@ -73,12 +73,13 @@ impl Default for SimulationBuilder {
             task_pool_options,
         });
 
+        builder.app.add_plugins(DestroyAtomsPlugin);
+        builder.app.add_plugins(GravityPlugin);
+        builder.app.add_plugins(InitiatePlugin);
         builder.app.add_plugins(IntegrationPlugin);
         builder.app.add_plugins(MagneticsPlugin);
+        builder.app.add_plugins(MarkerPlugin);
         builder.app.add_plugins(SimulationRegionPlugin);
-        builder.app.add_plugins(GravityPlugin);
-        builder.app.add_plugins(DestroyAtomsPlugin);
-        builder.app.add_plugins(InitiatePlugin);
         builder.app.add_systems(Update, console_output);
         builder
     }