Basics#

import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from astropy.utils.data import clear_download_cache
clear_download_cache()

from spextra import Spextrum, spextra_database, SpecLibrary, FilterSystem
/home/docs/checkouts/readthedocs.org/user_builds/spextrahb/envs/stable/lib/python3.11/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Examining the database#

sdb = spextra_database

Information of the database#

print(sdb)
Spextra Database:
  Remote URL: https://scopesim.univie.ac.at/spextra/database/
  Local path: /home/docs/.spextra_cache
Database contents:
├─libraries:
│ ├─ref: Library of reference stars
│ ├─kc96: Kinney-Calzetti Atlas
│ ├─pickles: Pickles Stellar Library
│ ├─dobos: SDSS galaxy composite spectra
│ ├─irtf: IRTF spectral library
│ ├─agn: AGN templates
│ ├─nebulae: Emission line nebulae
│ ├─brown: Galaxy SEDs from the UV to the Mid-IR
│ ├─kurucz: Subset of Kurucz 1993 Models
│ ├─sne: Supernova Legacy Survey
│ ├─moehler: flux/telluric standards with X-Shooter
│ ├─madden: High-Resolution Spectra of Habitable Zone Planets
│ ├─bosz/hr: BOSZ stellar atmosphere Grid - High Resolution
│ ├─bosz/mr: BOSZ stellar atmosphere Grid - Medium Resolution
│ ├─bosz/lr: BOSZ stellar atmosphere Grid - Low Resolution
│ ├─sky: Paranal sky background spectra
│ ├─shapley: Rest-Frame Ultraviolet Spectra of z ∼ 3 Lyman Break Galaxies
│ ├─etc/kurucz: ESO ETC subset of the Kurucz 1993 models
│ ├─etc/marcs/p: ESO ETC subset of the MARCS Stellar Models with Plane Parallel Geometry
│ ├─etc/marcs/s: ESO ETC subset of the MARCS Stellar Models with Spherical Geometry
│ └─etc/misc: Other templates, nubulae and qso
├─extinction_curves:
│ ├─gordon: LMC and SMC extinction laws
│ ├─cardelli: MW extinction laws
│ └─calzetti: extragalactic attenuation curves
└─filter_systems:
  ├─elt/micado: MICADO filters
  ├─elt/metis: METIS filters
  └─etc: ESO ETC standard filters
# Which templates are available?

print(SpecLibrary("sne"))
Spectral Library 'sne': Supernova spectral library
  spectral coverage: UV, VIS, NIR
  wave_unit: Angstrom
  flux_unit: FLAM
  Templates: sn1a, sn1b, sn1c, sn2l, sn2p, sn2n, hyper, pop3_3d, pop3_15d

Extinction curves and Filters#

print(sdb["extinction_curves"])
print(sdb["filter_systems"])
{'gordon': 'LMC and SMC extinction laws', 'cardelli': 'MW extinction laws', 'calzetti': 'extragalactic attenuation curves'}
{'elt/micado': 'MICADO filters', 'elt/metis': 'METIS filters', 'etc': 'ESO ETC standard filters'}
print(FilterSystem("micado"))
Filter system 'micado': <untitled>
  spectral coverage: 
  wave_unit: Angstrom
  filters: 

Retrieving the spectra#

sp1 = Spextrum("kc96/s0")
sp1.plot()
Downloading file 'libraries/kc96/s0.fits' from 'https://scopesim.univie.ac.at/spextra/database/libraries/kc96/s0.fits' to '/home/docs/.spextra_cache'.
  0%|                                              | 0.00/20.2k [00:00<?, ?B/s]
 36%|█████████████▏                       | 7.17k/20.2k [00:00<00:00, 63.0kB/s]
  0%|                                              | 0.00/20.2k [00:00<?, ?B/s]
100%|█████████████████████████████████████| 20.2k/20.2k [00:00<00:00, 19.9MB/s]

../_images/581f65f184fe94c989fae4931f2a3827690d35c5e382ad90fb12934441fc76ae.png
# another spectrum

