Add multi-wavelength spectroscopy tutorial#203
Add multi-wavelength spectroscopy tutorial#203aditya-pandey-dev wants to merge 14 commits intoJuliaAstro:mainfrom
Conversation
Adds an end-to-end spectroscopy tutorial covering: - Loading real SDSS spectra from FITS files (FITSIO.jl) - Physical units and Fλ↔Fν conversions (Unitful + UnitfulAstro) - Dust extinction with CCM89, OD94, CAL00 laws (DustExtinction.jl) - Spectral axis inspection via spectral_axis/flux_axis (Spectra.jl) - Automatic uncertainty propagation (Measurements.jl) - Synthetic blackbody spectra generation (Spectra.jl) - Spectral arithmetic: sky subtraction and scaling - Cosmological redshift and luminosity distances (Cosmology.jl) Includes test/test_tutorial.jl — standalone verification script. Updates docs/Project.toml with [compat] bounds for all new deps. Closes JuliaAstro#202
DustExtinction='0.6, 1' excluded v0.11.x (current release), creating an unsolvable constraint with FITSIO 0.17 in CI.
Tutorial uses plain julia code blocks (not @example), so docs build environment does not need to install tutorial dependencies. Extra packages were causing CI resolution failures in the JuliaAstroDocs test suite.
The CI test runner picks up all .jl files in test/ and tries to execute them. test_tutorial.jl requires packages not present in the test environment, causing JuliaAstroDocs to error. Users can run the tutorial locally by installing the packages listed in the tutorial's Packages section.
|
@icweaver can you check the PR. What changes should I do? |
cgarling
left a comment
There was a problem hiding this comment.
Interesting example, thanks for the PR. I made some specific comments/suggestions below, here is some high-level feedback:
- It would be nice to have this run on CI so the output can be displayed. Could you expand on why you can't use
@examplehere? I thought that if you setup the doc build environment correctly there shouldn't be a problem with importing Spectra.jl even if it isn't in the General registry yet (i.e., using [sources] to set the url in the Project.toml). - Mostly you are dealing with flux density units in this example so we should make the associated plot labels "Flux Density" rather than "Flux" for clarity
- I would try to be more consistent with using units when constructing the wavelength / flux density vectors, some places where you are constructing blackbody spectra you are not putting units on the wavelengths/temperatures and I think we want to encourage the use of units
- I think that
blackbodyreturns radiance units so you may need to calculate the luminosity density as L_lam(lambda) = 4 \pi^2 R^2 B(lambda,T) where R is the radius of the spherical blackbody (say 1 solar radius) and B(lambda, T) is what blackbody returns. It would be good to check this and amend theblackbodydocstring to be more descriptive. - We are mostly using Makie.jl for plotting in our documentation now, idk if anyone has strong feelings on this. I don't mind using Plots.jl but would probably prefer makie
|
|
||
| ### Dereddening the SDSS spectrum | ||
| ```julia | ||
| spec_dered = deredden(spec, 0.1) |
There was a problem hiding this comment.
It would be nice to use DustExtinction.SFD98Map here to get the E(B-V) value for the object here using its coordinates if they are saved in the file (then Av = R_v * E(B-V)) . Depending on the format of the coordinates, you might need to use SkyCoords.jl to convert them.
| ws = ustrip.(spectral_axis(spec)) | ||
| plot(ws, ustrip.(flux_axis(spec)), label = "Observed", lw = 0.5, color = :gray, alpha = 0.6) | ||
| plot!(ws, ustrip.(flux_axis(spec_dered)), label = "Dereddened", lw = 0.5, color = :steelblue) | ||
| plot!(xlabel = "Wavelength (Å)", ylabel = "Flux (erg s⁻¹ cm⁻² Å⁻¹)", |
There was a problem hiding this comment.
| plot!(xlabel = "Wavelength (Å)", ylabel = "Flux (erg s⁻¹ cm⁻² Å⁻¹)", | |
| plot!(xlabel = "Wavelength (Å)", ylabel = "Flux Density (erg s⁻¹ cm⁻² Å⁻¹)", |
|
|
||
| plot(ustrip.(waves), ustrip.(fluxes), | ||
| xlabel = "Wavelength (Å)", | ||
| ylabel = "Flux (erg s⁻¹ cm⁻² Å⁻¹)", |
There was a problem hiding this comment.
| ylabel = "Flux (erg s⁻¹ cm⁻² Å⁻¹)", | |
| ylabel = "Flux Density (erg s⁻¹ cm⁻² Å⁻¹)", |
| plot!(spectral_axis(bb), ustrip.(flux_axis(bb_05)), label = "Av = 0.5", lw = 2) | ||
| plot!(spectral_axis(bb), ustrip.(flux_axis(bb_15)), label = "Av = 1.5", lw = 2) | ||
| plot!(spectral_axis(bb), ustrip.(flux_axis(bb_30)), label = "Av = 3.0", lw = 2) | ||
| plot!(xlabel = "Wavelength (Å)", ylabel = "Flux", |
There was a problem hiding this comment.
| plot!(xlabel = "Wavelength (Å)", ylabel = "Flux", | |
| plot!(xlabel = "Wavelength (Å)", ylabel = "Flux Density", |
| # Zoom around Hα (6563 Å) | ||
| ha_mask = (ustrip.(waves) .> 6400) .& (ustrip.(waves) .< 6700) | ||
| plot(ustrip.(waves[ha_mask]), ustrip.(fluxes[ha_mask]), | ||
| xlabel = "Wavelength (Å)", ylabel = "Flux (erg s⁻¹ cm⁻² Å⁻¹)", |
There was a problem hiding this comment.
| xlabel = "Wavelength (Å)", ylabel = "Flux (erg s⁻¹ cm⁻² Å⁻¹)", | |
| xlabel = "Wavelength (Å)", ylabel = "Flux Density (erg s⁻¹ cm⁻² Å⁻¹)", |
| ## Part 5: Spectra with uncertainties | ||
| ```julia | ||
| using Measurements | ||
|
|
||
| wave_m = range(4000, 7000, length = 200) * u"angstrom" | ||
| flux_m = ((5e-17 .± 0.5e-17) .* ones(200)) * u"erg/s/cm^2/angstrom" | ||
| spec_m = spectrum(wave_m, flux_m) | ||
|
|
||
| spec_red = redden(spec_m, 0.5) | ||
|
|
||
| println("Original flux[1] : ", flux_axis(spec_m)[1]) | ||
| println("Reddened flux[1] : ", flux_axis(spec_red)[1]) | ||
| ``` |
There was a problem hiding this comment.
If this SDSS spectrum contains an error column (which I expect it does) I would show use of Measurements with that data as it's a more practical example.
| wave_bb = range(3000, 10000, length = 500) | ||
| bb = blackbody(wave_bb, 10_000) |
|
|
||
| for (z, c) in zip([0.0, 0.5, 1.0, 2.0], [:black, :blue, :green, :red]) | ||
| wave_obs = wave_rest .* (1 + z) | ||
| flux_obs = flux_axis(bb_rest) ./ (1 + z)^2 |
There was a problem hiding this comment.
Should this be (1 + z) / (4 \pi luminosity_dist(cosmo, z)^2)? Also note that I think the blackbody is in radiance units so you may need to calculate the luminosity density as L_lam(lambda) = 4 \pi^2 R^2 B(lambda,T) where R is the radius of the spherical blackbody (say 1 solar radius) and B(lambda, T) is what blackbody returns.
| wave_synth = range(1000, 20000, length = 1000) | ||
|
|
||
| bb_3k = blackbody(wave_synth, 3_000) | ||
| bb_6k = blackbody(wave_synth, 6_000) | ||
| bb_10k = blackbody(wave_synth, 10_000) | ||
| bb_30k = blackbody(wave_synth, 30_000) |
There was a problem hiding this comment.
I would put units on wave_synth and the temperature for blackbody(w, T)
- @example blocks for CI + Spectra via [sources] in Project.toml - Flux -> Flux Density in all plot ylabels - Units on wave_bb, wave_synth (.* u"angstrom"), temperatures (u"K") - L_lam = 4pi^2 * R^2 * B(lambda,T) correct luminosity formula - SDSS ivar column for real Measurements uncertainties - F_obs = L*(1+z)/(4pi*d_L^2) correct cosmological formula - Plots.jl -> CairoMakie throughout - SFD98Map with ICRSCoords -> GalCoords conversion
- Replace synthetic galaxy spectrum with real SDSS plate 1323 data - Use FITS ivar column for real per-pixel σ = 1/√ivar uncertainties - Addresses mentor comment: use Measurements with actual SDSS error column
|
@cgarling I've added Spectra.jl to docs/Project.toml via
|
|
I do Github Actions with network dependencies pretty frequently (e.g., through DataDeps.jl) and can't recall any issues. What exactly is the problem CI is creating, are you hitting a firewall (Github won't let you access external URLs)? |
cgarling
left a comment
There was a problem hiding this comment.
Thanks for the update, good progress! More comments below
| of 10⁻¹⁷ erg s⁻¹ cm⁻² Å⁻¹ and an `ivar` (inverse-variance) column | ||
| that encodes real per-pixel measurement uncertainties. | ||
| ```julia | ||
| using Downloads, FITSIO, Unitful, UnitfulAstro, Measurements |
There was a problem hiding this comment.
Some of these are duplicated in the ## Setup block above. I would probably just put all the imports in the first code block rather than doing the setup separately but it's a stylistic choice
| # Attach units + Measurements.jl uncertainties | ||
| flux_meas = ((flux_raw .± sigma_raw) .* 1e-17) * u"erg/s/cm^2/angstrom" | ||
|
|
||
| using Spectra |
There was a problem hiding this comment.
duplicated above
| using Spectra |
| wv = ustrip.(wave_aa) | ||
| fv = Measurements.value.(ustrip.(flux_meas)) | ||
| fe = Measurements.uncertainty.(ustrip.(flux_meas)) | ||
|
|
||
| fig, ax = lines(wv, fv; | ||
| axis = (xlabel = "Wavelength (Å)", | ||
| ylabel = "Flux Density (erg s⁻¹ cm⁻² Å⁻¹)", | ||
| title = "SDSS Galaxy — plate 1323, fiber 12")) | ||
| band!(ax, wv, fv .- fe, fv .+ fe; color = (:steelblue, 0.3)) | ||
| fig |
There was a problem hiding this comment.
CairoMakie should be able to take types from Unitful and Measurements directly, you might not have to do all the ustrip shenanigans, I would just try fig, ax = lines(wave_aa, flux_meas; ...).
| gal = convert(GalCoords, eq) | ||
|
|
||
| dustmap = SFD98Map() | ||
| ebv = dustmap(rad2deg(gal.l), rad2deg(gal.b)) |
There was a problem hiding this comment.
| ebv = dustmap(rad2deg(gal.l), rad2deg(gal.b)) | |
| ebv = dustmap(gal.l, gal.b) |
dustmap expects inputs in radians
| tau = [ustrip(ext(w * u"angstrom")) * Av for w in wave_raw] | ||
| corr = 10 .^ (tau ./ 2.5) | ||
|
|
||
| fv_raw = Measurements.value.(ustrip.(spec.flux)) | ||
| fv_ded = fv_raw .* corr | ||
|
|
||
| fig2, ax2 = lines(wv, fv_raw; label = "Raw", | ||
| axis = (xlabel = "Wavelength (Å)", | ||
| ylabel = "Flux Density (erg s⁻¹ cm⁻² Å⁻¹)", | ||
| title = "CCM89 Dust Dereddening")) | ||
| lines!(ax2, wv, fv_ded; label = "Dereddened", color = :orange) | ||
| axislegend(ax2) | ||
| fig2 |
There was a problem hiding this comment.
| tau = [ustrip(ext(w * u"angstrom")) * Av for w in wave_raw] | |
| corr = 10 .^ (tau ./ 2.5) | |
| fv_raw = Measurements.value.(ustrip.(spec.flux)) | |
| fv_ded = fv_raw .* corr | |
| fig2, ax2 = lines(wv, fv_raw; label = "Raw", | |
| axis = (xlabel = "Wavelength (Å)", | |
| ylabel = "Flux Density (erg s⁻¹ cm⁻² Å⁻¹)", | |
| title = "CCM89 Dust Dereddening")) | |
| lines!(ax2, wv, fv_ded; label = "Dereddened", color = :orange) | |
| axislegend(ax2) | |
| fig2 | |
| A_lambda = ext.(wave_raw) .* Av | |
| corr = exp10.(-A_lambda ./ 2.5) | |
| flux_dered = flux_meas ./ corr | |
| fig2, ax2 = lines(wave_aa, flux_meas; label = "Raw", | |
| axis = (title = "CCM89 Dust Dereddening",)) | |
| lines!(ax2, wave_aa, flux_dered; label = "Dereddened", color = :orange) | |
| axislegend(ax2) | |
| fig2 |
CCM89 returns A_\lambda / Av for input wavelength \lambda so you just multiply by Av to get the extinction in magnitudes at the input wavelength. Then corr here will be <1 by definition and you can deredden the spectrum by dividing the observed spectrum by the expected attenuation corr. When you deredden a spectrum you are trying to infer what its intrinsic luminosity is given an observed luminosity, so you must have flux_dered >= flux_raw
| # Compare multiple temperatures — wave_synth with units | ||
| wave_synth = collect(range(1000, 20000, length = 1000)) .* u"angstrom" | ||
| temps = [3_000u"K", 6_000u"K", 10_000u"K", 30_000u"K"] | ||
| labels_bb = ["3 000 K", "6 000 K", "10 000 K", "30 000 K"] | ||
| colors_bb = [:red, :orange, :steelblue, :purple] | ||
|
|
||
| fig3, ax3 = lines(ustrip.(wave_synth), | ||
| ustrip.(blackbody(wave_synth, temps[1]).flux); | ||
| label = labels_bb[1], color = colors_bb[1], | ||
| axis = (xlabel = "Wavelength (Å)", | ||
| ylabel = "Flux Density (erg s⁻¹ cm⁻² Å⁻¹ sr⁻¹)", | ||
| title = "Blackbody Spectra")) | ||
| for (T, lab, col) in zip(temps[2:end], labels_bb[2:end], colors_bb[2:end]) | ||
| lines!(ax3, ustrip.(wave_synth), | ||
| ustrip.(blackbody(wave_synth, T).flux); | ||
| label = lab, color = col) | ||
| end | ||
| axislegend(ax3, position = :rt) | ||
| fig3 |
There was a problem hiding this comment.
| # Compare multiple temperatures — wave_synth with units | |
| wave_synth = collect(range(1000, 20000, length = 1000)) .* u"angstrom" | |
| temps = [3_000u"K", 6_000u"K", 10_000u"K", 30_000u"K"] | |
| labels_bb = ["3 000 K", "6 000 K", "10 000 K", "30 000 K"] | |
| colors_bb = [:red, :orange, :steelblue, :purple] | |
| fig3, ax3 = lines(ustrip.(wave_synth), | |
| ustrip.(blackbody(wave_synth, temps[1]).flux); | |
| label = labels_bb[1], color = colors_bb[1], | |
| axis = (xlabel = "Wavelength (Å)", | |
| ylabel = "Flux Density (erg s⁻¹ cm⁻² Å⁻¹ sr⁻¹)", | |
| title = "Blackbody Spectra")) | |
| for (T, lab, col) in zip(temps[2:end], labels_bb[2:end], colors_bb[2:end]) | |
| lines!(ax3, ustrip.(wave_synth), | |
| ustrip.(blackbody(wave_synth, T).flux); | |
| label = lab, color = col) | |
| end | |
| axislegend(ax3, position = :rt) | |
| fig3 | |
| # Compare multiple temperatures — wave_synth with units | |
| wave_synth = collect(range(1000, 20000, length = 1000)) .* u"angstrom" | |
| temps = [3_000u"K", 6_000u"K", 10_000u"K", 30_000u"K"] | |
| labels_bb = ["3 000 K", "6 000 K", "10 000 K", "30 000 K"] | |
| colors_bb = [:red, :orange, :steelblue, :purple] | |
| fig3, ax3 = lines(wave_synth, | |
| blackbody(wave_synth, temps[1]).flux; | |
| label = labels_bb[1], color = colors_bb[1], | |
| axis = (title = "Blackbody Spectra", | |
| yscale = log10, yticks = LogTicks(LinearTicks(5)))) | |
| for (T, lab, col) in zip(temps[2:end], labels_bb[2:end], colors_bb[2:end]) | |
| lines!(ax3, wave_synth, | |
| blackbody(wave_synth, T).flux; | |
| label = lab, color = col) | |
| end | |
| axislegend(ax3, position = :rt) | |
| fig3 |
Here's an example showing how to let Makie.jl handle plotting quantities with units from Unitful.jl. I also made the y-axis a logscale to better distinguish between them.
| fig4, ax4 = lines(ustrip.(wave_rest), ustrip.(bb_rest.flux); | ||
| label = "z = 0", color = :black, | ||
| axis = (xlabel = "Wavelength (Å)", | ||
| ylabel = "Flux Density (erg s⁻¹ cm⁻² Å⁻¹)", | ||
| title = "Cosmological Surface-Brightness Dimming")) | ||
|
|
||
| for (z, col) in zip([0.5, 1.0, 2.0], [:blue, :green, :red]) | ||
| d_L = uconvert(u"cm", luminosity_dist(cosmo, z)) | ||
| w_obs = wave_rest .* (1 + z) | ||
| F_obs = L_lam .* (1 + z) ./ (4π .* d_L.^2) | ||
| lines!(ax4, ustrip.(w_obs), ustrip.(F_obs); label = "z = $z", color = col) | ||
| end | ||
| axislegend(ax4, position = :rt) | ||
| fig4 |
There was a problem hiding this comment.
| fig4, ax4 = lines(ustrip.(wave_rest), ustrip.(bb_rest.flux); | |
| label = "z = 0", color = :black, | |
| axis = (xlabel = "Wavelength (Å)", | |
| ylabel = "Flux Density (erg s⁻¹ cm⁻² Å⁻¹)", | |
| title = "Cosmological Surface-Brightness Dimming")) | |
| for (z, col) in zip([0.5, 1.0, 2.0], [:blue, :green, :red]) | |
| d_L = uconvert(u"cm", luminosity_dist(cosmo, z)) | |
| w_obs = wave_rest .* (1 + z) | |
| F_obs = L_lam .* (1 + z) ./ (4π .* d_L.^2) | |
| lines!(ax4, ustrip.(w_obs), ustrip.(F_obs); label = "z = $z", color = col) | |
| end | |
| axislegend(ax4, position = :rt) | |
| fig4 | |
| fig4, ax4 = lines(ustrip.(wave_rest), ustrip.(bb_rest.flux); | |
| label = "z = 0", color = :black, | |
| axis = (xlabel = "Wavelength (Å)", | |
| ylabel = "Flux Density (erg s⁻¹ cm⁻² Å⁻¹)", | |
| yscale = log10, yticks=LogTicks(LinearTicks(5)), | |
| title = "Cosmological Surface-Brightness Dimming")) | |
| for (z, col) in zip([0.1, 0.5, 1.0], [:blue, :green, :red]) | |
| d_L = uconvert(u"cm", luminosity_dist(cosmo, z)) | |
| w_obs = wave_rest .* (1 + z) | |
| F_obs = L_lam .* (1 + z) ./ (4π .* d_L.^2) | |
| lines!(ax4, ustrip.(w_obs), ustrip.(F_obs); label = "z = $z", color = col) | |
| end | |
| axislegend(ax4, position = :rt) | |
| fig4 |
Changed to display a little better
Closes #202
What's included
A new tutorial at
docs/src/tutorials/multi-wavelength-spectroscopy.mdcovering a realistic end-to-end spectroscopy workflow:spectral_axis/flux_axis(Spectra.jl)Testing
Includes
test/test_tutorial.jl— a standalone verification script:All 5 dependency tests pass on Julia 1.12.5
Packages exercised
Spectra.jl, FITSIO.jl, DustExtinction.jl, UnitfulAstro.jl, Measurements.jl, Cosmology.jl, Plots.jl
Notes
juliafences (not@example) since Spectra.jlis installed from GitHub and not yet in the General registry
continuum()is referenced as a planned future feature with a link toTowards
SpectraBase.jl[name tbd] Spectra.jl#41[compat]bounds added for all new registered dependencies