From 5fa660361c4717aaa0f809b4ba730f1387efcca6 Mon Sep 17 00:00:00 2001
From: jamesrafe <jamesrafe10@gmail.com>
Date: Fri, 7 Aug 2026 21:24:07 +0700
Subject: [PATCH 01/10] feat: init binned spectrum and add axis conversions

---
 src/SpectrumBase.jl              | 13 ++++++++----
 src/spectrum_binned.jl           | 28 +++++++++++++++++++++++++
 src/transforms/redshift.jl       |  3 ++-
 src/transforms/spectral_axis.jl  | 36 ++++++++++++++++++++++++++++++++
 src/transforms/transforms.jl     |  1 +
 test/runtests.jl                 |  2 +-
 test/spectrum.jl                 | 32 ++++++++++++++++++++++++++++
 test/transforms/spectral_axis.jl | 29 +++++++++++++++++++++++++
 8 files changed, 138 insertions(+), 6 deletions(-)
 create mode 100644 src/spectrum_binned.jl
 create mode 100644 src/transforms/spectral_axis.jl
 create mode 100644 test/transforms/spectral_axis.jl

diff --git a/src/SpectrumBase.jl b/src/SpectrumBase.jl
index e5f2795..87be181 100644
--- a/src/SpectrumBase.jl
+++ b/src/SpectrumBase.jl
@@ -4,11 +4,12 @@ module SpectrumBase
 export AbstractSpectrum, Spectrum, spectrum, spectral_axis, flux_axis
 
 # AbstractSpectrum types
-export SingleSpectrum, IFUSpectrum, EchelleSpectrum
+export SingleSpectrum, IFUSpectrum, EchelleSpectrum, BinnedSpectrum
 
 # Transforms
 export SpectrumResampler, redden, redden!, deredden, deredden!
 export redshift, redshift!, doppler_shift, doppler_shift!
+export to_frequency, to_wavelength
 
 # Utilities
 export blackbody #, line_flux, equivalent_width
