# zmatrix — Physics-Aware Linear Algebra in Rust
[](https://crates.io/crates/zmatrix)
[](https://codecov.io/gh/Treagzhao/zmatrix)
[](LICENSE)
A Rust math library for aerospace, robotics, and physics simulation. Features a compile-time dimensional analysis system and const-generic dense matrices.
> 📖 [中文 README](README.MD)
---
## Quick Look
```rust
use zmatrix::physics::basic::*;
use zmatrix::dense::Matrix;
use chrono::TimeDelta;
// Automatic unit conversion
let d = Distance::from_m(1000.0);
println!("{} km", d.as_km()); // 1 km
// Compile-time dimensional analysis
let v: Velocity = d / time_delta_from_secs_f64(10.0); // ✅ Correct
let a: Acceleration = v / time_delta_from_secs_f64(5.0); // ✅ Correct
// let wrong: Area = d / time_delta_from_secs_f64(10.0); // ❌ Compiler error: length / time is not area
```
> **Physical quantities cannot be matrix elements.**
> `T` in `Matrix<ROWS, COLS, T>` is bounded by `NumericElement` and thus restricted to
> primitive numbers (integers and floats). A single `T` forces the whole matrix to share one
> dimension, which cannot express the heterogeneous matrices engineers actually need
> (state-transition, covariance, Kalman gain, Jacobian). For dimensionless rotation matrices
> use
> [`CosMatrix`](https://docs.rs/zmatrix/latest/zmatrix/spatial_geometry/cos_matrix/struct.CosMatrix.html).
---
## 1. Physical Quantity System
32 physical quantities, each with multiple units. Operations between quantities are dimensionally checked at compile time.
### Supported Quantities
| `Distance` | m, km, ly |
| `Velocity` | m/s, km/h, km/s, c |
| `Acceleration` | m/s², km/h², g |
| `Angular` | rad, deg |
| `AngularVelocity` | rad/s, deg/s, rad/h, deg/h |
| `AngularAcceleration` | rad/s², deg/s² |
| `Area` | m², km² |
| `Volume` | m³, km³ |
| `Mass` | kg, g |
| `Inertia` | kg·m², kg·mm² |
| `Momentum` | kg·m/s, kg·km/s |
| `AngularMomentum` | kg·m²/s, kg·km²/s, N·m·s family |
| `Force` | N, mN, μN, nN, kN, MN |
| `Torque` | N·m, mN·m, μN·m, nN·m, kN·m, MN·m |
| `Energy` | J, eV, mJ, μJ, nJ, kJ, MJ |
| `Power` | W, hp, mW, μW, nW, kW, MW |
| `MagneticInduction` | T, G, mT, μT, nT |
| `MagneticMoment` | A·m², J/T |
| `MagneticAngularVelocity` | T·rad/s |
| `ElectricCurrent` | A, mA, μA, kA |
| `ElectricPotential` | V, mV, μV, kV |
| `ElectricResistance` | Ω, mΩ, kΩ, MΩ |
| `ElectricCharge` | C, mC, μC, nC |
| `ElectricCapacitance` | F, mF, μF, nF, pF |
| `ElectricConductance` | S, mS, μS |
| **`Temperature`** | K, °C, °F, °R |
| **`TemperatureDifference`** | ΔK, Δ°C, Δ°F, Δ°R |
| **`HeatFlux`** | W/m², kW/m², mW/cm² |
| **`ThermalConductivity`** | W/(m·K), kW/(m·K), mW/(cm·K) |
| **`HeatCapacity`** | J/K, kJ/K |
| **`ThermalConductance`** | W/K, mW/K |
| `Coef` | dimensionless |
### Time: `chrono::TimeDelta` only
Every time-taking operation accepts `chrono::TimeDelta`. `std::time::Duration` is **not supported**.
Rationale: `Duration` cannot represent negative spans, so retrograde burns and reversed
time stepping are inexpressible. The dual-track support was removed in `0.4.1`.
```rust
use chrono::TimeDelta;
let dt = TimeDelta::seconds(10);
let dt2 = TimeDelta::milliseconds(500);
let dt3 = time_delta_from_secs_f64(2.5); // float seconds
let dt4 = time_delta_from_secs_f64(-2.5); // negative — inexpressible with Duration
```
### Quantity Arithmetic
```rust
// Kinematics
let distance = Distance::from_m(1000.0);
let dt = time_delta_from_secs_f64(10.0);
let velocity: Velocity = distance / dt; // 100 m/s
let acceleration: Acceleration = velocity / dt; // 10 m/s²
// Dynamics
let mass = Mass::from_kg(2.0);
let force: Force = mass * acceleration; // F = m × a → 20 N
let energy: Energy = force * distance; // W = F × d → 20000 J
// Impulse–momentum theorem: F·Δt = Δp
let impulse: Momentum = force * dt; // 20 N × 10 s → 200 kg·m/s
// Rotational kinematics
let angle = Angular::from_deg(180.0);
let omega: AngularVelocity = angle / dt;
let alpha: AngularAcceleration = omega / dt;
// Angular impulse: M·Δt = ΔL
let torque = Torque::from_nm(5.0);
let angular_impulse: AngularMomentum = torque * dt;
// Electricity (Ohm's law & power)
let current = ElectricCurrent::from_a(2.0);
let resistance = ElectricResistance::from_ohm(5.0);
let voltage: ElectricPotential = current * resistance; // V = I × R → 10 V
let power: Power = voltage * current; // P = V × I → 20 W
// Electrical energy
let charge = ElectricCharge::from_coulomb(3.0);
let energy2: Energy = charge * voltage; // E = Q × V → 30 J
```
---
## 2. Matrix Operations
Generic matrices backed by const generics — rows and columns fixed at compile time.
### Basics
```rust
use zmatrix::dense::Matrix;
let mut m1 = Matrix::<2, 3, f64>::new([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]);
let m2 = Matrix::<2, 3, f64>::new([[1.0, 1.0, 1.0], [1.0, 1.0, 1.0]]);
let sum = m1 + m2;
let diff = m1 - m2;
let scaled = m1.scale(2.0);
let transposed = m1.T(); // transpose
// All-zero / all-one matrices
let zeros = Matrix::<3, 3, f64>::zeros();
let ones = Matrix::<3, 3, f64>::ones();
// Get & set: argument order is (column, row)
let val = m1.get(2, 1).unwrap(); // column 2, row 1
m1.set(2, 1, 99.0).unwrap();
```
### Matrix Product
```rust
let a = Matrix::<2, 3, f64>::new([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]);
let b = Matrix::<3, 2, f64>::new([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]]);
let c: Matrix<2, 2, f64> = a.product(&b).unwrap();
```
Both operands must share the same element type, which is restricted to primitive numbers
(see the note on physical quantities above).
> Note: `*` is **element-wise** (Hadamard) multiplication. Use `.product()` for matrix
> multiplication.
---
## 3. Vector Operations
3D vectors over multiple physical quantity types, with cross/dot products.
```rust
let displacement = Vector3::new(
Distance::from_m(10.0), Distance::from_m(20.0), Distance::from_m(30.0),
);
// Displacement ÷ time = velocity vector
let velocity = displacement / time_delta_from_secs_f64(10.0);
// Velocity × mass = momentum vector
let momentum = velocity * Mass::from_kg(1.0);
// Displacement × momentum = angular momentum
let angular_momentum = displacement * momentum;
// Dot product: takes a reference, returns the same type
let v1 = Vector3::new(Distance::from_m(1.0), Distance::from_m(0.0), Distance::from_m(0.0));
let v2 = Vector3::new(Distance::from_m(0.0), Distance::from_m(1.0), Distance::from_m(0.0));
let dot: Distance = v1.dot(&v2);
// Cross product: takes a reference, output type must implement VectorQuantity
let c1 = Vector3::new(Coef::new(1.0), Coef::new(0.0), Coef::new(0.0));
let c2 = Vector3::new(Coef::new(0.0), Coef::new(1.0), Coef::new(0.0));
let cross = c1.cross(&c2);
// Skew-symmetric matrices
let skew = v1.skew_symmetric_matrix(); // 3×3
let skew4 = v1.skew_symmetric_matrix_4(); // 4×4
```
> **On the cross-product type bound**: the output type of `cross` is the product of the two
> element types, and that type must implement `VectorQuantity`. Since
> `Distance × Distance = Area` and `Area` is not a vector quantity, `Vector3<Distance>` **cannot**
> be cross-multiplied; dimensionless `Vector3<Coef>` can.
### Vector time arithmetic
The six operations needed for force-domain stepping, all taking `chrono::TimeDelta`:
```rust
let dt = TimeDelta::seconds(2);
// Velocity ÷ time = acceleration vector
let velocity_vec = Vector3::new(
Velocity::from_m_per_sec(10.0),
Velocity::from_m_per_sec(4.0),
Velocity::from_m_per_sec(0.0),
);
let acc: Vector3<Acceleration> = velocity_vec / dt;
// Acceleration × time = velocity increment
let dv: Vector3<Velocity> = acc * dt;
// Force × time = impulse (vector form of the impulse–momentum theorem)
let force_vec = Vector3::new(
Force::from_newton(1.0),
Force::from_newton(2.0),
Force::from_newton(3.0),
);
let impulse_vec: Vector3<Momentum> = force_vec * dt;
// Torque × time = angular impulse
let torque_vec = Vector3::new(
Torque::from_nm(5.0),
Torque::from_nm(0.0),
Torque::from_nm(-1.0),
);
let ang_impulse: Vector3<AngularMomentum> = torque_vec * dt;
```
---
## 4. Rotational Inertia and Attitude Dynamics
Scalar rotational inertia `Inertia` (kg·m²), together with angular velocity, angular
acceleration, angular momentum and torque, closes the whole rigid-body rotational chain
inside the type system.
### Six mutually inverse operations
```rust
let j = Inertia::from_kg_m2(0.5); // flywheel rotor inertia
let omega = AngularVelocity::from_rad_per_second(10.0);
let torque = Torque::from_nm(3.0);
// Angular momentum: J·ω = L
let l: AngularMomentum = j * omega; // 5 kg·m²/s
// Solve back for spin rate: L / J = ω
let omega2: AngularVelocity = l / j;
// Angular acceleration: M / J = ω̇
let alpha: AngularAcceleration = torque / j;
// Torque: J·ω̇ = M
let torque2: Torque = j * alpha;
// Solve back for inertia: L / ω = J and M / ω̇ = J, two independent routes
let j1: Inertia = l / omega;
let j2: Inertia = torque / alpha;
```
### Two units
Engineering drawings and CAD usually report inertia in kg·mm² (10⁶ times larger than kg·m²):
```rust
let j_from_drawing = Inertia::from_kg_mm2(500_000.0);
assert_eq!(j_from_drawing.as_kg_m2(), 0.5); // converted automatically
```
### Full flywheel spin-up
```rust
// Flywheel with J = 0.5 kg·m², 3 N·m applied for 2 seconds
let j = Inertia::from_kg_m2(0.5);
let torque = Torque::from_nm(3.0);
let dt = TimeDelta::seconds(2);
let impulse: AngularMomentum = torque * dt; // angular impulse 6 N·m·s
let omega: AngularVelocity = impulse / j; // ω = 12 rad/s
let alpha: AngularAcceleration = torque / j; // ω̇ = 6 rad/s²
assert_eq!(alpha.as_rad_per_second2() * 2.0, omega.as_rad_per_second()); // both routes agree
```
> **Scalar vs. tensor**: `Inertia` describes rotational inertia about **one specified axis**.
> Full description of a rigid body about any axis needs the 3×3 tensor (with products of
> inertia) — see the next chapter.
---
## 5. Inertia Tensor and the `compound` Module
`compound` holds **structures assembled from physical quantities**, as opposed to
`physics::basic` (single quantities):
| `physics::basic` | single quantities (one value + unit) | — |
| `spatial_geometry` | pure geometry, dimensionless | all `f64` |
| **`compound`** | structures assembled from quantities | dimensioned |
Dependencies are one-way: `compound → physics::basic`, `compound → spatial_geometry`.
The reverse is not allowed.
### `InertiaTensor`
```rust
use zmatrix::compound::inertia_tensor::InertiaTensor;
// Six independent components. Symmetry is guaranteed by the interface —
// an asymmetric tensor cannot be constructed.
let j = InertiaTensor::new(
Inertia::from_kg_m2(120.0), // Ixx
Inertia::from_kg_m2(130.0), // Iyy
Inertia::from_kg_m2(45.0), // Izz
Inertia::from_kg_m2(0.0), // Ixy
Inertia::from_kg_m2(8.0), // Ixz
Inertia::from_kg_m2(0.0), // Iyz
);
// Principal axes / isotropic
let j2 = InertiaTensor::from_diagonal([
Inertia::from_kg_m2(1.0),
Inertia::from_kg_m2(2.0),
Inertia::from_kg_m2(3.0),
]);
let j3 = InertiaTensor::from_isotropic(Inertia::from_kg_m2(2.0));
// Angular momentum: J·ω = L. For an asymmetric tensor, L deviates from ω.
let omega = Vector3::new(
AngularVelocity::from_rad_per_second(0.0),
AngularVelocity::from_rad_per_second(0.0),
AngularVelocity::from_rad_per_second(1.0),
);
let l: Vector3<AngularMomentum> = j * omega;
// Torque: J·ω̇ = M
let alpha = Vector3::new(
AngularAcceleration::from_rad_per_second2(1.0),
AngularAcceleration::from_rad_per_second2(0.0),
AngularAcceleration::from_rad_per_second2(0.0),
);
let m: Vector3<Torque> = j * alpha;
// Inertia about an arbitrary axis: nᵀJn / nᵀn
// (normalized internally, n need not be a unit vector)
let dir = Vector3::new(Coef::new(1.0), Coef::new(1.0), Coef::new(0.0));
let about: Inertia = j.about_axis(dir)?;
// Isotropic equivalent: trace/3
let iso = j.isotropic_equivalent();
// Coordinate transform: R·J·Rᵀ
let rotated = j.rotate(&CosMatrix::unit());
```
> **`about_axis` cannot be replaced by a diagonal element.** For
> `[[2,-1,0],[-1,2,0],[0,0,3]]`, the z-axis value is 3 but the `(1,1,0)/√2` axis gives 1 —
> products of inertia pull it down from 2. Normalization happens inside the method, so a
> non-unit vector gives the same answer; a zero vector returns an error.
### Gyroscopic torque
The point of the tensor: `L` can deviate from `ω`, so `ω × L` is nonzero.
```text
ω = (0, 0, 1), L = (8, 0, 45) because Ixz ≠ 0
ω × L = (0, 8, 0) nonzero — why attitude dynamics needs the tensor
```
---
## 6. Temperature and Thermal Quantities
### Temperature and temperature difference must be separate types
Celsius and kelvin are related by an **affine** transform with an offset
(`T_K = T_C + 273.15`), while differences scale linearly with no offset:
```text
20 °C = 293.15 K temperature has an offset
Δ20 °C = 20 K difference does not
```
Collapsing them into one type has two consequences:
- `(t1 - t2).as_celsius()` returns nonsense like `-253.15`
- `G * T` and `G * ΔT` both typecheck (`W/K × K = W`) but differ by tens of times, with
**dimensional analysis unable to tell them apart**
So `Temperature` and `TemperatureDifference` are separate, and these operations are
deliberately absent:
| `T - T → ΔT` | ✅ | the only meaningful subtraction |
| `T + ΔT → T` | ✅ | heating / cooling (commutative) |
| `T / T → Coef` | ✅ | temperature ratio, e.g. Carnot efficiency |
| `T + T` | ❌ | sum of two temperatures is not a temperature |
| `T × T` | ❌ | output would be K², which has no use here |
| `T × f64` | ❌ | misleading from Celsius: 20 °C doubled gives 313 °C |
```rust
let t1 = Temperature::from_celsius(30.0);
let t2 = Temperature::from_celsius(20.0);
let dt: TemperatureDifference = t1 - t2; // 10 K
let heated = t2 + dt; // 30 °C
let ratio: Coef = t1 / t2; // based on kelvin, not 30/20
```
> Temperature is a scalar, so `Vector3<Temperature>` is unconstructable at compile time.
### Heat flux and Fourier's law
**Heat flux density** `HeatFlux` (W/m²) is thermal power per unit area. Distinct from two
easily confused quantities:
| Heat Q | J | accumulated energy | `Energy` |
| Thermal power Q̇ | W | energy per time | `Power` |
| Heat flux q | W/m² | energy per time per area | `HeatFlux` |
```rust
let q = HeatFlux::from_w_per_m2(1361.0); // solar constant
let p: Power = q * Area::from_m2(2.0); // 2722 W
// Datasheets often use mW/cm²: 1 mW/cm² = 10 W/m² (not 1e-3)
let q2 = HeatFlux::from_mw_per_cm2(1.0);
assert_eq!(q2.as_w_per_m2(), 10.0);
```
**Thermal conductivity** `ThermalConductivity` (W/(m·K)) is a material property.
Fourier's law:
```rust
let k = ThermalConductivity::from_w_per_m_k(237.0); // aluminium
let dt = TemperatureDifference::from_celsius(20.0);
let thickness = Distance::from_m(0.01);
let q: HeatFlux = k.heat_flux(dt, thickness); // 474 kW/m²
// Inverse — calibrate k from a measured heat flux
let k2 = ThermalConductivity::from_flux(q, dt, thickness);
```
> Takes `TemperatureDifference`, not `Temperature`. Passing an absolute temperature is a
> compile error — otherwise you get a heat flux measured from 0 K.
### The two elements of a lumped-parameter thermal network
$$C\frac{dT}{dt} = G\,\Delta T + \dot{Q}$$
| C | `HeatCapacity` | J/K | node storage |
| G | `ThermalConductance` | W/K | coupling between nodes |
| ΔT | `TemperatureDifference` | K | temperature difference |
| Q̇ | `Power` | W | external source |
```rust
// Storage: Q = C·ΔT
let c = HeatCapacity::from_j_per_k(9000.0);
let q: Energy = c * TemperatureDifference::from_celsius(20.0); // 180 kJ
// Inverse: ΔT = Q/C
let dt: TemperatureDifference = q / c;
// Transmission: Q̇ = G·ΔT
let g = ThermalConductance::from_w_per_k(12.0);
let p: Power = g * TemperatureDifference::from_celsius(20.0); // 240 W
```
> **Four similar names, two layers** (all distinct):
>
> | Layer | Storage | Transmission |
> |---|---|---|
> | Material property | specific heat capacity c, J/(kg·K) | thermal conductivity k, W/(m·K) |
> | Whole object | **heat capacity C, J/K** | **thermal conductance G, W/K** |
>
> Relations: `C = c·m`, `G = k·A/L`.
---
## 7. Spatial Geometry
Euler angles, direction cosine matrices (DCM), and quaternions with mutual conversions.
### Euler Angles
```rust
let euler: Vector3<Angular> = Vector3::new(
Angular::from_deg(10.0), Angular::from_deg(20.0), Angular::from_deg(30.0),
);
let quat = euler.to_quaternion();
let sin = euler.sin();
let cos = euler.cos();
```
### Direction Cosine Matrix
```rust
let cos = CosMatrix::unit();
let x_row = cos.get_x_row_vector();
let y_row = cos.get_y_row_vector();
let z_row = cos.get_z_row_vector();
let x_col = cos.get_x_col_vector();
let y_col = cos.get_y_col_vector();
let z_col = cos.get_z_col_vector();
let cos_t = cos.transfer(); // transpose
let q = cos.to_quaternion();
let euler = cos.to_pry(); // XZY sequence
let euler = cos.to_rpy(); // XYZ sequence
// Matrix product: takes references; the vector element type must implement
// VectorQuantity (e.g. Coef)
let cos_a = CosMatrix::unit();
let cos_b = CosMatrix::unit();
let cos_c = cos_a.product(&cos_b);
let dir = Vector3::new(Coef::new(1.0), Coef::new(0.0), Coef::new(0.0));
let vec = cos_c.product_vector(&dir);
```
### Quaternions
```rust
let q = Quaternion::new(1.0, 2.0, 3.0, 4.0);
let norm = q.norm();
let normalized = q.normalize();
let conjugate = q.conjugate();
let inverse = q.inverse();
let cos = q.to_cos_matrix();
let q1 = Quaternion::new(1.0, 0.0, 0.0, 0.0);
let q2 = Quaternion::new(0.0, 1.0, 0.0, 0.0);
let sum = q1 + q2;
let product = q1 * q2;
let quotient = q1 / q2;
```
---
## Not yet implemented
Ordered by current priority. **Deliberately out of scope** items are listed last so they do
not get re-litigated.
### Linear algebra gaps (blocking)
`dense::Matrix` currently has only add/subtract, scalar multiply, transpose, and matrix
product. Solving is missing entirely:
| `det()` | singularity checks, eigenvalues |
| `inverse()` | explicit inversion |
| `solve()` / LU | **steady-state thermal network `G·T = Q̇`**, any linear system |
This is a prerequisite for steady-state thermal analysis and for inverting `InertiaTensor`
(to solve `J·ω = L`). Note the dimensional problem: `J⁻¹` needs `(kg·m²)⁻¹` and `det` needs
`(kg·m²)³`, neither of which exists as a quantity. So even with a solver, it can only be
offered at the `f64` level, not as a dimensioned interface.
### Thermal quantities
| `SpecificHeatCapacity` | J/(kg·K) | material property, `C = c·m` |
| `HeatTransferCoefficient` | W/(m²·K) | surface convection, `q = h·ΔT`. Rare in vacuum, low priority |
| `ThermalResistance` | K/W | `R = 1/G`, clearer for series paths |
| Radiation helpers | — | `σT⁴` has no place in the type system; functions only: `radiative_flux(ε, T)` |
**Deliberately not planned**:
- `Entropy` (J/K) — unused in satellite thermal work; wait for a concrete need
- `ThermalDiffusivity` (m²/s) — derivable as `α = k/(ρc)`, redundant
- `Heat`, `InternalEnergy` — duplicate `Energy`
- `SpecificHeatRatio` γ — dimensionless quotient, `Coef` suffices; for solids `cp ≈ cv` and
`γ ≈ 1`, so it carries no information
- `Emissivity` — use `Coef` for now; no value-range constraints exist anywhere in this crate
### `compound` module
| `InertiaTensor` parallel-axis theorem | `J' = J + m(d²I - ddᵀ)`, needs `Mass` + `Distance` |
| `InertiaTensor` principal-axis diagonalization | Jacobi iteration to zero the products |
| `InertiaTensor` inverse / determinant | blocked by the linear algebra gap above |
| `ThermalNetwork` container | see below |
**Planned `ThermalNetwork` scope** (deliberately narrow):
```text
add_node(capacity, temperature) add a node
add_link(i, j, conductance) add a link
set_source(i, power) attach a heat source
step(dt) -> Result<(), _> advance one explicit step
temperature(i) -> Option<Temperature>
```
**Explicit stepping only**: implicit stepping and steady-state solving both need a linear
solver, blocked as above. `step()` will check the explicit stability criterion
`Δt < 2·Cᵢ/ΣGᵢⱼ` and return an error — without that check, a diverging solution still
produces plausible-looking curves early on, and by the time it blows up exponentially you
are already reading wrong conclusions.
### Other
- `Matrix` `det` / `inverse` / `solve` (see above)
- Quantity-to-quantity operators are still incomplete; additions land when a concrete case
needs them
- The unit system covers a practical engineering subset, not the full SI prefix hierarchy
---
## Dependencies
```toml
[dependencies]
zmatrix = "0.4"
chrono = "0.4"
```
> Time arithmetic uses `chrono::TimeDelta`, so `chrono` must be in scope.
> The helpers `time_delta_from_secs_f64` / `time_delta_to_secs_f64` let you avoid
> importing `chrono` directly.
---
## License
MIT