sp2 = Spextrum("agn/qso")
sp2.plot()
Downloading file 'libraries/agn/qso.fits' from 'https://scopesim.univie.ac.at/spextra/database/libraries/agn/qso.fits' to '/home/docs/.spextra_cache'.
  0%|                                              | 0.00/23.0k [00:00<?, ?B/s]
 31%|███████████▌                         | 7.17k/23.0k [00:00<00:00, 61.0kB/s]
  0%|                                              | 0.00/23.0k [00:00<?, ?B/s]
100%|█████████████████████████████████████| 23.0k/23.0k [00:00<00:00, 21.2MB/s]

../_images/e560c24cc74e502427e9157b70ad2df7a51ec238ee8a68e5e28ff1d1453c2eb9.png

Aritmetics#

simple arithmetics are possible

sp = sp1 + 3*sp2

sp.plot()
../_images/80f8a4e730887c5a3508cbf8d57276bb00f550cc473995eb024e84cca6882eb2.png

Adding emission lines#

It is possible to add emission lines, either individually or as a list. Parameters are center, flux and fwhm

sp3 = sp1.add_emi_lines(center=4000,flux=4e-13, fwhm=5*u.AA)
sp3.plot(left=3500, right=4500)
../_images/639f64ab722d481cf0aa7ea5ed5ea87064188eb44b3be6a6478615164f721d33.png

Redshifting spectra#

fig = plt.figure(figsize=(10,7))
sp4 = sp2.redshift(z=1)

wave = sp2.waveset
flux = sp2(wave, flux_unit="FLAM")

plt.plot(wave, flux)

plt.plot(sp4.waveset, 
         sp4(sp4.waveset, flux_unit="FLAM"))

plt.legend(['z=0', 'z=1'], loc='upper right')
<matplotlib.legend.Legend at 0x7adf0db74a50>
../_images/ca21185294139904f69523ae411b7958c1a49ccf986716a8faf89407704d3a5a.png

Or using velocity#

fig = plt.figure(figsize=(10,6))

sp1 = Spextrum("nebulae/orion")

vel = -1000 * u.km / u.s
sp2 = sp1.redshift(vel=vel)

plt.plot(sp1.waveset, 
         sp1(sp1.waveset, flux_unit="FLAM"))
plt.plot(sp2.waveset, 
         sp2(sp2.waveset, flux_unit="FLAM"))
plt.legend(['vel=0', 'vel=-1000 km/s'], loc='upper right')
plt.xlim(3000,5000)
Downloading file 'libraries/nebulae/orion.fits' from 'https://scopesim.univie.ac.at/spextra/database/libraries/nebulae/orion.fits' to '/home/docs/.spextra_cache'.
  0%|                                               | 0.00/170k [00:00<?, ?B/s]
  8%|███▎                                   | 14.3k/170k [00:00<00:01, 122kB/s]
 16%|██████                                 | 26.6k/170k [00:00<00:01, 111kB/s]
 27%|██████████▌                            | 46.1k/170k [00:00<00:00, 135kB/s]
 46%|█████████████████▊                     | 77.8k/170k [00:00<00:00, 186kB/s]
 74%|█████████████████████████████▍          | 125k/170k [00:00<00:00, 261kB/s]
  0%|                                               | 0.00/170k [00:00<?, ?B/s]
100%|████████████████████████████████████████| 170k/170k [00:00<00:00, 187MB/s]

(3000.0, 5000.0)
../_images/5b178d77ae644530f800f81a03c6be05cf5e56f7309d82dd60fd971b9bf7471e.png

Flat spectrum in any photometric system#

(aka reference spectrum if mag=0)

sp_vega = Spextrum.flat_spectrum(10*u.mag)
sp_ab = Spextrum.flat_spectrum(10*u.ABmag)
sp_st = Spextrum.flat_spectrum(10*u.STmag)


fig = plt.figure(figsize=(10,7))
wave = sp_vega.waveset
plt.plot(wave, sp_vega(wave, flux_unit="FLAM"), label="Vega")
plt.plot(wave, sp_ab(wave, flux_unit="FLAM"), label="AB")
plt.plot(wave, sp_st(wave, flux_unit="FLAM"), label="ST")

