.. _synphot-tutorials: Tutorials ========= This page contains tutorials of specific **synphot** functionality not explicitly covered in other sections. .. _tutorial_em_line: Emission Line ------------- This tutorial is adapted from `Exposure Time Calculator User's Guide on a similar topic `_. In this tutorial, you will learn how to manipulate and superimpose an emission line to a continuum spectrum. Create a continuum spectrum of a 5500 K blackbody with z=0.6:: >>> from synphot import SourceSpectrum >>> from synphot.models import BlackBodyNorm1D >>> bb = SourceSpectrum(BlackBodyNorm1D, temperature=5500, z=0.6) Create a Gaussian emission line with 8E-14 FLAM total flux, FWHM of 100 Angstrom, and centered at 7000 Angstrom:: >>> from astropy import units as u >>> from synphot import units >>> from synphot.models import GaussianFlux1D >>> em = SourceSpectrum( ... GaussianFlux1D, total_flux=8e-14*u.erg/(u.cm**2 * u.s), ... fwhm=100, mean=7000) Add emission line to continuum spectrum:: >>> sp = bb + em Apply extinction curve for LMC (average) with E(B-V)=1.3 to the composite spectrum:: >>> from synphot import ReddeningLaw >>> ext = ReddeningLaw.from_extinction_model( ... 'lmcavg').extinction_curve(1.3) # doctest: +REMOTE_DATA >>> my_spec = sp * ext # doctest: +REMOTE_DATA Plot the result:: >>> my_spec.plot(right=45000) # doctest: +SKIP .. image:: images/tutorial_em_line.png :width: 600px :alt: Emission line tutorial .. _tutorial_continuum_norm: Continuum-Normalized Spectrum ----------------------------- In this tutorial, you will learn how to create a composite spectrum with a noisy blackbody continuum, an emission line, and an absorption line. Then, you will divide it by a smooth continuum and plot the resultant continuum-normalized spectrum. .. plot:: :include-source: import matplotlib.pyplot as plt import numpy as np from astropy import units as u from synphot import SourceSpectrum, BaseUnitlessSpectrum from synphot.models import BlackBodyNorm1D, Empirical1D, GaussianFlux1D np.random.seed(1234) # For reproducibility # Create the smooth continuum that is a 5000 K blackbody. bb = SourceSpectrum(BlackBodyNorm1D, temperature=5000) # Then, add random noise to it. Since synphot spectrum object cannot # be multiplied with scalar array, this has to be done indirectly by # applying the noise as a unitless spectrum. wave = np.arange(100, 30001, 10) nse = 1 + np.random.normal(size=wave.size, scale=0.02) sp_nse = BaseUnitlessSpectrum(Empirical1D, points=wave, lookup_table=nse) bb_noisy = bb * sp_nse # Apply emission and absorption lines to the noisy continuum. tf_unit = u.erg / (u.cm**2 * u.s) g_em = SourceSpectrum( GaussianFlux1D, total_flux=2.65e-14*tf_unit, mean=15000, fwhm=500) g_ab = SourceSpectrum( GaussianFlux1D, total_flux=6.62e-14*tf_unit, mean=4500, fwhm=100) sp = bb_noisy + g_em - g_ab # Divide the noisy spectrum with lines with the original smooth continuum # to obtain the continuum-normalized spectrum. ratio = sp / bb with np.errstate(invalid='ignore'): ratio.plot(left=2500, right=17000, title='Continuum-normalized spectrum') plt.axhline(1, ls='--', color='k') .. _tutorial_dust_extinction: Using dust-extinction model --------------------------- In this tutorial, you will learn how to apply an extinction curve using a model from the :ref:`dust-extinction ` package: .. doctest-requires:: dust-extinction >>> import numpy as np >>> from astropy import units as u >>> from dust_extinction.parameter_averages import CCM89 >>> from synphot import SourceSpectrum, ReddeningLaw >>> from synphot.models import BlackBodyNorm1D >>> ccm89_model = CCM89(Rv=3.1) >>> wav = np.arange(0.1, 3, 0.001) * u.micron >>> ebv = 0.1 # E(B-V) >>> redlaw = ReddeningLaw(ccm89_model) >>> extcurve = redlaw.extinction_curve(ebv, wavelengths=wav) >>> bb = SourceSpectrum(BlackBodyNorm1D, temperature=5000 * u.K) >>> bb.integrate() # doctest: +FLOAT_CMP >>> sp = bb * extcurve >>> sp.integrate() # doctest: +FLOAT_CMP .. _tutorial_fit_ew: Fitting, Equivalent Width ------------------------- In this tutorial, you will learn how to fit a Gaussian model to some real data and calculate its equivalent width. This is not handled by **synphot** but it is included here for those who are interested to see how fitting in IRAF SYNPHOT is done in Python. See :ref:`astropy:astropy-modeling` for more information about fitting a model. Read in the real data. If your own data has a different format, you need to adjust the example accordingly:: >>> from astropy.io import fits >>> with fits.open('/path/to/combined_13330_G130M_v40_bin4.fits') as pf: # doctest: +SKIP ... dat = pf[1].data # doctest: +SKIP ... wave = dat.field('WAVELENGTH').flatten() # Angstrom # doctest: +SKIP ... flux = dat.field('FLUX').flatten() # FLAM # doctest: +SKIP For a good fit, only use data around the feature of interest. In this example, the feature is between 1202 and 1211 Angstrom:: >>> mask = (wave >= 1202) & (wave <= 1211) # doctest: +SKIP >>> x = wave[mask] # doctest: +SKIP >>> y = flux[mask] # doctest: +SKIP Create a composite model with some initial parameters close to the desired result (usually sufficient to guess from looking at a plot of the data) and fit it using some fitter that is best for the data (sometimes, several iterations are required for a good fit):: >>> from astropy.modeling import models, fitting >>> bg = models.Const1D(amplitude=3.5E-14) >>> gs = models.Gaussian1D(amplitude=3.5E-14, mean=1206, stddev=1) >>> init_model = bg - gs >>> fitter = fitting.LevMarLSQFitter() >>> fit_model = fitter(init_model, x, y) # doctest: +SKIP >>> y_fit = fit_model(x) # doctest: +SKIP >>> print(fit_model) # doctest: +SKIP Model: CompoundModel... Inputs: ('x',) Outputs: ('y',) Model set size: 1 Expression: [0] - [1] Components: [0]: [1]: Parameters: amplitude_0 amplitude_1 mean_1 stddev_1 ----------------- ----------------- ------------- ------------- 3.63064137361e-14 3.62623007738e-14 1206.27454371 0.23713207018 Plot the fitted model on top of input data:: >>> import matplotlib.pyplot as plt # doctest: +SKIP >>> from matplotlib import ticker # doctest: +SKIP >>> fig, ax = plt.subplots() # doctest: +SKIP >>> ax.plot(x, y, 'b', x, y_fit, 'r') # doctest: +SKIP >>> ax.get_xaxis().set_major_formatter( ... ticker.FuncFormatter(ticker.FormatStrFormatter('%.0f'))) # doctest: +SKIP >>> ax.set_xlabel('Wavelength (Angstrom)') # doctest: +SKIP >>> ax.set_ylabel('Flux (FLAM)') # doctest: +SKIP >>> ax.legend(['Data', 'Fit'], loc='lower right') # doctest: +SKIP .. image:: images/tutorial_fit_ab.png :width: 600px :alt: Fitting absorption line in data. Calculate equivalent width using the fitted model:: >>> import math >>> area = (math.sqrt(2 * math.pi) * fit_model.amplitude_1 * ... fit_model.stddev_1) # Area under curve # doctest: +SKIP >>> height = fit_model.amplitude_0 # Continuum level # doctest: +SKIP >>> print('EW = {:.4f} Angstrom'.format(area / height)) # doctest: +SKIP EW = 0.5937 Angstrom .. _tutorial_lyman_alpha: Lyman-Alpha Extinction ---------------------- In this tutorial, you will learn how to apply extinction curve due to Lyman-alpha forest (:ref:`Madau et al. 1995 `) to a source spectrum. For clarity, we will only use a flat source. .. plot:: :include-source: import matplotlib.pyplot as plt from synphot import SourceSpectrum, etau_madau from synphot.models import ConstFlux1D # Create a flat source sp = SourceSpectrum(ConstFlux1D, amplitude=1E-4) # Apply extinction for a given redshift z = 2 wave = range(2400, 4200) # Angstrom extcurve = etau_madau(wave, z) sp_ext = sp * extcurve # Compare the source with and without extinction plt.plot(wave, sp(wave), 'b--', wave, sp_ext(wave), 'r') plt.xlabel('Wavelength (Angstrom)') plt.ylabel('Flux (PHOTLAM)') plt.legend(['Original', 'Extincted'], loc='lower right') The chart below illustrates the Madau 1995 extinction curves for different redshift values. For clarity, they are plotted against rest wavelength, not the redshifted wavelength: .. plot:: :include-source: import matplotlib.pyplot as plt import numpy as np from synphot import etau_madau w_rest = np.arange(800, 1400) lc = ['k', 'navy', 'b', 'deepskyblue', 'mediumseagreen', 'lightgreen', 'y', 'orange', 'r'] for z in range(0, 9): wave = w_rest * (1 + z) extcurve = etau_madau(wave, z) plt.plot(w_rest, extcurve(wave), color=lc[z], label='z={}'.format(z)) plt.ylim(0, 1.1) plt.xlabel('Rest-Frame Wavelength (Angstrom)') plt.ylabel('Lyman-alpha Forest "Throughput"') plt.legend(loc='center right') plt.grid() .. _tutorial_xmm_area: Setting Area for XMM-OM ----------------------- Some telescopes (e.g., XMM-OM by ESA) provide throughput curves with effective area information embedded in them already. In the case of XMM-OM, its throughput curves are in ``effective_area * 100`` meter squared, for which we would divide each throughput curve loaded into ``SpectralElement`` by 100. We also would set the primary area, when it is needed, as follows: >>> from astropy import units as u >>> area = 1 * (u.m * u.m)