diff --git a/README.md b/README.md
index 5082a20..b9f421f 100644
--- a/README.md
+++ b/README.md
@@ -20,3 +20,33 @@ Currently this package can only be installed from github. To do so, either clone
     pkg> add https://github.com/JuliaAstro/Spectra.jl
 
 from the `pkg` command line (Press `]` from Julia REPL)
+
+## Developer documentation
+
+Below we show the commands to run from the package root level to develop the tests and documentation.
+
+### Tests
+
+```julia-repl
+julia --proj
+
+julia> import Pkg
+
+# List tests
+julia> Pkg.test("Spectra"; test_args = `--list`)
+
+# Run specific testsets by name. Will match with `startswith`
+julia> Pkg.test("Spectra"; test_args = `--verbose <testset name>`)
+```
+
+### Docs
+
+Assuming `LiveServer.jl` is in your global environment:
+
+```julia-repl
+julia --proj=docs/
+
+julia> using LiveServer
+
+julia> servedocs(; include_dirs = ["src/"])
+```
diff --git a/docs/Project.toml b/docs/Project.toml
index 4d4d0b2..2a6dfef 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -3,6 +3,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Spectra = "391af1a9-06f1-59d3-8d21-0be089654739"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAstro = "6112ee07-acf9-5e0f-b108-d242c714bf9f"
diff --git a/docs/make.jl b/docs/make.jl
index 950496e..e0cf339 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -2,12 +2,15 @@ using Documenter
 using Spectra
 using Unitful
 using Measurements
+using Revise
+
+Revise.revise()
 
 DocMeta.setdocmeta!(Spectra, :DocTestSetup, :(using Spectra); recursive = true)
 
 makedocs(sitename = "Spectra.jl",
     format = Documenter.HTML(;
-        prettyurls = get(ENV, "CI", nothing) == "true",
+        prettyurls = true,
         canonical = "https://juliaastro.org/Spectra/stable/",
     ),
     authors = "Miles Lucas and contributors.",
@@ -18,7 +21,7 @@ makedocs(sitename = "Spectra.jl",
         "spectrum.md",
         "transforms.md",
         "fitting.md",
-        "analysis.md",
+        "utils.md",
         "contrib.md",
     ],
     warnonly = [:missing_docs],
@@ -27,5 +30,7 @@ makedocs(sitename = "Spectra.jl",
 
 deploydocs(;
     repo = "github.com/JuliaAstro/Spectra.jl.git",
-    versions = ["stable" => "v^", "v#.#"] # Restrict to minor releases
+    devbranch = "main",
+    push_preview = true,
+    versions = ["stable" => "v^", "v#.#"], # Restrict to minor releases
 )
diff --git a/docs/src/index.md b/docs/src/index.md
index 92b6a7f..673683a 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -41,7 +41,10 @@ julia> wave = (10 .^ read(f[2], "loglam"))u"angstrom";
 julia> flux = (read(f[2], "flux") .* 1e-17)u"erg/s/cm^2/angstrom";
 
 julia> spec = spectrum(wave, flux)
-Spectrum(Quantity{Float32, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, erg, cm^-2, s^-1), 𝐌 𝐋^-1 𝐓^-3, nothing}})
+SingleSpectrum(Quantity{Float32, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, erg, cm^-2, s^-1), 𝐌 𝐋^-1 𝐓^-3, nothing}})
+  wave (3827,): 3815.0483f0 Å .. 9206.613f0 Å
+  flux (3827,): 2.182261505126953e-15 erg Å^-1 cm^-2 s^-1 .. 1.7559197998046877e-15 erg Å^-1 cm^-2 s^-1
+  meta: Dict{Symbol, Any}()
 
 julia> plot(spec);
 ```
@@ -50,9 +53,10 @@ julia> plot(spec);
 
 ```jldoctest guide
 julia> cont_fit = continuum(spec)
-Spectrum(Quantity{Float32, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, erg, cm^-2, s^-1), 𝐌 𝐋^-1 𝐓^-3, nothing}})
-  coeffs: Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, erg, cm^-2, s^-1), 𝐌 𝐋^-1 𝐓^-3, nothing}}[1.983152216046405e-15 erg Å^-1 cm^-2 s^-1, -1.8822245369267038e-16 erg Å^-1 cm^-2 s^-1, -1.0422750370065006e-16 erg Å^-1 cm^-2 s^-1, 4.8112282273206135e-17 erg Å^-1 cm^-2 s^-1]
-  normalized: true
+SingleSpectrum(Quantity{Float32, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, erg, cm^-2, s^-1), 𝐌 𝐋^-1 𝐓^-3, nothing}})
+  wave (3827,): 3815.0483f0 Å .. 9206.613f0 Å
+  flux (3827,): 1.0808438837160355 erg Å^-1 cm^-2 s^-1 .. 1.0098373106940344 erg Å^-1 cm^-2 s^-1
+  meta: Dict{Symbol, Any}(:coeffs => Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, erg, cm^-2, s^-1), 𝐌 𝐋^-1 𝐓^-3, nothing}}[1.983152216046405e-15 erg Å^-1 cm^-2 s^-1, -1.8822245369267038e-16 erg Å^-1 cm^-2 s^-1, -1.0422750370065006e-16 erg Å^-1 cm^-2 s^-1, 4.8112282273206135e-17 erg Å^-1 cm^-2 s^-1], :normalized => true)
 
 julia> plot(cont_fit, xlims=(6545, 6600));
 ```
diff --git a/docs/src/spectrum.md b/docs/src/spectrum.md
index 39995dd..119b217 100644
--- a/docs/src/spectrum.md
+++ b/docs/src/spectrum.md
@@ -1,10 +1,3 @@
-```@meta
-DocTestSetup = quote
-  using Spectra, Random
-  Random.seed!(11894)
-end
-```
-
 # Spectrum
 
 Here we will go over the different spectral types and how we use them.
@@ -14,7 +7,11 @@ Here we will go over the different spectral types and how we use them.
 Spectra are defined as possible subtypes of `AbstractSpectrum`. You can use these directly for construction, or use the catch-all [`spectrum`](@ref) function, which is preferred.
 
 ```@docs
+Spectra.AbstractSpectrum
 Spectra.Spectrum
