From 602e4ac81ecfde1d25385e27c2586516164dd18a Mon Sep 17 00:00:00 2001 From: icweaver Date: Sat, 1 Nov 2025 01:24:46 -0700 Subject: [PATCH 01/31] Spectrum: added indexing --- src/spectrum.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/spectrum.jl b/src/spectrum.jl index d636cd3..077c83a 100644 --- a/src/spectrum.jl +++ b/src/spectrum.jl @@ -11,6 +11,17 @@ end Spectrum(wave, flux, meta::Dict{Symbol,Any}) = Spectrum(collect(wave), collect(flux), meta) +Base.size(spec::Spectrum) = (length(spec.wave), ) +Base.IndexStyle(::Type{<:Spectrum}) = IndexLinear() + +function Base.getindex(spec::Spectrum, i::Int) + return Spectrum([spec.wave[i]], [spec.flux[i]], spec.meta) +end + +function Base.getindex(spec::Spectrum, inds) + return Spectrum(spec.wave[inds], spec.flux[inds], spec.meta) +end + function Base.show(io::IO, spec::Spectrum) print(io, "Spectrum($(eltype(spec.wave)), $(eltype(spec.flux)))") for (key, val) in spec.meta From aa7ef4849c03db272814ef1a11cad90b7270bf06 Mon Sep 17 00:00:00 2001 From: icweaver Date: Sat, 1 Nov 2025 01:42:17 -0700 Subject: [PATCH 02/31] quick index test --- test/spectrum.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/spectrum.jl b/test/spectrum.jl index 03afb3b..b13b0a2 100644 --- a/test/spectrum.jl +++ b/test/spectrum.jl @@ -10,12 +10,14 @@ 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 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,6 +27,7 @@ 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] From da7585eb985fc10541dbd10c3568096fc02ae521 Mon Sep 17 00:00:00 2001 From: icweaver Date: Sat, 1 Nov 2025 01:57:16 -0700 Subject: [PATCH 03/31] begin, end --- src/spectrum.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/spectrum.jl b/src/spectrum.jl index 077c83a..8b7a469 100644 --- a/src/spectrum.jl +++ b/src/spectrum.jl @@ -22,6 +22,9 @@ function Base.getindex(spec::Spectrum, inds) return Spectrum(spec.wave[inds], spec.flux[inds], spec.meta) end +Base.firstindex(spec::Spectrum) = firstindex(spec.wave) +Base.lastindex(spec::Spectrum) = lastindex(spec.wave) + function Base.show(io::IO, spec::Spectrum) print(io, "Spectrum($(eltype(spec.wave)), $(eltype(spec.flux)))") for (key, val) in spec.meta From a5f3cf0a1ccb385d9c09f856ade1e4e86bff9a6f Mon Sep 17 00:00:00 2001 From: icweaver Date: Sun, 2 Nov 2025 01:49:41 -0800 Subject: [PATCH 04/31] start moving EchelleSpectrum intro Spectrum --- src/EchelleSpectrum.jl | 31 ------------------------------- src/Spectra.jl | 1 - src/spectrum.jl | 32 ++++++++++++++++++++++++++++++++ 3 files changed, 32 insertions(+), 32 deletions(-) delete mode 100644 src/EchelleSpectrum.jl 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..4ee24fc 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -17,7 +17,6 @@ include("common.jl") # Spectrum types and basic arithmetic include("spectrum.jl") -include("EchelleSpectrum.jl") """ spectrum(wave, flux; kwds...) diff --git a/src/spectrum.jl b/src/spectrum.jl index 8b7a469..bec5f09 100644 --- a/src/spectrum.jl +++ b/src/spectrum.jl @@ -31,3 +31,35 @@ function Base.show(io::IO, spec::Spectrum) print(io, "\n $key: $val") end end + +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) From 36ea6b6252f595909eddc3dcb43b4b2521f772cc Mon Sep 17 00:00:00 2001 From: icweaver Date: Sun, 2 Nov 2025 03:35:19 -0800 Subject: [PATCH 05/31] not totally hating this --- src/Spectra.jl | 41 +++++++++++++------------- src/common.jl | 62 ++++++++++++++++++++++++++------------- src/spectrum.jl | 65 ----------------------------------------- src/spectrum_echelle.jl | 24 +++++++++++++++ src/spectrum_single.jl | 22 ++++++++++++++ 5 files changed, 108 insertions(+), 106 deletions(-) delete mode 100644 src/spectrum.jl create mode 100644 src/spectrum_echelle.jl create mode 100644 src/spectrum_single.jl diff --git a/src/Spectra.jl b/src/Spectra.jl index 4ee24fc..2217bf2 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -1,7 +1,7 @@ module Spectra # common.jl -export AbstractSpectrum, spectrum +export AbstractSpectrum, SingleSpectrum, EchelleSpectrum, spectrum # utils.jl export blackbody # transforms/redden.jl @@ -12,11 +12,12 @@ 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 +# AbstractSpectrum, Spectrum, and common functionality include("common.jl") # Spectrum types and basic arithmetic -include("spectrum.jl") +include("spectrum_single.jl") +include("spectrum_echelle.jl") """ spectrum(wave, flux; kwds...) @@ -73,27 +74,27 @@ function spectrum(wave::AbstractVector{<:Real}, flux::AbstractVector{<:Real}; kw 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" - 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" +# 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)) -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 -include("utils.jl") -include("transforms/transforms.jl") -include("plotting.jl") -include("fitting/fitting.jl") +#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)) +#end +# +## tools +#include("utils.jl") +#include("transforms/transforms.jl") +#include("plotting.jl") +#include("fitting/fitting.jl") end # module diff --git a/src/common.jl b/src/common.jl index 5c0d20c..1b6da2a 100644 --- a/src/common.jl +++ b/src/common.jl @@ -9,6 +9,19 @@ An abstract holder for astronomical spectra. All types inheriting from this must """ 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, N} <: AbstractSpectrum{W, F} + wave::AbstractArray{W, N} + flux::AbstractArray{F, N} + meta::Dict{Symbol,Any} +end + +Spectrum(wave, flux, meta::Dict{Symbol, Any}) = Spectrum(collect(wave), collect(flux), meta) + function Base.getproperty(spec::AbstractSpectrum, nm::Symbol) if nm in keys(getfield(spec, :meta)) return getfield(spec, :meta)[nm] @@ -19,7 +32,7 @@ end function Base.propertynames(spec::AbstractSpectrum) natural = (:wave, :flux, :meta) - meta = keys(spec.meta) + meta = keys(meta(spec)) return (natural..., meta...) end @@ -37,32 +50,39 @@ Return the flux of the spectrum. """ flux(spec::AbstractSpectrum) = spec.flux +""" + meta(::AbstractSpectrum) + +Return the meta of the spectrum. +""" +meta(spec::AbstractSpectrum) = spec.meta + # 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) +Base.eltype(spec::AbstractSpectrum) = eltype(flux(spec)) +Base.size(spec::AbstractSpectrum) = size(flux(spec)) +Base.size(spec::AbstractSpectrum, i) = size(flux(spec), i) +Base.length(spec::AbstractSpectrum) = length(flux(spec)) +Base.maximum(spec::AbstractSpectrum) = maximum(flux(spec)) +Base.minimum(spec::AbstractSpectrum) = minimum(flux(spec)) +Base.argmax(spec::AbstractSpectrum) = argmax(flux(spec)) +Base.argmin(spec::AbstractSpectrum) = argmin(flux(spec)) +Base.findmax(spec::AbstractSpectrum) = findmax(flux(spec)) +Base.findmin(spec::AbstractSpectrum) = findmin(flux(spec)) # 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::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(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) +Base.:+(s::T, o::T) where {T <: AbstractSpectrum} = T(wave(s), flux(s) .+ flux(s), meta(s)) +Base.:*(s::T, o::T) where {T <: AbstractSpectrum} = T(wave(s), flux(s) .* flux(s), meta(s)) +Base.:/(s::T, o::T) where {T <: AbstractSpectrum} = T(wave(s), flux(s) ./ flux(s) * unit(s)[2], meta(s)) +Base.:-(s::T, o::T) where {T <: AbstractSpectrum} = T(wave(s), flux(s) .- flux(s), meta(s)) """ Unitful.ustrip(::AbstractSpectrum) @@ -84,7 +104,7 @@ julia> ustrip(spec) Spectrum(Float64, Float64) ``` """ -Unitful.ustrip(spec::AbstractSpectrum) = spectrum(ustrip.(spec.wave), ustrip.(spec.flux); spec.meta...) +Unitful.ustrip(spec::AbstractSpectrum) = spectrum(ustrip.(wave(spec)), ustrip.(flux(spec)); meta(spec)...) """ Unitful.unit(::AbstractSpectrum) @@ -105,4 +125,4 @@ julia> w_unit, f_unit = unit(spec) (Å, W Å^-1 m^-2) ``` """ -Unitful.unit(spec::AbstractSpectrum) = unit(eltype(spec.wave)), unit(eltype(spec.flux)) +Unitful.unit(spec::AbstractSpectrum) = unit(eltype(wave(spec))), unit(eltype(flux(spec))) diff --git a/src/spectrum.jl b/src/spectrum.jl deleted file mode 100644 index bec5f09..0000000 --- a/src/spectrum.jl +++ /dev/null @@ -1,65 +0,0 @@ -""" - Spectrum <: AbstractSpectrum - -A 1-dimensional spectrum stored as vectors of real numbers. The wavelengths are assumed to be in angstrom. -""" -mutable 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) - -Base.size(spec::Spectrum) = (length(spec.wave), ) -Base.IndexStyle(::Type{<:Spectrum}) = IndexLinear() - -function Base.getindex(spec::Spectrum, i::Int) - return Spectrum([spec.wave[i]], [spec.flux[i]], spec.meta) -end - -function Base.getindex(spec::Spectrum, inds) - return Spectrum(spec.wave[inds], spec.flux[inds], spec.meta) -end - -Base.firstindex(spec::Spectrum) = firstindex(spec.wave) -Base.lastindex(spec::Spectrum) = lastindex(spec.wave) - -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 - -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/spectrum_echelle.jl b/src/spectrum_echelle.jl new file mode 100644 index 0000000..0557744 --- /dev/null +++ b/src/spectrum_echelle.jl @@ -0,0 +1,24 @@ +const EchelleSpectrum = Spectrum{W, F, 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) + println(io, "EchelleSpectrum($(eltype(wave(spec))), $(eltype(flux(spec))))") + println(io, " # orders: $(size(spec, 1))") + println(io, "meta: ", meta(spec)) +end diff --git a/src/spectrum_single.jl b/src/spectrum_single.jl new file mode 100644 index 0000000..4c62e57 --- /dev/null +++ b/src/spectrum_single.jl @@ -0,0 +1,22 @@ +const SingleSpectrum = Spectrum{W, F, 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) + println(io, "SingleSpectrum($(eltype(wave(spec))), $(eltype(flux(spec))))") + println(io, "wave: ", wave(spec)) + println(io, "flux: ", flux(spec)) + println(io, "meta: ", meta(spec)) +end From 337c79378271639fe65c6710f5b893da83c9d821 Mon Sep 17 00:00:00 2001 From: icweaver Date: Sun, 2 Nov 2025 12:33:34 -0800 Subject: [PATCH 06/31] current tests up --- src/Spectra.jl | 32 ++++++++++++++++---------------- src/common.jl | 40 ++++++++++++++++++++-------------------- src/plotting.jl | 2 +- src/spectrum_echelle.jl | 2 +- src/spectrum_single.jl | 6 +++--- test/spectrum.jl | 16 ++++++++++------ 6 files changed, 51 insertions(+), 47 deletions(-) diff --git a/src/Spectra.jl b/src/Spectra.jl index 2217bf2..dc0d8d1 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -74,27 +74,27 @@ function spectrum(wave::AbstractVector{<:Real}, flux::AbstractVector{<:Real}; kw 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" -# 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" + 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" 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)) -#end -# -## tools -#include("utils.jl") -#include("transforms/transforms.jl") -#include("plotting.jl") -#include("fitting/fitting.jl") +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" + Spectrum(wave, flux, Dict{Symbol,Any}(kwds)) +end + +# tools +include("utils.jl") +include("transforms/transforms.jl") +include("plotting.jl") +include("fitting/fitting.jl") end # module diff --git a/src/common.jl b/src/common.jl index 1b6da2a..d84d2fb 100644 --- a/src/common.jl +++ b/src/common.jl @@ -20,22 +20,6 @@ mutable struct Spectrum{W<:Number, F<:Number, N} <: AbstractSpectrum{W, F} meta::Dict{Symbol,Any} end -Spectrum(wave, flux, meta::Dict{Symbol, Any}) = Spectrum(collect(wave), collect(flux), 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) - meta = keys(meta(spec)) - return (natural..., meta...) -end - """ wave(::AbstractSpectrum) @@ -57,6 +41,22 @@ Return the meta of the spectrum. """ meta(spec::AbstractSpectrum) = spec.meta +Spectrum(wave, flux, meta::Dict{Symbol, Any}) = Spectrum(collect(wave), collect(flux), 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.eltype(spec::AbstractSpectrum) = eltype(flux(spec)) Base.size(spec::AbstractSpectrum) = size(flux(spec)) @@ -79,10 +79,10 @@ 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(s), meta(s)) -Base.:*(s::T, o::T) where {T <: AbstractSpectrum} = T(wave(s), flux(s) .* flux(s), meta(s)) -Base.:/(s::T, o::T) where {T <: AbstractSpectrum} = T(wave(s), flux(s) ./ flux(s) * unit(s)[2], meta(s)) -Base.:-(s::T, o::T) where {T <: AbstractSpectrum} = T(wave(s), flux(s) .- flux(s), 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), 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) 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_echelle.jl b/src/spectrum_echelle.jl index 0557744..9a491ee 100644 --- a/src/spectrum_echelle.jl +++ b/src/spectrum_echelle.jl @@ -20,5 +20,5 @@ Base.lastindex(spec::EchelleSpectrum) = lastindex(flux(spec), 1) function Base.show(io::IO, spec::EchelleSpectrum) println(io, "EchelleSpectrum($(eltype(wave(spec))), $(eltype(flux(spec))))") println(io, " # orders: $(size(spec, 1))") - println(io, "meta: ", meta(spec)) + print(io, " meta: ", meta(spec)) end diff --git a/src/spectrum_single.jl b/src/spectrum_single.jl index 4c62e57..7cf771c 100644 --- a/src/spectrum_single.jl +++ b/src/spectrum_single.jl @@ -16,7 +16,7 @@ Base.lastindex(spec::SingleSpectrum) = lastindex(wave(spec)) function Base.show(io::IO, spec::SingleSpectrum) println(io, "SingleSpectrum($(eltype(wave(spec))), $(eltype(flux(spec))))") - println(io, "wave: ", wave(spec)) - println(io, "flux: ", flux(spec)) - println(io, "meta: ", meta(spec)) + println(io, " wave: ", (extrema∘wave)(spec)) + println(io, " flux: ", (extrema∘flux)(spec)) + print(io, " meta: ", meta(spec)) end diff --git a/test/spectrum.jl b/test/spectrum.jl index b13b0a2..ef95e26 100644 --- a/test/spectrum.jl +++ b/test/spectrum.jl @@ -33,8 +33,10 @@ Random.seed!(8675309) flux_trimmed = flux[200:800] @test_throws AssertionError spectrum(wave, flux_trimmed) expected = """ - Spectrum(Float64, Measurements.Measurement{Float64}) - name: test spectrum""" + SingleSpectrum(Float64, Measurements.Measurement{Float64}) + wave: (10000.0, 50000.0) + flux: (1.0 ± 0.1, 1000.0 ± 1.0) + meta: Dict{Symbol, Any}(:name => "test spectrum")""" @test sprint(show, spec) == expected @test spec.name == "test spectrum" end @@ -85,7 +87,7 @@ end expected = """ EchelleSpectrum(Float64, Measurements.Measurement{Float64}) # orders: 3 - name: Test Echelle Spectrum""" + meta: Dict{Symbol, Any}(:name => "Test Echelle Spectrum")""" @test sprint(show, spec) == expected @test spec.name == "Test Echelle Spectrum" end @@ -128,8 +130,10 @@ 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: (10000.0 Å, 50000.0 Å) + flux: (1.0 ± 0.1 W Å^-1 m^-2, 1000.0 ± 1.0 W Å^-1 m^-2) + meta: Dict{Symbol, Any}(:name => "test")""" @test sprint(show, spec) == expected end @@ -183,7 +187,7 @@ end 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""" + meta: Dict{Symbol, Any}(:name => "test echelle")""" @test sprint(show, spec) == expected end From 8b9bf4d8bdaf128eacb68debb23e122baaa0fe3e Mon Sep 17 00:00:00 2001 From: icweaver Date: Sun, 2 Nov 2025 13:31:02 -0800 Subject: [PATCH 07/31] reorg utils, add more array interface + tests --- src/Spectra.jl | 4 +++- src/common.jl | 5 +++++ src/fitting/fitting.jl | 22 ---------------------- src/utils.jl | 21 +++++++++++++++++++++ test/fitting.jl | 14 ++++++++++++++ test/utils.jl | 12 +++++++++++- 6 files changed, 54 insertions(+), 24 deletions(-) create mode 100644 test/fitting.jl diff --git a/src/Spectra.jl b/src/Spectra.jl index dc0d8d1..4510c4d 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -3,7 +3,9 @@ module Spectra # common.jl export AbstractSpectrum, SingleSpectrum, EchelleSpectrum, spectrum # utils.jl -export blackbody +export blackbody, line_flux, equivalent_width +# fitting/fitting.jl +export continuum, continuum! # transforms/redden.jl export redden, redden!, deredden, deredden! diff --git a/src/common.jl b/src/common.jl index d84d2fb..70d830f 100644 --- a/src/common.jl +++ b/src/common.jl @@ -68,6 +68,11 @@ Base.argmax(spec::AbstractSpectrum) = argmax(flux(spec)) Base.argmin(spec::AbstractSpectrum) = argmin(flux(spec)) Base.findmax(spec::AbstractSpectrum) = findmax(flux(spec)) Base.findmin(spec::AbstractSpectrum) = findmin(flux(spec)) +function Base.iterate(spec::AbstractSpectrum, state=0) + state == length(spec) && return nothing + return spec[begin + state], state + 1 +end +Base.:(==)(s::AbstractSpectrum, o::AbstractSpectrum) = wave(s) == wave(o) && flux(s) == flux(o) && meta(s) == meta(o) # Arithmetic Base.:+(s::T, A) where {T <: AbstractSpectrum} = T(wave(s), flux(s) .+ A, meta(s)) 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/utils.jl b/src/utils.jl index 41a8598..984e132 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -56,3 +56,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/utils.jl b/test/utils.jl index 46b1ba3..ee217c2 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, 5, 4], [6, 7, 8]) + @test line_flux(spec) == 6 +end + +@testset "Equivalent width" begin + spec = spectrum([1, 2, 3], [1, -10, 1]) + @test equivalent_width(spec) == 11 +end From 03c3f346211c9e14673a18227f067bf455828961 Mon Sep 17 00:00:00 2001 From: icweaver Date: Sun, 2 Nov 2025 14:59:23 -0800 Subject: [PATCH 08/31] more docs + prettyurl --- docs/make.jl | 4 ++-- docs/src/index.md | 12 ++++++---- docs/src/spectrum.md | 35 ++++++++++++++++++++---------- docs/src/transforms.md | 15 ++++++++++--- docs/src/{analysis.md => utils.md} | 5 +++-- src/Spectra.jl | 17 +++++++++++---- src/common.jl | 24 ++++++++++++-------- src/utils.jl | 14 +++++++----- 8 files changed, 85 insertions(+), 41 deletions(-) rename docs/src/{analysis.md => utils.md} (60%) diff --git a/docs/make.jl b/docs/make.jl index 950496e..df28745 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -7,7 +7,7 @@ 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 +18,7 @@ makedocs(sitename = "Spectra.jl", "spectrum.md", "transforms.md", "fitting.md", - "analysis.md", + "utils.md", "contrib.md", ], warnonly = [:missing_docs], diff --git a/docs/src/index.md b/docs/src/index.md index 92b6a7f..21a28a3 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: (3815.0483f0 Å, 9206.613f0 Å) + flux: (9.416705322265625e-16 erg Å^-1 cm^-2 s^-1, 9.156107177734375e-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: (3815.0483f0 Å, 9206.613f0 Å) + flux: (0.4582525451059156 erg Å^-1 cm^-2 s^-1, 4.402153019328562 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..03431b4 100644 --- a/docs/src/spectrum.md +++ b/docs/src/spectrum.md @@ -14,6 +14,7 @@ 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 ``` @@ -28,25 +29,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 diff --git a/docs/src/transforms.md b/docs/src/transforms.md index a5a05c5..d33df7e 100644 --- a/docs/src/transforms.md +++ b/docs/src/transforms.md @@ -36,10 +36,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: (1.0 μm, 3.0 μm) + flux: (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: (1.0 μm, 3.0 μm) + flux: (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 +56,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: (1.0 μm, 3.0 μm) + flux: (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 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/Spectra.jl b/src/Spectra.jl index 4510c4d..470bfa4 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -33,11 +33,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: (10000.0, 40000.0) + flux: (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: (10000.0, 40000.0) + flux: (100.0, 100.0) + meta: Dict{Symbol, Any}(:name => "Just Noise") julia> spec.name "Just Noise" @@ -56,7 +61,10 @@ julia> sigma = randn(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: (1.0 μm, 4.0 μm) + flux: (100.0 ± 1.2 erg Å^-1 cm^-2 s^-1, 100.0 ± 1.1 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. @@ -69,6 +77,7 @@ julia> flux = ones(10, 100) .* collect(1:10); julia> spec = spectrum(wave, flux) EchelleSpectrum(Float64, Float64) # orders: 10 + meta: Dict{Symbol, Any}() ``` """ function spectrum(wave::AbstractVector{<:Real}, flux::AbstractVector{<:Real}; kwds...) diff --git a/src/common.jl b/src/common.jl index 70d830f..b8972ea 100644 --- a/src/common.jl +++ b/src/common.jl @@ -58,23 +58,23 @@ function Base.propertynames(spec::AbstractSpectrum) end # Collection -Base.eltype(spec::AbstractSpectrum) = eltype(flux(spec)) -Base.size(spec::AbstractSpectrum) = size(flux(spec)) -Base.size(spec::AbstractSpectrum, i) = size(flux(spec), i) -Base.length(spec::AbstractSpectrum) = length(flux(spec)) -Base.maximum(spec::AbstractSpectrum) = maximum(flux(spec)) -Base.minimum(spec::AbstractSpectrum) = minimum(flux(spec)) 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 -Base.:(==)(s::AbstractSpectrum, o::AbstractSpectrum) = wave(s) == wave(o) && flux(s) == flux(o) && meta(s) == meta(o) # 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)) @@ -103,10 +103,16 @@ 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}}) +SingleSpectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}) + wave: (10000.0 Å, 30000.0 Å) + flux: (99999.8952204731 W Å^-1 m^-2, 299999.8866277076 W Å^-1 m^-2) + meta: Dict{Symbol, Any}() julia> ustrip(spec) -Spectrum(Float64, Float64) +SingleSpectrum(Float64, Float64) + wave: (10000.0, 30000.0) + flux: (99999.8952204731, 299999.8866277076) + meta: Dict{Symbol, Any}() ``` """ Unitful.ustrip(spec::AbstractSpectrum) = spectrum(ustrip.(wave(spec)), ustrip.(flux(spec)); meta(spec)...) diff --git a/src/utils.jl b/src/utils.jl index 984e132..e5ad6d5 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: (1.0 μm, 3.0 μm) + flux: (49010.54557924032 W μm^-1 m^-2, 131058.60552995963 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: (10000.0, 30000.0) + flux: (40.04325690910415, 1190.9562575755397) + meta: Dict{Symbol, Any}(:T => 6000, :name => "Blackbody") julia> bb.wave[argmax(bb)] 1.4444444444444444 μm From dd6147567675f1acdb07783e590365b5770c74d5 Mon Sep 17 00:00:00 2001 From: icweaver Date: Sun, 2 Nov 2025 16:22:07 -0800 Subject: [PATCH 09/31] update make.jl --- docs/make.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index df28745..bb712bc 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -27,5 +27,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 ) From 314d6c2b6e09c1aa48b7632cf3376e9f04aee1f8 Mon Sep 17 00:00:00 2001 From: icweaver Date: Mon, 3 Nov 2025 12:36:55 -0800 Subject: [PATCH 10/31] add IFUSpectrum --- src/Spectra.jl | 14 ++++++++++---- src/common.jl | 4 ++-- src/spectrum_echelle.jl | 2 +- src/spectrum_ifu.jl | 33 +++++++++++++++++++++++++++++++++ src/spectrum_single.jl | 2 +- 5 files changed, 47 insertions(+), 8 deletions(-) create mode 100644 src/spectrum_ifu.jl diff --git a/src/Spectra.jl b/src/Spectra.jl index 470bfa4..f34ead2 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -1,7 +1,7 @@ module Spectra # common.jl -export AbstractSpectrum, SingleSpectrum, EchelleSpectrum, spectrum +export AbstractSpectrum, SingleSpectrum, IFUSpectrum, EchelleSpectrum, spectrum # utils.jl export blackbody, line_flux, equivalent_width # fitting/fitting.jl @@ -19,6 +19,7 @@ include("common.jl") # Spectrum types and basic arithmetic include("spectrum_single.jl") +include("spectrum_ifu.jl") include("spectrum_echelle.jl") """ @@ -85,9 +86,8 @@ function spectrum(wave::AbstractVector{<:Real}, flux::AbstractVector{<:Real}; kw 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::AbstractMatrix{<:Real}; kwds...) + @assert size(wave, 1) == size(flux, 2) "wave and flux in each order must have equal size" Spectrum(wave, flux, Dict{Symbol,Any}(kwds)) end @@ -96,6 +96,12 @@ function spectrum(wave::AbstractMatrix{<:Real}, flux::AbstractMatrix{<:Real}; kw 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" + 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" diff --git a/src/common.jl b/src/common.jl index b8972ea..901cc87 100644 --- a/src/common.jl +++ b/src/common.jl @@ -14,9 +14,9 @@ abstract type AbstractSpectrum{W,F} end 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, N} <: AbstractSpectrum{W, F} +mutable struct Spectrum{W<:Number, F<:Number, N, M} <: AbstractSpectrum{W, F} wave::AbstractArray{W, N} - flux::AbstractArray{F, N} + flux::AbstractArray{F, M} meta::Dict{Symbol,Any} end diff --git a/src/spectrum_echelle.jl b/src/spectrum_echelle.jl index 9a491ee..999faab 100644 --- a/src/spectrum_echelle.jl +++ b/src/spectrum_echelle.jl @@ -1,4 +1,4 @@ -const EchelleSpectrum = Spectrum{W, F, 2} where {W, F} +const EchelleSpectrum = Spectrum{W, F, 2, 2} where {W, F} function Base.getindex(spec::EchelleSpectrum, i::Int) w = wave(spec)[i, :] diff --git a/src/spectrum_ifu.jl b/src/spectrum_ifu.jl new file mode 100644 index 0000000..567448c --- /dev/null +++ b/src/spectrum_ifu.jl @@ -0,0 +1,33 @@ +const IFUSpectrum = Spectrum{W, F, 1, 2} where {W, F} + +function Base.getindex(spec::IFUSpectrum, i::Int) + w = wave(spec) + f = flux(spec)[i, :] + m = merge(Dict(:Order => i), meta(spec)) + return Spectrum(w, f, m) +end + +function Base.getindex(spec::IFUSpectrum, i::Int, J::AbstractVector) + w = wave(spec) + f = flux(spec)[i, J] + m = merge(Dict(:Order => i), meta(spec)) + return Spectrum(w, f, m) +end + +function Base.getindex(spec::IFUSpectrum, 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::IFUSpectrum) = firstindex(flux(spec), 1) +Base.firstindex(spec::IFUSpectrum, dim::Int) = firstindex(flux(spec), dim) +Base.lastindex(spec::IFUSpectrum) = lastindex(flux(spec), 1) + +function Base.show(io::IO, spec::IFUSpectrum) + println(io, "IFUSpectrum($(eltype(wave(spec))), $(eltype(flux(spec))))") + println(io, " # orders: $(size(spec, 1))") + println(io, " wave: ", (extrema∘wave)(spec)) + print(io, " meta: ", meta(spec)) +end diff --git a/src/spectrum_single.jl b/src/spectrum_single.jl index 7cf771c..1558da9 100644 --- a/src/spectrum_single.jl +++ b/src/spectrum_single.jl @@ -1,4 +1,4 @@ -const SingleSpectrum = Spectrum{W, F, 1} where {W, F} +const SingleSpectrum = Spectrum{W, F, 1, 1} where {W, F} #Base.size(spec::SingleSpectrum) = (length(wave(spec)), ) #Base.IndexStyle(::Type{<:SingleSpectrum}) = IndexLinear() From 76b6d0d7ffca6ad15f214e004768a578bdaa6357 Mon Sep 17 00:00:00 2001 From: icweaver Date: Tue, 4 Nov 2025 08:06:05 -0800 Subject: [PATCH 11/31] merge common.jl into Spectra.jl --- src/Spectra.jl | 147 +++++++++++++++++++++++++++++++++++++++++++++++-- src/common.jl | 139 ---------------------------------------------- 2 files changed, 143 insertions(+), 143 deletions(-) delete mode 100644 src/common.jl diff --git a/src/Spectra.jl b/src/Spectra.jl index f34ead2..bceb336 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -1,7 +1,9 @@ module Spectra -# common.jl -export AbstractSpectrum, SingleSpectrum, IFUSpectrum, EchelleSpectrum, spectrum +# Uniform API +export AbstractSpectrum, spectrum +# spectra_single.jl, spectra_ifu.jl, spectra_echelle.jl +export SingleSpectrum, IFUSpectrum, EchelleSpectrum # utils.jl export blackbody, line_flux, equivalent_width # fitting/fitting.jl @@ -14,8 +16,145 @@ using Measurements: Measurements, Measurement using Unitful: Unitful, Quantity, @u_str, ustrip, unit, dimension using PhysicalConstants.CODATA2018: h, c_0, k_B -# AbstractSpectrum, Spectrum, 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, N} +- flux::Array{F, N} +- meta::Dict{Symbol, Any} +""" +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, N, M} <: AbstractSpectrum{W, F} + wave::AbstractArray{W, N} + flux::AbstractArray{F, M} + meta::Dict{Symbol,Any} +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 + +""" + meta(::AbstractSpectrum) + +Return the meta of the spectrum. +""" +meta(spec::AbstractSpectrum) = spec.meta + +Spectrum(wave, flux, meta::Dict{Symbol, Any}) = Spectrum(collect(wave), collect(flux), 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 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") +SingleSpectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}) + wave: (10000.0 Å, 30000.0 Å) + flux: (99999.8952204731 W Å^-1 m^-2, 299999.8866277076 W Å^-1 m^-2) + meta: Dict{Symbol, Any}() + +julia> ustrip(spec) +SingleSpectrum(Float64, Float64) + wave: (10000.0, 30000.0) + flux: (99999.8952204731, 299999.8866277076) + 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 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(wave(spec))), unit(eltype(flux(spec))) # Spectrum types and basic arithmetic include("spectrum_single.jl") diff --git a/src/common.jl b/src/common.jl deleted file mode 100644 index 901cc87..0000000 --- a/src/common.jl +++ /dev/null @@ -1,139 +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 - -""" - 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, N, M} <: AbstractSpectrum{W, F} - wave::AbstractArray{W, N} - flux::AbstractArray{F, M} - meta::Dict{Symbol,Any} -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 - -""" - meta(::AbstractSpectrum) - -Return the meta of the spectrum. -""" -meta(spec::AbstractSpectrum) = spec.meta - -Spectrum(wave, flux, meta::Dict{Symbol, Any}) = Spectrum(collect(wave), collect(flux), 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 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") -SingleSpectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}) - wave: (10000.0 Å, 30000.0 Å) - flux: (99999.8952204731 W Å^-1 m^-2, 299999.8866277076 W Å^-1 m^-2) - meta: Dict{Symbol, Any}() - -julia> ustrip(spec) -SingleSpectrum(Float64, Float64) - wave: (10000.0, 30000.0) - flux: (99999.8952204731, 299999.8866277076) - 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 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(wave(spec))), unit(eltype(flux(spec))) From aef38db27ed7a0dcaf4b9ecb939ff51c2b14113c Mon Sep 17 00:00:00 2001 From: icweaver Date: Tue, 4 Nov 2025 09:23:40 -0800 Subject: [PATCH 12/31] use 3D flux for IFUSpectrum --- src/Spectra.jl | 4 ++-- src/spectrum_echelle.jl | 4 ++++ src/spectrum_ifu.jl | 38 +++++++++++++------------------------- src/spectrum_single.jl | 8 +++++--- 4 files changed, 24 insertions(+), 30 deletions(-) diff --git a/src/Spectra.jl b/src/Spectra.jl index bceb336..afa35a7 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -225,8 +225,8 @@ function spectrum(wave::AbstractVector{<:Real}, flux::AbstractVector{<:Real}; kw Spectrum(wave, flux, Dict{Symbol,Any}(kwds)) end -function spectrum(wave::AbstractVector{<:Real}, flux::AbstractMatrix{<:Real}; kwds...) - @assert size(wave, 1) == size(flux, 2) "wave and flux in each order must have equal size" +function spectrum(wave::AbstractVector{<:Real}, flux::AbstractArray{<:Real, 3}; kwds...) + @assert size(wave, 1) == size(flux, 1) "wave and flux in each order must have equal size" Spectrum(wave, flux, Dict{Symbol,Any}(kwds)) end diff --git a/src/spectrum_echelle.jl b/src/spectrum_echelle.jl index 999faab..c234d6b 100644 --- a/src/spectrum_echelle.jl +++ b/src/spectrum_echelle.jl @@ -18,7 +18,11 @@ 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, 1))") + 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 index 567448c..431baf7 100644 --- a/src/spectrum_ifu.jl +++ b/src/spectrum_ifu.jl @@ -1,33 +1,21 @@ -const IFUSpectrum = Spectrum{W, F, 1, 2} where {W, F} +const IFUSpectrum = Spectrum{W, F, 1, 3} where {W, F} -function Base.getindex(spec::IFUSpectrum, i::Int) - w = wave(spec) - f = flux(spec)[i, :] - m = merge(Dict(:Order => i), meta(spec)) - return Spectrum(w, f, m) -end - -function Base.getindex(spec::IFUSpectrum, i::Int, J::AbstractVector) - w = wave(spec) - f = flux(spec)[i, J] - m = merge(Dict(:Order => i), meta(spec)) - return Spectrum(w, f, m) -end +Base.getindex(spec::IFUSpectrum, i::Int, j, k) = flux(spec)[i, j, k] -function Base.getindex(spec::IFUSpectrum, 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) +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), 1) -Base.firstindex(spec::IFUSpectrum, dim::Int) = firstindex(flux(spec), dim) -Base.lastindex(spec::IFUSpectrum) = lastindex(flux(spec), 1) +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) - println(io, "IFUSpectrum($(eltype(wave(spec))), $(eltype(flux(spec))))") - println(io, " # orders: $(size(spec, 1))") - println(io, " wave: ", (extrema∘wave)(spec)) + 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 index 1558da9..5026e30 100644 --- a/src/spectrum_single.jl +++ b/src/spectrum_single.jl @@ -15,8 +15,10 @@ Base.firstindex(spec::SingleSpectrum) = firstindex(wave(spec)) Base.lastindex(spec::SingleSpectrum) = lastindex(wave(spec)) function Base.show(io::IO, spec::SingleSpectrum) - println(io, "SingleSpectrum($(eltype(wave(spec))), $(eltype(flux(spec))))") - println(io, " wave: ", (extrema∘wave)(spec)) - println(io, " flux: ", (extrema∘flux)(spec)) + 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 From 6ec0bf1fdacac61f614197f7858c216a586213f0 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Wed, 5 Nov 2025 13:55:15 -0800 Subject: [PATCH 13/31] get some color in there --- .github/workflows/CI.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index eddb8a5..dde27f2 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -29,6 +29,8 @@ jobs: - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 + - with: + julia_args: "--color=yes" - uses: julia-actions/julia-processcoverage@v1 - uses: codecov/codecov-action@v5 with: From c5b3e1f990b501c227addae3444b4294f9e17351 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Wed, 5 Nov 2025 14:13:17 -0800 Subject: [PATCH 14/31] typo --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index dde27f2..8bf184e 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -29,7 +29,7 @@ jobs: - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - - with: + with: julia_args: "--color=yes" - uses: julia-actions/julia-processcoverage@v1 - uses: codecov/codecov-action@v5 From 4fc12bab8f4eacb9ede5e18f143336f89251ee66 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Wed, 5 Nov 2025 14:55:46 -0800 Subject: [PATCH 15/31] nvm --- .github/workflows/CI.yml | 2 -- 1 file changed, 2 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 8bf184e..eddb8a5 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -29,8 +29,6 @@ jobs: - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - with: - julia_args: "--color=yes" - uses: julia-actions/julia-processcoverage@v1 - uses: codecov/codecov-action@v5 with: From 1f0d9b3803684278714dd3056a039bc297c65860 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Wed, 5 Nov 2025 15:41:06 -0800 Subject: [PATCH 16/31] update docstrings --- docs/src/index.md | 8 ++++---- docs/src/transforms.md | 12 ++++++------ src/Spectra.jl | 22 ++++++++++++---------- src/utils.jl | 8 ++++---- 4 files changed, 26 insertions(+), 24 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 21a28a3..a04374f 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -42,8 +42,8 @@ julia> flux = (read(f[2], "flux") .* 1e-17)u"erg/s/cm^2/angstrom"; julia> spec = spectrum(wave, flux) SingleSpectrum(Quantity{Float32, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, erg, cm^-2, s^-1), 𝐌 𝐋^-1 𝐓^-3, nothing}}) - wave: (3815.0483f0 Å, 9206.613f0 Å) - flux: (9.416705322265625e-16 erg Å^-1 cm^-2 s^-1, 9.156107177734375e-15 erg Å^-1 cm^-2 s^-1) + 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); @@ -54,8 +54,8 @@ julia> plot(spec); ```jldoctest guide julia> cont_fit = continuum(spec) SingleSpectrum(Quantity{Float32, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, erg, cm^-2, s^-1), 𝐌 𝐋^-1 𝐓^-3, nothing}}) - wave: (3815.0483f0 Å, 9206.613f0 Å) - flux: (0.4582525451059156 erg Å^-1 cm^-2 s^-1, 4.402153019328562 erg Å^-1 cm^-2 s^-1) + 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/transforms.md b/docs/src/transforms.md index d33df7e..9da8fb5 100644 --- a/docs/src/transforms.md +++ b/docs/src/transforms.md @@ -37,14 +37,14 @@ julia> flux = (100 .± sigma)u"W/m^2/μm" julia> spec = spectrum(wave, flux) SingleSpectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(μm,), 𝐋, nothing}}, Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(μm^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}) - wave: (1.0 μm, 3.0 μm) - flux: (100.0 ± 0.94 W μm^-1 m^-2, 100.0 ± -1.2 W μm^-1 m^-2) + 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) SingleSpectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(μm,), 𝐋, nothing}}, Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(μm^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}) - wave: (1.0 μm, 3.0 μm) - flux: (89.44 ± 0.84 W μm^-1 m^-2, 98.1 ± 1.2 W μm^-1 m^-2) + 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 @@ -57,8 +57,8 @@ julia> red.flux julia> deredden!(red, 0.3) SingleSpectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(μm,), 𝐋, nothing}}, Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(μm^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}) - wave: (1.0 μm, 3.0 μm) - flux: (100.0 ± 0.94 W μm^-1 m^-2, 100.0 ± 1.2 W μm^-1 m^-2) + 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 diff --git a/src/Spectra.jl b/src/Spectra.jl index afa35a7..a9cfcb7 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -122,14 +122,14 @@ julia> flux = wave .* 10 .+ randn(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: (10000.0 Å, 30000.0 Å) - flux: (99999.8952204731 W Å^-1 m^-2, 299999.8866277076 W Å^-1 m^-2) + wave ((1000,)): 10000.0 Å .. 30000.0 Å + flux ((1000,)): 99999.8952204731 W Å^-1 m^-2 .. 299999.8866277076 W Å^-1 m^-2 meta: Dict{Symbol, Any}() julia> ustrip(spec) SingleSpectrum(Float64, Float64) - wave: (10000.0, 30000.0) - flux: (99999.8952204731, 299999.8866277076) + wave ((1000,)): 10000.0 .. 30000.0 + flux ((1000,)): 99999.8952204731 .. 299999.8866277076 meta: Dict{Symbol, Any}() ``` """ @@ -174,14 +174,14 @@ julia> flux = 100 .* ones(size(wave)); julia> spec = spectrum(wave, flux) SingleSpectrum(Float64, Float64) - wave: (10000.0, 40000.0) - flux: (100.0, 100.0) + wave ((1000,)): 10000.0 .. 40000.0 + flux ((1000,)): 100.0 .. 100.0 meta: Dict{Symbol, Any}() julia> spec = spectrum(wave, flux, name="Just Noise") SingleSpectrum(Float64, Float64) - wave: (10000.0, 40000.0) - flux: (100.0, 100.0) + wave ((1000,)): 10000.0 .. 40000.0 + flux ((1000,)): 100.0 .. 100.0 meta: Dict{Symbol, Any}(:name => "Just Noise") julia> spec.name @@ -202,8 +202,8 @@ julia> flux = (100 .± sigma)u"erg/cm^2/s/angstrom"; julia> spec = spectrum(wave, flux) SingleSpectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(μm,), 𝐋, nothing}}, Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, erg, cm^-2, s^-1), 𝐌 𝐋^-1 𝐓^-3, nothing}}) - wave: (1.0 μm, 4.0 μm) - flux: (100.0 ± 1.2 erg Å^-1 cm^-2 s^-1, 100.0 ± 1.1 erg Å^-1 cm^-2 s^-1) + wave ((1000,)): 1.0 μm .. 4.0 μm + flux ((1000,)): 100.0 ± 1.2 erg Å^-1 cm^-2 s^-1 .. 100.0 ± 1.1 erg Å^-1 cm^-2 s^-1 meta: Dict{Symbol, Any}() ``` @@ -217,6 +217,8 @@ julia> flux = ones(10, 100) .* collect(1:10); julia> spec = spectrum(wave, flux) EchelleSpectrum(Float64, Float64) # orders: 10 + wave ((10, 100)): 100.0 .. 10000.0 + flux ((10, 100)): 1.0 .. 10.0 meta: Dict{Symbol, Any}() ``` """ diff --git a/src/utils.jl b/src/utils.jl index e5ad6d5..2a1dd1e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -22,14 +22,14 @@ julia> wave = range(1, 3, length=100)u"μm" julia> bb = blackbody(wave, 2000u"K") SingleSpectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(μm,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(μm^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}) - wave: (1.0 μm, 3.0 μm) - flux: (49010.54557924032 W μm^-1 m^-2, 131058.60552995963 W μm^-1 m^-2) + 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) SingleSpectrum(Float64, Float64) - wave: (10000.0, 30000.0) - flux: (40.04325690910415, 1190.9562575755397) + 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)] From 0174dd4a1e6e1e96dca85fc8288b96ec07dccf15 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Wed, 5 Nov 2025 16:04:21 -0800 Subject: [PATCH 17/31] update sprints --- test/spectrum.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/test/spectrum.jl b/test/spectrum.jl index ef95e26..07382f2 100644 --- a/test/spectrum.jl +++ b/test/spectrum.jl @@ -34,8 +34,8 @@ Random.seed!(8675309) @test_throws AssertionError spectrum(wave, flux_trimmed) expected = """ SingleSpectrum(Float64, Measurements.Measurement{Float64}) - wave: (10000.0, 50000.0) - flux: (1.0 ± 0.1, 1000.0 ± 1.0) + 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" @@ -87,6 +87,8 @@ end expected = """ EchelleSpectrum(Float64, Measurements.Measurement{Float64}) # orders: 3 + wave ((3, 1000)): 10000.0 .. 50000.0 + flux ((3, 1000)): 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" @@ -131,8 +133,8 @@ end @test strip_spec.meta == spec.meta expected = """ SingleSpectrum(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Unitful.Quantity{Measurements.Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}) - wave: (10000.0 Å, 50000.0 Å) - flux: (1.0 ± 0.1 W Å^-1 m^-2, 1000.0 ± 1.0 W Å^-1 m^-2) + 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 @@ -183,10 +185,11 @@ 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 + wave ((3, 1000)): 10000.0 Å .. 50000.0 Å + flux ((3, 1000)): 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 From 9327cb4860a92876d4d971b27f5f3e10c4c92d7c Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Wed, 5 Nov 2025 16:18:39 -0800 Subject: [PATCH 18/31] reduce number of parentheses used for display --- docs/src/index.md | 8 ++++---- docs/src/transforms.md | 12 ++++++------ src/Spectra.jl | 24 ++++++++++++------------ src/spectrum_echelle.jl | 4 ++-- src/spectrum_ifu.jl | 4 ++-- src/spectrum_single.jl | 4 ++-- src/utils.jl | 8 ++++---- test/spectrum.jl | 16 ++++++++-------- 8 files changed, 40 insertions(+), 40 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index a04374f..673683a 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -42,8 +42,8 @@ julia> flux = (read(f[2], "flux") .* 1e-17)u"erg/s/cm^2/angstrom"; julia> spec = spectrum(wave, flux) 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 + 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); @@ -54,8 +54,8 @@ julia> plot(spec); ```jldoctest guide julia> cont_fit = continuum(spec) 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 + 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/transforms.md b/docs/src/transforms.md index 9da8fb5..89b48be 100644 --- a/docs/src/transforms.md +++ b/docs/src/transforms.md @@ -37,14 +37,14 @@ julia> flux = (100 .± sigma)u"W/m^2/μm" julia> spec = spectrum(wave, flux) 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 + 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) 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 + 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 @@ -57,8 +57,8 @@ julia> red.flux julia> deredden!(red, 0.3) 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 + 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 diff --git a/src/Spectra.jl b/src/Spectra.jl index a9cfcb7..9f2bdf8 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -122,14 +122,14 @@ julia> flux = wave .* 10 .+ randn(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.8952204731 W Å^-1 m^-2 .. 299999.8866277076 W Å^-1 m^-2 + wave (1000,): 10000.0 Å .. 30000.0 Å + flux (1000,): 99999.8952204731 W Å^-1 m^-2 .. 299999.8866277076 W Å^-1 m^-2 meta: Dict{Symbol, Any}() julia> ustrip(spec) SingleSpectrum(Float64, Float64) - wave ((1000,)): 10000.0 .. 30000.0 - flux ((1000,)): 99999.8952204731 .. 299999.8866277076 + wave (1000,): 10000.0 .. 30000.0 + flux (1000,): 99999.8952204731 .. 299999.8866277076 meta: Dict{Symbol, Any}() ``` """ @@ -174,14 +174,14 @@ julia> flux = 100 .* ones(size(wave)); julia> spec = spectrum(wave, flux) SingleSpectrum(Float64, Float64) - wave ((1000,)): 10000.0 .. 40000.0 - flux ((1000,)): 100.0 .. 100.0 + wave (1000,): 10000.0 .. 40000.0 + flux (1000,): 100.0 .. 100.0 meta: Dict{Symbol, Any}() julia> spec = spectrum(wave, flux, name="Just Noise") SingleSpectrum(Float64, Float64) - wave ((1000,)): 10000.0 .. 40000.0 - flux ((1000,)): 100.0 .. 100.0 + wave (1000,): 10000.0 .. 40000.0 + flux (1000,): 100.0 .. 100.0 meta: Dict{Symbol, Any}(:name => "Just Noise") julia> spec.name @@ -202,8 +202,8 @@ julia> flux = (100 .± sigma)u"erg/cm^2/s/angstrom"; julia> spec = spectrum(wave, flux) 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 ± 1.2 erg Å^-1 cm^-2 s^-1 .. 100.0 ± 1.1 erg Å^-1 cm^-2 s^-1 + wave (1000,): 1.0 μm .. 4.0 μm + flux (1000,): 100.0 ± 1.2 erg Å^-1 cm^-2 s^-1 .. 100.0 ± 1.1 erg Å^-1 cm^-2 s^-1 meta: Dict{Symbol, Any}() ``` @@ -217,8 +217,8 @@ julia> flux = ones(10, 100) .* collect(1:10); julia> spec = spectrum(wave, flux) EchelleSpectrum(Float64, Float64) # orders: 10 - wave ((10, 100)): 100.0 .. 10000.0 - flux ((10, 100)): 1.0 .. 10.0 + wave (10, 100): 100.0 .. 10000.0 + flux (10, 100): 1.0 .. 10.0 meta: Dict{Symbol, Any}() ``` """ diff --git a/src/spectrum_echelle.jl b/src/spectrum_echelle.jl index c234d6b..8dbb8d2 100644 --- a/src/spectrum_echelle.jl +++ b/src/spectrum_echelle.jl @@ -22,7 +22,7 @@ function Base.show(io::IO, spec::EchelleSpectrum) f = flux(spec) println(io, "EchelleSpectrum($(eltype(wave(spec))), $(eltype(flux(spec))))") println(io, " # orders: $(size(spec, 1))") - println(io, " wave ($(size(w))): ", first(w), " .. ", last(w)) - println(io, " flux ($(size(f))): ", first(f), " .. ", last(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_ifu.jl b/src/spectrum_ifu.jl index 431baf7..91e2545 100644 --- a/src/spectrum_ifu.jl +++ b/src/spectrum_ifu.jl @@ -15,7 +15,7 @@ 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)) + 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 index 5026e30..6daa9f7 100644 --- a/src/spectrum_single.jl +++ b/src/spectrum_single.jl @@ -18,7 +18,7 @@ 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)) + 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 2a1dd1e..a6b9b59 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -22,14 +22,14 @@ julia> wave = range(1, 3, length=100)u"μm" julia> bb = blackbody(wave, 2000u"K") 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 + 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) SingleSpectrum(Float64, Float64) - wave ((100,)): 10000.0 .. 30000.0 - flux ((100,)): 1190.9562575755397 .. 40.04325690910415 + 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)] diff --git a/test/spectrum.jl b/test/spectrum.jl index 07382f2..a88a1e5 100644 --- a/test/spectrum.jl +++ b/test/spectrum.jl @@ -34,8 +34,8 @@ Random.seed!(8675309) @test_throws AssertionError spectrum(wave, flux_trimmed) expected = """ SingleSpectrum(Float64, Measurements.Measurement{Float64}) - wave ((1000,)): 10000.0 .. 50000.0 - flux ((1000,)): 100.0 ± -2.8 .. 100.0 ± 0.6 + 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" @@ -87,8 +87,8 @@ end expected = """ EchelleSpectrum(Float64, Measurements.Measurement{Float64}) # orders: 3 - wave ((3, 1000)): 10000.0 .. 50000.0 - flux ((3, 1000)): 100.0 ± -2.8 .. 100.0 ± 0.6 + wave (3, 1000): 10000.0 .. 50000.0 + flux (3, 1000): 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" @@ -133,8 +133,8 @@ end @test strip_spec.meta == spec.meta expected = """ 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 + 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 @@ -188,8 +188,8 @@ end 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 - wave ((3, 1000)): 10000.0 Å .. 50000.0 Å - flux ((3, 1000)): 100.0 ± -2.8 W Å^-1 m^-2 .. 100.0 ± 0.6 W Å^-1 m^-2 + wave (3, 1000): 10000.0 Å .. 50000.0 Å + flux (3, 1000): 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 From d355ea46f36c3a38bd47dec2781eaeb79a4039b2 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Thu, 6 Nov 2025 18:44:19 -0800 Subject: [PATCH 19/31] moar tests --- test/runtests.jl | 2 +- test/spectrum.jl | 25 +++++++++++++++++++++++++ 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 9d9b503..ed6cbdd 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, 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 a88a1e5..6d73314 100644 --- a/test/spectrum.jl +++ b/test/spectrum.jl @@ -94,6 +94,31 @@ end @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)) From 3fc80e9c542acc488badf78b10a9a21d74192a71 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Thu, 6 Nov 2025 19:28:36 -0800 Subject: [PATCH 20/31] more test coverage --- src/Spectra.jl | 5 +++-- test/runtests.jl | 2 +- test/spectrum.jl | 1 + 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/Spectra.jl b/src/Spectra.jl index 9f2bdf8..2a817ce 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -38,6 +38,9 @@ mutable struct Spectrum{W<:Number, F<:Number, N, M} <: AbstractSpectrum{W, F} meta::Dict{Symbol,Any} end +# Doesn't seem to be used atp +#Spectrum(wave, flux, meta::Dict{Symbol, Any}) = Spectrum(collect(wave), collect(flux), meta) + """ wave(::AbstractSpectrum) @@ -59,8 +62,6 @@ Return the meta of the spectrum. """ meta(spec::AbstractSpectrum) = spec.meta -Spectrum(wave, flux, meta::Dict{Symbol, Any}) = Spectrum(collect(wave), collect(flux), meta) - function Base.getproperty(spec::AbstractSpectrum, nm::Symbol) if nm in keys(getfield(spec, :meta)) return getfield(spec, :meta)[nm] diff --git a/test/runtests.jl b/test/runtests.jl index ed6cbdd..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, SingleSpectrum, EchelleSpectrum, IFUSpectrum, 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 6d73314..7859243 100644 --- a/test/spectrum.jl +++ b/test/spectrum.jl @@ -15,6 +15,7 @@ Random.seed!(8675309) @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 From 0a93d838e4d9afacd395e313728929a18e76fe29 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Thu, 6 Nov 2025 22:39:24 -0800 Subject: [PATCH 21/31] add size and monotonicity checks, switch from row --> col major representation --- src/Spectra.jl | 34 +++++++++++++++++++++++++++++++--- test/utils.jl | 4 ++-- 2 files changed, 33 insertions(+), 5 deletions(-) diff --git a/src/Spectra.jl b/src/Spectra.jl index 2a817ce..06372d7 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -32,10 +32,38 @@ abstract type AbstractSpectrum{W,F} end 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, N, M} <: AbstractSpectrum{W, F} - wave::AbstractArray{W, N} - flux::AbstractArray{F, M} +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 diff --git a/test/utils.jl b/test/utils.jl index ee217c2..27ddc62 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -19,8 +19,8 @@ using Spectra: spectrum, blackbody, line_flux, equivalent_width end @testset "Line flux" begin - spec = spectrum([3, 5, 4], [6, 7, 8]) - @test line_flux(spec) == 6 + spec = spectrum([3, 4, 5], [6, 7, 8]) + @test line_flux(spec) == 15 end @testset "Equivalent width" begin From 86ec55e2861031f6cd853216204f3e6fc8077ae8 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Thu, 6 Nov 2025 23:29:11 -0800 Subject: [PATCH 22/31] remove old asserts, finish updating tests --- src/Spectra.jl | 5 --- src/spectrum_echelle.jl | 10 +++--- test/spectrum.jl | 72 ++++++++++++++++++++++++++--------------- 3 files changed, 51 insertions(+), 36 deletions(-) diff --git a/src/Spectra.jl b/src/Spectra.jl index 06372d7..657e5d8 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -252,28 +252,23 @@ EchelleSpectrum(Float64, Float64) ``` """ 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{<:Real}, flux::AbstractArray{<:Real, 3}; kwds...) - @assert size(wave, 1) == size(flux, 1) "wave and flux in each order must have equal size" 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" 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" 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" Spectrum(wave, flux, Dict{Symbol,Any}(kwds)) end diff --git a/src/spectrum_echelle.jl b/src/spectrum_echelle.jl index 8dbb8d2..fe943da 100644 --- a/src/spectrum_echelle.jl +++ b/src/spectrum_echelle.jl @@ -1,15 +1,15 @@ 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, :] + 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, :] + w = wave(spec)[:, I] + f = flux(spec)[:, I] m = merge(Dict(:Orders => (first(I), last(I))), meta(spec)) return Spectrum(w, f, m) end @@ -21,7 +21,7 @@ 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, 1))") + 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)) diff --git a/test/spectrum.jl b/test/spectrum.jl index 7859243..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)) @@ -32,7 +52,7 @@ Random.seed!(8675309) @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 = """ SingleSpectrum(Float64, Measurements.Measurement{Float64}) wave (1000,): 10000.0 .. 50000.0 @@ -43,17 +63,17 @@ Random.seed!(8675309) 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") @@ -73,23 +93,23 @@ 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 - wave (3, 1000): 10000.0 .. 50000.0 - flux (3, 1000): 100.0 ± -2.8 .. 100.0 ± 0.6 + 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" @@ -166,17 +186,17 @@ end 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" @@ -192,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 @@ -214,8 +234,8 @@ end 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 - wave (3, 1000): 10000.0 Å .. 50000.0 Å - flux (3, 1000): 100.0 ± -2.8 W Å^-1 m^-2 .. 100.0 ± 0.6 W Å^-1 m^-2 + 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 From cb50530147a821172f8587d823e77fa2b2d041cd Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Fri, 7 Nov 2025 02:42:13 -0800 Subject: [PATCH 23/31] add some dev docs --- README.md | 30 ++++++++++++++++++++++++++++++ docs/Project.toml | 1 + docs/make.jl | 3 +++ 3 files changed, 34 insertions(+) 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 `) +``` + +### 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 bb712bc..e0cf339 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -2,6 +2,9 @@ using Documenter using Spectra using Unitful using Measurements +using Revise + +Revise.revise() DocMeta.setdocmeta!(Spectra, :DocTestSetup, :(using Spectra); recursive = true) From f8dc7b90c82785d474e1b74e865b43c597898081 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Fri, 7 Nov 2025 03:18:34 -0800 Subject: [PATCH 24/31] docs cleanup + more examples --- docs/src/spectrum.md | 14 +++------- docs/src/transforms.md | 24 +++-------------- src/Spectra.jl | 49 ++++++++++++++++++++++----------- src/spectrum_echelle.jl | 52 +++++++++++++++++++++++++++++++++++ src/spectrum_ifu.jl | 60 +++++++++++++++++++++++++++++++++++++++++ src/spectrum_single.jl | 19 +++++++++++++ 6 files changed, 170 insertions(+), 48 deletions(-) diff --git a/docs/src/spectrum.md b/docs/src/spectrum.md index 03431b4..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. @@ -16,6 +9,9 @@ Spectra are defined as possible subtypes of `AbstractSpectrum`. You can use thes ```@docs Spectra.AbstractSpectrum Spectra.Spectrum +Spectra.SingleSpectrum +Spectra.EchelleSpectrum +Spectra.IFUSpectrum ``` ## Constructors @@ -89,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 89b48be..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}}}: @@ -73,7 +59,3 @@ redden! deredden deredden! ``` - -```@meta -DocTestSetup = nothing -``` diff --git a/src/Spectra.jl b/src/Spectra.jl index 657e5d8..bc7665b 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -1,7 +1,7 @@ module Spectra # Uniform API -export AbstractSpectrum, spectrum +export AbstractSpectrum, Spectrum, spectrum # spectra_single.jl, spectra_ifu.jl, spectra_echelle.jl export SingleSpectrum, IFUSpectrum, EchelleSpectrum # utils.jl @@ -19,11 +19,13 @@ using PhysicalConstants.CODATA2018: h, c_0, k_B """ AbstractSpectrum{W<:Number, F<:Number} -An abstract holder for astronomical spectra. All types inheriting from this must have the following fields +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} +* `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 @@ -143,22 +145,27 @@ Remove the units from a spectrum. Useful for processing spectra in tools that do # 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(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.8952204731 W Å^-1 m^-2 .. 299999.8866277076 W Å^-1 m^-2 + 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.8952204731 .. 299999.8866277076 + flux (1000,): 99999.76809093042 .. 300000.2474309158 meta: Dict{Symbol, Any}() ``` """ @@ -171,11 +178,16 @@ Get the units of a spectrum. Returns a tuple of the wavelength units and flux/si # 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(1000); +julia> flux = wave .* 10 .+ randn(rng, 1000); julia> spec = spectrum(wave * u"angstrom", flux * u"W/m^2/angstrom"); @@ -221,33 +233,38 @@ 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) 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 ± 1.2 erg Å^-1 cm^-2 s^-1 .. 100.0 ± 1.1 erg Å^-1 cm^-2 s^-1 + 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 (10, 100): 100.0 .. 10000.0 - flux (10, 100): 1.0 .. 10.0 + wave (100, 10): 100.0 .. 10000.0 + flux (100, 10): 1.0 .. 10.0 meta: Dict{Symbol, Any}() ``` """ diff --git a/src/spectrum_echelle.jl b/src/spectrum_echelle.jl index fe943da..91b0e47 100644 --- a/src/spectrum_echelle.jl +++ b/src/spectrum_echelle.jl @@ -1,3 +1,55 @@ +""" + 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) diff --git a/src/spectrum_ifu.jl b/src/spectrum_ifu.jl index 91e2545..c1b90ee 100644 --- a/src/spectrum_ifu.jl +++ b/src/spectrum_ifu.jl @@ -1,5 +1,64 @@ +""" + 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) @@ -8,6 +67,7 @@ function Base.getindex(spec::IFUSpectrum, 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) diff --git a/src/spectrum_single.jl b/src/spectrum_single.jl index 6daa9f7..27ca085 100644 --- a/src/spectrum_single.jl +++ b/src/spectrum_single.jl @@ -1,3 +1,22 @@ +""" + 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)), ) From 216a40e07605a6a064a5db644ab1b5eca6c59615 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Wed, 26 Nov 2025 21:15:08 -0800 Subject: [PATCH 25/31] mention XSpecs in quickstart --- docs/src/index.md | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/docs/src/index.md b/docs/src/index.md index 673683a..cae3e18 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -18,7 +18,9 @@ julia> using Spectra ## Quick Start -Here is a quick demo of some of our features +Here is a quick demo of some of our features. + +### Spectrum construction ```jldoctest guide julia> using Spectra, FITSIO, Unitful, UnitfulAstro, Plots @@ -51,6 +53,10 @@ julia> plot(spec); ![](assets/sdss.svg) +For constructing higher dimensional spectra, e.g., for echelle or IFU spectra, see the docstrings for [EchelleSpectrum](@ref) and [IFUSpectrum](@ref), respectively. + +### Continuum fitting + ```jldoctest guide julia> cont_fit = continuum(spec) SingleSpectrum(Quantity{Float32, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, erg, cm^-2, s^-1), 𝐌 𝐋^-1 𝐓^-3, nothing}}) From 7d8ec3bedc418f4db09d0e4c86914eec74756b23 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Mon, 15 Dec 2025 17:40:29 -0800 Subject: [PATCH 26/31] generalized terms, shelved extra functionality --- src/Spectra.jl | 174 +++++++++++++++++++------------------- src/plotting.jl | 4 +- src/spectrum_echelle.jl | 34 ++++---- src/spectrum_ifu.jl | 42 ++++----- src/spectrum_single.jl | 26 +++--- src/transforms/redden.jl | 5 +- src/utils.jl | 8 +- test/fitting.jl | 14 --- test/old_tests/fitting.jl | 12 +++ test/old_tests/utils.jl | 29 +++++++ test/plotting.jl | 6 +- test/runtests.jl | 2 +- test/spectrum.jl | 158 +++++++++++++++++----------------- test/transforms/redden.jl | 20 ++--- test/utils.jl | 29 ------- 15 files changed, 282 insertions(+), 281 deletions(-) delete mode 100644 test/fitting.jl create mode 100644 test/old_tests/fitting.jl create mode 100644 test/old_tests/utils.jl delete mode 100644 test/utils.jl diff --git a/src/Spectra.jl b/src/Spectra.jl index bc7665b..bf3a53a 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -1,13 +1,13 @@ module Spectra # Uniform API -export AbstractSpectrum, Spectrum, spectrum -# spectra_single.jl, spectra_ifu.jl, spectra_echelle.jl -export SingleSpectrum, IFUSpectrum, EchelleSpectrum +export AbstractSpectrum, Spectrum, spectrum, spectral_axis, flux_axis +# spectra_single.jl, spectra_ifu.jl, spectra_echelle.jl, spectra_binned.jl +export SingleSpectrum, IFUSpectrum, EchelleSpectrum, spectral_axis, flux_axis # utils.jl -export blackbody, line_flux, equivalent_width +#export blackbody, line_flux, equivalent_width # fitting/fitting.jl -export continuum, continuum! +#export continuum, continuum! # transforms/redden.jl export redden, redden!, deredden, deredden! @@ -17,78 +17,81 @@ using Unitful: Unitful, Quantity, @u_str, ustrip, unit, dimension using PhysicalConstants.CODATA2018: h, c_0, k_B """ - AbstractSpectrum{W<:Number, F<:Number} + AbstractSpectrum{S <: 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}` +* `spectral_axis::Array{S, M}` +* `flux_axis::Array{F, N}` * `meta::Dict{Symbol, Any}` See [`SingleSpectrum`](@ref), [`EchelleSpectrum`](@ref), and [`IFUSpectrum`](@ref) for different subtypes. """ -abstract type AbstractSpectrum{W,F} end +abstract type AbstractSpectrum{S, F} end """ Spectrum <: AbstractSpectrum -A spectrum or spectra stored as arrays of real numbers. The wavelengths are assumed to be in angstrom. +A spectrum or spectra stored as arrays of real numbers. For UV/VIS/IR spectra, the `spectral_axis` is assumed to be wavelengths (in angstrom). For X-ray spectra, the `spectral_axis` is assumed to be energies (in keV). """ -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} +mutable struct Spectrum{S<:Number, F<:Number, M, N} <: AbstractSpectrum{S, F} + spectral_axis::AbstractArray{S, M} + flux_axis::AbstractArray{F, N} + meta::Dict{Symbol, Any} + function Spectrum{S, F, M, N}(s, f, meta) where {S<:Number, F<:Number, M, N} + # TODO: Investigate using Holy Traits to help with validation # 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) + if size(s, 1) != size(f, 1) + throw(ArgumentError( + """ + Spectral axis and flux axis are incompatible sizes. 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) + - TODO: BinnedSpectrum (final name(s) tbd): + - energy (M × 2 matrix), flux (N-length vector) + - others? + + See the documentation for each spectrum type for more. + """)) + end + + # Spectral axis monoticity check + spec_ax = eachcol(s) + if !(all(issorted, spec_ax) || all(x -> issorted(x; rev = true), spec_ax)) + throw(ArgumentError("Spectral axis must be strictly increasing or decreasing.")) + end + + return new{S, F, M, N}(s, f, meta) end end -function Spectrum(wave, flux, meta) - Spectrum{eltype(wave), eltype(flux), ndims(wave), ndims(flux)}(wave, flux, meta) +function Spectrum(s, f, meta) + Spectrum{eltype(s), eltype(f), ndims(s), ndims(f)}(s, f, meta) end # Doesn't seem to be used atp #Spectrum(wave, flux, meta::Dict{Symbol, Any}) = Spectrum(collect(wave), collect(flux), meta) """ - wave(::AbstractSpectrum) + spectral_axis(spec::AbstractSpectrum) -Return the wavelengths of the spectrum. +Return the spectral axis of `spec`. """ -wave(spec::AbstractSpectrum) = spec.wave +spectral_axis(spec::AbstractSpectrum) = spec.spectral_axis """ - flux(::AbstractSpectrum) + flux(spec::AbstractSpectrum) -Return the flux of the spectrum. +Return the flux axis of `spec`. """ -flux(spec::AbstractSpectrum) = spec.flux +flux_axis(spec::AbstractSpectrum) = spec.flux_axis """ - meta(::AbstractSpectrum) + meta(spec::AbstractSpectrum) -Return the meta of the spectrum. +Return the meta data associated with `spec`. """ meta(spec::AbstractSpectrum) = spec.meta @@ -101,42 +104,42 @@ function Base.getproperty(spec::AbstractSpectrum, nm::Symbol) end function Base.propertynames(spec::AbstractSpectrum) - natural = (:wave, :flux, :meta) + natural = (:spectral_axis, :flux_axis, :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) +Base.argmax(spec::AbstractSpectrum) = argmax(flux_axis(spec)) +Base.argmin(spec::AbstractSpectrum) = argmin(flux_axis(spec)) +Base.eltype(spec::AbstractSpectrum) = eltype(flux_axis(spec)) +Base.findmax(spec::AbstractSpectrum) = findmax(flux_axis(spec)) +Base.findmin(spec::AbstractSpectrum) = findmin(flux_axis(spec)) +Base.length(spec::AbstractSpectrum) = length(flux_axis(spec)) +Base.maximum(spec::AbstractSpectrum) = maximum(flux_axis(spec)) +Base.minimum(spec::AbstractSpectrum) = minimum(flux_axis(spec)) +Base.size(spec::AbstractSpectrum) = size(flux_axis(spec)) +Base.size(spec::AbstractSpectrum, i) = size(flux_axis(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, o::AbstractSpectrum) = spectral_axis(s) == spectral_axis(o) && flux_axis(s) == flux_axis(o) && meta(s) == meta(o) +Base.:+(s::T, A) where {T <: AbstractSpectrum} = T(spectral_axis(s), flux_axis(s) .+ A, meta(s)) +Base.:*(s::T, A::Union{Real, AbstractVector}) where {T <: AbstractSpectrum} = T(spectral_axis(s), flux_axis(s) .* A, meta(s)) +Base.:/(s::T, A) where {T <: AbstractSpectrum} = T(spectral_axis(s), flux_axis(s) ./ A, meta(s)) +Base.:-(s::T) where {T <: AbstractSpectrum} = T(spectral_axis(s), -flux_axis(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)) +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)) +Base.:-(s::T, o::T) where {T <: AbstractSpectrum} = T(spectral_axis(s), flux_axis(s) .- flux_axis(o), meta(s)) """ Unitful.ustrip(::AbstractSpectrum) @@ -169,12 +172,12 @@ SingleSpectrum(Float64, Float64) meta: Dict{Symbol, Any}() ``` """ -Unitful.ustrip(spec::AbstractSpectrum) = spectrum(ustrip.(wave(spec)), ustrip.(flux(spec)); meta(spec)...) +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 wavelength units and flux/sigma units +Get the units of a spectrum. Returns a tuple of the spectral axis units and flux/sigma units # Examples ```jldoctest @@ -195,17 +198,18 @@ julia> w_unit, f_unit = unit(spec) (Å, W Å^-1 m^-2) ``` """ -Unitful.unit(spec::AbstractSpectrum) = unit(eltype(wave(spec))), unit(eltype(flux(spec))) +Unitful.unit(spec::AbstractSpectrum) = unit(eltype(spectral_axis(spec))), unit(eltype(flux_axis(spec))) # Spectrum types and basic arithmetic include("spectrum_single.jl") -include("spectrum_ifu.jl") include("spectrum_echelle.jl") +include("spectrum_ifu.jl") +#include("spectrum_binned.jl") """ - spectrum(wave, flux; kwds...) + spectrum(spectral_axis, flux_axis, [meta]) -Construct a spectrum given the spectral wavelengths and fluxes. This will automatically dispatch the correct spectrum type given the shape and element type of the given flux. Any keyword arguments will be accessible from the spectrum as properties. +Construct a spectrum given the spectral axis and flux axis. This will automatically dispatch the correct spectrum type given the shape and element type of the given flux. Any keyword arguments will be accessible from the spectrum as properties. # Examples ```jldoctest @@ -268,32 +272,32 @@ EchelleSpectrum(Float64, Float64) meta: Dict{Symbol, Any}() ``` """ -function spectrum(wave::AbstractVector{<:Real}, flux::AbstractVector{<:Real}; kwds...) - Spectrum(wave, flux, Dict{Symbol,Any}(kwds)) +function spectrum(spectral_axis::AbstractVector{<:Real}, flux_axis::AbstractVector{<:Real}; kwds...) + Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds)) end -function spectrum(wave::AbstractVector{<:Real}, flux::AbstractArray{<:Real, 3}; kwds...) - Spectrum(wave, flux, Dict{Symbol,Any}(kwds)) +function spectrum(spectral_axis::AbstractVector{<:Real}, flux_axis::AbstractArray{<:Real, 3}; kwds...) + Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds)) end -function spectrum(wave::AbstractMatrix{<:Real}, flux::AbstractMatrix{<:Real}; kwds...) - Spectrum(wave, flux, Dict{Symbol,Any}(kwds)) +function spectrum(spectral_axis::AbstractMatrix{<:Real}, flux_axis::AbstractMatrix{<:Real}; kwds...) + Spectrum(spectral_axis, flux_axis, 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)) +function spectrum(spectral_axis::AbstractVector{<:Quantity}, flux_axis::AbstractVector{<: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(wave::AbstractMatrix{<:Quantity}, flux::AbstractMatrix{<:Quantity}; kwds...) - @assert dimension(eltype(wave)) == u"𝐋" "wave not recognized as having dimensions of wavelengths" - Spectrum(wave, flux, Dict{Symbol,Any}(kwds)) +function spectrum(spectral_axis::AbstractMatrix{<:Quantity}, flux_axis::AbstractMatrix{<: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 # tools include("utils.jl") include("transforms/transforms.jl") include("plotting.jl") -include("fitting/fitting.jl") +#include("fitting/fitting.jl") end # module diff --git a/src/plotting.jl b/src/plotting.jl index 427a704..9fef960 100644 --- a/src/plotting.jl +++ b/src/plotting.jl @@ -3,7 +3,7 @@ xlabel --> "wave" ylabel --> "flux density" label --> "" - ustrip.(spec.wave), Measurements.value.(ustrip.(spec.flux)) + ustrip.(spectral_axis(spec)), Measurements.value.(ustrip.(flux_axis(spec))) end @recipe function f(spec::EchelleSpectrum) @@ -11,5 +11,5 @@ end xlabel --> "wave" ylabel --> "flux density" label --> ["Order $i" for i in 1:size(spec, 1)] - spec.wave', spec.flux' + spectral_axis(spec)', flux_axis(spec)' end diff --git a/src/spectrum_echelle.jl b/src/spectrum_echelle.jl index 91b0e47..e3c4873 100644 --- a/src/spectrum_echelle.jl +++ b/src/spectrum_echelle.jl @@ -1,9 +1,9 @@ """ EchelleSpectrum <: AbstractSpectrum -An instance of [`Spectrum`](@ref) where the wavelength and flux are both 2D arrays, i.e., ``M = N = 2``. +An instance of [`Spectrum`](@ref) where the spectral and flux axes 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). +The spectral and flux matrices are both ``m`` rows in wavelength by ``n`` columns in [echelle order](https://en.wikipedia.org/wiki/Echelle_grating). # Examples @@ -37,14 +37,14 @@ julia> flux = repeat(1:4, 1, 10)' julia> spec = Spectrum(wave, flux, Dict()) EchelleSpectrum(Int64, Int64) # orders: 4 - wave (10, 4): 1 .. 40 - flux (10, 4): 1 .. 4 + spectral axis (10, 4): 1 .. 40 + flux axis (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 + spectral axis (10,): 1 .. 10 + flux axis (10,): 1 .. 1 meta: Dict{Symbol, Any}(:Order => 1) ``` @@ -53,28 +53,28 @@ See [`SingleSpectrum`](@ref) for a 1D variant, and [`IFUSpectrum`](@ref) for a 3 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] + w = spectral_axis(spec)[:, i] + f = flux_axis(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] + w = spectral_axis(spec)[:, I] + f = flux_axis(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) +Base.firstindex(spec::EchelleSpectrum) = firstindex(flux_axis(spec), 1) +Base.lastindex(spec::EchelleSpectrum) = lastindex(flux_axis(spec), 1) function Base.show(io::IO, spec::EchelleSpectrum) - w = wave(spec) - f = flux(spec) - println(io, "EchelleSpectrum($(eltype(wave(spec))), $(eltype(flux(spec))))") + w = spectral_axis(spec) + f = flux_axis(spec) + println(io, "EchelleSpectrum($(eltype(spectral_axis(spec))), $(eltype(flux_axis(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)) + println(io, " spectral axis $(size(w)): ", first(w), " .. ", last(w)) + println(io, " flux axis $(size(f)): ", first(f), " .. ", last(f)) print(io, " meta: ", meta(spec)) end diff --git a/src/spectrum_ifu.jl b/src/spectrum_ifu.jl index c1b90ee..4687eee 100644 --- a/src/spectrum_ifu.jl +++ b/src/spectrum_ifu.jl @@ -1,9 +1,9 @@ """ 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``. +An instance of [`Spectrum`](@ref) where the spectral axis is a 1D array and flux axis 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. +The spectral vector is of length ``m`` and flux 3D array has shape ``(m \\times n \\times k)`` where each entry of the flux along the spectral axis is an ``n \\times k`` matrix. # Examples @@ -17,8 +17,8 @@ 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 + spectral axis (5,): 20 .. 200 + flux axis (5, 10, 6): 0.4552384158732863 .. 0.11698905483599475 meta: Dict{Symbol, Any}() julia> spec[begin] # IFU image at first wavelength @@ -36,14 +36,14 @@ julia> spec[begin] # IFU image at first wavelength 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 + spectral axis (3,): 20 .. 120 + flux axis (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 + spectral axis (5,): 20 .. 200 + flux axis (5,): 0.4552384158732863 .. 0.02964765308691042 meta: Dict{Symbol, Any}() ``` @@ -51,31 +51,31 @@ See [`SingleSpectrum`](@ref) for a 1D variant, and [`EchelleSpectrum`](@ref) for """ const IFUSpectrum = Spectrum{W, F, 1, 3} where {W, F} -Base.getindex(spec::IFUSpectrum, i) = flux(spec)[i, :, :] +Base.getindex(spec::IFUSpectrum, i) = flux_axis(spec)[i, :, :] function Base.getindex(spec::IFUSpectrum, I::AbstractVector) - w = wave(spec)[I] - f = flux(spec)[I, :, :] + w = spectral_axis(spec)[I] + f = flux_axis(spec)[I, :, :] return Spectrum(w, f, meta(spec)) end -Base.getindex(spec::IFUSpectrum, i::Int, j, k) = flux(spec)[i, j, k] +Base.getindex(spec::IFUSpectrum, i::Int, j, k) = flux_axis(spec)[i, j, k] function Base.getindex(spec::IFUSpectrum, i, j, k) - w = wave(spec)[i] - f = flux(spec)[i, j, k] + w = spectral_axis(spec)[i] + f = flux_axis(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) +Base.firstindex(spec::IFUSpectrum) = firstindex(flux_axis(spec)) +Base.firstindex(spec::IFUSpectrum, i) = firstindex(flux_axis(spec), i) +Base.lastindex(spec::IFUSpectrum, i) = lastindex(flux_axis(spec), i) function Base.show(io::IO, spec::IFUSpectrum) - w = wave(spec) - f = flux(spec) + w = spectral_axis(spec) + f = flux_axis(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)) + println(io, " spectral axis $(size(w)): ", first(w), " .. ", last(w)) + println(io, " flux axis $(size(f)): ", first(f), " .. ", last(f)) print(io, " meta: ", meta(spec)) end diff --git a/src/spectrum_single.jl b/src/spectrum_single.jl index 27ca085..2615fd0 100644 --- a/src/spectrum_single.jl +++ b/src/spectrum_single.jl @@ -1,17 +1,17 @@ """ SingleSpectrum <: AbstractSpectrum -An instance of [`Spectrum`](@ref) where the wavelength and flux are both 1D arrays, i.e., ``M = N = 1``. +An instance of [`Spectrum`](@ref) where the spectral and flux axes are both 1D arrays, i.e., ``M = N = 1``. -Both wavelength and flux vectors must have the same length, ``m = n``, respectively. +Both spectral and flux axis 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 + spectral axis (4,): 6 .. 9 + flux axis (4,): 2 .. 5 meta: Dict{Symbol, Any}() ``` @@ -19,25 +19,25 @@ See [`EchelleSpectrum`](@ref) and [`IFUSpectrum`](@ref) for working with instanc """ const SingleSpectrum = Spectrum{W, F, 1, 1} where {W, F} -#Base.size(spec::SingleSpectrum) = (length(wave(spec)), ) +#Base.size(spec::SingleSpectrum) = (length(spectral_axis(spec)), ) #Base.IndexStyle(::Type{<:SingleSpectrum}) = IndexLinear() function Base.getindex(spec::SingleSpectrum, i::Int) - return Spectrum([wave(spec)[i]], [flux(spec)[i]], meta(spec)) + return Spectrum([spectral_axis(spec)[i]], [flux_axis(spec)[i]], meta(spec)) end function Base.getindex(spec::SingleSpectrum, inds) - return Spectrum(wave(spec)[inds], flux(spec)[inds], meta(spec)) + return Spectrum(spectral_axis(spec)[inds], flux_axis(spec)[inds], meta(spec)) end -Base.firstindex(spec::SingleSpectrum) = firstindex(wave(spec)) -Base.lastindex(spec::SingleSpectrum) = lastindex(wave(spec)) +Base.firstindex(spec::SingleSpectrum) = firstindex(spectral_axis(spec)) +Base.lastindex(spec::SingleSpectrum) = lastindex(spectral_axis(spec)) function Base.show(io::IO, spec::SingleSpectrum) - w = wave(spec) - f = flux(spec) + w = spectral_axis(spec) + f = flux_axis(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)) + println(io, " spectral axis $(size(w)): ", first(w), " .. ", last(w)) + println(io, " flux axis $(size(f)): ", first(f), " .. ", last(f)) print(io, " meta: ", meta(spec)) end diff --git a/src/transforms/redden.jl b/src/transforms/redden.jl index 9fcf9fe..f73dbe5 100644 --- a/src/transforms/redden.jl +++ b/src/transforms/redden.jl @@ -1,13 +1,12 @@ import DustExtinction - """ redden!(::AbstractSpectrum, Av; Rv = 3.1, law = DustExtinction.CCM89) In-place version of [`redden`](@ref) """ function redden!(spec::T, Av; Rv = 3.1, law = DustExtinction.CCM89) where {T <: AbstractSpectrum} - @. spec.flux = DustExtinction.redden(law, spec.wave, spec.flux; Rv, Av) + @. spec.flux_axis = DustExtinction.redden(law, spec.spectral_axis, spec.flux_axis; Rv, Av) return spec end @@ -31,7 +30,7 @@ end In-place version of [`deredden`](@ref) """ function deredden!(spec::AbstractSpectrum, Av; Rv = 3.1, law = DustExtinction.CCM89) - @. spec.flux = DustExtinction.deredden(law, spec.wave, spec.flux; Rv, Av) + @. spec.flux_axis = DustExtinction.deredden(law, spec.spectral_axis, spec.flux_axis; Rv, Av) return spec end diff --git a/src/utils.jl b/src/utils.jl index a6b9b59..f9df7d4 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -32,7 +32,7 @@ SingleSpectrum(Float64, Float64) flux (100,): 1190.9562575755397 .. 40.04325690910415 meta: Dict{Symbol, Any}(:T => 6000, :name => "Blackbody") -julia> bb.wave[argmax(bb)] +julia> spectral_axis(bb)[argmax(bb)] 1.4444444444444444 μm julia> 2898u"μm*K" / bb.T # See if it matches up with Wien's law @@ -65,7 +65,7 @@ blackbody(T::Quantity) = w->2h * c_0^2 / w^5 / (exp(h * c_0 / (w * k_B * T)) - 1 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] + dx = spectral_axis(spec)[end] - spectral_axis(spec)[1] flux = ustrip(line_flux(spec)) return dx - flux * unit(dx) end @@ -76,6 +76,6 @@ end 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) + avg_dx = diff(spectral_axis(spec)) + return sum(flux_axis(spec)[2:end] .* avg_dx) end diff --git a/test/fitting.jl b/test/fitting.jl deleted file mode 100644 index c09e9cc..0000000 --- a/test/fitting.jl +++ /dev/null @@ -1,14 +0,0 @@ -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/old_tests/fitting.jl b/test/old_tests/fitting.jl new file mode 100644 index 0000000..a7a98f3 --- /dev/null +++ b/test/old_tests/fitting.jl @@ -0,0 +1,12 @@ +#@testset "Continuum" begin +# spec = spectrum([1, 2, 3.], [1, -10, 1.]) +# spec_cont = continuum(spec) +# +# @test spectral_axis(spec_cont) == spectral_axis(spec) +# @test flux_axis(spec_cont) ≈ ones(eltype(flux_axis(spec)), length(flux_axis(spec))) +# @test meta(spec_cont)[:coeffs] == meta(spec_cont)[:coeffs] ≈ [-4.5, 0, 5.5, 0] +# @test meta(spec_cont)[:normalized] +# +# continuum!(spec) +# @test spec == spec_cont +#end diff --git a/test/old_tests/utils.jl b/test/old_tests/utils.jl new file mode 100644 index 0000000..ffebfab --- /dev/null +++ b/test/old_tests/utils.jl @@ -0,0 +1,29 @@ +#using Spectra: blackbody, line_flux, equivalent_width +# +#@testset "Blackbody T=$T" for T in [2000, 4000, 6000] +# wave = range(1e3, 5e4, length = 1000) +# b = 2.897771955185172e7 +# +# bb = @inferred blackbody(wave, T) +# @test typeof(bb) <: Spectra.Spectrum +# @test bb.T == T +# @test spectral_axis(bb)[argmax(bb)] ≈ b / T rtol = 0.01 +# +# wave *= u"angstrom" +# T *= u"K" +# bb = @inferred blackbody(wave, T) +# @test typeof(bb) <: Spectra.Spectrum +# @test unit(bb)[2] == u"W/m^2/angstrom" +# @test bb.T == T +# @test spectral_axis(bb)[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 diff --git a/test/plotting.jl b/test/plotting.jl index 3f1c12f..7ec70ab 100644 --- a/test/plotting.jl +++ b/test/plotting.jl @@ -14,7 +14,7 @@ using RecipesBase :xlabel => "wave", :ylabel => "flux density", :seriestype => :step) - @test rec[1].args == (ustrip.(spec.wave), Measurements.value.(ustrip.(spec.flux))) + @test rec[1].args == (ustrip.(spectral_axis(spec)), Measurements.value.(ustrip.(flux_axis(spec)))) strip_spec = ustrip(spec) @@ -24,7 +24,7 @@ using RecipesBase :xlabel => "wave", :ylabel => "flux density", :seriestype => :step) - @test rec[1].args == (strip_spec.wave, Measurements.value.(strip_spec.flux)) + @test rec[1].args == (spectral_axis(strip_spec), Measurements.value.(flux_axis(strip_spec))) end @testset "Plotting - Echelle" begin @@ -39,5 +39,5 @@ end :ylabel => "flux density", :seriestype => :step, ) - @test rec[1].args == (spec.wave', spec.flux') + @test rec[1].args == (spectral_axis(spec)', flux_axis(spec)') end diff --git a/test/runtests.jl b/test/runtests.jl index c26ce8b..595b212 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, SingleSpectrum, EchelleSpectrum, IFUSpectrum, spectrum + using Spectra: Spectra, Spectrum, SingleSpectrum, 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 1f135e6..b18af69 100644 --- a/test/spectrum.jl +++ b/test/spectrum.jl @@ -32,13 +32,13 @@ end 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 propertynames(spec) == (:spectral_axis, :flux_axis, :meta, :name) + @test spectral_axis(spec) == spectral_axis(spec) + @test flux_axis(spec) == flux_axis(spec) @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 eltype(spec) == eltype(flux_axis(spec)) + @test spectral_axis(spec) == wave + @test spectral_axis(spec_indexed) == wave @test size(spec) === (1000,) @test length(spec) == 1000 @test maximum(spec) == 1000 ± 1 @@ -47,16 +47,16 @@ end @test argmin(spec) == 134 @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 + @test flux_axis(spec) == flux + @test flux_axis(spec_indexed) == flux + @test Measurements.uncertainty.(flux_axis(spec)) ≈ sigma flux_trimmed = flux[200:800] @test_throws ArgumentError spectrum(wave, flux_trimmed) expected = """ SingleSpectrum(Float64, Measurements.Measurement{Float64}) - wave (1000,): 10000.0 .. 50000.0 - flux (1000,): 100.0 ± -2.8 .. 100.0 ± 0.6 + spectral axis (1000,): 10000.0 .. 50000.0 + flux axis (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" @@ -84,15 +84,15 @@ end spec_I = spec[I] spec_I_expected = spectrum(wave, flux, name=spec.name) - @test (spec_i.name, spec_i.wave, spec_i.flux) == (spec_i_expected.name, spec_i_expected.wave, spec_i_expected.flux) - @test (spec_I.name, spec_I.wave, spec_I.flux) == (spec_I_expected.name, spec_I_expected.wave, spec_I_expected.flux) - @test propertynames(spec) == (:wave, :flux, :meta, :name) - @test propertynames(spec_i) == (:wave, :flux, :meta, :Order, :name) - @test propertynames(spec_I) == (:wave, :flux, :meta, :name, :Orders) - @test Spectra.wave(spec) == spec.wave - @test Spectra.flux(spec) == spec.flux - @test eltype(spec) == eltype(spec.flux) - @test spec.wave == wave + @test (spec_i.name, spectral_axis(spec_i), flux_axis(spec_i)) == (spec_i_expected.name, spectral_axis(spec_i_expected), flux_axis(spec_i_expected)) + @test (spec_I.name, spectral_axis(spec_I), flux_axis(spec_I)) == (spec_I_expected.name, spectral_axis(spec_I_expected), flux_axis(spec_I_expected)) + @test propertynames(spec) == (:spectral_axis, :flux_axis, :meta, :name) + @test propertynames(spec_i) == (:spectral_axis, :flux_axis, :meta, :Order, :name) + @test propertynames(spec_I) == (:spectral_axis, :flux_axis, :meta, :name, :Orders) + @test spectral_axis(spec) == spectral_axis(spec) + @test flux_axis(spec) == flux_axis(spec) + @test eltype(spec) == eltype(flux_axis(spec)) + @test spectral_axis(spec) == wave @test size(spec) == (n_wavs, n_orders) @test length(spec) == n_wavs * n_orders @test maximum(spec) == 1000 ± 1 @@ -101,15 +101,15 @@ end @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) + @test eachcol(Measurements.uncertainty.(flux_axis(spec))) ≈ fill(sigma, n_orders) flux_trimmed = flux[200:800, :] @test_throws ArgumentError spectrum(wave, flux_trimmed) expected = """ EchelleSpectrum(Float64, Measurements.Measurement{Float64}) # orders: 3 - wave (1000, 3): 10000.0 .. 50000.0 - flux (1000, 3): 100.0 ± -2.8 .. 100.0 ± 0.6 + spectral axis (1000, 3): 10000.0 .. 50000.0 + flux axis (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" @@ -122,20 +122,20 @@ end expected = """ IFUSpectrum(Int64, Float64) - wave (5,): 20 .. 200 - flux (5, 10, 6): 0.9210599764489846 .. 0.47778429984485815 + spectral axis (5,): 20 .. 200 + flux axis (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 propertynames(spec) == (:spectral_axis, :flux_axis, :meta, :name) + @test spectral_axis(spec) == spectral_axis(spec) + @test flux_axis(spec) == flux_axis(spec) + @test eltype(spec) == eltype(flux_axis(spec)) + @test spectral_axis(spec) == wave @test size(spec) === (5, 10, 6) @test length(spec) == 300 - @test spec.flux == flux + @test flux_axis(spec) == flux @test spec[:, 1, 1] isa SingleSpectrum @test spec[:, begin:4, begin:3] isa IFUSpectrum end @@ -152,12 +152,12 @@ end funit = u"W/m^2/angstrom" spec = spectrum(wave * u"angstrom", flux * funit, name = "test") - @test spec.wave ≈ wave * u"angstrom" + @test spectral_axis(spec) ≈ wave * u"angstrom" - @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 propertynames(spec) == (:spectral_axis, :flux_axis, :meta, :name) + @test spectral_axis(spec) == spectral_axis(spec) + @test flux_axis(spec) == flux_axis(spec) + @test eltype(spec) == eltype(flux_axis(spec)) @test size(spec) === (1000,) @test length(spec) == 1000 @test maximum(spec) == (1000 ± 1) * funit @@ -174,13 +174,13 @@ end @test f_unit == u"W/m^2/angstrom" strip_spec = ustrip(spec) - @test strip_spec.wave == ustrip.(spec.wave) - @test strip_spec.flux == ustrip.(spec.flux) + @test spectral_axis(strip_spec) == ustrip.(spectral_axis(spec)) + @test flux_axis(strip_spec) == ustrip.(flux_axis(spec)) @test strip_spec.meta == spec.meta expected = """ 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 + spectral axis (1000,): 10000.0 Å .. 50000.0 Å + flux axis (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 @@ -206,12 +206,12 @@ end spec = spectrum(wave, flux, name = "test echelle") - @test spec.wave ≈ wave + @test spectral_axis(spec) ≈ wave - @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 propertynames(spec) == (:spectral_axis, :flux_axis, :meta, :name) + @test spectral_axis(spec) == spectral_axis(spec) + @test flux_axis(spec) == flux_axis(spec) + @test eltype(spec) == eltype(flux_axis(spec)) @test size(spec) === (n_wavs, n_orders) @test length(spec) == n_wavs * n_orders @test maximum(spec) == (1000 ± 1) * funit @@ -228,14 +228,14 @@ end @test f_unit == u"W/m^2/angstrom" strip_spec = ustrip(spec) - @test strip_spec.wave == ustrip.(spec.wave) - @test strip_spec.flux == ustrip.(spec.flux) + @test spectral_axis(strip_spec) == ustrip.(spectral_axis(spec)) + @test flux_axis(strip_spec) == ustrip.(flux_axis(spec)) @test strip_spec.meta == spec.meta 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 - 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 + spectral axis (1000, 3): 10000.0 Å .. 50000.0 Å + flux axis (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 @@ -250,82 +250,82 @@ end spec_2 = spectrum(wave, flux_2, name = "test spectrum") # negation - @test -spec.flux == -flux + @test -flux_axis(spec) == -flux # Scalars / vectors values = [10, randn(size(spec))] for A in values # addition s = spec + A - @test s.wave == spec.wave - @test s.flux ≈ spec.flux .+ A + @test spectral_axis(s) == spectral_axis(spec) + @test flux_axis(s) ≈ flux_axis(spec) .+ A # subtraction s = spec - A - @test s.wave == spec.wave - @test s.flux ≈ spec.flux .- A + @test spectral_axis(s) == spectral_axis(spec) + @test flux_axis(s) ≈ flux_axis(spec) .- A s = A - spec - @test s.wave == spec.wave - @test s.flux ≈ spec.flux .- A + @test spectral_axis(s) == spectral_axis(spec) + @test flux_axis(s) ≈ flux_axis(spec) .- A s = spec_2 - spec - @test s.wave == spec.wave - @test s.flux ≈ spec_2.flux .- spec.flux + @test spectral_axis(s) == spectral_axis(spec) + @test flux_axis(s) ≈ flux_axis(spec_2) .- flux_axis(spec) # multiplication s = spec * A - @test s.wave == spec.wave - @test s.flux ≈ spec.flux .* A + @test spectral_axis(s) == spectral_axis(spec) + @test flux_axis(s) ≈ flux_axis(spec) .* A # division s = spec / A - @test s.wave == spec.wave - @test s.flux ≈ spec.flux ./ A + @test spectral_axis(s) == spectral_axis(spec) + @test flux_axis(s) ≈ flux_axis(spec) ./ A end # Other spectra spec2 = deepcopy(spec) s = spec + spec2 - @test s.wave == spec.wave - @test s.flux ≈ 2 .* spec.flux + @test spectral_axis(s) == spectral_axis(spec) + @test flux_axis(s) ≈ 2 .* flux_axis(spec) s = spec - spec2 - @test s.wave == spec.wave - @test s.flux ≈ zeros(size(spec)) + @test spectral_axis(s) == spectral_axis(spec) + @test flux_axis(s) ≈ zeros(size(spec)) s = spec * spec2 - @test s.wave == spec.wave - @test s.flux ≈ spec.flux.^2 + @test spectral_axis(s) == spectral_axis(spec) + @test flux_axis(s) ≈ flux_axis(spec).^2 s = spec / spec2 - @test s.wave == spec.wave - @test s.flux ≈ ones(size(spec)) + @test spectral_axis(s) == spectral_axis(spec) + @test flux_axis(s) ≈ ones(size(spec)) - spec = spectrum(spec.wave * u"cm", spec.flux * u"W/m^2/cm", name = "test unitfulspectrum") + spec = spectrum(spectral_axis(spec) * u"cm", flux_axis(spec) * u"W/m^2/cm", name = "test unitfulspectrum") # Scalars / vectors for A in [10u"W/m^2/cm", randn(size(spec))u"W/m^2/cm"] # addition s = spec + A - @test s.wave == spec.wave - @test s.flux ≈ spec.flux .+ A + @test spectral_axis(s) == spectral_axis(spec) + @test flux_axis(s) ≈ flux_axis(spec) .+ A # subtraction s = spec - A - @test s.wave == spec.wave - @test s.flux ≈ spec.flux .- A + @test spectral_axis(s) == spectral_axis(spec) + @test flux_axis(s) ≈ flux_axis(spec) .- A end for A in [10, randn(size(spec))] # multiplication s = spec * A - @test s.wave == spec.wave - @test s.flux ≈ spec.flux .* A + @test spectral_axis(s) == spectral_axis(spec) + @test flux_axis(s) ≈ flux_axis(spec) .* A # division s = spec / 10 - @test s.wave == spec.wave - @test s.flux ≈ spec.flux ./ 10 + @test spectral_axis(s) == spectral_axis(spec) + @test flux_axis(s) ≈ flux_axis(spec) ./ 10 end end diff --git a/test/transforms/redden.jl b/test/transforms/redden.jl index 8377c93..3f17642 100644 --- a/test/transforms/redden.jl +++ b/test/transforms/redden.jl @@ -31,15 +31,15 @@ end spec2 = deepcopy(spec) reddened = @inferred redden(spec, Av) @inferred redden!(spec2, Av) - @test reddened.flux ≈ spec2.flux + @test flux_axis(reddened) ≈ flux_axis(spec2) @inferred deredden!(spec2, Av) - @test spec.flux ≈ spec2.flux + @test flux_axis(spec) ≈ flux_axis(spec2) dereddened = @inferred deredden(reddened, Av) - @test dereddened.flux ≈ spec.flux + @test flux_axis(dereddened) ≈ flux_axis(spec) # Custom law - expected = @. spec.flux * 10^(-0.4 * Av * CustomLaw(π)(spec.wave)) - @test expected ≈ redden(spec, Av; law=CustomLaw, Rv=π).flux + expected = @. spec.flux_axis * 10^(-0.4 * Av * CustomLaw(π)(spec.spectral_axis)) + #@test expected ≈ redden(spec, Av; law=CustomLaw, Rv=π) |> flux_axis # Bad law @test_throws MethodError redden(spec, Av, law = sin) @@ -49,11 +49,11 @@ end spec2 = deepcopy(spec) reddened = @inferred redden(spec, Av) @inferred redden!(spec2, Av) - @test reddened.flux ≈ spec2.flux + @test flux_axis(reddened) ≈ flux_axis(spec2) @inferred deredden!(spec2, Av) - @test spec.flux ≈ spec2.flux + @test flux_axis(spec) ≈ flux_axis(spec2) dereddened = @inferred deredden(reddened, Av) - @test dereddened.flux ≈ spec.flux + @test flux_axis(dereddened) ≈ flux_axis(spec) end @testset "Reddening Av=0" begin @@ -61,8 +61,8 @@ end Av = 0.0 spec = mock_spectrum() reddened = redden(spec, Av) - @test reddened.flux ≈ spec.flux + @test flux_axis(reddened) ≈ flux_axis(spec) dereddened = deredden(reddened, Av) - @test dereddened.flux ≈ spec.flux + @test flux_axis(dereddened) ≈ flux_axis(spec) end end diff --git a/test/utils.jl b/test/utils.jl deleted file mode 100644 index 27ddc62..0000000 --- a/test/utils.jl +++ /dev/null @@ -1,29 +0,0 @@ -using Spectra: spectrum, blackbody, line_flux, equivalent_width - -@testset "Blackbody T=$T" for T in [2000, 4000, 6000] - wave = range(1e3, 5e4, length = 1000) - b = 2.897771955185172e7 - - bb = @inferred blackbody(wave, T) - @test typeof(bb) <: Spectra.Spectrum - @test bb.T == T - @test bb.wave[argmax(bb)] ≈ b / T rtol = 0.01 - - wave *= u"angstrom" - T *= u"K" - bb = @inferred blackbody(wave, T) - @test typeof(bb) <: Spectra.Spectrum - @test unit(bb)[2] == u"W/m^2/angstrom" - @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 From 5daab57e47d5e514cabbe98ee78edfa4362dd6d9 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Mon, 15 Dec 2025 17:51:50 -0800 Subject: [PATCH 27/31] finish shelving old functionality --- src/Spectra.jl | 2 +- src/old_functionality/utils.jl | 81 ++++++++++++++++++++++++++++++++++ src/utils.jl | 81 ---------------------------------- 3 files changed, 82 insertions(+), 82 deletions(-) create mode 100644 src/old_functionality/utils.jl delete mode 100644 src/utils.jl diff --git a/src/Spectra.jl b/src/Spectra.jl index bf3a53a..4401bfc 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -295,7 +295,7 @@ function spectrum(spectral_axis::AbstractMatrix{<:Quantity}, flux_axis::Abstract end # tools -include("utils.jl") +#include("utils.jl") include("transforms/transforms.jl") include("plotting.jl") #include("fitting/fitting.jl") diff --git a/src/old_functionality/utils.jl b/src/old_functionality/utils.jl new file mode 100644 index 0000000..2225e75 --- /dev/null +++ b/src/old_functionality/utils.jl @@ -0,0 +1,81 @@ +#""" +# blackbody(wave::Vector{<:Quantity}, T::Quantity) +# blackbody(wave::Vector{<:Real}, T::Real) +# +#Create a blackbody spectrum using Planck's law. The curve follows the mathematical form +# +#``B_λ(T) = \\frac{2hc^2}{λ^5} \\frac{1}{e^{hc/λ k_B T} - 1}`` +# +#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. +# +## References +#[Planck's Law](https://en.wikipedia.org/wiki/Planck%27s_law) +# +## Examples +#```jldoctest +#julia> using Spectra, Unitful, UnitfulAstro +# +#julia> wave = range(1, 3, length=100)u"μm" +#(1.0:0.020202020202020204:3.0) μm +# +#julia> bb = blackbody(wave, 2000u"K") +#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) +#SingleSpectrum(Float64, Float64) +# wave (100,): 10000.0 .. 30000.0 +# flux (100,): 1190.9562575755397 .. 40.04325690910415 +# meta: Dict{Symbol, Any}(:T => 6000, :name => "Blackbody") +# +#julia> spectral_axis(bb)[argmax(bb)] +#1.4444444444444444 μm +# +#julia> 2898u"μm*K" / bb.T # See if it matches up with Wien's law +#1.449 μm +#``` +#""" +#function 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 +# +#function blackbody(wave::AbstractVector{<:Real}, T::Real) +# flux = ustrip.(u"W/m^2/angstrom", _blackbody(wave * u"angstrom", T * u"K")) +# return spectrum(wave, flux, name = "Blackbody", T = T) +#end +# +#_blackbody(wave::AbstractVector{<:Quantity}, T::Quantity) = blackbody(T).(wave) +# +#""" +# blackbody(T::Quantity) +# +#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 = spectral_axis(spec)[end] - spectral_axis(spec)[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(spectral_axis(spec)) +# return sum(flux_axis(spec)[2:end] .* avg_dx) +#end diff --git a/src/utils.jl b/src/utils.jl deleted file mode 100644 index f9df7d4..0000000 --- a/src/utils.jl +++ /dev/null @@ -1,81 +0,0 @@ -""" - blackbody(wave::Vector{<:Quantity}, T::Quantity) - blackbody(wave::Vector{<:Real}, T::Real) - -Create a blackbody spectrum using Planck's law. The curve follows the mathematical form - -``B_λ(T) = \\frac{2hc^2}{λ^5} \\frac{1}{e^{hc/λ k_B T} - 1}`` - -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. - -# References -[Planck's Law](https://en.wikipedia.org/wiki/Planck%27s_law) - -# Examples -```jldoctest -julia> using Spectra, Unitful, UnitfulAstro - -julia> wave = range(1, 3, length=100)u"μm" -(1.0:0.020202020202020204:3.0) μm - -julia> bb = blackbody(wave, 2000u"K") -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) -SingleSpectrum(Float64, Float64) - wave (100,): 10000.0 .. 30000.0 - flux (100,): 1190.9562575755397 .. 40.04325690910415 - meta: Dict{Symbol, Any}(:T => 6000, :name => "Blackbody") - -julia> spectral_axis(bb)[argmax(bb)] -1.4444444444444444 μm - -julia> 2898u"μm*K" / bb.T # See if it matches up with Wien's law -1.449 μm -``` -""" -function 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 - -function blackbody(wave::AbstractVector{<:Real}, T::Real) - flux = ustrip.(u"W/m^2/angstrom", _blackbody(wave * u"angstrom", T * u"K")) - return spectrum(wave, flux, name = "Blackbody", T = T) -end - -_blackbody(wave::AbstractVector{<:Quantity}, T::Quantity) = blackbody(T).(wave) - -""" - blackbody(T::Quantity) - -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 = spectral_axis(spec)[end] - spectral_axis(spec)[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(spectral_axis(spec)) - return sum(flux_axis(spec)[2:end] .* avg_dx) -end From 0e4d3ce80caa46148a23def422675171d1e9c9b9 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Mon, 15 Dec 2025 17:55:50 -0800 Subject: [PATCH 28/31] finish spectral axis rename --- src/spectrum_echelle.jl | 4 ++-- src/spectrum_ifu.jl | 2 +- src/spectrum_single.jl | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/spectrum_echelle.jl b/src/spectrum_echelle.jl index e3c4873..565f7ab 100644 --- a/src/spectrum_echelle.jl +++ b/src/spectrum_echelle.jl @@ -3,7 +3,7 @@ An instance of [`Spectrum`](@ref) where the spectral and flux axes are both 2D arrays, i.e., ``M = N = 2``. -The spectral and flux matrices are both ``m`` rows in wavelength by ``n`` columns in [echelle order](https://en.wikipedia.org/wiki/Echelle_grating). +The spectral and flux matrices are both ``m`` rows in spectral axis by ``n`` columns in [echelle order](https://en.wikipedia.org/wiki/Echelle_grating). # Examples @@ -50,7 +50,7 @@ SingleSpectrum(Int64, Int64) See [`SingleSpectrum`](@ref) for a 1D variant, and [`IFUSpectrum`](@ref) for a 3D variant. """ -const EchelleSpectrum = Spectrum{W, F, 2, 2} where {W, F} +const EchelleSpectrum = Spectrum{S, F, 2, 2} where {S, F} function Base.getindex(spec::EchelleSpectrum, i::Int) w = spectral_axis(spec)[:, i] diff --git a/src/spectrum_ifu.jl b/src/spectrum_ifu.jl index 4687eee..91a684b 100644 --- a/src/spectrum_ifu.jl +++ b/src/spectrum_ifu.jl @@ -49,7 +49,7 @@ SingleSpectrum(Int64, Float64) See [`SingleSpectrum`](@ref) for a 1D variant, and [`EchelleSpectrum`](@ref) for a 2D variant. """ -const IFUSpectrum = Spectrum{W, F, 1, 3} where {W, F} +const IFUSpectrum = Spectrum{S, F, 1, 3} where {S, F} Base.getindex(spec::IFUSpectrum, i) = flux_axis(spec)[i, :, :] diff --git a/src/spectrum_single.jl b/src/spectrum_single.jl index 2615fd0..aed1dba 100644 --- a/src/spectrum_single.jl +++ b/src/spectrum_single.jl @@ -17,7 +17,7 @@ SingleSpectrum(Int64, Int64) See [`EchelleSpectrum`](@ref) and [`IFUSpectrum`](@ref) for working with instances of higher dimensional spectra. """ -const SingleSpectrum = Spectrum{W, F, 1, 1} where {W, F} +const SingleSpectrum = Spectrum{S, F, 1, 1} where {S, F} #Base.size(spec::SingleSpectrum) = (length(spectral_axis(spec)), ) #Base.IndexStyle(::Type{<:SingleSpectrum}) = IndexLinear() From 6ae3a804f5cd7b796887e1e8ceed68ad4640cd14 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Mon, 15 Dec 2025 18:38:29 -0800 Subject: [PATCH 29/31] update docs --- docs/make.jl | 4 +- docs/{src => old_docs}/fitting.md | 0 docs/{src => old_docs}/utils.md | 0 docs/src/index.md | 18 +----- docs/src/spectrum.md | 4 +- docs/src/transforms.md | 16 ++--- src/Spectra.jl | 34 +++++----- src/fitting/fitting.jl | 46 -------------- src/old_functionality/fitting/fitting.jl | 46 ++++++++++++++ src/old_functionality/utils.jl | 81 ------------------------ src/utils.jl | 81 ++++++++++++++++++++++++ test/old_tests/utils.jl | 29 --------- test/utils.jl | 29 +++++++++ 13 files changed, 187 insertions(+), 201 deletions(-) rename docs/{src => old_docs}/fitting.md (100%) rename docs/{src => old_docs}/utils.md (100%) delete mode 100644 src/fitting/fitting.jl create mode 100644 src/old_functionality/fitting/fitting.jl delete mode 100644 src/old_functionality/utils.jl create mode 100644 src/utils.jl delete mode 100644 test/old_tests/utils.jl create mode 100644 test/utils.jl diff --git a/docs/make.jl b/docs/make.jl index e0cf339..459e71e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -20,8 +20,8 @@ makedocs(sitename = "Spectra.jl", "Home" => "index.md", "spectrum.md", "transforms.md", - "fitting.md", - "utils.md", + #"fitting.md", + #"utils.md", "contrib.md", ], warnonly = [:missing_docs], diff --git a/docs/src/fitting.md b/docs/old_docs/fitting.md similarity index 100% rename from docs/src/fitting.md rename to docs/old_docs/fitting.md diff --git a/docs/src/utils.md b/docs/old_docs/utils.md similarity index 100% rename from docs/src/utils.md rename to docs/old_docs/utils.md diff --git a/docs/src/index.md b/docs/src/index.md index cae3e18..176afb9 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -44,8 +44,8 @@ julia> flux = (read(f[2], "flux") .* 1e-17)u"erg/s/cm^2/angstrom"; julia> spec = spectrum(wave, flux) 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 + spectral axis (3827,): 3815.0483f0 Å .. 9206.613f0 Å + flux axis (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); @@ -55,20 +55,6 @@ julia> plot(spec); For constructing higher dimensional spectra, e.g., for echelle or IFU spectra, see the docstrings for [EchelleSpectrum](@ref) and [IFUSpectrum](@ref), respectively. -### Continuum fitting - -```jldoctest guide -julia> cont_fit = continuum(spec) -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)); -``` - -![](assets/sdss_cont.svg) - ## Citation If you found this software or any derivative work useful in your academic work, I ask that you please cite the code. diff --git a/docs/src/spectrum.md b/docs/src/spectrum.md index 119b217..c680315 100644 --- a/docs/src/spectrum.md +++ b/docs/src/spectrum.md @@ -27,8 +27,8 @@ For more advanced transformations, see [Transformations](@ref) ### Getters ```@docs -Spectra.wave(::AbstractSpectrum) -Spectra.flux(::AbstractSpectrum) +Spectra.spectral_axis(::AbstractSpectrum) +Spectra.flux_axis(::AbstractSpectrum) Spectra.meta(::AbstractSpectrum) ``` diff --git a/docs/src/transforms.md b/docs/src/transforms.md index 0f25a5a..57c5c31 100644 --- a/docs/src/transforms.md +++ b/docs/src/transforms.md @@ -23,17 +23,17 @@ julia> flux = (100 .± sigma)u"W/m^2/μm" julia> spec = spectrum(wave, flux) 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 + spectral axis (5,): 1.0 μm .. 3.0 μm + flux axis (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) 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 + spectral axis (5,): 1.0 μm .. 3.0 μm + flux axis (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 +julia> flux_axis(red) 5-element Vector{Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(μm^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}}: 89.44 ± 0.84 W μm^-1 m^-2 94.35 ± 0.13 W μm^-1 m^-2 @@ -43,11 +43,11 @@ julia> red.flux julia> deredden!(red, 0.3) 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 + spectral axis (5,): 1.0 μm .. 3.0 μm + flux axis (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 +julia> flux_axis(red) ≈ flux_axis(spec) true ``` diff --git a/src/Spectra.jl b/src/Spectra.jl index 4401bfc..e1ad5f9 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -5,7 +5,7 @@ export AbstractSpectrum, Spectrum, spectrum, spectral_axis, flux_axis # spectra_single.jl, spectra_ifu.jl, spectra_echelle.jl, spectra_binned.jl export SingleSpectrum, IFUSpectrum, EchelleSpectrum, spectral_axis, flux_axis # utils.jl -#export blackbody, line_flux, equivalent_width +export blackbody #, line_flux, equivalent_width # fitting/fitting.jl #export continuum, continuum! # transforms/redden.jl @@ -46,9 +46,9 @@ mutable struct Spectrum{S<:Number, F<:Number, M, N} <: AbstractSpectrum{S, F} """ Spectral axis and flux axis are incompatible sizes. 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) + - 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? @@ -161,14 +161,14 @@ 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 + 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) - wave (1000,): 10000.0 .. 30000.0 - flux (1000,): 99999.76809093042 .. 300000.2474309158 + spectral axis (1000,): 10000.0 .. 30000.0 + flux axis (1000,): 99999.76809093042 .. 300000.2474309158 meta: Dict{Symbol, Any}() ``` """ @@ -219,14 +219,14 @@ julia> flux = 100 .* ones(size(wave)); julia> spec = spectrum(wave, flux) SingleSpectrum(Float64, Float64) - wave (1000,): 10000.0 .. 40000.0 - flux (1000,): 100.0 .. 100.0 + spectral axis (1000,): 10000.0 .. 40000.0 + flux axis (1000,): 100.0 .. 100.0 meta: Dict{Symbol, Any}() julia> spec = spectrum(wave, flux, name="Just Noise") SingleSpectrum(Float64, Float64) - wave (1000,): 10000.0 .. 40000.0 - flux (1000,): 100.0 .. 100.0 + spectral axis (1000,): 10000.0 .. 40000.0 + flux axis (1000,): 100.0 .. 100.0 meta: Dict{Symbol, Any}(:name => "Just Noise") julia> spec.name @@ -252,8 +252,8 @@ julia> flux = (100 .± sigma)u"erg/cm^2/s/angstrom"; julia> spec = spectrum(wave, flux) 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 + spectral axis (1000,): 1.0 μm .. 4.0 μm + flux axis (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}() ``` @@ -267,8 +267,8 @@ 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 + spectral axis (100, 10): 100.0 .. 10000.0 + flux axis (100, 10): 1.0 .. 10.0 meta: Dict{Symbol, Any}() ``` """ @@ -295,7 +295,7 @@ function spectrum(spectral_axis::AbstractMatrix{<:Quantity}, flux_axis::Abstract end # tools -#include("utils.jl") +include("utils.jl") include("transforms/transforms.jl") include("plotting.jl") #include("fitting/fitting.jl") diff --git a/src/fitting/fitting.jl b/src/fitting/fitting.jl deleted file mode 100644 index 37f0154..0000000 --- a/src/fitting/fitting.jl +++ /dev/null @@ -1,46 +0,0 @@ -using LinearAlgebra: /, \, diagm, pinv - - -function chebvander(x::AbstractVector{T}, deg::Int) where {T <: Number} - v = Matrix{T}(undef, length(x), deg + 1) - v[:, 1] .= one(T) - x2 = 2 .* x - v[:, 2] .= x - @inbounds for i in 3:deg + 1 - v[:, i] .= v[:, i - 1] .* x2 .- v[:, i - 2] - end - return v -end - -function chebfit(wave, flux::AbstractVector{T}, deg) where {T} - # gotta be a simpler way to convert to -1 to 1 domain - x = (wave .- minimum(wave)) .* 2 ./ (maximum(wave) - minimum(wave)) .- 1 - vand = chebvander(x, deg) - if T <: Measurement - W = diagm(0 => Measurements.uncertainty.(flux)) - coeffs = (vand' * W * vand) \ (vand' * W * flux) - else - coeffs = pinv(vand) * flux - end - fit = vand * coeffs - return coeffs, fit -end - -""" - continuum!(::AbstractSpectrum, deg::Int=3) - -In-place version of [`continuum`](@ref) -""" -function continuum!(spec::AbstractSpectrum, deg::Int = 3) - coeffs, fit = chebfit(spec.wave, spec.flux, deg) - spec.flux .= spec.flux ./ ustrip.(fit) - merge!(spec.meta, Dict(:normalized => true, :coeffs => coeffs)) - return spec -end - -""" - continuum(::AbstractSpectrum, deg::Int=3) - -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) diff --git a/src/old_functionality/fitting/fitting.jl b/src/old_functionality/fitting/fitting.jl new file mode 100644 index 0000000..fc6106e --- /dev/null +++ b/src/old_functionality/fitting/fitting.jl @@ -0,0 +1,46 @@ +#using LinearAlgebra: /, \, diagm, pinv +# +# +#function chebvander(x::AbstractVector{T}, deg::Int) where {T <: Number} +# v = Matrix{T}(undef, length(x), deg + 1) +# v[:, 1] .= one(T) +# x2 = 2 .* x +# v[:, 2] .= x +# @inbounds for i in 3:deg + 1 +# v[:, i] .= v[:, i - 1] .* x2 .- v[:, i - 2] +# end +# return v +#end +# +#function chebfit(wave, flux::AbstractVector{T}, deg) where {T} +# # gotta be a simpler way to convert to -1 to 1 domain +# x = (wave .- minimum(wave)) .* 2 ./ (maximum(wave) - minimum(wave)) .- 1 +# vand = chebvander(x, deg) +# if T <: Measurement +# W = diagm(0 => Measurements.uncertainty.(flux)) +# coeffs = (vand' * W * vand) \ (vand' * W * flux) +# else +# coeffs = pinv(vand) * flux +# end +# fit = vand * coeffs +# return coeffs, fit +#end +# +#""" +# continuum!(::AbstractSpectrum, deg::Int=3) +# +#In-place version of [`continuum`](@ref) +#""" +#function continuum!(spec::AbstractSpectrum, deg::Int = 3) +# coeffs, fit = chebfit(spec.wave, spec.flux, deg) +# spec.flux .= spec.flux ./ ustrip.(fit) +# merge!(spec.meta, Dict(:normalized => true, :coeffs => coeffs)) +# return spec +#end +# +#""" +# continuum(::AbstractSpectrum, deg::Int=3) +# +#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) diff --git a/src/old_functionality/utils.jl b/src/old_functionality/utils.jl deleted file mode 100644 index 2225e75..0000000 --- a/src/old_functionality/utils.jl +++ /dev/null @@ -1,81 +0,0 @@ -#""" -# blackbody(wave::Vector{<:Quantity}, T::Quantity) -# blackbody(wave::Vector{<:Real}, T::Real) -# -#Create a blackbody spectrum using Planck's law. The curve follows the mathematical form -# -#``B_λ(T) = \\frac{2hc^2}{λ^5} \\frac{1}{e^{hc/λ k_B T} - 1}`` -# -#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. -# -## References -#[Planck's Law](https://en.wikipedia.org/wiki/Planck%27s_law) -# -## Examples -#```jldoctest -#julia> using Spectra, Unitful, UnitfulAstro -# -#julia> wave = range(1, 3, length=100)u"μm" -#(1.0:0.020202020202020204:3.0) μm -# -#julia> bb = blackbody(wave, 2000u"K") -#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) -#SingleSpectrum(Float64, Float64) -# wave (100,): 10000.0 .. 30000.0 -# flux (100,): 1190.9562575755397 .. 40.04325690910415 -# meta: Dict{Symbol, Any}(:T => 6000, :name => "Blackbody") -# -#julia> spectral_axis(bb)[argmax(bb)] -#1.4444444444444444 μm -# -#julia> 2898u"μm*K" / bb.T # See if it matches up with Wien's law -#1.449 μm -#``` -#""" -#function 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 -# -#function blackbody(wave::AbstractVector{<:Real}, T::Real) -# flux = ustrip.(u"W/m^2/angstrom", _blackbody(wave * u"angstrom", T * u"K")) -# return spectrum(wave, flux, name = "Blackbody", T = T) -#end -# -#_blackbody(wave::AbstractVector{<:Quantity}, T::Quantity) = blackbody(T).(wave) -# -#""" -# blackbody(T::Quantity) -# -#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 = spectral_axis(spec)[end] - spectral_axis(spec)[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(spectral_axis(spec)) -# return sum(flux_axis(spec)[2:end] .* avg_dx) -#end diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 0000000..919f921 --- /dev/null +++ b/src/utils.jl @@ -0,0 +1,81 @@ +""" + blackbody(wave::Vector{<:Quantity}, T::Quantity) + blackbody(wave::Vector{<:Real}, T::Real) + +Create a blackbody spectrum using Planck's law. The curve follows the mathematical form + +``B_λ(T) = \\frac{2hc^2}{λ^5} \\frac{1}{e^{hc/λ k_B T} - 1}`` + +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. + +# References +[Planck's Law](https://en.wikipedia.org/wiki/Planck%27s_law) + +# Examples +```jldoctest +julia> using Spectra, Unitful, UnitfulAstro + +julia> wave = range(1, 3, length=100)u"μm" +(1.0:0.020202020202020204:3.0) μm + +julia> bb = blackbody(wave, 2000u"K") +SingleSpectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(μm,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(μm^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}) + spectral axis (100,): 1.0 μm .. 3.0 μm + flux axis (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) +SingleSpectrum(Float64, Float64) + spectral axis (100,): 10000.0 .. 30000.0 + flux axis (100,): 1190.9562575755397 .. 40.04325690910415 + meta: Dict{Symbol, Any}(:T => 6000, :name => "Blackbody") + +julia> spectral_axis(bb)[argmax(bb)] +1.4444444444444444 μm + +julia> 2898u"μm*K" / bb.T # See if it matches up with Wien's law +1.449 μm +``` +""" +function 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 + +function blackbody(wave::AbstractVector{<:Real}, T::Real) + flux = ustrip.(u"W/m^2/angstrom", _blackbody(wave * u"angstrom", T * u"K")) + return spectrum(wave, flux, name = "Blackbody", T = T) +end + +_blackbody(wave::AbstractVector{<:Quantity}, T::Quantity) = blackbody(T).(wave) + +""" + blackbody(T::Quantity) + +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 = spectral_axis(spec)[end] - spectral_axis(spec)[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(spectral_axis(spec)) +# return sum(flux_axis(spec)[2:end] .* avg_dx) +#end diff --git a/test/old_tests/utils.jl b/test/old_tests/utils.jl deleted file mode 100644 index ffebfab..0000000 --- a/test/old_tests/utils.jl +++ /dev/null @@ -1,29 +0,0 @@ -#using Spectra: blackbody, line_flux, equivalent_width -# -#@testset "Blackbody T=$T" for T in [2000, 4000, 6000] -# wave = range(1e3, 5e4, length = 1000) -# b = 2.897771955185172e7 -# -# bb = @inferred blackbody(wave, T) -# @test typeof(bb) <: Spectra.Spectrum -# @test bb.T == T -# @test spectral_axis(bb)[argmax(bb)] ≈ b / T rtol = 0.01 -# -# wave *= u"angstrom" -# T *= u"K" -# bb = @inferred blackbody(wave, T) -# @test typeof(bb) <: Spectra.Spectrum -# @test unit(bb)[2] == u"W/m^2/angstrom" -# @test bb.T == T -# @test spectral_axis(bb)[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 diff --git a/test/utils.jl b/test/utils.jl new file mode 100644 index 0000000..fa2f7e3 --- /dev/null +++ b/test/utils.jl @@ -0,0 +1,29 @@ +using Spectra: blackbody#, line_flux, equivalent_width + +@testset "Blackbody T=$T" for T in [2000, 4000, 6000] + wave = range(1e3, 5e4, length = 1000) + b = 2.897771955185172e7 + + bb = @inferred blackbody(wave, T) + @test typeof(bb) <: Spectra.Spectrum + @test bb.T == T + @test spectral_axis(bb)[argmax(bb)] ≈ b / T rtol = 0.01 + + wave *= u"angstrom" + T *= u"K" + bb = @inferred blackbody(wave, T) + @test typeof(bb) <: Spectra.Spectrum + @test unit(bb)[2] == u"W/m^2/angstrom" + @test bb.T == T + @test spectral_axis(bb)[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 From c9b31f787e41b43c0470e262e7c0cf295908e691 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Tue, 16 Dec 2025 17:04:02 -0800 Subject: [PATCH 30/31] use accessor methods more --- src/transforms/redden.jl | 6 ++++-- test/transforms/redden.jl | 2 +- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/transforms/redden.jl b/src/transforms/redden.jl index f73dbe5..7930d6a 100644 --- a/src/transforms/redden.jl +++ b/src/transforms/redden.jl @@ -6,7 +6,8 @@ import DustExtinction In-place version of [`redden`](@ref) """ function redden!(spec::T, Av; Rv = 3.1, law = DustExtinction.CCM89) where {T <: AbstractSpectrum} - @. spec.flux_axis = DustExtinction.redden(law, spec.spectral_axis, spec.flux_axis; Rv, Av) + s, f = spectral_axis(spec), flux_axis(spec) + @. f = DustExtinction.redden(law, s, f; Rv, Av) return spec end @@ -30,7 +31,8 @@ end In-place version of [`deredden`](@ref) """ function deredden!(spec::AbstractSpectrum, Av; Rv = 3.1, law = DustExtinction.CCM89) - @. spec.flux_axis = DustExtinction.deredden(law, spec.spectral_axis, spec.flux_axis; Rv, Av) + s, f = spectral_axis(spec), flux_axis(spec) + @. f = DustExtinction.deredden(law, s, f; Rv, Av) return spec end diff --git a/test/transforms/redden.jl b/test/transforms/redden.jl index 3f17642..2975360 100644 --- a/test/transforms/redden.jl +++ b/test/transforms/redden.jl @@ -38,7 +38,7 @@ end @test flux_axis(dereddened) ≈ flux_axis(spec) # Custom law - expected = @. spec.flux_axis * 10^(-0.4 * Av * CustomLaw(π)(spec.spectral_axis)) + expected = @. $flux_axis(spec) * 10^(-0.4 * Av * CustomLaw(π)($spectral_axis(spec))) #@test expected ≈ redden(spec, Av; law=CustomLaw, Rv=π) |> flux_axis # Bad law From d147e6acce02e80bbcc70f9dfd78e729def11314 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Tue, 16 Dec 2025 17:13:30 -0800 Subject: [PATCH 31/31] typo --- test/transforms/redden.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/transforms/redden.jl b/test/transforms/redden.jl index 2975360..e899d48 100644 --- a/test/transforms/redden.jl +++ b/test/transforms/redden.jl @@ -38,7 +38,8 @@ end @test flux_axis(dereddened) ≈ flux_axis(spec) # Custom law - expected = @. $flux_axis(spec) * 10^(-0.4 * Av * CustomLaw(π)($spectral_axis(spec))) + s, f = spectral_axis(spec), flux_axis(spec) + expected = @. f * 10^(-0.4 * Av * CustomLaw(π)(s)) #@test expected ≈ redden(spec, Av; law=CustomLaw, Rv=π) |> flux_axis # Bad law