@@ -207,9 +208,9 @@ Unitful.unit(spec::AbstractSpectrum) = unit(eltype(spectral_axis(spec))), unit(e
 
 # Spectrum types and basic arithmetic
 include("spectrum_single.jl")
+include("spectrum_binned.jl")
 include("spectrum_echelle.jl")
 include("spectrum_ifu.jl")
-#include("spectrum_binned.jl")
 
 """
     spectrum(spectral_axis, flux_axis, [meta])
@@ -289,13 +290,17 @@ function spectrum(spectral_axis::AbstractMatrix{<:Real}, flux_axis::AbstractMatr
     Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds))
 end
 
+function spectrum(spectral_axis::AbstractMatrix{<:Real}, flux_axis::AbstractVector{<:Real}; kwds...)
+    Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds))
+end
+
 function spectrum(spectral_axis::AbstractVector{<:Unitful.Quantity}, flux_axis::AbstractVector{<:Unitful.Quantity}; kwds...)
-    @assert dimension(eltype(spectral_axis)) ∈ (u"𝐋", u"𝐋^2 * 𝐌 * 𝐓^-2") "spectral_axis not recognized as having dimensions of wavelength or energy."
+    @assert dimension(eltype(spectral_axis)) ∈ (u"𝐋", u"𝐓^-1", u"𝐋^2 * 𝐌 * 𝐓^-2") "spectral_axis not recognized as having dimensions of wavelength, frequency, or energy."
     Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds))
 end
 
 function spectrum(spectral_axis::AbstractMatrix{<:Unitful.Quantity}, flux_axis::AbstractMatrix{<:Unitful.Quantity}; kwds...)
-    @assert dimension(eltype(spectral_axis)) ∈ (u"𝐋", u"𝐋^2 * 𝐌 * 𝐓^-2") "spectral_axis not recognized as having dimensions of wavelength or energy."
+    @assert dimension(eltype(spectral_axis)) ∈ (u"𝐋", u"𝐓^-1", u"𝐋^2 * 𝐌 * 𝐓^-2") "spectral_axis not recognized as having dimensions of wavelength, frequency, or energy."
     Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds))
 end
 
diff --git a/src/spectrum_binned.jl b/src/spectrum_binned.jl
new file mode 100644
index 0000000..8b62ddf
--- /dev/null
+++ b/src/spectrum_binned.jl
@@ -0,0 +1,28 @@
+"""
+    BinnedSpectrum <: AbstractSpectrum
+
+An instance of [`Spectrum`](@ref) where the spectral axis is an ``m x 2`` matrix of bin edges and the flux axis is a 1D array.
+"""
+const BinnedSpectrum = Spectrum{S, F, 2, 1} where {S, F}
+
+function Base.getindex(spec::BinnedSpectrum, i::Int)
+    return Spectrum(spectral_axis(spec)[i:i, :], flux_axis(spec)[i:i], meta(spec))
+end
+
+function Base.getindex(spec::BinnedSpectrum, inds)
+    return Spectrum(spectral_axis(spec)[inds, :], flux_axis(spec)[inds], meta(spec))
+end
+
+Base.firstindex(spec::BinnedSpectrum) = firstindex(flux_axis(spec))
+Base.lastindex(spec::BinnedSpectrum) = lastindex(flux_axis(spec))
+
+function Base.show(io::IO, spec::BinnedSpectrum)
+    w = spectral_axis(spec)
+    f = flux_axis(spec)
+    println(io, "BinnedSpectrum($(eltype(w)), $(eltype(f)))")
+    println(io, "  spectral axis $(size(w)): ", w[1, 1], " .. ", w[end, 2])
+    println(io, "  flux axis $(size(f)): ", first(f), " .. ", last(f))
+    print(io, "  meta: ", meta(spec))
+end
+
+# rebinning -> important in the future, make sure rmf gets rebinned correctly or it could have different size
diff --git a/src/transforms/redshift.jl b/src/transforms/redshift.jl
index 762d3a5..98773d9 100644
--- a/src/transforms/redshift.jl
+++ b/src/transforms/redshift.jl
@@ -14,9 +14,10 @@ end
 # Resolve the (1 + z) stretch factor to the scalar actually multiplied onto a spectral axis of element type `S`.
 _axis_factor(::Type{<:Real}, factor) = factor # Unitless: assume wavelength
 _axis_factor(::Type{<:Unitful.Length}, factor) = factor # Wavelength: stretches
+_axis_factor(::Type{<:Unitful.Frequency}, factor) = inv(factor) # Frequency: compresses
 _axis_factor(::Type{<:Unitful.Energy}, factor) = inv(factor) # Energy: compresses
 _axis_factor(::Type{S}, factor) where {S <: Unitful.Quantity} =
-    throw(ArgumentError("cannot shift a spectral axis with dimension $(dimension(S)); expected a wavelength or energy"))
+    throw(ArgumentError("cannot shift a spectral axis with dimension $(dimension(S)); expected a wavelength, frequency, or energy"))
 _axis_factor(spec::AbstractSpectrum, factor) = _axis_factor(eltype(spectral_axis(spec)), factor)
 
 # Return a copy of `spec` whose spectral axis is scaled by `factor`. Rebuilding
diff --git a/src/transforms/spectral_axis.jl b/src/transforms/spectral_axis.jl
new file mode 100644
index 0000000..b12a22a
--- /dev/null
+++ b/src/transforms/spectral_axis.jl
@@ -0,0 +1,36 @@
+function _reciprocal_axis(spec::Spectrum, target_unit)
+    spec isa BinnedSpectrum && throw(ArgumentError("wavelength-frequency conversion is not implemented for binned spectra"))
+
+    axis = Unitful.uconvert.(target_unit, c_0 ./ spectral_axis(spec))
+    return Spectrum(axis, copy(flux_axis(spec)), deepcopy(meta(spec)))
+end
+
+"""
+    to_frequency(spec; unit=u"Hz")
+
+Convert the wavelength axis of a spectrum to frequency using ``ν = c / λ``.
+Flux values are preserved.
+"""
+function to_frequency(spec::Spectrum{S}; unit = u"Hz") where {S <: Unitful.Length}
+    dimension(unit) == u"𝐓^-1" || throw(ArgumentError("target unit must have frequency dimensions, got $(dimension(unit))"))
+    return _reciprocal_axis(spec, unit)
+end
+
+function to_frequency(spec::AbstractSpectrum; unit = u"Hz")
+    throw(ArgumentError("to_frequency requires a Unitful wavelength axis"))
+end
+
+"""
+    to_wavelength(spec; unit=u"m")
+
+Convert the frequency axis of a spectrum to wavelength using ``λ = c / ν``.
+Flux values are preserved.
+"""
+function to_wavelength(spec::Spectrum{S}; unit = u"m") where {S <: Unitful.Frequency}
+    dimension(unit) == u"𝐋" || throw(ArgumentError("target unit must have wavelength dimensions, got $(dimension(unit))"))
+    return _reciprocal_axis(spec, unit)
+end
+
+function to_wavelength(spec::AbstractSpectrum; unit = u"m")
+    throw(ArgumentError("to_wavelength requires a Unitful frequency axis"))
+end
diff --git a/src/transforms/transforms.jl b/src/transforms/transforms.jl
index fcd974a..fc6c02b 100644
--- a/src/transforms/transforms.jl
+++ b/src/transforms/transforms.jl
@@ -1,5 +1,6 @@
 include("redden.jl")
 include("redshift.jl")
+include("spectral_axis.jl")
 
 # ## Broadening ops
 # include("kernels.jl")
diff --git a/test/runtests.jl b/test/runtests.jl
index 1814ca0..b80f70d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,7 +2,7 @@ using ParallelTestRunner: runtests, find_tests, parse_args
 import SpectrumBase
 
 const init_code = quote
-    using SpectrumBase: SpectrumBase, Spectrum, SingleSpectrum, EchelleSpectrum, IFUSpectrum, spectrum, spectral_axis, flux_axis
+    using SpectrumBase: SpectrumBase, Spectrum, SingleSpectrum, BinnedSpectrum, EchelleSpectrum, IFUSpectrum, spectrum, spectral_axis, flux_axis
     using Measurements: Measurements, ±
     using Unitful: @u_str, unit, ustrip
     import Random
diff --git a/test/spectrum.jl b/test/spectrum.jl
index b18af69..592daa0 100644
--- a/test/spectrum.jl
+++ b/test/spectrum.jl
@@ -18,6 +18,9 @@ Random.seed!(8675309)
     wave[3] = 10
     @test_throws ArgumentError spectrum(wave, flux)
     @test_throws ArgumentError spectrum(wave[begin:end-1], flux)
+
+    wave, flux = [1 2; 2 3; 3 4; 4 5; 5 6], 1:5
+    @test spectrum(wave, flux) isa BinnedSpectrum
 end
 
 @testset "Spectrum - Single" begin
@@ -62,6 +65,35 @@ end
     @test spec.name == "test spectrum"
 end
 
+@testset "Spectrum - Binned" begin
+    wave = [1.0 2.0; 2.0 3.0; 3.0 4.0; 4.0 5.0]
+    flux = [10.0, 20.0, 30.0, 40.0]
+
+    spec = spectrum(wave, flux, name = "test binned spectrum")
+    spec_i = spec[2]
+    spec_I = spec[2:3]
+
+    @test spec isa BinnedSpectrum
+    @test spec_i isa BinnedSpectrum
+    @test spec_I isa BinnedSpectrum
+    @test spectral_axis(spec_i) == [2.0 3.0]
+    @test flux_axis(spec_i) == [20.0]
+    @test spectral_axis(spec_I) == [2.0 3.0; 3.0 4.0]
+    @test flux_axis(spec_I) == [20.0, 30.0]
+    @test firstindex(spec) == 1
+    @test lastindex(spec) == 4
+    @test spec[firstindex(spec)] == spec[1]
+    @test spec[lastindex(spec)] == spec[4]
+
+    expected = """
+    BinnedSpectrum(Float64, Float64)
+      spectral axis (4, 2): 1.0 .. 5.0
+      flux axis (4,): 10.0 .. 40.0
+      meta: Dict{Symbol, Any}(:name => "test binned spectrum")"""
+    @test sprint(show, spec) == expected
+    @test spec.name == "test binned spectrum"
+end
+
 @testset "Spectrum - Echelle" begin
     n_wavs = 1000
     n_orders = 3
diff --git a/test/transforms/spectral_axis.jl b/test/transforms/spectral_axis.jl
new file mode 100644
index 0000000..415cf96
--- /dev/null
+++ b/test/transforms/spectral_axis.jl
@@ -0,0 +1,29 @@
+using SpectrumBase: to_frequency, to_wavelength, redshift
+using PhysicalConstants.CODATA2018: c_0
+using Unitful: uconvert
+using UnitfulAstro
+
+@testset "Wavelength and frequency conversion" begin
+    wavelength = [1.0, 1.5, 2.0]u"μm"
+    flux = [1.0, 2.0, 3.0]u"Jy"
+    spec = spectrum(wavelength, flux, name = "example")
+
+    frequency_spec = to_frequency(spec; unit = u"THz")
+    expected_frequency = uconvert.(u"THz", c_0 ./ wavelength)
+
+    @test spectral_axis(frequency_spec) ≈ expected_frequency
+    @test flux_axis(frequency_spec) == flux
+    @test frequency_spec.name == "example"
+    @test spectral_axis(spec) == wavelength
+
+    wavelength_spec = to_wavelength(frequency_spec; unit = u"μm")
+    @test spectral_axis(wavelength_spec) ≈ wavelength
+    @test flux_axis(wavelength_spec) == flux
+
+    shifted = redshift(frequency_spec, 1.0)
+    @test spectral_axis(shifted) ≈ expected_frequency ./ 2
+
+    @test_throws ArgumentError to_frequency(spec; unit = u"μm")
+    @test_throws ArgumentError to_wavelength(frequency_spec; unit = u"GHz")
+    @test_throws ArgumentError to_frequency(spectrum([1.0, 2.0], [3.0, 4.0]))
+end

From 0ab78c9359a6e007f955f4c8b3d32dc6276e4046 Mon Sep 17 00:00:00 2001
From: jamesrafe <jamesrafe10@gmail.com>
Date: Fri, 14 Aug 2026 01:57:47 +0700
Subject: [PATCH 02/10] test: add test coverage

---
 test/transforms/spectral_axis.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/transforms/spectral_axis.jl b/test/transforms/spectral_axis.jl
index 45332e1..b68fcb8 100644
--- a/test/transforms/spectral_axis.jl
+++ b/test/transforms/spectral_axis.jl
@@ -23,4 +23,5 @@ using SpectrumBase: to_frequency, to_wavelength, redshift
     @test_throws ArgumentError to_frequency(spec; unit = u"μm")
     @test_throws ArgumentError to_wavelength(frequency_spec; unit = u"GHz")
     @test_throws ArgumentError to_frequency(spectrum([1.0, 2.0], [3.0, 4.0]))
+    @test_throws ArgumentError to_wavelength(spectrum([1.0, 2.0], [3.0, 4.0]))
 end

From dbdfadf908d1a26f0d3a275fecb19d3b43ca9b98 Mon Sep 17 00:00:00 2001
From: Ian Weaver <weaveric@gmail.com>
Date: Mon, 17 Aug 2026 05:18:39 -0400
Subject: [PATCH 03/10] update redshift

---
 src/transforms/redshift.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/transforms/redshift.jl b/src/transforms/redshift.jl
index e2265fe..dfc9d88 100644
--- a/src/transforms/redshift.jl
+++ b/src/transforms/redshift.jl
@@ -14,9 +14,10 @@ end
 # Resolve the (1 + z) stretch factor to the scalar actually multiplied onto a spectral axis of element type `S`.
 _spectral_axis_factor(::Type{<:Real}, factor) = factor # Unitless: assume wavelength
 _spectral_axis_factor(::Type{<:Unitful.Length}, factor) = factor # Wavelength: stretches
+_spectral_axis_factor(::Type{<:Unitful.Frequency}, factor) = inv(factor) # Frequency: compresses
 _spectral_axis_factor(::Type{<:Unitful.Energy}, factor) = inv(factor) # Energy: compresses
 _spectral_axis_factor(::Type{S}, factor) where {S <: Unitful.Quantity} =
-    throw(ArgumentError("cannot shift a spectral axis with dimension $(dimension(S)); expected a wavelength or energy"))
+    throw(ArgumentError("cannot shift a spectral axis with dimension $(dimension(S)); expected a wavelength, frequency, or energy"))
 _spectral_axis_factor(spec::AbstractSpectrum, factor) = _spectral_axis_factor(eltype(spectral_axis(spec)), factor)
 
 # Rebuilding through the constructor lets the element type promote and re-runs the invariant checks.

From bfb14068ed0efe1e9b99d192d99c0c710cf9496b Mon Sep 17 00:00:00 2001
From: Ian Weaver <weaveric@gmail.com>
Date: Mon, 17 Aug 2026 05:18:56 -0400
Subject: [PATCH 04/10] unitful conversions

---
 Project.toml                     |  5 ++
 docs/src/transforms.md           | 27 +++++++++
 src/SpectrumBase.jl              | 11 +++-
 src/transforms/spectral_axis.jl  | 81 +++++++++++++++++---------
 test/transforms/spectral_axis.jl | 98 +++++++++++++++++++++++++++++---
 5 files changed, 185 insertions(+), 37 deletions(-)

diff --git a/Project.toml b/Project.toml
index c7706bb..850422c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -14,6 +14,10 @@ PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAstro = "6112ee07-acf9-5e0f-b108-d242c714bf9f"
+UnitfulEquivalences = "da9c4bc3-91c8-4f02-8a40-6b990d2a7e0c"
+
+[sources]
+UnitfulEquivalences = {url = "https://github.com/icweaver/UnitfulEquivalences.jl", rev = "spectral-density"}
 
 [compat]
 DustExtinction = ">=0.11.7"
@@ -23,4 +27,5 @@ PhysicalConstants = "0.2"
 RecipesBase = "1"
 Unitful = "1"
 UnitfulAstro = "1"
+UnitfulEquivalences = "0.2"
 julia = "1.10"
diff --git a/docs/src/transforms.md b/docs/src/transforms.md
index dbcdd4c..ccb1956 100644
--- a/docs/src/transforms.md
+++ b/docs/src/transforms.md
@@ -66,6 +66,33 @@ doppler_shift
 doppler_shift!
 ```
 
+## Spectral axis conversion
+
+`uconvert` from [Unitful.jl](https://github.com/PainterQubits/Unitful.jl) is extended to convert the spectral axis between wavelength, frequency, and photon energy via the photon equivalence ``E = hν = hc/λ``, powered by [UnitfulEquivalences.jl](https://github.com/sostock/UnitfulEquivalences.jl). Passing an `(axis_unit, flux_unit)` tuple, i.e., the shape returned by `unit(spec)`, also converts flux density values:
+
+```jldoctest
+julia> using SpectrumBase, Unitful, UnitfulAstro
+
+julia> spec = spectrum([1.0, 1.5, 2.0]u"μm", [1.0, 2.0, 3.0]u"W/m^2/μm");
+
+julia> uconvert(u"THz", spec);
+
+julia> νspec = uconvert((u"THz", u"Jy"), spec);
+
+julia> issorted(spectral_axis(νspec); rev = true)
+true
+
+julia> flux_axis(νspec)[1] ≈ uconvert(u"Jy", 1.0u"W/m^2/μm" * (1.0u"μm")^2 / Unitful.c0)
+true
+```
+
+### API/Reference
+
+```@docs
+Unitful.uconvert
+SpectralDensity
+```
+
 ## Resampling
 
 External interpolators, e.g., from [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl) or [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl), can be used to resample spectra onto a given wavelength grid. Starting with a sample spectrum `spec`, we first create a [`SpectrumResampler`](@ref) object `resampler` which stores the initial spectrum and interpolator `interp` together. We then apply this object to the wavelength grid of our choice to produce the resampled spectrum. We show example usage in the docstring below:
diff --git a/src/SpectrumBase.jl b/src/SpectrumBase.jl
index 5721426..cc2b9c3 100644
--- a/src/SpectrumBase.jl
+++ b/src/SpectrumBase.jl
@@ -9,7 +9,9 @@ export SingleSpectrum, IFUSpectrum, EchelleSpectrum, BinnedSpectrum
 # Transforms
 export SpectrumResampler, redden, redden!, deredden, deredden!
 export redshift, redshift!, doppler_shift, doppler_shift!
-export to_frequency, to_wavelength
+
+# Unit equivalences (re-exported from UnitfulEquivalences.jl)
+export SpectralDensity
 
 # Utilities
 export blackbody #, line_flux, equivalent_width
@@ -295,13 +297,16 @@ function spectrum(spectral_axis::AbstractMatrix{<:Real}, flux_axis::AbstractVect
     Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds))
 end
 
+# Recognized spectral-axis dimensions: wavelength, frequency, photon energy.
+const SPECTRAL_DIMENSIONS = (u"𝐋", u"𝐓^-1", u"𝐋^2 * 𝐌 * 𝐓^-2")
+
 function spectrum(spectral_axis::AbstractVector{<:Unitful.Quantity}, flux_axis::AbstractVector{<:Unitful.Quantity}; kwds...)
-    @assert dimension(eltype(spectral_axis)) ∈ (u"𝐋", u"𝐓^-1", u"𝐋^2 * 𝐌 * 𝐓^-2") "spectral_axis not recognized as having dimensions of wavelength, frequency, or energy."
+    @assert dimension(eltype(spectral_axis)) ∈ SPECTRAL_DIMENSIONS "spectral_axis not recognized as having dimensions of wavelength, frequency, or energy."
     Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds))
 end
 
 function spectrum(spectral_axis::AbstractMatrix{<:Unitful.Quantity}, flux_axis::AbstractMatrix{<:Unitful.Quantity}; kwds...)
-    @assert dimension(eltype(spectral_axis)) ∈ (u"𝐋", u"𝐓^-1", u"𝐋^2 * 𝐌 * 𝐓^-2") "spectral_axis not recognized as having dimensions of wavelength, frequency, or energy."
+    @assert dimension(eltype(spectral_axis)) ∈ SPECTRAL_DIMENSIONS "spectral_axis not recognized as having dimensions of wavelength, frequency, or energy."
     Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds))
 end
 
diff --git a/src/transforms/spectral_axis.jl b/src/transforms/spectral_axis.jl
index b12a22a..e2f7883 100644
--- a/src/transforms/spectral_axis.jl
+++ b/src/transforms/spectral_axis.jl
@@ -1,36 +1,65 @@
-function _reciprocal_axis(spec::Spectrum, target_unit)
-    spec isa BinnedSpectrum && throw(ArgumentError("wavelength-frequency conversion is not implemented for binned spectra"))
+using UnitfulEquivalences: Spectral, SpectralDensity
 
-    axis = Unitful.uconvert.(target_unit, c_0 ./ spectral_axis(spec))
-    return Spectrum(axis, copy(flux_axis(spec)), deepcopy(meta(spec)))
-end
+# Same-dimension targets are a plain unit conversion, keeping conversions idempotent.
+_spectral_value(spectral_unit, x) =
+    dimension(spectral_unit) == dimension(x) ? Unitful.uconvert(spectral_unit, x) :
+    Unitful.uconvert(spectral_unit, x, Spectral())
 
-"""
-    to_frequency(spec; unit=u"Hz")
+_flux_value(flux_unit, x, at) =
+    dimension(flux_unit) == dimension(x) ? Unitful.uconvert(flux_unit, x) :
+    Unitful.uconvert(flux_unit, x, SpectralDensity(at))
 
-Convert the wavelength axis of a spectrum to frequency using ``ν = c / λ``.
-Flux values are preserved.
-"""
-function to_frequency(spec::Spectrum{S}; unit = u"Hz") where {S <: Unitful.Length}
-    dimension(unit) == u"𝐓^-1" || throw(ArgumentError("target unit must have frequency dimensions, got $(dimension(unit))"))
-    return _reciprocal_axis(spec, unit)
+_convert_spectral_axis(spec::Spectrum, spectral_unit) = _spectral_value.(spectral_unit, spectral_axis(spec))
+
+_convert_flux_axis(spec::Spectrum, ::Nothing) = copy(flux_axis(spec))
+function _convert_flux_axis(spec::Spectrum{<:Unitful.Quantity}, flux_unit::Unitful.Units)
+    spec isa BinnedSpectrum && throw(ArgumentError(
+        "flux density conversion is not defined for binned spectra: bin values are integrated quantities, so converting them requires a bin-integration convention"))
+    return _flux_value.(flux_unit, flux_axis(spec), spectral_axis(spec))
 end
 
-function to_frequency(spec::AbstractSpectrum; unit = u"Hz")
-    throw(ArgumentError("to_frequency requires a Unitful wavelength axis"))
+function _convert_spectrum(spec::Spectrum{<:Unitful.Quantity}, spectral_unit::Unitful.Units, flux_unit::Union{Nothing, Unitful.Units})
+    dimension(spectral_unit) ∈ SPECTRAL_DIMENSIONS || throw(ArgumentError(
+        "spectral axis unit must be a wavelength, frequency, or energy unit; got $spectral_unit"))
+    return Spectrum(_convert_spectral_axis(spec, spectral_unit), _convert_flux_axis(spec, flux_unit), deepcopy(meta(spec)))
 end
 
-"""
-    to_wavelength(spec; unit=u"m")
+@doc raw"""
+    Unitful.uconvert(spectral_unit::Unitful.Units, spec::Spectrum)
+    Unitful.uconvert(units::Tuple{Unitful.Units, Unitful.Units}, spec::Spectrum)
+
+Convert the spectral axis of `spec` to `spectral_unit`, which may be any wavelength, frequency, or
+photon-energy unit. Coordinates are related by the photon equivalence ``E = hν = hc/λ``.
+A target of the same dimension as the current axis is a plain unit conversion.
+
+The conversion is applied elementwise, so a reciprocal conversion (e.g., wavelength -->
+frequency) yields an axis with the opposite direction of monotonicity.
+
+Given a single unit, flux values are copied unchanged. Beware that a flux density then
+remains "per" its original coordinate (e.g., a per-wavelength density against a frequency
+axis). Pass a `(spectral_unit, flux_unit)` tuple, i.e., the shape returned by `unit(spec)`, to
+also convert flux density values into the matching convention via [`SpectralDensity`](@ref),
+which preserves integrals: ``F_ν = F_λ λ^2 / c``.
 
-Convert the frequency axis of a spectrum to wavelength using ``λ = c / ν``.
-Flux values are preserved.
+# Examples
+
+```jldoctest
+julia> using Unitful, UnitfulAstro
+
+julia> spec = spectrum([1.0, 1.5, 2.0]u"μm", [1.0, 2.0, 3.0]u"W/m^2/μm");
+
+julia> νspec = uconvert((u"THz", u"Jy"), spec);
+
+julia> issorted(spectral_axis(νspec); rev=true)
+true
+
+julia> flux_axis(νspec)[1] ≈ uconvert(u"Jy", 1.0u"W/m^2/μm" * (1.0u"μm")^2 / Unitful.c0)
+true
+```
 """
-function to_wavelength(spec::Spectrum{S}; unit = u"m") where {S <: Unitful.Frequency}
-    dimension(unit) == u"𝐋" || throw(ArgumentError("target unit must have wavelength dimensions, got $(dimension(unit))"))
-    return _reciprocal_axis(spec, unit)
-end
+Unitful.uconvert(spectral_unit::Unitful.Units, spec::Spectrum{<:Unitful.Quantity}) = _convert_spectrum(spec, spectral_unit, nothing)
+Unitful.uconvert(units::Tuple{Unitful.Units, Unitful.Units}, spec::Spectrum{<:Unitful.Quantity}) = _convert_spectrum(spec, units[1], units[2])
 
-function to_wavelength(spec::AbstractSpectrum; unit = u"m")
-    throw(ArgumentError("to_wavelength requires a Unitful frequency axis"))
-end
+const _UNITLESS_AXIS_MSG = "`uconvert` requires a spectral axis with Unitful units. This spectrum's axis is unitless. Attach units first, e.g., `spectrum(axis * u\"angstrom\", flux)`"
+Unitful.uconvert(::Unitful.Units, ::AbstractSpectrum) = throw(ArgumentError(_UNITLESS_AXIS_MSG))
+Unitful.uconvert(::Tuple{Unitful.Units, Unitful.Units}, ::AbstractSpectrum) = throw(ArgumentError(_UNITLESS_AXIS_MSG))
diff --git a/test/transforms/spectral_axis.jl b/test/transforms/spectral_axis.jl
index b68fcb8..d5a9e9e 100644
--- a/test/transforms/spectral_axis.jl
+++ b/test/transforms/spectral_axis.jl
@@ -1,27 +1,109 @@
-using SpectrumBase: to_frequency, to_wavelength, redshift
+using SpectrumBase: SpectralDensity, redshift
+
+const H_JS = 6.62607015e-34u"J*s" # Planck constant, exact by SI definition
 
 @testset "Wavelength and frequency conversion" begin
     wavelength = [1.0, 1.5, 2.0]u"μm"
     flux = [1.0, 2.0, 3.0]u"Jy"
     spec = spectrum(wavelength, flux, name = "example")
+    c = C_MPS * u"m/s"
 
-    frequency_spec = to_frequency(spec; unit = u"THz")
-    expected_frequency = uconvert.(u"THz", C_MPS * u"m/s" ./ wavelength)
+    frequency_spec = @inferred uconvert(u"THz", spec)
+    expected_frequency = uconvert.(u"THz", c ./ wavelength)
 
+    # The conversion is elementwise, so an ascending wavelength axis becomes descending
+    @test issorted(spectral_axis(frequency_spec); rev = true)
     @test spectral_axis(frequency_spec) ≈ expected_frequency
     @test flux_axis(frequency_spec) == flux
     @test frequency_spec.name == "example"
     @test spectral_axis(spec) == wavelength
 
-    wavelength_spec = to_wavelength(frequency_spec; unit = u"μm")
+    # Round trip restores the original axis and ordering
+    wavelength_spec = @inferred uconvert(u"μm", frequency_spec)
     @test spectral_axis(wavelength_spec) ≈ wavelength
     @test flux_axis(wavelength_spec) == flux
 
     shifted = redshift(frequency_spec, 1.0)
     @test spectral_axis(shifted) ≈ expected_frequency ./ 2
 
-    @test_throws ArgumentError to_frequency(spec; unit = u"μm")
-    @test_throws ArgumentError to_wavelength(frequency_spec; unit = u"GHz")
-    @test_throws ArgumentError to_frequency(spectrum([1.0, 2.0], [3.0, 4.0]))
-    @test_throws ArgumentError to_wavelength(spectrum([1.0, 2.0], [3.0, 4.0]))
+    # Same-dimension targets are plain unit conversions, so conversions are idempotent
+    @test spectral_axis(uconvert(u"GHz", frequency_spec)) ≈
+        uconvert.(u"GHz", spectral_axis(frequency_spec))
+    @test spectral_axis(uconvert(u"nm", spec)) ≈ uconvert.(u"nm", wavelength)
+
+    # A `unit(spec)` tuple converts one spectrum to another's units
+    matched = uconvert(unit(frequency_spec), spec)
+    @test spectral_axis(matched) == spectral_axis(frequency_spec)
+    @test flux_axis(matched) == flux_axis(frequency_spec)
+
+    @test_throws MethodError uconvert(1u"Hz", spec) # a quantity, not a unit
+    @test_throws ArgumentError uconvert(u"kg", spec)
+    @test_throws ArgumentError uconvert(u"THz", spectrum([1.0, 2.0], [3.0, 4.0]))
+end
+
+@testset "Energy axis conversion" begin
+    energies = [1.0, 2.0, 5.0]u"keV"
+    flux = [1.0, 2.0, 3.0]u"Jy"
+    xray = spectrum(energies, flux)
+    c = C_MPS * u"m/s"
+
+    λspec = @inferred uconvert(u"angstrom", xray)
+    @test spectral_axis(λspec) ≈ uconvert.(u"angstrom", H_JS * c ./ energies)
+    @test flux_axis(λspec) == flux
+
+    # E → ν is linear, so the axis direction is preserved
+    νspec = @inferred uconvert(u"Hz", xray)
+    @test spectral_axis(νspec) ≈ uconvert.(u"Hz", energies ./ H_JS)
+    @test flux_axis(νspec) == flux
+
+    back = @inferred uconvert(u"keV", λspec)
+    @test spectral_axis(back) ≈ energies
+    @test flux_axis(back) == flux
+end
+
+@testset "Flux density conversion" begin
+    wavelength = [1.0, 1.5, 2.0]u"μm"
+    F_λ = [1.0, 2.0, 3.0]u"W/m^2/μm"
+    spec = spectrum(wavelength, F_λ)
+    c = C_MPS * u"m/s"
+
+    # Pointwise, SpectralDensity relates F_λ and F_ν at a given spectral coordinate,
+    # which may be given in any of the equivalent coordinates
+    @test uconvert(u"Jy", 1.0u"W/m^2/μm", SpectralDensity(2.0u"μm")) ≈
+        uconvert(u"Jy", 1.0u"W/m^2/μm" * (2.0u"μm")^2 / c)
+    @test uconvert(u"Jy", 1.0u"W/m^2/μm", SpectralDensity(uconvert(u"THz", c / 2.0u"μm"))) ≈
+        uconvert(u"Jy", 1.0u"W/m^2/μm", SpectralDensity(2.0u"μm"))
+
+    converted = @inferred uconvert((u"THz", u"Jy"), spec)
+    expected_F_ν = uconvert.(u"Jy", F_λ .* wavelength .^ 2 ./ c)
+    @test spectral_axis(converted) ≈ uconvert.(u"THz", c ./ wavelength)
+    @test flux_axis(converted) ≈ expected_F_ν
+
+    # Round trip back to F_λ
+    back = uconvert((u"μm", u"W/m^2/μm"), converted)
+    @test spectral_axis(back) ≈ wavelength
+    @test flux_axis(back) ≈ F_λ
+
+    # Per-energy density, the X-ray convention
+    F_E_spec = @inferred uconvert((u"keV", u"W/m^2/keV"), spec)
+    @test flux_axis(F_E_spec) ≈ uconvert.(u"W/m^2/keV", expected_F_ν ./ H_JS)
+
+    # A same-dimension flux target is a plain unit conversion
+    same = uconvert((u"THz", u"W/m^2/μm"), spec)
+    @test flux_axis(same) == F_λ
+end
+
+@testset "Binned spectrum axis conversion" begin
+    edges = [1.0 2.0; 2.0 3.0; 3.0 4.0]u"μm"
+    flux = [10.0, 20.0, 30.0]u"Jy"
+    binned = SpectrumBase.Spectrum(edges, flux, Dict{Symbol, Any}())
+    c = C_MPS * u"m/s"
+
+    converted = uconvert(u"THz", binned)
+    @test spectral_axis(converted) ≈ uconvert.(u"THz", c ./ edges)
+    @test flux_axis(converted) == flux
+    @test all(col -> issorted(col; rev = true), eachcol(spectral_axis(converted)))
+
+    # Bin values are integrated quantities; flux density conversion is undefined for them
+    @test_throws ArgumentError uconvert((u"THz", u"W/m^2/μm"), binned)
 end

From 114bcc593a1bf23cf16ba5409e2a5f7952b2089f Mon Sep 17 00:00:00 2001
From: Ian Weaver <weaveric@gmail.com>
Date: Mon, 17 Aug 2026 05:55:58 -0400
Subject: [PATCH 05/10] docs: cleanup

---
 docs/src/spectrum.md            | 1 +
 docs/src/transforms.md          | 3 +--
 src/spectrum_binned.jl          | 2 +-
 src/transforms/spectral_axis.jl | 4 ++--
 4 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/docs/src/spectrum.md b/docs/src/spectrum.md
index 315b621..297ef38 100644
--- a/docs/src/spectrum.md
+++ b/docs/src/spectrum.md
@@ -12,6 +12,7 @@ SpectrumBase.Spectrum
 SpectrumBase.SingleSpectrum
 SpectrumBase.EchelleSpectrum
 SpectrumBase.IFUSpectrum
+SpectrumBase.BinnedSpectrum
 ```
 
 ## Constructors
diff --git a/docs/src/transforms.md b/docs/src/transforms.md
index ccb1956..574c6d9 100644
--- a/docs/src/transforms.md
+++ b/docs/src/transforms.md
@@ -68,7 +68,7 @@ doppler_shift!
 
 ## Spectral axis conversion
 
-`uconvert` from [Unitful.jl](https://github.com/PainterQubits/Unitful.jl) is extended to convert the spectral axis between wavelength, frequency, and photon energy via the photon equivalence ``E = hν = hc/λ``, powered by [UnitfulEquivalences.jl](https://github.com/sostock/UnitfulEquivalences.jl). Passing an `(axis_unit, flux_unit)` tuple, i.e., the shape returned by `unit(spec)`, also converts flux density values:
+`uconvert` from [Unitful.jl](https://github.com/JuliaPhysics/Unitful.jl) is extended to convert the spectral axis between wavelength, frequency, and photon energy via the photon equivalence ``E = hν = hc/λ``, powered by [UnitfulEquivalences.jl](https://github.com/sostock/UnitfulEquivalences.jl). Passing an `(axis_unit, flux_unit)` tuple, i.e., the shape returned by `unit(spec)`, also converts flux density values:
 
 ```jldoctest
 julia> using SpectrumBase, Unitful, UnitfulAstro
@@ -90,7 +90,6 @@ true
 
 ```@docs
 Unitful.uconvert
-SpectralDensity
 ```
 
 ## Resampling
diff --git a/src/spectrum_binned.jl b/src/spectrum_binned.jl
index 8b62ddf..6366b4d 100644
--- a/src/spectrum_binned.jl
+++ b/src/spectrum_binned.jl
@@ -1,7 +1,7 @@
 """
     BinnedSpectrum <: AbstractSpectrum
 
-An instance of [`Spectrum`](@ref) where the spectral axis is an ``m x 2`` matrix of bin edges and the flux axis is a 1D array.
+An instance of [`Spectrum`](@ref) where the spectral axis is an ``m \\times 2`` matrix of bin edges, and the flux axis is a 1D array.
 """
 const BinnedSpectrum = Spectrum{S, F, 2, 1} where {S, F}
 
diff --git a/src/transforms/spectral_axis.jl b/src/transforms/spectral_axis.jl
index e2f7883..89b6c96 100644
--- a/src/transforms/spectral_axis.jl
+++ b/src/transforms/spectral_axis.jl
@@ -38,7 +38,7 @@ frequency) yields an axis with the opposite direction of monotonicity.
 Given a single unit, flux values are copied unchanged. Beware that a flux density then
 remains "per" its original coordinate (e.g., a per-wavelength density against a frequency
 axis). Pass a `(spectral_unit, flux_unit)` tuple, i.e., the shape returned by `unit(spec)`, to