+Spectra.SingleSpectrum
+Spectra.EchelleSpectrum
+Spectra.IFUSpectrum
 ```
 
 ## Constructors
@@ -28,25 +25,37 @@ Spectra.spectrum
 
 For more advanced transformations, see [Transformations](@ref)
 
+### Getters
+```@docs
+Spectra.wave(::AbstractSpectrum)
+Spectra.flux(::AbstractSpectrum)
+Spectra.meta(::AbstractSpectrum)
+```
+
+### Array interface
+
 | Function                           |
 |:-----------------------------------|
+| `Base.argmax(::AbstractSpectrum)`  |
+| `Base.argmin(::AbstractSpectrum)`  |
+| `Base.eltype(::AbstractSpectrum)`  |
+| `Base.findmax(::AbstractSpectrum)` |
+| `Base.findmin(::AbstractSpectrum)` |
+| `Base.iterate(::AbstractSpectrum)` |
 | `Base.length(::AbstractSpectrum)`  |
-| `Base.size(::AbstractSpectrum)`    |
 | `Base.maximum(::AbstractSpectrum)` |
 | `Base.minimum(::AbstractSpectrum)` |
-| `Base.argmax(::AbstractSpectrum)`  |
-| `Base.argmin(::AbstractSpectrum)`  |
-| `Base.findmax(::AbstractSpectrum)`  |
-| `Base.findmin(::AbstractSpectrum)`  |
+| `Base.size(::AbstractSpectrum)`    |
 
 ### Arithmetic
 
-| Function                           |
-|:-----------------------------------|
-| `+(::AbstractSpectrum, A)`         |
-| `-(::AbstractSpectrum, A)`         |
-| `*(::AbstractSpectrum, A)`         |
-| `/(::AbstractSpectrum, A)`         |
+| Function                                            |
+|:----------------------------------------------------|
+| `+(::AbstractSpectrum, A)`                          |
+| `-(::AbstractSpectrum, A)`                          |
+| `*(::AbstractSpectrum, A)`                          |
+| `/(::AbstractSpectrum, A)`                          |
+| `Base.(==)(::AbstractSpectrum, ::AbstractSpectrum)` |
 
 ## Unitful helpers
 
@@ -76,7 +85,3 @@ savefig("spec-plot.svg"); nothing # hide
 ```@index
 Pages = ["spectrum.md"]
 ```
-
-```@meta
-DocTestSetup = nothing
-```
diff --git a/docs/src/transforms.md b/docs/src/transforms.md
index a5a05c5..0f25a5a 100644
--- a/docs/src/transforms.md
+++ b/docs/src/transforms.md
@@ -1,10 +1,3 @@
-```@meta
-DocTestSetup = quote
-  using Spectra, Random
-  Random.seed!(11894)
-end
-```
-
 # Transformations
 
 ## Extinction
@@ -12,20 +5,13 @@ end
 By levaraging [DustExtinction.jl](https://github.com/juliaastro/dustextinction.jl) we can apply common reddening laws to our spectra.
 
 ```jldoctest
-julia> using Unitful, Measurements, Random
+julia> using Spectra, Unitful, Measurements, Random
 
 julia> rng = Random.seed!(0);
 
-julia> wave = (1:0.5:3)u"μm"
-(1.0:0.5:3.0) μm
+julia> wave = (1:0.5:3)u"μm";
 
-julia> sigma = randn(rng, size(wave))
-5-element Vector{Float64}:
-  0.942970533446119
-  0.13392275765318448
-  1.5250689085124804
-  0.12390123120559722
- -1.205772284259936
+julia> sigma = randn(rng, size(wave));
 
 julia> flux = (100 .± sigma)u"W/m^2/μm"
 5-element Vector{Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(μm^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}}:
@@ -36,10 +22,16 @@ julia> flux = (100 .± sigma)u"W/m^2/μm"
  100.0 ± -1.2 W μm^-1 m^-2
 
 julia> spec = spectrum(wave, flux)
-Spectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(μm,), 𝐋, nothing}}, Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(μm^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}})
+SingleSpectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(μm,), 𝐋, nothing}}, Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(μm^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}})
+  wave (5,): 1.0 μm .. 3.0 μm
+  flux (5,): 100.0 ± 0.94 W μm^-1 m^-2 .. 100.0 ± -1.2 W μm^-1 m^-2
+  meta: Dict{Symbol, Any}()
 
 julia> red = redden(spec, 0.3)
-Spectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(μm,), 𝐋, nothing}}, Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(μm^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}})
+SingleSpectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(μm,), 𝐋, nothing}}, Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(μm^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}})
+  wave (5,): 1.0 μm .. 3.0 μm
+  flux (5,): 89.44 ± 0.84 W μm^-1 m^-2 .. 98.1 ± 1.2 W μm^-1 m^-2
+  meta: Dict{Symbol, Any}()
 
 julia> red.flux
 5-element Vector{Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(μm^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}}:
@@ -50,7 +42,10 @@ julia> red.flux
   98.1 ± 1.2 W μm^-1 m^-2
 
 julia> deredden!(red, 0.3)
-Spectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(μm,), 𝐋, nothing}}, Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(μm^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}})
+SingleSpectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(μm,), 𝐋, nothing}}, Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(μm^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}})
+  wave (5,): 1.0 μm .. 3.0 μm
+  flux (5,): 100.0 ± 0.94 W μm^-1 m^-2 .. 100.0 ± 1.2 W μm^-1 m^-2
+  meta: Dict{Symbol, Any}()
 
 julia> red.flux ≈ spec.flux
 true
@@ -64,7 +59,3 @@ redden!
 deredden
 deredden!
 ```
-
-```@meta
-DocTestSetup = nothing
-```
diff --git a/docs/src/analysis.md b/docs/src/utils.md
similarity index 60%
rename from docs/src/analysis.md
rename to docs/src/utils.md
index 6592083..8fb1048 100644
--- a/docs/src/analysis.md
+++ b/docs/src/utils.md
@@ -1,6 +1,7 @@
-# Analysis
+# Utilities
 
 ```@docs
+Spectra.blackbody
 Spectra.equivalent_width
 Spectra.line_flux
-```
\ No newline at end of file
+```
diff --git a/src/EchelleSpectrum.jl b/src/EchelleSpectrum.jl
deleted file mode 100644
index e14ba1e..0000000
--- a/src/EchelleSpectrum.jl
+++ /dev/null
@@ -1,31 +0,0 @@
-struct EchelleSpectrum{W <: Number,F <: Number} <: AbstractSpectrum{W,F}
-    wave::Matrix{W}
-    flux::Matrix{F}
-    meta::Dict{Symbol,Any}
-end
-
-EchelleSpectrum(wave, flux, meta::Dict{Symbol,Any}) = EchelleSpectrum(collect(wave), collect(flux), meta)
-
-function Base.show(io::IO, spec::EchelleSpectrum)
-    println(io, "EchelleSpectrum($(eltype(spec.wave)), $(eltype(spec.flux)))")
-    print(io, "  # orders: $(size(spec, 1))")
-    for (key, val) in spec.meta
-        print(io, "\n  $key: $val")
-    end
-end
-
-function Base.getindex(spec::EchelleSpectrum, i::Integer)
-    wave = spec.wave[i, :]
-    flux = spec.flux[i, :]
-    meta = merge(Dict(:Order => i), spec.meta)
-    return Spectrum(wave, flux, meta)
-end
-
-function Base.getindex(spec::EchelleSpectrum, I::AbstractVector)
-    waves = spec.wave[I, :]
-    flux = spec.flux[I, :]
-    meta = merge(Dict(:Orders => (first(I), last(I))), spec.meta)
-    return EchelleSpectrum(waves, flux, meta)
-end
-
-Base.lastindex(spec::EchelleSpectrum) = lastindex(spec.flux, 1)
diff --git a/src/Spectra.jl b/src/Spectra.jl
index f804c52..bc7665b 100644
--- a/src/Spectra.jl
+++ b/src/Spectra.jl
@@ -1,9 +1,13 @@
 module Spectra
 
-# common.jl
-export AbstractSpectrum, spectrum
+# Uniform API
+export AbstractSpectrum, Spectrum, spectrum
+# spectra_single.jl, spectra_ifu.jl, spectra_echelle.jl
+export SingleSpectrum, IFUSpectrum, EchelleSpectrum
 # utils.jl
-export blackbody
+export blackbody, line_flux, equivalent_width
+# fitting/fitting.jl
+export continuum, continuum!
 # transforms/redden.jl
 export redden, redden!, deredden, deredden!
 
@@ -12,12 +16,191 @@ using Measurements: Measurements, Measurement
 using Unitful: Unitful, Quantity, @u_str, ustrip, unit, dimension
 using PhysicalConstants.CODATA2018: h, c_0, k_B
 
