diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index aeb0be86..edebb624 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -33,6 +33,7 @@ jobs:
       - uses: julia-actions/cache@v1
       - uses: julia-actions/julia-buildpkg@v1
       - name: "Run tests"
+        shell: bash
         run: |
           julia --color=yes --project=. -e 'import Pkg; Pkg.add("Coverage")'
           julia --color=yes --threads=auto --check-bounds=yes --depwarn=yes --code-coverage=user --project=. -e 'import Pkg; Pkg.test(coverage=true)'
@@ -40,6 +41,7 @@ jobs:
           julia --color=yes --project=. coverage.jl
       - name: "Coveralls"
         uses: coverallsapp/github-action@v2
+        if: matrix.version == '1'
         with:
           github-token: ${{ secrets.GITHUB_TOKEN }}
           file: lcov.info
diff --git a/Project.toml b/Project.toml
index d5d9a3e3..d86ddd63 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "0.5.0"
 
 [deps]
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Tricks = "410a4b4d-49e4-4fbc-ab6d-cb71b17b3775"
 
 [weakdeps]
diff --git a/docs/make.jl b/docs/make.jl
index da2dcedf..3e6f7114 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -41,6 +41,8 @@ makedocs(;
         "Home" => "index.md",
         "Utilities" => "api.md",
         "Units" => "units.md",
+        "Constants" => "constants.md",
+        "Symbolic Units" => "symbolic_units.md",
         "Types" => "types.md",
     ]
 )
diff --git a/docs/src/constants.md b/docs/src/constants.md
new file mode 100644
index 00000000..5bdf7e81
--- /dev/null
+++ b/docs/src/constants.md
@@ -0,0 +1,46 @@
+
+# Units
+
+Many common physical constants are available as well:
+
+```@docs
+Constants.c
+Constants.h
+Constants.hbar
+Constants.e
+Constants.k_B
+Constants.N_A
+Constants.eV
+Constants.R
+Constants.F
+Constants.sigma_sb
+Constants.alpha
+Constants.u
+Constants.G
+Constants.mu_0
+Constants.eps_0
+Constants.m_e
+Constants.m_p
+Constants.m_n
+Constants.a_0
+Constants.k_e
+Constants.Ryd
+```
+
+## Astronomical constants
+
+```@docs
+Constants.M_earth
+Constants.M_sun
+Constants.M_jup
+Constants.R_earth
+Constants.R_jup
+Constants.R_sun
+Constants.L_sun
+Constants.L_bol0
+Constants.sigma_T
+Constants.au
+Constants.pc
+Constants.ly
+Constants.atm
+```
diff --git a/docs/src/symbolic_units.md b/docs/src/symbolic_units.md
new file mode 100644
index 00000000..6b17627a
--- /dev/null
+++ b/docs/src/symbolic_units.md
@@ -0,0 +1,21 @@
+# Symbolic Dimensions
+
+Whereas `u"..."` will automatically convert all units to the same
+base SI units, `us"..."` will not. This uses the `SymbolicDimensions`
+type, which is a subtype of `AbstractDimensions` that stores the
+dimensions symbolically. This is useful for keeping track of the
+original units and constants in a user-entered expression.
+
+The two main functions for working with symbolic
+units are `sym_uparse` and `us_str`:
+
+```@docs
+@us_str
+sym_uparse
+```
+
+To convert a quantity to its regular base SI units, use `expand_units`:
+
+```@docs
+expand_units
+```
diff --git a/docs/src/types.md b/docs/src/types.md
index 563e765d..1f8e2ff3 100644
--- a/docs/src/types.md
+++ b/docs/src/types.md
@@ -15,4 +15,10 @@ AbstractQuantity
 ```
 
 Note also that the `Quantity` object can take a custom `AbstractDimensions`
-as input, so there is often no need to subtype `AbstractQuantity` separately.
\ No newline at end of file
+as input, so there is often no need to subtype `AbstractQuantity` separately.
+
+Another type which subtypes `AbstractDimensions` is `SymbolicDimensions`:
+
+```@docs
+SymbolicDimensions
+```
diff --git a/docs/src/units.md b/docs/src/units.md
index 3e231f84..098bdbc5 100644
--- a/docs/src/units.md
+++ b/docs/src/units.md
@@ -16,7 +16,7 @@ in a namespace with all the units available.
 
 ```@docs
 Units.m
-Units.g
+Units.kg
 Units.s
 Units.A
 Units.K
@@ -24,6 +24,8 @@ Units.cd
 Units.mol
 ```
 
+### Derived units
+
 Several derived SI units are available as well:
 
 ```@docs
@@ -39,4 +41,4 @@ Units.Ω
 Units.T
 Units.L
 Units.bar
-```
\ No newline at end of file
+```
diff --git a/src/DynamicQuantities.jl b/src/DynamicQuantities.jl
index a636b171..2e82cabb 100644
--- a/src/DynamicQuantities.jl
+++ b/src/DynamicQuantities.jl
@@ -1,18 +1,25 @@
 module DynamicQuantities
 
+export Units, Constants
 export AbstractQuantity, AbstractDimensions
