Skip to content

Add multi-wavelength spectroscopy tutorial#203

Open
aditya-pandey-dev wants to merge 14 commits intoJuliaAstro:mainfrom
aditya-pandey-dev:add-spectroscopy-tutorial
Open

Add multi-wavelength spectroscopy tutorial#203
aditya-pandey-dev wants to merge 14 commits intoJuliaAstro:mainfrom
aditya-pandey-dev:add-spectroscopy-tutorial

Conversation

@aditya-pandey-dev
Copy link
Copy Markdown

@aditya-pandey-dev aditya-pandey-dev commented Mar 20, 2026

Closes #202

What's included

A new tutorial at docs/src/tutorials/multi-wavelength-spectroscopy.md covering a realistic end-to-end spectroscopy workflow:

  1. Loading a real SDSS DR14 spectrum from a FITS binary table (FITSIO.jl)
  2. Attaching physical units — Å, erg/s/cm²/Å (Unitful + UnitfulAstro)
  3. Fλ ↔ Fν conversion (Jansky)
  4. Dust extinction with CCM89, OD94, CAL00 laws (DustExtinction.jl)
  5. Spectral axis inspection via spectral_axis / flux_axis (Spectra.jl)
  6. Automatic uncertainty propagation (Measurements.jl)
  7. Synthetic blackbody spectra at 4 temperatures
  8. Spectral arithmetic — sky subtraction, scaling
  9. Cosmological redshift + luminosity distances (Cosmology.jl)

Testing

Includes test/test_tutorial.jl — a standalone verification script:

julia --project=docs/ test/test_tutorial.jl

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

  • Code blocks use plain julia fences (not @example) since Spectra.jl
    is installed from GitHub and not yet in the General registry
  • continuum() is referenced as a planned future feature with a link to
    Towards SpectraBase.jl [name tbd] Spectra.jl#41
  • [compat] bounds added for all new registered dependencies

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.
@aditya-pandey-dev
Copy link
Copy Markdown
Author

@icweaver can you check the PR. What changes should I do?

Copy link
Copy Markdown
Member

@cgarling cgarling left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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 @example here? 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 blackbody returns 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 the blackbody docstring 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)
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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⁻² Å⁻¹)",
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
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⁻² Å⁻¹)",
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
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",
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
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⁻² Å⁻¹)",
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
xlabel = "Wavelength (Å)", ylabel = "Flux (erg s⁻¹ cm⁻² Å⁻¹)",
xlabel = "Wavelength (Å)", ylabel = "Flux Density (erg s⁻¹ cm⁻² Å⁻¹)",

Comment on lines +201 to +213
## 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])
```
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Comment on lines +137 to +138
wave_bb = range(3000, 10000, length = 500)
bb = blackbody(wave_bb, 10_000)
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would add units to these


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
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Comment on lines +219 to +224
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)
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would put units on wave_synth and the temperature for blackbody(w, T)

aditya-pandey-dev added 7 commits March 31, 2026 14:51
- @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
@aditya-pandey-dev
Copy link
Copy Markdown
Author

@cgarling I've added Spectra.jl to docs/Project.toml via
[sources]. However, @example blocks still can't run on CI
because the tutorial makes network calls:

  1. Downloads a real SDSS FITS file (~50 MB)
  2. SFD98Map() downloads dust map tiles (~120 MB)
    So i think there are two suitable options : (a) Keep plain julia blocks with a note explaining network
    dependency (b) Split into @example (synthetic/offline parts) + plain
    julia (network parts)
    Which approach would you prefer?

@cgarling
Copy link
Copy Markdown
Member

cgarling commented Apr 8, 2026

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)?

Copy link
Copy Markdown
Member

@cgarling cgarling left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

duplicated above

Suggested change
using Spectra

Comment on lines +66 to +75
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
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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))
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
ebv = dustmap(rad2deg(gal.l), rad2deg(gal.b))
ebv = dustmap(gal.l, gal.b)

dustmap expects inputs in radians

Comment on lines +112 to +124
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
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
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

Comment on lines +148 to +166
# 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
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# 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.

Comment on lines +182 to +195
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
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

2 participants