Merge pull request #86 from karllark/stis_lsfs

STIS LSFs used in mocking of spectra

docs/model_standards.rst

@@ -57,3 +57,22 @@ observed.  The specific code is `utils/make_obsdata_from_model.py`.
 This code uses the 'merge_obsspec' functions to transform the model SEDs
 to the observed spectral formats.
 The `utils/make_all_tlusty_obsdata.py` runs on the `*.flux.gz` tlusty files.
+
+STIS Mocking
+^^^^^^^^^^^^
+
+The HST STIS observations are simulated by convolving the high spectral
+resolution model spectra to the STIS resolution using STIS line-spread-fuctions (LSFs)
+retrieved from
+`STScI <https://www.stsci.edu/hst/instrumentation/stis/performance/spectral-resolution>`_.
+These line-spread functions are provided at specific wavelengths and are linearly
+interpolated/extrapolated for other wavelengths (see `utils/mock_spectra_data.py`
+for details.)
+
+STIS LSFs (52x2 slit):
+
+.. image:: images/stis_lsfs.png
+
+Example of mocked STIS observations for all four low-resolution grating settings:
+
+.. image:: images/mock_stis_obs.png

measure_extinction/merge_obsspec.py

@@ -137,7 +137,7 @@ def merge_stis_obsspec(obstables, waveregion="UV", output_resolution=1000):
     if waveregion == "UV":
         wave_range = [1000.0, 3400.0] * u.angstrom
     elif waveregion == "Opt":
-        wave_range = [2900.0, 10250.0] * u.angstrom
+        wave_range = [2850.0, 10250.0] * u.angstrom
 
     iwave_range = wave_range.to(u.angstrom).value
     full_wave, full_wave_min, full_wave_max = _wavegrid(output_resolution, iwave_range)
@@ -150,19 +150,17 @@ def merge_stis_obsspec(obstables, waveregion="UV", output_resolution=1000):
         # may want to add in the SYS-ERROR, but need to be careful
         # to propagate it correctly, SYS-ERROR will not reduce with
         # multiple spectra or measurements in a wavelength bin
-        cuncs = ctable["STAT-ERROR"].data
+        cuncs = ctable["STAT-ERROR"]
         cwaves = ctable["WAVELENGTH"].data
-        cfluxes = ctable["FLUX"].data
+        cfluxes = ctable["FLUX"]
         cnpts = ctable["NPTS"].data
         for k in range(n_waves):
-            (indxs,) = np.where(
-                (cwaves >= full_wave_min[k]) & (cwaves < full_wave_max[k]) & (cnpts > 0)
-            )
-            if len(indxs) > 0:
-                weights = 1.0 / np.square(cuncs[indxs])
-                full_flux[k] += np.sum(weights * cfluxes[indxs])
+            gvals = (cwaves >= full_wave_min[k]) & (cwaves < full_wave_max[k]) & (cnpts > 0)
+            if np.sum(gvals) > 0:
+                weights = 1.0 / np.square(cuncs[gvals].value)
+                full_flux[k] += np.sum(weights * cfluxes[gvals].value)
                 full_unc[k] += np.sum(weights)
-                full_npts[k] += len(indxs)
+                full_npts[k] += np.sum(gvals)
 
     # divide by the net weights
     (indxs,) = np.where(full_npts > 0)

measure_extinction/utils/STIS_LSF/plot_lsfs.py