-export Quantity, Dimensions, DimensionError, ustrip, dimension, valid
+export Quantity, Dimensions, SymbolicDimensions, DimensionError
+export ustrip, dimension, valid
 export ulength, umass, utime, ucurrent, utemperature, uluminosity, uamount
-export uparse, @u_str
+export uparse, @u_str, sym_uparse, @us_str, expand_units
 
 include("fixed_rational.jl")
 include("types.jl")
 include("utils.jl")
 include("math.jl")
 include("units.jl")
+include("constants.jl")
+include("uparse.jl")
+include("symbolic_dimensions.jl")
 
 import Requires: @init, @require
-import .Units: uparse, @u_str
+import .Units
+import .Constants
+import .UnitsParse: uparse, @u_str
 
 if !isdefined(Base, :get_extension)
     @init @require Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" include("../ext/DynamicQuantitiesUnitfulExt.jl")
diff --git a/src/constants.jl b/src/constants.jl
new file mode 100644
index 00000000..736a71b4
--- /dev/null
+++ b/src/constants.jl
@@ -0,0 +1,224 @@
+module Constants
+
+import ..DEFAULT_QUANTITY_TYPE
+import ..Quantity
+import ..Units as U
+import ..Units: _add_prefixes
+
+const _CONSTANT_SYMBOLS = Symbol[]
+const _CONSTANT_VALUES = DEFAULT_QUANTITY_TYPE[]
+
+macro register_constant(name, value)
+    return esc(_register_constant(name, value))
+end
+
+macro add_prefixes(base_unit, prefixes)
+    @assert prefixes.head == :tuple
+    return esc(_add_prefixes(base_unit, prefixes.args, _register_constant))
+end
+
+function _register_constant(name::Symbol, value)
+    s = string(name)
+    return quote
+        const $name = $value
+        push!(_CONSTANT_SYMBOLS, Symbol($s))
+        push!(_CONSTANT_VALUES, $name)
+    end
+end
+
+# Source: http://physics.nist.gov/constants (2018)
+
+# Exact, base:
+@register_constant c 299792458 * U.m/U.s
+@register_constant h 6.62607015e-34 * U.J/U.Hz
+@register_constant hbar h / (2π)
+@register_constant e 1.602176634e-19 * U.C
+@register_constant k_B 1.380649e-23 * U.J/U.K
+@register_constant N_A 6.02214076e+23 / U.mol
+
+@doc(
+    "Speed of light in a vacuum. Standard.",
+    c,
+)
+@doc(
+    "Planck constant. Standard.",
+    h,
+)
+@doc(
+    "Reduced Planck constant (h/2π). Standard.",
+    hbar,
+)
+@doc(
+    "Elementary charge. Standard.",
+    e,
+)
+@doc(
+    "Boltzmann constant. Standard.",
+    k_B,
+)
+@doc(
+    "Avogadro constant. Standard.",
+    N_A,
+)
+
+# Exact, derived:
+@register_constant eV e * U.J/U.C
+@register_constant R N_A * k_B
+@register_constant F N_A * e
+@register_constant sigma_sb (π^2/60) * k_B^4/(hbar^3 * c^2)
+
+@add_prefixes eV (m, k, M, G, T)
+
+@doc(
+    "Electron volt. Standard.",
+    eV,
+)
+@doc(
+    "Molar gas constant. Standard.",
+    R,
+)
+@doc(
+    "Faraday constant. Standard.",
+    F,
+)
+@doc(
+    "Stefan-Boltzmann constant. Standard.",
+    sigma_sb,
+)
+
+# Measured
+@register_constant alpha DEFAULT_QUANTITY_TYPE(7.2973525693e-3)
+@register_constant u 1.66053906660e-27 * U.kg
+@register_constant G 6.67430e-11 * U.m^3 / (U.kg * U.s^2)
+@register_constant mu_0 4π * alpha * hbar / (e^2 * c)
+@register_constant eps_0 8.8541878128e-12 * U.F/U.m
+@register_constant m_e 9.1093837015e-31 * U.kg
+@register_constant m_p 1.67262192369e-27 * U.kg
+@register_constant m_n 1.67492749804e-27 * U.kg
+@register_constant a_0 hbar/(m_e * c * alpha)
+@register_constant k_e 1/(4π * eps_0)
+@register_constant Ryd alpha^2 * m_e * c^2 / (2 * h)
+
+@doc(
+    "Fine-structure constant. Measured.",
+    alpha,
+)
+@doc(
+    "Atomic mass unit (1/12th the mass of Carbon-12). Measured.",
+    u,
+)
+@doc(
+    "Newtonian constant of gravitation. Measured.",
+    G,
+)
+@doc(
+    "Vacuum magnetic permeability. Measured.",
+    mu_0,
+)
+@doc(
+    "Vacuum electric permittivity. Measured.",
+    eps_0,
+)
+@doc(
+    "Electron mass. Measured.",
+    m_e,
+)
+@doc(
+    "Proton mass. Measured.",
+    m_p,
+)
+@doc(
+    "Neutron mass. Measured.",
+    m_n,
+)
+@doc(
+    "Bohr radius. Measured.",
+    a_0,
+)
+@doc(
+    "Coulomb constant (Note: SI units only!). Measured.",
+    k_e,
+)
+@doc(
+    "Rydberg frequency. Measured.",
+    Ryd,
+)
+
+# Astro constants.
+# Source: https://arxiv.org/abs/1510.07674
+
+@register_constant M_earth 5.97216787e+24 * U.kg
+@register_constant M_sun 1.98840987e+30 * U.kg
+@register_constant M_jup 1.8981246e+27 * U.kg
+@register_constant R_earth 6.3781e+6 * U.m
+@register_constant R_jup 7.1492e+7 * U.m
+@register_constant R_sun 6.957e+8 * U.m
+@register_constant L_sun 3.828e+26 * U.W
+@register_constant L_bol0 3.0128e+28 * U.W
+@register_constant sigma_T 6.6524587321e-29 * U.m^2
+@register_constant au 149597870700 * U.m
+@register_constant pc (648000/π) * au
+@register_constant ly c * U.yr
+@register_constant atm 101325 * U.Pa
+
+@add_prefixes pc (k, M, G)
+
+@doc(
+    "Earth mass. Measured.",
+    M_earth,
+)
+@doc(
+    "Solar mass. Measured.",
+    M_sun,
+)
+@doc(
+    "Jupiter mass. Measured.",
+    M_jup,
+)
+@doc(
+    "Nominal Earth equatorial radius. Standard.",
+    R_earth,
+)
+@doc(
+    "Nominal Jupiter equatorial radius. Standard.",
+    R_jup,
+)
+@doc(
+    "Nominal solar radius. Standard.",
+    R_sun,
+)
+@doc(
+    "Nominal solar luminosity. Standard.",
+    L_sun,
+)
+@doc(
+    "Standard luminosity at absolute bolometric magnitude 0. Standard.",
+    L_bol0,
+)
+@doc(
+    "Thomson scattering cross-section. Measured.",
+    sigma_T,
+)
+@doc(
+    "Astronomical unit. Standard.",
+    au,
+)
+@doc(
+    "Parsec. Standard.",
+    pc,
+)
+@doc(
+    "Light year. Standard.",
+    ly,
+)
+@doc(
+    "Standard atmosphere. Standard.",
+    atm,
+)
+
+"""A tuple of all possible constants."""
+const CONSTANT_SYMBOLS = Tuple(_CONSTANT_SYMBOLS)
+const CONSTANT_VALUES = Tuple(_CONSTANT_VALUES)
+const CONSTANT_MAPPING = NamedTuple([s => i for (i, s) in enumerate(CONSTANT_SYMBOLS)])
+
+end
diff --git a/src/symbolic_dimensions.jl b/src/symbolic_dimensions.jl
new file mode 100644
index 00000000..50797ffe
--- /dev/null
+++ b/src/symbolic_dimensions.jl
@@ -0,0 +1,214 @@
+import .Units: UNIT_SYMBOLS, UNIT_MAPPING, UNIT_VALUES
+import .Constants: CONSTANT_SYMBOLS, CONSTANT_MAPPING, CONSTANT_VALUES
+import SparseArrays as SA
+
+const SYMBOL_CONFLICTS = intersect(UNIT_SYMBOLS, CONSTANT_SYMBOLS)
+
+# Prefer units over constants:
+# For example, this means we can't have a symbolic Planck's constant,
+# as it is just "hours" (h), which is more common.
+const ALL_SYMBOLS = (
+    UNIT_SYMBOLS...,
+    setdiff(CONSTANT_SYMBOLS, SYMBOL_CONFLICTS)...
+)
+const ALL_VALUES = vcat(
+    UNIT_VALUES...,
+    (
+        v
+        for (k, v) in zip(CONSTANT_SYMBOLS, CONSTANT_VALUES)
+        if k ∉ SYMBOL_CONFLICTS
+    )...
+)
+const ALL_MAPPING = NamedTuple([s => i for (i, s) in enumerate(ALL_SYMBOLS)])
+
+"""
+    SymbolicDimensions{R} <: AbstractDimensions{R}
+
+An `AbstractDimensions` with one dimension for every unit and constant symbol.
+This is to allow for lazily reducing to SI base units, whereas
+`Dimensions` is always in SI base units. Furthermore, `SymbolicDimensions`
+stores dimensions using a sparse vector for efficiency (since there
+are so many unit symbols).
+
+You can convert a quantity using `SymbolicDimensions` as its dimensions
+to one which uses `Dimensions` as its dimensions (i.e., base SI units)
+`expand_units`.
+"""
+struct SymbolicDimensions{R} <: AbstractDimensions{R}
+    _data::SA.SparseVector{R}
+
+    SymbolicDimensions(data::SA.SparseVector) = new{eltype(data)}(data)
+    SymbolicDimensions{_R}(data::SA.SparseVector) where {_R} = new{_R}(data)
+end
+
+static_fieldnames(::Type{<:SymbolicDimensions}) = ALL_SYMBOLS
+data(d::SymbolicDimensions) = getfield(d, :_data)
+Base.getproperty(d::SymbolicDimensions{R}, s::Symbol) where {R} = data(d)[ALL_MAPPING[s]]
+Base.getindex(d::SymbolicDimensions{R}, k::Symbol) where {R} = getproperty(d, k)
+constructor_of(::Type{<:SymbolicDimensions}) = SymbolicDimensions
+
+SymbolicDimensions{R}(d::SymbolicDimensions) where {R} = SymbolicDimensions{R}(data(d))
+(::Type{D})(; kws...) where {D<:SymbolicDimensions} = D(DEFAULT_DIM_BASE_TYPE; kws...)
+(::Type{D})(::Type{R}; kws...) where {R,D<:SymbolicDimensions} =
+    let constructor=constructor_of(D){R}
+        length(kws) == 0 && return constructor(SA.spzeros(R, length(ALL_SYMBOLS)))
+        I = [ALL_MAPPING[s] for s in keys(kws)]
+        V = [tryrationalize(R, v) for v in values(kws)]
+        data = SA.sparsevec(I, V, length(ALL_SYMBOLS))
+        return constructor(data)
+    end
+
+function Base.convert(::Type{Qout}, q::Quantity{<:Any,<:Dimensions}) where {T,D<:SymbolicDimensions,Qout<:Quantity{T,D}}
+    output = Qout(
+        convert(T, ustrip(q)),
+        D;
+        m=ulength(q),
+        kg=umass(q),
+        s=utime(q),
+        A=ucurrent(q),
+        K=utemperature(q),
+        cd=uluminosity(q),
+        mol=uamount(q),
+    )
+    SA.dropzeros!(data(dimension(output)))
+    return output
+end
+function Base.convert(::Type{Q}, q::Quantity{<:Any,<:SymbolicDimensions}) where {T,D<:Dimensions,Q<:Quantity{T,D}}
+    result = one(Q) * ustrip(q)
+    d = dimension(q)
+    for (idx, value) in zip(SA.findnz(data(d))...)
+        result = result * convert(Q, ALL_VALUES[idx]) ^ value
+    end
+    return result
+end
+
+"""
+    expand_units(q::Quantity{<:Any,<:SymbolicDimensions})
+
+Expand the symbolic units in a quantity to their base SI form.
+In other words, this converts a `Quantity` with `SymbolicDimensions`
+to one with `Dimensions`.
+"""
+function expand_units(q::Q) where {T,R,D<:SymbolicDimensions{R},Q<:Quantity{T,D}}
+    return convert(Quantity{T,Dimensions{R}}, q)
+end
+
+
+Base.copy(d::SymbolicDimensions) = SymbolicDimensions(copy(data(d)))
+Base.:(==)(l::SymbolicDimensions, r::SymbolicDimensions) = data(l) == data(r)
+Base.iszero(d::SymbolicDimensions) = iszero(data(d))
+Base.:*(l::SymbolicDimensions, r::SymbolicDimensions) = SymbolicDimensions(data(l) + data(r))
+Base.:/(l::SymbolicDimensions, r::SymbolicDimensions) = SymbolicDimensions(data(l) - data(r))
+Base.inv(d::SymbolicDimensions) = SymbolicDimensions(-data(d))
+_pow(l::SymbolicDimensions{R}, r::R) where {R} = SymbolicDimensions(data(l) * r)
+
+
+"""
+    SymbolicUnitsParse
+
+A separate module where each unit is treated as a separate dimension,
+to enable pretty-printing of units.
+"""
+module SymbolicUnitsParse
+
+    import ..UNIT_SYMBOLS
+    import ..CONSTANT_SYMBOLS
+    import ..SYMBOL_CONFLICTS
+    import ..SymbolicDimensions
+
+    import ...Quantity
+    import ...DEFAULT_VALUE_TYPE
+    import ...DEFAULT_DIM_BASE_TYPE
+
+    # Lazily create unit symbols (since there are so many)
+    module Constants
+        import ..CONSTANT_SYMBOLS
+        import ..SYMBOL_CONFLICTS
+        import ..SymbolicDimensions
+
+        import ..Quantity
+        import ..DEFAULT_VALUE_TYPE
+        import ..DEFAULT_DIM_BASE_TYPE
+
+        import ...Constants as EagerConstants
+
+        const CONSTANT_SYMBOLS_EXIST = Ref{Bool}(false)
+        const CONSTANT_SYMBOLS_LOCK = Threads.SpinLock()
+        function _generate_unit_symbols()
+            CONSTANT_SYMBOLS_EXIST[] || lock(CONSTANT_SYMBOLS_LOCK) do
+                CONSTANT_SYMBOLS_EXIST[] && return nothing
+                for unit in setdiff(CONSTANT_SYMBOLS, SYMBOL_CONFLICTS)
+                    @eval const $unit = Quantity(DEFAULT_VALUE_TYPE(1.0), SymbolicDimensions{DEFAULT_DIM_BASE_TYPE}; $(unit)=1)
+                end
+                # Evaluate conflicting symbols to non-symbolic form:
+                for unit in SYMBOL_CONFLICTS
+                    @eval const $unit = convert(Quantity{DEFAULT_VALUE_TYPE,SymbolicDimensions}, EagerConstants.$unit)
+                end
+                CONSTANT_SYMBOLS_EXIST[] = true
+            end
+            return nothing
+        end
+    end
+    import .Constants
+
+    const UNIT_SYMBOLS_EXIST = Ref{Bool}(false)
+    const UNIT_SYMBOLS_LOCK = Threads.SpinLock()
+    function _generate_unit_symbols()
+        UNIT_SYMBOLS_EXIST[] || lock(UNIT_SYMBOLS_LOCK) do
+            UNIT_SYMBOLS_EXIST[] && return nothing
+            for unit in UNIT_SYMBOLS
+                @eval const $unit = Quantity(DEFAULT_VALUE_TYPE(1.0), SymbolicDimensions{DEFAULT_DIM_BASE_TYPE}; $(unit)=1)
+            end
+            UNIT_SYMBOLS_EXIST[] = true
+        end
+        return nothing
+    end
+
+    """
+        sym_uparse(raw_string::AbstractString)
+
+    Parse a string containing an expression of units and return the
+    corresponding `Quantity` object with `Float64` value.
+    However, that unlike the regular `u"..."` macro, this macro uses
+    `SymbolicDimensions` for the dimension type, which means that all units and
+    constants are stored symbolically and will not automatically expand to SI
+    units. For example, `sym_uparse("km/s^2")` would be parsed to
+    `Quantity(1.0, SymbolicDimensions, km=1, s=-2)`.
+
+    Note that inside this expression, you also have access to the `Constants`
+    module. So, for example, `sym_uparse("Constants.c^2 * Hz^2")` would evaluate to
+    `Quantity(1.0, SymbolicDimensions, c=2, Hz=2)`. However, note that due to
+    namespace collisions, a few physical constants are not available.
+    """
+    function sym_uparse(raw_string::AbstractString)
+        _generate_unit_symbols()
+        Constants._generate_unit_symbols()
+        raw_result = eval(Meta.parse(raw_string))
+        return copy(as_quantity(raw_result))::Quantity{DEFAULT_VALUE_TYPE,SymbolicDimensions{DEFAULT_DIM_BASE_TYPE}}
+    end
+
+    as_quantity(q::Quantity) = q
+    as_quantity(x::Number) = Quantity(convert(DEFAULT_VALUE_TYPE, x), SymbolicDimensions{DEFAULT_DIM_BASE_TYPE})
+    as_quantity(x) = error("Unexpected type evaluated: $(typeof(x))")
+end
+
+import .SymbolicUnitsParse: sym_uparse
+
+"""
+    us"[unit expression]"
+
+Parse a string containing an expression of units and return the
+corresponding `Quantity` object with `Float64` value. However,
+unlike the regular `u"..."` macro, this macro uses `SymbolicDimensions`
+for the dimension type, which means that all units and constants
+are stored symbolically and will not automatically expand to SI units.
+For example, `us"km/s^2"` would be parsed to `Quantity(1.0, SymbolicDimensions, km=1, s=-2)`.
+
+Note that inside this expression, you also have access to the `Constants`
+module. So, for example, `us"Constants.c^2 * Hz^2"` would evaluate to
+`Quantity(1.0, SymbolicDimensions, c=2, Hz=2)`. However, note that due to
+namespace collisions, a few physical constants are not available.
+"""
+macro us_str(s)
+    return esc(SymbolicUnitsParse.sym_uparse(s))
+end
diff --git a/src/types.jl b/src/types.jl
index 78006836..7868d76c 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -121,6 +121,8 @@ end
 (::Type{Q})(q::AbstractQuantity) where {T,D<:AbstractDimensions,Q<:AbstractQuantity{T,D}} = constructor_of(Q)(convert(T, ustrip(q)), convert(D, dimension(q)))
 (::Type{Q})(q::AbstractQuantity) where {T,Q<:AbstractQuantity{T}} = constructor_of(Q)(convert(T, ustrip(q)), dimension(q))
 