-# AbstractSpectrum and common functionality
-include("common.jl")
+"""
+    AbstractSpectrum{W<:Number, F<:Number}
+
+An abstract holder for astronomical spectra. All types inheriting from this must have the following fields:
+
+* `wave::Array{W, M}`
+* `flux::Array{F, N}`
+* `meta::Dict{Symbol, Any}`
+
+See [`SingleSpectrum`](@ref), [`EchelleSpectrum`](@ref), and [`IFUSpectrum`](@ref) for different subtypes.
+"""
+abstract type AbstractSpectrum{W,F} end
+
+"""
+    Spectrum <: AbstractSpectrum
+
+A spectrum or spectra stored as arrays of real numbers. The wavelengths are assumed to be in angstrom.
+"""
+mutable struct Spectrum{W<:Number, F<:Number, M, N} <: AbstractSpectrum{W, F}
+    wave::AbstractArray{W, M}
+    flux::AbstractArray{F, N}
+    meta::Dict{Symbol,Any}
+    function Spectrum{W, F, M, N}(wave, flux, meta) where {W<:Number, F<:Number, M, N}
+        # Dimension compatibility check
+        size(wave, 1) != size(flux, 1) && throw(ArgumentError(
+        """
+        Wavelength and flux sizes are incompatible. Currently supported sizes are:
+
+        * SingleSpectrum: wave (M-length vector), flux (M-length vector)
+        * EchelleSpectrum: wave (M x N matrix), flux (M x N matrix)
+        * IFUSpectrum: wave (M-length vector), flux (M x N x K matrix)
+
+        See the documentation for each spectrum type for more.
+        """))
+
+        # Wavelength monoticity check
+        w = eachcol(wave)
+        !(
+            all(issorted, w) ||
+            all(x -> issorted(x; rev=true), w)
+        ) && throw(ArgumentError(
+        "Wavelengths must be strictly increasing or decreasing."
+        ))
+
+        return new{W, F, M, N}(wave, flux, meta)
+    end
+end
+
+function Spectrum(wave, flux, meta)
+    Spectrum{eltype(wave), eltype(flux), ndims(wave), ndims(flux)}(wave, flux, meta)
+end
+
+# Doesn't seem to be used atp
+#Spectrum(wave, flux, meta::Dict{Symbol, Any}) = Spectrum(collect(wave), collect(flux), meta)
+
+"""
+    wave(::AbstractSpectrum)
+
+Return the wavelengths of the spectrum.
+"""
+wave(spec::AbstractSpectrum) = spec.wave
+
+"""
+    flux(::AbstractSpectrum)
+
+Return the flux of the spectrum.
+"""
+flux(spec::AbstractSpectrum) = spec.flux
+
+"""
+    meta(::AbstractSpectrum)
+
+Return the meta of the spectrum.
+"""
+meta(spec::AbstractSpectrum) = spec.meta
+
+function Base.getproperty(spec::AbstractSpectrum, nm::Symbol)
+    if nm in keys(getfield(spec, :meta))
+        return getfield(spec, :meta)[nm]
+    else
+        return getfield(spec, nm)
+    end
+end
+
+function Base.propertynames(spec::AbstractSpectrum)
+    natural = (:wave, :flux, :meta)
+    m = keys(meta(spec))
+    return (natural..., m...)
+end
+
+# Collection
+Base.argmax(spec::AbstractSpectrum) = argmax(flux(spec))
+Base.argmin(spec::AbstractSpectrum) = argmin(flux(spec))
+Base.eltype(spec::AbstractSpectrum) = eltype(flux(spec))
+Base.findmax(spec::AbstractSpectrum) = findmax(flux(spec))
+Base.findmin(spec::AbstractSpectrum) = findmin(flux(spec))
+Base.length(spec::AbstractSpectrum) = length(flux(spec))
+Base.maximum(spec::AbstractSpectrum) = maximum(flux(spec))
+Base.minimum(spec::AbstractSpectrum) = minimum(flux(spec))
+Base.size(spec::AbstractSpectrum) = size(flux(spec))
+Base.size(spec::AbstractSpectrum, i) = size(flux(spec), i)
+function Base.iterate(spec::AbstractSpectrum, state=0)
+    state == length(spec) && return nothing
+    return spec[begin + state], state + 1
+end
+
+# Arithmetic
+Base.:(==)(s::AbstractSpectrum, o::AbstractSpectrum) = wave(s) == wave(o) && flux(s) == flux(o) && meta(s) == meta(o)
+Base.:+(s::T, A) where {T <: AbstractSpectrum} = T(wave(s), flux(s) .+ A, meta(s))
+Base.:*(s::T, A::Union{Real, AbstractVector}) where {T <: AbstractSpectrum} = T(wave(s), flux(s) .* A, meta(s))
+Base.:/(s::T, A) where {T <: AbstractSpectrum} = T(wave(s), flux(s) ./ A, meta(s))
+Base.:-(s::T) where {T <: AbstractSpectrum} = T(wave(s), -flux(s), meta(s))
+Base.:-(s::AbstractSpectrum, A) = s + -A
+Base.:-(A, s::AbstractSpectrum) = s - A
+Base.:-(s::AbstractSpectrum, o::AbstractSpectrum) = s - o # Satisfy Aqua
+
+# Multi-Spectrum
+Base.:+(s::T, o::T) where {T <: AbstractSpectrum} = T(wave(s), flux(s) .+ flux(o), meta(s))
+Base.:*(s::T, o::T) where {T <: AbstractSpectrum} = T(wave(s), flux(s) .* flux(o), meta(s))
+Base.:/(s::T, o::T) where {T <: AbstractSpectrum} = T(wave(s), flux(s) ./ flux(o) * unit(s)[2], meta(s))
+Base.:-(s::T, o::T) where {T <: AbstractSpectrum} = T(wave(s), flux(s) .- flux(o), meta(s))
+
+"""
+    Unitful.ustrip(::AbstractSpectrum)
+
+Remove the units from a spectrum. Useful for processing spectra in tools that don't play nicely with `Unitful.jl`
+
+# Examples
+```jldoctest
+julia> using Random
+
+julia> rng = Random.seed!(0)
+TaskLocalRNG()
+
+julia> using Unitful, UnitfulAstro
+
+julia> wave = range(1e4, 3e4, length=1000);
+
+julia> flux = wave .* 10 .+ randn(rng, 1000);
+
+julia> spec = spectrum(wave*u"angstrom", flux*u"W/m^2/angstrom")
+SingleSpectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}})
+  wave (1000,): 10000.0 Å .. 30000.0 Å
+  flux (1000,): 99999.76809093042 W Å^-1 m^-2 .. 300000.2474309158 W Å^-1 m^-2
+  meta: Dict{Symbol, Any}()
+
+julia> ustrip(spec)
+SingleSpectrum(Float64, Float64)
+  wave (1000,): 10000.0 .. 30000.0
+  flux (1000,): 99999.76809093042 .. 300000.2474309158
+  meta: Dict{Symbol, Any}()
+```
+"""
+Unitful.ustrip(spec::AbstractSpectrum) = spectrum(ustrip.(wave(spec)), ustrip.(flux(spec)); meta(spec)...)
+
+"""
+    Unitful.unit(::AbstractSpectrum)
+
+Get the units of a spectrum. Returns a tuple of the wavelength units and flux/sigma units
+
+# Examples
+```jldoctest
+julia> using Random
+
+julia> rng = Random.seed!(0)
+TaskLocalRNG()
+
+julia> using Unitful, UnitfulAstro
+
+julia> wave = range(1e4, 3e4, length=1000);
+
+julia> flux = wave .* 10 .+ randn(rng, 1000);
+
+julia> spec = spectrum(wave * u"angstrom", flux * u"W/m^2/angstrom");
+
+julia> w_unit, f_unit = unit(spec)
+(Å, W Å^-1 m^-2)
+```
+"""
+Unitful.unit(spec::AbstractSpectrum) = unit(eltype(wave(spec))), unit(eltype(flux(spec)))
 
 # Spectrum types and basic arithmetic