-also convert flux density values into the matching convention via [`SpectralDensity`](@ref),
+also convert flux density values into the matching convention via [`UnitfulEquivalences.SpectralDensity`](https://github.com/sostock/UnitfulEquivalences.jl/pull/31),
 which preserves integrals: ``F_ν = F_λ λ^2 / c``.
 
 # Examples
@@ -50,7 +50,7 @@ julia> spec = spectrum([1.0, 1.5, 2.0]u"μm", [1.0, 2.0, 3.0]u"W/m^2/μm");
 
 julia> νspec = uconvert((u"THz", u"Jy"), spec);
 
-julia> issorted(spectral_axis(νspec); rev=true)
+julia> issorted(spectral_axis(νspec); rev = true)
 true
 
 julia> flux_axis(νspec)[1] ≈ uconvert(u"Jy", 1.0u"W/m^2/μm" * (1.0u"μm")^2 / Unitful.c0)

From ad249a959fc0647a16b1e827920c1ded9409e39b Mon Sep 17 00:00:00 2001
From: Ian Weaver <weaveric@gmail.com>
Date: Mon, 17 Aug 2026 06:22:26 -0400
Subject: [PATCH 06/10] docs: cleanup

---
 docs/src/transforms.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/src/transforms.md b/docs/src/transforms.md
index 574c6d9..83d9950 100644
--- a/docs/src/transforms.md
+++ b/docs/src/transforms.md
@@ -66,9 +66,9 @@ doppler_shift
 doppler_shift!
 ```
 
-## Spectral axis conversion
+## Axis conversions
 
-`uconvert` from [Unitful.jl](https://github.com/JuliaPhysics/Unitful.jl) is extended to convert the spectral axis between wavelength, frequency, and photon energy via the photon equivalence ``E = hν = hc/λ``, powered by [UnitfulEquivalences.jl](https://github.com/sostock/UnitfulEquivalences.jl). Passing an `(axis_unit, flux_unit)` tuple, i.e., the shape returned by `unit(spec)`, also converts flux density values:
+Powered by [UnitfulEquivalences.jl](https://github.com/sostock/UnitfulEquivalences.jl), [`uconvert`](@ref) from [Unitful.jl](https://github.com/JuliaPhysics/Unitful.jl) is extended to convert the spectral axis between wavelength, frequency, and photon energy, along with its associated flux density:
 
 ```jldoctest
 julia> using SpectrumBase, Unitful, UnitfulAstro

From 862f9dac8f1cad99235a2750021b5446c1e6fe9c Mon Sep 17 00:00:00 2001
From: Ian Weaver <weaveric@gmail.com>
Date: Mon, 17 Aug 2026 06:25:00 -0400
Subject: [PATCH 07/10] extra BinnedSpectrum logic up

---
 src/SpectrumBase.jl              | 22 +++++++++++++++++++---
 src/spectrum_binned.jl           | 13 ++++++++-----
 src/transforms/redden.jl         |  5 +++++
 test/spectrum.jl                 | 31 +++++++++++++++++++++++++++++++
 test/transforms/redden.jl        |  8 ++++++++
 test/transforms/spectral_axis.jl |  2 +-
 6 files changed, 72 insertions(+), 9 deletions(-)

diff --git a/src/SpectrumBase.jl b/src/SpectrumBase.jl
index cc2b9c3..f067708 100644
--- a/src/SpectrumBase.jl
+++ b/src/SpectrumBase.jl
@@ -57,9 +57,7 @@ mutable struct Spectrum{S<:Number, F<:Number, M, N, A<:AbstractArray{S, M}, B<:A
             - SingleSpectrum: spectral axis (M-length vector), flux axis (M-length vector)
             - EchelleSpectrum: spectral axis (M x N matrix), flux axis (M x N matrix)
             - IFUSpectrum: spectral axis (M-length vector), flux axis (M x N x K matrix)
-            - TODO: BinnedSpectrum (final name(s) tbd):
-                - energy (M × 2 matrix), flux (N-length vector)
-                - others?
+            - BinnedSpectrum: spectral axis (M x 2 matrix of bin edges), flux axis (M-length vector)
 
             See the documentation for each spectrum type for more.
             """))
@@ -71,6 +69,19 @@ mutable struct Spectrum{S<:Number, F<:Number, M, N, A<:AbstractArray{S, M}, B<:A
             throw(ArgumentError("Spectral axis must be strictly increasing or decreasing."))
         end
 
+        # BinnedSpectrum invariants: m × 2 bin edges, each row strictly ordered in the axis direction
+        if M == 2 && N == 1
+            if size(s, 2) != 2
+                throw(ArgumentError("BinnedSpectrum spectral axis must be an m × 2 matrix of bin edges; got size $(size(s))"))
+            end
+            rows = eachrow(s)
+            ordered = all(issorted, spec_ax) ? all(r -> first(r) < last(r), rows) :
+                all(r -> first(r) > last(r), rows)
+            if !ordered
+                throw(ArgumentError("BinnedSpectrum bin edges must be strictly ordered within each row, in the direction of the spectral axis."))
+            end
+        end
+
         return new{S, F, M, N, A, B}(s, f, meta)
     end
 end
@@ -310,6 +321,11 @@ function spectrum(spectral_axis::AbstractMatrix{<:Unitful.Quantity}, flux_axis::
     Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds))
 end
 
+function spectrum(spectral_axis::AbstractMatrix{<:Unitful.Quantity}, flux_axis::AbstractVector{<:Unitful.Quantity}; kwds...)
+    @assert dimension(eltype(spectral_axis)) ∈ SPECTRAL_DIMENSIONS "spectral_axis not recognized as having dimensions of wavelength, frequency, or energy."
+    Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds))
+end
+
 include("utils.jl")
 include("transforms/resampler.jl")
 include("transforms/transforms.jl")
diff --git a/src/spectrum_binned.jl b/src/spectrum_binned.jl
index 6366b4d..17978ea 100644
--- a/src/spectrum_binned.jl
+++ b/src/spectrum_binned.jl
@@ -5,9 +5,7 @@ An instance of [`Spectrum`](@ref) where the spectral axis is an ``m \\times 2``
 """
 const BinnedSpectrum = Spectrum{S, F, 2, 1} where {S, F}
 
-function Base.getindex(spec::BinnedSpectrum, i::Int)
-    return Spectrum(spectral_axis(spec)[i:i, :], flux_axis(spec)[i:i], meta(spec))
-end
+Base.getindex(spec::BinnedSpectrum, i::Integer) = spec[i:i]
 
 function Base.getindex(spec::BinnedSpectrum, inds)
     return Spectrum(spectral_axis(spec)[inds, :], flux_axis(spec)[inds], meta(spec))
@@ -20,8 +18,13 @@ function Base.show(io::IO, spec::BinnedSpectrum)
     w = spectral_axis(spec)
     f = flux_axis(spec)
     println(io, "BinnedSpectrum($(eltype(w)), $(eltype(f)))")
-    println(io, "  spectral axis $(size(w)): ", w[1, 1], " .. ", w[end, 2])
-    println(io, "  flux axis $(size(f)): ", first(f), " .. ", last(f))
+    if isempty(f)
+        println(io, "  spectral axis $(size(w)): empty")
+        println(io, "  flux axis $(size(f)): empty")
+    else
+        println(io, "  spectral axis $(size(w)): ", first(w), " .. ", last(w))
+        println(io, "  flux axis $(size(f)): ", first(f), " .. ", last(f))
+    end
     print(io, "  meta: ", meta(spec))
 end
 
diff --git a/src/transforms/redden.jl b/src/transforms/redden.jl
index a5cf031..b34f910 100644
--- a/src/transforms/redden.jl
+++ b/src/transforms/redden.jl
@@ -1,5 +1,10 @@
 import DustExtinction
 
+# Extinction laws evaluate at a single wavelength per sample, which bin edges cannot provide.
+const _BINNED_REDDEN_MSG = "reddening is not defined for binned spectra: applying an extinction law to bin edges requires a bin-integration convention"
+redden!(::BinnedSpectrum, Av; kwargs...) = throw(ArgumentError(_BINNED_REDDEN_MSG))
+deredden!(::BinnedSpectrum, Av; kwargs...) = throw(ArgumentError(_BINNED_REDDEN_MSG))
+
 """
     redden!(::AbstractSpectrum, Av; Rv = 3.1, law = DustExtinction.CCM89)
 
diff --git a/test/spectrum.jl b/test/spectrum.jl
index 490dc1f..303b24f 100644
--- a/test/spectrum.jl
+++ b/test/spectrum.jl
@@ -85,6 +85,22 @@ end
     @test spec[firstindex(spec)] == spec[1]
     @test spec[lastindex(spec)] == spec[4]
 
+    # Any integer index selects a single bin
+    @test spec[Int32(2)] == spec[2] == spec[2:2]
+
+    # m × 2 shape and row-order invariants
+    @test_throws ArgumentError spectrum([1.0 2.0 3.0; 2.0 3.0 4.0], [1.0, 2.0]) # not two columns
+    @test_throws ArgumentError spectrum(reshape([1.0, 2.0], 2, 1), [1.0, 2.0]) # not two columns
+    @test_throws ArgumentError spectrum([2.0 1.0; 3.0 2.0], [1.0, 2.0]) # edges swapped within rows
+    @test_throws ArgumentError spectrum([1.0 2.0; 3.0 2.5], [1.0, 2.0]) # inconsistent row order
+    @test_throws ArgumentError spectrum([1.0 1.0; 2.0 2.0], [1.0, 2.0]) # zero-width bins
+    @test spectrum([4.0 3.0; 3.0 2.0], [1.0, 2.0]) isa BinnedSpectrum # descending edges are valid
+
+    # Empty selections construct and display
+    empty_spec = spec[falses(4)]
+    @test length(empty_spec) == 0
+    @test occursin("empty", sprint(show, empty_spec))
+
     expected = """
     BinnedSpectrum(Float64, Float64)
       spectral axis (4, 2): 1.0 .. 5.0
@@ -353,3 +369,18 @@ end
         @test flux_axis(s) ≈ flux_axis(spec) ./ 10
     end
 end
+
+@testset "Unitful Spectrum - Binned" begin
+    edges = [1.0 2.0; 2.0 3.0; 3.0 4.0]u"keV"
+    flux = [1.0, 2.0, 3.0]u"Jy"
+
+    spec = spectrum(edges, flux, name = "test xray spectrum")
+
+    @test spec isa BinnedSpectrum
+    @test unit(spec) == (u"keV", u"Jy")
+    @test spectral_axis(spec) == edges
+    @test flux_axis(spec) == flux
+    @test spec.name == "test xray spectrum"
+
+    @test_throws AssertionError spectrum([1.0 2.0; 2.0 3.0]u"kg", flux)
+end
diff --git a/test/transforms/redden.jl b/test/transforms/redden.jl
index 88cd372..71fc2b6 100644
--- a/test/transforms/redden.jl
+++ b/test/transforms/redden.jl
@@ -67,3 +67,11 @@ end
         @test flux_axis(dereddened) ≈ flux_axis(spec)
     end
 end
+
+@testset "Binned spectra are rejected" begin
+    binned = spectrum([4000.0 4010.0; 4010.0 4020.0], [1.0, 2.0])
+    @test_throws ArgumentError redden(binned, 0.3)
+    @test_throws ArgumentError redden!(binned, 0.3)
+    @test_throws ArgumentError deredden(binned, 0.3)
+    @test_throws ArgumentError deredden!(binned, 0.3)
+end
diff --git a/test/transforms/spectral_axis.jl b/test/transforms/spectral_axis.jl
index d5a9e9e..76f9574 100644
--- a/test/transforms/spectral_axis.jl
+++ b/test/transforms/spectral_axis.jl
@@ -96,7 +96,7 @@ end
 @testset "Binned spectrum axis conversion" begin
     edges = [1.0 2.0; 2.0 3.0; 3.0 4.0]u"μm"
     flux = [10.0, 20.0, 30.0]u"Jy"
-    binned = SpectrumBase.Spectrum(edges, flux, Dict{Symbol, Any}())
+    binned = spectrum(edges, flux)
     c = C_MPS * u"m/s"
 
     converted = uconvert(u"THz", binned)

From 43526897866da2521892221af7825581cc380888 Mon Sep 17 00:00:00 2001
From: Ian Weaver <weaveric@gmail.com>
Date: Fri, 28 Aug 2026 14:11:44 -0700
Subject: [PATCH 08/10] additional redshift notes + checks

---
 docs/src/transforms.md      | 18 +++++++++++++++++-
 src/transforms/redshift.jl  | 14 +++++++++++---
 test/transforms/redshift.jl | 31 +++++++++++++++++++++++++++++++
 3 files changed, 59 insertions(+), 4 deletions(-)

diff --git a/docs/src/transforms.md b/docs/src/transforms.md
index 83d9950..ee1c248 100644
--- a/docs/src/transforms.md
+++ b/docs/src/transforms.md
@@ -43,7 +43,7 @@ julia> spectral_axis(shifted) ≈ spectral_axis(spec) .* 1.5
 true
 ```
 
-**Doppler shift** shifts by a radial velocity `v`. Pass a `Unitful` velocity or a plain number (interpreted as m/s). Set `relativistic=true` for the full relativistic formula:
+**Doppler shift** shifts by a radial velocity `v`. Pass a `Unitful` velocity or a plain number (interpreted as km/s). Set `relativistic=true` for the full relativistic formula:
 
 ```jldoctest
 julia> using SpectrumBase, Unitful
@@ -57,6 +57,22 @@ julia> shifted_rel = doppler_shift(spec, 100u"km/s"; relativistic=true);
 
 Both `redshift` and `doppler_shift` return a new spectrum. In-place variants `redshift!` and `doppler_shift!` are also available.
 
+Binned spectra are supported as well: both bin edges are shifted and the bin values are carried unchanged. An energy axis is compressed rather than stretched:
+
+```jldoctest
+julia> using SpectrumBase, Unitful
+
+julia> spec = spectrum([1.0 2.0; 2.0 3.0; 3.0 4.0]u"keV", [10.0, 20.0, 30.0]u"erg/s/cm^2");
+
+julia> shifted = redshift(spec, 0.5);
+
+julia> spectral_axis(shifted) ≈ spectral_axis(spec) ./ 1.5
+true
+
+julia> flux_axis(shifted) == flux_axis(spec)
+true
+```
+
 ### API/Reference
 
 ```@docs
diff --git a/src/transforms/redshift.jl b/src/transforms/redshift.jl
index dfc9d88..24629a2 100644
--- a/src/transforms/redshift.jl
+++ b/src/transforms/redshift.jl
@@ -11,6 +11,12 @@ function _doppler_factor(v; relativistic)
     return sqrt((1 + β) / (1 - β))
 end
 
+# Cosmological redshift
+function _redshift_factor(z)
+    z > -1 || throw(DomainError(z, "Cosmological redshift requires z > -1."))
+    return 1 + z
+end
+
 # Resolve the (1 + z) stretch factor to the scalar actually multiplied onto a spectral axis of element type `S`.
 _spectral_axis_factor(::Type{<:Real}, factor) = factor # Unitless: assume wavelength
 _spectral_axis_factor(::Type{<:Unitful.Length}, factor) = factor # Wavelength: stretches
@@ -36,7 +42,7 @@ Reassigns the spectral axis, so its element and array type cannot change.
 Use [`redshift`](@ref) for spectra with an integer spectral axis.
 """
 function redshift!(spec::AbstractSpectrum, z::Real)
-    spec.spectral_axis = spectral_axis(spec) .* _spectral_axis_factor(spec, 1 + z)
+    spec.spectral_axis = spectral_axis(spec) .* _spectral_axis_factor(spec, _redshift_factor(z))
     return spec
 end
 
@@ -53,6 +59,7 @@ The observed wavelength is related to the rest-frame wavelength by:
 
 where ``z`` is the cosmological redshift parameter. Only the spectral axis is transformed.
 Flux density values are not corrected for the stretching of the wavelength bins.
+For a [`BinnedSpectrum`](@ref), both bin edges are shifted and the bin values are carried unchanged.
 
 A wavelength axis is stretched by ``1 + z``. An energy axis (distinguished by its `Unitful` dimension)
 is compressed by the same factor. A unitless axis is assumed to be a wavelength.
@@ -60,7 +67,7 @@ is compressed by the same factor. A unitless axis is assumed to be a wavelength.
 # Arguments
 
 - `spec`: The input spectrum.
-- `z`: Redshift parameter. Positive values redshift (longer wavelengths),
+- `z`: Redshift parameter, must satisfy `z > -1`. Positive values redshift (longer wavelengths),
   negative values blueshift (shorter wavelengths).
 
 # Examples
@@ -76,7 +83,7 @@ true
 
 See also [`doppler_shift`](@ref) for velocity-based Doppler shifting.
 """
-redshift(spec::AbstractSpectrum, z::Real) = _scale_spectral_axis(spec, 1 + z)
+redshift(spec::AbstractSpectrum, z::Real) = _scale_spectral_axis(spec, _redshift_factor(z))
 
 """
     doppler_shift!(spec::AbstractSpectrum, v; relativistic=false)
@@ -109,6 +116,7 @@ Apply a Doppler shift to a spectrum, returning a new spectrum with shifted wavel
 ```
 
 Only the spectral axis is transformed; flux density values are preserved as-is.
+For a [`BinnedSpectrum`](@ref), both bin edges are shifted and the bin values are carried unchanged.
 
 # Arguments
 
diff --git a/test/transforms/redshift.jl b/test/transforms/redshift.jl
index f1c7e05..d4bedf8 100644
--- a/test/transforms/redshift.jl
+++ b/test/transforms/redshift.jl
@@ -1,4 +1,5 @@
 using SpectrumBase:
+    BinnedSpectrum,
     redshift,
     redshift!,
     doppler_shift,
@@ -81,6 +82,36 @@ end
         spec = SpectrumBase.Spectrum(axis, ones(10) * u"erg/s/cm^2/angstrom", Dict{Symbol, Any}())
         @test_throws ArgumentError redshift(spec, 0.5)
     end
+
+    @testset "Unphysical redshift rejected" begin
+        spec = mock_spectrum_redshift()
+        @test_throws DomainError redshift(spec, -1.0) # Exactly z = -1
+        @test_throws DomainError redshift(spec, -1.5)
+        @test_throws DomainError redshift!(spec, -1.0)
+        # Blueshifts above the bound remain valid
+        @test redshift(spec, -0.5) isa typeof(spec)
+    end
+
+    @testset "Binned spectra" begin
+        edges = [1.0 2.0; 2.0 3.0; 3.0 4.0]
+        flux = [10.0, 20.0, 30.0]
+        binned = spectrum(edges, flux)
+
+        shifted = @inferred redshift(binned, 0.5)
+        @test shifted isa BinnedSpectrum
+        @test spectral_axis(shifted) ≈ edges .* 1.5
+        @test flux_axis(shifted) == flux
+        @test spectral_axis(binned) == edges
+
+        @test redshift!(binned, 0.5) === binned
+        @test spectral_axis(binned) ≈ edges .* 1.5
+
+        # Energy axis compresses
+        kev = spectrum(edges * u"keV", flux * u"Jy")
+        @test spectral_axis(@inferred redshift(kev, 0.5)) ≈ edges * u"keV" ./ 1.5
+        @test spectral_axis(doppler_shift(kev, 1000.0)) ≈ edges * u"keV" ./ (1 + 1000 / C_KMPS)
+        @test flux_axis(redshift(kev, 0.5)) == flux_axis(kev)
+    end
 end
 
 @testset "doppler_shift" begin

From ae382babe2616e6f30bd99f9c451184f6539ae60 Mon Sep 17 00:00:00 2001
From: Ian Weaver <weaveric@gmail.com>
Date: Wed, 16 Sep 2026 17:33:05 -0700
Subject: [PATCH 09/10] refactor: `SpectralDensity` now upstreamed

---
 Project.toml                     | 5 +----
 src/SpectrumBase.jl              | 3 ---
 src/transforms/spectral_axis.jl  | 2 +-
 test/transforms/spectral_axis.jl | 9 +--------
 4 files changed, 3 insertions(+), 16 deletions(-)

diff --git a/Project.toml b/Project.toml
index 850422c..aabfe9c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -16,9 +16,6 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAstro = "6112ee07-acf9-5e0f-b108-d242c714bf9f"
 UnitfulEquivalences = "da9c4bc3-91c8-4f02-8a40-6b990d2a7e0c"
 
-[sources]
-UnitfulEquivalences = {url = "https://github.com/icweaver/UnitfulEquivalences.jl", rev = "spectral-density"}
-
 [compat]
 DustExtinction = ">=0.11.7"
 LinearAlgebra = "1.5"
@@ -27,5 +24,5 @@ PhysicalConstants = "0.2"
 RecipesBase = "1"
 Unitful = "1"
 UnitfulAstro = "1"
-UnitfulEquivalences = "0.2"
+UnitfulEquivalences = "0.2.1"
 julia = "1.10"
diff --git a/src/SpectrumBase.jl b/src/SpectrumBase.jl
index f067708..94777d6 100644
--- a/src/SpectrumBase.jl
+++ b/src/SpectrumBase.jl
@@ -10,9 +10,6 @@ export SingleSpectrum, IFUSpectrum, EchelleSpectrum, BinnedSpectrum
 export SpectrumResampler, redden, redden!, deredden, deredden!
 export redshift, redshift!, doppler_shift, doppler_shift!
 
-# Unit equivalences (re-exported from UnitfulEquivalences.jl)
-export SpectralDensity
-
 # Utilities
 export blackbody #, line_flux, equivalent_width
 
diff --git a/src/transforms/spectral_axis.jl b/src/transforms/spectral_axis.jl
index 89b6c96..f24d00f 100644
--- a/src/transforms/spectral_axis.jl
+++ b/src/transforms/spectral_axis.jl
@@ -38,7 +38,7 @@ frequency) yields an axis with the opposite direction of monotonicity.
 Given a single unit, flux values are copied unchanged. Beware that a flux density then
 remains "per" its original coordinate (e.g., a per-wavelength density against a frequency
 axis). Pass a `(spectral_unit, flux_unit)` tuple, i.e., the shape returned by `unit(spec)`, to
-also convert flux density values into the matching convention via [`UnitfulEquivalences.SpectralDensity`](https://github.com/sostock/UnitfulEquivalences.jl/pull/31),
+also convert flux density values into the matching convention via [`UnitfulEquivalences.SpectralDensity`](https://github.com/sostock/UnitfulEquivalences.jl),
 which preserves integrals: ``F_ν = F_λ λ^2 / c``.
 
 # Examples
diff --git a/test/transforms/spectral_axis.jl b/test/transforms/spectral_axis.jl
index 76f9574..2d231b0 100644
--- a/test/transforms/spectral_axis.jl
+++ b/test/transforms/spectral_axis.jl
@@ -1,4 +1,4 @@
-using SpectrumBase: SpectralDensity, redshift
+using SpectrumBase: redshift
 
 const H_JS = 6.62607015e-34u"J*s" # Planck constant, exact by SI definition
 
@@ -67,13 +67,6 @@ end
     spec = spectrum(wavelength, F_λ)
     c = C_MPS * u"m/s"
 
-    # Pointwise, SpectralDensity relates F_λ and F_ν at a given spectral coordinate,
-    # which may be given in any of the equivalent coordinates
-    @test uconvert(u"Jy", 1.0u"W/m^2/μm", SpectralDensity(2.0u"μm")) ≈
-        uconvert(u"Jy", 1.0u"W/m^2/μm" * (2.0u"μm")^2 / c)
-    @test uconvert(u"Jy", 1.0u"W/m^2/μm", SpectralDensity(uconvert(u"THz", c / 2.0u"μm"))) ≈
-        uconvert(u"Jy", 1.0u"W/m^2/μm", SpectralDensity(2.0u"μm"))
-
     converted = @inferred uconvert((u"THz", u"Jy"), spec)
     expected_F_ν = uconvert.(u"Jy", F_λ .* wavelength .^ 2 ./ c)
     @test spectral_axis(converted) ≈ uconvert.(u"THz", c ./ wavelength)

From 277d0199c6891e49bc8779532ef143f9e156aa36 Mon Sep 17 00:00:00 2001
From: Ian Weaver <weaveric@gmail.com>
Date: Wed, 16 Sep 2026 18:40:03 -0700
Subject: [PATCH 10/10] typo

---
 docs/src/index.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/index.md b/docs/src/index.md
index 1df9a08..f5bfb19 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -37,7 +37,7 @@ spec = spectrum(wave, flux)
 plot(spec)
 ```
 
-For constructing higher dimensional spectra, e.g., for echelle or IFU spectra, see the docstrings for [EchelleSpectrum](@ref) and [IFUSpectrum](@ref), respectively.
+For constructing higher dimensional spectra, e.g., for echelle or IFU spectra, see the docstrings for [`EchelleSpectrum`](@ref) and [`IFUSpectrum`](@ref), respectively.
 
 ## Citation