+const DEFAULT_QUANTITY_TYPE = Quantity{DEFAULT_VALUE_TYPE, DEFAULT_DIM_TYPE}
+
 new_dimensions(::Type{D}, dims...) where {D<:AbstractDimensions} = constructor_of(D)(dims...)
 new_quantity(::Type{Q}, l, r) where {Q<:AbstractQuantity} = constructor_of(Q)(l, r)
 
diff --git a/src/units.jl b/src/units.jl
index bbba5ffa..644d7448 100644
--- a/src/units.jl
+++ b/src/units.jl
@@ -1,48 +1,55 @@
 module Units
 
-export uparse, @u_str
-
 import ..DEFAULT_DIM_TYPE
 import ..DEFAULT_VALUE_TYPE
+import ..DEFAULT_QUANTITY_TYPE
 import ..Quantity
 
 @assert DEFAULT_VALUE_TYPE == Float64 "`units.jl` must be updated to support a different default value type."
 
+const _UNIT_SYMBOLS = Symbol[]
+const _UNIT_VALUES = DEFAULT_QUANTITY_TYPE[]
+
+macro register_unit(name, value)
+    return esc(_register_unit(name, value))
+end
+
 macro add_prefixes(base_unit, prefixes)
     @assert prefixes.head == :tuple