-include("spectrum.jl")
-include("EchelleSpectrum.jl")
+include("spectrum_single.jl")
+include("spectrum_ifu.jl")
+include("spectrum_echelle.jl")
 
 """
     spectrum(wave, flux; kwds...)
@@ -31,11 +214,16 @@ julia> wave = range(1e4, 4e4, length=1000);
 julia> flux = 100 .* ones(size(wave));
 
 julia> spec = spectrum(wave, flux)
-Spectrum(Float64, Float64)
+SingleSpectrum(Float64, Float64)
+  wave (1000,): 10000.0 .. 40000.0
+  flux (1000,): 100.0 .. 100.0
+  meta: Dict{Symbol, Any}()
 
 julia> spec = spectrum(wave, flux, name="Just Noise")
-Spectrum(Float64, Float64)
-  name: Just Noise
+SingleSpectrum(Float64, Float64)
+  wave (1000,): 10000.0 .. 40000.0
+  flux (1000,): 100.0 .. 100.0
+  meta: Dict{Symbol, Any}(:name => "Just Noise")
 
 julia> spec.name
 "Just Noise"
@@ -45,50 +233,61 @@ There is easy integration with [Unitful.jl](https://github.com/JuliaPhysics/Unit
 and its sub-projects and [Measurements.jl](https://github.com/juliaphysics/measurements.jl)
 
 ```jldoctest
+julia> using Random
+
+julia> rng = Random.seed!(0)
+TaskLocalRNG()
+
 julia> using Unitful, UnitfulAstro, Measurements
 
 julia> wave = range(1, 4, length=1000)u"μm";
 
-julia> sigma = randn(size(wave));
+julia> sigma = randn(rng, size(wave));
 
 julia> flux = (100 .± sigma)u"erg/cm^2/s/angstrom";
 
 julia> spec = spectrum(wave, flux)
-Spectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(μm,), 𝐋, nothing}}, Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, erg, cm^-2, s^-1), 𝐌 𝐋^-1 𝐓^-3, nothing}})
+SingleSpectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(μm,), 𝐋, nothing}}, Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, erg, cm^-2, s^-1), 𝐌 𝐋^-1 𝐓^-3, nothing}})
+  wave (1000,): 1.0 μm .. 4.0 μm
+  flux (1000,): 100.0 ± -0.23 erg Å^-1 cm^-2 s^-1 .. 100.0 ± 0.25 erg Å^-1 cm^-2 s^-1
+  meta: Dict{Symbol, Any}()
 ```
 
-For a multi-order spectrum, all orders must have the same length, so be sure to pad any ragged orders with NaN.
+For an echelle spectrum, all orders must have the same length, so be sure to pad any ragged orders with NaN.
 
 ```jldoctest
-julia> wave = reshape(range(100, 1e4, length=1000), 100, 10)';
+julia> wave = reshape(range(100, 1e4, length=1000), 100, 10);
 
-julia> flux = ones(10, 100) .* collect(1:10);
+julia> flux = repeat(1:10.0, 1, 100)';
 
 julia> spec = spectrum(wave, flux)
 EchelleSpectrum(Float64, Float64)
   # orders: 10
+  wave (100, 10): 100.0 .. 10000.0
+  flux (100, 10): 1.0 .. 10.0
+  meta: Dict{Symbol, Any}()
 ```
 """
 function spectrum(wave::AbstractVector{<:Real}, flux::AbstractVector{<:Real}; kwds...)
-    @assert size(wave) == size(flux) "wave and flux must have equal size"
     Spectrum(wave, flux, Dict{Symbol,Any}(kwds))
 end
 
-function spectrum(wave::AbstractVector{<:Quantity}, flux::AbstractVector{<:Quantity}; kwds...)
-    @assert size(wave) == size(flux) "wave and flux must have equal size"
-    @assert dimension(eltype(wave)) == u"𝐋" "wave not recognized as having dimensions of wavelengths"
+function spectrum(wave::AbstractVector{<:Real}, flux::AbstractArray{<:Real, 3}; kwds...)
     Spectrum(wave, flux, Dict{Symbol,Any}(kwds))
 end
 
 function spectrum(wave::AbstractMatrix{<:Real}, flux::AbstractMatrix{<:Real}; kwds...)
-    @assert size(wave) == size(flux) "wave and flux must have equal size"
-    EchelleSpectrum(wave, flux, Dict{Symbol,Any}(kwds))
+    Spectrum(wave, flux, Dict{Symbol,Any}(kwds))
+end
+
+function spectrum(wave::AbstractVector{<:Quantity}, flux::AbstractVector{<:Quantity}; kwds...)
+    @assert dimension(eltype(wave)) == u"𝐋" "wave not recognized as having dimensions of wavelengths"
+    Spectrum(wave, flux, Dict{Symbol,Any}(kwds))
 end
 
 function spectrum(wave::AbstractMatrix{<:Quantity}, flux::AbstractMatrix{<:Quantity}; kwds...)
-    @assert size(wave) == size(flux) "wave and flux must have equal size"
     @assert dimension(eltype(wave)) == u"𝐋" "wave not recognized as having dimensions of wavelengths"
-    EchelleSpectrum(wave, flux, Dict{Symbol,Any}(kwds))
+    Spectrum(wave, flux, Dict{Symbol,Any}(kwds))
 end
 
 # tools