@@ -0,0 +1,93 @@
+import argparse
+
+import matplotlib.pyplot as plt
+from astropy.table import Table
+
+
+if __name__ == "__main__":
+
+    # commandline parser
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--png", help="save figure as a png file", action="store_true")
+    parser.add_argument("--pdf", help="save figure as a pdf file", action="store_true")
+    args = parser.parse_args()
+
+    fig, ax = plt.subplots(ncols=4, figsize=(16, 8))
+
+    # setup the plots
+    fontsize = 12
+    font = {"size": fontsize}
+    plt.rc("font", **font)
+    plt.rc("lines", linewidth=2)
+    plt.rc("axes", linewidth=2)
+    plt.rc("xtick.major", width=2)
+    plt.rc("ytick.major", width=2)
+
+    # resolutions in A/pixel
+    g140l_res = 0.5831004076162168
+    g230l_res = 1.548239955089159
+    g430l_res = 2.7463714795193273
+    g750l_res = 4.879600079462369
+
+    tags = ["G140L_1200", "G140L_1500"]
+    for i, ctag in enumerate(tags):
+        a = Table.read(
+            f"data/LSF_{ctag}.txt", format="ascii.commented_header", header_start=-1
+        )
+
+        # slits = ["52x0.1", "52x0.2", "52x0.5", "52x2.0"]
+        slits = ["52x2.0"]
+        for cslit in slits:
+            ax[0].plot(g140l_res * a["Rel_pixel"], a[cslit], label=ctag)
+        ax[0].legend()
+    ax[0].set_xlabel(r"$\Delta\lambda$ [$\AA$]")
+    ax[0].set_xlim(-20 * g140l_res, 20 * g140l_res)
+
+    tags = ["G230L_1700", "G230L_2400"]
+    for i, ctag in enumerate(tags):
+        a = Table.read(
+            f"data/LSF_{ctag}.txt", format="ascii.commented_header", header_start=-1
+        )
+
+        # slits = ["52x0.1", "52x0.2", "52x0.5", "52x2.0"]
+        slits = ["52x2.0"]
+        for cslit in slits:
+            ax[1].plot(g230l_res * a["Rel_pixel"], a[cslit], label=ctag)
+        ax[1].legend()
+    ax[1].set_xlabel(r"$\Delta\lambda$ [$\AA$]")
+    ax[1].set_xlim(-20 * g230l_res, 20 * g230l_res)
+
+    tags = ["G430L_3200", "G430L_5500"]
+    for i, ctag in enumerate(tags):
+        a = Table.read(
+            f"data/LSF_{ctag}.txt", format="ascii.commented_header", header_start=-1
+        )
+
+        # slits = ["52x0.1", "52x0.2", "52x0.5", "52x2.0"]
+        slits = ["52x2.0"]
+        for cslit in slits:
+            ax[2].plot(g430l_res * a["Rel_pixel"], a[cslit], label=ctag)
+        ax[2].legend()
+    ax[2].set_xlabel(r"$\Delta\lambda$ [$\AA$]")
+    ax[2].set_xlim(-20 * g430l_res, 20 * g430l_res)
+
+    tags = ["G750L_7000"]
+    for i, ctag in enumerate(tags):
+        a = Table.read(
+            f"data/LSF_{ctag}.txt", format="ascii.commented_header", header_start=-1
+        )
+
+        # slits = ["52x0.1", "52x0.2", "52x0.5", "52x2.0"]
+        slits = ["52x2.0"]
+        for cslit in slits:
+            ax[3].plot(g750l_res * a["Rel_pixel"], a[cslit], label=ctag)
+        ax[3].legend()
+    ax[3].set_xlabel(r"$\Delta\lambda$ [$\AA$]")
+    ax[3].set_xlim(-20 * g750l_res, 20 * g750l_res)
+
+    fig.tight_layout()
+
+    if args.png:
+        fig.savefig("stis_lsfs.png")
+    else:
+        plt.show()

measure_extinction/utils/make_obsdata_from_model.py

@@ -1,22 +1,21 @@
 #!/usr/bin/env python
 
-from __future__ import absolute_import, division, print_function, unicode_literals
-
 import pkg_resources
 
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.interpolate import interp1d
 
 from astropy.io import ascii
-from astropy.table import Table, Column
+from astropy.table import QTable, Column
 import astropy.units as u
 from astropy.convolution import Gaussian1DKernel, convolve
 from synphot import SpectralElement
 import stsynphot as STS
 
 from measure_extinction.stardata import BandData
 from measure_extinction.merge_obsspec import merge_stis_obsspec, merge_irs_obsspec
+from measure_extinction.utils.mock_spectra_data import mock_stis_data
 
 __all__ = ["make_obsdata_from_model"]
 