-    expr = _add_prefixes(base_unit, prefixes.args)
-    return expr |> esc
+    return esc(_add_prefixes(base_unit, prefixes.args, _register_unit))
 end
 
-function _add_prefixes(base_unit::Symbol, prefixes)
+function _register_unit(name::Symbol, value)
+    s = string(name)
+    return quote
+        const $name = $value
+        push!(_UNIT_SYMBOLS, Symbol($s))
+        push!(_UNIT_VALUES, $name)
+    end
+end
+
+function _add_prefixes(base_unit::Symbol, prefixes, register_function)
     all_prefixes = (
         f=1e-15, p=1e-12, n=1e-9, μ=1e-6, u=1e-6, m=1e-3, c=1e-2, d=1e-1,
-        k=1e3, M=1e6, G=1e9
+        k=1e3, M=1e6, G=1e9, T=1e12
     )
     expr = Expr(:block)
     for (prefix, value) in zip(keys(all_prefixes), values(all_prefixes))
         prefix in prefixes || continue
         new_unit = Symbol(prefix, base_unit)
-        push!(expr.args, :(const $new_unit = $value * $base_unit))
+        push!(expr.args, register_function(new_unit, :($value * $base_unit)))
     end
     return expr
 end
 
 # SI base units
-"Length in meters. Available variants: `fm`, `pm`, `nm`, `μm` (/`um`), `cm`, `dm`, `mm`, `km`, `Mm`, `Gm`."
-const m = Quantity(1.0, length=1)
-"Mass in grams. Available variants: `μg` (/`ug`), `mg`, `kg`."
-const g = Quantity(1e-3, mass=1)
-"Time in seconds. Available variants: `fs`, `ps`, `ns`, `μs` (/`us`), `ms`, `min`, `h` (/`hr`), `day`, `yr`, `kyr`, `Myr`, `Gyr`."
-const s = Quantity(1.0, time=1)
-"Current in Amperes. Available variants: `nA`, `μA` (/`uA`), `mA`, `kA`."
-const A = Quantity(1.0, current=1)
-"Temperature in Kelvin. Available variant: `mK`."
-const K = Quantity(1.0, temperature=1)
-"Luminosity in candela. Available variant: `mcd`."
-const cd = Quantity(1.0, luminosity=1)
-"Amount in moles. Available variant: `mmol`."
-const mol = Quantity(1.0, amount=1)
+@register_unit m DEFAULT_QUANTITY_TYPE(1.0, length=1)
+@register_unit g DEFAULT_QUANTITY_TYPE(1e-3, mass=1)
+@register_unit s DEFAULT_QUANTITY_TYPE(1.0, time=1)
+@register_unit A DEFAULT_QUANTITY_TYPE(1.0, current=1)
+@register_unit K DEFAULT_QUANTITY_TYPE(1.0, temperature=1)
+@register_unit cd DEFAULT_QUANTITY_TYPE(1.0, luminosity=1)
+@register_unit mol DEFAULT_QUANTITY_TYPE(1.0, amount=1)
 
 @add_prefixes m (f, p, n, μ, u, c, d, m, k, M, G)
 @add_prefixes g (μ, u, m, k)
