|
| 1 | +# [Multi-wavelength spectroscopy](@id tutorial-multiwavelength-spectroscopy) |
| 2 | + |
| 3 | +Spectroscopy is one of the most powerful tools in an astronomer's toolkit. By splitting |
| 4 | +light into its component wavelengths, we can learn about a source's chemical composition, |
| 5 | +temperature, velocity, and much more. |
| 6 | + |
| 7 | +In this tutorial we walk through a realistic workflow: downloading a real SDSS optical |
| 8 | +spectrum, loading it with [`FITSIO.jl`](https://juliaastro.org/FITSIO/stable/), wrapping |
| 9 | +it in a `Spectrum` object with physical units via |
| 10 | +[`UnitfulAstro.jl`](https://juliaastro.org/UnitfulAstro/stable/), applying dust corrections |
| 11 | +with [`DustExtinction.jl`](https://juliaastro.org/DustExtinction/stable/), propagating |
| 12 | +uncertainties with [`Measurements.jl`](https://juliaphysics.github.io/Measurements.jl/stable/), |
| 13 | +and computing cosmological distances with |
| 14 | +[`Cosmology.jl`](https://juliaastro.org/Cosmology/stable/). |
| 15 | + |
| 16 | +!!! note |
| 17 | + This tutorial uses Julia 1.9+. `Spectra.jl` is under active development and installed |
| 18 | + from GitHub. Code blocks show the intended API; some call signatures may evolve. |
| 19 | + |
| 20 | +## Packages |
| 21 | +``` |
| 22 | +pkg> add FITSIO Unitful UnitfulAstro DustExtinction Measurements Cosmology Plots |
| 23 | +pkg> add https://github.com/JuliaAstro/Spectra.jl |
| 24 | +``` |
| 25 | + |
| 26 | +| Package | Role | |
| 27 | +|---|---| |
| 28 | +| `FITSIO` | Read spectral data from FITS binary tables | |
| 29 | +| `Spectra` | Core `Spectrum` type, blackbody, reddening | |
| 30 | +| `DustExtinction` | Empirical extinction laws (CCM89, OD94, CAL00) | |
| 31 | +| `Unitful` / `UnitfulAstro` | Physical units — Å, nm, Jy, erg/s/cm²/Å | |
| 32 | +| `Measurements` | Automatic uncertainty propagation | |
| 33 | +| `Cosmology` | ΛCDM luminosity distances | |
| 34 | +| `Plots` | Visualisation | |
| 35 | + |
| 36 | +--- |
| 37 | + |
| 38 | +## Part 1: Loading a real SDSS spectrum |
| 39 | + |
| 40 | +We work with a publicly available spectrum from SDSS DR14 — a galaxy observed on |
| 41 | +plate 1323, MJD 52797, fiber 12. The FITS file stores the coadded spectrum as a |
| 42 | +binary table with log₁₀-spaced wavelengths and calibrated flux in units of |
| 43 | +10⁻¹⁷ erg/s/cm²/Å. |
| 44 | + |
| 45 | +### Downloading the data |
| 46 | +```julia |
| 47 | +using Downloads |
| 48 | + |
| 49 | +sdss_url = "https://dr14.sdss.org/optical/spectrum/view/data/format=fits/spec=lite" * |
| 50 | + "?plateid=1323&mjd=52797&fiberid=12" |
| 51 | +sdss_file = joinpath(@__DIR__, "sdss_example.fits") |
| 52 | + |
| 53 | +if !isfile(sdss_file) |
| 54 | + Downloads.download(sdss_url, sdss_file) |
| 55 | +end |
| 56 | +``` |
| 57 | + |
| 58 | +### Reading the FITS file |
| 59 | +```julia |
| 60 | +using FITSIO, Unitful, UnitfulAstro |
| 61 | + |
| 62 | +f = FITS(sdss_file) |
| 63 | +loglam = read(f[2], "loglam") |
| 64 | +flux_raw = read(f[2], "flux") |
| 65 | +ivar = read(f[2], "ivar") |
| 66 | + |
| 67 | +wave = (10 .^ loglam) * u"angstrom" |
| 68 | +flux = (flux_raw .* 1e-17) * u"erg/s/cm^2/angstrom" |
| 69 | +close(f) |
| 70 | + |
| 71 | +println("Wavelength range : ", minimum(wave), " — ", maximum(wave)) |
| 72 | +println("Number of pixels : ", length(wave)) |
| 73 | +``` |
| 74 | + |
| 75 | +### Creating a Spectrum object |
| 76 | +```julia |
| 77 | +using Spectra |
| 78 | + |
| 79 | +spec = spectrum(wave, flux) |
| 80 | +println(spec) |
| 81 | +``` |
| 82 | + |
| 83 | +### Plotting the raw spectrum |
| 84 | +```julia |
| 85 | +using Plots |
| 86 | + |
| 87 | +waves = spectral_axis(spec) |
| 88 | +fluxes = flux_axis(spec) |
| 89 | + |
| 90 | +plot(ustrip.(waves), ustrip.(fluxes), |
| 91 | + xlabel = "Wavelength (Å)", |
| 92 | + ylabel = "Flux (erg s⁻¹ cm⁻² Å⁻¹)", |
| 93 | + title = "SDSS Galaxy — plate 1323, fiber 12", |
| 94 | + legend = false, lw = 0.5, color = :steelblue, size = (800, 400)) |
| 95 | +``` |
| 96 | + |
| 97 | +--- |
| 98 | + |
| 99 | +## Part 2: Working with physical units |
| 100 | + |
| 101 | +### Wavelength conversions |
| 102 | +```julia |
| 103 | +ha_wave = 6563.0u"angstrom" |
| 104 | + |
| 105 | +ha_nm = uconvert(u"nm", ha_wave) |
| 106 | +ha_um = uconvert(u"μm", ha_wave) |
| 107 | + |
| 108 | +println("Hα = $ha_wave = $ha_nm = $ha_um") |
| 109 | +``` |
| 110 | + |
| 111 | +### Converting between Fλ and Fν |
| 112 | + |
| 113 | +Radio and infrared data are usually stored as Fν (flux per unit frequency, often in |
| 114 | +Jansky). The conversion is Fν = Fλ × λ² / c. |
| 115 | +```julia |
| 116 | +using Unitful: c0 |
| 117 | + |
| 118 | +Flambda = 2.5e-17u"erg/s/cm^2/angstrom" |
| 119 | +Fnu = Flambda * ha_wave^2 / c0 |> u"erg/s/cm^2/Hz" |
| 120 | +Fnu_jy = uconvert(u"Jy", Fnu) |
| 121 | + |
| 122 | +println("Fλ = $Flambda") |
| 123 | +println("Fν = $Fnu = $Fnu_jy") |
| 124 | +``` |
| 125 | + |
| 126 | +--- |
| 127 | + |
| 128 | +## Part 3: Dust extinction and dereddening |
| 129 | + |
| 130 | +Interstellar dust makes sources appear fainter and redder — it preferentially absorbs |
| 131 | +blue photons. `Spectra.jl` integrates with `DustExtinction.jl` for one-line corrections. |
| 132 | + |
| 133 | +### Reddening a blackbody |
| 134 | +```julia |
| 135 | +using DustExtinction |
| 136 | + |
| 137 | +wave_bb = range(3000, 10000, length = 500) |
| 138 | +bb = blackbody(wave_bb, 10_000) |
| 139 | + |
| 140 | +bb_05 = redden(bb, 0.5) |
| 141 | +bb_15 = redden(bb, 1.5) |
| 142 | +bb_30 = redden(bb, 3.0) |
| 143 | + |
| 144 | +plot(spectral_axis(bb), ustrip.(flux_axis(bb)), label = "Av = 0", lw = 2) |
| 145 | +plot!(spectral_axis(bb), ustrip.(flux_axis(bb_05)), label = "Av = 0.5", lw = 2) |
| 146 | +plot!(spectral_axis(bb), ustrip.(flux_axis(bb_15)), label = "Av = 1.5", lw = 2) |
| 147 | +plot!(spectral_axis(bb), ustrip.(flux_axis(bb_30)), label = "Av = 3.0", lw = 2) |
| 148 | +plot!(xlabel = "Wavelength (Å)", ylabel = "Flux", |
| 149 | + title = "Effect of dust extinction (T = 10 000 K blackbody)", |
| 150 | + size = (800, 400)) |
| 151 | +``` |
| 152 | + |
| 153 | +### Comparing extinction laws |
| 154 | +```julia |
| 155 | +bb_ccm = redden(bb, 1.0, law = CCM89) |
| 156 | +bb_od94 = redden(bb, 1.0, law = OD94) |
| 157 | +bb_cal = redden(bb, 1.0, law = CAL00) |
| 158 | +``` |
| 159 | + |
| 160 | +### Dereddening the SDSS spectrum |
| 161 | +```julia |
| 162 | +spec_dered = deredden(spec, 0.1) |
| 163 | + |
| 164 | +ws = ustrip.(spectral_axis(spec)) |
| 165 | +plot(ws, ustrip.(flux_axis(spec)), label = "Observed", lw = 0.5, color = :gray, alpha = 0.6) |
| 166 | +plot!(ws, ustrip.(flux_axis(spec_dered)), label = "Dereddened", lw = 0.5, color = :steelblue) |
| 167 | +plot!(xlabel = "Wavelength (Å)", ylabel = "Flux (erg s⁻¹ cm⁻² Å⁻¹)", |
| 168 | + title = "Dereddening the SDSS spectrum", size = (800, 400)) |
| 169 | +``` |
| 170 | + |
| 171 | +--- |
| 172 | + |
| 173 | +## Part 4: Inspecting spectral axes |
| 174 | + |
| 175 | +`Spectra.jl` provides `spectral_axis` and `flux_axis` helpers to extract the wavelength |
| 176 | +and flux arrays from any `Spectrum` object. This is useful for passing data to external |
| 177 | +fitting routines. |
| 178 | +```julia |
| 179 | +waves = spectral_axis(spec_dered) |
| 180 | +fluxes = flux_axis(spec_dered) |
| 181 | + |
| 182 | +println("First wavelength : ", waves[1]) |
| 183 | +println("Last wavelength : ", waves[end]) |
| 184 | + |
| 185 | +# Zoom around Hα (6563 Å) |
| 186 | +ha_mask = (ustrip.(waves) .> 6400) .& (ustrip.(waves) .< 6700) |
| 187 | +plot(ustrip.(waves[ha_mask]), ustrip.(fluxes[ha_mask]), |
| 188 | + xlabel = "Wavelength (Å)", ylabel = "Flux (erg s⁻¹ cm⁻² Å⁻¹)", |
| 189 | + title = "SDSS spectrum near Hα (6563 Å)", |
| 190 | + lw = 1.5, color = :steelblue, legend = false, size = (800, 400)) |
| 191 | +``` |
| 192 | + |
| 193 | +!!! note |
| 194 | + Chebyshev continuum fitting (`continuum()`) is planned for a future `Spectra.jl` |
| 195 | + release — see [JuliaAstro/Spectra.jl#41](https://github.com/JuliaAstro/Spectra.jl/issues/41). |
| 196 | + For now, `SpectralFitting.jl` provides full continuum and line fitting with |
| 197 | + OGIP-compatible models. |
| 198 | + |
| 199 | +--- |
| 200 | + |
| 201 | +## Part 5: Spectra with uncertainties |
| 202 | +```julia |
| 203 | +using Measurements |
| 204 | + |
| 205 | +wave_m = range(4000, 7000, length = 200) * u"angstrom" |
| 206 | +flux_m = ((5e-17 .± 0.5e-17) .* ones(200)) * u"erg/s/cm^2/angstrom" |
| 207 | +spec_m = spectrum(wave_m, flux_m) |
| 208 | + |
| 209 | +spec_red = redden(spec_m, 0.5) |
| 210 | + |
| 211 | +println("Original flux[1] : ", flux_axis(spec_m)[1]) |
| 212 | +println("Reddened flux[1] : ", flux_axis(spec_red)[1]) |
| 213 | +``` |
| 214 | + |
| 215 | +--- |
| 216 | + |
| 217 | +## Part 6: Synthetic blackbody spectra |
| 218 | +```julia |
| 219 | +wave_synth = range(1000, 20000, length = 1000) |
| 220 | + |
| 221 | +bb_3k = blackbody(wave_synth, 3_000) |
| 222 | +bb_6k = blackbody(wave_synth, 6_000) |
| 223 | +bb_10k = blackbody(wave_synth, 10_000) |
| 224 | +bb_30k = blackbody(wave_synth, 30_000) |
| 225 | + |
| 226 | +plot(spectral_axis(bb_3k), ustrip.(flux_axis(bb_3k)), label = "3 000 K", lw = 2, color = :red) |
| 227 | +plot!(spectral_axis(bb_6k), ustrip.(flux_axis(bb_6k)), label = "6 000 K", lw = 2, color = :orange) |
| 228 | +plot!(spectral_axis(bb_10k), ustrip.(flux_axis(bb_10k)), label = "10 000 K", lw = 2, color = :steelblue) |
| 229 | +plot!(spectral_axis(bb_30k), ustrip.(flux_axis(bb_30k)), label = "30 000 K", lw = 2, color = :purple) |
| 230 | +plot!(xlabel = "Wavelength (Å)", ylabel = "Flux", |
| 231 | + title = "Planck blackbodies at four temperatures", size = (800, 400)) |
| 232 | +``` |
| 233 | + |
| 234 | +--- |
| 235 | + |
| 236 | +## Part 7: Spectral arithmetic |
| 237 | +```julia |
| 238 | +wave_a = range(4000, 8000, length = 500) |
| 239 | +source = blackbody(wave_a, 8_000) |
| 240 | +sky = blackbody(wave_a, 300) * 0.1 |
| 241 | + |
| 242 | +science = source - sky |
| 243 | + |
| 244 | +plot(spectral_axis(source), ustrip.(flux_axis(source)), label = "Source", lw = 2) |
| 245 | +plot!(spectral_axis(sky), ustrip.(flux_axis(sky)), label = "Sky", lw = 2) |
| 246 | +plot!(spectral_axis(science),ustrip.(flux_axis(science)), label = "Source − Sky", lw = 2, ls = :dash) |
| 247 | +plot!(xlabel = "Wavelength (Å)", ylabel = "Flux", |
| 248 | + title = "Spectral arithmetic: sky subtraction", size = (800, 400)) |
| 249 | +``` |
| 250 | + |
| 251 | +--- |
| 252 | + |
| 253 | +## Part 8: Cosmological redshift |
| 254 | +```julia |
| 255 | +using Cosmology |
| 256 | + |
| 257 | +cosmo = cosmology() |
| 258 | + |
| 259 | +for z in [0.1, 0.5, 1.0, 2.0] |
| 260 | + d = luminosity_dist(cosmo, z) |
| 261 | + println("z = $z → dL = ", round(typeof(d), d, digits = 0)) |
| 262 | +end |
| 263 | + |
| 264 | +wave_rest = range(3000, 8000, length = 500) |
| 265 | +bb_rest = blackbody(wave_rest, 6_000) |
| 266 | + |
| 267 | +p = plot(xlabel = "Observed wavelength (Å)", ylabel = "Flux", |
| 268 | + title = "Cosmological redshift", size = (800, 400)) |
| 269 | + |
| 270 | +for (z, c) in zip([0.0, 0.5, 1.0, 2.0], [:black, :blue, :green, :red]) |
| 271 | + wave_obs = wave_rest .* (1 + z) |
| 272 | + flux_obs = flux_axis(bb_rest) ./ (1 + z)^2 |
| 273 | + spec_z = spectrum(wave_obs, flux_obs) |
| 274 | + lbl = z == 0 ? "z = 0 (rest frame)" : "z = $z" |
| 275 | + plot!(p, spectral_axis(spec_z), ustrip.(flux_axis(spec_z)), label = lbl, lw = 2, color = c) |
| 276 | +end |
| 277 | +p |
| 278 | +``` |
| 279 | + |
| 280 | +--- |
| 281 | + |
| 282 | +## Summary |
| 283 | + |
| 284 | +1. **Loading real spectral data** from SDSS FITS files with `FITSIO.jl` |
| 285 | +2. **Physical units** with `Unitful.jl` + `UnitfulAstro.jl` |
| 286 | +3. **Dust extinction** with CCM89, OD94, CAL00 laws via `DustExtinction.jl` |
| 287 | +4. **Spectral axis inspection** using `spectral_axis` and `flux_axis` |
| 288 | +5. **Uncertainty propagation** via `Measurements.jl` |
| 289 | +6. **Synthetic blackbody spectra** from `Spectra.jl` |
| 290 | +7. **Spectral arithmetic** — sky subtraction, scaling |
| 291 | +8. **Cosmological redshift** and luminosity distances via `Cosmology.jl` |
| 292 | + |
| 293 | +### Where to go next |
| 294 | + |
| 295 | +| Package | What it adds | |
| 296 | +|---|---| |
| 297 | +| [`SpectralFitting.jl`](https://juliaastro.org/SpectralFitting/stable/) | X-ray spectral fitting with OGIP PHA/RMF/ARF | |
| 298 | +| [`Korg.jl`](https://ajwheeler.github.io/Korg.jl/stable/) | Synthetic stellar spectra from model atmospheres | |
| 299 | +| [`AstroImages.jl`](https://juliaastro.org/AstroImages/stable/) | IFU data-cubes — spectra at every spaxel | |
0 commit comments