From 8ba1b6584d407858396c5ef57a6192f6e22c4a4c Mon Sep 17 00:00:00 2001 From: Kyle Lastine Date: Mon, 6 Oct 2025 00:47:35 -0500 Subject: [PATCH 1/4] Add basic GNC state and control --- src/gnc.rs | 87 ++++++++++++++++++++++++++++++ src/lib.rs | 13 +++-- src/orbital.rs | 141 +++++++++++++++++++++++++++++++------------------ src/utils.rs | 82 ++++++++++++++++++++++++++++ 4 files changed, 270 insertions(+), 53 deletions(-) create mode 100644 src/gnc.rs diff --git a/src/gnc.rs b/src/gnc.rs new file mode 100644 index 0000000..46a9c63 --- /dev/null +++ b/src/gnc.rs @@ -0,0 +1,87 @@ +use crate::utils::{Vec3, EARTH_MU, EARTH_RADIUS}; + +#[derive(Clone, Copy, Debug)] +pub struct State { + pub r: Vec3, // position (km) + pub v: Vec3, // velocity (km/s) + pub m: f64, // mass (kg) - include if you want fuel usage; otherwise keep constant +} + +#[derive(Clone, Copy, Debug, Default)] +pub struct Control { + /// Thrust vector in ECI (N). If you prefer body frame + attitude, add attitude and rotate. + pub thrust_eci: Vec3, + /// Specific impulse (s) used for mass flow (optional). + pub isp: f64, +} + +#[inline] +fn two_body_accel(mu: f64, r: Vec3) -> Vec3 { + let rnorm = r.magnitude(); + let inv = 1.0 / (rnorm * rnorm * rnorm); + r * (-mu * inv) +} + +fn j2_accel(r: Vec3) -> Vec3 { + const J2: f64 = 1.08263e-3; + const RE: f64 = EARTH_RADIUS; // km + let z2 = r.z * r.z; + let r2 = r.dot(r); + let r1 = r2.sqrt(); + let rx = r.x; let ry = r.y; let rz = r.z; + let k = 1.5 * J2 * EARTH_MU * (RE * RE) / (r1.powi(5)); + let factor = 5.0 * z2 / r2; + Vec3::new( + rx * (factor - 1.0), + ry * (factor - 1.0), + rz * (factor - 3.0), + ) * k +} + +/// Mass flow from thrust (Tsiolkovsky; here as continuous mdot = -T/(Isp*g0)) +fn mdot_from_thrust(thrust: f64, isp: f64) -> f64 { + const G0: f64 = 9.80665; // m/s^2 + if isp > 0.0 { -(thrust) / (isp * G0) } else { 0.0 } // N / (s*m/s^2) -> kg/s +} + +/// Continuous-time EOM: xdot = f(x, u) +// time included but not used for future simulations +pub fn eom(_t: f64, x: &State, u: &Control) -> State { + // Forces + let a_grav = two_body_accel(EARTH_MU, x.r); + let a_j2 = j2_accel(x.r); + + // Thrust acceleration: a = F/m, convert N to (km/s^2): 1 N = 1 kg·m/s^2; 1 km = 1000 m + let thrust_km = u.thrust_eci / 1000.0; + let a_thrust = if x.m > 0.0 { thrust_km / x.m } else { Vec3::ZERO }; + + // Mass flow + let thrust_mag_n = u.thrust_eci.magnitude(); + let mdot = mdot_from_thrust(thrust_mag_n, u.isp); + + State { + r: x.v, // rdot = v + v: a_grav + a_j2 + a_thrust, // vdot = a_total + m: mdot, // mdot + } +} + +/// RK4 step +pub fn rk4_step(f: F, t: f64, x: &State, u: &Control, dt: f64) -> State +where + F: Fn(f64, &State, &Control) -> State, +{ + let k1 = f(t, x, u); + let x2 = State { r: x.r + k1.r*(dt*0.5), v: x.v + k1.v*(dt*0.5), m: x.m + k1.m*(dt*0.5) }; + let k2 = f(t + 0.5*dt, &x2, u); + let x3 = State { r: x.r + k2.r*(dt*0.5), v: x.v + k2.v*(dt*0.5), m: x.m + k2.m*(dt*0.5) }; + let k3 = f(t + 0.5*dt, &x3, u); + let x4 = State { r: x.r + k3.r*dt, v: x.v + k3.v*dt, m: x.m + k3.m*dt }; + let k4 = f(t + dt, &x4, u); + + State { + r: x.r + (k1.r + k2.r*2.0 + k3.r*2.0 + k4.r) * (dt / 6.0), + v: x.v + (k1.v + k2.v*2.0 + k3.v*2.0 + k4.v) * (dt / 6.0), + m: x.m + (k1.m + k2.m*2.0 + k3.m*2.0 + k4.m) * (dt / 6.0), + } +} diff --git a/src/lib.rs b/src/lib.rs index 6c2c987..63cfdfa 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,5 +1,6 @@ mod utils; mod orbital; +mod gnc; use wasm_bindgen::prelude::*; @@ -36,6 +37,9 @@ impl Satellite { raan: f64, arg_periapsis: f64, true_anomaly: f64, + mass: f64, + isp: f64, + thrust_limit_n: f64, ) -> Satellite { Satellite { inner: orbital::Satellite::new( @@ -45,6 +49,9 @@ impl Satellite { raan, arg_periapsis, true_anomaly, + mass, + isp, + thrust_limit_n, ), } } @@ -56,17 +63,17 @@ impl Satellite { /// Get current X position in kilometers pub fn get_x(&self) -> f64 { - self.inner.position.x + self.inner.state.r.x } /// Get current Y position in kilometers pub fn get_y(&self) -> f64 { - self.inner.position.y + self.inner.state.r.y } /// Get current Z position in kilometers pub fn get_z(&self) -> f64 { - self.inner.position.z + self.inner.state.r.z } /// Get orbital period in seconds diff --git a/src/orbital.rs b/src/orbital.rs index 7f81ae2..3135904 100644 --- a/src/orbital.rs +++ b/src/orbital.rs @@ -1,34 +1,7 @@ use crate::utils; +use crate::utils::{Vec3, EARTH_MU, EARTH_RADIUS}; +use crate::gnc; -// Physical constants -const EARTH_MU: f64 = 398600.4418; -const EARTH_RADIUS: f64 = 6371.0; - -#[derive(Debug, Clone, Copy)] -pub struct Vec3 { - pub x: f64, - pub y: f64, - pub z: f64, -} - -impl Vec3 { - pub fn new(x: f64, y: f64, z: f64) -> Self { - Vec3 { x, y, z } - } - - pub fn magnitude(&self) -> f64 { - libm::sqrt(self.x * self.x + self.y * self.y + self.z * self.z) - } - - pub fn normalize(&self) -> Vec3 { - let mag = self.magnitude(); - if mag > 0.0 { - Vec3::new(self.x / mag, self.y / mag, self.z / mag) - } else { - Vec3::new(0.0, 0.0, 0.0) - } - } -} /// Orbital elements (Keplerian elements) #[derive(Debug, Clone)] @@ -65,9 +38,11 @@ impl OrbitalElements { #[derive(Debug, Clone)] pub struct Satellite { pub elements: OrbitalElements, - pub position: Vec3, - pub velocity: Vec3, + pub state: gnc::State, + pub control: gnc::Control, pub time: f64, + pub thrust_limit_n: f64, + pub j2_on: bool, } impl Satellite { @@ -78,6 +53,9 @@ impl Satellite { raan: f64, arg_periapsis: f64, true_anomaly: f64, + mass: f64, + isp: f64, + thrust_limit_n: f64, ) -> Self { let elements = OrbitalElements::new( semi_major_axis, @@ -92,12 +70,15 @@ impl Satellite { Satellite { elements, - position, - velocity, + state: gnc::State { r: position, v: velocity, m: mass }, + control: gnc::Control { thrust_eci: Vec3::ZERO, isp: isp }, time: 0.0, + thrust_limit_n, + j2_on: true, } } + // https://en.wikipedia.org/wiki/Orbital_elements fn elements_to_state_vectors(elements: &OrbitalElements) -> (Vec3, Vec3) { let a = elements.semi_major_axis; let e = elements.eccentricity; @@ -140,23 +121,83 @@ impl Satellite { (Vec3::new(x, y, z), Vec3::new(vx, vy, vz)) } + + // https://en.wikipedia.org/wiki/Orbital_elements + pub fn state_vectors_to_elements(r: Vec3, v: Vec3) -> OrbitalElements { + // Constants & helpers + const EPS: f64 = 1e-10; - /// Update satellite position using simplified orbital propagation - /// This uses mean motion for circular/elliptical orbits - pub fn update(&mut self, delta_time: f64) { - self.time += delta_time; - - // Calculate mean motion (radians per second) - let n = libm::sqrt(EARTH_MU / libm::pow(self.elements.semi_major_axis, 3.0)); - - // Update true anomaly (simplified - assumes near-circular orbit) - self.elements.true_anomaly += n * delta_time; - self.elements.true_anomaly %= 2.0 * std::f64::consts::PI; - - // Recalculate position and velocity - let (position, velocity) = Self::elements_to_state_vectors(&self.elements); - self.position = position; - self.velocity = velocity; + let r_norm = r.magnitude(); + let v_norm2 = v.dot(v); + let h_vec = r.cross(v); + let h_norm = h_vec.magnitude(); + + // Inclination + let i = (h_vec.z / h_norm).acos(); + + // Node vector (pointing to ascending node): n = k × h + let n_vec = Vec3::new(-h_vec.y, h_vec.x, 0.0); + let n_norm = n_vec.magnitude(); + + let e_vec = (v.cross(h_vec) / EARTH_MU) - (r / r_norm); + let e = e_vec.magnitude(); + + let a = 1.0 / (2.0 / r_norm - v_norm2 / EARTH_MU); + + let raan = if n_norm > EPS { + utils::normalize_angle(libm::atan2(n_vec.y, n_vec.x)) + } else { + 0.0 + }; + + // Argument of periapsis (ω) + let arg_periapsis = if e > EPS && n_norm > EPS { + let cos_w = (n_vec.dot(e_vec)) / (n_norm * e); + let sin_w = (n_vec.cross(e_vec)).z / (n_norm * e); + utils::normalize_angle(libm::atan2(sin_w, cos_w)) + } else if e > EPS && n_norm <= EPS { + // Equatorial but eccentric: measure from x-axis + utils::normalize_angle(libm::atan2(e_vec.y, e_vec.x)) + } else { + 0.0 + }; + + // True anomaly (ν) + let true_anomaly = if e > EPS { + // cosν = (e·r)/(e|r|) + let cos_nu = (e_vec.dot(r)) / (e * r_norm); + let cos_nu = cos_nu.clamp(-1.0, 1.0); + let mut nu = libm::acos(cos_nu); + // Resolve quadrant using r·v + if r.dot(v) < 0.0 { + nu = 2.0 * std::f64::consts::PI - nu; + } + utils::normalize_angle(nu) + } else if n_norm > EPS { + let cos_u = (n_vec.dot(r)) / (n_norm * r_norm); + let sin_u = (n_vec.cross(r)).z / (n_norm * r_norm); + utils::normalize_angle(libm::atan2(sin_u, cos_u)) + } else { + utils::normalize_angle(libm::atan2(r.y, r.x)) + }; + + OrbitalElements::new(a, e, i, raan, arg_periapsis, true_anomaly) + } + + pub fn set_thrust_eci(&mut self, thrust_n_eci: Vec3) { + // clamp to actuator capability + let mag = thrust_n_eci.magnitude(); + self.control.thrust_eci = if mag > self.thrust_limit_n { + thrust_n_eci * (self.thrust_limit_n / mag) + } else { thrust_n_eci }; + } + + pub fn update(&mut self, dt: f64) { + // integrate one step + self.state = gnc::rk4_step(gnc::eom, self.time, &self.state, &self.control, dt); + self.time += dt; + + self.elements = Self::state_vectors_to_elements(self.state.r, self.state.v); } /// Get orbital period in seconds @@ -168,7 +209,7 @@ impl Satellite { /// Get current velocity magnitude in km/s pub fn get_velocity(&self) -> f64 { - self.velocity.magnitude() + self.state.v.magnitude() } /// Get apoapsis altitude (highest point) in kilometers diff --git a/src/utils.rs b/src/utils.rs index 0fbb98e..90921ef 100644 --- a/src/utils.rs +++ b/src/utils.rs @@ -1,3 +1,7 @@ +// Physical constants +pub const EARTH_MU: f64 = 398600.4418; +pub const EARTH_RADIUS: f64 = 6371.0; + /// Set panic hook for better error messages pub fn set_panic_hook() { #[cfg(feature = "console_error_panic_hook")] @@ -25,3 +29,81 @@ pub fn clamp(value: f64, min: f64, max: f64) -> f64 { } } +pub fn normalize_angle(mut ang: f64) -> f64 { + let two_pi = 2.0 * std::f64::consts::PI; + ang = ang % two_pi; + if ang < 0.0 { ang += two_pi; } + ang +} + +#[derive(Debug, Clone, Copy, Default)] +pub struct Vec3 { + pub x: f64, + pub y: f64, + pub z: f64, +} + +impl Vec3 { + pub const ZERO: Self = Self { x: 0.0, y: 0.0, z: 0.0 }; + + pub fn new(x: f64, y: f64, z: f64) -> Self { + Self { x, y, z } + } + + pub fn dot(&self, other: Self) -> f64 { + self.x * other.x + self.y * other.y + self.z * other.z + } + + pub fn magnitude(&self) -> f64 { + (self.x * self.x + self.y * self.y + self.z * self.z).sqrt() + } + + pub fn normalize(&self) -> Self { + let mag = self.magnitude(); + if mag == 0.0 { + *self + } else { + *self / mag + } + } + + pub fn cross(&self, other: Self) -> Self { + Self { + x: self.y * other.z - self.z * other.y, + y: self.z * other.x - self.x * other.z, + z: self.x * other.y - self.y * other.x, + } + } +} + +impl std::ops::Add for Vec3 { + type Output = Vec3; + + fn add(self, other: Vec3) -> Vec3 { + Vec3::new(self.x + other.x, self.y + other.y, self.z + other.z) + } +} + +impl std::ops::Sub for Vec3 { + type Output = Vec3; + + fn sub(self, other: Vec3) -> Vec3 { + Vec3::new(self.x - other.x, self.y - other.y, self.z - other.z) + } +} + +impl std::ops::Mul for Vec3 { + type Output = Vec3; + + fn mul(self, scalar: f64) -> Vec3 { + Vec3::new(self.x * scalar, self.y * scalar, self.z * scalar) + } +} + +impl std::ops::Div for Vec3 { + type Output = Vec3; + + fn div(self, scalar: f64) -> Vec3 { + Vec3::new(self.x / scalar, self.y / scalar, self.z / scalar) + } +} \ No newline at end of file From 8c047445e0bbb53df07e559b580202d6fa3c48f3 Mon Sep 17 00:00:00 2001 From: Kyle Lastine Date: Mon, 6 Oct 2025 01:33:03 -0500 Subject: [PATCH 2/4] ci rework --- .github/workflows/deploy.yml | 64 +++++++++++++++++++++++++++++------- 1 file changed, 53 insertions(+), 11 deletions(-) diff --git a/.github/workflows/deploy.yml b/.github/workflows/deploy.yml index 7ede708..d54b668 100644 --- a/.github/workflows/deploy.yml +++ b/.github/workflows/deploy.yml @@ -12,11 +12,63 @@ permissions: id-token: write jobs: + build: + runs-on: ubuntu-latest + if: github.event_name == 'pull_request' + + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Setup Rust + uses: actions-rs/toolchain@v1 + with: + toolchain: stable + components: rustfmt, clippy + + - name: Cache Rust dependencies + uses: actions/cache@v3 + with: + path: | + ~/.cargo/registry + ~/.cargo/git + target + key: ${{ runner.os }}-cargo-${{ hashFiles('**/Cargo.lock') }} + + - name: Install wasm-pack + run: curl https://rustwasm.github.io/wasm-pack/installer/init.sh -sSf | sh + + - name: Build WASM + run: wasm-pack build --target web --out-dir pkg + + - name: Setup Node.js + uses: actions/setup-node@v4 + with: + node-version: '18' + cache: 'npm' + cache-dependency-path: www/package-lock.json + + - name: Install frontend dependencies + run: | + cd www + npm ci + + - name: Copy WASM to frontend + run: | + cp -r pkg www/ + + - name: Build frontend + run: | + cd www + npm run build + ls -la dist/ + build-and-deploy: runs-on: ubuntu-latest environment: name: github-pages url: ${{ steps.deployment.outputs.page_url }} + if: github.ref == 'refs/heads/main' steps: - name: Checkout @@ -66,24 +118,14 @@ jobs: ls -la dist/ - name: Setup Pages - if: github.ref == 'refs/heads/main' uses: actions/configure-pages@v4 - name: Upload artifact - if: github.ref == 'refs/heads/main' uses: actions/upload-pages-artifact@v3 with: path: www/dist - name: Deploy to GitHub Pages - if: github.ref == 'refs/heads/main' id: deployment uses: actions/deploy-pages@v4 - - - name: Deploy PR preview - if: github.event_name == 'pull_request' - uses: peaceiris/actions-gh-pages@v3 - with: - github_token: ${{ secrets.GITHUB_TOKEN }} - publish_dir: www/dist - destination_dir: preview/${{ github.event.number }} + From 85c625716d74f5d6419fead4a9385cee2a06b72f Mon Sep 17 00:00:00 2001 From: Kyle Lastine Date: Mon, 6 Oct 2025 01:35:21 -0500 Subject: [PATCH 3/4] fix ci --- www/src/components/OrbitCanvas.tsx | 10 ++++++++-- www/src/wasm.ts | 9 ++++++--- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/www/src/components/OrbitCanvas.tsx b/www/src/components/OrbitCanvas.tsx index 841f266..ed05bd9 100644 --- a/www/src/components/OrbitCanvas.tsx +++ b/www/src/components/OrbitCanvas.tsx @@ -205,7 +205,10 @@ const OrbitCanvas: React.FC = () => { 51.6, // inclination (degrees) 0.0, // RAAN (degrees) 0.0, // argument of periapsis (degrees) - 0.0 // initial true anomaly (degrees) + 0.0, // initial true anomaly (degrees) + 1000.0, // mass (kg) + 300.0, // isp (specific impulse in seconds) + 1000.0 // thrust limit (N) ) console.log('Satellite created successfully:', satelliteRef.current) @@ -277,7 +280,10 @@ const OrbitCanvas: React.FC = () => { inclinationDeg, 0, 0, - 0 + 0, + 1000.0, // mass (kg) + 300.0, // isp (specific impulse in seconds) + 1000.0 // thrust limit (N) ) } // reset integrator clock to avoid jumps diff --git a/www/src/wasm.ts b/www/src/wasm.ts index b954a2d..12a9f4a 100644 --- a/www/src/wasm.ts +++ b/www/src/wasm.ts @@ -17,9 +17,12 @@ export interface WasmModule { inclination: number, raan: number, argPeriapsis: number, - trueAnomaly: number + trueAnomaly: number, + mass: number, + isp: number, + thrustLimitN: number ) => Satellite; - default: () => Promise; // WASM initialization function + init: () => void; // WASM initialization function } // Load WASM module @@ -27,7 +30,7 @@ export const loadWasm = async (): Promise => { try { // Import the WASM module from the local pkg directory const wasm = await import('../pkg/orbit_crab.js') - return wasm as WasmModule + return wasm as unknown as WasmModule } catch (error) { console.error('Failed to load WASM module:', error) throw error From 0fcd3d99ee46987929126ca74b4c25e36824f666 Mon Sep 17 00:00:00 2001 From: Kyle Lastine Date: Mon, 6 Oct 2025 01:36:19 -0500 Subject: [PATCH 4/4] fix ci --- www/src/components/OrbitCanvas.tsx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/www/src/components/OrbitCanvas.tsx b/www/src/components/OrbitCanvas.tsx index ed05bd9..58fe232 100644 --- a/www/src/components/OrbitCanvas.tsx +++ b/www/src/components/OrbitCanvas.tsx @@ -191,7 +191,7 @@ const OrbitCanvas: React.FC = () => { console.log('WASM loaded successfully:', wasm) console.log('Initializing WASM...') - await wasm.default() + wasm.init() console.log('WASM initialized successfully') // Create ISS-like satellite (408 km altitude, circular orbit, 51.6° inclination)