@@ -52,28 +59,47 @@ const mol = Quantity(1.0, amount=1)
 @add_prefixes cd (m,)
 @add_prefixes mol (m,)
 
+@doc(
+    "Length in meters. Available variants: `fm`, `pm`, `nm`, `μm` (/`um`), `cm`, `dm`, `mm`, `km`, `Mm`, `Gm`.",
+    m,
+)
+@doc(
+    "Mass in kilograms. Available variants: `μg` (/`ug`), `mg`, `g`.",
+    kg,
+)
+@doc(
+    "Time in seconds. Available variants: `fs`, `ps`, `ns`, `μs` (/`us`), `ms`, `min`, `h` (/`hr`), `day`, `yr`, `kyr`, `Myr`, `Gyr`.",
+    s,
+)
+@doc(
+    "Current in Amperes. Available variants: `nA`, `μA` (/`uA`), `mA`, `kA`.",
+    A,
+)
+@doc(
+    "Temperature in Kelvin. Available variant: `mK`.",
+    K,
+)
+@doc(
+    "Luminosity in candela. Available variant: `mcd`.",
+    cd,
+)
+@doc(
+    "Amount in moles. Available variant: `mmol`.",
+    mol,
+)
+
 # SI derived units
-"Frequency in Hertz. Available variants: `kHz`, `MHz`, `GHz`."
-const Hz = inv(s)
-"Force in Newtons."
-const N = kg * m / s^2
-"Pressure in Pascals. Available variant: `kPa`."
-const Pa = N / m^2
-"Energy in Joules. Available variant: `kJ`."
-const J = N * m
-"Power in Watts. Available variants: `kW`, `MW`, `GW`."
-const W = J / s
-"Charge in Coulombs."
-const C = A * s
-"Voltage in Volts. Available variants: `kV`, `MV`, `GV`."
-const V = W / A
-"Capacitance in Farads."
-const F = C / V
-"Resistance in Ohms. Available variant: `mΩ`. Also available is ASCII `ohm` (with variant `mohm`)."
-const Ω = V / A
-const ohm = Ω
-"Magnetic flux density in Teslas."
-const T = N / (A * m)
+@register_unit Hz inv(s)
+@register_unit N kg * m / s^2
+@register_unit Pa N / m^2
+@register_unit J N * m
+@register_unit W J / s
+@register_unit C A * s
+@register_unit V W / A
+@register_unit F C / V
+@register_unit Ω V / A
+@register_unit ohm Ω
+@register_unit T N / (A * m)
 
 @add_prefixes Hz (k, M, G)
 @add_prefixes N ()