plt.xlim(3000,1e4)
plt.ylim(0,0.2e-11)
plt.xlabel("wavelength")
plt.ylabel("flux (FLAM)")
plt.legend()
Downloading file 'libraries/ref/vega.fits' from 'https://scopesim.univie.ac.at/spextra/database/libraries/ref/vega.fits' to '/home/docs/.spextra_cache'.
  0%|                                               | 0.00/276k [00:00<?, ?B/s]
  5%|██                                     | 14.3k/276k [00:00<00:02, 124kB/s]
 17%|██████▌                                | 46.1k/276k [00:00<00:01, 210kB/s]
 29%|███████████▎                           | 79.9k/276k [00:00<00:00, 245kB/s]
 51%|████████████████████▎                   | 140k/276k [00:00<00:00, 350kB/s]
 69%|███████████████████████████▍            | 189k/276k [00:00<00:00, 374kB/s]
 92%|████████████████████████████████████▉   | 255k/276k [00:00<00:00, 435kB/s]
  0%|                                               | 0.00/276k [00:00<?, ?B/s]
100%|████████████████████████████████████████| 276k/276k [00:00<00:00, 241MB/s]

<matplotlib.legend.Legend at 0x7adf0d91c710>
../_images/8357cc5a888c9851776cbac64e381a27d6aca7f72103851b6b78afaf42fb694a.png

Scaling to a magnitude#

sp1 = Spextrum("kc96/s0").scale_to_magnitude(amplitude=13 * u.ABmag, filter_curve="g")
sp2 = sp1.scale_to_magnitude(amplitude=15 * u.ABmag, filter_curve="g")

sig = plt.figure(figsize=(10,7))
plt.plot(sp1.waveset, 
         sp1(sp1.waveset, flux_unit="FLAM"))
plt.plot(sp2.waveset, 
         sp2(sp2.waveset, flux_unit="FLAM"))
plt.legend(['mag=13', 'mag=15'], loc='upper right')
plt.xlim(4000,7000)
(4000.0, 7000.0)
../_images/917c67de4cfc3b937031065730ff902f7b4792cff9879be5adc1cc65eaacce0a.png

Obtaining magnitudes from spectra#

print("Magnitude spectra 1:", sp1.get_magnitude(filter_curve="g"), 
      sp1.get_magnitude(filter_curve="g", system_name="Vega"), "Vega")
print("Magnitude spectra 2:", sp2.get_magnitude(filter_curve="g"), 
      sp2.get_magnitude(filter_curve="g", system_name="Vega"), "Vega")
Magnitude spectra 1: 13.0 mag(AB) 13.095325163298817 mag Vega
Magnitude spectra 2: 15.0 mag(AB) 15.095325163298817 mag Vega

Rebin spectra#

a new wavelength array must be passed

sp1 = Spextrum("agn/qso")
new_waves = np.linspace(np.min(sp1.waveset),
                        np.max(sp1.waveset),
                        100)
sp2 = sp1.rebin_spectra(new_waves=new_waves)

sig = plt.figure(figsize=(10,7))
plt.plot(sp1.waveset, 
         sp1(sp1.waveset, flux_unit="FLAM"))
plt.plot(sp2.waveset, 
         sp2(sp2.waveset, flux_unit="FLAM"))
plt.xlim(1000,4000)
(1000.0, 4000.0)
../_images/4654f142fe1a82f1a4ea7d2495945bbfeef9bb312c65fe79665b1e380869385f.png

Smooth the spectral with a velocity kernel#

sp1 = Spextrum("nebulae/pn")

sigma = 500*(u.km / u.s)
sp2 = sp1.smooth(sigma=sigma)

fig = plt.figure(figsize=(10,7))
plt.plot(sp1.waveset, 
         sp1(sp1.waveset, flux_unit="FLAM"), label="original")
plt.plot(sp2.waveset, 
         sp2(sp2.waveset, flux_unit="FLAM"), label="broadened with 500 km/s")

