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 950496e..459e71e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -2,12 +2,15 @@ using Documenter using Spectra using Unitful using Measurements +using Revise + +Revise.revise() DocMeta.setdocmeta!(Spectra, :DocTestSetup, :(using Spectra); recursive = true) makedocs(sitename = "Spectra.jl", format = Documenter.HTML(; - prettyurls = get(ENV, "CI", nothing) == "true", + prettyurls = true, canonical = "https://juliaastro.org/Spectra/stable/", ), authors = "Miles Lucas and contributors.", @@ -17,8 +20,8 @@ makedocs(sitename = "Spectra.jl", "Home" => "index.md", "spectrum.md", "transforms.md", - "fitting.md", - "analysis.md", + #"fitting.md", + #"utils.md", "contrib.md", ], warnonly = [:missing_docs], @@ -27,5 +30,7 @@ makedocs(sitename = "Spectra.jl", deploydocs(; repo = "github.com/JuliaAstro/Spectra.jl.git", - versions = ["stable" => "v^", "v#.#"] # Restrict to minor releases + devbranch = "main", + push_preview = true, + versions = ["stable" => "v^", "v#.#"], # Restrict to minor releases ) diff --git a/docs/src/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/analysis.md b/docs/old_docs/utils.md similarity index 60% rename from docs/src/analysis.md rename to docs/old_docs/utils.md index 6592083..8fb1048 100644 --- a/docs/src/analysis.md +++ b/docs/old_docs/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/docs/src/index.md b/docs/src/index.md index 92b6a7f..176afb9 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 @@ -41,23 +43,17 @@ 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}}) + 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); ``` ![](assets/sdss.svg) -```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 - -julia> plot(cont_fit, xlims=(6545, 6600)); -``` - -![](assets/sdss_cont.svg) +For constructing higher dimensional spectra, e.g., for echelle or IFU spectra, see the docstrings for [EchelleSpectrum](@ref) and [IFUSpectrum](@ref), respectively. ## Citation diff --git a/docs/src/spectrum.md b/docs/src/spectrum.md index 39995dd..c680315 100644 --- a/docs/src/spectrum.md +++ b/docs/src/spectrum.md @@ -1,10 +1,3 @@ -```@meta -DocTestSetup = quote - using Spectra, Random - Random.seed!(11894) -end -``` - # Spectrum Here we will go over the different spectral types and how we use them. @@ -14,7 +7,11 @@ Here we will go over the different spectral types and how we use them. Spectra are defined as possible subtypes of `AbstractSpectrum`. You can use these directly for construction, or use the catch-all [`spectrum`](@ref) function, which is preferred. ```@docs +Spectra.AbstractSpectrum Spectra.Spectrum +Spectra.SingleSpectrum +Spectra.EchelleSpectrum +Spectra.IFUSpectrum ``` ## Constructors @@ -28,25 +25,37 @@ Spectra.spectrum For more advanced transformations, see [Transformations](@ref) +### Getters +```@docs +Spectra.spectral_axis(::AbstractSpectrum) +Spectra.flux_axis(::AbstractSpectrum) +Spectra.meta(::AbstractSpectrum) +``` + +### Array interface + | Function | |:-----------------------------------| +| `Base.argmax(::AbstractSpectrum)` | +| `Base.argmin(::AbstractSpectrum)` | +| `Base.eltype(::AbstractSpectrum)` | +| `Base.findmax(::AbstractSpectrum)` | +| `Base.findmin(::AbstractSpectrum)` | +| `Base.iterate(::AbstractSpectrum)` | | `Base.length(::AbstractSpectrum)` | -| `Base.size(::AbstractSpectrum)` | | `Base.maximum(::AbstractSpectrum)` | | `Base.minimum(::AbstractSpectrum)` | -| `Base.argmax(::AbstractSpectrum)` | -| `Base.argmin(::AbstractSpectrum)` | -| `Base.findmax(::AbstractSpectrum)` | -| `Base.findmin(::AbstractSpectrum)` | +| `Base.size(::AbstractSpectrum)` | ### Arithmetic -| Function | -|:-----------------------------------| -| `+(::AbstractSpectrum, A)` | -| `-(::AbstractSpectrum, A)` | -| `*(::AbstractSpectrum, A)` | -| `/(::AbstractSpectrum, A)` | +| Function | +|:----------------------------------------------------| +| `+(::AbstractSpectrum, A)` | +| `-(::AbstractSpectrum, A)` | +| `*(::AbstractSpectrum, A)` | +| `/(::AbstractSpectrum, A)` | +| `Base.(==)(::AbstractSpectrum, ::AbstractSpectrum)` | ## Unitful helpers @@ -76,7 +85,3 @@ savefig("spec-plot.svg"); nothing # hide ```@index Pages = ["spectrum.md"] ``` - -```@meta -DocTestSetup = nothing -``` diff --git a/docs/src/transforms.md b/docs/src/transforms.md index a5a05c5..57c5c31 100644 --- a/docs/src/transforms.md +++ b/docs/src/transforms.md @@ -1,10 +1,3 @@ -```@meta -DocTestSetup = quote - using Spectra, Random - Random.seed!(11894) -end -``` - # Transformations ## Extinction @@ -12,20 +5,13 @@ end By levaraging [DustExtinction.jl](https://github.com/juliaastro/dustextinction.jl) we can apply common reddening laws to our spectra. ```jldoctest -julia> using Unitful, Measurements, Random +julia> using Spectra, Unitful, Measurements, Random julia> rng = Random.seed!(0); -julia> wave = (1:0.5:3)u"Ξm" -(1.0:0.5:3.0) Ξm +julia> wave = (1:0.5:3)u"Ξm"; -julia> sigma = randn(rng, size(wave)) -5-element Vector{Float64}: - 0.942970533446119 - 0.13392275765318448 - 1.5250689085124804 - 0.12390123120559722 - -1.205772284259936 +julia> sigma = randn(rng, size(wave)); julia> flux = (100 .Âą sigma)u"W/m^2/Ξm" 5-element Vector{Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Ξm^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}}: @@ -36,12 +22,18 @@ 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}}) + 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) -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}}) + 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 @@ -50,9 +42,12 @@ 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}}) + 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 ``` @@ -64,7 +59,3 @@ redden! deredden deredden! ``` - -```@meta -DocTestSetup = nothing -``` 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..e1ad5f9 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -1,9 +1,13 @@ module Spectra -# common.jl -export AbstractSpectrum, spectrum +# Uniform API +export AbstractSpectrum, Spectrum, spectrum, 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 +export blackbody #, line_flux, equivalent_width +# fitting/fitting.jl +#export continuum, continuum! # transforms/redden.jl export redden, redden!, deredden, deredden! @@ -12,17 +16,200 @@ using Measurements: Measurements, Measurement using Unitful: Unitful, Quantity, @u_str, ustrip, unit, dimension using PhysicalConstants.CODATA2018: h, c_0, k_B -# AbstractSpectrum and common functionality -include("common.jl") +""" + AbstractSpectrum{S <: Number, F <: Number} + +An abstract holder for astronomical spectra. All types inheriting from this must have the following fields: + +* `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{S, F} end + +""" + Spectrum <: AbstractSpectrum + +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{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 + if size(s, 1) != size(f, 1) + throw(ArgumentError( + """ + Spectral axis and flux axis are incompatible sizes. Currently supported sizes are: + + - 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? + + 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(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) + +""" + spectral_axis(spec::AbstractSpectrum) + +Return the spectral axis of `spec`. +""" +spectral_axis(spec::AbstractSpectrum) = spec.spectral_axis + +""" + flux(spec::AbstractSpectrum) + +Return the flux axis of `spec`. +""" +flux_axis(spec::AbstractSpectrum) = spec.flux_axis + +""" + meta(spec::AbstractSpectrum) + +Return the meta data associated with `spec`. +""" +meta(spec::AbstractSpectrum) = spec.meta + +function Base.getproperty(spec::AbstractSpectrum, nm::Symbol) + if nm in keys(getfield(spec, :meta)) + return getfield(spec, :meta)[nm] + else + return getfield(spec, nm) + end +end + +function Base.propertynames(spec::AbstractSpectrum) + natural = (:spectral_axis, :flux_axis, :meta) + m = keys(meta(spec)) + return (natural..., m...) +end + +# Collection +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) = 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(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) + +Remove the units from a spectrum. Useful for processing spectra in tools that don't play nicely with `Unitful.jl` + +# Examples +```jldoctest +julia> using Random + +julia> rng = Random.seed!(0) +TaskLocalRNG() + +julia> using Unitful, UnitfulAstro + +julia> wave = range(1e4, 3e4, length=1000); + +julia> flux = wave .* 10 .+ randn(rng, 1000); + +julia> spec = spectrum(wave*u"angstrom", flux*u"W/m^2/angstrom") +SingleSpectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}) + spectral axis (1000,): 10000.0 Å .. 30000.0 Å + flux axis (1000,): 99999.76809093042 W Å^-1 m^-2 .. 300000.2474309158 W Å^-1 m^-2 + meta: Dict{Symbol, Any}() + +julia> ustrip(spec) +SingleSpectrum(Float64, Float64) + spectral axis (1000,): 10000.0 .. 30000.0 + flux axis (1000,): 99999.76809093042 .. 300000.2474309158 + meta: Dict{Symbol, Any}() +``` +""" +Unitful.ustrip(spec::AbstractSpectrum) = spectrum(ustrip.(spectral_axis(spec)), ustrip.(flux_axis(spec)); meta(spec)...) + +""" + Unitful.unit(::AbstractSpectrum) + +Get the units of a spectrum. Returns a tuple of the spectral axis units and flux/sigma units + +# Examples +```jldoctest +julia> using Random + +julia> rng = Random.seed!(0) +TaskLocalRNG() + +julia> using Unitful, UnitfulAstro + +julia> wave = range(1e4, 3e4, length=1000); + +julia> flux = wave .* 10 .+ randn(rng, 1000); + +julia> spec = spectrum(wave * u"angstrom", flux * u"W/m^2/angstrom"); + +julia> w_unit, f_unit = unit(spec) +(Å, W Å^-1 m^-2) +``` +""" +Unitful.unit(spec::AbstractSpectrum) = unit(eltype(spectral_axis(spec))), unit(eltype(flux_axis(spec))) # Spectrum types and basic arithmetic -include("spectrum.jl") -include("EchelleSpectrum.jl") +include("spectrum_single.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 @@ -31,11 +218,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) + 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") -Spectrum(Float64, Float64) - name: Just Noise +SingleSpectrum(Float64, Float64) + spectral axis (1000,): 10000.0 .. 40000.0 + flux axis (1000,): 100.0 .. 100.0 + meta: Dict{Symbol, Any}(:name => "Just Noise") julia> spec.name "Just Noise" @@ -45,56 +237,67 @@ There is easy integration with [Unitful.jl](https://github.com/JuliaPhysics/Unit and its sub-projects and [Measurements.jl](https://github.com/juliaphysics/measurements.jl) ```jldoctest +julia> using Random + +julia> rng = Random.seed!(0) +TaskLocalRNG() + julia> using Unitful, UnitfulAstro, Measurements julia> wave = range(1, 4, length=1000)u"Ξm"; -julia> sigma = randn(size(wave)); +julia> sigma = randn(rng, size(wave)); julia> flux = (100 .Âą sigma)u"erg/cm^2/s/angstrom"; julia> spec = spectrum(wave, flux) -Spectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(Ξm,), 𝐋, nothing}}, Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, erg, cm^-2, s^-1), 𝐌 𝐋^-1 𝐓^-3, nothing}}) +SingleSpectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(Ξm,), 𝐋, nothing}}, Quantity{Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, erg, cm^-2, s^-1), 𝐌 𝐋^-1 𝐓^-3, nothing}}) + 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}() ``` -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 + spectral axis (100, 10): 100.0 .. 10000.0 + flux axis (100, 10): 1.0 .. 10.0 + meta: Dict{Symbol, Any}() ``` """ -function spectrum(wave::AbstractVector{<:Real}, flux::AbstractVector{<:Real}; kwds...) - @assert size(wave) == size(flux) "wave and flux must have equal size" - Spectrum(wave, flux, Dict{Symbol,Any}(kwds)) +function spectrum(spectral_axis::AbstractVector{<:Real}, flux_axis::AbstractVector{<:Real}; kwds...) + Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds)) +end + +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::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)) +function spectrum(spectral_axis::AbstractMatrix{<:Real}, flux_axis::AbstractMatrix{<:Real}; kwds...) + Spectrum(spectral_axis, flux_axis, 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)) +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 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)) +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/common.jl b/src/common.jl deleted file mode 100644 index 5c0d20c..0000000 --- a/src/common.jl +++ /dev/null @@ -1,108 +0,0 @@ -""" - AbstractSpectrum{W<:Number, F<:Number} - -An abstract holder for astronomical spectra. All types inheriting from this must have the following fields - -- wave::Array{W, N} -- flux::Array{F, N} -- meta::Dict{Symbol, Any} -""" -abstract type AbstractSpectrum{W,F} end - -function Base.getproperty(spec::AbstractSpectrum, nm::Symbol) - if nm in keys(getfield(spec, :meta)) - return getfield(spec, :meta)[nm] - else - return getfield(spec, nm) - end -end - -function Base.propertynames(spec::AbstractSpectrum) - natural = (:wave, :flux, :meta) - meta = keys(spec.meta) - return (natural..., meta...) -end - -""" - wave(::AbstractSpectrum) - -Return the wavelengths of the spectrum. -""" -wave(spec::AbstractSpectrum) = spec.wave - -""" - flux(::AbstractSpectrum) - -Return the flux of the spectrum. -""" -flux(spec::AbstractSpectrum) = spec.flux - -# Collection -Base.eltype(spec::AbstractSpectrum) = eltype(spec.flux) -Base.size(spec::AbstractSpectrum) = size(spec.flux) -Base.size(spec::AbstractSpectrum, i) = size(spec.flux, i) -Base.length(spec::AbstractSpectrum) = length(spec.flux) -Base.maximum(spec::AbstractSpectrum) = maximum(spec.flux) -Base.minimum(spec::AbstractSpectrum) = minimum(spec.flux) -Base.argmax(spec::AbstractSpectrum) = argmax(spec.flux) -Base.argmin(spec::AbstractSpectrum) = argmin(spec.flux) -Base.findmax(spec::AbstractSpectrum) = findmax(spec.flux) -Base.findmin(spec::AbstractSpectrum) = findmin(spec.flux) - -# Arithmetic -Base.:+(s::T, A) where {T <: AbstractSpectrum} = T(s.wave, s.flux .+ A, s.meta) -Base.:*(s::T, A::Union{Real, AbstractVector}) where {T <: AbstractSpectrum} = T(s.wave, s.flux .* A, s.meta) -Base.:/(s::T, A) where {T <: AbstractSpectrum} = T(s.wave, s.flux ./ A, s.meta) -Base.:-(s::T) where {T <: AbstractSpectrum} = T(s.wave, -s.flux, s.meta) -Base.:-(s::AbstractSpectrum, A) = s + -A -Base.:-(A, s::AbstractSpectrum) = s - A -Base.:-(s::AbstractSpectrum, o::AbstractSpectrum) = s - o # Satisfy Aqua - -# Multi-Spectrum -Base.:+(s::T, o::T) where {T <: AbstractSpectrum} = T(s.wave, s.flux .+ o.flux, s.meta) -Base.:*(s::T, o::T) where {T <: AbstractSpectrum} = T(s.wave, s.flux .* o.flux, s.meta) -Base.:/(s::T, o::T) where {T <: AbstractSpectrum} = T(s.wave, s.flux ./ o.flux * unit(s)[2], s.meta) -Base.:-(s::T, o::T) where {T <: AbstractSpectrum} = T(s.wave, s.flux .- o.flux, s.meta) - -""" - Unitful.ustrip(::AbstractSpectrum) - -Remove the units from a spectrum. Useful for processing spectra in tools that don't play nicely with `Unitful.jl` - -# Examples -```jldoctest -julia> using Unitful, UnitfulAstro - -julia> wave = range(1e4, 3e4, length=1000); - -julia> flux = wave .* 10 .+ randn(1000); - -julia> spec = spectrum(wave*u"angstrom", flux*u"W/m^2/angstrom") -Spectrum(Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}) - -julia> ustrip(spec) -Spectrum(Float64, Float64) -``` -""" -Unitful.ustrip(spec::AbstractSpectrum) = spectrum(ustrip.(spec.wave), ustrip.(spec.flux); spec.meta...) - -""" - Unitful.unit(::AbstractSpectrum) - -Get the units of a spectrum. Returns a tuple of the wavelength units and flux/sigma units - -# Examples -```jldoctest -julia> using Unitful, UnitfulAstro - -julia> wave = range(1e4, 3e4, length=1000); - -julia> flux = wave .* 10 .+ randn(1000); - -julia> spec = spectrum(wave * u"angstrom", flux * u"W/m^2/angstrom"); - -julia> w_unit, f_unit = unit(spec) -(Å, W Å^-1 m^-2) -``` -""" -Unitful.unit(spec::AbstractSpectrum) = unit(eltype(spec.wave)), unit(eltype(spec.flux)) diff --git a/src/fitting/fitting.jl b/src/fitting/fitting.jl deleted file mode 100644 index d2e3343..0000000 --- a/src/fitting/fitting.jl +++ /dev/null @@ -1,68 +0,0 @@ -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) - 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) - -""" - 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/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/plotting.jl b/src/plotting.jl index 3515c2d..9fef960 100644 --- a/src/plotting.jl +++ b/src/plotting.jl @@ -1,9 +1,9 @@ -@recipe function f(spec::Spectrum) +@recipe function f(spec::SingleSpectrum) seriestype --> :step 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.jl b/src/spectrum.jl deleted file mode 100644 index c1a20f9..0000000 --- a/src/spectrum.jl +++ /dev/null @@ -1,19 +0,0 @@ -""" - Spectrum <: AbstractSpectrum - -A 1-dimensional spectrum stored as vectors of real numbers. The wavelengths are assumed to be in angstrom. -""" -struct Spectrum{W <: Number,F <: Number} <: AbstractSpectrum{W,F} - wave::Vector{W} - flux::Vector{F} - meta::Dict{Symbol,Any} -end - -Spectrum(wave, flux, meta::Dict{Symbol,Any}) = Spectrum(collect(wave), collect(flux), meta) - -function Base.show(io::IO, spec::Spectrum) - print(io, "Spectrum($(eltype(spec.wave)), $(eltype(spec.flux)))") - for (key, val) in spec.meta - print(io, "\n $key: $val") - end -end diff --git a/src/spectrum_echelle.jl b/src/spectrum_echelle.jl new file mode 100644 index 0000000..565f7ab --- /dev/null +++ b/src/spectrum_echelle.jl @@ -0,0 +1,80 @@ +""" + EchelleSpectrum <: AbstractSpectrum + +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 spectral axis 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 + 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) + spectral axis (10,): 1 .. 10 + flux axis (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{S, F, 2, 2} where {S, F} + +function Base.getindex(spec::EchelleSpectrum, i::Int) + 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 = 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_axis(spec), 1) +Base.lastindex(spec::EchelleSpectrum) = lastindex(flux_axis(spec), 1) + +function Base.show(io::IO, spec::EchelleSpectrum) + 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, " 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 new file mode 100644 index 0000000..91a684b --- /dev/null +++ b/src/spectrum_ifu.jl @@ -0,0 +1,81 @@ +""" + IFUSpectrum <: AbstractSpectrum + +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 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 + +```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) + 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 +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) + 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) + spectral axis (5,): 20 .. 200 + flux axis (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{S, F, 1, 3} where {S, F} + +Base.getindex(spec::IFUSpectrum, i) = flux_axis(spec)[i, :, :] + +function Base.getindex(spec::IFUSpectrum, I::AbstractVector) + 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_axis(spec)[i, j, k] + +function Base.getindex(spec::IFUSpectrum, 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_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 = spectral_axis(spec) + f = flux_axis(spec) + println(io, "IFUSpectrum($(eltype(w)), $(eltype(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 new file mode 100644 index 0000000..aed1dba --- /dev/null +++ b/src/spectrum_single.jl @@ -0,0 +1,43 @@ +""" + SingleSpectrum <: AbstractSpectrum + +An instance of [`Spectrum`](@ref) where the spectral and flux axes are both 1D arrays, i.e., ``M = N = 1``. + +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) + spectral axis (4,): 6 .. 9 + flux axis (4,): 2 .. 5 + meta: Dict{Symbol, Any}() +``` + +See [`EchelleSpectrum`](@ref) and [`IFUSpectrum`](@ref) for working with instances of higher dimensional spectra. +""" +const SingleSpectrum = Spectrum{S, F, 1, 1} where {S, F} + +#Base.size(spec::SingleSpectrum) = (length(spectral_axis(spec)), ) +#Base.IndexStyle(::Type{<:SingleSpectrum}) = IndexLinear() + +function Base.getindex(spec::SingleSpectrum, i::Int) + return Spectrum([spectral_axis(spec)[i]], [flux_axis(spec)[i]], meta(spec)) +end + +function Base.getindex(spec::SingleSpectrum, inds) + return Spectrum(spectral_axis(spec)[inds], flux_axis(spec)[inds], meta(spec)) +end + +Base.firstindex(spec::SingleSpectrum) = firstindex(spectral_axis(spec)) +Base.lastindex(spec::SingleSpectrum) = lastindex(spectral_axis(spec)) + +function Base.show(io::IO, spec::SingleSpectrum) + w = spectral_axis(spec) + f = flux_axis(spec) + println(io, "SingleSpectrum($(eltype(w)), $(eltype(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..7930d6a 100644 --- a/src/transforms/redden.jl +++ b/src/transforms/redden.jl @@ -1,13 +1,13 @@ 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) + s, f = spectral_axis(spec), flux_axis(spec) + @. f = DustExtinction.redden(law, s, f; Rv, Av) return spec end @@ -31,7 +31,8 @@ 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) + s, f = spectral_axis(spec), flux_axis(spec) + @. f = DustExtinction.deredden(law, s, f; Rv, Av) return spec end diff --git a/src/utils.jl b/src/utils.jl index 41a8598..919f921 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -21,16 +21,18 @@ 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}}) + 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) -Spectrum(Float64, Float64) - T: 6000 - name: Blackbody +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> 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 @@ -56,3 +58,24 @@ _blackbody(wave::AbstractVector{<:Quantity}, T::Quantity) = blackbody(T).(wave) Returns a function for calculating blackbody curves. """ blackbody(T::Quantity) = w->2h * c_0^2 / w^5 / (exp(h * c_0 / (w * k_B * T)) - 1) + +#""" +# equivalent_width(::AbstractSpectrum) +# +#Calculate the equivalent width of the given continuum-normalized spectrum. Return value has units equal to wavelengths. +#""" +#function equivalent_width(spec::AbstractSpectrum) +# dx = 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/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/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 9d9b503..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 + 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 03afb3b..b18af69 100644 --- a/test/spectrum.jl +++ b/test/spectrum.jl @@ -1,5 +1,25 @@ Random.seed!(8675309) +@testset "Spectrum types" begin + wave, flux = collect(1:5), 1:5 + @test spectrum(wave, flux) isa SingleSpectrum + wave[3] = 10 + @test_throws ArgumentError spectrum(wave, flux) + @test_throws ArgumentError spectrum(wave[begin:end-1], flux) + + wave, flux = repeat(1:5, 1, 3), rand(5, 3) + @test spectrum(wave, flux) isa EchelleSpectrum + wave[3, 3] = 10 + @test_throws ArgumentError spectrum(wave, flux) + @test_throws ArgumentError spectrum(wave[begin:end-1, :], flux) + + wave, flux = collect(1:5), rand(5, 4, 3) + @test spectrum(wave, flux) isa IFUSpectrum + wave[3] = 10 + @test_throws ArgumentError spectrum(wave, flux) + @test_throws ArgumentError spectrum(wave[begin:end-1], flux) +end + @testset "Spectrum - Single" begin wave = range(1e4, 5e4, length = 1000) sigma = randn(size(wave)) @@ -10,12 +30,15 @@ Random.seed!(8675309) flux[134] = 1 Âą 0.1 spec = spectrum(wave, flux, name = "test spectrum") - - @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 + spec_indexed = spec[begin:end] + + @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(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 @@ -24,30 +47,33 @@ Random.seed!(8675309) @test argmin(spec) == 134 @test findmax(spec) == (1000 Âą 1, 7) @test findmin(spec) == (1 Âą 0.1, 134) - @test spec.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 AssertionError spectrum(wave, flux_trimmed) + @test_throws ArgumentError spectrum(wave, flux_trimmed) expected = """ - Spectrum(Float64, Measurements.Measurement{Float64}) - name: test spectrum""" + SingleSpectrum(Float64, Measurements.Measurement{Float64}) + 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" 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") @@ -58,35 +84,62 @@ 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 size(spec) == (n_orders, n_wavs) - @test length(spec) == n_orders * n_wavs + @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 @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.(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 - name: Test Echelle Spectrum""" + 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" 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) + 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) == (: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 flux_axis(spec) == 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)) @@ -99,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 @@ -121,27 +174,29 @@ 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 = """ - 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}}) + 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 @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" @@ -151,20 +206,20 @@ 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 size(spec) === (n_orders, n_wavs) + @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 @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 @@ -173,14 +228,15 @@ 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 - sprint(show, spec) expected = """ EchelleSpectrum(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Unitful.Quantity{Measurements.Measurement{Float64}, 𝐌 𝐋^-1 𝐓^-3, Unitful.FreeUnits{(Å^-1, m^-2, W), 𝐌 𝐋^-1 𝐓^-3, nothing}}) # orders: 3 - name: test echelle""" + 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 @@ -194,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..e899d48 100644 --- a/test/transforms/redden.jl +++ b/test/transforms/redden.jl @@ -31,15 +31,16 @@ 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 + 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 @test_throws MethodError redden(spec, Av, law = sin) @@ -49,11 +50,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 +62,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 index 46b1ba3..fa2f7e3 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -1,4 +1,4 @@ -using Spectra: blackbody +using Spectra: blackbody#, line_flux, equivalent_width @testset "Blackbody T=$T" for T in [2000, 4000, 6000] wave = range(1e3, 5e4, length = 1000) @@ -7,7 +7,7 @@ using Spectra: blackbody bb = @inferred blackbody(wave, T) @test typeof(bb) <: Spectra.Spectrum @test bb.T == T - @test bb.wave[argmax(bb)] ≈ b / T rtol = 0.01 + @test spectral_axis(bb)[argmax(bb)] ≈ b / T rtol = 0.01 wave *= u"angstrom" T *= u"K" @@ -15,5 +15,15 @@ using Spectra: blackbody @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 + @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