diff --git a/src/common.jl b/src/common.jl
deleted file mode 100644
index 5c0d20c..0000000
--- a/src/common.jl
+++ /dev/null
@@ -1,108 +0,0 @@
-"""
-    AbstractSpectrum{W<:Number, F<:Number}
-
-An abstract holder for astronomical spectra. All types inheriting from this must have the following fields
-
-- wave::Array{W, N}
-- flux::Array{F, N}
-- meta::Dict{Symbol, Any}
-"""
-abstract type AbstractSpectrum{W,F} end
-
-function Base.getproperty(spec::AbstractSpectrum, nm::Symbol)
-    if nm in keys(getfield(spec, :meta))
-        return getfield(spec, :meta)[nm]
-    else
-        return getfield(spec, nm)
-    end
-end
-
-function Base.propertynames(spec::AbstractSpectrum)
-    natural = (:wave, :flux, :meta)
-    meta = keys(spec.meta)
-    return (natural..., meta...)
-end
-
-"""
-    wave(::AbstractSpectrum)
-
-Return the wavelengths of the spectrum.
-"""
-wave(spec::AbstractSpectrum) = spec.wave
-
-"""
-    flux(::AbstractSpectrum)
-
-Return the flux of the spectrum.
-"""
-flux(spec::AbstractSpectrum) = spec.flux
-
-# Collection
-Base.eltype(spec::AbstractSpectrum) = eltype(spec.flux)
-Base.size(spec::AbstractSpectrum) = size(spec.flux)
-Base.size(spec::AbstractSpectrum, i) = size(spec.flux, i)
-Base.length(spec::AbstractSpectrum) = length(spec.flux)
-Base.maximum(spec::AbstractSpectrum) = maximum(spec.flux)
-Base.minimum(spec::AbstractSpectrum) = minimum(spec.flux)
-Base.argmax(spec::AbstractSpectrum) = argmax(spec.flux)
-Base.argmin(spec::AbstractSpectrum) = argmin(spec.flux)
-Base.findmax(spec::AbstractSpectrum) = findmax(spec.flux)
-Base.findmin(spec::AbstractSpectrum) = findmin(spec.flux)
-
-# Arithmetic
-Base.:+(s::T, A) where {T <: AbstractSpectrum} = T(s.wave, s.flux .+ A, s.meta)
-Base.:*(s::T, A::Union{Real, AbstractVector}) where {T <: AbstractSpectrum} = T(s.wave, s.flux .* A, s.meta)
-Base.:/(s::T, A) where {T <: AbstractSpectrum} = T(s.wave, s.flux ./ A, s.meta)
-Base.:-(s::T) where {T <: AbstractSpectrum} = T(s.wave, -s.flux, s.meta)
-Base.:-(s::AbstractSpectrum, A) = s + -A
-Base.:-(A, s::AbstractSpectrum) = s - A
-Base.:-(s::AbstractSpectrum, o::AbstractSpectrum) = s - o # Satisfy Aqua
-
-# Multi-Spectrum
-Base.:+(s::T, o::T) where {T <: AbstractSpectrum} = T(s.wave, s.flux .+ o.flux, s.meta)
-Base.:*(s::T, o::T) where {T <: AbstractSpectrum} = T(s.wave, s.flux .* o.flux, s.meta)
-Base.:/(s::T, o::T) where {T <: AbstractSpectrum} = T(s.wave, s.flux ./ o.flux * unit(s)[2], s.meta)
-Base.:-(s::T, o::T) where {T <: AbstractSpectrum} = T(s.wave, s.flux .- o.flux, s.meta)
-
-"""
-    Unitful.ustrip(::AbstractSpectrum)
-
-Remove the units from a spectrum. Useful for processing spectra in tools that don't play nicely with `Unitful.jl`
-
-# Examples
-```jldoctest
-julia> using Unitful, UnitfulAstro
-
-julia> wave = range(1e4, 3e4, length=1000);
-
-julia> flux = wave .* 10 .+ randn(1000);
-
-julia> spec = spectrum(wave*u"angstrom", flux*u"W/m^2/angstrom")
-Spectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}})
-
-julia> ustrip(spec)
-Spectrum(Float64, Float64)
-```
-"""
-Unitful.ustrip(spec::AbstractSpectrum) = spectrum(ustrip.(spec.wave), ustrip.(spec.flux); spec.meta...)
-
-"""
-    Unitful.unit(::AbstractSpectrum)
-
-Get the units of a spectrum. Returns a tuple of the wavelength units and flux/sigma units
-
-# Examples
-```jldoctest
-julia> using Unitful, UnitfulAstro
-
-julia> wave = range(1e4, 3e4, length=1000);
-
-julia> flux = wave .* 10 .+ randn(1000);
-
-julia> spec = spectrum(wave * u"angstrom", flux * u"W/m^2/angstrom");
-
-julia> w_unit, f_unit = unit(spec)
-(Å, W Å^-1 m^-2)
-```
-"""
-Unitful.unit(spec::AbstractSpectrum) = unit(eltype(spec.wave)), unit(eltype(spec.flux))
diff --git a/src/fitting/fitting.jl b/src/fitting/fitting.jl
index d2e3343..37f0154 100644
--- a/src/fitting/fitting.jl
+++ b/src/fitting/fitting.jl
@@ -1,6 +1,5 @@
 using LinearAlgebra: /, \, diagm, pinv
 
-export continuum, line_flux, equivalent_width
 
 function chebvander(x::AbstractVector{T}, deg::Int) where {T <: Number}
     v = Matrix{T}(undef, length(x), deg + 1)
