From 1be6e03452726d6b8a8203d4b56c721aa4414bac Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Fri, 10 Jul 2026 00:31:27 -0700 Subject: [PATCH] Initial sketch up --- Project.toml | 12 ++- docs/Project.toml | 1 + docs/make.jl | 8 +- docs/src/spectrum.md | 14 +++ ext/SpectrumBaseDynamicQuantitiesExt.jl | 89 ++++++++++++++++++ ext/SpectrumBaseUnitfulExt.jl | 119 ++++++++++++++++++++++++ src/SpectrumBase.jl | 95 +++++-------------- src/plotting.jl | 2 +- src/transforms/redshift.jl | 37 +++++--- src/utils.jl | 33 ++++--- test/Project.toml | 1 + test/ext/dynamicquantities.jl | 100 ++++++++++++++++++++ 12 files changed, 405 insertions(+), 106 deletions(-) create mode 100644 ext/SpectrumBaseDynamicQuantitiesExt.jl create mode 100644 ext/SpectrumBaseUnitfulExt.jl create mode 100644 test/ext/dynamicquantities.jl diff --git a/Project.toml b/Project.toml index c7706bb..f392b45 100644 --- a/Project.toml +++ b/Project.toml @@ -10,17 +10,21 @@ projects = ["test", "docs"] DustExtinction = "fb44c06c-c62f-5397-83f5-69249e0a3c8e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" -PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" + +[weakdeps] +DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" -UnitfulAstro = "6112ee07-acf9-5e0f-b108-d242c714bf9f" + +[extensions] +SpectrumBaseDynamicQuantitiesExt = "DynamicQuantities" +SpectrumBaseUnitfulExt = "Unitful" [compat] DustExtinction = ">=0.11.7" +DynamicQuantities = "1" LinearAlgebra = "1.5" Measurements = "2" -PhysicalConstants = "0.2" RecipesBase = "1" Unitful = "1" -UnitfulAstro = "1" julia = "1.10" diff --git a/docs/Project.toml b/docs/Project.toml index 4ae2119..9693e45 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -2,6 +2,7 @@ DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821" FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" diff --git a/docs/make.jl b/docs/make.jl index da47359..5b2bd03 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,8 +1,14 @@ using Documenter using SpectrumBase using Unitful +using DynamicQuantities using Measurements +# The Unitful and DynamicQuantities integrations live in package extensions, so +# their docstrings must be pulled from those modules explicitly. +const unitful_ext = Base.get_extension(SpectrumBase, :SpectrumBaseUnitfulExt) +const dynamicquantities_ext = Base.get_extension(SpectrumBase, :SpectrumBaseDynamicQuantitiesExt) + makedocs(sitename = "SpectrumBase.jl", format = Documenter.HTML(; prettyurls = true, @@ -10,7 +16,7 @@ makedocs(sitename = "SpectrumBase.jl", ), authors = "Miles Lucas and contributors.", linkcheck = !("skiplinks" in ARGS), - modules = [SpectrumBase], + modules = [SpectrumBase, unitful_ext, dynamicquantities_ext], pages = [ "Home" => "index.md", "spectrum.md", diff --git a/docs/src/spectrum.md b/docs/src/spectrum.md index 315b621..c615f50 100644 --- a/docs/src/spectrum.md +++ b/docs/src/spectrum.md @@ -59,11 +59,25 @@ SpectrumBase.meta(::AbstractSpectrum) ## Unitful helpers +Loading [Unitful.jl](https://github.com/PainterQubits/Unitful.jl) activates an +extension that lets spectra carry units and adds the following helpers. + ```@docs Unitful.unit Unitful.ustrip ``` +## DynamicQuantities helpers + +Loading [DynamicQuantities.jl](https://github.com/SymbolicML/DynamicQuantities.jl) +activates an analogous extension. Since DynamicQuantities has no unit objects, +the `unit` helper is replaced by a `dimension` one. + +```@docs +DynamicQuantities.ustrip +DynamicQuantities.dimension +``` + ## Utilities ```@docs diff --git a/ext/SpectrumBaseDynamicQuantitiesExt.jl b/ext/SpectrumBaseDynamicQuantitiesExt.jl new file mode 100644 index 0000000..9d2ae1d --- /dev/null +++ b/ext/SpectrumBaseDynamicQuantitiesExt.jl @@ -0,0 +1,89 @@ +module SpectrumBaseDynamicQuantitiesExt + +using SpectrumBase +using SpectrumBase: AbstractSpectrum, Spectrum, spectrum, spectral_axis, flux_axis, meta, blackbody +import SpectrumBase: _ustrip, _velocity_ratio, _axis_factor +using DynamicQuantities: DynamicQuantities, @u_str, ustrip, dimension, uexpand, + UnionAbstractQuantity, Constants + +# DynamicQuantities does not encode a quantity's dimension in its type (all +# `Quantity{Float64,Dimensions}` share one concrete type), so — unlike the +# Unitful backend — every check here inspects the *value*. `uexpand` normalizes +# symbolic units (`us"..."`) to base `Dimensions` so the comparisons below work +# regardless of how the axis was constructed. +const _LENGTH_DIM = dimension(u"m") +const _ENERGY_DIM = dimension(u"J") + +_dimension(q) = dimension(uexpand(q)) +_is_wavelength(q) = _dimension(q) == _LENGTH_DIM +_is_energy(q) = _dimension(q) == _ENERGY_DIM + +# --- units-agnostic hooks from the core ------------------------------------- + +_ustrip(x::UnionAbstractQuantity) = ustrip(x) + +_velocity_ratio(v::UnionAbstractQuantity) = ustrip(uexpand(v / Constants.c)) + +# Resolved from the axis value rather than the element type (see note above): a +# wavelength axis stretches, an energy axis compresses. +function _axis_factor(spec::AbstractSpectrum{<:UnionAbstractQuantity}, factor) + q = first(spectral_axis(spec)) + _is_wavelength(q) && return factor + _is_energy(q) && return inv(factor) + throw(ArgumentError("cannot shift a spectral axis with dimension $(_dimension(q)); expected a wavelength or energy")) +end + +# --- spectrum constructors for DynamicQuantities axes ----------------------- + +function _assert_spectral_dimension(axis) + q = first(axis) + @assert (_is_wavelength(q) || _is_energy(q)) "spectral_axis not recognized as having dimensions of wavelength or energy." +end + +function SpectrumBase.spectrum(spectral_axis::AbstractVector{<:UnionAbstractQuantity}, flux_axis::AbstractVector{<:UnionAbstractQuantity}; kwds...) + _assert_spectral_dimension(spectral_axis) + Spectrum(spectral_axis, flux_axis, Dict{Symbol, Any}(kwds)) +end + +function SpectrumBase.spectrum(spectral_axis::AbstractMatrix{<:UnionAbstractQuantity}, flux_axis::AbstractMatrix{<:UnionAbstractQuantity}; kwds...) + _assert_spectral_dimension(spectral_axis) + Spectrum(spectral_axis, flux_axis, Dict{Symbol, Any}(kwds)) +end + +# --- ustrip / dimension on whole spectra ------------------------------------ + +""" + DynamicQuantities.ustrip(::AbstractSpectrum) + +Remove the units from a spectrum, returning a plain-number spectrum with the axis +values expressed in SI base units. The DynamicQuantities analog of the Unitful +`ustrip(::AbstractSpectrum)` method. +""" +DynamicQuantities.ustrip(spec::AbstractSpectrum) = spectrum(ustrip.(spectral_axis(spec)), ustrip.(flux_axis(spec)); meta(spec)...) + +""" + DynamicQuantities.dimension(::AbstractSpectrum) + +Get the dimensions of a spectrum as a tuple of the spectral-axis and flux-axis +dimensions. The DynamicQuantities analog of the Unitful `unit(::AbstractSpectrum)` +method (DynamicQuantities has no unit objects, only dimensions). +""" +DynamicQuantities.dimension(spec::AbstractSpectrum) = _dimension(first(spectral_axis(spec))), _dimension(first(flux_axis(spec))) + +# --- blackbody with DynamicQuantities arguments ----------------------------- + +function SpectrumBase.blackbody(wave::AbstractVector{<:UnionAbstractQuantity}, T::UnionAbstractQuantity) + flux = SpectrumBase.blackbody(T).(wave) + return spectrum(wave, flux, name = "Blackbody", T = T) +end + +""" + blackbody(T::DynamicQuantities.UnionAbstractQuantity) + +Returns a function for calculating blackbody curves. The flux is returned in SI +base units (W m⁻² m⁻¹). +""" +SpectrumBase.blackbody(T::UnionAbstractQuantity) = + w -> 2 * Constants.h * Constants.c^2 / w^5 / (exp(Constants.h * Constants.c / (w * Constants.k_B * T)) - 1) + +end # module diff --git a/ext/SpectrumBaseUnitfulExt.jl b/ext/SpectrumBaseUnitfulExt.jl new file mode 100644 index 0000000..2855920 --- /dev/null +++ b/ext/SpectrumBaseUnitfulExt.jl @@ -0,0 +1,119 @@ +module SpectrumBaseUnitfulExt + +using SpectrumBase +using SpectrumBase: AbstractSpectrum, Spectrum, spectrum, spectral_axis, flux_axis, meta, blackbody +import SpectrumBase: _ustrip, _velocity_ratio, _axis_factor +using Unitful: Unitful, @u_str, ustrip, unit, dimension, Quantity, Length, Energy, NoUnits + +# Unit-ful physical constants, rebuilt from the SI floats that live in +# `SpectrumBase` so the Unitful backend agrees with the numeric core to the last +# digit. +const h = SpectrumBase.PLANCK_CONSTANT * u"J*s" +const c_0 = SpectrumBase.SPEED_OF_LIGHT * u"m/s" +const k_B = SpectrumBase.BOLTZMANN_CONSTANT * u"J/K" + +# Wavelength (𝐋) or energy (𝐋²𝐌𝐓⁻²) — the only spectral axes we accept. +const _WAVELENGTH_OR_ENERGY = (u"𝐋", u"𝐋^2 * 𝐌 * 𝐓^-2") + +# --- units-agnostic hooks from the core ------------------------------------- + +_ustrip(x::Quantity) = ustrip(x) + +_velocity_ratio(v::Quantity) = ustrip(NoUnits, v / c_0) + +_axis_factor(::Type{<:Length}, factor) = factor # Wavelength: stretches +_axis_factor(::Type{<:Energy}, factor) = inv(factor) # Energy: compresses +_axis_factor(::Type{S}, factor) where {S <: Quantity} = + throw(ArgumentError("cannot shift a spectral axis with dimension $(dimension(S)); expected a wavelength or energy")) + +# --- spectrum constructors for Unitful axes --------------------------------- + +function SpectrumBase.spectrum(spectral_axis::AbstractVector{<:Quantity}, flux_axis::AbstractVector{<:Quantity}; kwds...) + @assert dimension(eltype(spectral_axis)) ∈ _WAVELENGTH_OR_ENERGY "spectral_axis not recognized as having dimensions of wavelength or energy." + Spectrum(spectral_axis, flux_axis, Dict{Symbol, Any}(kwds)) +end + +function SpectrumBase.spectrum(spectral_axis::AbstractMatrix{<:Quantity}, flux_axis::AbstractMatrix{<:Quantity}; kwds...) + @assert dimension(eltype(spectral_axis)) ∈ _WAVELENGTH_OR_ENERGY "spectral_axis not recognized as having dimensions of wavelength or energy." + Spectrum(spectral_axis, flux_axis, Dict{Symbol, Any}(kwds)) +end + +# --- ustrip / unit on whole spectra ----------------------------------------- + +""" + 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(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Unitful.Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}) + spectral axis (1000,): 10000.0 Å .. 30000.0 Å + flux axis (1000,): 99999.76809093042 W Å^-1 m^-2 .. 300000.2474309158 W Å^-1 m^-2 + meta: Dict{Symbol, Any}() + +julia> ustrip(spec) +SingleSpectrum(Float64, Float64) + spectral axis (1000,): 10000.0 .. 30000.0 + flux axis (1000,): 99999.76809093042 .. 300000.2474309158 + meta: Dict{Symbol, Any}() +``` +""" +Unitful.ustrip(spec::AbstractSpectrum) = spectrum(ustrip.(spectral_axis(spec)), ustrip.(flux_axis(spec)); meta(spec)...) + +""" + Unitful.unit(::AbstractSpectrum) + +Get the units of a spectrum. Returns a tuple of the spectral axis 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(spectral_axis(spec))), unit(eltype(flux_axis(spec))) + +# --- blackbody with Unitful arguments --------------------------------------- + +function SpectrumBase.blackbody(wave::AbstractVector{<:Quantity}, T::Quantity) + out_unit = u"W/m^2" / unit(eltype(wave)) + flux = _blackbody(wave, T) .|> out_unit + return spectrum(wave, flux, name = "Blackbody", T = T) +end + +""" + blackbody(T::Unitful.Quantity) + +Returns a function for calculating blackbody curves. +""" +SpectrumBase.blackbody(T::Quantity) = w -> 2h * c_0^2 / w^5 / (exp(h * c_0 / (w * k_B * T)) - 1) + +_blackbody(wave::AbstractVector{<:Quantity}, T::Quantity) = SpectrumBase.blackbody(T).(wave) + +end # module diff --git a/src/SpectrumBase.jl b/src/SpectrumBase.jl index e5f2795..686f403 100644 --- a/src/SpectrumBase.jl +++ b/src/SpectrumBase.jl @@ -18,8 +18,19 @@ export blackbody #, line_flux, equivalent_width using RecipesBase: @recipe using Measurements: Measurements -using Unitful: Unitful, @u_str, ustrip, unit, dimension -using PhysicalConstants.CODATA2018: h, c_0, k_B + +# Physical constants, CODATA2018 / SI-2019 exact values in SI base units. These are +# the single source of truth for the numeric core and for both units extensions +# (`SpectrumBaseUnitfulExt`, `SpectrumBaseDynamicQuantitiesExt`), which re-attach +# units to these same numbers so every backend agrees to the last digit. +const PLANCK_CONSTANT = 6.62607015e-34 # h, J s +const SPEED_OF_LIGHT = 2.99792458e8 # c_0, m s^-1 +const BOLTZMANN_CONSTANT = 1.380649e-23 # k_B, J K^-1 + +# Strip units from a scalar. The unitless core is the identity; the Unitful and +# DynamicQuantities extensions specialize this on their quantity types so that +# plotting and other numeric consumers see plain numbers. +_ustrip(x) = x """ AbstractSpectrum{S <: Number, F <: Number} @@ -143,67 +154,16 @@ Base.:-(s::AbstractSpectrum, o::AbstractSpectrum) = s - o # Satisfy Aqua # Multi-Spectrum Base.:+(s::T, o::T) where {T <: AbstractSpectrum} = T(spectral_axis(s), flux_axis(s) .+ flux_axis(o), meta(s)) Base.:*(s::T, o::T) where {T <: AbstractSpectrum} = T(spectral_axis(s), flux_axis(s) .* flux_axis(o), meta(s)) -Base.:/(s::T, o::T) where {T <: AbstractSpectrum} = T(spectral_axis(s), flux_axis(s) ./ flux_axis(o) * unit(s)[2], meta(s)) +# Dividing two spectra cancels the flux units; multiply back a unit-valued one so +# the ratio keeps flux units. `oneunit` is units-agnostic (it is `1` for plain +# reals), which keeps this method in the unitless core. +Base.:/(s::T, o::T) where {T <: AbstractSpectrum} = T(spectral_axis(s), flux_axis(s) ./ flux_axis(o) * oneunit(eltype(flux_axis(s))), meta(s)) Base.:-(s::T, o::T) where {T <: AbstractSpectrum} = T(spectral_axis(s), flux_axis(s) .- flux_axis(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(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Unitful.Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}) - spectral axis (1000,): 10000.0 Å .. 30000.0 Å - flux axis (1000,): 99999.76809093042 W Å^-1 m^-2 .. 300000.2474309158 W Å^-1 m^-2 - meta: Dict{Symbol, Any}() - -julia> ustrip(spec) -SingleSpectrum(Float64, Float64) - spectral axis (1000,): 10000.0 .. 30000.0 - flux axis (1000,): 99999.76809093042 .. 300000.2474309158 - meta: Dict{Symbol, Any}() -``` -""" -Unitful.ustrip(spec::AbstractSpectrum) = spectrum(ustrip.(spectral_axis(spec)), ustrip.(flux_axis(spec)); meta(spec)...) - -""" - Unitful.unit(::AbstractSpectrum) - -Get the units of a spectrum. Returns a tuple of the spectral axis 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(spectral_axis(spec))), unit(eltype(flux_axis(spec))) +# `ustrip(spec)` and `unit(spec)` (and the DynamicQuantities `dimension(spec)` +# analog) live in the units extensions, since they extend functions owned by the +# respective units package. See `ext/SpectrumBaseUnitfulExt.jl` and +# `ext/SpectrumBaseDynamicQuantitiesExt.jl`. # Spectrum types and basic arithmetic include("spectrum_single.jl") @@ -289,15 +249,10 @@ function spectrum(spectral_axis::AbstractMatrix{<:Real}, flux_axis::AbstractMatr 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." - 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." - Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds)) -end +# `spectrum(::AbstractArray{<:Quantity}, ...)` methods, which validate that the +# spectral axis carries wavelength or energy dimensions, live in the units +# extensions (`ext/SpectrumBaseUnitfulExt.jl`, +# `ext/SpectrumBaseDynamicQuantitiesExt.jl`). include("utils.jl") include("transforms/resampler.jl") diff --git a/src/plotting.jl b/src/plotting.jl index 9fef960..4a9b16c 100644 --- a/src/plotting.jl +++ b/src/plotting.jl @@ -3,7 +3,7 @@ xlabel --> "wave" ylabel --> "flux density" label --> "" - ustrip.(spectral_axis(spec)), Measurements.value.(ustrip.(flux_axis(spec))) + _ustrip.(spectral_axis(spec)), Measurements.value.(_ustrip.(flux_axis(spec))) end @recipe function f(spec::EchelleSpectrum) diff --git a/src/transforms/redshift.jl b/src/transforms/redshift.jl index 762d3a5..3897f82 100644 --- a/src/transforms/redshift.jl +++ b/src/transforms/redshift.jl @@ -1,7 +1,9 @@ # Dimensionless velocity β ≡ v / c against the speed of light in vacuum. -# A plain number is taken as km/s, the usual unit for radial velocities. -_velocity_ratio(v::Unitful.Quantity) = ustrip(Unitful.NoUnits, v / c_0) -_velocity_ratio(v::Real) = _velocity_ratio(v * u"km/s") +# A plain number is taken as km/s, the usual unit for radial velocities; the +# `Quantity` methods (which honor the velocity's own units) live in the units +# extensions. +const _SPEED_OF_LIGHT_KMPS = SPEED_OF_LIGHT / 1000 +_velocity_ratio(v::Real) = v / _SPEED_OF_LIGHT_KMPS # Doppler shift function _doppler_factor(v; relativistic) @@ -11,22 +13,31 @@ function _doppler_factor(v; relativistic) return sqrt((1 + β) / (1 - β)) 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.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")) -_axis_factor(spec::AbstractSpectrum, factor) = _axis_factor(eltype(spectral_axis(spec)), factor) +# Resolve the (1 + z) stretch factor to the scalar actually multiplied onto a +# spectrum's spectral axis. A unitless axis is assumed to be a wavelength +# (stretches). The dimension-aware behavior for `Quantity` axes (wavelength +# stretches, energy compresses) lives in the units extensions: the Unitful +# backend dispatches on the axis element type, while the DynamicQuantities +# backend inspects the axis value (its dimension is not encoded in the type), so +# both hook in through the `spec`-based method below. +_axis_factor(::Type{<:Real}, factor) = factor +# Dispatch on the spectral-axis element type parameter `S` (known at compile time) +# rather than the runtime axis value, so the factor's type stays inferable. The +# DynamicQuantities extension, whose element type does not carry its dimension, +# instead specializes the more-specific `AbstractSpectrum{<:AbstractQuantity}` +# method and resolves the factor from the value. +_axis_factor(spec::AbstractSpectrum{S}, factor) where {S} = _axis_factor(S, factor) # Return a copy of `spec` whose spectral axis is scaled by `factor`. Rebuilding # through the constructor (rather than reassigning the field) lets the element # type promote, e.g., an integer axis becomes floating point, and re-runs the # `Spectrum` invariant checks. `R` and the type assertions recover the inference -# that is otherwise lost through the `getproperty` overload. +# that is otherwise lost through the `getproperty` overload. `R` is derived with +# `promote_op` rather than `oneunit(S)`, since a DynamicQuantities element type +# does not carry its dimension in the type and so has no `oneunit(::Type)`. function _scale_spectral_axis(spec::Spectrum{S, F, M, N}, factor) where {S, F, M, N} - f = _axis_factor(S, factor) - R = typeof(oneunit(S) * f) + f = _axis_factor(spec, factor) + R = Base.promote_op(*, S, typeof(f)) axis = (spectral_axis(spec) .* f)::AbstractArray{R, M} flux = copy(flux_axis(spec))::AbstractArray{F, N} return Spectrum(axis, flux, deepcopy(meta(spec))) diff --git a/src/utils.jl b/src/utils.jl index 3b69a16..fc23952 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -10,7 +10,10 @@ Create a blackbody spectrum using Planck's law. The curve follows the mathematic If `wave` and `T` are not `Unitful.Quantity`, they are assumed to be in angstrom and Kelvin, and the returned flux will be in units `W m^-2 Å^-1`. -The physical constants are calculated using [PhysicalConstants.jl](https://github.com/juliaphysics/physicalconstants.jl), specifically the CODATA2018 measurement set. +The physical constants are the CODATA2018 / SI-2019 exact values. + +When `wave` and `T` are `Unitful` or `DynamicQuantities` quantities, the corresponding +units extension handles the calculation and the flux carries the appropriate units. # References @@ -33,7 +36,7 @@ SingleSpectrum(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(μm,), 𝐋, n julia> blackbody(ustrip.(u"angstrom", wave), 6000) SingleSpectrum(Float64, Float64) spectral axis (100,): 10000.0 .. 30000.0 - flux axis (100,): 1190.9562575755397 .. 40.04325690910415 + flux axis (100,): 1190.9562575755401 .. 40.043256909104144 meta: Dict{Symbol, Any}(:T => 6000, :name => "Blackbody") julia> spectral_axis(bb)[argmax(bb)] @@ -43,25 +46,21 @@ julia> 2898u"μm*K" / bb.T # See if it matches up with Wien's law 1.449 μm ``` """ -function blackbody(wave::AbstractVector{<:Unitful.Quantity}, T::Unitful.Quantity) - out_unit = u"W/m^2" / unit(eltype(wave)) - flux = _blackbody(wave, T) .|> out_unit - return spectrum(wave, flux, name = "Blackbody", T = T) -end - function blackbody(wave::AbstractVector{<:Real}, T::Real) - flux = ustrip.(u"W/m^2/angstrom", _blackbody(wave * u"angstrom", T * u"K")) + flux = _planck_per_angstrom.(wave, T) return spectrum(wave, flux, name = "Blackbody", T = T) end -_blackbody(wave::AbstractVector{<:Unitful.Quantity}, T::Unitful.Quantity) = blackbody(T).(wave) - -""" - blackbody(T::Unitful.Quantity) - -Returns a function for calculating blackbody curves. -""" -blackbody(T::Unitful.Quantity) = w->2h * c_0^2 / w^5 / (exp(h * c_0 / (w * k_B * T)) - 1) +# Planck's law Bλ(T) for wavelength `λ` in angstrom and temperature `T` in Kelvin, +# returned in W m^-2 Å^-1. Uses the SI constants defined in `SpectrumBase`; the +# unit-aware `blackbody` methods (Unitful/DynamicQuantities `Quantity` arguments) +# live in the units extensions. +function _planck_per_angstrom(λ, T) + λ_m = λ * 1e-10 # angstrom -> m + spectral_radiance = 2 * PLANCK_CONSTANT * SPEED_OF_LIGHT^2 / λ_m^5 / + expm1(PLANCK_CONSTANT * SPEED_OF_LIGHT / (λ_m * BOLTZMANN_CONSTANT * T)) + return spectral_radiance * 1e-10 # W m^-2 m^-1 -> W m^-2 Å^-1 +end #""" # equivalent_width(::AbstractSpectrum) diff --git a/test/Project.toml b/test/Project.toml index 7b080ef..f3d3a7c 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -3,6 +3,7 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DustExtinction = "fb44c06c-c62f-5397-83f5-69249e0a3c8e" +DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc" diff --git a/test/ext/dynamicquantities.jl b/test/ext/dynamicquantities.jl new file mode 100644 index 0000000..bb8e430 --- /dev/null +++ b/test/ext/dynamicquantities.jl @@ -0,0 +1,100 @@ +import DynamicQuantities as DQ +using SpectrumBase: redshift, redshift!, doppler_shift, doppler_shift!, blackbody + +const C_MPS = 299792458.0 +const C_KMPS = C_MPS / 1000 + +# DynamicQuantities normalizes to SI base units, so `angstrom` and `keV` are +# stored as their SI-base values; `DQ.ustrip` returns those base-unit numbers. +Å = DQ.u"angstrom" +keV = DQ.Constants.keV +Wm2 = DQ.u"W/m^2" + +@testset "DynamicQuantities extension" begin + @testset "spectrum constructor" begin + wave = collect(range(4000.0, 8000.0, length = 10)) .* Å + flux = ones(10) .* Wm2 + spec = spectrum(wave, flux) + @test spec isa SingleSpectrum + + # energy axis is also accepted + energy = collect(range(1.0, 5.0, length = 10)) .* keV + @test spectrum(energy, flux) isa SingleSpectrum + + # a non-wavelength/energy axis is rejected + bad = collect(range(1.0, 10.0, length = 10)) .* DQ.u"s" + @test_throws AssertionError spectrum(bad, flux) + end + + @testset "Wavelength axis stretches" begin + wave = collect(range(4000.0, 8000.0, length = 10)) .* Å + spec = spectrum(wave, ones(10) .* Wm2) + + shifted = @inferred redshift(spec, 0.5) + @test all(spectral_axis(shifted) .≈ wave .* 1.5) + @test flux_axis(shifted) == flux_axis(spec) + + # roundtrip + @test all(spectral_axis(redshift(shifted, -1 / 3)) .≈ wave) + end + + @testset "Energy axis compresses" begin + energy = collect(range(1.0, 5.0, length = 10)) .* keV + spec = spectrum(energy, ones(10) .* Wm2) + + shifted = @inferred redshift(spec, 0.5) + @test all(spectral_axis(shifted) .≈ energy ./ 1.5) + @test all(spectral_axis(redshift(shifted, -1 / 3)) .≈ energy) + end + + @testset "redshift in-place" begin + wave = collect(range(4000.0, 8000.0, length = 10)) .* Å + spec = spectrum(wave, ones(10) .* Wm2) + expected = spectral_axis(spec) .* 1.5 + result = @inferred redshift!(spec, 0.5) + @test result === spec + @test all(spectral_axis(spec) .≈ expected) + end + + @testset "Doppler shift" begin + wave = collect(range(4000.0, 8000.0, length = 10)) .* Å + spec = spectrum(wave, ones(10) .* Wm2) + + # bare Real velocity is km/s, as in the numeric core + v = 100.0 + @test all(spectral_axis(doppler_shift(spec, v)) .≈ wave .* (1 + v / C_KMPS)) + + # DynamicQuantities velocity honors its own units + vq = 100.0 * DQ.u"km/s" + β = 1e5 / C_MPS + @test all(spectral_axis(doppler_shift(spec, vq)) .≈ wave .* (1 + β)) + + shifted_rel = doppler_shift(spec, vq; relativistic = true) + @test all(spectral_axis(shifted_rel) .≈ wave .* sqrt((1 + β) / (1 - β))) + + # superluminal is rejected in the relativistic case + @test_throws DomainError doppler_shift(spec, DQ.Constants.c; relativistic = true) + end + + @testset "ustrip and dimension" begin + wave = collect(range(4000.0, 8000.0, length = 10)) .* Å + spec = spectrum(wave, ones(10) .* Wm2) + + stripped = DQ.ustrip(spec) + @test eltype(spectral_axis(stripped)) <: Real + @test spectral_axis(stripped) ≈ DQ.ustrip.(wave) + + w_dim, f_dim = DQ.dimension(spec) + @test w_dim == DQ.dimension(DQ.u"m") + @test f_dim == DQ.dimension(Wm2) + end + + @testset "blackbody" begin + wave = collect(range(1.0, 3.0, length = 10)) .* DQ.u"μm" + bb = blackbody(wave, 2000.0 * DQ.u"K") + @test bb isa SingleSpectrum + @test spectral_axis(bb) == wave + @test all(DQ.ustrip.(flux_axis(bb)) .> 0) + @test bb.T == 2000.0 * DQ.u"K" + end +end