# Tutorials¶

This page contains tutorials of specific stsynphot functionality not explicitly covered in other sections.

## Count Rates for Multiple Apertures¶

In this tutorial, you will learn how to calculate count rates for observations of the same source and bandpass, but with different apertures. Note that this feature is only available for observing modes that allow encircled energy (EE) radius specification (see Appendix B: OBSMODE Keywords).

Create two observations of Vega (renormalized to 20 STMAG in Johnson V) with ACS/WFC1 F555W bandpass, with 0.3 and 1.0 arcsec EE radii, respectively:

```>>> import stsynphot as stsyn
>>> from astropy import units as u
>>> from synphot import Observation
>>> sp = stsyn.Vega.normalize(20 * u.STmag, stsyn.band('johnson,v'))
>>> obs03 = Observation(sp, stsyn.band('acs,wfc1,f555w,aper#0.3'))
>>> obs10 = Observation(sp, stsyn.band('acs,wfc1,f555w,aper#1.0'))
```

Calculate the count rates for both and display the results:

```>>> c03 = obs03.countrate(stsyn.conf.area)
>>> c10 = obs10.countrate(stsyn.conf.area)
>>> print(f'Count rate for 0.3" is {c03:.3f}\n'
...       f'Count rate for 1.0" is {c10:.3f}')
Count rate for 0.3" is 174.801 ct / s
Count rate for 1.0" is 186.521 ct / s
```

## Bandpass STMAG Zeropoint¶

HST bandpasses store their Bandpass Unit Response values under the `PHOTFLAM` keyword in image headers. This keyword is then used to compute STMAG zeropoint for the respective bandpass (e.g., ACS and WFC3).

In this tutorial, you will learn how to calculate the STMAG zeropoint for the ACS/WFC1 F555W bandpass, which happens to be time-dependent:

```>>> import numpy as np
>>> import stsynphot as stsyn
>>> from astropy.time import Time
>>> obsdate = Time('2017-05-30').mjd
>>> obsmode = f'acs,wfc1,f555w,mjd#{obsdate}'
>>> bp = stsyn.band(obsmode)
>>> photflam = bp.unit_response(stsyn.conf.area)
>>> photflam
<Quantity 1.9647813651514673e-19 FLAM>
```
```>>> st_zpt = -2.5 * np.log10(photflam.value) - 21.1
>>> print(f'STmag zeropoint for {obsmode} is {st_zpt:.5f}')
STmag zeropoint for acs,wfc1,f555w,mjd#57903.0 is 25.66671
```

## Bandpass ABMAG Zeropoint¶

For ABMAG zeropoint, it extends from Bandpass STMAG Zeropoint by also using `PHOTPLAM` keyword in image headers.

In this tutorial, you will learn how to calculate the ABMAG zeropoint for the ACS/WFC1 F555W bandpass, which happens to be time-dependent:

```>>> import numpy as np
>>> import stsynphot as stsyn
>>> from astropy.time import Time
>>> obsdate = Time('2017-05-30').mjd
>>> obsmode = f'acs,wfc1,f555w,mjd#{obsdate}'
>>> bp = stsyn.band(obsmode)
>>> photflam = bp.unit_response(stsyn.conf.area)
>>> photplam = bp.pivot()
>>> photplam
<Quantity 5360.938362432486 Angstrom>
```
```>>> ab_zpt = (-2.5 * np.log10(photflam.value) - 21.1 -
...           5 * np.log10(photplam.value) + 18.6921)
>>> print(f'ABmag zeropoint for {obsmode} is {ab_zpt:.5f}')
ABmag zeropoint for acs,wfc1,f555w,mjd#57903.0 is 25.71261
```

## Bandpass VEGAMAG Zeropoint¶

In addition to Bandpass STMAG Zeropoint and Bandpass ABMAG Zeropoint, HST bandpasses also provide zeropoints in `VEGAMAG`, which is a magnitude system where Vega has magnitude 0 at all wavelengths. Note that this zeropoint strongly depends on the actual Vega spectrum used; Therefore, VEGAMAG zeropoint values for the same filter might vary in literature as the authors use their favorite Vega spectra.

In this tutorial, you will learn how to calculate the VEGAMAG zeropoint for the ACS/WFC1 F555W bandpass, which happens to be time-dependent:

```>>> import numpy as np
>>> import stsynphot as stsyn
>>> from astropy.time import Time
>>> from synphot import Observation
>>> obsdate = Time('2017-05-30').mjd
>>> bp = stsyn.band(f'acs,wfc1,f555w,mjd#{obsdate}')
>>> obs = Observation(stsyn.Vega, bp, binset=bp.binset)
>>> vega_zpt = -obs.effstim(flux_unit='obmag', area=stsyn.conf.area)
>>> print(f'VEGAMAG zeropoint for {bp.obsmode} is {vega_zpt:.5f}')
VEGAMAG zeropoint for acs,wfc1,f555w,mjd#57903.0 is 25.71235 OBMAG
```

## Sun’s Abs. Mag. in HST Filters¶

In this tutorial, you will learn how to calculate the absolute magnitude of the Sun for three different HST filters. Sun’s spectrum can be obtained from HST Calibration Spectra but needs to be normalized to literature value (e.g., http://www.astronomynotes.com/starprop/s4.htm):

```>>> import stsynphot as stsyn
>>> from synphot import units, SourceSpectrum, SpectralElement, Observation
>>> v_band = SpectralElement.from_filter('johnson_v')
>>> sun_file = 'https://archive.stsci.edu/hlsps/reference-atlases/cdbs/calspec/sun_reference_stis_002.fits'
>>> sun_raw = SourceSpectrum.from_file(sun_file)
>>> sun = sun_raw.normalize(4.83 * units.VEGAMAG, v_band, vegaspec=stsyn.Vega)
>>> for obsmode in ['acs,wfc1,f555w', 'wfc3,uvis2,f336w', 'wfc3,ir,f160w']:
...     bp = stsyn.band(obsmode)
...     obs = Observation(sun, bp, binset=bp.binset)
...     m = obs.effstim('vegamag', vegaspec=stsyn.Vega)
...     print(f"Sun's abs mag in {bp.obsmode} is {m:.4f}")
Sun's abs mag in acs,wfc1,f555w is 4.8395 VEGAMAG
Sun's abs mag in wfc3,uvis2,f336w is 5.4864 VEGAMAG
Sun's abs mag in wfc3,ir,f160w is 3.4127 VEGAMAG
```

## Custom Wavelength Table¶

In this tutorial, you will learn how to create a custom wavelength array and save it to a FITS table using `astropy.io.fits`. Then, you will read the array back in from file, and use it to define detector binning for an observation.

Suppose we want a wavelength set that ranges from 2000 to 8000 Angstrom, with 1 Angstrom spacing over most of the range, but 0.1 Angstrom spacing around the [O III] forbidden lines at 4959 and 5007 Angstrom.

Create the 3 regions separately, concatenate them, and display the result:

```>>> import numpy as np
>>> lowave = np.arange(2000, 4950)
>>> mdwave = np.arange(4950, 5010, 0.1)  # [O III]
>>> hiwave = np.arange(5010, 8000)
>>> wave = np.concatenate([lowave, mdwave, hiwave])
>>> wave
array([2000., 2001., 2002., ..., 7997., 7998., 7999.])
```

Create an Astropy table from the concatenated array above and save it out as a FITS table:

```>>> from astropy.io import fits
>>> col = fits.Column(
...     name='wavelength', unit='angstrom', format='E', array=wave)
>>> tabhdu = fits.BinTableHDU.from_columns([col])
>>> tabhdu.writeto('mywaveset.fits')
```

Read the custom wavelength set back in from file using Astropy table:

```>>> from astropy.table import QTable
>>> tab = QTable.read('mywaveset.fits')  # Ignore the UnitsWarning
WARNING: UnitsWarning: The unit 'angstrom' has been deprecated...
```
```>>> wave = tab['wavelength']
>>> wave
<Quantity [2000., 2001., 2002., ..., 7997., 7998., 7999.] Angstrom>
```

Create an observation of Vega with ACS/WFC1 F555W bandpass, using the custom wavelength binning above, and then check that the binned wavelength set is indeed the given one:

```>>> import stsynphot as stsyn
>>> from synphot import Observation
>>> obs = Observation(stsyn.Vega, stsyn.band('acs,wfc1,f555w'), binset=wave)
>>> obs.binset
<Quantity [ 2000., 2001., 2002.,...,  7997., 7998., 7999.] Angstrom>
```