@@ -262,14 +261,14 @@ def make_obsdata_from_model(
     #   means that SFlux is read in as a string
     # solution is to remove the rows with the problem and replace
     #   the fortran 'D' with an 'E' and then convert to floats
-    if mspec["SFlux"].dtype != np.float:
+    if mspec["SFlux"].dtype != float:
         indxs = [k for k in range(len(mspec)) if "D" not in mspec["SFlux"][k]]
         if len(indxs) > 0:
             indxs = [k for k in range(len(mspec)) if "D" in mspec["SFlux"][k]]
             mspec = mspec[indxs]
             new_strs = [cval.replace("D", "E") for cval in mspec["SFlux"].data]
             mspec["SFlux"] = new_strs
-            mspec["SFlux"] = mspec["SFlux"].astype(np.float)
+            mspec["SFlux"] = mspec["SFlux"].astype(float)
 
     # set the units
     mspec["Freq"].unit = u.Hz
@@ -290,11 +289,11 @@ def make_obsdata_from_model(
     )
 
     # save the full spectrum to a binary FITS table
-    otable = Table()
+    otable = QTable()
     otable["WAVELENGTH"] = Column(wave_r5000, unit=u.angstrom)
     otable["FLUX"] = Column(flux_r5000, unit=u.erg / (u.s * u.cm * u.cm * u.angstrom))
     otable["SIGMA"] = Column(
-        flux_r5000 * 0.0, unit=u.erg / (u.s * u.cm * u.cm * u.angstrom)
+        flux_r5000 * 0.01, unit=u.erg / (u.s * u.cm * u.cm * u.angstrom)
     )
     otable["NPTS"] = Column(npts_r5000)
     otable.write(
@@ -305,28 +304,16 @@ def make_obsdata_from_model(
     specinfo = {}
 
     # create the ultraviolet HST/STIS mock observation
-    # first create the spectrum convolved to the STIS low resolution
-    # Resolution approximately 1000
-    stis_fwhm_pix = 5000.0 / 1000.0
-    g = Gaussian1DKernel(stddev=stis_fwhm_pix / 2.355)
-
-    # Convolve data
-    nflux = convolve(otable["FLUX"].data, g)
+    stis_table = mock_stis_data(otable)
 
-    stis_table = Table()
-    stis_table["WAVELENGTH"] = otable["WAVELENGTH"]
-    stis_table["FLUX"] = nflux
-    stis_table["NPTS"] = otable["NPTS"]
-    stis_table["STAT-ERROR"] = Column(np.full((len(stis_table)), 1.0))
-    stis_table["SYS-ERROR"] = otable["SIGMA"]
     # UV STIS obs
-    rb_stis_uv = merge_stis_obsspec([stis_table], waveregion="UV")
+    rb_stis_uv = merge_stis_obsspec(stis_table[0:2], waveregion="UV")
     rb_stis_uv["SIGMA"] = rb_stis_uv["FLUX"] * 0.0
     stis_uv_file = "%s_stis_uv.fits" % (output_filebase)
     rb_stis_uv.write("%s/Models/%s" % (output_path, stis_uv_file), overwrite=True)
     specinfo["STIS"] = stis_uv_file
     # Optical STIS obs
-    rb_stis_opt = merge_stis_obsspec([stis_table], waveregion="Opt")
+    rb_stis_opt = merge_stis_obsspec(stis_table[2:4], waveregion="Opt")
     rb_stis_opt["SIGMA"] = rb_stis_opt["FLUX"] * 0.0
     stis_opt_file = "%s_stis_opt.fits" % (output_filebase)
     rb_stis_opt.write("%s/Models/%s" % (output_path, stis_opt_file), overwrite=True)
@@ -340,7 +327,7 @@ def make_obsdata_from_model(
     # Convolve data
     nflux = convolve(otable["FLUX"].data, g)
 
-    lrs_table = Table()
+    lrs_table = QTable()
     lrs_table["WAVELENGTH"] = otable["WAVELENGTH"]
     lrs_table["FLUX"] = nflux
     lrs_table["NPTS"] = otable["NPTS"]