diff --git a/.github/workflows/cargo.yml b/.github/workflows/cargo.yml index d643668..40da11b 100644 --- a/.github/workflows/cargo.yml +++ b/.github/workflows/cargo.yml @@ -18,42 +18,24 @@ jobs: uses: actions/checkout@v4 with: fetch-depth: 0 + - name: Setup Rust toolchain uses: actions-rust-lang/setup-rust-toolchain@v1 - - name: Restore cargo cache - uses: actions/cache@v3 + + - name: Populate target directory from cache + uses: Leafwing-Studios/cargo-cache@v2 with: - path: | - ~/.cargo/registry/index - ~/.cargo/registry/cache - target - key: ${{ runner.os }}-cargo-${{ hashFiles('**/Cargo.lock') }} - restore-keys: | - ${{ runner.os }}-cargo- + sweep-cache: true + + - name: Install Bevy dependencies + run: sudo apt-get update; sudo apt-get install --no-install-recommends libasound2-dev libudev-dev libwayland-dev libxkbcommon-dev + - name: Cargo Check run: RUSTFLAGS="$RUSTFLAGS -A dead_code" cargo check + - name: Cargo Build run: cargo build --release --all-features - test: - needs: build - runs-on: ubuntu-latest - steps: - - name: Checkout - uses: actions/checkout@v4 - with: - fetch-depth: 0 - - name: Setup Rust toolchain - uses: actions-rust-lang/setup-rust-toolchain@v1 - - name: Restore cargo cache - uses: actions/cache@v3 - with: - path: | - ~/.cargo/registry/index - ~/.cargo/registry/cache - target - key: ${{ runner.os }}-cargo-${{ hashFiles('**/Cargo.lock') }} - restore-keys: | - ${{ runner.os }}-cargo- + doc: runs-on: ubuntu-latest needs: @@ -63,24 +45,26 @@ jobs: uses: actions/checkout@v4 with: fetch-depth: 0 + - name: Setup Rust toolchain uses: actions-rust-lang/setup-rust-toolchain@v1 - - name: Restore cargo cache - uses: actions/cache@v3 + + - name: Install Bevy dependencies + run: sudo apt-get update; sudo apt-get install --no-install-recommends libasound2-dev libudev-dev libwayland-dev libxkbcommon-dev + + - name: Populate target directory from cache + uses: Leafwing-Studios/cargo-cache@v2 with: - path: | - ~/.cargo/registry/index - ~/.cargo/registry/cache - target - key: ${{ runner.os }}-cargo-${{ hashFiles('**/Cargo.lock') }} - restore-keys: | - ${{ runner.os }}-cargo- + sweep-cache: true + - name: Cargo Doc run: cargo doc --no-deps + - uses: actions/upload-artifact@v4 with: name: doc path: target/doc + deploy: if: github.ref == 'refs/heads/main' needs: doc @@ -97,14 +81,17 @@ jobs: - uses: actions/checkout@v3 with: fetch-depth: 0 + - uses: actions/download-artifact@v4 with: name: doc path: ${{ github.workspace }} + - name: Upload Pages artifact uses: actions/upload-pages-artifact@v3 with: path: target/doc + - name: Deploy to GitHub Pages id: deployment uses: actions/deploy-pages@v4 diff --git a/Cargo.toml b/Cargo.toml index b58b982..7434818 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -2,49 +2,101 @@ name = "yahs" description = "Yet Another HAB Simulator" authors = ["Philip Linden "] -version = "0.2.0" +version = "0.4.0" edition = "2021" readme = "README.md" license-file = "LICENSE" [features] -default = ["gui"] -gui = [ - "egui", - "egui_extras", - "egui_plot", - "eframe", - "emath", - "rfd", +default = [ + # Default to a native dev build. + "dev_native", ] +dev = [ + # Improve compile times for dev builds by linking Bevy as a dynamic + # library. + "bevy/dynamic_linking", + "bevy/bevy_dev_tools", +] +dev_native = [ + "dev", + # Enable system information plugin for native dev builds. + "bevy/sysinfo_plugin", + "iyes_perf_ui/sysinfo", +] +config-files = ["ron", "bevy_common_assets", "serde"] +inspect = ["bevy-inspector-egui", "bevy_panorbit_camera/bevy_egui"] [dependencies] -pretty_env_logger = "0.5.0" -libm = "0.2.1" -toml = "0.8.10" -clap = { version = "4.1", default-features = false, features = [ - "derive", - "std", - "help", - "usage", - "error-context", - "suggestions", -] } -csv = "1.2.1" -serde = { version = "1.0.196", features = ["derive"] } -log = { version = "0.4.20", features = ["release_max_level_debug"] } -egui = { version = "0.26.2", features = ["log", "serde"], optional = true } -egui_plot = { version = "0.26.2", features = ["serde"], optional = true } -eframe = { version = "0.26.2", features = ["persistence"], optional = true } -emath = { version = "0.26.2", optional = true } -egui_extras = { version = "0.26.2", features = [ - "chrono", - "datepicker", - "file", -], optional = true } -rfd = { version = "0.14.0", optional = true } -ultraviolet = { version = "0.9.2", features = ["serde"] } +# core dependencies +bevy = "0.14.2" +bevy-trait-query = "0.6.0" +# physics dependencies +avian3d = { version = "0.1.2", features = ["debug-plugin"] } +# ui dependencies +bevy_panorbit_camera = { version = "0.20.0" } +bevy-inspector-egui = { version = "0.27.0", features = ["highlight_changes"], optional = true } +iyes_perf_ui = "0.3.0" +# file io dependencies +bevy_common_assets = { version = "0.11.0", features = ["ron"], optional = true } +ron = { version = "0.8.1", optional = true } +serde = { version = "1.0", features = ["derive"], optional = true } +parry3d = { version = "0.17.2", features = ["parallel"] } [[bin]] name = "yahs" path = "src/main.rs" + +# ----------------------------------------------------------------------------- +# Some Bevy optimizations +# ----------------------------------------------------------------------------- + +# Idiomatic Bevy code often triggers these lints, and the CI workflow treats +# them as errors. In some cases they may still signal poor code quality however, +# so consider commenting out these lines. +[lints.clippy] +# Bevy supplies arguments to systems via dependency injection, so it's +# natural for systems to request more than 7 arguments -- which triggers +# this lint. +too_many_arguments = "allow" +# Queries that access many components may trigger this lint. +type_complexity = "allow" + + +# Compile with Performance Optimizations: +# https://bevyengine.org/learn/quick-start/getting-started/setup/#compile-with-performance-optimizations + +# Enable a small amount of optimization in the dev profile. +[profile.dev] +opt-level = 1 + +# Enable a large amount of optimization in the dev profile for dependencies. +[profile.dev.package."*"] +opt-level = 3 + +# The default profile is optimized for Wasm builds because that's what [Trunk +# reads](https://github.com/trunk-rs/trunk/issues/605). Optimize for size in the +# wasm-release profile to reduce load times and bandwidth usage on web. +[profile.release] +# Compile the entire crate as one unit. Slows compile times, marginal +# improvements. +codegen-units = 1 +# Do a second optimization pass over the entire program, including +# dependencies. Slows compile times, marginal improvements. +lto = "thin" +# Optimize with size in mind (also try "z", sometimes it is better). +# Slightly slows compile times, great improvements to file size and runtime +# performance. +opt-level = "s" +# Strip all debugging information from the binary to slightly reduce file +# size. +strip = "debuginfo" + +# Override some settings for native builds. +[profile.release-native] +# Default to release profile values. +inherits = "release" +# Optimize with performance in mind. +opt-level = 3 +# Keep debug information in the binary. +strip = "none" diff --git a/README.md b/README.md index 9786b1f..c62f4de 100644 --- a/README.md +++ b/README.md @@ -1,25 +1,12 @@ # yet another HAB simulator -A high altitude balloon flight simulator based on [tkschuler/EarthSHAB](https://github.com/tkschuler/EarthSHAB) and [brickworks/mfc-apps](https://github.com/Brickworks/mfc-apps), built on Rust. +A high altitude balloon flight simulator based on +[tkschuler/EarthSHAB](https://github.com/tkschuler/EarthSHAB) and +[brickworks/mfc-apps](https://github.com/Brickworks/mfc-apps), built on Rust. ## Installation + ```sh cargo install --git https://github.com/Brickworks/yahs.git --features gui yahs --help ``` - -## Usage (CLI) -Set the simulation configuration in `config/default.toml` - -Then use the `start` command to start a flight simulation. Use the `RUST_LOG` -environment variable to specify the log level to report. - -```sh -yahs start -``` - -## Usage (GUI) - -```sh -yahs gui -``` diff --git a/assets/configs/properties.ron b/assets/configs/properties.ron new file mode 100644 index 0000000..518cd3c --- /dev/null +++ b/assets/configs/properties.ron @@ -0,0 +1,103 @@ +( + gases: [ + // Species of gas with a known molar mass (kg/mol) + ( + name: "Air", + abbreviation: "Air", + molar_mass: 0.02897, + ), + ( + name: "Helium", + abbreviation: "He", + molar_mass: 0.0040026, + ), + ( + name: "Hydrogen", + abbreviation: "H2", + molar_mass: 0.00201594, + ), + ( + name: "Nitrogen", + abbreviation: "N2", + molar_mass: 0.0280134, + ), + ( + name: "Oxygen", + abbreviation: "O2", + molar_mass: 0.0319988, + ), + ( + name: "Argon", + abbreviation: "Ar", + molar_mass: 0.039948, + ), + ( + name: "Carbon Dioxide", + abbreviation: "CO2", + molar_mass: 0.04400995, + ), + ( + name: "Neon", + abbreviation: "Ne", + molar_mass: 0.020183, + ), + ( + name: "Krypton", + abbreviation: "Kr", + molar_mass: 0.08380, + ), + ( + name: "Xenon", + abbreviation: "Xe", + molar_mass: 0.13130, + ), + ( + name: "Methane", + abbreviation: "CH4", + molar_mass: 0.01604303, + ), + ], + // Properties, attributes and functions related to the material properties. + // Source: https://www.matweb.com/ + materials: [ + ( + name: "Nothing", + max_temperature: 1e30, + density: 0.0, + emissivity: 1.0, + absorptivity: 0.0, + thermal_conductivity: 1e30, + specific_heat: 0.0, + poissons_ratio: 0.5, + elasticity: 1e30, + max_strain: 1e30, + max_stress: 1e30, + ), + ( + name: "Rubber", + max_temperature: 385.0, + density: 1000.0, + emissivity: 0.86, + absorptivity: 0.86, + thermal_conductivity: 0.25, + specific_heat: 1490.0, + poissons_ratio: 0.5, + elasticity: 4000000.0, + max_strain: 8.0, + max_stress: 25000000.0, + ), + ( + name: "LowDensityPolyethylene", + max_temperature: 348.0, + density: 919.0, + emissivity: 0.94, + absorptivity: 0.94, + thermal_conductivity: 0.3175, + specific_heat: 2600.0, + poissons_ratio: 0.5, + elasticity: 300_000_000.0, + max_strain: 6.25, + max_stress: 10_000_000.0, + ), + ] +) diff --git a/assets/textures/balloon.png b/assets/textures/balloon.png new file mode 100644 index 0000000..029adff Binary files /dev/null and b/assets/textures/balloon.png differ diff --git a/assets/textures/splash.png b/assets/textures/splash.png new file mode 100644 index 0000000..d984dbb Binary files /dev/null and b/assets/textures/splash.png differ diff --git a/config/default.toml b/config/default.toml deleted file mode 100644 index 464dc57..0000000 --- a/config/default.toml +++ /dev/null @@ -1,29 +0,0 @@ -[environment] -real_time = false -tick_rate_hz = 10.0 -max_elapsed_time_s = 10_000.0 -initial_altitude_m = 0.0 -initial_velocity_m_s = 0.0 - -[balloon] -material = "Rubber" -thickness_m = 0.0001 # thicker balloons have less stress -barely_inflated_diameter_m = 1.5 # larger balloons have less strain - -[balloon.lift_gas] -species = "Helium" -mass_kg = 0.3 - -[bus.body] -mass_kg = 1.427 -drag_area_m2 = 0.04 -drag_coeff = 2.1 - -[bus.parachute] -total_mass_kg = 0.2 -drogue_area_m2 = 0.2 -drogue_drag_coeff = 1.5 -main_area_m2 = 1.1 -main_drag_coeff = 0.8 -deploy_force_n = 100 -deploy_time_s = 5.0 diff --git a/docs/Drag-Force-in-External-Flows-Lesson-2-Handout.pdf b/docs/Drag-Force-in-External-Flows-Lesson-2-Handout.pdf new file mode 100644 index 0000000..67efca4 Binary files /dev/null and b/docs/Drag-Force-in-External-Flows-Lesson-2-Handout.pdf differ diff --git a/docs/Fluids-lecture.pdf b/docs/Fluids-lecture.pdf new file mode 100644 index 0000000..3080575 Binary files /dev/null and b/docs/Fluids-lecture.pdf differ diff --git a/docs/Newtonian_Aerodynamics_for_General_Body_Shapes.pdf b/docs/Newtonian_Aerodynamics_for_General_Body_Shapes.pdf new file mode 100644 index 0000000..a1517a4 Binary files /dev/null and b/docs/Newtonian_Aerodynamics_for_General_Body_Shapes.pdf differ diff --git a/docs/devlog.md b/docs/devlog.md new file mode 100644 index 0000000..4824b0d --- /dev/null +++ b/docs/devlog.md @@ -0,0 +1,330 @@ +# development log + +## 2024-11-18 again + +I think I was a bit naive to install `bevy-trait-query`. It works for now but in +the future we should really move away from it. It is currently a crutch. + +For some reason when the balloon bounces, it accelerates upward to oblivion. I +added some errors to stop the simulation when the balloon goes out of bounds, +but right now it panics. Not great but better than nothing. Best I can do is +bring in the "out of bounds" level to be less than the true bounds. + +The drag force is the one that is causing the bad acceleration. Weight and +buoyancy don't change with time since the balloon has constant size right now. +Annoying since the drag is supposed to be the placeholder. + +I found the bug. The drag's "flow velocity" variable was being set to the +_balloon_ velocity, rather than the opposing force. After making that change, +the balloon gently floats to the ground when filled with air, and gently rises +when filled with helium. + +- Added pause/play controls. Default key is `Space`. +- Added a new ui that displays the simulation state. +- Added a new ui that displays forces in real time, assuming only one balloon. +- Added an `Anomaly` state to the sim that is supposed to freeze. + +## 2024-11-18 + +Iterating with AI on the drag calculations. It is set up to perform a raycast +along the velocity vector and compute the drag as a function of the angle of +attack. Each raycast hits a point on the surface of the collider and the normal +and differential area is used to compute the drag force for each one. + +Turns out air viscosity also changes with altitude for the standard atmosphere. +I learned today that there is also _dynamic_ viscosity and _kinematic_ +viscosity. The latter is the former divided by the density. So as density +changes, so too does viscosity. These two types are really just different ways +to express the same thing. + +In any case, I will wire up drag to be calculated from the bounding-sphere of +the collider. That way I can get a simple drag calculation working before +wiring up the more complex stuff and ignore things like shear and asymmetry for +now. + +I reorganized the module hierarchy again shut up I'm not the problem you are. + +- Split the `forces` module into `body` and `aero`. The base `forces` module + contains the common code for all forces. +- Added `AeroPlugin` for computing drag on solid bodies. +- Added `BodyForcesPlugin` that deals with forces related to the mass and + volume of the rigid bodies. +- Added `ForceVisualizationPlugin` for visualizing the forces as vectors in the + scene. +- Moved serde features behind the `config-files` feature flag. I'm hoping to + pare down the default dependencies for the project. +- Added `intro_to_drag.md` as I learn how to compute drag on arbitrary shapes. +- Added `camera.rs` as a simple pan-orbit camera controller. +- Added `controls.rs` for managing keybindings, and `KeyBindingsConfig` resource + for storing the keybindings. I use a resource instead of constants because + it makes the keybindings easier to change at runtime and feels more natural + alongside how Bevy handles inputs anyway. +- Moved `scene` and `ui` modules to a `graphics` module. Hopefully this will + make it easier to separate the concerns of physics simulation and graphics, + like if I ever want to add a TUI or CLI instead of a 3D graphics UI. + +## 2024-11-17 + +Time to see what's going on with the forces through some new UIs instead of +relying on the inspector plugin. It would be nice to have some gizmos that show +the forces as vectors in the scene. + +I added a new plugin that lets me query components by a trait. This greatly +simplifies the logic for collecting forces, and makes the coupling between the +force components and the equations that define them more direct. + +I used the [Newtonian Aerodynamics for General Body Shapes](https://ntrs.nasa.gov/citations/19780014285) +paper to help me understand the equations for computing drag. o1-mini helped +interpret the equations and implement them. A summary of the paper is saved in +[drag.md](drag.md). I overcomplicated the drag calculation by sampling points on +the surface of the collider and then projecting them onto the surface but it was +fun to implement. + +## 2024-11-16 + +I figured out how to get the forces to update correctly setting `WeightForce`, +`BuoyantForce`, and `DragForce` to wrap an `ExternalForce` component was a mistake. +Instead, I had to change them to wrap a `Vec3`, which is the input to the +`set_force` method. What I'm doing now is having systems that update the forces, +then after that run the system that adds the forces into a single vector. That +is what gets applied to the rigid body. + +Remember that `WeightForce`, `BuoyantForce`, and `DragForce` are _components_, +not forces themselves. We name them so we can query the _kinds_ of forces, +otherwise generic force vectors would not be distinguishable from each other. +The alternative (when we can't tell which vector is which) is to add the vectors +to a "blind" array and then update the `ExternalForce` component as before, +except that when we sum them we don't know which vector is which. Now that I +think about it, that might be easier and simpler in the long run. If we need to +know each force's value and identity, we can add it to a resource or something +that the UI can query. + +Update: it was not better. + +I fell in love with [iyes_perf_ui](https://github.com/IyesGames/iyes_perf_ui). +It is simple and easy to work with. Egui was adding a lot of overhead, more than +I want right now. This is perfect until I'm done with the basic features of the +simulation. + +- Fixed systems for computing the weight, buoyancy, and drag forces. +- Dropped egui. +- Added `iyes_perf_ui` for performance monitoring and other UI. + +## 2024-11-14 + +My focus today is to spend some time getting the balloon to "look right". Things +like the balloon's shape, skin material, and thickness. Also things like the +tether and a payload box. My plan for balloon material stress is to calculate +stress from the change in surface area. The gas volume can change shape all it +wants, as long as the density follows the ideal gas law. The balloon skin will +stretch and deform to follow the gas volume shape, and the stress will be +calculated from the change in surface area of the skin's mesh. + +My goal at the end of today is to see a balloon with a little payload hanging +from it that slowly rises into the air. + +In Avian, forces persist across frames by default. The way the forces are set up +now ends up applying the same force every frame, which is accumulating and not +correct. I think there are also some issues with applying forces in separate +systems. I'm not entirely sure what the best way to handle this is. I think if +we remove the persistence assumption then every frame the only forces applied +are the ones that _we_ explicitly apply in systems. This situation is more +complex than the simple case Avian is pre-configured for, like 1G of gravity +down and maybe some friction. + +I'm second-guessing whether I should use a new `ExternalForce` component for +each force instead of mutating the default component that spawns in when the +rigid body is created. The default component is persistent by default, and +this is probably why the forces are accumulating. It feels gross to have to +clear the force every frame or apply an impulse every frame instead of just +updating the force. + +- Added a bunch of handy geometry functions to `properties.rs`. +- Moved the `PbrBundle` into the `BalloonBundle` so the balloon is more tightly + associated with its mesh, material, and transform. Here the PBR (physics based + rendering) is not relevant to the physics simulation, it just happens to be + the bundle that contains these components. +- Added `PayloadPlugin` for spawning a payload and tether. Placeholder for now. +- Removed `#![allow(dead_code)]` from a few modules as I remove placeholders. + +## 2024-11-13 + +Avian has its own `Mass` component but it is not as convenient for this project +as something custom that can be updated with the other thermodynamics properties +like temperature and pressure. The Avian component is not reflected in the +inspector plugin so it's hard to debug, too. I can make a custom `Mass` that +hooks up to the Avian one and then otherwise is whatever I want. + +I realized I was using `ExternalForce` incorrectly. Instead of adding more +force components, I should be applying additional vectors to the existing +`ExternalForce` component using the `apply_force` method. + +- Added names to the all entities in the test scene. +- Moved core properties (Temperature, Pressure, Volume, Density, Mass) to their + own module called `properties.rs` and added a plugin for registering them. To + make things more ergonomic, these properties are exposed at the top level of + the simulator plugin so that they can be used in other systems. +- Removed the `thermodynamics` module since it became redundant. +- Added a `CorePhysicsPlugin` for registering the physics plugins. +- Added a `ScenePlugin` for setting up the scene with the HAB and atmosphere. +- Added systems to update the density, weight, and buoyant forces every frame. + This wasn't working until I added the default Avian physics plugins, + obviously. +- Got buoyancy and weight forces working on the collider. When the gas species + is set to air, the balloon slowly sinks to the ground. When the gas species + is set to helium, the balloon floats. + +## 2024-11-11 + +More tightening up of the systems using `bevy-inspector-egui` and making sure +the numbers update correctly. + +I am considering reducing the `IdealGas` component to be paired with a rigid +body. The rigid body would have a `Mass` and a `Collider` (sphere, ellipsoid, +etc.) and the `IdealGas` would have a `Volume` and a `GasSpecies`. This would +allow for more flexibility in how the gas volume is computed and updated. For +example, if the gas volume is not attached to the rigid body it may be +desirable to compute the volume based on the mesh instead. It would also allow +for more complex gas interactions with the rigid body, like inside a parachute. + +Converted forces to use the new `ExternalForce` component from Avian. + +- Renamed `mechanics` module to `forces` +- Corrected the buoyancy force calculation to use only the ambient density of + the atmosphere and not the density of the gas inside the balloon. The weight + of the gas is already accounted for in the `mass` of the gas. +- Added the `Component` trait to `Temperature`, `Pressure`, `Density`, and + `Volume`. +- Removed `config.rs` since its duties are fulfilled by `assets.rs` + +## 2024-11-10 + +A little bit of cleanup to get the thermodynamics plugin functional and neat. +Also getting the balloon to appear as a simple sphere in the simulation but it's +annoyingly full of boilerplate and other setup unrelated to the physics. + +It's a lot harder than I thought to compute the volume of a mesh, mostly because +meshes are surfaces and not volumes. +[Possible, but not simple.](https://github.com/dimforge/parry/blob/68abfc1c22c0beb6d8eba11d57acbb29b4837577/src/mass_properties/mass_properties_trimesh3d.rs#L155) +Instead of having everything be _directly_ coupled, it is suffiecient to have +all of these systems (volume, mesh, collider, etc.) update together based on +the same state. + +Instead, I'll make the debug view show the current atmospheric conditions at a +given point in the simulation. + +- Added an About window. +- Added `bevy-inspector-egui` plugin for debugging. +- Removed `VolumetricBody` from `IdealGas` to make things more direct. +- Removed `compute_volume_from_mesh` function. + +## 2024-11-09 + +Pushing forward toward bringing the basic features into Bevy ECS. Today I am +focusing on the ideal gas law and dynamics functions. + +- Moved basic equations around, like into a new thermodynamics module +- Added custom structs for `Temperature`, `Pressure`, `Volume`, `Density`, and + `MolarMass` to abstract unit conversions and such. It is safer to work with + structs than bare `f32`s for physics equations. + +## 2024-11-08 + +Focusing on getting the app running as a simple splash screen followed by a +loading screen followed by the running application. Loading the configuration +file and setting up the simulation from there. + +I am having trouble getting the Bevy asset server to load the configuration +file. It seems to be advancing to the running state before the configuration +file is loaded. I added a splash screen to practice with state transitions. + +I also added dev tools provided by the Bevy Quickstart Template +([link](https://github.com/TheBevyFlock/bevy_new_2d/blob/main/src/dev_tools.rs)) +to help with debugging. This is toggled with the `F3` key by default, and only +added when the `dev` feature is enabled (it is not enabled by default when +building with cargo and omitted from the release build). I borrowed some other +patterns from the Bevy Quickstart Template for the asset tracking and a few +small things. + +I'm having trouble integrating simple things like loading in values from the +config files. I can't tell what's getting in my way. Time to move things out of +configs and into hard-coded statics. + +Another note about migrating to Bevy is that the engine assumes we are doing +graphics outside of a terminal. We could implement a terminal UI or continue to +support a CLI, but that's 2xs or 3x the work to maintain multiples. Since I plan +to leverage the 3D physics (with things like raycasting) we will drop CLI +support for now. + +- Added splash screen to the application. +- Changed the generic asset loader to a configuration loader for now. +- Added asset tracking plugin. +- Added dev tools provided by the Bevy Quickstart Template. +- Added `dev` feature flag and Bevy build optimiztions to `Cargo.toml`. +- Added `lib.rs` and moved some things around to clean up the root directory. +- Replaced all logging with bevy's built-in logging plugin. + +## 2024-11-07 + +I am switching to Bevy for the simulation. Bevy is a game engine built upon an +ECS framework. It allows for high performance, multi-threaded, dynamic +simulations with 3D graphics and rich interactive GUIs. + +The first reason is that I want to be able to spend more time on the +interactions between the HAB components and less time on the fundamental physics +and simulation scheduling loop. Bevy has a very nice schedule system that allows +for easy parallelization of systems. It also has a component system that allows +me to keep all the logic for the physics systems close to the objects that they +act on. For example, all of the solid bodies that will need to have drag applied +will have a `Body` component, and the logic to calculate the drag on those +bodies will be computed from the their mesh using custom colliders and forces on +top of the physics engine, [Avian](https://github.com/Jondolf/avian), that takes +care of equations of motion, collisions, and constraints. + +The second reason is that I want to be able to run the simulation in a web +browser. Bevy has a web backend that allows for this and very nice tools for +visualizing the simulation state. It also has first-class support for +[Egui](https://github.com/emilk/egui) which is a library for building +interactive GUIs with [Bevy Egui](https://github.com/mvlabat/bevy_egui), and +first-class support for loading assets like configs, 3D models, and textures. + +The first thing I want to do is to get a simple ballistic flight working in +Bevy so that I can validate the fundamental physics assumptions that I have +made. To do this I'll need to start duplicating some of the functionality that +I had in the previous Rust simulation. Namely, I need to: + +1. Initialize the atmosphere. +2. Create the solid bodies (balloon, parachute, payload). +3. Set up a schedule for computing the physics updates every tick/frame. + +The Bevy schedule system will completely replace the threaded event loop that I +had in the previous simulation, including things like `SimCommands`, +`SimOutput`, `SimInstant`, and the `AsyncSim` struct. + +- **Main Application (`src/main.rs`):** + - Integrated `simulator::SimulatorPlugins` to replace the previous UI plugins + setup. + +- **Atmosphere Module (`src/simulator/atmosphere.rs`):** + - Introduced a new `Atmosphere` struct to model atmospheric conditions using + the US Standard Atmosphere, 1976. + - Added functions for calculating temperature and pressure based on altitude. + +- **Balloon Module (`src/simulator/balloon.rs`):** + - Added `BalloonPlugin` struct implementing the `Plugin` trait for Bevy + integration. + +- **Forces Module (`src/simulator/forces.rs`):** + - Updated to use the new `atmosphere::Atmosphere` for buoyancy and drag + calculations. + +- **Gas Module (`src/simulator/gas.rs`):** + - Moved `Atmosphere` struct to its own module. + - Added utility functions for calculating gas volume and density. + +- **Simulator Module (`src/simulator/mod.rs`):** + - Added new modules: `atmosphere` and `units`. + - Converted `SimulatorPlugin` to `SimulatorPlugins` as a `PluginGroup`. + +- **Units Module (`src/simulator/units.rs`):** + - Added a utility function to convert Celsius to Kelvin. diff --git a/docs/newtonian_aero.md b/docs/newtonian_aero.md new file mode 100644 index 0000000..ef9e184 --- /dev/null +++ b/docs/newtonian_aero.md @@ -0,0 +1,127 @@ +# Introduction to Newtonian Aerodynamics + +This is a summary of the paper +[Newtonian Aerodynamics for General Body Shapes](https://ntrs.nasa.gov/citations/19660012440) +and its implementation in our simulation. It is generated using ChatGPT +`o1-mini`. + +## 1. Introduction to Newtonian Aerodynamics + +Newtonian Aerodynamics applies Newtonian mechanics to predict aerodynamic forces +such as drag and lift on objects moving through a fluid (like air). This +approach is particularly useful for objects in supersonic or hypersonic flows. + +## 2. Fundamental Concepts + +- *Drag Force (Fₓ)* is the resistance force acting opposite to the direction of + an object's velocity. + +- *Lift Force (Fᵧ)* is the force perpendicular to the direction of motion, often + associated with winged objects. + +- *Flow Regimes* + + - **Subsonic**: Flow speeds less than the speed of sound. + - **Supersonic**: Flow speeds greater than the speed of sound. + - **Hypersonic**: Flow speeds significantly higher than supersonic speeds. + +## 3. Newtonian Impact Theory + +The Newtonian Impact Theory estimates aerodynamic forces based on the assumption +that fluid particles impacting the body either stick or rebound elastically. + +### Key Assumptions + +1. The fluid is treated as a collection of particles. +2. Particles impacting the body either stick or reflect off elastically. +3. The flow is steady and uniform far from the body. +4. Viscous effects are neglected. + +### Mathematical Formulation + +#### a. Pressure Coefficient (Cₚ) + +\[ C_p = \frac{P - P_{\infty}}{\frac{1}{2} \rho V_{\infty}^2} \] Where: + +- \( P \) = Local pressure on the body. +- \( P_{\infty} \) = Free-stream pressure. +- \( \rho \) = Fluid density. +- \( V_{\infty} \) = Free-stream velocity. + +#### b. Newtonian Drag Coefficient (C_d) + +\[ C_d = \int \left(1 - \cos(\theta)\right) dA \] Where: + +- \( \theta \) = Angle between the local surface normal and the free-stream + velocity vector. +- \( dA \) = Differential area element. + +## 4. Calculating Aerodynamic Forces + +### a. Drag Force (Fₓ) + +\[ F_x = \frac{1}{2} \rho V_{\infty}^2 C_d A \] Where: + +- \( A \) = Reference area (e.g., frontal area). + +### b. Lift Force (Fᵧ) + +\[ F_y = \frac{1}{2} \rho V_{\infty}^2 C_l A \] Where: + +- \( C_l \) = Lift coefficient. + +## 5. Application to General Body Shapes + +The paper extends Newtonian Impact Theory to arbitrary body shapes by +decomposing the body's surface into differential elements, calculating local +impact angles, and integrating contributions to overall aerodynamic forces. + +### a. Surface Sampling + +To handle complex geometries, the body's surface is sampled at numerous points. +Each point's normal vector and its angle relative to the flow direction are +determined. + +### b. Numerical Integration + +Due to the complexity of analytical solutions for irregular shapes, numerical +techniques are employed to perform surface integrals for drag and lift. + +## 6. Enhancements and Corrections + +While Newtonian Impact Theory provides foundational estimates, the paper +discusses enhancements for improved accuracy, including: + +- Correction factors based on empirical data. +- Considering how the body's orientation relative to the flow affects the + distribution of impact angles. + +## 7. Practical Implementation Considerations + +### a. Computational Efficiency + +Sampling the surface at too many points can be computationally intensive. +Optimal sampling densities should balance accuracy and performance. + +### b. Integration with Physics Engines + +Aerodynamic forces must be integrated with the physics engine's update loops, +ensuring forces are applied correctly based on the body's current state and +orientation. + +### c. Validation + +Comparing simulation results with experimental data to validate and calibrate +the implemented models. + +## 8. Summary of Mathematical Steps for Implementation + +1. **Surface Sampling**: Discretize the body's surface into numerous points. +2. **Local Impact Angle Calculation**: Compute the angle between each normal + vector and the free-stream velocity vector. +3. **Pressure Coefficient Determination**: Use the local impact angles to + estimate the pressure coefficient \( C_p \). +4. **Drag and Lift Force Integration**: Numerically integrate the contributions + of all sampled points to compute total drag and lift forces. +5. **Normalization and Scaling**: Adjust the computed forces based on fluid + density and relative velocity squared. diff --git a/src/app3d/mod.rs b/src/app3d/mod.rs new file mode 100644 index 0000000..40c0245 --- /dev/null +++ b/src/app3d/mod.rs @@ -0,0 +1,6 @@ +mod scene; +mod ui; + +// Re-export the plugins so they can be added to the app with `app.add_plugins`. +pub use scene::ScenePlugin; +pub use ui::InterfacePlugins; diff --git a/src/app3d/scene/camera.rs b/src/app3d/scene/camera.rs new file mode 100644 index 0000000..8372f55 --- /dev/null +++ b/src/app3d/scene/camera.rs @@ -0,0 +1,65 @@ +use bevy::prelude::*; +use bevy_panorbit_camera::{PanOrbitCamera, PanOrbitCameraPlugin}; +use avian3d::math::TAU; + +use crate::controls::KeyBindingsConfig; + +pub struct CameraPlugin; + +impl Plugin for CameraPlugin { + fn build(&self, app: &mut App) { + app.add_plugins(PanOrbitCameraPlugin); + app.add_systems(Startup, setup); + app.add_systems(Update, toggle_camera_controls_system); + } +} + +fn setup(mut commands: Commands, key_bindings: Res) { + commands.spawn(( + // Note we're setting the initial position below with yaw, pitch, and radius, hence + // we don't set transform on the camera. + Camera3dBundle::default(), + PanOrbitCamera { + // Set focal point (what the camera should look at) + focus: Vec3::new(0.0, 1.0, 0.0), + // Set the starting position, relative to focus (overrides camera's transform). + yaw: Some(TAU / 8.0), + pitch: Some(TAU / 8.0), + radius: Some(5.0), + // Set limits on rotation and zoom + yaw_upper_limit: Some(TAU / 4.0), + yaw_lower_limit: Some(-TAU / 4.0), + pitch_upper_limit: Some(TAU / 3.0), + pitch_lower_limit: Some(-TAU / 3.0), + zoom_upper_limit: Some(100.0), + zoom_lower_limit: 1.0, + // Adjust sensitivity of controls + orbit_sensitivity: 1.5, + pan_sensitivity: 0.5, + zoom_sensitivity: 0.5, + // Allow the camera to go upside down + allow_upside_down: true, + // Change the controls (these match Blender) + button_orbit: key_bindings.camera_controls.button_orbit, + button_pan: key_bindings.camera_controls.button_pan, + modifier_pan: key_bindings.camera_controls.modifier_pan, + // Reverse the zoom direction + reversed_zoom: false, + ..default() + }, + )); +} + +// This is how you can change config at runtime. +// Press 'T' to toggle the camera zoom direction. +fn toggle_camera_controls_system( + key_input: Res>, + key_bindings: Res, + mut pan_orbit_query: Query<&mut PanOrbitCamera>, +) { + if key_input.just_pressed(key_bindings.camera_controls.toggle_zoom_direction) { + for mut pan_orbit in pan_orbit_query.iter_mut() { + pan_orbit.reversed_zoom = !pan_orbit.reversed_zoom; + } + } +} diff --git a/src/app3d/scene/mod.rs b/src/app3d/scene/mod.rs new file mode 100644 index 0000000..cb7c6da --- /dev/null +++ b/src/app3d/scene/mod.rs @@ -0,0 +1,57 @@ +mod camera; + +use avian3d::prelude::*; +use bevy::prelude::*; + +pub struct ScenePlugin; + +impl Plugin for ScenePlugin { + fn build(&self, app: &mut App) { + app.add_plugins(camera::CameraPlugin); + app.add_systems(Startup, simple_scene); + } +} + +/// set up a simple 3D scene +fn simple_scene( + mut commands: Commands, + mut meshes: ResMut>, + mut materials: ResMut>, +) { + // light + commands.spawn(( + Name::new("Light"), + PointLightBundle { + point_light: PointLight { + shadows_enabled: true, + ..default() + }, + transform: Transform::from_xyz(4.0, 8.0, 4.0), + ..default() + }, + )); + // ground + let ground_size = Vec3::new(4.0, 0.1, 4.0); + commands.spawn(( + Name::new("Ground"), + PbrBundle { + mesh: meshes.add(Cuboid::new(ground_size.x, ground_size.y, ground_size.z)), + material: materials.add(Color::srgba(0.75, 0.75, 0.75, 0.1)), + transform: Transform::from_translation(Vec3::new(0.0, -2.0, 0.0)), + ..default() + }, + RigidBody::Static, + Collider::cuboid(ground_size.x, ground_size.y, ground_size.z), + )); + commands.spawn(( + Name::new("Ceiling"), + PbrBundle { + mesh: meshes.add(Cuboid::new(ground_size.x, ground_size.y, ground_size.z)), + material: materials.add(Color::srgba(0.75, 0.75, 0.75, 0.1)), + transform: Transform::from_translation(Vec3::new(0.0, 3.0, 0.0)), + ..default() + }, + RigidBody::Static, + Collider::cuboid(ground_size.x, ground_size.y, ground_size.z), + )); +} diff --git a/src/app3d/scene/vectors.rs b/src/app3d/scene/vectors.rs new file mode 100644 index 0000000..8f21a06 --- /dev/null +++ b/src/app3d/scene/vectors.rs @@ -0,0 +1,35 @@ +use bevy::{prelude::*, render::primitives::Aabb}; +use avian3d::math::PI; +use crate::simulator::forces::Force; + +pub struct ForceVectorPlugin; + +impl Plugin for ForceVectorPlugin { + fn build(&self, app: &mut App) { + app.add_systems(Update, draw_axes); + } +} + +/// System to visualize all force vectors +// fn visualize_forces( +// gizmos: Res, +// query: Query<(&Transform, &dyn Force)>, +// ) { +// for (transform, all_forces) in query.iter() { +// let origin = transform.translation; +// let segments: Vec<(Vec3, Vec3)> = all_forces.iter().map(|force| { +// let force_vector = force.force(); +// (origin, origin + force_vector) +// }).collect(); + +// } +// } + +// This system draws the axes based on the cube's transform, with length based on the size of +// the entity's axis-aligned bounding box (AABB). +fn draw_axes(mut gizmos: Gizmos, query: Query<(&Transform, &Aabb)>) { + for (&transform, &aabb) in &query { + let length = aabb.half_extents.length() * 0.1; + gizmos.axes(transform, length); + } +} diff --git a/src/app3d/ui/core.rs b/src/app3d/ui/core.rs new file mode 100644 index 0000000..b1b8d9b --- /dev/null +++ b/src/app3d/ui/core.rs @@ -0,0 +1,56 @@ +use bevy::{app::PluginGroupBuilder, prelude::*}; +use iyes_perf_ui::prelude::*; + +use crate::controls::KeyBindingsConfig; +use crate::simulator::SimState; + +use super::*; + +/// A plugin group that includes all interface-related plugins +pub struct InterfacePlugins; + +impl PluginGroup for InterfacePlugins { + fn build(self) -> PluginGroupBuilder { + PluginGroupBuilder::start::() + .add(CoreUiPlugin) + .add(PausePlayPlugin) + .add(monitors::MonitorsPlugin) + } +} + +/// Base UI plugin. This sets up the base plugins that all other ui plugins +/// need. .Placeholder for now +pub struct CoreUiPlugin; + +impl Plugin for CoreUiPlugin { + fn build(&self, app: &mut App) { + app.add_plugins(( + PerfUiPlugin, + #[cfg(feature = "dev")] + dev_tools::plugin, + )); + } +} + +pub struct PausePlayPlugin; + +impl Plugin for PausePlayPlugin { + fn build(&self, app: &mut App) { + app.add_systems(Update, toggle_pause); + } +} + +fn toggle_pause( + sim_state: Res>, + mut next_state: ResMut>, + key_input: Res>, + key_bindings: Res, +) { + if key_input.just_pressed(key_bindings.time_controls.toggle_pause) { + match sim_state.as_ref().get() { + SimState::Stopped => next_state.set(SimState::Running), + SimState::Running => next_state.set(SimState::Stopped), + _ => () + } + } +} diff --git a/src/app3d/ui/dev_tools.rs b/src/app3d/ui/dev_tools.rs new file mode 100644 index 0000000..dd4dae6 --- /dev/null +++ b/src/app3d/ui/dev_tools.rs @@ -0,0 +1,108 @@ +//! Development tools for the game. This plugin is only enabled in dev builds. +#[allow(unused_imports)] +use bevy::{ + // dev_tools::states::log_transitions, + diagnostic::{ + FrameTimeDiagnosticsPlugin, EntityCountDiagnosticsPlugin, + SystemInformationDiagnosticsPlugin, + }, + input::common_conditions::input_just_pressed, + prelude::*, +}; +use iyes_perf_ui::prelude::*; + +#[cfg(not(target_arch = "wasm32"))] +use bevy::pbr::wireframe::{WireframeConfig, WireframePlugin}; + +use avian3d::debug_render::PhysicsDebugPlugin; + +use crate::controls::KeyBindingsConfig; + +pub(super) fn plugin(app: &mut App) { + // Toggle the debug overlay for UI. + app.add_plugins(( + // physics + PhysicsDebugPlugin::default(), + // performance + FrameTimeDiagnosticsPlugin, + EntityCountDiagnosticsPlugin, + SystemInformationDiagnosticsPlugin, + // rendering + #[cfg(not(target_arch = "wasm32"))] + WireframePlugin, + )); + app.add_systems( + Update, + toggle_debug_ui + .before(iyes_perf_ui::PerfUiSet::Setup), + ); + + // Wireframe doesn't work on WASM + #[cfg(not(target_arch = "wasm32"))] + app.add_systems( + Update, + toggle_wireframe, + ); + + #[cfg(feature = "inspect")] + { + use bevy_inspector_egui::quick::WorldInspectorPlugin; + app.add_plugins(WorldInspectorPlugin::new()); + } +} + +/// Toggle the debug overlay +fn toggle_debug_ui( + mut commands: Commands, + q_root: Query>, + key_input: Res>, + key_bindings: Res, +) { + if key_input.just_pressed(key_bindings.debug_controls.toggle_perf_ui) { + if let Ok(e) = q_root.get_single() { + // despawn the existing Perf UI + commands.entity(e).despawn_recursive(); + } else { + // create a simple Perf UI with default settings + // and all entries provided by the crate: + commands.spawn(( + PerfUiRoot { + // set a fixed width to make all the bars line up + values_col_width: Some(160.0), + ..Default::default() + }, + // when we have lots of entries, we have to group them + // into tuples, because of Bevy Rust syntax limitations + ( + PerfUiWidgetBar::new(PerfUiEntryFPS::default()), + PerfUiWidgetBar::new(PerfUiEntryFPSWorst::default()), + PerfUiWidgetBar::new(PerfUiEntryFrameTime::default()), + PerfUiWidgetBar::new(PerfUiEntryFrameTimeWorst::default()), + PerfUiWidgetBar::new(PerfUiEntryEntityCount::default()), + ), + ( + PerfUiEntryRunningTime::default(), + PerfUiEntryClock::default(), + ), + ( + PerfUiEntryCursorPosition::default(), + // PerfUiEntryWindowResolution::default(), + // PerfUiEntryWindowScaleFactor::default(), + // PerfUiEntryWindowMode::default(), + // PerfUiEntryWindowPresentMode::default(), + ), + )); + } + } +} + +#[cfg(not(target_arch = "wasm32"))] +fn toggle_wireframe( + mut wireframe_config: ResMut, + key_input: Res>, + key_bindings: Res, +) { + if key_input.just_pressed(key_bindings.debug_controls.toggle_wireframe) { + wireframe_config.global = !wireframe_config.global; + } +} diff --git a/src/app3d/ui/mod.rs b/src/app3d/ui/mod.rs new file mode 100644 index 0000000..af959cc --- /dev/null +++ b/src/app3d/ui/mod.rs @@ -0,0 +1,7 @@ +mod core; +mod monitors; + +#[cfg(feature = "dev")] +mod dev_tools; + +pub use core::InterfacePlugins; diff --git a/src/app3d/ui/monitors.rs b/src/app3d/ui/monitors.rs new file mode 100644 index 0000000..a344b19 --- /dev/null +++ b/src/app3d/ui/monitors.rs @@ -0,0 +1,238 @@ +//! UI for monitoring the simulation. +#![allow(unused_imports)] +use bevy::{ + ecs::system::{lifetimeless::SRes, SystemParam}, + prelude::*, +}; +use iyes_perf_ui::{entry::PerfUiEntry, prelude::*, utils::format_pretty_float}; + +use crate::simulator::{ + forces::{Buoyancy, Drag, Force, Weight}, + SimState, SimulatedBody, +}; + +pub struct MonitorsPlugin; + +impl Plugin for MonitorsPlugin { + fn build(&self, app: &mut App) { + app.add_perf_ui_simple_entry::(); + app.add_perf_ui_simple_entry::(); + app.add_systems(Startup, spawn_monitors); + app.add_systems(Update, update_force_monitor_values); + app.init_resource::(); + } +} + +fn spawn_monitors(mut commands: Commands) { + commands.spawn(( + PerfUiRoot { + position: PerfUiPosition::BottomRight, + ..default() + }, + SimStateMonitor::default(), + ForceMonitor::default(), + )); +} + +#[derive(Component)] +struct SimStateMonitor { + /// The label text to display, to allow customization + pub label: String, + /// Should we display units? + pub display_units: bool, + /// Highlight the value if it goes above this threshold + pub threshold_highlight: Option, + /// Support color gradients! + pub color_gradient: ColorGradient, + /// Width for formatting the string + pub digits: u8, + /// Precision for formatting the string + pub precision: u8, + /// Required to ensure the entry appears in the correct place in the Perf UI + pub sort_key: i32, +} + +impl Default for SimStateMonitor { + fn default() -> Self { + SimStateMonitor { + label: String::new(), + display_units: false, + threshold_highlight: Some(10.0), + color_gradient: ColorGradient::new_preset_gyr(0.0, 10.0, 100.0).unwrap(), + digits: 7, + precision: 0, + sort_key: iyes_perf_ui::utils::next_sort_key(), + } + } +} + +impl PerfUiEntry for SimStateMonitor { + type Value = SimState; + type SystemParam = SRes>; + + fn label(&self) -> &str { + if self.label.is_empty() { + "Simulation State" + } else { + &self.label + } + } + + fn sort_key(&self) -> i32 { + self.sort_key + } + + fn format_value(&self, value: &Self::Value) -> String { + match value { + SimState::Running => String::from("Running"), + SimState::Stopped => String::from("Stopped"), + SimState::Anomaly => String::from("ANOMALY") + } + } + + fn update_value( + &self, + state: &mut ::Item<'_, '_>, + ) -> Option { + Some(*state.as_ref().get()) + } + + // (optional) We should add a width hint, so that the displayed + // strings in the UI can be correctly aligned. + // This value represents the largest length the formatted string + // is expected to have. + fn width_hint(&self) -> usize { + // there is a helper we can use, since we use `format_pretty_float` + let w = iyes_perf_ui::utils::width_hint_pretty_float(self.digits, self.precision); + if self.display_units { + w + 2 + } else { + w + } + } + + // (optional) Called every frame to determine if a custom color should be used for the value + fn value_color(&self, value: &Self::Value) -> Option { + match *value { + SimState::Running => self.color_gradient.get_color_for_value(0.0), + SimState::Stopped => self.color_gradient.get_color_for_value(10.0), + _ => self.color_gradient.get_color_for_value(100.0), + } + } + + // (optional) Called every frame to determine if the value should be highlighted + fn value_highlight(&self, value: &Self::Value) -> bool { + self.threshold_highlight + .map(|_| value == &SimState::Stopped) + .unwrap_or(false) + } +} + +#[derive(Resource, Reflect, Default)] +struct ForceMonitorResource { + pub weight: Vec3, + pub buoyancy: Vec3, + pub drag: Vec3, +} + +fn update_force_monitor_values( + mut force_resource: ResMut, + forces: Query<(&Weight, &Buoyancy, &Drag), With>, +) { + for (weight, bouyancy, drag) in forces.iter() { + // assume there's only one simulated body for now + force_resource.weight = weight.force(); + force_resource.buoyancy = bouyancy.force(); + force_resource.drag = drag.force(); + } +} + +#[derive(Component)] +struct ForceMonitor { + /// The label text to display, to allow customization + pub label: String, + /// Should we display units? + pub display_units: bool, + /// Highlight the value if it goes above this threshold + #[allow(dead_code)] + pub threshold_highlight: Option, + /// Support color gradients! + #[allow(dead_code)] + pub color_gradient: ColorGradient, + /// Width for formatting the string + pub digits: u8, + /// Precision for formatting the string + pub precision: u8, + /// Required to ensure the entry appears in the correct place in the Perf UI + pub sort_key: i32, +} + +impl Default for ForceMonitor { + fn default() -> Self { + ForceMonitor { + label: String::new(), + display_units: true, + threshold_highlight: Some(10.0), + color_gradient: ColorGradient::new_preset_gyr(0.0, 10.0, 100.0).unwrap(), + digits: 5, + precision: 2, + sort_key: iyes_perf_ui::utils::next_sort_key(), + } + } +} + +impl PerfUiEntry for ForceMonitor { + type Value = (f32, f32, f32); + type SystemParam = SRes; + + fn label(&self) -> &str { + if self.label.is_empty() { + "Forces" + } else { + &self.label + } + } + + fn sort_key(&self) -> i32 { + self.sort_key + } + + fn format_value(&self, value: &Self::Value) -> String { + // we can use a premade helper function for nice-looking formatting + let mut f_g = format_pretty_float(self.digits, self.precision, value.0 as f64); + let mut f_b = format_pretty_float(self.digits, self.precision, value.1 as f64); + let mut f_d = format_pretty_float(self.digits, self.precision, value.2 as f64); + // (and append units to it) + if self.display_units { + f_g.push_str(" N"); + f_b.push_str(" N"); + f_d.push_str(" N"); + } + format!("Fg {:} Fb {:} Fd {:}", f_g, f_b, f_d) + } + + fn update_value( + &self, + force_resource: &mut ::Item<'_, '_>, + ) -> Option { + Some(( + force_resource.weight.length(), + force_resource.buoyancy.length(), + force_resource.drag.length(), + )) + } + + // (optional) We should add a width hint, so that the displayed + // strings in the UI can be correctly aligned. + // This value represents the largest length the formatted string + // is expected to have. + fn width_hint(&self) -> usize { + // there is a helper we can use, since we use `format_pretty_float` + let w = iyes_perf_ui::utils::width_hint_pretty_float(self.digits, self.precision); + if self.display_units { + w + 2 + } else { + w + } + } +} diff --git a/src/cli.rs b/src/cli.rs deleted file mode 100644 index 8b28a77..0000000 --- a/src/cli.rs +++ /dev/null @@ -1,97 +0,0 @@ -use std::path::PathBuf; - -use clap::{Parser, Subcommand}; -use log::error; - -use crate::simulator::schedule::{AsyncSim, Rate}; - -#[derive(Parser)] -#[clap(author, version, about)] -struct Cli { - #[clap(subcommand)] - command: Commands, -} - -#[derive(Subcommand)] -enum Commands { - /// Start a new simulation process - /// - /// Configure an asynchronous physics simulation in the background. This - /// simulation runs on the MFC with flight software code running in the - /// loop and logs the simulation output to a CSV file. - Start { - /// Sets a custom simulation config file - #[clap( - short, - long, - value_name = "TOML", - default_value = "config/default.toml" - )] - config: PathBuf, - - /// Sets a custom output file - #[clap(short, long, value_name = "CSV", default_value = "./out.csv")] - outpath: PathBuf, - }, - - /// Inspect a physics parameter in an existing simulation - Get { - /// Parameter to be inspect - param: String, - }, - - /// Modify a physics parameter in an existing simulation - Set { - /// Parameter to be modified - param: String, - /// New value to set - value: String, - }, - - /// Open a graphical user interface instead of the terminal interface - Gui, -} - -pub fn parse_inputs() { - // parse CLI input args and options - let cli = Cli::parse(); - match &cli.command { - Commands::Start { config, outpath } => { - start_sim(config, outpath); - } - Commands::Gui => { - #[cfg(feature = "gui")] - start_gui(); - - #[cfg(not(feature = "gui"))] - error!("GUI feature not enabled. Reinstall with `--features gui`") - } - _ => { - error!("Command not implemented yet!") - } - } -} - -pub fn start_sim(config: &PathBuf, outpath: &PathBuf) { - // initialize the simulation - let mut sim = AsyncSim::new(config, outpath.clone()); - let mut rate_sleeper = Rate::new(1.0); - - // start the sim - sim.start(); - loop { - sim.get_sim_output(); - rate_sleeper.sleep(); - } -} - -#[cfg(feature = "gui")] -pub fn start_gui() { - use crate::gui; - let native_options = eframe::NativeOptions::default(); - let _ = eframe::run_native( - "Mission Control", - native_options, - Box::new(|cc| Box::new(gui::Shell::new(cc))), - ); -} diff --git a/src/controls.rs b/src/controls.rs new file mode 100644 index 0000000..68ea9ae --- /dev/null +++ b/src/controls.rs @@ -0,0 +1,68 @@ +use bevy::prelude::*; + +pub struct ControlsPlugin; + +impl Plugin for ControlsPlugin { + fn build(&self, app: &mut App) { + app.init_resource::(); + } +} + +#[derive(Resource, Default)] +pub struct KeyBindingsConfig { + pub camera_controls: CameraControls, + pub debug_controls: DebugControls, + pub time_controls: TimeControls, +} + +#[derive(Reflect)] +pub struct CameraControls { + pub modifier_pan: Option, + pub button_pan: MouseButton, + pub button_orbit: MouseButton, + pub toggle_zoom_direction: KeyCode, +} + +#[derive(Reflect)] +pub struct DebugControls { + pub toggle_wireframe: KeyCode, + pub toggle_inspector: KeyCode, + pub toggle_perf_ui: KeyCode, +} + +#[derive(Reflect)] +pub struct TimeControls { + pub toggle_pause: KeyCode, +} + +// ============================ DEFAULT KEYBINDINGS ============================ + +/// Defaults follow Blender conventions +impl Default for CameraControls { + fn default() -> Self { + Self { + modifier_pan: Some(KeyCode::ShiftLeft), + button_pan: MouseButton::Middle, + button_orbit: MouseButton::Middle, + toggle_zoom_direction: KeyCode::KeyZ, + } + } +} + +impl Default for DebugControls { + fn default() -> Self { + Self { + toggle_wireframe: KeyCode::KeyW, + toggle_inspector: KeyCode::KeyI, + toggle_perf_ui: KeyCode::KeyP, + } + } +} + +impl Default for TimeControls { + fn default() -> Self { + Self { + toggle_pause: KeyCode::Space, + } + } +} diff --git a/src/gui/mod.rs b/src/gui/mod.rs deleted file mode 100644 index f228dc4..0000000 --- a/src/gui/mod.rs +++ /dev/null @@ -1,23 +0,0 @@ -mod shell; -mod views; - -pub use shell::Shell; - -/// Something to view in the monitor windows -pub trait View { - fn ui(&mut self, ui: &mut egui::Ui); -} - -/// Something to view -pub trait UiPanel { - /// Is the monitor enabled for this integration? - fn is_enabled(&self, _ctx: &egui::Context) -> bool { - true - } - - /// `&'static` so we can also use it as a key to store open/close state. - fn name(&self) -> &'static str; - - /// Show windows, etc - fn show(&mut self, ctx: &egui::Context, open: &mut bool); -} diff --git a/src/gui/shell.rs b/src/gui/shell.rs deleted file mode 100644 index 8e5dd88..0000000 --- a/src/gui/shell.rs +++ /dev/null @@ -1,264 +0,0 @@ -use std::{ - collections::BTreeSet, - path::PathBuf, - sync::{Arc, Mutex}, - thread::JoinHandle, -}; -use egui::{Context, Modifiers, ScrollArea, Ui}; - -use super::UiPanel; -use crate::simulator::{ - io::SimOutput, - config::{self, Config}, - schedule::AsyncSim, -}; - -/// A menu bar in which you can select different info windows to show. -#[cfg_attr(feature = "serde", derive(serde::Deserialize, serde::Serialize))] -#[cfg_attr(feature = "serde", serde(default))] -pub struct Shell { - screens: Screens, - config: Config, - output: Option, - run_handle: Option>, -} - -impl Default for Shell { - fn default() -> Self { - let default_config_path = PathBuf::from("config/default.toml"); - let config = config::parse_from_file(&default_config_path); - Self { - screens: Screens::default(), - config, - output: None, - run_handle: None, - } - } -} - -impl Shell { - /// Called once before the first frame. - pub fn new(_cc: &eframe::CreationContext<'_>) -> Self { - // This is also where you can customize the look and feel of egui using - // `cc.egui_ctx.set_visuals` and `cc.egui_ctx.set_fonts`. - Default::default() - } - - /// Show the app ui (menu bar and windows). - pub fn ui(&mut self, ctx: &Context) { - egui::SidePanel::left("mission_control_panel") - .resizable(false) - .default_width(150.0) - .show(ctx, |ui| { - ui.vertical_centered(|ui| { - ui.heading("yahs"); - }); - - ui.separator(); - - use egui::special_emojis::GITHUB; - ui.hyperlink_to( - format!("{GITHUB} yahs on GitHub"), - "https://github.com/brickworks/yahs", - ); - ui.hyperlink_to( - format!("{GITHUB} @philiplinden"), - "https://github.com/philiplinden", - ); - - ui.separator(); - egui::widgets::global_dark_light_mode_buttons(ui); - ui.separator(); - self.screen_list_ui(ui); - ui.separator(); - self.sim_control_buttons(ui); - ui.separator(); - }); - - egui::TopBottomPanel::top("menu_bar").show(ctx, |ui| { - egui::menu::bar(ui, |ui| { - file_menu_button(ui); - }); - }); - - self.show_windows(ctx); - } - - /// Show the open windows. - fn show_windows(&mut self, ctx: &Context) { - self.screens.windows(ctx); - } - - fn screen_list_ui(&mut self, ui: &mut egui::Ui) { - ScrollArea::vertical().show(ui, |ui| { - ui.with_layout(egui::Layout::top_down_justified(egui::Align::LEFT), |ui| { - self.screens.checkboxes(ui); - }); - }); - } - - fn sim_control_buttons(&mut self, ui: &mut egui::Ui) { - if ui.button("Simulate").clicked() { - - if self.run_handle.is_some() { - panic!("Can't start again, sim already ran. Need to stop.") - } - let config = self.config.clone(); - let outpath = PathBuf::from("out.csv"); - let init_state = Arc::new(Mutex::new(SimOutput::default())); - AsyncSim::run_sim(config, init_state, outpath) - } - } -} - -impl eframe::App for Shell { - /// Called each time the UI needs repainting, which may be many times per second. - fn update(&mut self, ctx: &egui::Context, _frame: &mut eframe::Frame) { - // Put your widgets into a `SidePanel`, `TopBottomPanel`, `CentralPanel`, `Window` or `Area`. - // For inspiration and more examples, go to https://emilk.github.io/egui - egui::CentralPanel::default().show(ctx, |_ui| { - // The central panel the region left after adding TopPanel's and SidePanel's - self.ui(ctx); - }); - egui::TopBottomPanel::bottom("powered_by_eframe").show(ctx, |ui| { - ui.with_layout(egui::Layout::bottom_up(egui::Align::LEFT), |ui| { - powered_by_egui_and_eframe(ui); - egui::warn_if_debug_build(ui); - }); - }); - } -} - -fn powered_by_egui_and_eframe(ui: &mut egui::Ui) { - ui.horizontal(|ui| { - ui.spacing_mut().item_spacing.x = 0.0; - ui.label("Powered by "); - ui.hyperlink_to("egui", "https://github.com/emilk/egui"); - ui.label(" and "); - ui.hyperlink_to( - "eframe", - "https://github.com/emilk/egui/tree/master/crates/eframe", - ); - ui.label("."); - }); -} - -fn file_menu_button(ui: &mut Ui) { - let organize_shortcut = - egui::KeyboardShortcut::new(Modifiers::CTRL | Modifiers::SHIFT, egui::Key::O); - let reset_shortcut = - egui::KeyboardShortcut::new(Modifiers::CTRL | Modifiers::SHIFT, egui::Key::R); - - // NOTE: we must check the shortcuts OUTSIDE of the actual "File" menu, - // or else they would only be checked if the "File" menu was actually open! - - if ui.input_mut(|i| i.consume_shortcut(&organize_shortcut)) { - ui.ctx().memory_mut(|mem| mem.reset_areas()); - } - - if ui.input_mut(|i| i.consume_shortcut(&reset_shortcut)) { - ui.ctx().memory_mut(|mem| *mem = Default::default()); - } - - ui.menu_button("View", |ui| { - ui.set_min_width(220.0); - ui.style_mut().wrap = Some(false); - - // On the web the browser controls the zoom - #[cfg(not(target_arch = "wasm32"))] - { - egui::gui_zoom::zoom_menu_buttons(ui); - ui.weak(format!( - "Current zoom: {:.0}%", - 100.0 * ui.ctx().zoom_factor() - )) - .on_hover_text("The UI zoom level, on top of the operating system's default value"); - ui.separator(); - } - - if ui - .add( - egui::Button::new("Organize Windows") - .shortcut_text(ui.ctx().format_shortcut(&organize_shortcut)), - ) - .clicked() - { - ui.ctx().memory_mut(|mem| mem.reset_areas()); - ui.close_menu(); - } - - if ui - .add( - egui::Button::new("Reset app memory") - .shortcut_text(ui.ctx().format_shortcut(&reset_shortcut)), - ) - .on_hover_text("Forget scroll, positions, sizes etc") - .clicked() - { - ui.ctx().memory_mut(|mem| *mem = Default::default()); - ui.close_menu(); - } - }); -} - -// ---------------------------------------------------------------------------- - -#[cfg_attr(feature = "serde", derive(serde::Deserialize, serde::Serialize))] -#[cfg_attr(feature = "serde", serde(default))] -struct Screens { - #[cfg_attr(feature = "serde", serde(skip))] - screens: Vec>, - - open: BTreeSet, -} - -impl Default for Screens { - fn default() -> Self { - Self::from_demos(vec![ - Box::::default(), - Box::::default(), - Box::::default(), - ]) - } -} - -impl Screens { - pub fn from_demos(screens: Vec>) -> Self { - let mut open = BTreeSet::new(); - open.insert(super::views::ConfigView::default().name().to_owned()); - - Self { screens, open } - } - - pub fn checkboxes(&mut self, ui: &mut Ui) { - let Self { screens, open } = self; - for screen in screens { - if screen.is_enabled(ui.ctx()) { - let mut is_open = open.contains(screen.name()); - ui.toggle_value(&mut is_open, screen.name()); - set_open(open, screen.name(), is_open); - } - } - } - - pub fn windows(&mut self, ctx: &Context) { - let Self { screens, open } = self; - for screen in screens { - let mut is_open = open.contains(screen.name()); - screen.show(ctx, &mut is_open); - set_open(open, screen.name(), is_open); - } - } -} - -#[cfg_attr(feature = "serde", derive(serde::Deserialize, serde::Serialize))] -#[cfg_attr(feature = "serde", serde(default))] -fn set_open(open: &mut BTreeSet, key: &'static str, is_open: bool) { - if is_open { - if !open.contains(key) { - open.insert(key.to_owned()); - } - } else { - open.remove(key); - } -} diff --git a/src/gui/views.rs b/src/gui/views.rs deleted file mode 100644 index 09ac3ea..0000000 --- a/src/gui/views.rs +++ /dev/null @@ -1,308 +0,0 @@ -use egui::*; -use egui_plot::{ - Bar, BarChart, BoxElem, BoxPlot, BoxSpread, Legend, Line, Plot, -}; - -use crate::gui::View; -use crate::simulator::config::{self, Config}; - -// ---------------------------------------------------------------------------- - -#[derive(PartialEq)] -pub struct ConfigView { - config: Config, - picked_path: Option, -} - -impl Default for ConfigView { - fn default() -> Self { - Self { - config: config::parse_from_file("config/default.toml"), - picked_path: None, - } - } -} -impl super::UiPanel for ConfigView { - fn name(&self) -> &'static str { - "🗠 Config" - } - - fn show(&mut self, ctx: &Context, open: &mut bool) { - Window::new(self.name()) - .open(open) - .vscroll(false) - .default_size(vec2(400.0, 400.0)) - .show(ctx, |ui| self.ui(ui)); - } -} - -impl super::View for ConfigView { - fn ui(&mut self, ui: &mut Ui) { - ui.horizontal(|ui| { - - if ui.button("Import").clicked() { - if let Some(path) = rfd::FileDialog::new().pick_file() { - self.picked_path = Some(path.display().to_string()); - } - } - if let Some(picked_path) = &self.picked_path { - self.config = config::parse_from_file(picked_path); - ui.horizontal(|ui| { - ui.monospace(picked_path); - }); - } - }); - } -} - -// ---------------------------------------------------------------------------- - -#[derive(PartialEq, Default)] -pub struct FlightView; - -impl super::UiPanel for FlightView { - fn name(&self) -> &'static str { - "🗠 Flight" - } - - fn show(&mut self, ctx: &Context, open: &mut bool) { - Window::new(self.name()) - .open(open) - .vscroll(false) - .default_size(vec2(400.0, 400.0)) - .show(ctx, |ui| self.ui(ui)); - } -} - -impl super::View for FlightView { - fn ui(&mut self, ui: &mut Ui) { - ui.horizontal(|ui| { - egui::reset_button(ui, self); - ui.collapsing("Instructions", |ui| { - ui.label("Pan by dragging, or scroll (+ shift = horizontal)."); - ui.label("Box zooming: Right click to zoom in and zoom out using a selection."); - if cfg!(target_arch = "wasm32") { - ui.label("Zoom with ctrl / ⌘ + pointer wheel, or with pinch gesture."); - } else if cfg!(target_os = "macos") { - ui.label("Zoom with ctrl / ⌘ + scroll."); - } else { - ui.label("Zoom with ctrl + scroll."); - } - ui.label("Reset view with double-click."); - }); - }); - ui.separator(); - Plot::new("left-bottom") - .data_aspect(0.5) - .x_axis_label("flight time (s)") - .show(ui, Self::configure_plot); - } -} - -impl FlightView { - fn cos() -> Line { - Line::new(egui_plot::PlotPoints::from_explicit_callback( - move |x| x.cos(), - .., - 100, - )) - } - - fn configure_plot(plot_ui: &mut egui_plot::PlotUi) { - plot_ui.line(Self::cos()); - } -} - -// ---------------------------------------------------------------------------- - -#[derive(PartialEq, Eq)] -enum Chart { - GaussBars, - StackedBars, - BoxPlot, -} - -impl Default for Chart { - fn default() -> Self { - Self::GaussBars - } -} - -#[derive(PartialEq, Default)] -pub struct StatsView { - chart: Chart, -} - -impl super::UiPanel for StatsView { - fn name(&self) -> &'static str { - "🗠 Stats" - } - - fn show(&mut self, ctx: &Context, open: &mut bool) { - use super::View as _; - Window::new(self.name()) - .open(open) - .default_size(vec2(400.0, 400.0)) - .vscroll(false) - .show(ctx, |ui| self.ui(ui)); - } -} - -impl super::View for StatsView { - fn ui(&mut self, ui: &mut Ui) { - ui.horizontal(|ui| { - egui::reset_button(ui, self); - ui.collapsing("Instructions", |ui| { - ui.label("Pan by dragging, or scroll (+ shift = horizontal)."); - ui.label("Box zooming: Right click to zoom in and zoom out using a selection."); - if cfg!(target_arch = "wasm32") { - ui.label("Zoom with ctrl / ⌘ + pointer wheel, or with pinch gesture."); - } else if cfg!(target_os = "macos") { - ui.label("Zoom with ctrl / ⌘ + scroll."); - } else { - ui.label("Zoom with ctrl + scroll."); - } - ui.label("Reset view with double-click."); - }); - }); - ui.separator(); - ui.horizontal(|ui| { - ui.selectable_value(&mut self.chart, Chart::GaussBars, "Histogram"); - ui.selectable_value(&mut self.chart, Chart::StackedBars, "Stacked Bar Chart"); - ui.selectable_value(&mut self.chart, Chart::BoxPlot, "Box Plot"); - }); - ui.separator(); - let _ = match self.chart { - Chart::GaussBars => self.bar_gauss(ui), - Chart::StackedBars => self.bar_stacked(ui), - Chart::BoxPlot => self.box_plot(ui), - }; - } -} - -impl StatsView { - fn bar_gauss(&self, ui: &mut Ui) -> Response { - let chart = BarChart::new( - (-395..=395) - .step_by(10) - .map(|x| x as f64 * 0.01) - .map(|x| { - ( - x, - (-x * x / 2.0).exp() / (2.0 * std::f64::consts::PI).sqrt(), - ) - }) - // The 10 factor here is purely for a nice 1:1 aspect ratio - .map(|(x, f)| Bar::new(x, f * 10.0).width(0.095)) - .collect(), - ) - .color(Color32::LIGHT_BLUE) - .name("Normal Distribution"); - - Plot::new("Normal Distribution Demo") - .legend(Legend::default()) - .clamp_grid(true) - .y_axis_width(3) - .allow_zoom(true) - .allow_drag(true) - .show(ui, |plot_ui| plot_ui.bar_chart(chart)) - .response - } - - fn bar_stacked(&self, ui: &mut Ui) -> Response { - let chart1 = BarChart::new(vec![ - Bar::new(0.5, 1.0).name("Day 1"), - Bar::new(1.5, 3.0).name("Day 2"), - Bar::new(2.5, 1.0).name("Day 3"), - Bar::new(3.5, 2.0).name("Day 4"), - Bar::new(4.5, 4.0).name("Day 5"), - ]) - .width(0.7) - .name("Set 1"); - - let chart2 = BarChart::new(vec![ - Bar::new(0.5, 1.0), - Bar::new(1.5, 1.5), - Bar::new(2.5, 0.1), - Bar::new(3.5, 0.7), - Bar::new(4.5, 0.8), - ]) - .width(0.7) - .name("Set 2") - .stack_on(&[&chart1]); - - let chart3 = BarChart::new(vec![ - Bar::new(0.5, -0.5), - Bar::new(1.5, 1.0), - Bar::new(2.5, 0.5), - Bar::new(3.5, -1.0), - Bar::new(4.5, 0.3), - ]) - .width(0.7) - .name("Set 3") - .stack_on(&[&chart1, &chart2]); - - let chart4 = BarChart::new(vec![ - Bar::new(0.5, 0.5), - Bar::new(1.5, 1.0), - Bar::new(2.5, 0.5), - Bar::new(3.5, -0.5), - Bar::new(4.5, -0.5), - ]) - .width(0.7) - .name("Set 4") - .stack_on(&[&chart1, &chart2, &chart3]); - - Plot::new("Stacked Bar Chart Demo") - .legend(Legend::default()) - .data_aspect(1.0) - .allow_drag(true) - .show(ui, |plot_ui| { - plot_ui.bar_chart(chart1); - plot_ui.bar_chart(chart2); - plot_ui.bar_chart(chart3); - plot_ui.bar_chart(chart4); - }) - .response - } - - fn box_plot(&self, ui: &mut Ui) -> Response { - let yellow = Color32::from_rgb(248, 252, 168); - let box1 = BoxPlot::new(vec![ - BoxElem::new(0.5, BoxSpread::new(1.5, 2.2, 2.5, 2.6, 3.1)).name("Day 1"), - BoxElem::new(2.5, BoxSpread::new(0.4, 1.0, 1.1, 1.4, 2.1)).name("Day 2"), - BoxElem::new(4.5, BoxSpread::new(1.7, 2.0, 2.2, 2.5, 2.9)).name("Day 3"), - ]) - .name("Experiment A"); - - let box2 = BoxPlot::new(vec![ - BoxElem::new(1.0, BoxSpread::new(0.2, 0.5, 1.0, 2.0, 2.7)).name("Day 1"), - BoxElem::new(3.0, BoxSpread::new(1.5, 1.7, 2.1, 2.9, 3.3)) - .name("Day 2: interesting") - .stroke(Stroke::new(1.5, yellow)) - .fill(yellow.linear_multiply(0.2)), - BoxElem::new(5.0, BoxSpread::new(1.3, 2.0, 2.3, 2.9, 4.0)).name("Day 3"), - ]) - .name("Experiment B"); - - let box3 = BoxPlot::new(vec![ - BoxElem::new(1.5, BoxSpread::new(2.1, 2.2, 2.6, 2.8, 3.0)).name("Day 1"), - BoxElem::new(3.5, BoxSpread::new(1.3, 1.5, 1.9, 2.2, 2.4)).name("Day 2"), - BoxElem::new(5.5, BoxSpread::new(0.2, 0.4, 1.0, 1.3, 1.5)).name("Day 3"), - ]) - .name("Experiment C"); - - Plot::new("Box Plot Demo") - .legend(Legend::default()) - .allow_zoom(true) - .allow_drag(true) - .show(ui, |plot_ui| { - plot_ui.box_plot(box1); - plot_ui.box_plot(box2); - plot_ui.box_plot(box3); - }) - .response - } -} - -// ---------------------------------------------------------------------------- diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..80f997b --- /dev/null +++ b/src/lib.rs @@ -0,0 +1,47 @@ +mod app3d; +mod controls; +mod simulator; + +use bevy::{asset::AssetMetaCheck, prelude::*}; + +pub struct YahsPlugin; + +impl Plugin for YahsPlugin { + fn build(&self, app: &mut App) { + // Add Bevy plugins. + app.add_plugins( + DefaultPlugins + .set(AssetPlugin { + // Wasm builds will check for meta files (that don't exist) + // if this isn't set. This causes errors and even panics on + // web build on itch. See + // https://github.com/bevyengine/bevy_github_ci_template/issues/48. + meta_check: AssetMetaCheck::Never, + ..default() + }) + .set(WindowPlugin { + primary_window: Window { + title: "🎈".to_string(), + canvas: Some("#bevy".to_string()), + fit_canvas_to_parent: true, + prevent_default_event_handling: true, + ..default() + } + .into(), + ..default() + }), + ); + + // Add the simulator plugins that don't deal with graphics. + app.add_plugins(( + simulator::SimulatorPlugins, + controls::ControlsPlugin, + )); + + // Add the 3D application plugins. + app.add_plugins(( + app3d::ScenePlugin, + app3d::InterfacePlugins, + )); + } +} diff --git a/src/main.rs b/src/main.rs index 976aacd..7d22685 100644 --- a/src/main.rs +++ b/src/main.rs @@ -1,15 +1,9 @@ -mod cli; -mod simulator; +// Disable console on Windows for non-dev builds. +#![cfg_attr(not(feature = "dev"), windows_subsystem = "windows")] -#[cfg(feature = "gui")] -mod gui; +use bevy::prelude::*; +use yahs::YahsPlugin; -fn main() { - // look for the RUST_LOG env var or default to "info" - let rust_log = std::env::var("RUST_LOG").unwrap_or_else(|_| "info".to_owned()); - std::env::set_var("RUST_LOG", rust_log); - // initialize pretty print logger - pretty_env_logger::init(); - // parse the commands, arguments, and options - cli::parse_inputs(); +fn main() -> AppExit { + App::new().add_plugins(YahsPlugin).run() } diff --git a/src/simulator/atmosphere.rs b/src/simulator/atmosphere.rs new file mode 100644 index 0000000..bc65af3 --- /dev/null +++ b/src/simulator/atmosphere.rs @@ -0,0 +1,123 @@ +//! Atmosphere model based on the US Standard Atmosphere, 1976. +//! +//! Reference: +//! - https://apps.dtic.mil/dtic/tr/fulltext/u2/a035728.pdf +//! - https://www.translatorscafe.com/unit-converter/en-US/calculator/altitude +//! - https://www.grc.nasa.gov/WWW/K-12/airplane/atmosmet.html + +use bevy::prelude::*; + +use super::{ + ideal_gas::{ideal_gas_density, GasSpecies}, + Density, Position, Pressure, SimState, SimulatedBody, Temperature, +}; + +pub struct AtmospherePlugin; +impl Plugin for AtmospherePlugin { + fn build(&self, app: &mut App) { + app.insert_resource(Atmosphere); + app.add_systems(Update, fault_if_out_of_bounds); + } +} + +/// US Standard Atmosphere, 1976 +#[derive(Resource)] +pub struct Atmosphere; + +impl Atmosphere { + pub const MAX_ALTITUDE: f32 = 84999.0; // small margin to avoid panics + pub const MIN_ALTITUDE: f32 = -56.0; // small margin to avoid panics + + pub fn out_of_bounds(&self, position: Vec3) -> bool { + match position.y { + y if (y > Atmosphere::MAX_ALTITUDE) => true, + y if (y < Atmosphere::MIN_ALTITUDE) => true, + _ => false, + } + } + + /// Temperature (K) of the atmosphere at a position. + pub fn temperature(&self, position: Vec3) -> Temperature { + // TODO: Look up temperature based on latitude, longitude, not just altitude + coesa_temperature(position.y).unwrap() // we should handle this better + } + + /// Pressure (Pa) of the atmosphere at a position. + pub fn pressure(&self, position: Vec3) -> Pressure { + // TODO: Look up pressure based on latitude, longitude, not just altitude + coesa_pressure(position.y).unwrap() // we should handle this better + } + + /// Density (kg/m³) of the atmosphere at a position. + pub fn density(&self, position: Vec3) -> Density { + ideal_gas_density( + self.temperature(position), + self.pressure(position), + &GasSpecies::air(), + ) + } +} + +/// If any of the simulated bodies are out of bounds, set the app state to anomaly +/// TODO: we should use an event for this +fn fault_if_out_of_bounds( + atmosphere: Res, + bodies: Query<(Entity, &Position), With>, + mut next_state: ResMut>, +) { + for (_, position) in bodies.iter() { + if atmosphere.out_of_bounds(position.0) { + next_state.set(SimState::Anomaly) + }; + } +} + +/// Temperature (K) of the atmosphere at a given altitude (m). +/// Only valid for altitudes below 85,000 meters. +/// Based on the US Standard Atmosphere, 1976. (aka COESA) +fn coesa_temperature(altitude: f32) -> Result { + if (-57.0..11000.0).contains(&altitude) { + Ok(Temperature::from_celsius(15.04 - 0.00649 * altitude)) + } else if (11000.0..25000.0).contains(&altitude) { + Ok(Temperature::from_celsius(-56.46)) + } else if (25000.0..85000.0).contains(&altitude) { + Ok(Temperature::from_celsius(-131.21 + 0.00299 * altitude)) + } else { + Err(format!( + "Altitude {:}m is outside of the accepted range! Must be 0-85000m", + altitude + )) + } +} + +/// Pressure (Pa) of the atmosphere at a given altitude (m). +/// Only valid for altitudes below 85,000 meters. +/// Based on the US Standard Atmosphere, 1976. (aka COESA) +fn coesa_pressure(altitude: f32) -> Result { + if (-57.0..11000.0).contains(&altitude) { + Ok(Pressure::from_kilopascal( + 101.29 + * f32::powf( + coesa_temperature(altitude).unwrap_or_default().kelvin() / 288.08, + 5.256, + ), + )) + } else if (11000.0..25000.0).contains(&altitude) { + Ok(Pressure::from_kilopascal( + 22.65 * f32::exp(1.73 - 0.000157 * altitude), + )) + } else if (25000.0..85000.0).contains(&altitude) { + Ok(Pressure::from_kilopascal( + 2.488 + * f32::powf( + coesa_temperature(altitude).unwrap_or_default().kelvin() / 216.6, + -11.388, + ), + )) + } else { + Err(format!( + "Altitude {:}m is outside of the accepted range! Must be 0-85000m", + altitude + )) + } +} diff --git a/src/simulator/balloon.rs b/src/simulator/balloon.rs index b5c211d..b94ab97 100644 --- a/src/simulator/balloon.rs +++ b/src/simulator/balloon.rs @@ -1,243 +1,41 @@ -// ---------------------------------------------------------------------------- -// Balloon -// ------- -// Properties, attributes and functions related to the balloon. -// ---------------------------------------------------------------------------- +//! Properties, attributes and functions related to the balloon. -#![allow(dead_code)] +use avian3d::prelude::*; +use bevy::prelude::*; +#[cfg(feature = "config-files")] +use serde::{Deserialize, Serialize}; -extern crate libm; - -use log::debug; -use serde::Deserialize; -use std::f32::consts::PI; -use std::fmt; - -use super::{gas, SolidBody}; - -#[derive(Copy, Clone)] -pub struct Balloon { - pub intact: bool, // whether or not it has burst - pub lift_gas: gas::GasVolume, // gas inside the balloon - pub material: Material, // what the balloon is made of - pub skin_thickness: f32, // thickness of the skin of the balloon (m) - mass: f32, // balloon mass (kg) - drag_coeff: f32, // drag coefficient - unstretched_thickness: f32, // thickness of the skin of the balloon without stretch (m) - unstretched_radius: f32, // radius of balloon without stretch (m) - pub temperature: f32, // fail if surface temperature exceeds this (K) - stress: f32, - strain: f32, -} - -impl Balloon { - pub fn new( - material: Material, // material of balloon skin - skin_thickness: f32, // balloon skin thickness (m) at zero pressure - barely_inflated_diameter: f32, // internal diameter (m) at zero pressure - lift_gas: gas::GasVolume, // species of gas inside balloon - ) -> Self { - let unstretched_radius = barely_inflated_diameter / 2.0; - let mass = shell_volume(unstretched_radius, skin_thickness) * material.density; - Balloon { - intact: true, - mass, - temperature: 293.0, - drag_coeff: 0.3, - lift_gas, - material, - skin_thickness, - unstretched_thickness: skin_thickness, - unstretched_radius, - stress: 0.0, - strain: 1.0, - } - } - - pub fn surface_area(self) -> f32 { - sphere_surface_area(sphere_radius_from_volume(self.lift_gas.volume())) - } - - pub fn radius(self) -> f32 { - sphere_radius_from_volume(self.volume()) - } - - pub fn volume(self) -> f32 { - self.lift_gas.volume() - } - - fn set_volume(&mut self, new_volume: f32) { - self.lift_gas.set_volume(new_volume) - } - - pub fn pressure(&mut self) -> f32 { - self.lift_gas.pressure() - } - - fn set_pressure(&mut self, new_pressure: f32) { - self.lift_gas.set_pressure(new_pressure) - } - - fn set_thickness(&mut self, new_thickness: f32) { - self.skin_thickness = new_thickness - } - - fn gage_pressure(self, external_pressure: f32) -> f32 { - self.lift_gas.pressure() - external_pressure - } - - pub fn stress(self) -> f32 { - self.stress - } - - fn set_stress(&mut self, external_pressure: f32) { - // hoop stress (Pa) of thin-walled hollow sphere from internal pressure - // https://en.wikipedia.org/wiki/Pressure_vessel#Stress_in_thin-walled_pressure_vessels - // https://pkel015.connect.amazon.auckland.ac.nz/SolidMechanicsBooks/Part_I/BookSM_Part_I/07_ElasticityApplications/07_Elasticity_Applications_03_Presure_Vessels.pdf - self.stress = self.gage_pressure(external_pressure) * self.radius() / (2.0 * self.skin_thickness); - if self.stress > self.material.max_stress { - self.burst(format!( - "Hoop stress ({:?} Pa) exceeded maximum stress ({:?} Pa)", - self.stress, self.material.max_stress - )); - } - } - - pub fn strain(self) -> f32 { - self.strain - } - - fn set_strain(&mut self) { - // strain (%) of thin-walled hollow sphere from internal pressure - // https://en.wikipedia.org/wiki/Pressure_vessel#Stress_in_thin-walled_pressure_vessels - // https://pkel015.connect.amazon.auckland.ac.nz/SolidMechanicsBooks/Part_I/BookSM_Part_I/07_ElasticityApplications/07_Elasticity_Applications_03_Presure_Vessels.pdf - self.strain = self.radius() / self.unstretched_radius; - if self.strain > self.material.max_strain { - self.burst(format!( - "Tangential strain ({:?} %) exceeded maximum strain ({:?} %)", - self.strain * 100.0, - self.material.max_strain * 100.0 - )); - } - } - - fn radial_displacement(self, external_pressure: f32) -> f32 { - // https://pkel015.connect.amazon.auckland.ac.nz/SolidMechanicsBooks/Part_I/BookSM_Part_I/07_ElasticityApplications/07_Elasticity_Applications_03_Presure_Vessels.pdf - ((1.0 - self.material.poissons_ratio) / self.material.elasticity) - * ((self.gage_pressure(external_pressure) * libm::powf(self.radius(), 2.0)) / 2.0 - * self.skin_thickness) - } - - fn rebound(&mut self, radial_displacement: f32) -> f32 { - // https://physics.stackexchange.com/questions/10372/inflating-a-balloon-expansion-resistance - self.set_thickness( - self.unstretched_thickness * libm::powf(self.unstretched_radius / self.radius(), 2.0), - ); - 2.0 * self.material.elasticity - * radial_displacement - * self.unstretched_thickness - * self.unstretched_radius - / libm::powf(self.radius(), 3.0) - } - - pub fn stretch(&mut self, external_pressure: f32) { - // stretch the balloon and/or compress the gas inside. - // - the gas wants to be at the same pressure as ambient - // - the balloon will stretch in response to the pressure difference - // - the balloon will likely not stretch enough to reach equilibrium - // - the difference between the ideal gas volume and the deformed - // balloon volume is the new pressure difference - // - the balloon fails when it starts to plasticly deform, in other - // words the balloon stretches as long as tangential stress is less - // than the material's yield stress - debug!( - "current gage pressure: {:?}", - self.gage_pressure(external_pressure) - ); - - self.set_stress(external_pressure); - self.set_strain(); - - if self.intact { - let delta_r = self.radial_displacement(external_pressure); - debug!( - "radius before stretch: {:?} delta_r: {:?}", - self.radius(), - delta_r - ); - let internal_pressure = self.rebound(delta_r); - self.set_pressure(internal_pressure + external_pressure); - debug!("radius after stretch: {:?}", self.radius()); - debug!( - "gage pressure after stretch: {:?}", - self.gage_pressure(external_pressure) - ); - } - - } - - fn burst(&mut self, reason: String) { - // Assert new balloon attributes to reflect that it has burst - self.intact = false; - self.set_volume(0.0); - self.lift_gas.set_mass(0.0); - log::warn!("The balloon has burst! Reason: {:?}", reason) - } -} +use super::{ + SimulatedBody, + ideal_gas::{GasSpecies, IdealGasBundle}, + properties::*, +}; -impl SolidBody for Balloon { - fn drag_area(&self) -> f32 { - if self.intact { - sphere_radius_from_volume(self.volume()) - } else { - 0.0 - } - } +pub struct BalloonPlugin; - fn drag_coeff(&self) -> f32 { - self.drag_coeff - } +impl Plugin for BalloonPlugin { + fn build(&self, app: &mut App) { + app.add_systems(Startup, spawn_balloon); - fn total_mass(&self) -> f32 { - self.mass + // Register types for reflection + app.register_type::(); + app.register_type::(); } } -fn sphere_volume(radius: f32) -> f32 { - (4.0 / 3.0) * PI * libm::powf(radius, 3.0) +#[derive(Bundle)] +pub struct BalloonBundle { + pub balloon: Balloon, + pub gas: IdealGasBundle, + pub pbr: PbrBundle, } -fn shell_volume(internal_radius: f32, thickness: f32) -> f32 { - let external_radius = internal_radius + thickness; - let internal_volume = sphere_volume(internal_radius); - let external_volume = sphere_volume(external_radius); - external_volume - internal_volume -} - -fn sphere_radius_from_volume(volume: f32) -> f32 { - libm::powf(volume, 1.0 / 3.0) / (4.0 / 3.0) * PI -} - -fn sphere_surface_area(radius: f32) -> f32 { - 4.0 * PI * libm::powf(radius, 2.0) -} - -pub fn projected_spherical_area(volume: f32) -> f32 { - // Get the projected area (m^2) of a sphere with a given volume (m^3) - libm::powf(sphere_radius_from_volume(volume), 2.0) * PI -} - -// ---------------------------------------------------------------------------- -// Materials -// --------- -// Properties, attributes and functions related to the material properties. -// Source: https://www.matweb.com/ -// ---------------------------------------------------------------------------- - -#[derive(Copy, Clone, PartialEq)] -pub struct Material { +#[derive(Debug, Clone, PartialEq, Reflect)] +#[cfg_attr(feature = "config-files", derive(Serialize, Deserialize))] +pub struct BalloonMaterial { + pub name: String, pub max_temperature: f32, // temperature (K) where the given material fails - pub density: f32, // density (kg/m^3) + pub density: f32, // density (kg/m³) pub emissivity: f32, // how much thermal radiation is emitted pub absorptivity: f32, // how much thermal radiation is absorbed pub thermal_conductivity: f32, // thermal conductivity (W/mK) of the material at room temperature @@ -248,77 +46,167 @@ pub struct Material { pub max_stress: f32, // tangential stress at failure (Pa) } -impl Material { - pub fn new(material_type: MaterialType) -> Self { - match material_type { - MaterialType::Rubber => RUBBER, - MaterialType::LDPE | MaterialType::LowDensityPolyethylene => LOW_DENSITY_POLYETHYLENE, - _ => NOTHING, +impl Default for BalloonMaterial { + fn default() -> Self { + BalloonMaterial { + name: "Latex".to_string(), + max_temperature: 373.0, // Example value in Kelvin + density: 920.0, // Example density in kg/m³ + emissivity: 0.9, // Example emissivity + absorptivity: 0.9, // Example absorptivity + thermal_conductivity: 0.13, // Example thermal conductivity in W/mK + specific_heat: 2000.0, // Example specific heat in J/kgK + poissons_ratio: 0.5, // Example Poisson's ratio + elasticity: 0.01e9, // Example Young's Modulus in Pa + max_strain: 0.8, // Example max strain (unitless) + max_stress: 0.5e6, // Example max stress in Pa } } } -#[derive(Copy, Clone, PartialEq, Deserialize)] -pub enum MaterialType { - // Species of gas with a known molar mass (kg/mol) - Nothing, - Rubber, - LDPE, - LowDensityPolyethylene, +/// Balloon properties. The balloon always conforms to the surface of a +/// collider. It does not have its own rigid body. +#[derive(Component, Reflect)] +#[cfg_attr(feature = "config-files", derive(Serialize, Deserialize))] +pub struct Balloon { + /// Balloon material type + pub skin_material: BalloonMaterial, + /// Thickness of balloon membrane in meters. For use in calculating stress. + pub unstretched_thickness: f32, + /// surface area of balloon without stretch (m²). For use in calculating stress. + pub unstretched_area: f32, } -impl fmt::Display for MaterialType { - fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - match *self { - MaterialType::Nothing => write!(f, "nothing"), - MaterialType::Rubber => write!(f, "rubber"), - MaterialType::LDPE | MaterialType::LowDensityPolyethylene => { - write!(f, "low-density polyethylene (LDPE)") - } +impl Default for Balloon { + fn default() -> Self { + Balloon { + skin_material: BalloonMaterial::default(), + unstretched_thickness: 0.001, + unstretched_area: 4.0 * std::f32::consts::PI, } } } +fn spawn_balloon( + mut commands: Commands, + mut meshes: ResMut>, + mut materials: ResMut>, +) { + let radius = 0.3; + commands.spawn(( + Name::new("BalloonBundle"), + SimulatedBody, + BalloonBundle { + balloon: Balloon::default(), + gas: IdealGasBundle::new( + Collider::sphere(radius), + GasSpecies::helium(), + Temperature::STANDARD, + Pressure::STANDARD, + ), + pbr: PbrBundle { + mesh: meshes.add(Sphere::new(radius)), + material: materials.add(Color::srgb_u8(124, 144, 255)), + transform: Transform::from_xyz(0.0, 0.0, 0.0), + ..default() + }, + }, + RigidBody::Dynamic, + )); +} -pub const NOTHING: Material = Material { - // nothing - max_temperature: f32::INFINITY, - density: 0.0, - emissivity: 1.0, - absorptivity: 0.0, - thermal_conductivity: f32::INFINITY, - specific_heat: 0.0, - poissons_ratio: 0.5, - elasticity: f32::INFINITY, - max_strain: f32::INFINITY, - max_stress: f32::INFINITY, -}; - -pub const RUBBER: Material = Material { - // Nitrile Butadiene Rubber - // https://designerdata.nl/materials/plastics/rubbers/nitrile-butadiene-rubber - max_temperature: 385.0, - density: 1000.0, - emissivity: 0.86, - absorptivity: 0.86, - thermal_conductivity: 0.25, - specific_heat: 1490.0, - poissons_ratio: 0.5, - elasticity: 4_000_000.0, - max_strain: 8.0, - max_stress: 25_000_000.0, -}; - -pub const LOW_DENSITY_POLYETHYLENE: Material = Material { - // Low Density Polyethylene (LDPE) - // https://designerdata.nl/materials/plastics/thermo-plastics/low-density-polyethylene - max_temperature: 348.0, - density: 919.0, - emissivity: 0.94, - absorptivity: 0.94, - thermal_conductivity: 0.3175, - specific_heat: 2600.0, - poissons_ratio: 0.5, - elasticity: 300_000_000.0, - max_strain: 6.25, - max_stress: 10_000_000.0, -}; +// impl Balloon { + +// pub fn gage_pressure(&self, external_pressure: f32) -> f32 { +// self.lift_gas.pressure() - external_pressure +// } + +// fn set_stress(&mut self, external_pressure: f32) { +// // hoop stress (Pa) of thin-walled hollow sphere from internal pressure +// // https://en.wikipedia.org/wiki/Pressure_vessel#Stress_in_thin-walled_pressure_vessels +// // https://pkel015.connect.amazon.auckland.ac.nz/SolidMechanicsBooks/Part_I/BookSM_Part_I/07_ElasticityApplications/07_Elasticity_Applications_03_Presure_Vessels.pdf +// self.stress = +// self.gage_pressure(external_pressure) * self.radius() / (2.0 * self.skin_thickness); +// if self.stress > self.material.max_stress { +// self.burst(format!( +// "Hoop stress ({:?} Pa) exceeded maximum stress ({:?} Pa)", +// self.stress, self.material.max_stress +// )); +// } +// } + +// fn set_strain(&mut self) { +// // strain (%) of thin-walled hollow sphere from internal pressure +// // https://en.wikipedia.org/wiki/Pressure_vessel#Stress_in_thin-walled_pressure_vessels +// // https://pkel015.connect.amazon.auckland.ac.nz/SolidMechanicsBooks/Part_I/BookSM_Part_I/07_ElasticityApplications/07_Elasticity_Applications_03_Presure_Vessels.pdf +// self.strain = self.radius() / self.unstretched_radius; +// if self.strain > self.material.max_strain { +// self.burst(format!( +// "Tangential strain ({:?} %) exceeded maximum strain ({:?} %)", +// self.strain * 100.0, +// self.material.max_strain * 100.0 +// )); +// } +// } + +// pub fn radial_displacement(&self, external_pressure: f32) -> f32 { +// // https://pkel015.connect.amazon.auckland.ac.nz/SolidMechanicsBooks/Part_I/BookSM_Part_I/07_ElasticityApplications/07_Elasticity_Applications_03_Presure_Vessels.pdf +// ((1.0 - self.material.poissons_ratio) / self.material.elasticity) +// * ((self.gage_pressure(external_pressure) * f32::powf(self.radius(), 2.0)) / 2.0 +// * self.skin_thickness) +// } + +// fn rebound(&mut self, radial_displacement: f32) -> f32 { +// // https://physics.stackexchange.com/questions/10372/inflating-a-balloon-expansion-resistance +// self.set_thickness( +// self.unstretched_thickness * f32::powf(self.unstretched_radius / self.radius(), 2.0), +// ); +// 2.0 * self.material.elasticity +// * radial_displacement +// * self.unstretched_thickness +// * self.unstretched_radius +// / f32::powf(self.radius(), 3.0) +// } + +// pub fn stretch(&mut self, external_pressure: f32) { +// // stretch the balloon and/or compress the gas inside. +// // - the gas wants to be at the same pressure as ambient +// // - the balloon will stretch in response to the pressure difference +// // - the balloon will likely not stretch enough to reach equilibrium +// // - the difference between the ideal gas volume and the deformed +// // balloon volume is the new pressure difference +// // - the balloon fails when it starts to plasticly deform, in other +// // words the balloon stretches as long as tangential stress is less +// // than the material's yield stress +// debug!( +// "current gage pressure: {:?}", +// self.gage_pressure(external_pressure) +// ); + +// self.set_stress(external_pressure); +// self.set_strain(); + +// if self.intact { +// let delta_r = self.radial_displacement(external_pressure); +// debug!( +// "radius before stretch: {:?} delta_r: {:?}", +// self.radius(), +// delta_r +// ); +// let internal_pressure = self.rebound(delta_r); +// self.set_pressure(internal_pressure + external_pressure); +// debug!("radius after stretch: {:?}", self.radius()); +// debug!( +// "gage pressure after stretch: {:?}", +// self.gage_pressure(external_pressure) +// ); +// } +// } + +// fn burst(&mut self, reason: String) { +// // Assert new balloon attributes to reflect that it has burst +// self.intact = false; +// self.set_volume(0.0); +// self.lift_gas.set_mass(0.0); +// warn!("The balloon has burst! Reason: {:?}", reason) +// } +// } diff --git a/src/simulator/bus.rs b/src/simulator/bus.rs deleted file mode 100644 index dfa6631..0000000 --- a/src/simulator/bus.rs +++ /dev/null @@ -1,140 +0,0 @@ -// ---------------------------------------------------------------------------- -// Bus -// --- -// Properties and functions of the balloon's structure and avionics buses -// ---------------------------------------------------------------------------- - -use super::{ - config::{BodyConfig, ParachuteConfig}, - forces, - gas::Atmosphere, - SolidBody, -}; - - -#[derive(Copy, Clone)] -pub struct Body { - mass: f32, // kg - drag_coeff: f32, // drag coefficient - drag_area: f32, // effective area used for drag, m^2 -} - -impl Body { - pub fn new(config: BodyConfig) -> Self { - Self { - mass: config.mass_kg, - drag_coeff: config.drag_coeff, - drag_area: config.drag_area_m2, - } - } -} - -impl SolidBody for Body { - fn drag_area(&self) -> f32 { - self.drag_area - } - fn drag_coeff(&self) -> f32 { - self.drag_coeff - } - fn total_mass(&self) -> f32 { - self.mass - } -} - - -#[derive(Copy, Clone)] -enum ParachuteState { - Stowed, - PartiallyOpen, - FullyOpen, -} - - -#[derive(Copy, Clone)] -pub struct Parachute { - state: ParachuteState, - mass: f32, //kg - drag_coeff: f32, // drag coefficient - deploy_force: f32, // force needed to go from partial to full open, N - drag_area: f32, // effective area used for drag, m^2 - deploy_progress: f32, // percentage of how open the parachute is - deploy_rate: f32, // progress increment per time-step -} - -impl SolidBody for Parachute { - fn drag_area(&self) -> f32 { - match self.state { - ParachuteState::Stowed => 0.0, - ParachuteState::PartiallyOpen => self.drag_area * self.deploy_progress, - ParachuteState::FullyOpen => self.drag_area, - } - } - fn drag_coeff(&self) -> f32 { - self.drag_coeff - } - fn total_mass(&self) -> f32 { - self.mass - } -} - -impl Parachute { - fn deploy(&mut self, force: f32) { - match self.state { - ParachuteState::Stowed => { - if force >= self.deploy_force { - self.state = ParachuteState::PartiallyOpen - } - } - ParachuteState::PartiallyOpen => self.continue_deploying(), - ParachuteState::FullyOpen => {} - } - } - - fn continue_deploying(&mut self) { - self.deploy_progress += &self.deploy_rate; - if self.deploy_progress >= 1.0 { - self.state = ParachuteState::FullyOpen - } - } -} - - -#[derive(Copy, Clone)] -pub struct ParachuteSystem { - pub main: Parachute, - pub drogue: Parachute, -} - -impl ParachuteSystem { - pub fn new(config: ParachuteConfig, timestep_size: f32) -> Self { - let main = Parachute { - mass: config.total_mass_kg, - drag_coeff: config.main_drag_coeff, - drag_area: config.main_area_m2, - deploy_force: config.deploy_force_n, - deploy_progress: 0.0, - deploy_rate: config.deploy_time_s / timestep_size, - state: ParachuteState::Stowed, - }; - let drogue = Parachute { - mass: 0.0, - drag_coeff: config.drogue_drag_coeff, - drag_area: config.drogue_area_m2, - deploy_force: 0.0, - deploy_progress: 0.0, - deploy_rate: 1.0, - state: ParachuteState::Stowed, - }; - Self { main, drogue } - } - - pub fn total_mass(&self) -> f32 { - self.main.mass + self.drogue.mass - } - - pub fn deploy(&mut self, atmo: Atmosphere, velocity: f32) { - self.drogue.deploy(100.0); - let drogue_drag = forces::drag(atmo, velocity, self.drogue); - self.main.deploy(drogue_drag); - } -} diff --git a/src/simulator/config.rs b/src/simulator/config.rs deleted file mode 100644 index 573b118..0000000 --- a/src/simulator/config.rs +++ /dev/null @@ -1,82 +0,0 @@ -use serde::Deserialize; - -use super::balloon::MaterialType; -use super::gas::GasSpecies; -use std::{fs, path::Path}; - -pub fn parse_from_file>(path: P) -> Config { - // Read the contents of the configuration file as string - let contents = match fs::read_to_string(path) { - // If successful return the files text as `contents`. - // `c` is a local variable. - Ok(c) => c, - // Handle the `error` case. - Err(_) => { - // Write `msg` to `stderr`. - eprintln!("Could not read file."); - // Exit the program with exit code `1`. - std::process::exit(1); - } - }; - parse(&contents) -} - -pub fn parse(contents: &String) -> Config { - // unpack the config TOML from string - toml::from_str(contents).unwrap() -} - -#[derive(Clone, Deserialize, PartialEq)] -pub struct Config { - pub environment: EnvConfig, - pub balloon: BalloonConfig, - pub bus: BusConfig, -} - -#[derive(Clone, Deserialize, PartialEq)] -pub struct EnvConfig { - pub real_time: bool, - pub tick_rate_hz: f32, - pub max_elapsed_time_s: f32, - pub initial_altitude_m: f32, - pub initial_velocity_m_s: f32, -} - -#[derive(Clone, Deserialize, PartialEq)] -pub struct BalloonConfig { - pub material: MaterialType, // balloon material - pub thickness_m: f32, // thickness of balloon membrane - pub barely_inflated_diameter_m: f32, // assuming balloon is a sphere, diameter of "unstressed" balloon membrane when filled - pub lift_gas: GasConfig, -} - -#[derive(Clone, Deserialize, PartialEq)] -pub struct GasConfig { - pub species: GasSpecies, - pub mass_kg: f32, -} - -#[derive(Clone, Deserialize, PartialEq)] -pub struct BusConfig { - pub body: BodyConfig, - pub parachute: ParachuteConfig, - // pub control: ControlConfig, -} - -#[derive(Copy, Clone, Deserialize, PartialEq)] -pub struct BodyConfig { - pub mass_kg: f32, // mass of all components less ballast material - pub drag_area_m2: f32, // effective area used for drag calculations during freefall - pub drag_coeff: f32, // drag coefficient of the payload during freefall -} - -#[derive(Copy, Clone, Deserialize, PartialEq)] -pub struct ParachuteConfig { - pub total_mass_kg: f32, // mass of the parachute system (main + drogue) - pub drogue_area_m2: f32, // drogue parachute effective area used for drag calculations - pub drogue_drag_coeff: f32, // drogue parachute drag coefficient - pub main_area_m2: f32, // main parachute effective area used for drag calculations - pub main_drag_coeff: f32, // main parachute drag coefficient when fully deployed - pub deploy_force_n: f32, // force needed for the drogue to deploy the main chute - pub deploy_time_s: f32, // how long the main chute stays in the partially open state -} diff --git a/src/simulator/constants.rs b/src/simulator/constants.rs deleted file mode 100644 index 1138fdc..0000000 --- a/src/simulator/constants.rs +++ /dev/null @@ -1,14 +0,0 @@ -// ---------------------------------------------------------------------------- -// Cosntants -// --------- -// Physical constants (read-only). -// ---------------------------------------------------------------------------- -pub const STANDARD_TEMPERATURE: f32 = 273.15; // [K] -pub const STANDARD_PRESSURE: f32 = 101325.0; // [Pa] -pub const BOLTZMANN_CONSTANT: f32 = 1.38e-23_f32; // [J/K] -pub const AVOGADRO_CONSTANT: f32 = 6.022e+23_f32; // [1/mol] -pub const R: f32 = BOLTZMANN_CONSTANT * AVOGADRO_CONSTANT; // [J/K-mol] Ideal gas constant - -pub const STANDARD_G: f32 = 9.80665; // [m/s^2] standard gravitational acceleration -pub const EARTH_RADIUS_M: f32 = 6371007.2; // [m] mean radius of Earth - diff --git a/src/simulator/forces.rs b/src/simulator/forces.rs deleted file mode 100644 index b39f390..0000000 --- a/src/simulator/forces.rs +++ /dev/null @@ -1,53 +0,0 @@ -// ---------------------------------------------------------------------------- -// Forces -// ------- -// Forces that act in the vertical axis. All forces assume a positive-up -// coordinate frame and aR_E R_Eported in Newtons. -// ---------------------------------------------------------------------------- - -#![allow(dead_code)] - -extern crate libm; - -use super::constants::{EARTH_RADIUS_M, STANDARD_G}; -use super::{gas, SolidBody}; - -fn g(altitude: f32) -> f32 { - // Acceleration (m/s^2) from gravity at an altitude (m) above mean sea level. - -STANDARD_G * (EARTH_RADIUS_M / (EARTH_RADIUS_M + altitude)) // [m/s^2] -} - -pub fn weight(altitude: f32, mass: f32) -> f32 { - // Weight (N) as a function of altitude (m) and mass (kg). - g(altitude) * mass // [N] -} - -pub fn buoyancy(altitude: f32, atmo: gas::Atmosphere, lift_gas: gas::GasVolume) -> f32 { - // Force (N) due to air displaced by the given gas volume. - let v = lift_gas.volume(); - if v > 0.0 { - let rho_atmo = atmo.density(); - let rho_lift = lift_gas.density(); - return lift_gas.volume() * (rho_lift - rho_atmo) * g(altitude) - } else { - return 0.0 - } -} - -pub fn drag(atmo: gas::Atmosphere, velocity: f32, body: T) -> f32 { - // Force (N) due to drag against the balloon - let direction = -libm::copysignf(1.0, velocity); - direction * body.drag_coeff() / 2.0 * atmo.density() * libm::powf(velocity, 2.0) * body.drag_area() -} - -pub fn gross_lift(atmo: gas::Atmosphere, lift_gas: gas::GasVolume) -> f32 { - // [kg] - let rho_atmo = atmo.density(); - let rho_lift = lift_gas.density(); - lift_gas.volume() * (rho_lift - rho_atmo) -} - -pub fn free_lift(atmo: gas::Atmosphere, lift_gas: gas::GasVolume, total_dry_mass: f32) -> f32 { - // [kg] - gross_lift(atmo, lift_gas) - total_dry_mass -} diff --git a/src/simulator/forces/aero.rs b/src/simulator/forces/aero.rs new file mode 100644 index 0000000..9c5497f --- /dev/null +++ b/src/simulator/forces/aero.rs @@ -0,0 +1,106 @@ +//! Forces applied to rigid bodies due to aerodynamic drag. + +use avian3d::{math::PI, prelude::*}; +use parry3d::shape::{ShapeType, Shape, Ball}; +use bevy::prelude::*; +use bevy_trait_query::{self, RegisterExt}; + +use super::{Atmosphere, Density, ForceUpdateOrder, Force, SimulatedBody}; + +pub struct AeroForcesPlugin; + +impl Plugin for AeroForcesPlugin { + fn build(&self, app: &mut App) { + app.register_type::(); + app.register_component_as::(); + app.add_systems(Update, update_drag_parameters.in_set(ForceUpdateOrder::Prepare)); + } +} + +/// Force (N) due to drag as a solid body moves through a fluid. +#[derive(Component, Reflect)] +pub struct Drag { + flow_velocity: Vec3, + ambient_density: Density, + drag_area: f32, + drag_coeff: f32, +} +impl Default for Drag { + fn default() -> Self { + Self { + flow_velocity: Vec3::ZERO, + ambient_density: Density::ZERO, + drag_area: 0.0, + drag_coeff: 1.0, + } + } +} +impl Drag { + pub fn update( + &mut self, + flow_velocity: Vec3, + ambient_density: Density, + drag_area: f32, + drag_coeff: f32, + ) { + self.flow_velocity = flow_velocity; + self.ambient_density = ambient_density; + self.drag_area = drag_area; + self.drag_coeff = drag_coeff; + } +} +impl Force for Drag { + fn name(&self) -> String { + String::from("Drag") + } + fn force(&self) -> Vec3 { + drag( + self.flow_velocity, + self.ambient_density.kg_per_m3(), + self.drag_area, + self.drag_coeff, + ) + } + fn point_of_application(&self) -> Vec3 { + Vec3::ZERO + } +} + +fn update_drag_parameters( + atmosphere: Res, + mut bodies: Query<(&mut Drag, &Position, &LinearVelocity, &Collider), With>, +) { + for (mut drag, position, velocity, collider) in bodies.iter_mut() { + let bounding_sphere = collider.shape().compute_bounding_sphere(&position.0.into()); + drag.update( + velocity.0, + atmosphere.density(position.0), + projected_spherical_area(bounding_sphere.radius()), + drag_coefficient(&Ball::new(bounding_sphere.radius()), &atmosphere), + ); + } +} + +/// Force (N) due to drag as a solid body moves through a fluid. +pub fn drag(velocity: Vec3, ambient_density: f32, drag_area: f32, drag_coeff: f32) -> Vec3 { + let drag_direction = -velocity.normalize_or_zero(); // oppose the object's velocity + let drag_magnitude = drag_coeff / 2.0 + * ambient_density + * velocity.length_squared() + * drag_area; + drag_direction * drag_magnitude +} + +/// Get the projected area (m^2) of a sphere with a given radius (m) +fn projected_spherical_area(radius: f32) -> f32 { + f32::powf(radius, 2.0) * PI +} + +/// Get the drag coefficient for a given shape and ambient conditions. +fn drag_coefficient(shape: &dyn Shape, _atmosphere: &Atmosphere) -> f32 { + match shape.shape_type() { + ShapeType::Ball => 1.17, + ShapeType::Cuboid => 2.05, + _ => 1.0, + } +} diff --git a/src/simulator/forces/body.rs b/src/simulator/forces/body.rs new file mode 100644 index 0000000..7cb0a0d --- /dev/null +++ b/src/simulator/forces/body.rs @@ -0,0 +1,129 @@ +//! Forces applied to rigid bodies due to gravity and buoyancy. + +use avian3d::prelude::*; +use bevy::prelude::*; +use bevy_trait_query::{self, RegisterExt}; + +use super::{Atmosphere, Density, Force, ForceUpdateOrder, Mass, SimulatedBody, Volume}; +use crate::simulator::properties::{EARTH_RADIUS_M, STANDARD_G}; + +pub struct BodyForcesPlugin; + +impl Plugin for BodyForcesPlugin { + fn build(&self, app: &mut App) { + app.register_type::(); + app.register_type::(); + + app.register_component_as::(); + app.register_component_as::(); + + app.add_systems( + Update, + (update_weight_parameters, update_buoyant_parameters).in_set(ForceUpdateOrder::Prepare), + ); + } +} + +/// Downward force (N) vector due to gravity as a function of altitude (m) and +/// mass (kg). The direction of this force is always world-space down. +#[derive(Component, Reflect)] +pub struct Weight { + position: Vec3, + mass: f32, +} +impl Default for Weight { + fn default() -> Self { + Self { + position: Vec3::ZERO, + mass: 0.0, + } + } +} +impl Weight { + pub fn update(&mut self, position: Vec3, mass: f32) { + self.position = position; + self.mass = mass; + } +} +impl Force for Weight { + fn name(&self) -> String { + String::from("Weight") + } + fn force(&self) -> Vec3 { + weight(self.position, self.mass) + } + fn point_of_application(&self) -> Vec3 { + self.position + } +} + +/// Force (N) from gravity at an altitude (m) above mean sea level. +fn g(position: Vec3) -> f32 { + let altitude = position.y; // [m] + STANDARD_G * (EARTH_RADIUS_M / (EARTH_RADIUS_M + altitude)) +} + +/// Downward force (N) vector due to gravity as a function of altitude (m) and +/// mass (kg). The direction of this force is always world-space down. +pub fn weight(position: Vec3, mass: f32) -> Vec3 { + Vec3::NEG_Y * g(position) * mass // [N] +} + +fn update_weight_parameters( + mut bodies: Query<(&mut Weight, &Position, &Mass), With>, +) { + for (mut weight, position, mass) in bodies.iter_mut() { + weight.update(position.0, mass.kg()); + } +} + +/// Upward force (N) vector due to atmosphere displaced by the given gas volume. +#[derive(Component, Reflect)] +pub struct Buoyancy { + position: Vec3, + displaced_volume: Volume, + ambient_density: Density, +} +impl Default for Buoyancy { + fn default() -> Self { + Self { + position: Vec3::ZERO, + displaced_volume: Volume::ZERO, + ambient_density: Density::ZERO, + } + } +} +impl Buoyancy { + pub fn update(&mut self, position: Vec3, displaced_volume: Volume, ambient_density: Density) { + self.position = position; + self.displaced_volume = displaced_volume; + self.ambient_density = ambient_density; + } +} +impl Force for Buoyancy { + fn name(&self) -> String { + String::from("Buoyancy") + } + fn force(&self) -> Vec3 { + buoyancy(self.position, self.displaced_volume, self.ambient_density) + } + fn point_of_application(&self) -> Vec3 { + self.position + } +} + +/// Upward force (N) vector due to atmosphere displaced by the given gas volume. +/// The direction of this force is always world-space up (it opposes gravity). +pub fn buoyancy(position: Vec3, displaced_volume: Volume, ambient_density: Density) -> Vec3 { + Vec3::Y * (displaced_volume.cubic_meters() * ambient_density.kg_per_m3() * g(position)) +} + +fn update_buoyant_parameters( + atmosphere: Res, + mut bodies: Query<(&mut Buoyancy, &Position, &Volume), With>, +) { + for (mut buoyancy, position, volume) in bodies.iter_mut() { + let density = atmosphere.density(position.0); + buoyancy.update(position.0, *volume, density); + } +} diff --git a/src/simulator/forces/mod.rs b/src/simulator/forces/mod.rs new file mode 100644 index 0000000..a5bb436 --- /dev/null +++ b/src/simulator/forces/mod.rs @@ -0,0 +1,111 @@ +//! Forces applied to rigid bodies. +pub mod aero; +pub mod body; + +use avian3d::prelude::*; +use bevy::prelude::*; +use bevy_trait_query; + +// Re-expert common forces +#[allow(unused_imports)] +pub use body::{Weight, Buoyancy}; +#[allow(unused_imports)] +pub use aero::Drag; + +use super::{Atmosphere, Density, Mass, Volume, SimulatedBody}; +pub struct ForcesPlugin; + +impl Plugin for ForcesPlugin { + fn build(&self, app: &mut App) { + // Disable the default forces since we apply our own. + app.insert_resource(Gravity(Vec3::ZERO)); + + // Update force vectors before solving physics. + app.configure_sets( + Update, + ( + ForceUpdateOrder::First, + ForceUpdateOrder::Prepare.after(ForceUpdateOrder::First), + ForceUpdateOrder::Apply + .after(ForceUpdateOrder::Prepare) + .before(PhysicsStepSet::First), + ), + ); + app.add_systems( + Update, + on_simulated_body_added.in_set(ForceUpdateOrder::First), + ); + app.add_systems( + Update, + update_total_external_force.in_set(ForceUpdateOrder::Apply), + ); + + app.add_plugins((aero::AeroForcesPlugin, body::BodyForcesPlugin)); + } +} + +#[derive(SystemSet, Debug, Clone, PartialEq, Eq, Hash)] +enum ForceUpdateOrder { + First, + Prepare, + Apply, +} + +/// A bundle of force components to be added to entities with a `RigidBody`. The +/// body can only have one of each type of force component. +#[derive(Bundle, Default)] +pub struct ForceBundle { + weight: body::Weight, + buoyancy: body::Buoyancy, + drag: aero::Drag, +} + +/// Add a `ForceBundle` to entities with a `RigidBody` when they are added. +fn on_simulated_body_added(mut commands: Commands, query: Query>) { + for entity in &query { + commands.entity(entity).insert(ForceBundle::default()); + } +} + +/// This trait is used to identify a force vector component. All forces are +/// collected and summed to determine the net force acting on a rigid body. All +/// forces assume a right-handed Y-up coordinate frame and are reported in +/// Newtons. +#[bevy_trait_query::queryable] +pub trait Force { + fn name(&self) -> String { + String::from("Force") + } + fn force(&self) -> Vec3; + fn direction(&self) -> Vec3 { + self.force().normalize() + } + fn magnitude(&self) -> f32 { + self.force().length() + } + fn point_of_application(&self) -> Vec3; +} + +/// Set the `ExternalForce` to the sum of all forces in the `Forces` collection. +/// This effectively applies all the calculated force vectors to the physics +/// rigid body without regard to where the forces came from. +/// +/// TODO: preserve the position of the total force vector and apply it at that +/// point instead of the center of mass. +fn update_total_external_force( + mut body_forces: Query<(&mut ExternalForce, &dyn Force, &RigidBody), With>, +) { + // Iterate over each entity that has force vector components. + for (mut physics_force_component, acting_forces, rigid_body) in body_forces.iter_mut() { + // Forces only act on dynamic bodies. Don't bother with other kinds. + if rigid_body.is_dynamic() { + let mut net_force = Vec3::ZERO; // reset the net force to zero + + // Iterate over each force vector component and compute its value. + for force in acting_forces.iter() { + net_force += force.force(); + } + physics_force_component.set_force(net_force); + } + } +} diff --git a/src/simulator/forces/scratch.rs b/src/simulator/forces/scratch.rs new file mode 100644 index 0000000..f9fdbae --- /dev/null +++ b/src/simulator/forces/scratch.rs @@ -0,0 +1,232 @@ + +use std::f64::consts::PI; +use nalgebra as na; + +// Constants +const AIR_DENSITY: f64 = 1.225; // kg/m³ at sea level, 15°C +const AIR_VISCOSITY: f64 = 1.81e-5; // kg/(m·s) at 15°C + +#[derive(Debug, Clone)] +struct Vertex(na::Point3); + +#[derive(Debug, Clone)] +struct Face { + vertices: Vec, // Indices into vertex array + normal: na::Vector3, + area: f64, +} + +#[derive(Debug)] +struct IrregularPolyhedron { + vertices: Vec, + faces: Vec, + center_of_mass: na::Point3, + orientation: na::Rotation3, +} + +impl IrregularPolyhedron { + fn new(vertices: Vec<[f64; 3]>, faces: Vec>) -> Self { + let vertices: Vec = vertices + .into_iter() + .map(|v| Vertex(na::Point3::new(v[0], v[1], v[2]))) + .collect(); + + let mut polyhedron_faces = Vec::new(); + for face_vertices in faces { + let v0 = vertices[face_vertices[0]].0; + let v1 = vertices[face_vertices[1]].0; + let v2 = vertices[face_vertices[2]].0; + + // Calculate face normal and area + let edge1 = v1 - v0; + let edge2 = v2 - v0; + let normal = edge1.cross(&edge2).normalize(); + + // Calculate area using Newell's method for arbitrary polygons + let mut area = 0.0; + for i in 0..face_vertices.len() { + let j = (i + 1) % face_vertices.len(); + let vi = vertices[face_vertices[i]].0; + let vj = vertices[face_vertices[j]].0; + area += (vj - vi).norm() * (vi + vj).norm() / 2.0; + } + + polyhedron_faces.push(Face { + vertices: face_vertices, + normal: normal, + area, + }); + } + + // Calculate center of mass (assuming uniform density) + let com = vertices.iter().fold(na::Point3::origin(), |acc, v| acc + v.0.coords) + / vertices.len() as f64; + + IrregularPolyhedron { + vertices, + faces: polyhedron_faces, + center_of_mass: com, + orientation: na::Rotation3::identity(), + } + } + + fn rotate(&mut self, axis: &na::Vector3, angle: f64) { + self.orientation *= na::Rotation3::from_axis_angle(&na::Unit::new_normalize(*axis), angle); + } + + fn get_projected_area(&self, flow_direction: &na::Vector3) -> f64 { + let flow_dir = self.orientation * flow_direction; + + // Calculate projected area by summing contributions from each face + self.faces.iter().map(|face| { + let cos_angle = face.normal.dot(&flow_dir).abs(); + face.area * cos_angle + }).sum() + } + + fn get_characteristic_length(&self) -> f64 { + // Use maximum distance between any two vertices as characteristic length + let mut max_distance = 0.0; + for (i, v1) in self.vertices.iter().enumerate() { + for v2 in self.vertices.iter().skip(i + 1) { + let distance = (v2.0 - v1.0).norm(); + max_distance = max_distance.max(distance); + } + } + max_distance + } + + fn get_surface_area(&self) -> f64 { + self.faces.iter().map(|face| face.area).sum() + } +} + +struct FlowConditions { + velocity: na::Vector3, + density: f64, + viscosity: f64, +} + +impl Default for FlowConditions { + fn default() -> Self { + FlowConditions { + velocity: na::Vector3::new(10.0, 0.0, 0.0), // 10 m/s in x direction + density: AIR_DENSITY, + viscosity: AIR_VISCOSITY, + } + } +} + +fn calculate_reynolds_number(length: f64, velocity: f64, conditions: &FlowConditions) -> f64 { + (conditions.density * velocity * length) / conditions.viscosity +} + +fn calculate_drag_coefficients(reynolds: f64, shape_complexity: f64) -> (f64, f64) { + // Pressure drag coefficient + let cd_pressure = match reynolds { + re if re < 1.0 => 24.0 / re, + re if re < 1000.0 => 24.0 / re + 6.0 / (1.0 + re.sqrt()), + re if re < 2e5 => 1.1 * shape_complexity, + re if re < 7e5 => 0.8 * shape_complexity, + _ => 0.6 * shape_complexity, + }; + + // Friction drag coefficient + let cd_friction = match reynolds { + re if re < 1e5 => 1.328 / re.sqrt(), + _ => 0.074 / reynolds.powf(0.2), + }; + + (cd_pressure, cd_friction) +} + +struct DragResult { + pressure_drag: na::Vector3, + friction_drag: na::Vector3, + reynolds_number: f64, + cd_pressure: f64, + cd_friction: f64, +} + +fn calculate_drag(polyhedron: &IrregularPolyhedron, conditions: &FlowConditions) -> DragResult { + let velocity_magnitude = conditions.velocity.norm(); + let flow_direction = conditions.velocity.normalize(); + + // Get characteristic length and Reynolds number + let char_length = polyhedron.get_characteristic_length(); + let reynolds = calculate_reynolds_number(char_length, velocity_magnitude, conditions); + + // Calculate shape complexity factor (ratio of actual surface area to minimum possible surface area) + let actual_surface_area = polyhedron.get_surface_area(); + let min_surface_area = 4.0 * PI * (char_length / 2.0).powi(2); // sphere surface area + let shape_complexity = actual_surface_area / min_surface_area; + + // Get drag coefficients + let (cd_pressure, cd_friction) = calculate_drag_coefficients(reynolds, shape_complexity); + + // Calculate projected area + let projected_area = polyhedron.get_projected_area(&flow_direction); + + // Calculate drag forces + let dynamic_pressure = 0.5 * conditions.density * velocity_magnitude.powi(2); + + let pressure_drag = flow_direction * (cd_pressure * dynamic_pressure * projected_area); + let friction_drag = flow_direction * (cd_friction * dynamic_pressure * actual_surface_area); + + DragResult { + pressure_drag, + friction_drag, + reynolds_number: reynolds, + cd_pressure, + cd_friction, + } +} + +fn main() { + // Example usage with an irregular tetrahedron + let vertices = vec![ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.3, 1.2, 0.0], + [0.5, 0.4, 1.5], + ]; + + let faces = vec![ + vec![0, 1, 2], + vec![0, 1, 3], + vec![1, 2, 3], + vec![0, 2, 3], + ]; + + let mut polyhedron = IrregularPolyhedron::new(vertices, faces); + + // Rotate the polyhedron 30 degrees around the y-axis + polyhedron.rotate(&na::Vector3::new(0.0, 1.0, 0.0), PI / 6.0); + + let conditions = FlowConditions::default(); + + let result = calculate_drag(&polyhedron, &conditions); + + println!("Irregular Polyhedron Analysis:"); + println!("Characteristic Length: {:.3} m", polyhedron.get_characteristic_length()); + println!("Surface Area: {:.3} m²", polyhedron.get_surface_area()); + println!("Projected Area: {:.3} m²", + polyhedron.get_projected_area(&conditions.velocity.normalize())); + + println!("\nFlow Conditions:"); + println!("Velocity: [{:.1}, {:.1}, {:.1}] m/s", + conditions.velocity.x, conditions.velocity.y, conditions.velocity.z); + + println!("\nResults:"); + println!("Reynolds Number: {:.1e}", result.reynolds_number); + println!("Pressure Drag Coefficient: {:.3}", result.cd_pressure); + println!("Friction Drag Coefficient: {:.3}", result.cd_friction); + println!("Pressure Drag: [{:.3}, {:.3}, {:.3}] N", + result.pressure_drag.x, result.pressure_drag.y, result.pressure_drag.z); + println!("Friction Drag: [{:.3}, {:.3}, {:.3}] N", + result.friction_drag.x, result.friction_drag.y, result.friction_drag.z); + println!("Total Drag: [{:.3}, {:.3}, {:.3}] N", + result.pressure_drag.x + result.friction_drag.x, + result.pressure_drag.y + result.friction_drag.y, + result.pressure_drag.z + result.friction_drag.z); +} diff --git a/src/simulator/gas.rs b/src/simulator/gas.rs deleted file mode 100644 index af8c7c9..0000000 --- a/src/simulator/gas.rs +++ /dev/null @@ -1,313 +0,0 @@ -// ---------------------------------------------------------------------------- -// Gas -// --- -// Ideal gas equations. This model assumes that temperature and pressure are -// known, and that volume is not constrained. As such, only a gas's species, -// mass, temperature, and pressure can be set explicitly; volume and density -// are read-only derived attributes. -// ---------------------------------------------------------------------------- -// Atmosphere -// ---------- -// Approximate atmospheric temperature and pressure as a function of altitude. -// Based on the US Standard Atmosphere, 1976. (aka COESA) -// Reference: -// https://apps.dtic.mil/dtic/tr/fulltext/u2/a035728.pdf -// https://www.translatorscafe.com/unit-converter/en-US/calculator/altitude -// https://www.grc.nasa.gov/WWW/K-12/airplane/atmosmet.html -// ---------------------------------------------------------------------------- - -#![allow(dead_code)] - -use super::constants::{R, STANDARD_PRESSURE, STANDARD_TEMPERATURE}; -use log::error; -use serde::Deserialize; -use std::fmt; - -#[derive(Copy, Clone, Deserialize, PartialEq)] -pub enum GasSpecies { - // Species of gas with a known molar mass (kg/mol) - Air, - He, - Helium, - H2, - Hydrogen, - N2, - Nitrogen, - O2, - Oxygen, - Ar, - Argon, - CO2, - CarbonDioxide, - Ne, - Neon, - Kr, - Krypton, - Xe, - Xenon, - CH4, - Methane, -} - -impl fmt::Display for GasSpecies { - fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - match *self { - GasSpecies::Air => write!(f, "Air"), - GasSpecies::He => write!(f, "Helium"), - GasSpecies::Helium => write!(f, "Helium"), - GasSpecies::H2 => write!(f, "Hydrogen"), - GasSpecies::Hydrogen => write!(f, "Hydrogen"), - GasSpecies::N2 => write!(f, "Nitrogen"), - GasSpecies::Nitrogen => write!(f, "Nitrogen"), - GasSpecies::O2 => write!(f, "Oxygen"), - GasSpecies::Oxygen => write!(f, "Oxygen"), - GasSpecies::Ar => write!(f, "Argon"), - GasSpecies::Argon => write!(f, "Argon"), - GasSpecies::CO2 => write!(f, "Carbon Dioxide"), - GasSpecies::CarbonDioxide => write!(f, "Carbon Dioxide"), - GasSpecies::Ne => write!(f, "Neon"), - GasSpecies::Neon => write!(f, "Neon"), - GasSpecies::Kr => write!(f, "Krypton"), - GasSpecies::Krypton => write!(f, "Krypton"), - GasSpecies::Xe => write!(f, "Xenon"), - GasSpecies::Xenon => write!(f, "Xenon"), - GasSpecies::CH4 => write!(f, "Methane"), - GasSpecies::Methane => write!(f, "Methane"), - } - } -} - -#[derive(Copy, Clone)] -pub struct GasVolume { - // A finite amount of a particular gas - species: GasSpecies, // type of gas - mass: f32, // [kg] amount of gas in the volume - temperature: f32, // [K] temperature - pressure: f32, // [Pa] pressure - molar_mass: f32, // [kg/mol] molar mass a.k.a. molecular weight - volume: f32, // [m^3] volume -} - -impl fmt::Display for GasVolume { - fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - write!( - f, - "{:}: {:} kg | {:} K | {:} Pa | {:} m^3", - self.species, self.mass, self.temperature, self.pressure, self.volume, - ) - } -} - -impl GasVolume { - pub fn new(species: GasSpecies, mass: f32) -> Self { - // Create a new gas volume as a finite amount of mass (kg) of a - // particular species of gas. Gas is initialized at standard - // temperature and pressure. - GasVolume { - species, - mass, - temperature: STANDARD_TEMPERATURE, - pressure: STANDARD_PRESSURE, - molar_mass: molar_mass(species), - volume: volume( - STANDARD_TEMPERATURE, - STANDARD_PRESSURE, - mass, - molar_mass(species), - ), - } - } - - pub fn update_from_ambient(&mut self, atmo: Atmosphere) { - self.temperature = atmo.temperature(); - self.pressure = atmo.pressure(); - } - - pub fn set_temperature(&mut self, new_temperature: f32) { - // set the temperature (K) of the GasVolume - self.temperature = new_temperature; - } - - pub fn set_pressure(&mut self, new_pressure: f32) { - // set the pressure (Pa) of the GasVolume - self.pressure = new_pressure; - } - - pub fn set_volume(&mut self, new_volume: f32) { - // set the volume (m^3) of the GasVolume and update the pressure - // according to the ideal gas law. - self.volume = new_volume; - let new_pressure = ((self.mass / self.molar_mass) * R * self.temperature) / self.volume; - self.set_pressure(new_pressure); - } - - pub fn set_mass(&mut self, new_mass: f32) { - // set the mass (kg) of the GasVolume - if new_mass >= 0.0 { - self.mass = new_mass; - } else { - error!("Cannot set mass below zero!") - } - } - - pub fn set_mass_from_volume(&mut self) { - // set the mass (kg) based on the current volume (m^3), - // density (kg/m^3), and molar mass (kg/mol) - self.mass = self.volume * (self.molar_mass / R) * (self.pressure / self.temperature); - } - - pub fn temperature(self) -> f32 { - // temperature (K) - self.temperature - } - - pub fn pressure(self) -> f32 { - // pressure (Pa) - self.pressure - } - - pub fn mass(self) -> f32 { - // mass (kg) - self.mass - } - pub fn density(self) -> f32 { - // density (kg/m^3) - density(self.temperature, self.pressure, self.molar_mass) - } - - pub fn volume(&self) -> f32 { - // volume (m^3) - volume(self.temperature, self.pressure, self.mass, self.molar_mass) - } -} - -#[derive(Copy, Clone)] -pub struct Atmosphere { - // US Standard Atmosphere, 1976 - altitude: f32, // [m] altitude (which determines the other attributes) - temperature: f32, // [K] temperature - pressure: f32, // [Pa] pressure - density: f32, // [kg/m^3] density - molar_mass: f32, // [kg/mol] molar mass a.k.a. molecular weight -} - -impl fmt::Display for Atmosphere { - fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - write!( - f, - "{:} K | {:} Pa | {:} kg/m^3", - self.temperature, self.pressure, self.density, - ) - } -} - -impl Atmosphere { - pub fn new(altitude: f32) -> Self { - Atmosphere { - altitude, - temperature: coesa_temperature(altitude), - pressure: coesa_pressure(altitude), - density: density( - coesa_temperature(altitude), - coesa_pressure(altitude), - molar_mass(GasSpecies::Air), - ), - molar_mass: molar_mass(GasSpecies::Air), - } - } - pub fn set_altitude(&mut self, new_altitude: f32) { - self.altitude = new_altitude; - // update all params - self.temperature = coesa_temperature(new_altitude); - self.pressure = coesa_pressure(new_altitude); - self.density = density(self.temperature, self.pressure, self.molar_mass); - } - - pub fn temperature(self) -> f32 { - // Temperature (K) - self.temperature - } - - pub fn pressure(self) -> f32 { - // Pressure (Pa) - self.pressure - } - - pub fn density(self) -> f32 { - // Density (kg/m^3) - self.density - } -} - -fn volume(temperature: f32, pressure: f32, mass: f32, molar_mass: f32) -> f32 { - // Volume (m^3) of an ideal gas from its temperature (K), pressure (Pa), - // mass (kg) and molar mass (kg/mol). - (mass / molar_mass) * R * temperature / pressure // [m^3] -} - -fn density(temperature: f32, pressure: f32, molar_mass: f32) -> f32 { - // Density (kg/m^3) of an ideal gas frorm its temperature (K), pressure (Pa), - // and molar mass (kg/mol) - (molar_mass * pressure) / (R * temperature) // [kg/m^3] -} - -fn molar_mass(species: GasSpecies) -> f32 { - // Get the molecular weight (kg/mol) of a dry gas at sea level. - // Source: US Standard Atmosphere, 1976 - match species { - GasSpecies::Air => 0.02897, - GasSpecies::He | GasSpecies::Helium => 0.0040026, - GasSpecies::H2 | GasSpecies::Hydrogen => 0.00201594, - GasSpecies::N2 | GasSpecies::Nitrogen => 0.0280134, - GasSpecies::O2 | GasSpecies::Oxygen => 0.0319988, - GasSpecies::Ar | GasSpecies::Argon => 0.039948, - GasSpecies::CO2 | GasSpecies::CarbonDioxide => 0.04400995, - GasSpecies::Ne | GasSpecies::Neon => 0.020183, - GasSpecies::Kr | GasSpecies::Krypton => 0.08380, - GasSpecies::Xe | GasSpecies::Xenon => 0.13130, - GasSpecies::CH4 | GasSpecies::Methane => 0.01604303, - } -} - -fn coesa_temperature(altitude: f32) -> f32 { - // Temperature (K) of the atmosphere at a given altitude (m). - // Only valid for altitudes below 85,000 meters. - // Based on the US Standard Atmosphere, 1976. (aka COESA) - if (-57.0..11000.0).contains(&altitude) { - celsius2kelvin(15.04 - 0.00649 * altitude) - } else if (11000.0..25000.0).contains(&altitude) { - celsius2kelvin(-56.46) - } else if (25000.0..85000.0).contains(&altitude) { - celsius2kelvin(-131.21 + 0.00299 * altitude) - } else { - error!( - "Altitude {:}m is outside of the accepted range! Must be 0-85000m", - altitude - ); - 0.0 - } -} - -fn coesa_pressure(altitude: f32) -> f32 { - // Pressure (Pa) of the atmosphere at a given altitude (m). - // Only valid for altitudes below 85,000 meters. - // Based on the US Standard Atmosphere, 1976. (aka COESA) - if (-57.0..11000.0).contains(&altitude) { - (101.29 * libm::powf(coesa_temperature(altitude) / 288.08, 5.256)) * 1000.0 - } else if (11000.0..25000.0).contains(&altitude) { - (22.65 * libm::expf(1.73 - 0.000157 * altitude)) * 1000.0 - } else if (25000.0..85000.0).contains(&altitude) { - (2.488 * libm::powf(coesa_temperature(altitude) / 216.6, -11.388)) * 1000.0 - } else { - error!( - "Altitude {:}m is outside of the accepted range! Must be 0-85000m", - altitude - ); - 0.0 - } -} - -fn celsius2kelvin(deg_celsius: f32) -> f32 { - // Convert degrees C to Kelvin - deg_celsius + 273.15 -} diff --git a/src/simulator/heat.rs b/src/simulator/heat.rs index f82b24b..91b7207 100644 --- a/src/simulator/heat.rs +++ b/src/simulator/heat.rs @@ -1,30 +1,24 @@ -// ---------------------------------------------------------------------------- -// Heat -// ---- -// Heat transferred between and stored in materials -// -// Assumptions: -// 1. convection heat transfer is linear -// This is reasonable for airflows with velocity under 5 m/s -// 2. boundary nodes do not change temperature -// 3. thermal resistances are constant between nodes -// -// https://www.sciencedirect.com/science/article/pii/S1000936118301018 -// https://materion.com/-/media/files/alloy/newsletters/technical-tidbits/issue-no-114-thermal-emissivity-and-radiative-heat-transfer.pdf -// ---------------------------------------------------------------------------- +//! Heat transferred between and stored in materials +//! +//! Assumptions: +//! - convection heat transfer is linear +//! This is reasonable for airflows with velocity under 5 m/s +//! - boundary nodes do not change temperature +//! - thermal resistances are constant between nodes +//! +//! References: +//! - https://www.sciencedirect.com/science/article/pii/S1000936118301018 +//! - https://materion.com/-/media/files/alloy/newsletters/technical-tidbits/issue-no-114-thermal-emissivity-and-radiative-heat-transfer.pdf #![allow(dead_code)] #![allow(unused_imports)] -use super::balloon::Balloon; -use super::gas::Atmosphere; - -fn absorbed(balloon: Balloon, heat_flux: f32) -> f32 { - // absorbed - // A * q'' - let surface_area: f32 = balloon.surface_area(); - surface_area * heat_flux -} +// fn absorbed(balloon: Balloon, heat_flux: f32) -> f32 { +// // absorbed +// // A * q'' +// let surface_area: f32 = balloon.surface_area(); +// surface_area * heat_flux +// } // convected // A * h * (T - Tatmo) @@ -33,4 +27,4 @@ fn absorbed(balloon: Balloon, heat_flux: f32) -> f32 { // eps * sigma * A * (T^4 - Tatmo^4) // stored -// rho * V * cp * dT/dt \ No newline at end of file +// rho * V * cp * dT/dt diff --git a/src/simulator/ideal_gas.rs b/src/simulator/ideal_gas.rs new file mode 100644 index 0000000..789d3f2 --- /dev/null +++ b/src/simulator/ideal_gas.rs @@ -0,0 +1,161 @@ +//! Ideal gas equations. +#![allow(dead_code)] + +use avian3d::prelude::*; +use bevy::prelude::*; +#[cfg(feature = "config-files")] +use serde::{Deserialize, Serialize}; + +use crate::simulator::properties::{Mass as SimMass, *}; + +pub const R: f32 = BOLTZMANN_CONSTANT * AVOGADRO_CONSTANT; // [J/K-mol] Ideal gas constant + +pub struct IdealGasPlugin; + +impl Plugin for IdealGasPlugin { + fn build(&self, app: &mut App) { + app.register_type::(); + app.add_systems(Update, ( + update_ideal_gas_volume_from_pressure, + update_ideal_gas_density_from_volume, + )); + } +} + +#[derive(Component, Debug, Clone, PartialEq, Reflect)] +#[cfg_attr(feature = "config-files", derive(Serialize, Deserialize))] +pub struct GasSpecies { + pub name: String, + pub abbreviation: String, + pub molar_mass: MolarMass, // [kg/mol] molar mass a.k.a. molecular weight +} + +impl GasSpecies { + /// Dry air. + pub fn air() -> Self { + GasSpecies { + name: "Air".to_string(), + abbreviation: "AIR".to_string(), + molar_mass: MolarMass(0.0289647), + } + } + + pub fn helium() -> Self { + GasSpecies { + name: "Helium".to_string(), + abbreviation: "He".to_string(), + molar_mass: MolarMass(0.0040026), + } + } +} + +impl Default for GasSpecies { + fn default() -> Self { + GasSpecies::air() + } +} + +impl GasSpecies { + pub fn new(name: String, abbreviation: String, molar_mass: MolarMass) -> Self { + GasSpecies { + name, + abbreviation, + molar_mass, + } + } +} + +/// Volume (m³) of an ideal gas from its temperature (K), pressure (Pa), +/// mass (kg) and molar mass (kg/mol). +pub fn ideal_gas_volume( + temperature: Temperature, + pressure: Pressure, + mass: SimMass, + species: &GasSpecies, +) -> Volume { + Volume( + (mass.kilograms() / species.molar_mass.kilograms_per_mole()) * R * temperature.kelvin() + / pressure.pascal(), + ) +} + +/// Density (kg/m³) of an ideal gas from its temperature (K), pressure (Pa), +/// and molar mass (kg/mol) +pub fn ideal_gas_density( + temperature: Temperature, + pressure: Pressure, + species: &GasSpecies, +) -> Density { + Density( + (species.molar_mass.kilograms_per_mole() * pressure.pascal()) / (R * temperature.kelvin()), + ) +} + +/// Gage pressure (Pa) of an ideal gas. This is the relative pressure compared +/// to the ambient pressure. Use `Atmosphere::pressure()` to get ambient +/// conditions. +pub fn gage_pressure(pressure: Pressure, ambient_pressure: Pressure) -> Pressure { + pressure - ambient_pressure +} + +/// A finite amount of a particular ideal gas +#[derive(Component, Debug)] +pub struct IdealGas; + +#[derive(Bundle, Debug)] +pub struct IdealGasBundle { + pub collider: Collider, + pub species: GasSpecies, + pub temperature: Temperature, + pub pressure: Pressure, + pub volume: Volume, + pub mass: SimMass, +} + +impl IdealGasBundle { + pub fn new( + collider: Collider, + species: GasSpecies, + temperature: Temperature, + pressure: Pressure, + ) -> Self { + let density = ideal_gas_density(temperature, pressure, &species); + let mass_props = collider.mass_properties(density.kg_per_m3()); + let mass = SimMass::from_mass_properties(mass_props); + Self { + collider, + species: species.clone(), + temperature, + pressure, + volume: ideal_gas_volume(temperature, pressure, mass, &species), + mass, + } + } +} + +impl Default for IdealGasBundle { + fn default() -> Self { + IdealGasBundle::new( + Collider::sphere(1.0), + GasSpecies::default(), + Temperature::STANDARD, + Pressure::STANDARD, + ) + } +} + +fn update_ideal_gas_volume_from_pressure( + mut query: Query<(&mut Volume, &Temperature, &Pressure, &SimMass, &GasSpecies), With>, +) { + for (mut volume, temperature, pressure, mass, species) in query.iter_mut() { + *volume = ideal_gas_volume(*temperature, *pressure, *mass, species); + } +} + +fn update_ideal_gas_density_from_volume( + mut query: Query<(&mut Density, &Volume, &SimMass), With>, +) { + for (mut density, volume, mass) in query.iter_mut() { + *density = *mass / *volume; + } +} diff --git a/src/simulator/io.rs b/src/simulator/io.rs deleted file mode 100644 index ab7ac5f..0000000 --- a/src/simulator/io.rs +++ /dev/null @@ -1,64 +0,0 @@ -use std::{fs, path::PathBuf}; -fn init_log_file(outpath: PathBuf) -> csv::Writer { - let mut writer = csv::Writer::from_path(outpath).unwrap(); - writer - .write_record(&[ - "time_s", - "altitude_m", - "ascent_rate_m_s", - "acceleration_m_s2", - "atmo_temp_K", - "atmo_pres_Pa", - "balloon_pres_Pa", - "balloon_radius_m", - "balloon_stress_Pa", - "balloon_strain_pct", - "balloon_thickness_m", - "drogue_parachute_area_m2", - "main_parachute_area_m2", - ]) - .unwrap(); - writer -} - -pub fn log_to_file(sim_output: &SimOutput, writer: &mut csv::Writer) { - writer - .write_record(&[ - sim_output.time_s.to_string(), - sim_output.altitude.to_string(), - sim_output.ascent_rate.to_string(), - sim_output.acceleration.to_string(), - sim_output.atmo_temp.to_string(), - sim_output.atmo_pres.to_string(), - sim_output.balloon_pres.to_string(), - sim_output.balloon_radius.to_string(), - sim_output.balloon_stress.to_string(), - sim_output.balloon_strain.to_string(), - sim_output.balloon_thickness.to_string(), - sim_output.drogue_parachute_area.to_string(), - sim_output.main_parachute_area.to_string(), - ]) - .unwrap(); - writer.flush().unwrap(); -} - -#[derive(Default, Copy, Clone)] -pub struct SimOutput { - pub time_s: f32, - pub altitude: f32, - pub ascent_rate: f32, - pub acceleration: f32, - pub atmo_temp: f32, - pub atmo_pres: f32, - pub balloon_pres: f32, - pub balloon_radius: f32, - pub balloon_stress: f32, - pub balloon_strain: f32, - pub balloon_thickness: f32, - pub drogue_parachute_area: f32, - pub main_parachute_area: f32, -} - -pub struct SimCommands { - //TODO: add ability to inject commands to logic controllers -} diff --git a/src/simulator/mod.rs b/src/simulator/mod.rs index ae1ea6b..c24cd24 100644 --- a/src/simulator/mod.rs +++ b/src/simulator/mod.rs @@ -1,15 +1,63 @@ -mod balloon; -mod bus; -pub mod config; -mod constants; -mod forces; -mod gas; -mod heat; -pub mod io; -pub mod schedule; - -pub trait SolidBody { - fn drag_area(&self) -> f32; - fn drag_coeff(&self) -> f32; - fn total_mass(&self) -> f32; +pub mod atmosphere; +pub mod balloon; +pub mod forces; +pub mod ideal_gas; +pub mod payload; +pub mod properties; + +use bevy::app::PluginGroupBuilder; +use bevy::prelude::*; +use avian3d::prelude::*; + +// Re-export the properties module at the top level. +pub use properties::{Temperature, Pressure, Volume, Density, Mass}; +pub use atmosphere::Atmosphere; + +/// A marker component for entities that are simulated. +#[derive(Component, Default)] +pub struct SimulatedBody; + + +pub struct SimulatorPlugins; + +impl PluginGroup for SimulatorPlugins { + fn build(self) -> PluginGroupBuilder { + PluginGroupBuilder::start::() + .add(CorePhysicsPlugin) + .add(atmosphere::AtmospherePlugin) + .add(balloon::BalloonPlugin) + } +} + +struct CorePhysicsPlugin; + +impl Plugin for CorePhysicsPlugin { + fn build(&self, app: &mut App) { + app.add_plugins(( + PhysicsPlugins::default(), + properties::CorePropertiesPlugin, + ideal_gas::IdealGasPlugin, + forces::ForcesPlugin, + )); + app.init_state::(); + app.add_systems(Update, pause_physics_time); + } +} + +#[derive(States, Debug, Default, Clone, Copy, Hash, PartialEq, Eq)] +pub enum SimState { + #[default] + Running, + Stopped, + Anomaly, +} + +fn pause_physics_time( + sim_state: Res>, + mut physics_time: ResMut>) { + match sim_state.as_ref().get() { + SimState::Running => physics_time.unpause(), + SimState::Stopped => physics_time.pause(), + SimState::Anomaly => physics_time.pause(), + } } diff --git a/src/simulator/payload.rs b/src/simulator/payload.rs new file mode 100644 index 0000000..3c107c5 --- /dev/null +++ b/src/simulator/payload.rs @@ -0,0 +1,21 @@ +//! The things that are carried by the balloon. +#![allow(dead_code)] + +use bevy::prelude::*; + +pub struct PayloadPlugin; + +impl Plugin for PayloadPlugin { + fn build(&self, _app: &mut App) { + // app.add_systems(Startup, spawn_payload); + } +} + +/// A thing carried by the balloon. +#[derive(Component)] +pub struct Payload; + + +/// A tether that connects the balloon to the payload. +#[derive(Component)] +pub struct Tether; diff --git a/src/simulator/properties.rs b/src/simulator/properties.rs new file mode 100644 index 0000000..6615665 --- /dev/null +++ b/src/simulator/properties.rs @@ -0,0 +1,376 @@ +//! Basic physical properties. + +#![allow(dead_code)] + +use std::ops::{Add, Div, Mul, Sub}; + +use avian3d::{ + math::{Scalar, PI}, + prelude::{ColliderDensity, ColliderMassProperties, PhysicsSet}, +}; +use bevy::{prelude::*, reflect::Reflect}; +#[cfg(feature = "config-files")] +use serde::{Deserialize, Serialize}; + +pub const BOLTZMANN_CONSTANT: f32 = 1.38e-23_f32; // [J/K] +pub const AVOGADRO_CONSTANT: f32 = 6.022e+23_f32; // [1/mol] + +pub const STANDARD_G: f32 = 9.80665; // [m/s^2] standard gravitational acceleration +pub const EARTH_RADIUS_M: f32 = 6371007.2; // [m] mean radius of Earth + +fn sphere_volume(radius: f32) -> f32 { + (4.0 / 3.0) * PI * f32::powf(radius, 3.0) +} + +fn shell_volume(internal_radius: f32, thickness: f32) -> f32 { + let external_radius = internal_radius + thickness; + let internal_volume = sphere_volume(internal_radius); + let external_volume = sphere_volume(external_radius); + external_volume - internal_volume +} + +fn sphere_radius_from_volume(volume: f32) -> f32 { + f32::powf(volume, 1.0 / 3.0) / (4.0 / 3.0) * PI +} + +fn sphere_surface_area(radius: f32) -> f32 { + 4.0 * PI * f32::powf(radius, 2.0) +} + +pub struct CorePropertiesPlugin; + +impl Plugin for CorePropertiesPlugin { + fn build(&self, app: &mut App) { + app.register_type::(); + app.register_type::(); + app.register_type::(); + app.register_type::(); + app.register_type::(); + app.register_type::(); + + // Ensure that the Avian density matches our computed mass and density + // before it starts solving physics. + app.add_systems( + Update, + (sync_avian_mass, sync_avian_density) + .chain() + .in_set(PhysicsSet::Prepare), + ); + } +} + +/// Temperature (K) +#[derive(Component, Debug, Default, Clone, Copy, PartialEq, Reflect)] +#[cfg_attr(feature = "config-files", derive(Serialize, Deserialize))] +pub struct Temperature(pub Scalar); + +impl Temperature { + pub const STANDARD: Self = Temperature(273.15); + + pub fn new(kelvin: f32) -> Self { + Temperature(kelvin) + } + + pub fn from_celsius(deg_celsius: f32) -> Self { + Temperature(deg_celsius + 273.15) + } + + pub fn kelvin(&self) -> f32 { + self.0 + } + + pub fn celsius(&self) -> f32 { + self.kelvin() - 273.15 + } +} + +impl Add for Temperature { + type Output = Temperature; + + fn add(self, rhs: Temperature) -> Self::Output { + Temperature(self.0 + rhs.0) + } +} + +impl Sub for Temperature { + type Output = Temperature; + + fn sub(self, rhs: Temperature) -> Self::Output { + Temperature(self.0 - rhs.0) + } +} + +impl Mul for Temperature { + type Output = Temperature; + + fn mul(self, rhs: Scalar) -> Self::Output { + Temperature(self.0 * rhs) + } +} + +impl Div for Temperature { + type Output = Temperature; + + fn div(self, rhs: Scalar) -> Self::Output { + Temperature(self.0 / rhs) + } +} + +/// Pressure (Pa) +#[derive(Component, Debug, Default, Clone, Copy, PartialEq, Reflect)] +#[cfg_attr(feature = "config-files", derive(Serialize, Deserialize))] +pub struct Pressure(pub Scalar); + +impl Pressure { + pub const STANDARD: Self = Pressure(101325.0); + + pub fn new(pascal: f32) -> Self { + Pressure(pascal) + } + + pub fn from_kilopascal(kilopascal: f32) -> Self { + Pressure(kilopascal * 1000.0) + } + + pub fn pascal(&self) -> f32 { + self.0 + } + + pub fn kilopascal(&self) -> f32 { + self.pascal() / 1000.0 + } +} + +impl Add for Pressure { + type Output = Pressure; + + fn add(self, rhs: Pressure) -> Self::Output { + Pressure(self.0 + rhs.0) + } +} + +impl Sub for Pressure { + type Output = Pressure; + + fn sub(self, rhs: Pressure) -> Self::Output { + Pressure(self.0 - rhs.0) + } +} + +impl Mul for Pressure { + type Output = Pressure; + + fn mul(self, rhs: Scalar) -> Self::Output { + Pressure(self.0 * rhs) + } +} + +impl Div for Pressure { + type Output = Pressure; + + fn div(self, rhs: Scalar) -> Self::Output { + Pressure(self.0 / rhs) + } +} + +/// The volume of a body in cubic meters. +#[derive(Component, Debug, Default, Clone, Copy, PartialEq, Reflect)] +#[cfg_attr(feature = "config-files", derive(Serialize, Deserialize))] +pub struct Volume(pub Scalar); + +impl Volume { + /// Zero volume. + pub const ZERO: Self = Self(0.0); + + pub fn cubic_meters(&self) -> f32 { + self.0 + } + + pub fn m3(&self) -> f32 { + self.0 + } +} + +impl Add for Volume { + type Output = Volume; + + fn add(self, rhs: Volume) -> Self::Output { + Volume(self.0 + rhs.0) + } +} + +impl Sub for Volume { + type Output = Volume; + + fn sub(self, rhs: Volume) -> Self::Output { + Volume(self.0 - rhs.0) + } +} + +impl Mul for Volume { + type Output = Volume; + + fn mul(self, rhs: Scalar) -> Self::Output { + Volume(self.0 * rhs) + } +} + +impl Div for Volume { + type Output = Volume; + + fn div(self, rhs: Scalar) -> Self::Output { + Volume(self.0 / rhs) + } +} + +/// Density (kg/m³) +#[derive(Component, Debug, Default, Clone, Copy, PartialEq, Reflect)] +#[cfg_attr(feature = "config-files", derive(Serialize, Deserialize))] +pub struct Density(pub Scalar); + +impl Density { + pub const ZERO: Self = Density(0.0); + + pub fn new(kilograms: Mass, volume: Volume) -> Self { + Density(kilograms.0 / volume.0) + } + + pub fn kilograms_per_cubic_meter(&self) -> f32 { + self.0 + } + + pub fn kg_per_m3(&self) -> f32 { + self.0 + } +} + +impl Add for Density { + type Output = Density; + + fn add(self, rhs: Density) -> Self::Output { + Density(self.0 + rhs.0) + } +} + +impl Sub for Density { + type Output = Density; + + fn sub(self, rhs: Density) -> Self::Output { + Density(self.0 - rhs.0) + } +} + +impl Mul for Density { + type Output = Density; + + fn mul(self, rhs: Scalar) -> Self::Output { + Density(self.0 * rhs) + } +} + +impl Div for Density { + type Output = Density; + + fn div(self, rhs: Scalar) -> Self::Output { + Density(self.0 / rhs) + } +} + +fn sync_avian_density(mut densities: Query<(&mut ColliderDensity, &Volume, &Mass)>) { + for (mut density, volume, mass) in densities.iter_mut() { + let our_density = mass.kg() / volume.m3(); + density.0 = our_density; + } +} + +/// Mass (kg) +#[derive(Component, Debug, Default, Clone, Copy, PartialEq, Reflect)] +#[cfg_attr(feature = "config-files", derive(Serialize, Deserialize))] +pub struct Mass(pub Scalar); + +impl Mass { + pub fn kilograms(&self) -> f32 { + self.0 + } + + pub fn kg(&self) -> f32 { + self.0 + } + + pub fn from_mass_properties(mass_props: ColliderMassProperties) -> Self { + Mass(mass_props.mass.0) + } +} + +impl Add for Mass { + type Output = Mass; + + fn add(self, rhs: Mass) -> Self::Output { + Mass(self.0 + rhs.0) + } +} + +impl Sub for Mass { + type Output = Mass; + + fn sub(self, rhs: Mass) -> Self::Output { + Mass(self.0 - rhs.0) + } +} + +impl Mul for Mass { + type Output = Mass; + + fn mul(self, rhs: Scalar) -> Self::Output { + Mass(self.0 * rhs) + } +} + +impl Div for Mass { + type Output = Density; + + fn div(self, rhs: Volume) -> Self::Output { + Density(self.0 / rhs.0) + } +} + +impl Div for Mass { + type Output = Mass; + + fn div(self, rhs: Scalar) -> Self::Output { + Mass(self.0 / rhs) + } +} + +fn sync_avian_mass(mut bodies: Query<(&mut ColliderMassProperties, &Mass)>) { + for (mut mass_props, mass) in bodies.iter_mut() { + mass_props.mass.0 = mass.0; + } +} + +/// Molar mass (kg/mol) of a substance. +#[derive(Component, Debug, Default, Clone, Copy, PartialEq, Reflect)] +#[cfg_attr(feature = "config-files", derive(Serialize, Deserialize))] +pub struct MolarMass(pub Scalar); + +impl MolarMass { + pub fn kilograms_per_mole(&self) -> f32 { + self.0 + } +} + +impl Mul for MolarMass { + type Output = MolarMass; + + fn mul(self, rhs: Scalar) -> Self::Output { + MolarMass(self.0 * rhs) + } +} + +impl Div for MolarMass { + type Output = MolarMass; + + fn div(self, rhs: Scalar) -> Self::Output { + MolarMass(self.0 / rhs) + } +} diff --git a/src/simulator/schedule.rs b/src/simulator/schedule.rs deleted file mode 100644 index f1e3803..0000000 --- a/src/simulator/schedule.rs +++ /dev/null @@ -1,283 +0,0 @@ -use log::{debug, info, error}; -use std::{ - path::PathBuf, - process::exit, - sync::{ - mpsc::Sender, - Arc, Mutex, - }, - thread::JoinHandle, - time::{Duration, Instant}, -}; - -use crate::simulator::{ - balloon::{Balloon, Material}, - bus::{Body, ParachuteSystem}, - config::{self, Config}, - forces, - gas::{Atmosphere, GasVolume}, - io::{SimCommands, SimOutput}, - SolidBody, -}; - -/// The simulation state at a given instant -pub struct SimInstant { - pub time: f32, - pub altitude: f32, - pub ascent_rate: f32, - pub acceleration: f32, - pub atmosphere: Atmosphere, - pub balloon: Balloon, - pub body: Body, - pub parachute: ParachuteSystem, -} - -fn initialize(config: &Config) -> SimInstant { - // create an initial time step based on the config - let atmo = Atmosphere::new(config.environment.initial_altitude_m); - let material = Material::new(config.balloon.material); - let mut lift_gas = GasVolume::new( - config.balloon.lift_gas.species, - config.balloon.lift_gas.mass_kg, - ); - lift_gas.update_from_ambient(atmo); - let body = Body::new(config.bus.body); - let parachute = - ParachuteSystem::new(config.bus.parachute, 1.0 / config.environment.tick_rate_hz); - - SimInstant { - time: 0.0, - altitude: config.environment.initial_altitude_m, - ascent_rate: config.environment.initial_velocity_m_s, - acceleration: 0.0, - atmosphere: atmo, - balloon: Balloon::new( - material, - config.balloon.thickness_m, - config.balloon.barely_inflated_diameter_m, // ballon diameter (m) - lift_gas, - ), - body, - parachute, - } -} - -pub fn step(input: SimInstant, config: &Config) -> SimInstant { - // propagate the closed loop simulation forward by one time step - let delta_t = 1.0 / config.environment.tick_rate_hz; - let time = input.time + delta_t; - let mut atmosphere = input.atmosphere; - let mut balloon = input.balloon; - let body = input.body; - let mut parachute = input.parachute; - - if balloon.intact { - // balloon is intact - balloon.stretch(atmosphere.pressure()); - } else { - parachute.deploy(atmosphere, input.ascent_rate); - }; - let total_dry_mass = body.total_mass() + parachute.total_mass(); - let weight_force = forces::weight(input.altitude, total_dry_mass); - let buoyancy_force = forces::buoyancy(input.altitude, atmosphere, balloon.lift_gas); - - let total_drag_force = forces::drag(atmosphere, input.ascent_rate, balloon) - + forces::drag(atmosphere, input.ascent_rate, body) - + forces::drag(atmosphere, input.ascent_rate, parachute.main) - + forces::drag(atmosphere, input.ascent_rate, parachute.drogue); - debug!( - "weight: {:?} buoyancy: {:?} drag: {:?}", - weight_force, buoyancy_force, total_drag_force - ); - - // calculate the net force - let net_force = weight_force + buoyancy_force + total_drag_force; - let acceleration = net_force / total_dry_mass; - let ascent_rate = input.ascent_rate + acceleration * delta_t; - let altitude = input.altitude + ascent_rate * delta_t; - - atmosphere.set_altitude(altitude); - - SimInstant { - time, - altitude, - ascent_rate, - acceleration, - atmosphere, - balloon, - body, - parachute, - } -} - -pub struct Rate { - cycle_time: Duration, - end_of_last_sleep: Option, -} - -impl Rate { - pub fn new(rate_hz: f32) -> Self { - Self { - cycle_time: Duration::from_secs_f32(1.0 / rate_hz), - end_of_last_sleep: None, - } - } - - pub fn sleep(&mut self) { - let now = Instant::now(); - - let sleep_duration = match self.end_of_last_sleep { - Some(v) => self - .cycle_time - .checked_sub(now.checked_duration_since(v).expect( - "Rate sleep experienced a last sleep with time ahead of the current time", - )) - .expect("Rate sleep detected a blown cycle"), - None => self.cycle_time, - }; - - std::thread::sleep(sleep_duration); - - self.end_of_last_sleep = Some(Instant::now()); - } -} - -pub struct AsyncSim { - config: Config, - sim_output: Arc>, - outpath: PathBuf, - command_sender: Option>, - /// keep track of - run_handle: Option>, -} - -impl AsyncSim { - pub fn new(config_path: &PathBuf, outpath: PathBuf) -> Self { - Self { - config: config::parse_from_file(config_path), - sim_output: Arc::new(Mutex::new(SimOutput::default())), - outpath, - command_sender: None, - run_handle: None, - } - } - - pub fn get_sim_output(&self) -> SimOutput { - *self.sim_output.lock().unwrap() - } - - /// Start a thread to run the sim - pub fn start(&mut self) { - if self.run_handle.is_some() { - panic!("Can't start again, sim already ran. Need to stop.") - } - let config = self.config.clone(); - let output = self.sim_output.clone(); - let outpath = self.outpath.clone(); - - debug!("Creating simulation handler..."); - self.run_handle = Some(std::thread::spawn(move || { - debug!("Simulation handler created. Initializing run..."); - AsyncSim::run_sim(config, output, outpath) - })); - } - - pub fn run_sim( - config: Config, - _sim_output: Arc>, - _outpath: PathBuf, - ) { - let mut sim_state = initialize(&config); - // configure simulation - let physics_rate = config.environment.tick_rate_hz; - let max_sim_time = config.environment.max_elapsed_time_s; - let real_time = config.environment.real_time; - let mut rate_sleeper = Rate::new(physics_rate); - - // set up data logger - // let mut writer = init_log_file(outpath); - - info!("Simulation run initialized. Starting loop..."); - loop { - if real_time { - rate_sleeper.sleep(); - } - sim_state = step(sim_state, &config); - - //log output - - // Print log to terminal - debug!( - "[{:.3} s] | Atmosphere @ {:} m: {:} K, {:} Pa", - sim_state.time, - sim_state.altitude, - sim_state.atmosphere.temperature(), - sim_state.atmosphere.temperature() - ); - debug!( - "[{:.3} s] | HAB @ {:.2} m, {:.3} m/s, {:.3} m/s^2 | {:.2} m radius, {:.2} Pa stress, {:.2} % strain", - sim_state.time, - sim_state.altitude, - sim_state.ascent_rate, - sim_state.acceleration, - sim_state.balloon.radius(), - sim_state.balloon.stress(), - sim_state.balloon.strain() * 100.0, - ); - // Stop if there is a problem - if sim_state.altitude.is_nan() - | sim_state.ascent_rate.is_nan() - | sim_state.acceleration.is_nan() - { - let status = format!("Something went wrong, a physical value is NaN!"); - #[cfg(feature = "gui")] - { - error!("{}", status); - break - } - #[cfg(not(feature = "gui"))] - { - Self::terminate(1, status); - } - } - // Run for a certain amount of sim time or to a certain altitude - if sim_state.time >= max_sim_time { - let status = format!("Reached maximum time step ({:?} s)", sim_state.time); - #[cfg(feature = "gui")] - { - info!("{}", status); - break - } - #[cfg(not(feature = "gui"))] - { - Self::terminate( - 0, status, - ); - } - } - if sim_state.altitude < 0.0 { - let status = format!("Altitude at or below zero."); - #[cfg(feature = "gui")] - { - info!("{}", status); - break - } - #[cfg(not(feature = "gui"))] - { - Self::terminate(0, status); - } - } - } - } - fn terminate(code: i32, reason: String) { - if code > 0 { - error!( - "Simulation terminated abnormally with code {:?}. Reason: {:?}", - code, reason - ); - } else { - info!("Simulation terminated normally. Reason: {:?}", reason); - } - exit(code); - } -}