@@ -45,24 +44,3 @@ end
 Return a continuum-normalized spectrum by fitting the continuum with a Chebyshev polynomial of degree `deg`.
 """
 continuum(spec::AbstractSpectrum, deg::Int = 3) = continuum!(deepcopy(spec), deg)
-
-"""
-    equivalent_width(::AbstractSpectrum)
-
-Calculate the equivalent width of the given continuum-normalized spectrum. Return value has units equal to wavelengths.
-"""
-function equivalent_width(spec::AbstractSpectrum)
-    dx = spec.wave[end] - spec.wave[1]
-    flux = ustrip(line_flux(spec))
-    return dx - flux * unit(dx)
-end
-
-"""
-    line_flux(::AbstractSpectrum)
-
-Calculate the line flux of the given continuum-normalized spectrum. Return value has units equal to flux.
-"""
-function line_flux(spec::AbstractSpectrum)
-    avg_dx = diff(spec.wave)
-    return sum(spec.flux[2:end] .* avg_dx)
-end
diff --git a/src/plotting.jl b/src/plotting.jl
index 3515c2d..427a704 100644
--- a/src/plotting.jl
+++ b/src/plotting.jl
@@ -1,4 +1,4 @@
-@recipe function f(spec::Spectrum)
+@recipe function f(spec::SingleSpectrum)
     seriestype --> :step
     xlabel --> "wave"
     ylabel --> "flux density"
diff --git a/src/spectrum.jl b/src/spectrum.jl
deleted file mode 100644
index c1a20f9..0000000
--- a/src/spectrum.jl
+++ /dev/null
@@ -1,19 +0,0 @@
-"""
-    Spectrum <: AbstractSpectrum
-
-A 1-dimensional spectrum stored as vectors of real numbers. The wavelengths are assumed to be in angstrom.
-"""
-struct Spectrum{W <: Number,F <: Number} <: AbstractSpectrum{W,F}
-    wave::Vector{W}
-    flux::Vector{F}
-    meta::Dict{Symbol,Any}
-end
-
-Spectrum(wave, flux, meta::Dict{Symbol,Any}) = Spectrum(collect(wave), collect(flux), meta)
-
-function Base.show(io::IO, spec::Spectrum)
-    print(io, "Spectrum($(eltype(spec.wave)), $(eltype(spec.flux)))")
-    for (key, val) in spec.meta
-        print(io, "\n  $key: $val")
-    end
-end
diff --git a/src/spectrum_echelle.jl b/src/spectrum_echelle.jl
new file mode 100644
index 0000000..91b0e47
--- /dev/null
+++ b/src/spectrum_echelle.jl
@@ -0,0 +1,80 @@
+"""
+    EchelleSpectrum <: AbstractSpectrum
+
+An instance of [`Spectrum`](@ref) where the wavelength and flux are both 2D arrays, i.e., ``M = N = 2``.
+
+The wavelength and flux matrices are both ``m`` rows in wavelength by ``n`` columns in [echelle order](https://en.wikipedia.org/wiki/Echelle_grating).
+
+# Examples
+
+```jldoctest
+julia> wave = reshape(1:40, 10, 4)
+10×4 reshape(::UnitRange{Int64}, 10, 4) with eltype Int64:
+  1  11  21  31
+  2  12  22  32
+  3  13  23  33
+  4  14  24  34
+  5  15  25  35
+  6  16  26  36
+  7  17  27  37
+  8  18  28  38
+  9  19  29  39
+ 10  20  30  40
+
+julia> flux = repeat(1:4, 1, 10)'
+10×4 adjoint(::Matrix{Int64}) with eltype Int64:
+ 1  2  3  4
+ 1  2  3  4
+ 1  2  3  4
+ 1  2  3  4
+ 1  2  3  4
+ 1  2  3  4
+ 1  2  3  4
+ 1  2  3  4
+ 1  2  3  4
+ 1  2  3  4
+
+julia> spec = Spectrum(wave, flux, Dict())
+EchelleSpectrum(Int64, Int64)
+  # orders: 4
+  wave (10, 4): 1 .. 40
+  flux (10, 4): 1 .. 4
+  meta: Dict{Symbol, Any}()
+
+julia> spec[1] # Indexing returns a `SingleSpectrum`
+SingleSpectrum(Int64, Int64)
+  wave (10,): 1 .. 10
+  flux (10,): 1 .. 1
+  meta: Dict{Symbol, Any}(:Order => 1)
+```
+
+See [`SingleSpectrum`](@ref) for a 1D variant, and [`IFUSpectrum`](@ref) for a 3D variant.
+"""
+const EchelleSpectrum = Spectrum{W, F, 2, 2} where {W, F}
+
+function Base.getindex(spec::EchelleSpectrum, i::Int)
+    w = wave(spec)[:, i]
+    f = flux(spec)[:, i]
+    m = merge(Dict(:Order => i), meta(spec))
+    return Spectrum(w, f, m)
+end
+
+function Base.getindex(spec::EchelleSpectrum, I::AbstractVector)
+    w = wave(spec)[:, I]
+    f = flux(spec)[:, I]
+    m = merge(Dict(:Orders => (first(I), last(I))), meta(spec))
+    return Spectrum(w, f, m)
+end
+
+Base.firstindex(spec::EchelleSpectrum) = firstindex(flux(spec), 1)
+Base.lastindex(spec::EchelleSpectrum) = lastindex(flux(spec), 1)
+
+function Base.show(io::IO, spec::EchelleSpectrum)
+    w = wave(spec)
+    f = flux(spec)
+    println(io, "EchelleSpectrum($(eltype(wave(spec))), $(eltype(flux(spec))))")
+    println(io, "  # orders: $(size(spec, 2))")
+    println(io, "  wave $(size(w)): ", first(w), " .. ", last(w))
+    println(io, "  flux $(size(f)): ", first(f), " .. ", last(f))
+    print(io, "  meta: ", meta(spec))
+end
diff --git a/src/spectrum_ifu.jl b/src/spectrum_ifu.jl
new file mode 100644
index 0000000..c1b90ee
--- /dev/null
+++ b/src/spectrum_ifu.jl
@@ -0,0 +1,81 @@
+"""
+    IFUSpectrum <: AbstractSpectrum
+
+An instance of [`Spectrum`](@ref) where the wavelength is a 1D array and flux is a 3D array, i.e., ``M = 1`` and ``N = 3``.
+
+The wavelength vector is of length ``m`` and flux 3D array has shape ``(m \\times n \\times k)`` where each entry of the flux along the wavelength axis is an ``n \\times k`` matrix.
+
+# Examples
+
+```jldoctest
+julia> using Random
+
+julia> rng = Random.seed!(0)
+TaskLocalRNG()
+
+julia> wave, flux = [20, 40, 120, 160, 200], rand(rng, 5, 10, 6);
+
+julia> spec = Spectrum(wave, flux, Dict())
+IFUSpectrum(Int64, Float64)
+  wave (5,): 20 .. 200
+  flux (5, 10, 6): 0.4552384158732863 .. 0.11698905483599475
+  meta: Dict{Symbol, Any}()
+
+julia> spec[begin] # IFU image at first wavelength
+10×6 Matrix{Float64}:
+ 0.455238   0.828104   0.735106   0.042069  0.554894  0.0715802
+ 0.746943   0.149248   0.864755   0.116243  0.519913  0.310438
+ 0.0997382  0.523732   0.315933   0.935547  0.274589  0.250664
+ 0.470257   0.654557   0.351769   0.812597  0.158201  0.617466
+ 0.678779   0.312182   0.0568161  0.622296  0.61899   0.191777
+ 0.385452   0.345902   0.448835   0.041962  0.458694  0.791756
+ 0.908402   0.609104   0.108874   0.430905  0.91365   0.430885
+ 0.0256413  0.0831649  0.179467   0.799997  0.982336  0.721449
+ 0.408092   0.361884   0.849442   0.527004  0.341892  0.499461
+ 0.239929   0.3754     0.247219   0.92438   0.733984  0.432918
+
+julia> spec[begin:3] # IFU spectrum at first three wavelengths
+IFUSpectrum(Int64, Float64)
+  wave (3,): 20 .. 120
+  flux (3, 10, 6): 0.4552384158732863 .. 0.35149138733595564
+  meta: Dict{Symbol, Any}()
+
+julia> spec[:, begin, begin] # 1D spectrum at spaxel (1, 1)
+SingleSpectrum(Int64, Float64)
+  wave (5,): 20 .. 200
+  flux (5,): 0.4552384158732863 .. 0.02964765308691042
+  meta: Dict{Symbol, Any}()
+```
+
+See [`SingleSpectrum`](@ref) for a 1D variant, and [`EchelleSpectrum`](@ref) for a 2D variant.
+"""
+const IFUSpectrum = Spectrum{W, F, 1, 3} where {W, F}
+
+Base.getindex(spec::IFUSpectrum, i) = flux(spec)[i, :, :]
+
+function Base.getindex(spec::IFUSpectrum, I::AbstractVector)
+    w = wave(spec)[I]
+    f = flux(spec)[I, :, :]
+    return Spectrum(w, f, meta(spec))
+end
+
+Base.getindex(spec::IFUSpectrum, i::Int, j, k) = flux(spec)[i, j, k]
+
+function Base.getindex(spec::IFUSpectrum, i, j, k)
+    w = wave(spec)[i]
+    f = flux(spec)[i, j, k]
+    return Spectrum(w, f, meta(spec))
+end
+
+Base.firstindex(spec::IFUSpectrum) = firstindex(flux(spec))
+Base.firstindex(spec::IFUSpectrum, i) = firstindex(flux(spec), i)
+Base.lastindex(spec::IFUSpectrum, i) = lastindex(flux(spec), i)
+
+function Base.show(io::IO, spec::IFUSpectrum)
+    w = wave(spec)
+    f = flux(spec)
+    println(io, "IFUSpectrum($(eltype(w)), $(eltype(f)))")
+    println(io, "  wave $(size(w)): ", first(w), " .. ", last(w))
+    println(io, "  flux $(size(f)): ", first(f), " .. ", last(f))
+    print(io, "  meta: ", meta(spec))
+end
diff --git a/src/spectrum_single.jl b/src/spectrum_single.jl
new file mode 100644
index 0000000..27ca085
--- /dev/null
+++ b/src/spectrum_single.jl
@@ -0,0 +1,43 @@
+"""
+    SingleSpectrum <: AbstractSpectrum
+
+An instance of [`Spectrum`](@ref) where the wavelength and flux are both 1D arrays, i.e., ``M = N = 1``.
+
+Both wavelength and flux vectors must have the same length, ``m = n``, respectively.
+
+# Examples
+
+```jldoctest
+julia> spec = Spectrum([6, 7, 8, 9], [2, 3, 4, 5], Dict())
+SingleSpectrum(Int64, Int64)
+  wave (4,): 6 .. 9
+  flux (4,): 2 .. 5
+  meta: Dict{Symbol, Any}()
+```
+
+See [`EchelleSpectrum`](@ref) and [`IFUSpectrum`](@ref) for working with instances of higher dimensional spectra.
+"""
+const SingleSpectrum = Spectrum{W, F, 1, 1} where {W, F}
+
+#Base.size(spec::SingleSpectrum) = (length(wave(spec)), )
+#Base.IndexStyle(::Type{<:SingleSpectrum}) = IndexLinear()
+
+function Base.getindex(spec::SingleSpectrum, i::Int)
+    return Spectrum([wave(spec)[i]], [flux(spec)[i]], meta(spec))
+end
+
+function Base.getindex(spec::SingleSpectrum, inds)
+    return Spectrum(wave(spec)[inds], flux(spec)[inds], meta(spec))
+end
+
+Base.firstindex(spec::SingleSpectrum) = firstindex(wave(spec))
+Base.lastindex(spec::SingleSpectrum) = lastindex(wave(spec))
+
+function Base.show(io::IO, spec::SingleSpectrum)
+    w = wave(spec)
+    f = flux(spec)
+    println(io, "SingleSpectrum($(eltype(w)), $(eltype(f)))")
+    println(io, "  wave $(size(w)): ", first(w), " .. ", last(w))
+    println(io, "  flux $(size(f)): ", first(f), " .. ", last(f))
+    print(io, "  meta: ", meta(spec))
+end
diff --git a/src/utils.jl b/src/utils.jl
index 41a8598..a6b9b59 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -21,14 +21,16 @@ julia> wave = range(1, 3, length=100)u"μm"
 (1.0:0.020202020202020204:3.0) μm
 
 julia> bb = blackbody(wave, 2000u"K")
-Spectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(μm,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(μm^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}})
-  T: 2000 K
-  name: Blackbody
+SingleSpectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(μm,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(μm^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}})
+  wave (100,): 1.0 μm .. 3.0 μm
+  flux (100,): 89534.30930426194 W μm^-1 m^-2 .. 49010.54557924032 W μm^-1 m^-2
+  meta: Dict{Symbol, Any}(:T => 2000 K, :name => "Blackbody")
 
 julia> blackbody(ustrip.(u"angstrom", wave), 6000)
-Spectrum(Float64, Float64)
-  T: 6000
-  name: Blackbody
+SingleSpectrum(Float64, Float64)
+  wave (100,): 10000.0 .. 30000.0
+  flux (100,): 1190.9562575755397 .. 40.04325690910415
+  meta: Dict{Symbol, Any}(:T => 6000, :name => "Blackbody")
 
 julia> bb.wave[argmax(bb)]
 1.4444444444444444 μm
@@ -56,3 +58,24 @@ _blackbody(wave::AbstractVector{<:Quantity}, T::Quantity) = blackbody(T).(wave)
 Returns a function for calculating blackbody curves.
 """
 blackbody(T::Quantity) = w->2h * c_0^2 / w^5 / (exp(h * c_0 / (w * k_B * T)) - 1)
+
+"""
+    equivalent_width(::AbstractSpectrum)
+
+Calculate the equivalent width of the given continuum-normalized spectrum. Return value has units equal to wavelengths.
+"""
+function equivalent_width(spec::AbstractSpectrum)
+    dx = spec.wave[end] - spec.wave[1]
+    flux = ustrip(line_flux(spec))
+    return dx - flux * unit(dx)
+end
+
+"""
+    line_flux(::AbstractSpectrum)
+
+Calculate the line flux of the given continuum-normalized spectrum. Return value has units equal to flux.
+"""
+function line_flux(spec::AbstractSpectrum)
+    avg_dx = diff(spec.wave)
+    return sum(spec.flux[2:end] .* avg_dx)
+end
diff --git a/test/fitting.jl b/test/fitting.jl
new file mode 100644
index 0000000..c09e9cc
--- /dev/null
+++ b/test/fitting.jl
@@ -0,0 +1,14 @@
+using Spectra: spectrum, continuum, continuum!
+
+@testset "Continuum" begin
+    spec = spectrum([1, 2, 3.], [1, -10, 1.])
+    spec_cont = continuum(spec)
+
+    @test spec_cont.wave == spec.wave
+    @test spec_cont.flux ≈ ones(eltype(spec.flux), length(spec.flux))
+    @test spec_cont.meta[:coeffs] == spec_cont.meta[:coeffs] ≈ [-4.5, 0, 5.5, 0]
+    @test spec_cont.meta[:normalized]
+
+    continuum!(spec)
+    @test spec == spec_cont
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 9d9b503..c26ce8b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,7 +2,7 @@ using ParallelTestRunner: runtests, find_tests, parse_args
 import Spectra
 
 const init_code = quote
-    using Spectra: Spectra, spectrum
+    using Spectra: Spectra, Spectrum, SingleSpectrum, EchelleSpectrum, IFUSpectrum, spectrum
     using Measurements: Measurements, ±
     using Unitful: @u_str, unit, ustrip
     import Random
diff --git a/test/spectrum.jl b/test/spectrum.jl
index 03afb3b..1f135e6 100644
--- a/test/spectrum.jl
+++ b/test/spectrum.jl
@@ -1,5 +1,25 @@
 Random.seed!(8675309)
 
+@testset "Spectrum types" begin
+    wave, flux = collect(1:5), 1:5
+    @test spectrum(wave, flux) isa SingleSpectrum
+    wave[3] = 10
+    @test_throws ArgumentError spectrum(wave, flux)
+    @test_throws ArgumentError spectrum(wave[begin:end-1], flux)
+
+    wave, flux = repeat(1:5, 1, 3), rand(5, 3)
+    @test spectrum(wave, flux) isa EchelleSpectrum
+    wave[3, 3] = 10
+    @test_throws ArgumentError spectrum(wave, flux)
+    @test_throws ArgumentError spectrum(wave[begin:end-1, :], flux)
+
+    wave, flux = collect(1:5), rand(5, 4, 3)
+    @test spectrum(wave, flux) isa IFUSpectrum
+    wave[3] = 10
+    @test_throws ArgumentError spectrum(wave, flux)
+    @test_throws ArgumentError spectrum(wave[begin:end-1], flux)
+end
+
 @testset "Spectrum - Single" begin
     wave = range(1e4, 5e4, length = 1000)
     sigma = randn(size(wave))
@@ -10,12 +30,15 @@ Random.seed!(8675309)
     flux[134] = 1 ± 0.1
 
     spec = spectrum(wave, flux, name = "test spectrum")
+    spec_indexed = spec[begin:end]
 
     @test propertynames(spec) == (:wave, :flux, :meta, :name)
     @test Spectra.wave(spec) == spec.wave
     @test Spectra.flux(spec) == spec.flux
+    @test [s for s in spec] isa Vector{SingleSpectrum{Float64, Measurements.Measurement{Float64}}}
     @test eltype(spec) == eltype(spec.flux)
     @test spec.wave == wave
+    @test spec_indexed.wave == wave
     @test size(spec) === (1000,)
     @test length(spec) == 1000
     @test maximum(spec) == 1000 ± 1
@@ -25,29 +48,32 @@ Random.seed!(8675309)
     @test findmax(spec) == (1000 ± 1, 7)
     @test findmin(spec) == (1 ± 0.1, 134)
     @test spec.flux == flux
+    @test spec_indexed.flux == flux
     @test Measurements.uncertainty.(spec.flux) ≈ sigma
 
     flux_trimmed = flux[200:800]
-    @test_throws AssertionError spectrum(wave, flux_trimmed)
+    @test_throws ArgumentError spectrum(wave, flux_trimmed)
     expected = """
-    Spectrum(Float64, Measurements.Measurement{Float64})
-      name: test spectrum"""
+    SingleSpectrum(Float64, Measurements.Measurement{Float64})
+      wave (1000,): 10000.0 .. 50000.0
+      flux (1000,): 100.0 ± -2.8 .. 100.0 ± 0.6
+      meta: Dict{Symbol, Any}(:name => "test spectrum")"""
     @test sprint(show, spec) == expected
     @test spec.name == "test spectrum"
 end
 
 @testset "Spectrum - Echelle" begin
-    n_orders = 3
     n_wavs = 1000
+    n_orders = 3
     wave_1 = range(1e4, 5e4, length=n_wavs)
-    wave = repeat(wave_1, 1, n_orders)'
+    wave = repeat(wave_1, 1, n_orders)
     sigma = randn(size(wave_1))
     sigma[7] = 1
     sigma[134] = 0.1
     flux_1 = 100 .± sigma
     flux_1[7] = 1000 ± 1
     flux_1[134] = 1 ± 0.1
-    flux = repeat(flux_1, 1, n_orders)'
+    flux = repeat(flux_1, 1, n_orders)
 
     spec = spectrum(wave, flux, name = "Test Echelle Spectrum")
 
@@ -67,26 +93,53 @@ end
     @test Spectra.flux(spec) == spec.flux
     @test eltype(spec) == eltype(spec.flux)
     @test spec.wave == wave
-    @test size(spec) == (n_orders, n_wavs)
-    @test length(spec) == n_orders * n_wavs
+    @test size(spec) == (n_wavs, n_orders)
+    @test length(spec) == n_wavs * n_orders
     @test maximum(spec) == 1000 ± 1
     @test minimum(spec) == 1 ± 0.1
-    @test argmax(spec) == CartesianIndex(1, 7)
-    @test argmin(spec) == CartesianIndex(1, 134)
-    @test findmax(spec) == (1000 ± 1, CartesianIndex(1, 7))
-    @test findmin(spec) == (1 ± 0.1, CartesianIndex(1, 134))
-    @test eachrow(Measurements.uncertainty.(spec.flux)) ≈ fill(sigma, n_orders)
-
-    flux_trimmed = flux[:, 200:800]
-    @test_throws AssertionError spectrum(wave, flux_trimmed)
+    @test argmax(spec) == CartesianIndex(7, 1)
+    @test argmin(spec) == CartesianIndex(134, 1)
+    @test findmax(spec) == (1000 ± 1, CartesianIndex(7, 1))
+    @test findmin(spec) == (1 ± 0.1, CartesianIndex(134, 1))
+    @test eachcol(Measurements.uncertainty.(spec.flux)) ≈ fill(sigma, n_orders)
+
+    flux_trimmed = flux[200:800, :]
+    @test_throws ArgumentError spectrum(wave, flux_trimmed)
     expected = """
     EchelleSpectrum(Float64, Measurements.Measurement{Float64})
       # orders: 3
-      name: Test Echelle Spectrum"""
+      wave (1000, 3): 10000.0 .. 50000.0
+      flux (1000, 3): 100.0 ± -2.8 .. 100.0 ± 0.6
+      meta: Dict{Symbol, Any}(:name => "Test Echelle Spectrum")"""
     @test sprint(show, spec) == expected
     @test spec.name == "Test Echelle Spectrum"
 end
 
+@testset "Spectrum - IFU" begin
+    wave, flux = [20, 40, 120, 160, 200], rand(5, 10, 6)
+
+    spec = spectrum(wave, flux, name = "test spectrum")
+
+    expected = """
+    IFUSpectrum(Int64, Float64)
+      wave (5,): 20 .. 200
+      flux (5, 10, 6): 0.9210599764489846 .. 0.47778429984485815
+      meta: Dict{Symbol, Any}(:name => "test spectrum")"""
+
+    @test sprint(show, spec) == expected
+    @test spec.name == "test spectrum"
+    @test propertynames(spec) == (:wave, :flux, :meta, :name)
+    @test Spectra.wave(spec) == spec.wave
+    @test Spectra.flux(spec) == spec.flux
+    @test eltype(spec) == eltype(spec.flux)
+    @test spec.wave == wave
+    @test size(spec) === (5, 10, 6)
+    @test length(spec) == 300
+    @test spec.flux == flux
+    @test spec[:, 1, 1] isa SingleSpectrum
+    @test spec[:, begin:4, begin:3] isa IFUSpectrum
+end
+
 @testset "Unitful Spectrum - Single" begin
     wave = range(1e4, 5e4, length = 1000)
     sigma = randn(size(wave))
@@ -125,23 +178,25 @@ end
     @test strip_spec.flux == ustrip.(spec.flux)
     @test strip_spec.meta == spec.meta
     expected = """
-    Spectrum(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Unitful.Quantity{Measurements.Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}})
-      name: test"""
+    SingleSpectrum(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Unitful.Quantity{Measurements.Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}})
+      wave (1000,): 10000.0 Å .. 50000.0 Å
+      flux (1000,): 100.0 ± -2.8 W Å^-1 m^-2 .. 100.0 ± 0.6 W Å^-1 m^-2
+      meta: Dict{Symbol, Any}(:name => "test")"""
     @test sprint(show, spec) == expected
 end
 
 @testset "Unitful Spectrum - Echelle" begin
-    n_orders = 3
     n_wavs = 1000
+    n_orders = 3
     wave_1 = range(1e4, 5e4, length=n_wavs)
-    wave = repeat(wave_1, 1, n_orders)'
+    wave = repeat(wave_1, 1, n_orders)
     sigma = randn(size(wave_1))
     sigma[7] = 1
     sigma[134] = 0.1
     flux_1 = 100 .± sigma
     flux_1[7] = 1000 ± 1
     flux_1[134] = 1 ± 0.1
-    flux = repeat(flux_1, 1, n_orders)'
+    flux = repeat(flux_1, 1, n_orders)
 
     wunit = u"angstrom"
     funit = u"W/m^2/angstrom"
@@ -157,14 +212,14 @@ end
     @test Spectra.wave(spec) == spec.wave
     @test Spectra.flux(spec) == spec.flux
     @test eltype(spec) == eltype(spec.flux)
-    @test size(spec) === (n_orders, n_wavs)
+    @test size(spec) === (n_wavs, n_orders)
     @test length(spec) == n_wavs * n_orders
     @test maximum(spec) == (1000 ± 1) * funit
     @test minimum(spec) == (1 ± 0.1) * funit
-    @test argmax(spec) == CartesianIndex(1, 7)
-    @test argmin(spec) == CartesianIndex(1, 134)
-    @test findmax(spec) == ((1000 ± 1) * funit, CartesianIndex(1, 7))
-    @test findmin(spec) == ((1 ± 0.1) * funit, CartesianIndex(1, 134))
+    @test argmax(spec) == CartesianIndex(7, 1)
+    @test argmin(spec) == CartesianIndex(134, 1)
+    @test findmax(spec) == ((1000 ± 1) * funit, CartesianIndex(7, 1))
+    @test findmin(spec) == ((1 ± 0.1) * funit, CartesianIndex(134, 1))
     @test spec.name == "test echelle"
 
     # Test stripping
@@ -176,11 +231,12 @@ end
     @test strip_spec.wave == ustrip.(spec.wave)
     @test strip_spec.flux == ustrip.(spec.flux)
     @test strip_spec.meta == spec.meta
-    sprint(show, spec)
     expected = """
     EchelleSpectrum(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Unitful.Quantity{Measurements.Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}})
       # orders: 3
-      name: test echelle"""
+      wave (1000, 3): 10000.0 Å .. 50000.0 Å
+      flux (1000, 3): 100.0 ± -2.8 W Å^-1 m^-2 .. 100.0 ± 0.6 W Å^-1 m^-2
+      meta: Dict{Symbol, Any}(:name => "test echelle")"""
     @test sprint(show, spec) == expected
 end
 
diff --git a/test/utils.jl b/test/utils.jl
index 46b1ba3..27ddc62 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -1,4 +1,4 @@
-using Spectra: blackbody
+using Spectra: spectrum, blackbody, line_flux, equivalent_width
 
 @testset "Blackbody T=$T" for T in [2000, 4000, 6000]
     wave = range(1e3, 5e4, length = 1000)
@@ -17,3 +17,13 @@ using Spectra: blackbody
     @test bb.T == T
     @test bb.wave[argmax(bb)] ≈ b * u"angstrom*K" / T rtol = 0.01
 end
+
+@testset "Line flux" begin
+    spec = spectrum([3, 4, 5], [6, 7, 8])
+    @test line_flux(spec) == 15
+end
+
+@testset "Equivalent width" begin
+    spec = spectrum([1, 2, 3], [1, -10, 1])
+    @test equivalent_width(spec) == 11
+end