@@ -87,13 +113,55 @@ const T = N / (A * m)
 @add_prefixes ohm (m,)
 @add_prefixes T ()
 
+# SI derived units
+@doc(
+    "Frequency in Hertz. Available variants: `kHz`, `MHz`, `GHz`.",
+    Hz,
+)
+@doc(
+    "Force in Newtons.",
+    N,
+)
+@doc(
+    "Pressure in Pascals. Available variant: `kPa`.",
+    Pa,
+)
+@doc(
+    "Energy in Joules. Available variant: `kJ`.",
+    J,
+)
+@doc(
+    "Power in Watts. Available variants: `kW`, `MW`, `GW`.",
+    W,
+)
+@doc(
+    "Charge in Coulombs.",
+    C,
+)
+@doc(
+    "Voltage in Volts. Available variants: `kV`, `MV`, `GV`.",
+    V,
+)
+@doc(
+    "Capacitance in Farads.",
+    F,
+)
+@doc(
+    "Resistance in Ohms. Available variant: `mΩ`. Also available is ASCII `ohm` (with variant `mohm`).",
+    Ω,
+)
+@doc(
+    "Magnetic flux density in Teslas.",
+    T,
+)
+
 # Common assorted units
 ## Time
-const min = 60 * s
-const h = 60 * min
-const hr = h
-const day = 24 * h
-const yr = 365.25 * day
+@register_unit min 60 * s
+@register_unit h 60 * min
+@register_unit hr h
+@register_unit day 24 * h
+@register_unit yr 365.25 * day
 
 @add_prefixes min ()
 @add_prefixes h ()
