diff --git a/CHANGELOG.md b/CHANGELOG.md
index b922c92..5c1cbb6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,26 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
## [Unreleased]
+## [0.9.3] - 2026-02-26
+
+### Added
+
+- Natural units support with `ConstantAwareBasisTransform` (#206)
+ - `NATURAL` basis where c = ℏ = k_B = 1 (all quantities reduce to powers of energy)
+ - `SI_TO_NATURAL` transform maps SI dimensions to natural units (e.g., velocity → dimensionless)
+ - `NATURAL_TO_SI` inverse transform with constant bindings for reconstruction
+ - `ConstantBinding` records relationships between dimensions via physical constants
+ - `LossyProjection` exception for dimensions without natural unit representation (e.g., current)
+ - `allow_projection=True` option to drop unrepresentable dimensions
+- Example demos for alternative unit systems (`examples/basis/`)
+ - `natural_units_demo.py` — particle physics natural units
+ - `geometrized_units_demo.py` — general relativity units (c = G = 1)
+ - `elemental_units_demo.py` — custom "elemental" basis
+- Natural units guide in `docs/guides/natural-units.md`
+- Natural units section in API reference `docs/reference/api.md`
+- Custom constants documentation in README and API reference
+- Comprehensive test suite for natural units (`tests/ucon/test_natural_units.py`)
+
### Fixed
- MCP session state persistence across tool calls (#209)
@@ -418,7 +438,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
- Initial commit
-[Unreleased]: https://github.com/withtwoemms/ucon/compare/0.9.2...HEAD
+[Unreleased]: https://github.com/withtwoemms/ucon/compare/0.9.3...HEAD
+[0.9.3]: https://github.com/withtwoemms/ucon/compare/0.9.2...0.9.3
[0.9.2]: https://github.com/withtwoemms/ucon/compare/0.9.1...0.9.2
[0.9.1]: https://github.com/withtwoemms/ucon/compare/0.9.0...0.9.1
[0.9.0]: https://github.com/withtwoemms/ucon/compare/0.8.5...0.9.0
diff --git a/README.md b/README.md
index a42c4fd..1f335bf 100644
--- a/README.md
+++ b/README.md
@@ -128,11 +128,14 @@ AI agents can then convert units, check dimensions, and perform factor-label cal
## Features
- **Physical constants** — CODATA 2022 values with uncertainty propagation (`E = m * c**2`)
+- **Custom constants** — Define domain-specific constants with uncertainty propagation via MCP or Python API
- **String parsing** — `parse("9.81 m/s^2")` with uncertainty support (`1.234 ± 0.005 m`)
- **Dimensional algebra** — Units combine through multiplication/division with automatic dimension tracking
- **Scale prefixes** — Full SI (kilo, milli, micro, etc.) and binary (kibi, mebi) prefix support
- **Uncertainty propagation** — Errors propagate through arithmetic and conversions
- **Pseudo-dimensions** — Semantically isolated handling of angles, ratios, and counts
+- **Natural units** — Custom dimensional bases where c=ℏ=k_B=1 for particle physics
+- **Logarithmic units** — dB, pH, and neper conversions with uncertainty propagation
- **Pydantic v2** — Type-safe API validation and JSON serialization
- **MCP server** — AI agent integration with Claude, Cursor, and other MCP clients
- **ConversionGraph** — Extensible conversion registry with custom unit support
@@ -149,7 +152,7 @@ AI agents can then convert units, check dimensions, and perform factor-label cal
| **0.6.x** | Pydantic + MCP Server | Complete |
| **0.7.x** | Compute Tool + Extension API | Complete |
| **0.8.x** | Basis Abstraction + String Parsing | Complete |
-| **0.9.x** | Constants + Logarithmic Units | In Progress |
+| **0.9.x** | Constants + Natural Units | Complete |
| **0.10.x** | NumPy/Polars Integration | Planned |
| **1.0.0** | API Stability | Planned |
diff --git a/ROADMAP.md b/ROADMAP.md
index 32a9956..c10b3cb 100644
--- a/ROADMAP.md
+++ b/ROADMAP.md
@@ -46,23 +46,24 @@ ucon is a dimensional analysis library for engineers building systems where unit
| v0.9.0 | Physical Constants | Complete |
| v0.9.1 | Logarithmic Units | Complete |
| v0.9.2 | MCP Constants Tools | Complete |
-| v0.9.3 | Natural Units | Planned |
+| v0.9.3 | Natural Units + MCP Session Fixes | Complete |
| v0.10.0 | Scientific Computing | Planned |
| v1.0.0 | API Stability | Planned |
---
-## Current Version: **v0.9.2** (complete)
+## Current Version: **v0.9.3** (complete)
-Building on v0.9.1 baseline:
-- `ucon.basis` (`Basis`, `BasisComponent`, `Vector`, `BasisTransform`, `BasisGraph`)
-- `ucon.bases` (standard bases: `SI`, `CGS`, `CGS_ESU`; standard transforms)
+Building on v0.9.2 baseline:
+- `ucon.basis` (`Basis`, `BasisComponent`, `Vector`, `BasisTransform`, `BasisGraph`, `ConstantAwareBasisTransform`)
+- `ucon.bases` (standard bases: `SI`, `CGS`, `CGS_ESU`, `NATURAL`; standard transforms including `SI_TO_NATURAL`)
- `ucon.dimension` (`Dimension` as frozen dataclass backed by basis-aware `Vector`)
- `ucon.core` (`Scale`, `Unit`, `UnitFactor`, `UnitProduct`, `Number`, `Ratio`, `UnitSystem`, `RebasedUnit`, `Exponent`)
- `ucon.maps` (`Map`, `LinearMap`, `AffineMap`, `ComposedMap`, `LogMap`, `ExpMap`)
- `ucon.graph` (`ConversionGraph`, default graph, `get_default_graph()`, `using_graph()`, cross-basis conversion)
- `ucon.units` (SI + imperial + information + angle + ratio units, callable syntax, `si` and `imperial` systems, `get_unit_by_name()`)
- `ucon.pydantic` (`Number` type for Pydantic v2 models)
+- `ucon.mcp` (`SessionState`, `DefaultSessionState` for injectable session management)
- Callable unit API: `meter(5)`, `(mile / hour)(60)`
- `Number.simplify()` for base-scale normalization
- Pseudo-dimensions: `ANGLE`, `SOLID_ANGLE`, `RATIO`, `COUNT` with semantic isolation
@@ -77,6 +78,8 @@ Building on v0.9.1 baseline:
- Physical constants: `Constant` class with CODATA 2022 values and uncertainty propagation
- Logarithmic units: pH with concentration dimension, dBm, dBW, dBV, dBSPL
- MCP constants tools: `list_constants()`, `define_constant()` for AI agent access
+- Natural units: `NATURAL` basis with c=ℏ=k_B=1, `ConstantAwareBasisTransform` for non-square transforms
+- MCP session persistence: lifespan-scoped session state across tool calls
---
@@ -628,18 +631,27 @@ Prerequisite for factor-label chains with countable items (tablets, doses).
---
-## v0.9.3 — Natural Units
+## v0.9.3 — Natural Units + MCP Session Fixes (Complete)
-**Theme:** Custom dimensional bases where c=ℏ=1.
+**Theme:** Custom dimensional bases where c=ℏ=1, plus MCP reliability fixes.
-- [ ] `ConstantAwareBasisTransform` with `inverse()` for non-square transforms
-- [ ] `NATURAL` basis with single energy dimension
-- [ ] SI ↔ NATURAL transform using constant bindings
-- [ ] Documentation and examples
+### Natural Units
+
+- [x] `ConstantAwareBasisTransform` with `inverse()` for non-square transforms
+- [x] `ConstantBinding` dataclass for tracking which constants absorb dimensions
+- [x] `NATURAL` basis with single energy dimension
+- [x] SI → NATURAL transform: T→E⁻¹, L→E⁻¹, M→E, Θ→E (via ℏ, c, k_B)
+- [x] NATURAL → SI inverse transform
+- [x] Documentation and examples (`examples/basis/`)
+
+### MCP Safety Improvements (Feedback Issues)
+
+- [x] All session-dependent tools use FastMCP `Context` injection
**Outcomes:**
- Natural units leverage custom basis infrastructure
- Foundation for particle physics and quantum field theory domains
+- MCP tools are more reliable for multi-call agent workflows
---
diff --git a/docs/architecture/kind-of-quantity.md b/docs/architecture/kind-of-quantity.md
new file mode 100644
index 0000000..60a066f
--- /dev/null
+++ b/docs/architecture/kind-of-quantity.md
@@ -0,0 +1,344 @@
+# The Kind-of-Quantity Problem
+
+Core conceptual foundation for understanding dimensional ambiguity and how ucon addresses it.
+
+---
+
+## The Problem
+
+When working with physical quantities, a fundamental limitation of dimensional
+analysis becomes apparent: **units and dimensions do not uniquely identify
+the kind of quantity being measured.**
+
+### Examples of Ambiguity
+
+| Unit | Possible Quantities |
+|------|---------------------|
+| `kg·m²·s⁻²` | Energy, Torque, Work, Heat |
+| `s⁻¹` | Frequency, Angular velocity, Radioactive activity |
+| `1` (dimensionless) | Angle, Refractive index, Efficiency, Count |
+
+This is known as the **Kind-of-Quantity (KOQ) problem** in metrology literature.
+
+### Why It Matters
+
+- **Silent errors**: Adding energy to torque produces nonsense but passes dimensional checks
+- **Ambiguous conversions**: `s⁻¹ → Hz` is valid for frequency but not for activity
+- **Lost semantics**: Dimensionless quantities lose all physical meaning
+
+---
+
+## ucon's Two Solutions
+
+ucon addresses KOQ through two complementary mechanisms:
+
+1. **Pseudo-dimensions** — Semantic tags for quantities dimensionless within a basis
+2. **Basis abstraction** — Coordinate transforms across dimensional systems
+
+Both solutions build on the relationship between `BasisVector` and `Dimension`.
+
+---
+
+## BasisVector and Dimension
+
+Before diving into solutions, it's important to understand how ucon represents
+dimensions internally.
+
+### BasisVector: The Coordinate Representation
+
+A `BasisVector` is the raw numerical representation of a dimension—a tuple of
+exponents over a basis. Think of it as coordinates in dimensional space:
+
+```python
+from ucon import BasisVector, SI
+from fractions import Fraction
+
+# Velocity = L/T = L¹·T⁻¹ in SI
+# SI order: (T, L, M, I, Θ, J, N, B)
+velocity_vector = BasisVector(SI, (
+ Fraction(-1), # T⁻¹
+ Fraction(1), # L¹
+ Fraction(0), # M⁰
+ Fraction(0), # I⁰
+ Fraction(0), # Θ⁰
+ Fraction(0), # J⁰
+ Fraction(0), # N⁰
+ Fraction(0), # B⁰
+))
+
+velocity_vector.is_dimensionless() # False
+```
+
+`BasisVector` supports algebraic operations (multiplication, division, exponentiation)
+that add, subtract, or scale exponents:
+
+```python
+# length_vec * length_vec = area_vec (exponents add)
+# length_vec / time_vec = velocity_vec (exponents subtract)
+```
+
+### Dimension: The Named, Semantic Wrapper
+
+A `Dimension` wraps a `BasisVector` with additional metadata:
+
+- **name**: Human-readable identifier (e.g., "velocity", "energy")
+- **symbol**: Short notation (e.g., "L", "T")
+- **tag**: For pseudo-dimensions (e.g., "angle", "count")
+
+```python
+from ucon import Dimension
+
+# Access predefined dimensions
+Dimension.velocity.name # "velocity"
+Dimension.velocity.vector # The underlying BasisVector
+
+# Dimensions resolve from vectors automatically
+from ucon import resolve_dimension
+resolved = resolve_dimension(velocity_vector)
+resolved == Dimension.velocity # True
+```
+
+### The Key Relationship
+
+```mermaid
+flowchart LR
+ subgraph BasisVector["BasisVector"]
+ bv_basis["basis"]
+ bv_components["components"]
+ bv_algebra["algebra"]
+ end
+
+ subgraph Dimension["Dimension"]
+ d_vector["vector"]
+ d_name["name"]
+ d_symbol["symbol"]
+ d_tag["tag (pseudo)"]
+ end
+
+ Dimension -- "wraps" --> BasisVector
+
+ BasisVector --> coords["(−1, 1, 0, 0, 0, 0, 0, 0)
raw coordinates in dimensional space"]
+ Dimension --> semantic["Dimension(velocity)
semantic identity for human reasoning"]
+```
+
+**Why does this matter for KOQ?**
+
+- `BasisVector` handles the algebra—it knows nothing about "torque" vs "energy"
+- `Dimension` provides the semantic layer—it can distinguish quantities via names or tags
+- The KOQ problem arises when the same `BasisVector` maps to multiple physical meanings
+- ucon's solutions work by either:
+ - Adding a **tag** to the `Dimension` (pseudo-dimensions)
+ - Choosing a different **basis** so the `BasisVector` differs (basis abstraction)
+
+---
+
+## Solution 1: Pseudo-Dimensions
+
+For quantities that are mathematically dimensionless but physically distinct,
+ucon provides **pseudo-dimensions**:
+
+```python
+from ucon import Dimension, units
+
+# These are all dimensionless (zero vector) but semantically isolated
+Dimension.angle # radians, degrees
+Dimension.solid_angle # steradians
+Dimension.ratio # percent, ppm, ppb
+Dimension.count # discrete items
+
+# Semantic isolation prevents nonsensical operations
+angle = units.radian(3.14)
+count = units.each(5)
+
+# This raises TypeError — incompatible pseudo-dimensions
+result = angle + count # Error!
+
+# But algebra within the same pseudo-dimension works
+total_angle = units.radian(1.57) + units.radian(1.57) # OK
+```
+
+### How Pseudo-Dimensions Work
+
+| Aspect | Behavior |
+|--------|----------|
+| Algebraic | Zero vector — behaves as dimensionless in multiplication |
+| Semantic | Tagged — prevents mixing incompatible types |
+| Cancellation | `angle / angle = dimensionless` (not angle) |
+
+### When to Use
+
+- Angles (radian, degree, steradian)
+- Counts of discrete items
+- Ratios and proportions (percent, ppm)
+- Any quantity dimensionless in all reasonable unit systems
+
+---
+
+## Solution 2: Basis Abstraction
+
+The key insight is that **dimensional ambiguity is basis-dependent**. Two quantities
+that share the same dimension in one basis may be distinguished in another.
+
+### Example: Torque vs Energy
+
+In standard SI, torque and energy are dimensionally identical:
+
+| Quantity | SI Dimension | Problem |
+|----------|--------------|---------|
+| Energy | M·L²·T⁻² | Same as torque! |
+| Torque | M·L²·T⁻² | Same as energy! |
+
+But physically, torque = force × lever arm × **angle**. If we choose a basis where
+angle is NOT dimensionless, the ambiguity resolves:
+
+```python
+from ucon import Basis, BasisComponent, BasisVector
+from fractions import Fraction
+
+# Define a basis where angle (θ) is a base dimension
+MECHANICAL_WITH_ANGLE = Basis(
+ "mechanical_angle",
+ [
+ BasisComponent("time", "T"),
+ BasisComponent("length", "L"),
+ BasisComponent("mass", "M"),
+ BasisComponent("angle", "θ"), # Angle as base dimension!
+ ],
+)
+
+# In this basis:
+# Energy = M·L²·T⁻²·θ⁰ (no angle dependence)
+# Torque = M·L²·T⁻²·θ¹ (depends on angle!)
+
+energy_vec = BasisVector(MECHANICAL_WITH_ANGLE, (
+ Fraction(-2), Fraction(2), Fraction(1), Fraction(0) # T⁻² L² M¹ θ⁰
+))
+
+torque_vec = BasisVector(MECHANICAL_WITH_ANGLE, (
+ Fraction(-2), Fraction(2), Fraction(1), Fraction(1) # T⁻² L² M¹ θ¹
+))
+
+energy_vec == torque_vec # False! Now distinguishable
+```
+
+### Example: Gray (Gy) vs Sievert (Sv)
+
+In radiation dosimetry, absorbed dose (Gy) and dose equivalent (Sv) have the
+same SI dimension (m²·s⁻²) but measure different things:
+
+- **Gray**: Physical energy absorbed per unit mass
+- **Sievert**: Biological effect (weighted by radiation type)
+
+A domain-specific basis can distinguish them:
+
+```python
+from ucon import Basis, BasisComponent, BasisVector
+from fractions import Fraction
+
+DOSIMETRY = Basis(
+ "dosimetry",
+ [
+ BasisComponent("time", "T"),
+ BasisComponent("length", "L"),
+ BasisComponent("mass", "M"),
+ BasisComponent("biological_weight", "Q"), # Quality factor dimension
+ ],
+)
+
+# Absorbed dose (Gy): L²·T⁻²·Q⁰
+# Dose equivalent (Sv): L²·T⁻²·Q¹
+
+absorbed_dose = BasisVector(DOSIMETRY, (
+ Fraction(-2), Fraction(2), Fraction(0), Fraction(0)
+))
+
+dose_equivalent = BasisVector(DOSIMETRY, (
+ Fraction(-2), Fraction(2), Fraction(0), Fraction(1) # Q¹ factor!
+))
+
+# Now they're dimensionally distinct
+absorbed_dose == dose_equivalent # False
+```
+
+### Natural Units (Advanced)
+
+For particle physics, where c = ℏ = k_B = 1, dimensions collapse further:
+
+| Quantity | SI Dimension | Natural Units |
+|----------|--------------|---------------|
+| Velocity | L·T⁻¹ | E⁰ (dimensionless!) |
+| Length | L | E⁻¹ |
+| Mass | M | E |
+
+```python
+from ucon import SI, NATURAL, SI_TO_NATURAL, BasisVector
+from fractions import Fraction
+
+# Velocity becomes dimensionless in natural units
+si_velocity = BasisVector(SI, (
+ Fraction(-1), Fraction(1), Fraction(0), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0),
+))
+
+natural_velocity = SI_TO_NATURAL(si_velocity)
+natural_velocity.is_dimensionless() # True! (c = 1)
+```
+
+### When to Use Basis Abstraction
+
+| Scenario | Basis Strategy |
+|----------|----------------|
+| Torque vs Energy | Add angle as base dimension |
+| Gy vs Sv | Add biological weighting dimension |
+| SI ↔ CGS | Standard basis transforms |
+| Particle physics | Natural units (E only) |
+| Electromagnetism | CGS-ESU (charge as base) |
+
+---
+
+## Comparing the Two Approaches
+
+| Aspect | Pseudo-Dimensions | Basis Abstraction |
+|--------|-------------------|-------------------|
+| **Problem solved** | Dimensionless quantities that differ semantically | Quantities with same dimension that differ physically |
+| **Mechanism** | Semantic tag on zero-vector | Choose a basis where dimensions differ |
+| **Scope** | Works within any single basis | Requires defining/choosing appropriate basis |
+| **Examples** | angle vs count vs ratio | torque vs energy, Gy vs Sv |
+| **Declaration** | `Dimension.pseudo("tag")` | `Basis(name, components)` |
+| **Trade-off** | Simple but limited to "truly" dimensionless | Powerful but requires domain modeling |
+
+---
+
+## The Unified Picture
+
+Both solutions preserve **kind-of-quantity information** that pure dimensional
+analysis loses. The choice depends on the nature of the ambiguity:
+
+```mermaid
+flowchart TB
+ KOQ["KOQ Ambiguity"]
+
+ KOQ --> dimensionless["Dimensionless but
semantically distinct"]
+ KOQ --> same_dim["Same dimension but
physically distinct"]
+
+ dimensionless --> pseudo["PSEUDO-DIMENSION
angle
count
ratio"]
+ same_dim --> basis["BASIS ABSTRACTION
Expand basis to
include hidden
dimension"]
+
+ pseudo --> tag["Tag on zero-vector
(0,0,0,...,'angle')"]
+ basis --> vectors["Different vectors
in richer basis"]
+```
+
+**Key insight**: The SI's 7-dimensional basis is a pragmatic choice, not a
+physical necessity. Some KOQ problems arise because the SI basis is too
+coarse to distinguish certain quantities. ucon lets you:
+
+1. **Add semantic tags** for quantities the SI treats as dimensionless
+2. **Define richer bases** for quantities the SI conflates
+
+---
+
+## Further Reading
+
+- Hall, B.D. (2022). "Representing quantities and units in digital systems." *Measurement: Sensors*
+- Hall, B.D. (2023). "Modelling expressions of physical quantities." *IC3K 2023*
+- Hall, B.D. (2025). "Interpreting SI notation for units and dimensions." *Measurement: Sensors*
diff --git a/docs/guides/index.md b/docs/guides/index.md
index f052496..54c4d70 100644
--- a/docs/guides/index.md
+++ b/docs/guides/index.md
@@ -17,6 +17,7 @@ Task-oriented guides for common use cases.
- **[Dimensional Analysis](dimensional-analysis.md)** - Step-by-step factor-label calculations
- **[Custom Units & Graphs](custom-units-and-graphs.md)** - Define domain-specific units
+- **[Natural Units](natural-units.md)** - Particle physics natural units (c = h_bar = k_B = 1)
## Domain Walkthroughs
diff --git a/docs/guides/natural-units.md b/docs/guides/natural-units.md
new file mode 100644
index 0000000..f3b4375
--- /dev/null
+++ b/docs/guides/natural-units.md
@@ -0,0 +1,271 @@
+# Natural Units Guide
+
+This guide explains ucon's support for natural units, commonly used in particle physics and theoretical physics.
+
+---
+
+## What Are Natural Units?
+
+Natural units simplify physical equations by setting certain fundamental constants to 1:
+
+| Convention | Constants Set to 1 | Common Use |
+|------------|-------------------|------------|
+| **Particle Physics** | c, h_bar, k_B | High-energy physics |
+| **Planck Units** | c, h_bar, G, k_B | Quantum gravity |
+| **Stoney Units** | c, e, G, k_B | Historical |
+
+In ucon's `NATURAL` basis (particle physics convention), setting **c = h_bar = k_B = 1** collapses length, time, mass, and temperature into expressions of a single dimension: **energy**.
+
+---
+
+## Key Consequences
+
+### Velocity is Dimensionless
+
+With c = 1, the speed of light is just the number 1. All velocities are expressed as fractions of c:
+
+```python
+from ucon import SI, NATURAL, SI_TO_NATURAL, BasisVector
+from fractions import Fraction
+
+# Velocity in SI: L^1 T^-1
+si_velocity = BasisVector(SI, (
+ Fraction(-1), # T^-1
+ Fraction(1), # L^1
+ Fraction(0), Fraction(0), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0),
+))
+
+natural_velocity = SI_TO_NATURAL(si_velocity)
+natural_velocity.is_dimensionless() # True!
+```
+
+### Mass and Energy Share Dimension
+
+With E = mc^2 and c = 1, mass **is** energy:
+
+```python
+# Mass in SI: M^1
+si_mass = BasisVector(SI, (
+ Fraction(0), Fraction(0), Fraction(1), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0),
+))
+
+natural_mass = SI_TO_NATURAL(si_mass)
+natural_mass["E"] # Fraction(1) - pure energy!
+```
+
+Particle physicists express masses in energy units: "The electron mass is 0.511 MeV."
+
+### Action is Dimensionless
+
+With h_bar = 1, action (energy x time) becomes dimensionless:
+
+```python
+# Action in SI: M^1 L^2 T^-1 (same as h_bar)
+si_action = BasisVector(SI, (
+ Fraction(-1), # T^-1
+ Fraction(2), # L^2
+ Fraction(1), # M^1
+ Fraction(0), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0),
+))
+
+natural_action = SI_TO_NATURAL(si_action)
+natural_action.is_dimensionless() # True!
+```
+
+### Length and Time are Inverse Energy
+
+Since L = h_bar*c/E and T = h_bar/E, both length and time scale as E^-1:
+
+```python
+# Length: L -> E^-1
+# Time: T -> E^-1
+
+# This means frequency (T^-1) has dimension E:
+si_frequency = BasisVector(SI, (
+ Fraction(-1), Fraction(0), Fraction(0), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0),
+))
+SI_TO_NATURAL(si_frequency)["E"] # Fraction(1)
+```
+
+---
+
+## The NATURAL Basis
+
+ucon provides the `NATURAL` basis with a single component:
+
+```python
+from ucon import NATURAL
+
+len(NATURAL) # 1
+NATURAL[0].name # "energy"
+NATURAL[0].symbol # "E"
+```
+
+### Dimension Mapping Table
+
+| SI Dimension | Natural | Physical Origin |
+|--------------|---------|-----------------|
+| Length (L) | E^-1 | L = h_bar*c / E |
+| Time (T) | E^-1 | T = h_bar / E |
+| Mass (M) | E^1 | E = mc^2 |
+| Temperature (Theta) | E^1 | E = k_B*T |
+| Current (I) | -- | Not representable |
+| Luminous Int. (J) | -- | Not representable |
+| Amount (N) | -- | Not representable |
+| Information (B) | -- | Not representable |
+
+---
+
+## Using SI_TO_NATURAL
+
+Transform SI dimensional vectors to natural units:
+
+```python
+from ucon import SI, SI_TO_NATURAL, BasisVector
+from fractions import Fraction
+
+# Energy in SI: M^1 L^2 T^-2
+si_energy = BasisVector(SI, (
+ Fraction(-2), Fraction(2), Fraction(1), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0),
+))
+
+natural_energy = SI_TO_NATURAL(si_energy)
+print(f"Energy: E^{natural_energy['E']}") # "Energy: E^1"
+```
+
+### Handling Electromagnetic Dimensions
+
+Electromagnetic quantities (current, charge, etc.) are **not representable** in pure natural units. Attempting to transform them raises `LossyProjection`:
+
+```python
+from ucon import SI, SI_TO_NATURAL, BasisVector, LossyProjection
+from fractions import Fraction
+
+# Current: I^1
+si_current = BasisVector(SI, (
+ Fraction(0), Fraction(0), Fraction(0), Fraction(1),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0),
+))
+
+# This raises LossyProjection
+try:
+ SI_TO_NATURAL(si_current)
+except LossyProjection as e:
+ print("Cannot transform current to natural units")
+```
+
+If you want to proceed anyway (projecting to zero):
+
+```python
+result = SI_TO_NATURAL(si_current, allow_projection=True)
+result.is_dimensionless() # True (current dimension is lost)
+```
+
+---
+
+## Going Back: NATURAL_TO_SI
+
+The inverse transform uses constant bindings to reverse the mapping:
+
+```python
+from ucon import NATURAL, NATURAL_TO_SI, BasisVector
+from fractions import Fraction
+
+# Energy in natural units
+natural_energy = BasisVector(NATURAL, (Fraction(1),))
+
+# Transform to SI representation
+si_result = NATURAL_TO_SI(natural_energy)
+print(si_result.basis.name) # "SI"
+```
+
+Note: The inverse transform recovers the **primary** SI representation. Since mass, temperature, and energy all have E^1, the inverse picks one canonical form.
+
+---
+
+## Particle Physics Examples
+
+### Cross-Section (Area)
+
+Cross-sections are measured in units of area (L^2), which becomes E^-2 in natural units:
+
+```python
+si_cross_section = BasisVector(SI, (
+ Fraction(0), Fraction(2), Fraction(0), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0),
+))
+
+natural_xs = SI_TO_NATURAL(si_cross_section)
+print(f"Cross-section: E^{natural_xs['E']}") # "E^-2"
+```
+
+In particle physics, cross-sections are often given in inverse GeV^2 (or barns).
+
+### Decay Width
+
+Decay widths have dimension of inverse time (T^-1), which is E in natural units:
+
+```python
+si_decay_width = BasisVector(SI, (
+ Fraction(-1), Fraction(0), Fraction(0), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0),
+))
+
+natural_width = SI_TO_NATURAL(si_decay_width)
+print(f"Decay width: E^{natural_width['E']}") # "E^1"
+```
+
+Particle physicists express decay widths directly in energy units (MeV, GeV).
+
+### Momentum
+
+Momentum (M*L/T) has dimension E in natural units:
+
+```python
+si_momentum = BasisVector(SI, (
+ Fraction(-1), Fraction(1), Fraction(1), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0),
+))
+
+natural_p = SI_TO_NATURAL(si_momentum)
+print(f"Momentum: E^{natural_p['E']}") # "E^1"
+```
+
+Hence "the proton momentum is 100 GeV" -- momentum and energy share units.
+
+---
+
+## Advanced: ConstantBinding
+
+Under the hood, `SI_TO_NATURAL` uses `ConstantBinding` to record which constants relate each dimension:
+
+```python
+from ucon.bases import SI_TO_NATURAL
+
+# The transform has 4 bindings (L, T, M, Theta)
+len(SI_TO_NATURAL.bindings) # 4
+
+for binding in SI_TO_NATURAL.bindings:
+ print(f"{binding.source_component.symbol} via {binding.constant_symbol}^{binding.exponent}")
+```
+
+These bindings enable the inverse transform to work despite the non-square (8x1) matrix.
+
+---
+
+## Summary
+
+| Feature | Description |
+|---------|-------------|
+| `NATURAL` | Single-dimensional basis (energy only) |
+| `SI_TO_NATURAL` | Transform SI dimensions to natural units |
+| `NATURAL_TO_SI` | Inverse transform using constant bindings |
+| `ConstantBinding` | Records dimension-constant relationships |
+| `LossyProjection` | Raised for non-representable dimensions |
+
+Natural units in ucon enable dimensional analysis for particle physics, where c = h_bar = k_B = 1 simplifies calculations and all quantities reduce to powers of energy.
diff --git a/docs/reference/api.md b/docs/reference/api.md
index 876fa67..a148790 100644
--- a/docs/reference/api.md
+++ b/docs/reference/api.md
@@ -422,6 +422,28 @@ c_num = c.as_number()
type(c_num) # Number
```
+### Creating Custom Constants
+
+```python
+from ucon.constants import Constant
+from ucon import units
+
+# Domain-specific constant
+custom_k = Constant(
+ symbol="k_custom",
+ name="custom spring constant",
+ value=250.0,
+ unit=units.newton / units.meter,
+ uncertainty=0.5, # or None for exact
+ source="Lab measurement",
+ category="session",
+)
+
+# Use in calculations
+x = units.meter(0.1)
+F = custom_k.as_number() * x # Hooke's law
+```
+
---
## ConversionGraph
@@ -734,3 +756,190 @@ reset_default_basis_graph()
| `using_basis_graph(graph)` | Scoped BasisGraph override |
| `set_default_basis_graph(graph)` | Replace module default |
| `reset_default_basis_graph()` | Restore standard graph |
+
+---
+
+## Natural Units
+
+Support for natural units where physical constants c = h_bar = k_B = 1.
+
+```python
+from ucon import NATURAL, SI_TO_NATURAL, NATURAL_TO_SI
+from ucon.basis import ConstantBinding, ConstantAwareBasisTransform
+```
+
+### NATURAL Basis
+
+Single-dimensional basis where all quantities reduce to powers of energy:
+
+```python
+from ucon import NATURAL, BasisVector
+from fractions import Fraction
+
+# NATURAL has one component: energy (E)
+len(NATURAL) # 1
+NATURAL[0].name # "energy"
+NATURAL[0].symbol # "E"
+
+# Create natural unit vectors
+energy = BasisVector(NATURAL, (Fraction(1),)) # E^1
+inverse_energy = BasisVector(NATURAL, (Fraction(-1),)) # E^-1
+dimensionless = BasisVector(NATURAL, (Fraction(0),)) # E^0
+```
+
+### Dimension Mappings
+
+With c = h_bar = k_B = 1:
+
+| SI Dimension | Natural Units | Physical Meaning |
+|--------------|---------------|------------------|
+| Length (L) | E^-1 | L = h_bar*c / E |
+| Time (T) | E^-1 | T = h_bar / E |
+| Mass (M) | E | E = mc^2 |
+| Temperature (Theta) | E | E = k_B * T |
+| Velocity (L/T) | E^0 | Dimensionless (c=1) |
+| Action (ML^2/T) | E^0 | Dimensionless (h_bar=1) |
+| Energy (ML^2/T^2) | E | Energy is fundamental |
+
+### SI_TO_NATURAL Transform
+
+Transform SI dimensions to natural units:
+
+```python
+from ucon import SI, NATURAL, SI_TO_NATURAL, BasisVector
+from fractions import Fraction
+
+# Velocity: L^1 T^-1 in SI
+si_velocity = BasisVector(SI, (
+ Fraction(-1), # T^-1
+ Fraction(1), # L^1
+ Fraction(0), Fraction(0), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0),
+))
+
+natural_velocity = SI_TO_NATURAL(si_velocity)
+natural_velocity.is_dimensionless() # True! (c = 1)
+
+# Energy: M^1 L^2 T^-2 in SI
+si_energy = BasisVector(SI, (
+ Fraction(-2), # T^-2
+ Fraction(2), # L^2
+ Fraction(1), # M^1
+ Fraction(0), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0),
+))
+
+natural_energy = SI_TO_NATURAL(si_energy)
+natural_energy["E"] # Fraction(1) - pure energy
+```
+
+### Lossy Projections
+
+Electromagnetic dimensions (current, luminous intensity, etc.) are not representable in pure natural units:
+
+```python
+from ucon import SI, SI_TO_NATURAL, BasisVector, LossyProjection
+from fractions import Fraction
+
+# Current has no representation in natural units
+si_current = BasisVector(SI, (
+ Fraction(0), Fraction(0), Fraction(0), Fraction(1),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0),
+))
+
+# Raises LossyProjection
+try:
+ SI_TO_NATURAL(si_current)
+except LossyProjection as e:
+ print(e) # "current has no representation in natural"
+
+# Use allow_projection=True to drop
+result = SI_TO_NATURAL(si_current, allow_projection=True)
+result.is_dimensionless() # True (projected to zero)
+```
+
+### NATURAL_TO_SI (Inverse)
+
+Transform natural unit dimensions back to SI using constant bindings:
+
+```python
+from ucon import NATURAL, NATURAL_TO_SI, BasisVector
+from fractions import Fraction
+
+# Energy in natural units
+natural_energy = BasisVector(NATURAL, (Fraction(1),))
+
+# Transform back to SI
+si_result = NATURAL_TO_SI(natural_energy)
+si_result.basis.name # "SI"
+```
+
+### ConstantBinding
+
+Records the relationship between source and target dimensions via physical constants:
+
+```python
+from ucon import ConstantBinding, NATURAL, BasisVector
+from ucon.basis import BasisComponent
+from fractions import Fraction
+
+# Length relates to inverse energy via h_bar*c
+binding = ConstantBinding(
+ source_component=BasisComponent("length", "L"),
+ target_expression=BasisVector(NATURAL, (Fraction(-1),)),
+ constant_symbol="h_bar_c", # h_bar * c
+ exponent=Fraction(1),
+)
+
+binding.constant_symbol # "h_bar_c"
+binding.exponent # Fraction(1)
+```
+
+### ConstantAwareBasisTransform
+
+Basis transform with bindings that enable inversion of non-square matrices:
+
+```python
+from ucon import ConstantAwareBasisTransform
+
+# SI_TO_NATURAL is a ConstantAwareBasisTransform (8x1 matrix)
+isinstance(SI_TO_NATURAL, ConstantAwareBasisTransform) # True
+
+# Inverse works despite non-square matrix
+NATURAL_TO_SI = SI_TO_NATURAL.inverse()
+NATURAL_TO_SI.source.name # "natural"
+NATURAL_TO_SI.target.name # "SI"
+
+# Convert to plain BasisTransform (loses binding info)
+plain = SI_TO_NATURAL.as_basis_transform()
+```
+
+### Particle Physics Applications
+
+Common particle physics dimensions:
+
+```python
+from ucon import SI, SI_TO_NATURAL, BasisVector
+from fractions import Fraction
+
+# Cross-section (area): L^2 -> E^-2
+si_cross_section = BasisVector(SI, (
+ Fraction(0), Fraction(2), Fraction(0), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0),
+))
+SI_TO_NATURAL(si_cross_section)["E"] # Fraction(-2)
+
+# Decay width (inverse time): T^-1 -> E
+si_decay_width = BasisVector(SI, (
+ Fraction(-1), Fraction(0), Fraction(0), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0),
+))
+SI_TO_NATURAL(si_decay_width)["E"] # Fraction(1)
+
+# Momentum: M*L/T -> E
+si_momentum = BasisVector(SI, (
+ Fraction(-1), Fraction(1), Fraction(1), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0),
+))
+SI_TO_NATURAL(si_momentum)["E"] # Fraction(1)
+```
diff --git a/examples/basis/__init__.py b/examples/basis/__init__.py
new file mode 100644
index 0000000..cfee5fc
--- /dev/null
+++ b/examples/basis/__init__.py
@@ -0,0 +1 @@
+# Basis transform examples
diff --git a/examples/basis/elemental_units_demo.py b/examples/basis/elemental_units_demo.py
new file mode 100644
index 0000000..4a93a5b
--- /dev/null
+++ b/examples/basis/elemental_units_demo.py
@@ -0,0 +1,463 @@
+#!/usr/bin/env python3
+"""
+ConstantAwareBasisTransform Demo: Classical Elemental Units
+
+A fantasy unit system based on classical elements where dimensions are
+FUNDAMENTALLY different from SI - not just scaled versions.
+
+The Elemental Basis:
+--------------------
+In the world of Alchemia, natural philosophers discovered that reality
+is composed of five fundamental essences:
+
+ - Fire (Fi): The principle of transformation and energy release
+ - Water (Wa): The principle of flow and continuity
+ - Earth (Ea): The principle of substance and permanence
+ - Air (Ai): The principle of motion and change
+ - Aether (Ae): The quintessence - purely magical, has NO physical analog
+
+Key Features:
+-------------
+1. CONVOLUTIONS: Elemental dimensions map to COMBINATIONS of SI dimensions:
+ - Fire is not just energy, it's POWER (energy/time) - transformation rate
+ - Water is not volume, it's FLOW (volume/time) - continuous motion
+ - Air is not length, it's VELOCITY (length/time) - change itself
+
+2. NON-SI DIMENSION: Aether has no physical representation and raises
+ LossyProjection when transformed to SI - it's purely magical!
+
+3. NON-SQUARE MATRIX: 5 elemental dimensions -> 4 SI dimensions
+
+This demonstrates ConstantAwareBasisTransform with true dimensional complexity.
+"""
+
+from fractions import Fraction
+
+from ucon.basis import (
+ Basis,
+ BasisComponent,
+ ConstantAwareBasisTransform,
+ ConstantBinding,
+ LossyProjection,
+ Vector,
+)
+
+
+def main():
+ print("=" * 70)
+ print(" THE ELEMENTAL BASIS: Classical Alchemy Meets Dimensional Analysis")
+ print("=" * 70)
+
+ # =========================================================================
+ # PART 1: Define the Elemental and SI Bases
+ # =========================================================================
+ print("\n" + "=" * 70)
+ print("PART 1: The Five Elements")
+ print("=" * 70)
+
+ ELEMENTAL = Basis(
+ "Elemental",
+ [
+ BasisComponent("fire", "Fi"), # Transformation/power
+ BasisComponent("water", "Wa"), # Flow
+ BasisComponent("earth", "Ea"), # Substance
+ BasisComponent("air", "Ai"), # Motion
+ BasisComponent("aether", "Ae"), # The quintessence (magical)
+ ],
+ )
+
+ SI_MECH = Basis(
+ "SI-Mechanical",
+ [
+ BasisComponent("length", "L"),
+ BasisComponent("mass", "M"),
+ BasisComponent("time", "T"),
+ BasisComponent("temperature", "K"), # For heat/fire concepts
+ ],
+ )
+
+ print(f"""
+The Elemental Basis ({len(ELEMENTAL)} dimensions):
+
+ Fire (Fi): "That which transforms"
+ Not mere energy, but the RATE of transformation
+ Maps to: POWER = M * L^2 * T^-3
+
+ Water (Wa): "That which flows"
+ Not volume, but continuous FLOW
+ Maps to: VOLUMETRIC FLOW = L^3 * T^-1
+
+ Earth (Ea): "That which endures"
+ Pure substance, mass itself
+ Maps to: MASS = M
+
+ Air (Ai): "That which moves"
+ Not distance, but VELOCITY - motion itself
+ Maps to: VELOCITY = L * T^-1
+
+ Aether (Ae): "The quintessence"
+ Purely magical, exists beyond physical reality
+ Maps to: NOTHING (raises LossyProjection!)
+
+SI-Mechanical Basis ({len(SI_MECH)} dimensions): L, M, T, K
+Matrix shape: {len(ELEMENTAL)} x {len(SI_MECH)} = 5 x 4 (non-square!)
+""")
+
+ # =========================================================================
+ # PART 2: Define the Alchemical Constants
+ # =========================================================================
+ print("=" * 70)
+ print("PART 2: The Alchemical Constants")
+ print("=" * 70)
+
+ print("""
+The great alchemist Paracelsus discovered the Universal Constants that
+bridge elemental essences to mundane physics:
+
+ Phi (phi_fire): The Promethean constant
+ Relates Fire to Power (M*L^2*T^-3)
+ "One unit of Fire equals phi watts of transformation"
+
+ Psi (psi_water): The Heraclitan constant
+ Relates Water to Flow (L^3*T^-1)
+ "All things flow at rate psi"
+
+ Sigma (sigma_earth): The Stoic constant
+ Relates Earth to Mass (M)
+ "Substance is conserved at ratio sigma"
+
+ Omega (omega_air): The Zephyrian constant
+ Relates Air to Velocity (L*T^-1)
+ "Motion itself has measure omega"
+
+ Aether has NO constant - it transcends physical measurement!
+""")
+
+ # =========================================================================
+ # PART 3: Build the ConstantAwareBasisTransform
+ # =========================================================================
+ print("=" * 70)
+ print("PART 3: The Elemental Transform Matrix")
+ print("=" * 70)
+
+ # SI order: L, M, T, K
+ # Fire -> Power = M^1 * L^2 * T^-3 = (2, 1, -3, 0)
+ # Water -> Flow = L^3 * T^-1 = (3, 0, -1, 0)
+ # Earth -> Mass = M^1 = (0, 1, 0, 0)
+ # Air -> Velocity = L^1 * T^-1 = (1, 0, -1, 0)
+ # Aether -> Nothing = 0 = (0, 0, 0, 0)
+
+ bindings = (
+ ConstantBinding(
+ source_component=ELEMENTAL[0], # Fire
+ target_expression=Vector(SI_MECH, (Fraction(2), Fraction(1), Fraction(-3), Fraction(0))),
+ constant_symbol="phi", # Promethean constant
+ exponent=Fraction(1),
+ ),
+ ConstantBinding(
+ source_component=ELEMENTAL[1], # Water
+ target_expression=Vector(SI_MECH, (Fraction(3), Fraction(0), Fraction(-1), Fraction(0))),
+ constant_symbol="psi", # Heraclitan constant
+ exponent=Fraction(1),
+ ),
+ ConstantBinding(
+ source_component=ELEMENTAL[2], # Earth
+ target_expression=Vector(SI_MECH, (Fraction(0), Fraction(1), Fraction(0), Fraction(0))),
+ constant_symbol="sigma", # Stoic constant
+ exponent=Fraction(1),
+ ),
+ ConstantBinding(
+ source_component=ELEMENTAL[3], # Air
+ target_expression=Vector(SI_MECH, (Fraction(1), Fraction(0), Fraction(-1), Fraction(0))),
+ constant_symbol="omega", # Zephyrian constant
+ exponent=Fraction(1),
+ ),
+ # NOTE: No binding for Aether - it has no physical representation!
+ )
+
+ ELEMENTAL_TO_SI = ConstantAwareBasisTransform(
+ source=ELEMENTAL,
+ target=SI_MECH,
+ matrix=(
+ # L, M, T, K
+ (Fraction(2), Fraction(1), Fraction(-3), Fraction(0)), # Fi -> M*L^2*T^-3 (power)
+ (Fraction(3), Fraction(0), Fraction(-1), Fraction(0)), # Wa -> L^3*T^-1 (flow)
+ (Fraction(0), Fraction(1), Fraction(0), Fraction(0)), # Ea -> M (mass)
+ (Fraction(1), Fraction(0), Fraction(-1), Fraction(0)), # Ai -> L*T^-1 (velocity)
+ (Fraction(0), Fraction(0), Fraction(0), Fraction(0)), # Ae -> 0 (nothing!)
+ ),
+ bindings=bindings,
+ )
+
+ print(f"Transform: {ELEMENTAL_TO_SI}")
+ print(f"\nTransformation matrix (5 x 4, non-square!):")
+ print(" L M T K")
+ elem_names = ["Fire ", "Water", "Earth", "Air ", "Aethr"]
+ for i, name in enumerate(elem_names):
+ row = ELEMENTAL_TO_SI.matrix[i]
+ formatted = []
+ for v in row:
+ if v == 0:
+ formatted.append(" .")
+ elif v > 0:
+ formatted.append(f" +{v}")
+ else:
+ formatted.append(f" {v}")
+ print(f" {name} {''.join(formatted)}")
+
+ print("""
+Notice:
+ - Fire's row (2, 1, -3, 0) encodes POWER: L^2 * M * T^-3
+ - Water's row (3, 0, -1, 0) encodes FLOW: L^3 * T^-1
+ - Air's row (1, 0, -1, 0) encodes VELOCITY: L * T^-1
+ - Aether's row is all zeros - it has NO physical representation!
+""")
+
+ # =========================================================================
+ # PART 4: Transform Base Elemental Dimensions
+ # =========================================================================
+ print("=" * 70)
+ print("PART 4: Base Element Transformations")
+ print("=" * 70)
+
+ elements = [
+ ("Fire", (1, 0, 0, 0, 0)),
+ ("Water", (0, 1, 0, 0, 0)),
+ ("Earth", (0, 0, 1, 0, 0)),
+ ("Air", (0, 0, 0, 1, 0)),
+ ]
+
+ for name, comps in elements:
+ vec = Vector(ELEMENTAL, tuple(Fraction(c) for c in comps))
+ result = ELEMENTAL_TO_SI(vec)
+
+ # Format SI result
+ si_parts = []
+ for c in SI_MECH:
+ exp = result[c.symbol]
+ if exp == 1:
+ si_parts.append(c.symbol)
+ elif exp == -1:
+ si_parts.append(f"{c.symbol}^-1")
+ elif exp != 0:
+ si_parts.append(f"{c.symbol}^{exp}")
+
+ print(f" {name:6} -> {' * '.join(si_parts)}")
+
+ # =========================================================================
+ # PART 5: Aether - The Unphysical Dimension
+ # =========================================================================
+ print("\n" + "=" * 70)
+ print("PART 5: Aether - Beyond Physical Reality")
+ print("=" * 70)
+
+ aether = Vector(ELEMENTAL, (Fraction(0), Fraction(0), Fraction(0), Fraction(0), Fraction(1)))
+
+ print(f"""
+Aether (Ae^1) is the quintessence - it exists in the elemental realm
+but has NO representation in mundane physics.
+
+Attempting to transform pure Aether to SI:
+ {aether}
+""")
+
+ try:
+ ELEMENTAL_TO_SI(aether)
+ print(" (This should not happen!)")
+ except LossyProjection as e:
+ print(f" Raises LossyProjection!")
+ print(f" '{e.component.name}' cannot be expressed in {e.target.name}")
+
+ print("\n With allow_projection=True (the Aether is... lost):")
+ result = ELEMENTAL_TO_SI(aether, allow_projection=True)
+ print(f" {aether} -> {result}")
+ print(f" Dimensionless? {result.is_dimensionless()} (the magic dissipates)")
+
+ # =========================================================================
+ # PART 6: Compound Elemental Quantities
+ # =========================================================================
+ print("\n" + "=" * 70)
+ print("PART 6: Compound Elemental Quantities")
+ print("=" * 70)
+
+ print("""
+Alchemists work with compound essences. Let's see how they map to physics:
+""")
+
+ compounds = [
+ ("Steam (Fire + Water)", (1, 1, 0, 0, 0),
+ "Hot flowing substance"),
+ ("Earthquake (Earth * Air)", (0, 0, 1, 1, 0),
+ "Moving mass = momentum!"),
+ ("Lightning (Fire / Air)", (1, 0, 0, -1, 0),
+ "Power per velocity = force!"),
+ ("Flood (Water / Air)", (0, 1, 0, -1, 0),
+ "Flow per velocity = area!"),
+ ("Volcano (Fire * Earth / Water)", (1, -1, 1, 0, 0),
+ "Complex transformation"),
+ ("Enchantment (Fire + Aether)", (1, 0, 0, 0, 1),
+ "Magical power - LOSSY!"),
+ ]
+
+ for name, comps, desc in compounds:
+ vec = Vector(ELEMENTAL, tuple(Fraction(c) for c in comps))
+
+ # Check if it has Aether
+ has_aether = comps[4] != 0
+
+ if has_aether:
+ print(f"\n {name}")
+ print(f" Description: {desc}")
+ print(f" Elemental: {vec}")
+ try:
+ ELEMENTAL_TO_SI(vec)
+ except LossyProjection:
+ print(f" SI: LOSSY! (contains Aether)")
+ projected = ELEMENTAL_TO_SI(vec, allow_projection=True)
+ si_parts = []
+ for c in SI_MECH:
+ exp = projected[c.symbol]
+ if exp == 1:
+ si_parts.append(c.symbol)
+ elif exp != 0:
+ si_parts.append(f"{c.symbol}^{exp}")
+ print(f" Projected: {' * '.join(si_parts) if si_parts else '(dimensionless)'}")
+ else:
+ result = ELEMENTAL_TO_SI(vec)
+ si_parts = []
+ for c in SI_MECH:
+ exp = result[c.symbol]
+ if exp == 1:
+ si_parts.append(c.symbol)
+ elif exp != 0:
+ si_parts.append(f"{c.symbol}^{exp}")
+
+ print(f"\n {name}")
+ print(f" Description: {desc}")
+ print(f" Elemental: {vec}")
+ print(f" SI: {' * '.join(si_parts) if si_parts else '(dimensionless)'}")
+
+ # =========================================================================
+ # PART 7: The Inverse Transform
+ # =========================================================================
+ print("\n" + "=" * 70)
+ print("PART 7: The Inverse Transform (SI to Elemental)")
+ print("=" * 70)
+
+ SI_TO_ELEMENTAL = ELEMENTAL_TO_SI.inverse()
+
+ print(f"""
+The inverse transform lets us express physical quantities in elemental terms.
+Note: The inverse is a 4x5 matrix!
+
+{SI_TO_ELEMENTAL}
+
+Inverse bindings (constants inverted):""")
+
+ for b in SI_TO_ELEMENTAL.bindings:
+ src = b.source_component.symbol
+ # Find non-zero target components
+ tgt_parts = []
+ for c in ELEMENTAL:
+ exp = b.target_expression[c.symbol]
+ if exp == 1:
+ tgt_parts.append(c.symbol)
+ elif exp != 0:
+ tgt_parts.append(f"{c.symbol}^{exp}")
+ print(f" {src} -> {' * '.join(tgt_parts)} via {b.constant_symbol}^{b.exponent}")
+
+ # =========================================================================
+ # PART 8: Physical Quantities in Elemental Terms
+ # =========================================================================
+ print("\n" + "=" * 70)
+ print("PART 8: Expressing Physics in Elemental Terms")
+ print("=" * 70)
+
+ print("""
+How do familiar physical quantities appear in the elemental worldview?
+""")
+
+ # SI quantities: (L, M, T, K)
+ physical = [
+ ("Length (L)", (1, 0, 0, 0)),
+ ("Mass (M)", (0, 1, 0, 0)),
+ ("Time (T)", (0, 0, 1, 0)),
+ ("Velocity (L/T)", (1, 0, -1, 0)),
+ ("Force (M*L/T^2)", (1, 1, -2, 0)),
+ ("Energy (M*L^2/T^2)", (2, 1, -2, 0)),
+ ("Power (M*L^2/T^3)", (2, 1, -3, 0)),
+ ]
+
+ for name, comps in physical:
+ si_vec = Vector(SI_MECH, tuple(Fraction(c) for c in comps))
+ elem_result = SI_TO_ELEMENTAL(si_vec)
+
+ # Format elemental result
+ elem_parts = []
+ for c in ELEMENTAL:
+ exp = elem_result[c.symbol]
+ if exp == 1:
+ elem_parts.append(c.symbol)
+ elif exp == -1:
+ elem_parts.append(f"{c.symbol}^-1")
+ elif exp != 0:
+ elem_parts.append(f"{c.symbol}^{exp}")
+
+ print(f" {name:22} -> {' * '.join(elem_parts) if elem_parts else '(pure number)'}")
+
+ print("""
+Observations:
+ - Velocity IS Air (Ai) - motion itself!
+ - Power IS Fire (Fi) - transformation rate!
+ - Energy is Fire * Air^-1 - transformation without motion (stored potential)
+ - Force is Fire * Air^-2 - transformation resisting motion
+""")
+
+ # =========================================================================
+ # PART 9: The Philosopher's Stone Problem
+ # =========================================================================
+ print("=" * 70)
+ print("PART 9: The Philosopher's Stone Problem")
+ print("=" * 70)
+
+ print("""
+The legendary Philosopher's Stone is said to combine ALL five elements
+in perfect unity. What happens when we try to measure it physically?
+
+ Philosopher's Stone: Fi^1 * Wa^1 * Ea^1 * Ai^1 * Ae^1
+""")
+
+ stone = Vector(ELEMENTAL, (Fraction(1), Fraction(1), Fraction(1), Fraction(1), Fraction(1)))
+ print(f" Elemental dimension: {stone}")
+
+ try:
+ ELEMENTAL_TO_SI(stone)
+ except LossyProjection:
+ print(" Physical dimension: UNMEASURABLE (contains Aether)")
+ projected = ELEMENTAL_TO_SI(stone, allow_projection=True)
+
+ si_parts = []
+ for c in SI_MECH:
+ exp = projected[c.symbol]
+ if exp != 0:
+ si_parts.append(f"{c.symbol}^{exp}")
+
+ print(f" Projected (losing Ae): {' * '.join(si_parts)}")
+ print("""
+ The Stone's physical projection is L^6 * M^2 * T^-5:
+ = (Power) * (Flow) * (Mass) * (Velocity)
+ = (L^2*M*T^-3) * (L^3*T^-1) * (M) * (L*T^-1)
+ = L^6 * M^2 * T^-5
+
+ But without Aether, it's just... a very energetic fluid.
+ The TRUE Philosopher's Stone transcends physical measurement!
+""")
+
+ print("=" * 70)
+ print(" 'The wise alchemist knows: some things cannot be weighed.'")
+ print(" - Attributed to Hermes Trismegistus")
+ print("=" * 70 + "\n")
+
+
+if __name__ == "__main__":
+ main()
diff --git a/examples/basis/geometrized_units_demo.py b/examples/basis/geometrized_units_demo.py
new file mode 100644
index 0000000..21be350
--- /dev/null
+++ b/examples/basis/geometrized_units_demo.py
@@ -0,0 +1,261 @@
+#!/usr/bin/env python3
+"""
+ConstantAwareBasisTransform Demo: Geometrized Units (General Relativity)
+
+Demonstrates creating a custom ConstantAwareBasisTransform for geometrized
+units, commonly used in general relativity where G = c = 1.
+
+Background:
+-----------
+In geometrized units:
+ - c = 1: Length and time have the same dimension
+ - G = 1: Mass becomes a length (Schwarzschild radius ~ GM/c^2)
+
+This gives:
+ - Length [L] -> L (unchanged, our reference)
+ - Time [T] -> L (via c)
+ - Mass [M] -> L (via G/c^2)
+
+Key result: Everything is measured in lengths!
+ - 1 solar mass ~ 1.48 km (Schwarzschild radius / 2)
+ - 1 second ~ 3e8 meters (light-second)
+
+This example shows:
+ 1. Creating a custom geometrized basis
+ 2. Building a ConstantAwareBasisTransform from SI to geometrized
+ 3. Transforming physical quantities
+ 4. Computing the inverse transform
+ 5. Verifying that the Schwarzschild radius formula simplifies
+"""
+
+from fractions import Fraction
+
+from ucon.basis import (
+ Basis,
+ BasisComponent,
+ ConstantAwareBasisTransform,
+ ConstantBinding,
+ Vector,
+)
+
+
+def main():
+ print("=" * 70)
+ print("Geometrized Units (General Relativity): G = c = 1")
+ print("=" * 70)
+ print("""
+In general relativity, setting G = c = 1 simplifies equations dramatically.
+The Schwarzschild radius formula:
+
+ r_s = 2GM/c^2 (SI) --> r_s = 2M (geometrized)
+
+This example builds a ConstantAwareBasisTransform for this system.
+""")
+
+ # -------------------------------------------------------------------------
+ # Define the bases
+ # -------------------------------------------------------------------------
+ # Simplified SI with just mechanical dimensions
+ SI_MECH = Basis(
+ "SI-Mechanical",
+ [
+ BasisComponent("length", "L"),
+ BasisComponent("mass", "M"),
+ BasisComponent("time", "T"),
+ ],
+ )
+
+ # Geometrized: everything is length
+ GEOMETRIZED = Basis(
+ "Geometrized",
+ [BasisComponent("length", "L")],
+ )
+
+ print(f"Source: {SI_MECH.name} ({len(SI_MECH)} components: L, M, T)")
+ print(f"Target: {GEOMETRIZED.name} ({len(GEOMETRIZED)} component: L)")
+
+ # -------------------------------------------------------------------------
+ # Define the dimensional mappings with constant bindings
+ # -------------------------------------------------------------------------
+ # c = 1: T -> L (1 second = c meters ~ 3e8 m)
+ # G = 1: M -> L (1 kg = G/c^2 meters ~ 7.4e-28 m)
+ #
+ # More precisely:
+ # [L] = L (reference)
+ # [T] = L / c so T -> L via c^(-1)
+ # [M] = L * c^2 / G so M -> L via (G/c^2)^(+1)
+
+ bindings = (
+ # Length stays length (coefficient 1, no constant needed)
+ # We don't need a binding for clean 1:1 mappings
+
+ # Time -> Length via c
+ ConstantBinding(
+ source_component=SI_MECH[2], # time
+ target_expression=Vector(GEOMETRIZED, (Fraction(1),)), # L^1
+ constant_symbol="c",
+ exponent=Fraction(1), # multiply by c to get length
+ ),
+ # Mass -> Length via G/c^2
+ ConstantBinding(
+ source_component=SI_MECH[1], # mass
+ target_expression=Vector(GEOMETRIZED, (Fraction(1),)), # L^1
+ constant_symbol="G/c^2",
+ exponent=Fraction(1), # multiply by G/c^2 to get length
+ ),
+ )
+
+ print("\nConstant bindings:")
+ print(" L -> L (identity, no binding needed)")
+ for b in bindings:
+ print(f" {b.source_component.symbol} -> L via {b.constant_symbol}^{b.exponent}")
+
+ # -------------------------------------------------------------------------
+ # Build the transform
+ # -------------------------------------------------------------------------
+ # Matrix: each SI dimension maps to L^1 in geometrized
+ SI_TO_GEOMETRIZED = ConstantAwareBasisTransform(
+ source=SI_MECH,
+ target=GEOMETRIZED,
+ matrix=(
+ (Fraction(1),), # L -> L
+ (Fraction(1),), # M -> L (mass is a length!)
+ (Fraction(1),), # T -> L (time is a length!)
+ ),
+ bindings=bindings,
+ )
+
+ print(f"\nTransform: {SI_TO_GEOMETRIZED}")
+ print(f"Matrix: 3x1 (non-square)")
+
+ # -------------------------------------------------------------------------
+ # Transform example dimensions
+ # -------------------------------------------------------------------------
+ print("\n" + "-" * 50)
+ print("Dimensional transformations:")
+ print("-" * 50)
+
+ examples = [
+ ("Length", (1, 0, 0)),
+ ("Mass", (0, 1, 0)),
+ ("Time", (0, 0, 1)),
+ ("Velocity (L/T)", (1, 0, -1)),
+ ("Acceleration (L/T^2)", (1, 0, -2)),
+ ("Force (ML/T^2)", (1, 1, -2)),
+ ("Energy (ML^2/T^2)", (2, 1, -2)),
+ ("Schwarzschild factor (M/L)", (-1, 1, 0)),
+ ]
+
+ for name, (l, m, t) in examples:
+ vec = Vector(SI_MECH, (Fraction(l), Fraction(m), Fraction(t)))
+ result = SI_TO_GEOMETRIZED(vec)
+ exp = result["L"]
+
+ if result.is_dimensionless():
+ interpretation = "dimensionless"
+ elif exp == 1:
+ interpretation = "length"
+ elif exp == -1:
+ interpretation = "inverse length"
+ elif exp == 2:
+ interpretation = "area"
+ else:
+ interpretation = f"L^{exp}"
+
+ print(f" {name:25} -> L^{exp:+2} ({interpretation})")
+
+ # -------------------------------------------------------------------------
+ # Physical interpretation
+ # -------------------------------------------------------------------------
+ print("\n" + "-" * 50)
+ print("Physical interpretations:")
+ print("-" * 50)
+
+ print("""
+In geometrized units:
+ - Velocity is dimensionless (like c = 1 in natural units)
+ - Mass has dimension of length (Schwarzschild radius)
+ - Energy has dimension of length (via E = Mc^2 = M in geom.)
+ - The Schwarzschild radius r_s = 2GM/c^2 becomes r_s = 2M
+
+Example masses in geometrized units:
+ - Sun: M_sun ~ 1.48 km
+ - Earth: M_earth ~ 4.4 mm
+ - Human: M_human ~ 5e-26 m (unmeasurably small!)
+""")
+
+ # -------------------------------------------------------------------------
+ # Compute inverse transform
+ # -------------------------------------------------------------------------
+ print("-" * 50)
+ print("Inverse transform (GEOMETRIZED -> SI):")
+ print("-" * 50)
+
+ GEOMETRIZED_TO_SI = SI_TO_GEOMETRIZED.inverse()
+ print(f"\n{GEOMETRIZED_TO_SI}")
+ print(f"Matrix: 1x3")
+
+ print("\nInverse bindings (exponents negated):")
+ for b in GEOMETRIZED_TO_SI.bindings:
+ print(f" {b.source_component.symbol} -> {b.target_expression}")
+ print(f" via {b.constant_symbol}^{b.exponent}")
+
+ # -------------------------------------------------------------------------
+ # Verify round-trip
+ # -------------------------------------------------------------------------
+ print("\n" + "-" * 50)
+ print("Round-trip verification:")
+ print("-" * 50)
+
+ # L^1 in geometrized
+ geom_length = Vector(GEOMETRIZED, (Fraction(1),))
+ si_result = GEOMETRIZED_TO_SI(geom_length)
+ print(f"\nL^1 (geometrized) -> SI: {si_result}")
+
+ # L^2 in geometrized (area)
+ geom_area = Vector(GEOMETRIZED, (Fraction(2),))
+ si_area = GEOMETRIZED_TO_SI(geom_area)
+ print(f"L^2 (geometrized) -> SI: {si_area}")
+
+ # -------------------------------------------------------------------------
+ # The Schwarzschild radius simplification
+ # -------------------------------------------------------------------------
+ print("\n" + "-" * 50)
+ print("Schwarzschild radius formula:")
+ print("-" * 50)
+ print("""
+SI formula: r_s = 2GM/c^2
+ [r_s] = L
+
+Dimensional analysis in SI:
+ [G] = L^3 M^-1 T^-2
+ [M] = M
+ [c^2] = L^2 T^-2
+
+ [GM/c^2] = (L^3 M^-1 T^-2)(M) / (L^2 T^-2)
+ = L^3 T^-2 / L^2 T^-2
+ = L (correct!)
+
+In geometrized units (G = c = 1):
+ r_s = 2M
+
+ Since [M] = L in geometrized units,
+ [r_s] = L (still correct, but the formula is simpler!)
+""")
+
+ # Verify dimensionally
+ print("Verification via transform:")
+ # GM/c^2 in SI has dimension [L^3 M^-1 T^-2][M] / [L^2 T^-2] = L
+ # Which equals just [M] when transformed to geometrized
+ si_mass_dim = Vector(SI_MECH, (Fraction(0), Fraction(1), Fraction(0)))
+ geom_mass = SI_TO_GEOMETRIZED(si_mass_dim)
+ print(f" [M]_SI in geometrized: {geom_mass}")
+ print(f" Confirms: mass dimension = length dimension in geometrized units")
+
+ print("\n" + "=" * 70)
+ print("Demo complete!")
+ print("=" * 70 + "\n")
+
+
+if __name__ == "__main__":
+ main()
diff --git a/examples/basis/natural_units_demo.py b/examples/basis/natural_units_demo.py
new file mode 100644
index 0000000..0185d9f
--- /dev/null
+++ b/examples/basis/natural_units_demo.py
@@ -0,0 +1,284 @@
+#!/usr/bin/env python3
+"""
+ConstantAwareBasisTransform Demo: Natural Units
+
+Demonstrates how ConstantAwareBasisTransform enables inversion of non-square
+transformation matrices by tracking physical constant bindings.
+
+Background:
+-----------
+In natural units (particle physics), physical constants c, h-bar, and k_B
+are set to 1, collapsing length, time, mass, and temperature into a single
+energy dimension. This creates a non-square transform (8 SI dims -> 1 energy).
+
+A regular BasisTransform cannot invert a non-square matrix. But by recording
+which constants define each dimensional relationship (via ConstantBinding),
+ConstantAwareBasisTransform can compute the inverse.
+
+Key concepts:
+- Non-square transformation matrices (e.g., 8x1)
+- ConstantBinding records constant-dimension relationships
+- inverse() works despite non-square matrix
+- Velocity becomes dimensionless when c=1
+- Mass and energy share the same dimension (E=mc^2)
+"""
+
+from fractions import Fraction
+
+from ucon.basis import (
+ Basis,
+ BasisComponent,
+ ConstantAwareBasisTransform,
+ ConstantBinding,
+ LossyProjection,
+ Vector,
+)
+
+
+def custom_example():
+ """Demonstrate building a custom ConstantAwareBasisTransform."""
+ print("=" * 70)
+ print("PART 1: Custom ConstantAwareBasisTransform")
+ print("=" * 70)
+
+ # -------------------------------------------------------------------------
+ # Define two bases: Mechanics (3D) and Energy (1D)
+ # -------------------------------------------------------------------------
+ MECHANICS = Basis(
+ "Mechanics",
+ [
+ BasisComponent("length", "L"),
+ BasisComponent("mass", "M"),
+ BasisComponent("time", "T"),
+ ],
+ )
+
+ ENERGY = Basis(
+ "Energy",
+ [BasisComponent("energy", "E")],
+ )
+
+ print(f"\nSource basis: {MECHANICS.name} ({len(MECHANICS)} components)")
+ print(f"Target basis: {ENERGY.name} ({len(ENERGY)} components)")
+ print(f"Matrix shape: {len(MECHANICS)} x {len(ENERGY)} (non-square)")
+
+ # -------------------------------------------------------------------------
+ # Define how each mechanics dimension maps to energy
+ # -------------------------------------------------------------------------
+ # In natural units with c = h-bar = 1:
+ # Length: L = h-bar*c / E => L ~ E^(-1)
+ # Mass: M = E / c^2 => M ~ E^(+1)
+ # Time: T = h-bar / E => T ~ E^(-1)
+
+ bindings = (
+ ConstantBinding(
+ source_component=MECHANICS[0], # length
+ target_expression=Vector(ENERGY, (Fraction(-1),)), # E^(-1)
+ constant_symbol="h-bar*c",
+ exponent=Fraction(1),
+ ),
+ ConstantBinding(
+ source_component=MECHANICS[1], # mass
+ target_expression=Vector(ENERGY, (Fraction(1),)), # E^(+1)
+ constant_symbol="c",
+ exponent=Fraction(-2),
+ ),
+ ConstantBinding(
+ source_component=MECHANICS[2], # time
+ target_expression=Vector(ENERGY, (Fraction(-1),)), # E^(-1)
+ constant_symbol="h-bar",
+ exponent=Fraction(1),
+ ),
+ )
+
+ print("\nConstant bindings (how each dimension relates to energy):")
+ for b in bindings:
+ sign = "+" if b.target_expression["E"] > 0 else ""
+ print(f" {b.source_component.symbol} -> E^({sign}{b.target_expression['E']}) "
+ f"via {b.constant_symbol}^{b.exponent}")
+
+ # -------------------------------------------------------------------------
+ # Create the ConstantAwareBasisTransform
+ # -------------------------------------------------------------------------
+ transform = ConstantAwareBasisTransform(
+ source=MECHANICS,
+ target=ENERGY,
+ matrix=(
+ (Fraction(-1),), # L -> E^(-1)
+ (Fraction(1),), # M -> E^(+1)
+ (Fraction(-1),), # T -> E^(-1)
+ ),
+ bindings=bindings,
+ )
+
+ print(f"\nCreated transform: {transform}")
+
+ # -------------------------------------------------------------------------
+ # Transform dimensional vectors
+ # -------------------------------------------------------------------------
+ print("\n" + "-" * 50)
+ print("Transforming dimensional vectors:")
+ print("-" * 50)
+
+ # Length
+ length = Vector(MECHANICS, (Fraction(1), Fraction(0), Fraction(0)))
+ length_natural = transform(length)
+ print(f"\nLength (L^1):")
+ print(f" Mechanics: {length}")
+ print(f" Energy: {length_natural} (inverse energy)")
+
+ # Mass
+ mass = Vector(MECHANICS, (Fraction(0), Fraction(1), Fraction(0)))
+ mass_natural = transform(mass)
+ print(f"\nMass (M^1):")
+ print(f" Mechanics: {mass}")
+ print(f" Energy: {mass_natural} (energy!)")
+
+ # Velocity (L/T)
+ velocity = Vector(MECHANICS, (Fraction(1), Fraction(0), Fraction(-1)))
+ velocity_natural = transform(velocity)
+ print(f"\nVelocity (L^1 T^-1):")
+ print(f" Mechanics: {velocity}")
+ print(f" Energy: {velocity_natural}")
+ print(f" Dimensionless? {velocity_natural.is_dimensionless()} (c = 1!)")
+
+ # Energy (M L^2 T^-2)
+ energy = Vector(MECHANICS, (Fraction(2), Fraction(1), Fraction(-2)))
+ energy_natural = transform(energy)
+ print(f"\nEnergy (M^1 L^2 T^-2):")
+ print(f" Mechanics: {energy}")
+ print(f" Energy: {energy_natural} (pure energy, as expected)")
+
+ # -------------------------------------------------------------------------
+ # Compute the inverse transform
+ # -------------------------------------------------------------------------
+ print("\n" + "-" * 50)
+ print("Computing inverse (non-square matrix!):")
+ print("-" * 50)
+
+ inverse = transform.inverse()
+ print(f"\nInverse transform: {inverse}")
+ print(f"Matrix shape: {len(inverse.source)} x {len(inverse.target)}")
+
+ print("\nInverse bindings (exponents negated):")
+ for b in inverse.bindings:
+ print(f" {b.source_component.symbol} -> {b.target_expression.basis.name} "
+ f"via {b.constant_symbol}^{b.exponent}")
+
+ # -------------------------------------------------------------------------
+ # Verify round-trip
+ # -------------------------------------------------------------------------
+ print("\n" + "-" * 50)
+ print("Round-trip verification:")
+ print("-" * 50)
+
+ # E^1 in energy basis
+ e1 = Vector(ENERGY, (Fraction(1),))
+ back_to_mechanics = inverse(e1)
+ print(f"\nE^1 in Energy basis -> Mechanics:")
+ print(f" {e1} -> {back_to_mechanics}")
+
+
+def builtin_natural_units():
+ """Demonstrate the built-in SI_TO_NATURAL transform."""
+ print("\n")
+ print("=" * 70)
+ print("PART 2: Built-in SI -> NATURAL Transform")
+ print("=" * 70)
+
+ from ucon import NATURAL, SI, SI_TO_NATURAL, NATURAL_TO_SI, BasisVector
+
+ print(f"\nSI basis: {len(SI)} components")
+ print(f"NATURAL basis: {len(NATURAL)} component (energy only)")
+ print(f"\nSI_TO_NATURAL is a ConstantAwareBasisTransform:")
+ print(f" Type: {type(SI_TO_NATURAL).__name__}")
+ print(f" Bindings: {len(SI_TO_NATURAL.bindings)}")
+
+ print("\n" + "-" * 50)
+ print("SI dimension mappings to natural units:")
+ print("-" * 50)
+
+ # SI order: T, L, M, I, Theta, J, N, B
+ dims = [
+ ("Time", (1, 0, 0, 0, 0, 0, 0, 0)),
+ ("Length", (0, 1, 0, 0, 0, 0, 0, 0)),
+ ("Mass", (0, 0, 1, 0, 0, 0, 0, 0)),
+ ("Temperature", (0, 0, 0, 0, 1, 0, 0, 0)),
+ ("Velocity (L/T)", (-1, 1, 0, 0, 0, 0, 0, 0)),
+ ("Energy (ML^2/T^2)", (-2, 2, 1, 0, 0, 0, 0, 0)),
+ ("Action (ML^2/T)", (-1, 2, 1, 0, 0, 0, 0, 0)),
+ ("Momentum (ML/T)", (-1, 1, 1, 0, 0, 0, 0, 0)),
+ ]
+
+ for name, comps in dims:
+ vec = BasisVector(SI, tuple(Fraction(c) for c in comps))
+ result = SI_TO_NATURAL(vec)
+ dimless = " (dimensionless!)" if result.is_dimensionless() else ""
+ print(f" {name:20} -> E^{result['E']}{dimless}")
+
+ # -------------------------------------------------------------------------
+ # Show lossy projection for electromagnetic dimensions
+ # -------------------------------------------------------------------------
+ print("\n" + "-" * 50)
+ print("Electromagnetic dimensions (not representable):")
+ print("-" * 50)
+
+ si_current = BasisVector(SI, (
+ Fraction(0), Fraction(0), Fraction(0), Fraction(1),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0),
+ ))
+
+ print(f"\nSI Current (I^1): {si_current}")
+ try:
+ SI_TO_NATURAL(si_current)
+ except LossyProjection as e:
+ print(f" Raises LossyProjection: current has no representation")
+ print(f" (Pure natural units only cover mechanical + thermal dimensions)")
+
+ # With allow_projection=True
+ result = SI_TO_NATURAL(si_current, allow_projection=True)
+ print(f"\n With allow_projection=True: {result}")
+
+ # -------------------------------------------------------------------------
+ # Inverse transform
+ # -------------------------------------------------------------------------
+ print("\n" + "-" * 50)
+ print("Inverse transform (NATURAL_TO_SI):")
+ print("-" * 50)
+
+ print(f"\nNATURAL_TO_SI computed via inverse():")
+ print(f" Type: {type(NATURAL_TO_SI).__name__}")
+ print(f" Source: {NATURAL_TO_SI.source.name}")
+ print(f" Target: {NATURAL_TO_SI.target.name}")
+
+ # Transform E^1 back to SI
+ natural_energy = BasisVector(NATURAL, (Fraction(1),))
+ si_result = NATURAL_TO_SI(natural_energy)
+ print(f"\n E^1 in natural -> {si_result}")
+
+
+def main():
+ """Run all demonstrations."""
+ print("\n" + "=" * 70)
+ print("ConstantAwareBasisTransform Demonstration")
+ print("=" * 70)
+ print("""
+This demo shows how ConstantAwareBasisTransform enables inversion of
+non-square transformation matrices by recording constant bindings.
+
+In particle physics natural units (c = h-bar = k_B = 1):
+ - Velocity becomes dimensionless (c = 1)
+ - Mass and energy share dimension (E = mc^2)
+ - Length and time are inverse energy (L = h-bar*c/E, T = h-bar/E)
+""")
+
+ custom_example()
+ builtin_natural_units()
+
+ print("\n" + "=" * 70)
+ print("Demo complete!")
+ print("=" * 70 + "\n")
+
+
+if __name__ == "__main__":
+ main()
diff --git a/mkdocs.yml b/mkdocs.yml
index de885c6..dcb8f1d 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -53,6 +53,7 @@ nav:
- Config Safety: guides/dimensional-safety-config.md
- Dimensional Analysis: guides/dimensional-analysis.md
- Custom Units & Graphs: guides/custom-units-and-graphs.md
+ - Natural Units: guides/natural-units.md
- MCP Server:
- guides/mcp-server/index.md
- Custom Units: guides/mcp-server/custom-units.md
@@ -69,6 +70,7 @@ nav:
- Architecture:
- architecture/index.md
- Design Principles: architecture/design-principles.md
+ - Kind-of-Quantity Problem: architecture/kind-of-quantity.md
- ConversionGraph: architecture/conversion-graph.md
- Dual-Graph Architecture: architecture/dual-graph-architecture.md
- Suggestions & Recovery: architecture/suggestions-and-recovery.md
diff --git a/tests/ucon/mcp/test_constants.py b/tests/ucon/mcp/test_constants.py
new file mode 100644
index 0000000..cef7de3
--- /dev/null
+++ b/tests/ucon/mcp/test_constants.py
@@ -0,0 +1,267 @@
+# Copyright 2026 The Radiativity Company
+# Licensed under the Apache License, Version 2.0
+
+"""Tests for MCP constants tools (v0.9.2)."""
+
+import sys
+import pytest
+
+
+# Skip MCP tests on Python < 3.10 (FastMCP requires 3.10+)
+pytestmark = pytest.mark.skipif(
+ sys.version_info < (3, 10),
+ reason="MCP requires Python 3.10+",
+)
+
+
+class TestListConstants:
+ """Tests for list_constants tool."""
+
+ def test_returns_list(self):
+ from ucon.mcp.server import list_constants
+ result = list_constants()
+ assert isinstance(result, list)
+ assert len(result) >= 17 # Built-in constants
+
+ def test_filter_exact(self):
+ from ucon.mcp.server import list_constants
+ result = list_constants(category="exact")
+ assert isinstance(result, list)
+ for const in result:
+ assert const.is_exact
+ assert const.category == "exact"
+ assert len(result) == 7
+
+ def test_filter_derived(self):
+ from ucon.mcp.server import list_constants
+ result = list_constants(category="derived")
+ assert isinstance(result, list)
+ for const in result:
+ assert const.category == "derived"
+ assert len(result) == 3
+
+ def test_filter_measured(self):
+ from ucon.mcp.server import list_constants
+ result = list_constants(category="measured")
+ assert isinstance(result, list)
+ for const in result:
+ assert not const.is_exact
+ assert const.category == "measured"
+ assert len(result) == 7
+
+ def test_filter_session_empty_initially(self):
+ from ucon.mcp.server import list_constants, reset_session
+ reset_session()
+ result = list_constants(category="session")
+ assert isinstance(result, list)
+ assert len(result) == 0
+
+ def test_filter_invalid_category(self):
+ from ucon.mcp.server import list_constants, ConstantError
+ result = list_constants(category="invalid")
+ assert isinstance(result, ConstantError)
+ assert result.error_type == "invalid_input"
+
+ def test_includes_speed_of_light(self):
+ from ucon.mcp.server import list_constants
+ result = list_constants()
+ symbols = [c.symbol for c in result]
+ assert "c" in symbols
+
+ def test_includes_gravitational_constant(self):
+ from ucon.mcp.server import list_constants
+ result = list_constants()
+ symbols = [c.symbol for c in result]
+ assert "G" in symbols
+
+ def test_constant_info_has_all_fields(self):
+ from ucon.mcp.server import list_constants
+ result = list_constants()
+ const = result[0]
+ assert hasattr(const, 'symbol')
+ assert hasattr(const, 'name')
+ assert hasattr(const, 'value')
+ assert hasattr(const, 'unit')
+ assert hasattr(const, 'dimension')
+ assert hasattr(const, 'uncertainty')
+ assert hasattr(const, 'is_exact')
+ assert hasattr(const, 'source')
+ assert hasattr(const, 'category')
+
+
+class TestDefineConstant:
+ """Tests for define_constant tool."""
+
+ def test_define_success(self):
+ from ucon.mcp.server import define_constant, reset_session, ConstantDefinitionResult
+ reset_session()
+ result = define_constant(
+ symbol="test_vs",
+ name="speed of sound",
+ value=343,
+ unit="m/s",
+ )
+ assert isinstance(result, ConstantDefinitionResult)
+ assert result.success
+ assert result.symbol == "test_vs"
+
+ def test_duplicate_builtin_symbol_fails(self):
+ from ucon.mcp.server import define_constant, ConstantError
+ result = define_constant(
+ symbol="c",
+ name="my speed",
+ value=300000000,
+ unit="m/s",
+ )
+ assert isinstance(result, ConstantError)
+ assert result.error_type == "duplicate_symbol"
+
+ def test_duplicate_session_symbol_fails(self):
+ from ucon.mcp.server import define_constant, reset_session, ConstantError
+ reset_session()
+ # Define first time
+ define_constant(
+ symbol="test_dup",
+ name="first",
+ value=1,
+ unit="m",
+ )
+ # Try to define again
+ result = define_constant(
+ symbol="test_dup",
+ name="second",
+ value=2,
+ unit="m",
+ )
+ assert isinstance(result, ConstantError)
+ assert result.error_type == "duplicate_symbol"
+
+ def test_invalid_unit_fails(self):
+ from ucon.mcp.server import define_constant, reset_session, ConstantError
+ reset_session()
+ result = define_constant(
+ symbol="test_Y",
+ name="Y",
+ value=1,
+ unit="invalid_unit_xyz",
+ )
+ assert isinstance(result, ConstantError)
+ assert result.error_type == "invalid_unit"
+
+ def test_nan_value_fails(self):
+ import math
+ from ucon.mcp.server import define_constant, reset_session, ConstantError
+ reset_session()
+ result = define_constant(
+ symbol="test_nan",
+ name="NaN constant",
+ value=math.nan,
+ unit="m",
+ )
+ assert isinstance(result, ConstantError)
+ assert result.error_type == "invalid_value"
+
+ def test_inf_value_fails(self):
+ import math
+ from ucon.mcp.server import define_constant, reset_session, ConstantError
+ reset_session()
+ result = define_constant(
+ symbol="test_inf",
+ name="Inf constant",
+ value=math.inf,
+ unit="m",
+ )
+ assert isinstance(result, ConstantError)
+ assert result.error_type == "invalid_value"
+
+ def test_negative_uncertainty_fails(self):
+ from ucon.mcp.server import define_constant, reset_session, ConstantError
+ reset_session()
+ result = define_constant(
+ symbol="test_neg_unc",
+ name="negative uncertainty",
+ value=1.0,
+ unit="m",
+ uncertainty=-0.1,
+ )
+ assert isinstance(result, ConstantError)
+ assert result.error_type == "invalid_value"
+
+ def test_with_uncertainty(self):
+ from ucon.mcp.server import define_constant, reset_session, ConstantDefinitionResult
+ reset_session()
+ result = define_constant(
+ symbol="test_unc",
+ name="with uncertainty",
+ value=1.0,
+ unit="m",
+ uncertainty=0.01,
+ )
+ assert isinstance(result, ConstantDefinitionResult)
+ assert result.success
+ assert result.uncertainty == 0.01
+
+ def test_defined_constant_appears_in_session_list(self):
+ from ucon.mcp.server import define_constant, list_constants, reset_session
+ reset_session()
+ define_constant(
+ symbol="test_sess",
+ name="session constant",
+ value=42,
+ unit="kg",
+ )
+ result = list_constants(category="session")
+ assert len(result) == 1
+ assert result[0].symbol == "test_sess"
+ assert result[0].category == "session"
+
+
+class TestResetSession:
+ """Tests for reset_session clearing constants."""
+
+ def test_reset_clears_session_constants(self):
+ from ucon.mcp.server import define_constant, list_constants, reset_session
+ reset_session()
+ # Define a constant
+ define_constant(
+ symbol="test_clear",
+ name="to be cleared",
+ value=1,
+ unit="m",
+ )
+ # Verify it exists
+ result = list_constants(category="session")
+ assert len(result) == 1
+
+ # Reset
+ reset_session()
+
+ # Verify it's gone
+ result = list_constants(category="session")
+ assert len(result) == 0
+
+
+class TestConstantsCategoryCounts:
+ """Tests for correct counts of constants by category."""
+
+ def test_total_builtin_constants(self):
+ from ucon.mcp.server import list_constants
+ exact = list_constants(category="exact")
+ derived = list_constants(category="derived")
+ measured = list_constants(category="measured")
+ total = len(exact) + len(derived) + len(measured)
+ assert total == 17
+
+ def test_exact_constants_are_exact(self):
+ from ucon.mcp.server import list_constants
+ exact = list_constants(category="exact")
+ for const in exact:
+ assert const.uncertainty is None
+ assert const.is_exact is True
+
+ def test_measured_constants_have_uncertainty(self):
+ from ucon.mcp.server import list_constants
+ measured = list_constants(category="measured")
+ for const in measured:
+ assert const.uncertainty is not None
+ assert const.is_exact is False
diff --git a/tests/ucon/test_natural_units.py b/tests/ucon/test_natural_units.py
new file mode 100644
index 0000000..b62634e
--- /dev/null
+++ b/tests/ucon/test_natural_units.py
@@ -0,0 +1,449 @@
+# Copyright 2026 The Radiativity Company
+# Licensed under the Apache License, Version 2.0
+
+"""Tests for natural units support (v0.9.3).
+
+Natural units set c = h_bar = k_B = 1, collapsing length, time, mass, and
+temperature into expressions of a single energy dimension.
+"""
+
+from fractions import Fraction
+
+import pytest
+
+from ucon.basis import (
+ BasisComponent,
+ ConstantAwareBasisTransform,
+ ConstantBinding,
+ LossyProjection,
+ Vector,
+)
+from ucon.bases import (
+ NATURAL,
+ NATURAL_TO_SI,
+ SI,
+ SI_TO_NATURAL,
+)
+
+
+# -----------------------------------------------------------------------------
+# NATURAL Basis Structure Tests
+# -----------------------------------------------------------------------------
+
+
+class TestNaturalBasis:
+ """Tests for the NATURAL basis structure."""
+
+ def test_has_single_component(self):
+ """GIVEN NATURAL basis, THEN it has exactly 1 component (energy)."""
+ assert len(NATURAL) == 1
+
+ def test_component_is_energy(self):
+ """GIVEN NATURAL basis, THEN its component is named 'energy'."""
+ assert NATURAL[0].name == "energy"
+ assert NATURAL[0].symbol == "E"
+
+ def test_energy_in_basis(self):
+ """GIVEN NATURAL basis, THEN 'energy' and 'E' are in the basis."""
+ assert "energy" in NATURAL
+ assert "E" in NATURAL
+
+ def test_basis_name(self):
+ """GIVEN NATURAL basis, THEN its name is 'natural'."""
+ assert NATURAL.name == "natural"
+
+
+# -----------------------------------------------------------------------------
+# ConstantBinding Tests
+# -----------------------------------------------------------------------------
+
+
+class TestConstantBinding:
+ """Tests for ConstantBinding dataclass."""
+
+ def test_frozen_immutable(self):
+ """GIVEN a ConstantBinding, THEN it cannot be mutated."""
+ binding = ConstantBinding(
+ source_component=BasisComponent("length", "L"),
+ target_expression=Vector(NATURAL, (Fraction(-1),)),
+ constant_symbol="c",
+ exponent=Fraction(1),
+ )
+ with pytest.raises(AttributeError):
+ binding.constant_symbol = "h"
+
+ def test_hashable(self):
+ """GIVEN a ConstantBinding, THEN it is hashable for dict/set usage."""
+ binding = ConstantBinding(
+ source_component=BasisComponent("length", "L"),
+ target_expression=Vector(NATURAL, (Fraction(-1),)),
+ constant_symbol="c",
+ exponent=Fraction(1),
+ )
+ d = {binding: "value"}
+ assert d[binding] == "value"
+
+ def test_default_exponent(self):
+ """GIVEN a ConstantBinding without exponent, THEN default is 1."""
+ binding = ConstantBinding(
+ source_component=BasisComponent("length", "L"),
+ target_expression=Vector(NATURAL, (Fraction(-1),)),
+ constant_symbol="c",
+ )
+ assert binding.exponent == Fraction(1)
+
+ def test_fractional_exponent(self):
+ """GIVEN a ConstantBinding with fractional exponent, THEN it is stored."""
+ binding = ConstantBinding(
+ source_component=BasisComponent("current", "I"),
+ target_expression=Vector(NATURAL, (Fraction(1),)),
+ constant_symbol="e",
+ exponent=Fraction(1, 2),
+ )
+ assert binding.exponent == Fraction(1, 2)
+
+
+# -----------------------------------------------------------------------------
+# SI_TO_NATURAL Transform Tests
+# -----------------------------------------------------------------------------
+
+
+class TestSIToNatural:
+ """Tests for SI_TO_NATURAL transform."""
+
+ def test_is_constant_aware_transform(self):
+ """GIVEN SI_TO_NATURAL, THEN it is a ConstantAwareBasisTransform."""
+ assert isinstance(SI_TO_NATURAL, ConstantAwareBasisTransform)
+
+ def test_source_is_si(self):
+ """GIVEN SI_TO_NATURAL, THEN source is SI basis."""
+ assert SI_TO_NATURAL.source == SI
+
+ def test_target_is_natural(self):
+ """GIVEN SI_TO_NATURAL, THEN target is NATURAL basis."""
+ assert SI_TO_NATURAL.target == NATURAL
+
+ def test_has_bindings(self):
+ """GIVEN SI_TO_NATURAL, THEN it has constant bindings."""
+ assert len(SI_TO_NATURAL.bindings) == 4 # L, T, M, Θ
+
+ def test_length_to_inverse_energy(self):
+ """GIVEN length dimension, THEN it transforms to E⁻¹."""
+ # SI: T=0, L=1, M=0, I=0, Θ=0, J=0, N=0, B=0
+ si_length = Vector(
+ SI,
+ (Fraction(0), Fraction(1), Fraction(0), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0)),
+ )
+ natural_result = SI_TO_NATURAL(si_length)
+
+ assert natural_result.basis == NATURAL
+ assert natural_result["E"] == Fraction(-1)
+
+ def test_time_to_inverse_energy(self):
+ """GIVEN time dimension, THEN it transforms to E⁻¹."""
+ # SI: T=1, L=0, M=0, I=0, Θ=0, J=0, N=0, B=0
+ si_time = Vector(
+ SI,
+ (Fraction(1), Fraction(0), Fraction(0), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0)),
+ )
+ natural_result = SI_TO_NATURAL(si_time)
+
+ assert natural_result["E"] == Fraction(-1)
+
+ def test_mass_to_energy(self):
+ """GIVEN mass dimension, THEN it transforms to E."""
+ # SI: T=0, L=0, M=1, I=0, Θ=0, J=0, N=0, B=0
+ si_mass = Vector(
+ SI,
+ (Fraction(0), Fraction(0), Fraction(1), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0)),
+ )
+ natural_result = SI_TO_NATURAL(si_mass)
+
+ assert natural_result["E"] == Fraction(1)
+
+ def test_temperature_to_energy(self):
+ """GIVEN temperature dimension, THEN it transforms to E."""
+ # SI: T=0, L=0, M=0, I=0, Θ=1, J=0, N=0, B=0
+ si_temp = Vector(
+ SI,
+ (Fraction(0), Fraction(0), Fraction(0), Fraction(0),
+ Fraction(1), Fraction(0), Fraction(0), Fraction(0)),
+ )
+ natural_result = SI_TO_NATURAL(si_temp)
+
+ assert natural_result["E"] == Fraction(1)
+
+ def test_velocity_is_dimensionless(self):
+ """GIVEN velocity (L/T), THEN it transforms to E⁰ (dimensionless).
+
+ This is the key consequence of c = 1: velocity is dimensionless.
+ """
+ # Velocity: L¹T⁻¹
+ # SI: T=-1, L=1, M=0, I=0, Θ=0, J=0, N=0, B=0
+ si_velocity = Vector(
+ SI,
+ (Fraction(-1), Fraction(1), Fraction(0), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0)),
+ )
+ natural_result = SI_TO_NATURAL(si_velocity)
+
+ # L → E⁻¹, T⁻¹ → E¹, so L¹T⁻¹ → E⁻¹·E¹ = E⁰
+ assert natural_result["E"] == Fraction(0)
+ assert natural_result.is_dimensionless()
+
+ def test_energy_dimension_preserved(self):
+ """GIVEN energy dimension (ML²T⁻²), THEN it transforms to E¹."""
+ # Energy: M¹L²T⁻²
+ # SI: T=-2, L=2, M=1, I=0, Θ=0, J=0, N=0, B=0
+ si_energy = Vector(
+ SI,
+ (Fraction(-2), Fraction(2), Fraction(1), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0)),
+ )
+ natural_result = SI_TO_NATURAL(si_energy)
+
+ # M → E, L² → E⁻², T⁻² → E², so M¹L²T⁻² → E¹·E⁻²·E² = E¹
+ assert natural_result["E"] == Fraction(1)
+
+ def test_action_is_dimensionless(self):
+ """GIVEN action (ML²T⁻¹), THEN it transforms to E⁰ (dimensionless).
+
+ This is the consequence of ℏ = 1: action is dimensionless.
+ """
+ # Action: M¹L²T⁻¹ (same as ℏ)
+ # SI: T=-1, L=2, M=1, I=0, Θ=0, J=0, N=0, B=0
+ si_action = Vector(
+ SI,
+ (Fraction(-1), Fraction(2), Fraction(1), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0)),
+ )
+ natural_result = SI_TO_NATURAL(si_action)
+
+ # M → E, L² → E⁻², T⁻¹ → E¹, so M¹L²T⁻¹ → E¹·E⁻²·E¹ = E⁰
+ assert natural_result["E"] == Fraction(0)
+ assert natural_result.is_dimensionless()
+
+ def test_current_raises_lossy_projection(self):
+ """GIVEN current dimension (I), THEN LossyProjection is raised."""
+ # SI: T=0, L=0, M=0, I=1, Θ=0, J=0, N=0, B=0
+ si_current = Vector(
+ SI,
+ (Fraction(0), Fraction(0), Fraction(0), Fraction(1),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0)),
+ )
+
+ with pytest.raises(LossyProjection) as exc_info:
+ SI_TO_NATURAL(si_current)
+
+ assert "current" in str(exc_info.value)
+
+ def test_current_with_allow_projection(self):
+ """GIVEN current with allow_projection=True, THEN projected to zero."""
+ si_current = Vector(
+ SI,
+ (Fraction(0), Fraction(0), Fraction(0), Fraction(1),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0)),
+ )
+
+ result = SI_TO_NATURAL(si_current, allow_projection=True)
+ assert result.is_dimensionless()
+
+ def test_luminous_intensity_raises_lossy_projection(self):
+ """GIVEN luminous_intensity (J), THEN LossyProjection is raised."""
+ # SI: T=0, L=0, M=0, I=0, Θ=0, J=1, N=0, B=0
+ si_luminosity = Vector(
+ SI,
+ (Fraction(0), Fraction(0), Fraction(0), Fraction(0),
+ Fraction(0), Fraction(1), Fraction(0), Fraction(0)),
+ )
+
+ with pytest.raises(LossyProjection):
+ SI_TO_NATURAL(si_luminosity)
+
+
+# -----------------------------------------------------------------------------
+# NATURAL_TO_SI (Inverse Transform) Tests
+# -----------------------------------------------------------------------------
+
+
+class TestNaturalToSI:
+ """Tests for NATURAL_TO_SI inverse transform."""
+
+ def test_is_constant_aware_transform(self):
+ """GIVEN NATURAL_TO_SI, THEN it is a ConstantAwareBasisTransform."""
+ assert isinstance(NATURAL_TO_SI, ConstantAwareBasisTransform)
+
+ def test_source_is_natural(self):
+ """GIVEN NATURAL_TO_SI, THEN source is NATURAL basis."""
+ assert NATURAL_TO_SI.source == NATURAL
+
+ def test_target_is_si(self):
+ """GIVEN NATURAL_TO_SI, THEN target is SI basis."""
+ assert NATURAL_TO_SI.target == SI
+
+ def test_bindings_have_negated_exponents(self):
+ """GIVEN NATURAL_TO_SI, THEN binding exponents are negated from forward."""
+ # The inverse should have bindings with negated exponents
+ for inv_binding in NATURAL_TO_SI.bindings:
+ # Find corresponding forward binding
+ for fwd_binding in SI_TO_NATURAL.bindings:
+ if inv_binding.constant_symbol == fwd_binding.constant_symbol:
+ assert inv_binding.exponent == -fwd_binding.exponent
+ break
+
+ def test_energy_transforms_to_si(self):
+ """GIVEN energy in natural units, THEN it can transform to SI."""
+ natural_energy = Vector(NATURAL, (Fraction(1),))
+ si_result = NATURAL_TO_SI(natural_energy)
+
+ assert si_result.basis == SI
+ # The result should have the SI energy dimension signature
+ # Note: The exact mapping depends on which binding is primary
+
+
+# -----------------------------------------------------------------------------
+# Round-Trip Tests
+# -----------------------------------------------------------------------------
+
+
+class TestRoundTrip:
+ """Tests for SI → NATURAL → SI round-trip consistency."""
+
+ def test_velocity_round_trip(self):
+ """GIVEN velocity, WHEN round-tripped, THEN dimensionless preserved."""
+ # Velocity L¹T⁻¹ → dimensionless (c=1)
+ si_velocity = Vector(
+ SI,
+ (Fraction(-1), Fraction(1), Fraction(0), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0)),
+ )
+
+ natural = SI_TO_NATURAL(si_velocity)
+ assert natural.is_dimensionless()
+
+ # Going back: dimensionless stays dimensionless
+ si_back = NATURAL_TO_SI(natural)
+ assert si_back.is_dimensionless()
+
+ def test_energy_round_trip(self):
+ """GIVEN energy dimension, WHEN round-tripped via natural units, THEN consistent."""
+ # Energy M¹L²T⁻² → E¹
+ si_energy = Vector(
+ SI,
+ (Fraction(-2), Fraction(2), Fraction(1), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0)),
+ )
+
+ natural = SI_TO_NATURAL(si_energy)
+ assert natural["E"] == Fraction(1)
+
+ # Going back gives the primary SI representation of energy
+ si_back = NATURAL_TO_SI(natural)
+ assert si_back.basis == SI
+
+
+# -----------------------------------------------------------------------------
+# Particle Physics Dimensions Tests
+# -----------------------------------------------------------------------------
+
+
+class TestParticlePhysicsDimensions:
+ """Tests for particle physics dimensional analysis."""
+
+ def test_cross_section_dimension(self):
+ """GIVEN cross-section (L²), THEN it transforms to E⁻²."""
+ # Cross-section: L²
+ si_cross_section = Vector(
+ SI,
+ (Fraction(0), Fraction(2), Fraction(0), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0)),
+ )
+ natural_result = SI_TO_NATURAL(si_cross_section)
+
+ # L² → (E⁻¹)² = E⁻²
+ assert natural_result["E"] == Fraction(-2)
+
+ def test_decay_width_dimension(self):
+ """GIVEN decay width (T⁻¹), THEN it transforms to E."""
+ # Decay width: T⁻¹ (same dimension as frequency)
+ si_decay_width = Vector(
+ SI,
+ (Fraction(-1), Fraction(0), Fraction(0), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0)),
+ )
+ natural_result = SI_TO_NATURAL(si_decay_width)
+
+ # T⁻¹ → (E⁻¹)⁻¹ = E
+ assert natural_result["E"] == Fraction(1)
+
+ def test_momentum_dimension(self):
+ """GIVEN momentum (MLT⁻¹), THEN it transforms to E."""
+ # Momentum: M¹L¹T⁻¹
+ si_momentum = Vector(
+ SI,
+ (Fraction(-1), Fraction(1), Fraction(1), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0)),
+ )
+ natural_result = SI_TO_NATURAL(si_momentum)
+
+ # M → E, L → E⁻¹, T⁻¹ → E, so M¹L¹T⁻¹ → E¹·E⁻¹·E¹ = E¹
+ assert natural_result["E"] == Fraction(1)
+
+ def test_wavelength_dimension(self):
+ """GIVEN wavelength (L), THEN it transforms to E⁻¹."""
+ # Wavelength is just length
+ si_wavelength = Vector(
+ SI,
+ (Fraction(0), Fraction(1), Fraction(0), Fraction(0),
+ Fraction(0), Fraction(0), Fraction(0), Fraction(0)),
+ )
+ natural_result = SI_TO_NATURAL(si_wavelength)
+
+ assert natural_result["E"] == Fraction(-1)
+
+
+# -----------------------------------------------------------------------------
+# ConstantAwareBasisTransform General Tests
+# -----------------------------------------------------------------------------
+
+
+class TestConstantAwareBasisTransform:
+ """General tests for ConstantAwareBasisTransform."""
+
+ def test_wrong_basis_raises(self):
+ """GIVEN a vector in wrong basis, THEN ValueError raised."""
+ natural_vec = Vector(NATURAL, (Fraction(1),))
+
+ with pytest.raises(ValueError, match="expects basis 'SI'"):
+ SI_TO_NATURAL(natural_vec)
+
+ def test_matrix_dimension_validation(self):
+ """GIVEN wrong matrix dimensions, THEN ValueError raised."""
+ with pytest.raises(ValueError, match="Matrix has 2 rows"):
+ ConstantAwareBasisTransform(
+ source=SI, # 8 components
+ target=NATURAL, # 1 component
+ matrix=(
+ (Fraction(1),),
+ (Fraction(1),),
+ # Missing 6 rows
+ ),
+ )
+
+ def test_as_basis_transform(self):
+ """GIVEN ConstantAwareBasisTransform, THEN as_basis_transform works."""
+ from ucon.basis import BasisTransform
+
+ plain = SI_TO_NATURAL.as_basis_transform()
+ assert isinstance(plain, BasisTransform)
+ assert plain.source == SI
+ assert plain.target == NATURAL
+ assert plain.matrix == SI_TO_NATURAL.matrix
+
+ def test_repr(self):
+ """GIVEN ConstantAwareBasisTransform, THEN repr is informative."""
+ assert "SI" in repr(SI_TO_NATURAL)
+ assert "natural" in repr(SI_TO_NATURAL)
diff --git a/ucon/__init__.py b/ucon/__init__.py
index 109c584..f708d80 100644
--- a/ucon/__init__.py
+++ b/ucon/__init__.py
@@ -44,6 +44,8 @@
BasisComponent,
BasisGraph,
BasisTransform,
+ ConstantAwareBasisTransform,
+ ConstantBinding,
LossyProjection,
NoTransformPath,
Vector as BasisVector,
@@ -58,9 +60,12 @@
CGS,
CGS_ESU,
CGS_TO_SI,
+ NATURAL,
+ NATURAL_TO_SI,
SI,
SI_TO_CGS,
SI_TO_CGS_ESU,
+ SI_TO_NATURAL,
)
from ucon.core import (
DimConstraint,
@@ -93,6 +98,8 @@
'BasisGraph',
'BasisTransform',
'BasisVector',
+ 'ConstantAwareBasisTransform',
+ 'ConstantBinding',
'LossyProjection',
'NoTransformPath',
# Basis context scoping
@@ -105,11 +112,14 @@
# Standard bases
'CGS',
'CGS_ESU',
+ 'NATURAL',
'SI',
# Standard transforms
'CGS_TO_SI',
+ 'NATURAL_TO_SI',
'SI_TO_CGS',
'SI_TO_CGS_ESU',
+ 'SI_TO_NATURAL',
# Core types
'Constant',
'DimConstraint',
diff --git a/ucon/bases.py b/ucon/bases.py
index 2ee6f15..1cfef33 100644
--- a/ucon/bases.py
+++ b/ucon/bases.py
@@ -22,7 +22,14 @@
from fractions import Fraction
-from ucon.basis import Basis, BasisComponent, BasisTransform
+from ucon.basis import (
+ Basis,
+ BasisComponent,
+ BasisTransform,
+ ConstantAwareBasisTransform,
+ ConstantBinding,
+ Vector,
+)
# -----------------------------------------------------------------------------
@@ -153,3 +160,112 @@
# is not a clean projection — current (I) maps to a complex derived dimension
# L^(3/2) M^(1/2) T^(-2), not a simple 1:1 mapping. Users needing CGS-ESU -> SI
# conversion should construct the transform manually based on their use case.
+
+
+# -----------------------------------------------------------------------------
+# Natural Units Basis
+# -----------------------------------------------------------------------------
+
+NATURAL = Basis(
+ "natural",
+ [BasisComponent("energy", "E")],
+)
+"""Natural units basis with single energy dimension.
+
+In natural units (particle physics conventions), physical constants c, ℏ, and k_B
+are set to 1. This collapses length, time, mass, and temperature into expressions
+of energy:
+
+- c = 1 → Length and time have same dimension
+- ℏ = 1 → Energy × time is dimensionless
+- k_B = 1 → Temperature is energy
+
+As a result:
+- Mass [M] → E (via E = mc²)
+- Length [L] → E⁻¹ (via ℏc/E)
+- Time [T] → E⁻¹ (via ℏ/E)
+- Temperature [Θ] → E (via k_B T)
+- Velocity is dimensionless (c = 1)
+
+Electromagnetic quantities (current, etc.) are not representable in pure
+natural units and will raise LossyProjection.
+"""
+
+# Bindings for SI → NATURAL (c = ℏ = k_B = 1)
+_NATURAL_BINDINGS = (
+ # Length: L = ℏc/E → L ~ E⁻¹ via ℏc
+ ConstantBinding(
+ source_component=SI[1], # length (index 1 in SI)
+ target_expression=Vector(NATURAL, (Fraction(-1),)),
+ constant_symbol="ℏc",
+ exponent=Fraction(1),
+ ),
+ # Time: T = ℏ/E → T ~ E⁻¹ via ℏ
+ ConstantBinding(
+ source_component=SI[0], # time (index 0 in SI)
+ target_expression=Vector(NATURAL, (Fraction(-1),)),
+ constant_symbol="ℏ",
+ exponent=Fraction(1),
+ ),
+ # Mass: M = E/c² → M ~ E via c⁻²
+ ConstantBinding(
+ source_component=SI[2], # mass (index 2 in SI)
+ target_expression=Vector(NATURAL, (Fraction(1),)),
+ constant_symbol="c",
+ exponent=Fraction(-2),
+ ),
+ # Temperature: Θ = E/k_B → Θ ~ E via k_B⁻¹
+ ConstantBinding(
+ source_component=SI[4], # temperature (index 4 in SI)
+ target_expression=Vector(NATURAL, (Fraction(1),)),
+ constant_symbol="k_B",
+ exponent=Fraction(-1),
+ ),
+)
+
+SI_TO_NATURAL = ConstantAwareBasisTransform(
+ source=SI,
+ target=NATURAL,
+ matrix=(
+ # SI order: T, L, M, I, Θ, J, N, B
+ # NATURAL order: E
+ (Fraction(-1),), # T → E⁻¹
+ (Fraction(-1),), # L → E⁻¹
+ (Fraction(1),), # M → E
+ (Fraction(0),), # I → 0 (not representable)
+ (Fraction(1),), # Θ → E
+ (Fraction(0),), # J → 0 (not representable)
+ (Fraction(0),), # N → 0 (not representable)
+ (Fraction(0),), # B → 0 (not representable)
+ ),
+ bindings=_NATURAL_BINDINGS,
+)
+"""Transform from SI to natural units.
+
+Maps SI dimensions to the single energy dimension in natural units:
+- Time (T) → E⁻¹ (via ℏ)
+- Length (L) → E⁻¹ (via ℏc)
+- Mass (M) → E (via c²)
+- Temperature (Θ) → E (via k_B)
+
+Current (I), luminous_intensity (J), amount_of_substance (N), and
+information (B) are not representable in natural units and will raise
+LossyProjection if non-zero (unless allow_projection=True).
+
+Key consequences:
+- Velocity (L/T) → E⁰ (dimensionless, since c = 1)
+- Energy (ML²T⁻²) → E (as expected)
+- Action (ML²T⁻¹) → E⁰ (dimensionless, since ℏ = 1)
+"""
+
+NATURAL_TO_SI = SI_TO_NATURAL.inverse()
+"""Transform from natural units back to SI.
+
+This is the inverse of SI_TO_NATURAL, computed using the constant bindings.
+Allows converting natural unit dimensions back to their SI representation.
+
+Note: Information about which specific combination of L, T, M, Θ a given
+E dimension originated from is tracked via the constant bindings. However,
+the numeric conversion factors require the actual constant values from
+ucon.constants.
+"""
diff --git a/ucon/basis.py b/ucon/basis.py
index 36d8118..9d3be7c 100644
--- a/ucon/basis.py
+++ b/ucon/basis.py
@@ -11,7 +11,7 @@
from __future__ import annotations
-from dataclasses import dataclass
+from dataclasses import dataclass, field
from fractions import Fraction
from typing import TYPE_CHECKING, Iterator, Sequence
@@ -621,6 +621,285 @@ def is_identity(self) -> bool:
return True
+# -----------------------------------------------------------------------------
+# ConstantBinding
+# -----------------------------------------------------------------------------
+
+
+@dataclass(frozen=True)
+class ConstantBinding:
+ """Binds a source dimension to a target expression via a physical constant.
+
+ Records that `source_component` in the source basis becomes
+ `target_expression` in the target basis, with the relationship
+ defined by a physical constant raised to `exponent`.
+
+ Parameters
+ ----------
+ source_component : BasisComponent
+ The fundamental component being transformed.
+ target_expression : Vector
+ How it expresses in the target basis.
+ constant_symbol : str
+ Symbol of the constant defining this relationship (e.g., "c", "ℏ").
+ We use a string rather than a Constant object to avoid circular imports
+ between ucon.basis and ucon.constants.
+ exponent : Fraction
+ Power of the constant in the relationship (usually ±1/2 or ±1).
+
+ Examples
+ --------
+ >>> # In natural units: length = ℏc/E, so L → E⁻¹ via ℏc
+ >>> from fractions import Fraction
+ >>> binding = ConstantBinding(
+ ... source_component=BasisComponent("length", "L"),
+ ... target_expression=Vector(NATURAL, (Fraction(-1),)),
+ ... constant_symbol="ℏc",
+ ... exponent=Fraction(1),
+ ... )
+ """
+
+ source_component: BasisComponent
+ target_expression: "Vector"
+ constant_symbol: str
+ exponent: Fraction = field(default_factory=lambda: Fraction(1))
+
+
+# -----------------------------------------------------------------------------
+# ConstantAwareBasisTransform
+# -----------------------------------------------------------------------------
+
+
+@dataclass(frozen=True)
+class ConstantAwareBasisTransform:
+ """A basis transform with constants that enable inversion of non-square matrices.
+
+ Extends BasisTransform with explicit constant bindings that record which
+ constants define each non-trivial mapping. This enables `inverse()` to work
+ on non-square transforms by providing the information needed to reverse
+ derived mappings.
+
+ Parameters
+ ----------
+ source : Basis
+ Source basis.
+ target : Basis
+ Target basis.
+ matrix : tuple[tuple[Fraction, ...], ...]
+ Dimensional transformation matrix. Shape is (len(source), len(target)).
+ bindings : tuple[ConstantBinding, ...]
+ Constants that define derived relationships.
+
+ Examples
+ --------
+ >>> # SI (8 dimensions) → NATURAL (1 dimension) transform
+ >>> SI_TO_NATURAL = ConstantAwareBasisTransform(
+ ... source=SI,
+ ... target=NATURAL,
+ ... matrix=(...), # 8×1 matrix
+ ... bindings=(...), # Bindings for L, T, M, Θ → E
+ ... )
+ >>> # This works because bindings record how to reverse!
+ >>> NATURAL_TO_SI = SI_TO_NATURAL.inverse()
+ """
+
+ source: Basis
+ target: Basis
+ matrix: tuple[tuple[Fraction, ...], ...]
+ bindings: tuple[ConstantBinding, ...] = ()
+
+ def __post_init__(self) -> None:
+ if len(self.matrix) != len(self.source):
+ raise ValueError(
+ f"Matrix has {len(self.matrix)} rows but source basis "
+ f"'{self.source.name}' has {len(self.source)} components"
+ )
+ for i, row in enumerate(self.matrix):
+ if len(row) != len(self.target):
+ raise ValueError(
+ f"Matrix row {i} has {len(row)} columns but target basis "
+ f"'{self.target.name}' has {len(self.target)} components"
+ )
+
+ def __repr__(self) -> str:
+ return f"ConstantAwareBasisTransform({self.source.name} -> {self.target.name})"
+
+ def __call__(self, vector: Vector, *, allow_projection: bool = False) -> Vector:
+ """Transform a vector from source basis to target basis.
+
+ Args:
+ vector: The vector to transform (must be in source basis).
+ allow_projection: If False (default), raise LossyProjection when
+ a non-zero component would be discarded. If True, silently
+ project (drop) unrepresentable components.
+
+ Returns:
+ A new Vector in the target basis.
+
+ Raises:
+ ValueError: If vector is not in the source basis.
+ LossyProjection: If allow_projection=False and a non-zero
+ component maps entirely to zeros in the target basis.
+ """
+ if vector.basis != self.source:
+ raise ValueError(
+ f"Vector is in basis '{vector.basis.name}' but transform "
+ f"expects basis '{self.source.name}'"
+ )
+
+ # Check for lossy projection (unless explicitly allowed)
+ if not allow_projection:
+ for i, src_exp in enumerate(vector.components):
+ if src_exp != 0:
+ # Does this component have ANY representation in target?
+ if all(self.matrix[i][j] == 0 for j in range(len(self.target))):
+ raise LossyProjection(
+ self.source[i],
+ self.source,
+ self.target,
+ )
+
+ # Matrix-vector multiplication
+ result = [Fraction(0)] * len(self.target)
+ for i, src_exp in enumerate(vector.components):
+ for j, transform_coeff in enumerate(self.matrix[i]):
+ result[j] += src_exp * transform_coeff
+
+ return Vector(self.target, tuple(result))
+
+ def inverse(self) -> "ConstantAwareBasisTransform":
+ """Compute the inverse transform using constant bindings.
+
+ For each binding (src_component → target_expression via constant):
+ - The inverse maps target_expression → src_component
+ - The constant exponent is negated (constant^(-exponent))
+
+ Returns
+ -------
+ ConstantAwareBasisTransform
+ The inverse transform from target to source.
+
+ Raises
+ ------
+ ValueError
+ If a source component has no binding and cannot be inverted.
+
+ Notes
+ -----
+ Components that map to zero (truly dropped) cannot be recovered
+ and will map to zero in the inverse. Only components with bindings
+ are invertible for non-square matrices.
+ """
+ inv_matrix = self._compute_inverse_matrix()
+ inv_bindings = self._invert_bindings()
+
+ return ConstantAwareBasisTransform(
+ source=self.target,
+ target=self.source,
+ matrix=inv_matrix,
+ bindings=inv_bindings,
+ )
+
+ def _compute_inverse_matrix(self) -> tuple[tuple[Fraction, ...], ...]:
+ """Compute inverse matrix using constant bindings.
+
+ For each binding, we record how to map from target back to source.
+ The inverse matrix has shape (len(target), len(source)).
+ """
+ n_src = len(self.source)
+ n_tgt = len(self.target)
+
+ # Start with zero matrix of shape (n_tgt, n_src)
+ inv: list[list[Fraction]] = [
+ [Fraction(0)] * n_src for _ in range(n_tgt)
+ ]
+
+ # Build a set of source components that have bindings
+ bound_src_indices: set[int] = set()
+
+ # For each binding, record the reverse mapping
+ for binding in self.bindings:
+ src_idx = self.source.index(binding.source_component.name)
+ bound_src_indices.add(src_idx)
+
+ # The target expression tells us which target components
+ # contribute to this source component
+ for tgt_idx, coeff in enumerate(binding.target_expression.components):
+ if coeff != 0:
+ # For the inverse: we need to invert the coefficient
+ # E.g., if L → E⁻¹ (coeff = -1), then E → L⁻¹ means
+ # multiplying by -1 to get back the L exponent
+ inv[tgt_idx][src_idx] = Fraction(1) / coeff
+
+ # For source components without bindings, check for clean 1:1 mappings
+ for i in range(n_src):
+ if i in bound_src_indices:
+ continue
+
+ # Find which target component(s) this source maps to
+ non_zero = [
+ (j, self.matrix[i][j])
+ for j in range(n_tgt)
+ if self.matrix[i][j] != 0
+ ]
+
+ if len(non_zero) == 1 and non_zero[0][1] == 1:
+ # Clean 1-to-1 mapping: source[i] -> target[j] with coeff 1
+ j = non_zero[0][0]
+ inv[j][i] = Fraction(1)
+ elif len(non_zero) == 0:
+ # Source component is dropped (projected out) - leave as zero
+ pass
+ # For more complex mappings, we rely on bindings
+
+ return tuple(tuple(row) for row in inv)
+
+ def _invert_bindings(self) -> tuple[ConstantBinding, ...]:
+ """Invert the constant bindings (negate exponents, swap source/target)."""
+ inverted: list[ConstantBinding] = []
+
+ for binding in self.bindings:
+ src_idx = self.source.index(binding.source_component.name)
+
+ # Find the primary target component (first non-zero)
+ primary_tgt_idx = None
+ for tgt_idx, coeff in enumerate(binding.target_expression.components):
+ if coeff != 0:
+ primary_tgt_idx = tgt_idx
+ break
+
+ if primary_tgt_idx is None:
+ continue
+
+ # Create inverse binding: target component → source component
+ inv_binding = ConstantBinding(
+ source_component=self.target[primary_tgt_idx],
+ target_expression=Vector(
+ self.source,
+ tuple(
+ Fraction(1) if k == src_idx else Fraction(0)
+ for k in range(len(self.source))
+ ),
+ ),
+ constant_symbol=binding.constant_symbol,
+ exponent=-binding.exponent, # Negate the exponent
+ )
+ inverted.append(inv_binding)
+
+ return tuple(inverted)
+
+ def as_basis_transform(self) -> BasisTransform:
+ """Return as plain BasisTransform (loses binding information).
+
+ Useful for interoperability with code expecting BasisTransform.
+ """
+ return BasisTransform(
+ source=self.source,
+ target=self.target,
+ matrix=self.matrix,
+ )
+
+
# -----------------------------------------------------------------------------
# BasisGraph
# -----------------------------------------------------------------------------