plt.xlim(4800,5200)
plt.legend()
Downloading file 'libraries/nebulae/pn.fits' from 'https://scopesim.univie.ac.at/spextra/database/libraries/nebulae/pn.fits' to '/home/docs/.spextra_cache'.
  0%|                                               | 0.00/170k [00:00<?, ?B/s]
  4%|█▌                                    | 7.17k/170k [00:00<00:02, 62.5kB/s]
  8%|███▏                                  | 14.3k/170k [00:00<00:02, 61.8kB/s]
 23%|████████▉                              | 38.9k/170k [00:00<00:01, 130kB/s]
 37%|██████████████▌                        | 63.5k/170k [00:00<00:00, 161kB/s]
 55%|█████████████████████▌                 | 94.2k/170k [00:00<00:00, 198kB/s]
 74%|█████████████████████████████▍          | 125k/170k [00:00<00:00, 220kB/s]
  0%|                                               | 0.00/170k [00:00<?, ?B/s]
100%|████████████████████████████████████████| 170k/170k [00:00<00:00, 190MB/s]

<matplotlib.legend.Legend at 0x7adf0db34110>
../_images/8949bd34178219703f57820a987f4dd1e927886a5d52790ffc2b83638de30b91.png

Blackbody spectrum and extinction curves#

sp1 = Spextrum.black_body_spectrum(temperature=5500, 
                                   amplitude=10 * u.ABmag, 
                                   filter_curve="r")
sp2 = sp1.redden("gordon/smc_bar", Ebv=0.15)


fig = plt.figure(figsize=(10,7))
plt.plot(sp1.waveset, 
         sp1(sp1.waveset, flux_unit="FLAM"), label="original")
plt.plot(sp2.waveset, 
         sp2(sp2.waveset, flux_unit="FLAM"), label="attenuated")

plt.xlim(1800,15200)
plt.legend()
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/spextrahb/envs/stable/lib/python3.11/site-packages/astropy/units/decorators.py:69, in _validate_arg_value(param_name, func_name, arg, targets, equivalencies, strict_dimensionless)
     68 try:
---> 69     if arg.unit.is_equivalent(allowed_unit, equivalencies=equivalencies):
     70         break

AttributeError: 'int' object has no attribute 'unit'

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
Cell In[18], line 1
----> 1 sp1 = Spextrum.black_body_spectrum(temperature=5500, 
      2                                    amplitude=10 * u.ABmag, 
      3                                    filter_curve="r")
      4 sp2 = sp1.redden("gordon/smc_bar", Ebv=0.15)
      7 fig = plt.figure(figsize=(10,7))

File ~/checkouts/readthedocs.org/user_builds/spextrahb/envs/stable/lib/python3.11/site-packages/astropy/units/decorators.py:298, in QuantityInput.__call__.<locals>.wrapper(*func_args, **func_kwargs)
    290         valid_targets = [
    291             t
    292             for t in valid_targets
    293             if isinstance(t, (str, UnitBase, PhysicalType))
    294         ]
    296     # Now we loop over the allowed units/physical types and validate
    297     #   the value of the argument:
--> 298     _validate_arg_value(
    299         param.name,
    300         wrapped_function.__name__,
    301         arg,
    302         valid_targets,
    303         self.equivalencies,
    304         self.strict_dimensionless,
    305     )
    307 if self.equivalencies:
    308     equiv_context = add_enabled_equivalencies(self.equivalencies)

File ~/checkouts/readthedocs.org/user_builds/spextrahb/envs/stable/lib/python3.11/site-packages/astropy/units/decorators.py:78, in _validate_arg_value(param_name, func_name, arg, targets, equivalencies, strict_dimensionless)
     75         else:
     76             error_msg = "no 'unit' attribute"
---> 78         raise TypeError(
     79             f"Argument '{param_name}' to function '{func_name}'"
     80             f" has {error_msg}. You should pass in an astropy "
     81             "Quantity instead."
     82         )
     84 else:
     85     error_msg = (
     86         f"Argument '{param_name}' to function '{func_name}' must "
     87         "be in units convertible to"
     88     )

TypeError: Argument 'temperature' to function 'black_body_spectrum' has no 'unit' attribute. You should pass in an astropy Quantity instead.

Photons within a filter#

(or between wmin or wmax)

n_photons = sp2.photons_in_range(area=2*u.m**2,
                                 filter_curve="V")

print(n_photons)