@@ -102,47 +170,34 @@ const yr = 365.25 * day
 @add_prefixes yr (k, M, G)
 
 ## Volume
-"Volume in liters. Available variants: `mL`, `dL`."
-const L = dm^3
+@register_unit L dm^3
 
 @add_prefixes L (m, d)
 
+@doc(
+    "Volume in liters. Available variants: `mL`, `dL`.",
+    L,
+)
+
 ## Pressure
-"Pressure in bars."
-const bar = 100 * kPa
+@register_unit bar 100 * kPa
 
 @add_prefixes bar ()
 
+@doc(
+    "Pressure in bars.",
+    bar,
+)
+
 # Do not wish to define Gaussian units, as it changes
 # some formulas. Safer to force user to work exclusively in one unit system.
 
 # Do not wish to define physical constants, as the number of symbols might lead to ambiguity.
 # The user should define these instead.
 
-"""
-    uparse(s::AbstractString)
-
-Parse a string containing an expression of units and return the
-corresponding `Quantity` object with `Float64` value. For example,
-`uparse("m/s")` would be parsed to `Quantity(1.0, length=1, time=-1)`.
-"""
-function uparse(s::AbstractString)
-    return as_quantity(eval(Meta.parse(s)))::Quantity{DEFAULT_VALUE_TYPE,DEFAULT_DIM_TYPE}
-end
-
-as_quantity(q::Quantity) = q
-as_quantity(x::Number) = Quantity(convert(DEFAULT_VALUE_TYPE, x), DEFAULT_DIM_TYPE)
-as_quantity(x) = error("Unexpected type evaluated: $(typeof(x))")
-
-"""
-    u"[unit expression]"
-
-Parse a string containing an expression of units and return the
-corresponding `Quantity` object with `Float64` value. For example,
-`u"km/s^2"` would be parsed to `Quantity(1000.0, length=1, time=-2)`.
-"""
-macro u_str(s)
-    return esc(uparse(s))
-end
+"""A tuple of all possible unit symbols."""
+const UNIT_SYMBOLS = Tuple(_UNIT_SYMBOLS)
+const UNIT_VALUES = Tuple(_UNIT_VALUES)
+const UNIT_MAPPING = NamedTuple([s => i for (i, s) in enumerate(UNIT_SYMBOLS)])
 
 end
diff --git a/src/uparse.jl b/src/uparse.jl
new file mode 100644
index 00000000..e6f22854
--- /dev/null
+++ b/src/uparse.jl
@@ -0,0 +1,56 @@
+module UnitsParse
+
+import ..Quantity
+import ..DEFAULT_DIM_TYPE
+import ..DEFAULT_VALUE_TYPE
+import ..Units: UNIT_SYMBOLS
+import ..Constants
+
+function _generate_units_import()
+    import_expr = :(import ..Units: _)
+    deleteat!(first(import_expr.args).args, 2)
+    for symb in UNIT_SYMBOLS
+        push!(first(import_expr.args).args, Expr(:., symb))
+    end
+    return import_expr
+end
+
+eval(_generate_units_import())
+
+"""
+    uparse(s::AbstractString)
+
+Parse a string containing an expression of units and return the
+corresponding `Quantity` object with `Float64` value. For example,
+`uparse("m/s")` would be parsed to `Quantity(1.0, length=1, time=-1)`.
+
+Note that inside this expression, you also have access to the `Constants`
+module. So, for example, `uparse("Constants.c^2 * Hz^2")` would evaluate to
+the quantity corresponding to the speed of light multiplied by Hertz,
+squared.
+"""
+function uparse(s::AbstractString)
+    return as_quantity(eval(Meta.parse(s)))::Quantity{DEFAULT_VALUE_TYPE,DEFAULT_DIM_TYPE}
+end
+
+as_quantity(q::Quantity) = q
+as_quantity(x::Number) = Quantity(convert(DEFAULT_VALUE_TYPE, x), DEFAULT_DIM_TYPE)
+as_quantity(x) = error("Unexpected type evaluated: $(typeof(x))")
+
+"""
+    u"[unit expression]"
+
+Parse a string containing an expression of units and return the
+corresponding `Quantity` object with `Float64` value. For example,
+`u"km/s^2"` would be parsed to `Quantity(1000.0, length=1, time=-2)`.
+
+Note that inside this expression, you also have access to the `Constants`
+module. So, for example, `u"Constants.c^2 * Hz^2"` would evaluate to
+the quantity corresponding to the speed of light multiplied by Hertz,
+squared.
+"""
+macro u_str(s)
+    return esc(uparse(s))
+end
+
+end
diff --git a/src/utils.jl b/src/utils.jl
index 0e89e8f4..d779f140 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -115,6 +115,8 @@ Base.convert(::Type{Q}, q::AbstractQuantity) where {T,D,Q<:AbstractQuantity{T,D}
 Base.convert(::Type{D}, d::AbstractDimensions) where {D<:AbstractDimensions} = d
 Base.convert(::Type{D}, d::AbstractDimensions) where {R,D<:AbstractDimensions{R}} = D(d)
 
+Base.copy(q::Q) where {Q<:AbstractQuantity} = new_quantity(Q, copy(ustrip(q)), copy(dimension(q)))
+
 """
     ustrip(q::AbstractQuantity)
 
diff --git a/test/unittests.jl b/test/unittests.jl
index cea8a5e5..b3799309 100644
--- a/test/unittests.jl
+++ b/test/unittests.jl
@@ -346,6 +346,14 @@ end
     @test_throws LoadError eval(:(u":x"))
 end
 
+@testset "Constants" begin
+    @test Constants.h * Constants.c / (1000.0u"nm") ≈ 1.9864458571489284e-19u"J"
+
+    # Compute period of Earth based on solar mass and semi-major axis:
+    a = u"Constants.au"
+    @test isapprox(sqrt(4π^2 * a^3 / (Constants.G * Constants.M_sun)), 1u"yr"; rtol=1e-3)
+end
+
 @testset "Additional tests of FixedRational" begin
     @test convert(Int64, FixedRational{Int64,1000}(2 // 1)) == 2
     @test convert(Int32, FixedRational{Int64,1000}(3 // 1)) == 3
@@ -406,3 +414,38 @@ end
     # But, we always need to use a quantity when mixing with mathematical operations:
     @test_throws ErrorException MyQuantity(0.1) + 0.1 * MyDimensions()
 end
+
+@testset "Symbolic dimensions" begin
+    q = 1.5us"km/s"
+    @test q == 1.5 * us"km" / us"s"
+    @test typeof(q) <: Quantity{Float64,<:SymbolicDimensions}
+    @test string(dimension(q)) == "s⁻¹ km"
+    @test expand_units(q) == 1.5u"km/s"
+    @test string(dimension(us"Constants.au^1.5")) == "au³ᐟ²"
+    @test string(dimension(expand_units(us"Constants.au^1.5"))) == "m³ᐟ²"
+    @test expand_units(2.3us"Constants.au^1.5") ≈ 2.3u"Constants.au^1.5"
+    @test iszero(dimension(us"1.0")) == true
+    @test expand_units(inv(us"Constants.au")) ≈ 1/u"Constants.au"
+    @test dimension(inv(us"s") * us"km") == dimension(us"km/s")
+    @test dimension(inv(us"s") * us"m") != dimension(us"km/s")
+    @test dimension(expand_units(inv(us"s") * us"m")) == dimension(expand_units(us"km/s"))
+
+    @test_throws ErrorException sym_uparse("'c'")
+
+    # For constants which have a namespace collision, the numerical expansion is used:
+    @test dimension(us"Constants.au")[:au] == 1
+    @test dimension(us"Constants.h")[:h] == 0
+    @test dimension(us"h")[:h] == 1
+
+    @test us"Constants.h" != us"h"
+    @test expand_units(us"Constants.h") == u"Constants.h"
+
+    # Actually expands to:
+    @test dimension(us"Constants.h")[:m] == 2
+    @test dimension(us"Constants.h")[:s] == -1
+    @test dimension(us"Constants.h")[:kg] == 1
+
+    # So the numerical value is different from other constants:
+    @test ustrip(us"Constants.h") == ustrip(u"Constants.h")
+    @test ustrip(us"Constants.au") != ustrip(u"Constants.au")
+end