Finalize NIRCam spectra support, updates for STIS, etc.

measure_extinction/extdata.py

@@ -27,6 +27,7 @@
     "SpeX_SXD",
     "SpeX_LXD",
     "IRS",
+    "NIRCam_SS",
     "MIRI_IFU",
 ]
 
@@ -263,7 +264,7 @@ class ExtData:
     """
     Extinction for a single line-of-sight
 
-    Atributes:
+    Attributes:
 
     type : string
         extinction curve type (e.g., elx or alax)

measure_extinction/merge_obsspec.py

@@ -285,7 +285,7 @@ def merge_gen_obsspec(obstables, wave_range, output_resolution=100):
     Parameters
     ----------
     obstables : list of astropy Table objects
-        list of tables containing the observed IRS spectra
+        list of tables containing the observed spectra
         usually the result of reading tables
 
     wave_range : 2 element float
@@ -367,7 +367,7 @@ def merge_irs_obsspec(obstables, output_resolution=150):
     """
     wave_range = [5.0, 40.0] * u.micron
     otable = merge_gen_obsspec(
-        obstables, wave_range, output_resolution=output_resolution
+        obstables, wave_range, output_resolution=output_resolution,
     )
     return otable
 
@@ -421,6 +421,7 @@ def merge_miri_ifu_obsspec(obstables, output_resolution=3000):
     """
     wave_range = [4.8, 29.0] * u.micron
     otable = merge_gen_obsspec(
-        obstables, wave_range, output_resolution=output_resolution
+        obstables, wave_range, output_resolution=output_resolution,
     )
+
     return otable

measure_extinction/modeldata.py

@@ -209,7 +209,7 @@ def stellar_sed(self, params, velocity=None):
 
         return sed
 
-    def dust_extinguished_sed_FM90_G23(self, params, sed, velocity=0.0):
+    def dust_extinguished_sed(self, params, sed, velocity=0.0):
         """
         Dust extinguished sed given the extinction parameters
 
@@ -254,7 +254,7 @@ def dust_extinguished_sed_FM90_G23(self, params, sed, velocity=0.0):
                 params[6],
                 optnir_axav_x.value,
                 optnir_axav_y,
-                [0.2, 11.0],
+                [0.033, 11.0],
                 "FM90_G23_measure_extinction",
             )
             ext_sed[cspec] = sed[cspec] * (10 ** (-0.4 * axav * params[0]))
@@ -270,7 +270,7 @@ def dust_extinguished_sed_FM90_G23(self, params, sed, velocity=0.0):
 
         return ext_sed
 
-    def dust_extinguished_sed(self, params, sed, velocity=0.0):
+    def dust_extinguished_sed_FM04(self, params, sed, velocity=0.0):
         """
         Dust extinguished sed given the extinction parameters
 

measure_extinction/stardata.py

@@ -81,7 +81,7 @@
         self.band_waves = OrderedDict()
         self.band_fluxes = OrderedDict()
 
-    def read_bands(self, lines):
+    def read_bands(self, lines, only_bands=None):
         """
         Read the photometric band data from a DAT file
         and upate class variables.
@@ -92,6 +92,8 @@
         ----------
         lines : list of string
             lines from a DAT formatted file
+        only_bands : list
+            Only read in the bands given
 
         Returns
         -------
@@ -110,15 +112,23 @@
                 elif line.find(";") != -1 and line.find("mJy") != -1:
                     colpos = min(line.find(";"), line.find("mJy"))
                 band_name = line[0:eqpos].strip()
-                self.bands[band_name] = (
-                    float(line[eqpos + 1 : pmpos]),
-                    float(line[pmpos + 3 : colpos]),
-                )
-                # units
-                if line.find("mJy") >= 0:
-                    self.band_units[band_name] = "mJy"
+
+                save_band = False
+                if only_bands is None:
+                    save_band = True
                 else:
-                    self.band_units[band_name] = "mag"
+                    if band_name in only_bands:
+                        save_band = True
+                if save_band:
+                    self.bands[band_name] = (
+                        float(line[eqpos + 1 : pmpos]),
+                        float(line[pmpos + 3 : colpos]),
+                    )
+                    # units
+                    if line.find("mJy") >= 0:
+                        self.band_units[band_name] = "mJy"
+                    else:
+                        self.band_units[band_name] = "mag"
 
         self.n_bands = len(self.bands)
 
@@ -815,6 +825,30 @@
         self.fluxes = self.fluxes.value * (u.Jy)
         self.uncs = self.uncs.value * (u.Jy)
 
+    def read_nircam_ss(self, line, path="./"):
+        """
+        Read in Webb/NIRCam slitless spectra
+
+        Parameters
+        ----------
+        line : string
+            formatted line from DAT file
+            example: 'NIRCam_SS = hd029647_irs.fits'
+
+        path : string, optional
+            location of the FITS files path
+
+        Returns
+        -------
+        Updates self.(file, wave_range, waves, flux, uncs, npts, n_waves)
+        """
+        self.read_spectra(line, path)
+
+        self.fluxes = self.fluxes.to(
+            fluxunit, equivalencies=u.spectral_density(self.waves)
+        )
+        self.uncs = self.uncs.to(fluxunit, equivalencies=u.spectral_density(self.waves))
+
     def read_miri_ifu(self, line, path="./"):
         """
         Read in Webb/MRS IFU spectra
@@ -823,7 +857,7 @@
         ----------
         line : string
             formatted line from DAT file
-            example: 'IRS = hd029647_irs.fits'
+            example: 'MIRI_IFU = hd029647_irs.fits'
 
         path : string, optional
             location of the FITS files path
@@ -834,9 +868,6 @@
         """
         self.read_spectra(line, path)
 
-        # add units
-        self.fluxes = self.fluxes.value * u.Jy
-        self.uncs = self.uncs.value * u.Jy
         self.fluxes = self.fluxes.to(
             fluxunit, equivalencies=u.spectral_density(self.waves)
         )
@@ -934,7 +965,13 @@
     """
 
     def __init__(
-        self, datfile, path="", photonly=False, use_corfac=True, deredden=False
+        self,
+        datfile,
+        path="",
+        photonly=False,
+        use_corfac=True,
+        deredden=False,
+        only_bands=None,
     ):
         """
         Parameters
@@ -955,6 +992,9 @@
         deredden : boolean [default=False]
            Deredden the data based on dereddening parameters given in the DAT file.
            Generally used to deredden standards.
+
+        only_bands : list
+            Only read in the bands given
         """
         self.file = datfile
         self.path = path
@@ -969,9 +1009,9 @@
         self.dereddened = deredden
 
         if self.file is not None:
-            self.read(deredden=deredden)
+            self.read(deredden=deredden, only_bands=only_bands)
 
-    def read(self, deredden=False):
+    def read(self, deredden=False, only_bands=None):
         """
         Populate the object from a DAT file + spectral files
 
@@ -980,6 +1020,8 @@
         deredden : boolean [default=False]
            Deredden the data based on dereddening parameters given in the DAT file.
            Generally used to deredden standards.
+        only_bands : list
+            Only read in the bands given
         """
 
         # open and read all the lines in the file
@@ -988,7 +1030,7 @@
         f.close()
         # get the photometric band data
         self.data["BAND"] = BandData("BAND")
-        self.data["BAND"].read_bands(self.datfile_lines)
+        self.data["BAND"].read_bands(self.datfile_lines, only_bands=only_bands)
 
         # covert the photoemtric band data to fluxes in all possible bands
         self.data["BAND"].get_band_fluxes()
@@ -1093,6 +1135,16 @@
                         )
                     else:
                         warnings.warn(f"{fname} does not exist", UserWarning)
+                elif line.find("NIRCam_SS") == 0:
+                    fname = _getspecfilename(line, self.path)
+                    if os.path.isfile(fname):
+                        self.data["NIRCam_SS"] = SpecData("NIRCam_SS")
+                        self.data["NIRCam_SS"].read_miri_ifu(
+                            line,
+                            path=self.path,
+                        )
+                    else:
+                        warnings.warn(f"{fname} does not exist", UserWarning)
                 elif line.find("MIRI_IFU") == 0:
                     fname = _getspecfilename(line, self.path)
                     if os.path.isfile(fname):

measure_extinction/utils/merge_miri_ifu_spec.py

@@ -4,12 +4,17 @@
 import argparse
 import numpy as np
 import pkg_resources
+import matplotlib.pyplot as plt
 
 from astropy.table import QTable
+import astropy.units as u
 
 from measure_extinction.merge_obsspec import merge_miri_ifu_obsspec
 
 
+fluxunit = u.erg / (u.cm * u.cm * u.s * u.angstrom)
+
+
 if __name__ == "__main__":
 
     # commandline parser
@@ -31,15 +36,15 @@
     parser.add_argument("--pdf", help="save figure as a pdf file", action="store_true")
     args = parser.parse_args()
 
-    sfilename = f"{args.inpath}{args.starname}*.fits"
-    print(sfilename)
+    sfilename = f"{args.inpath}{args.starname}*order*.fits"
     sfiles = glob.glob(sfilename)
-    print(sfiles)
     stable = []
     for cfile in sfiles:
-        print(cfile)
         cdata = QTable.read(cfile)
-        cdata.rename_column("FLUX_ERROR", "ERROR")
+        cdata.rename_column("wavelength", "WAVELENGTH")
+        cdata.rename_column("flux", "FLUX")
+        cdata.rename_column("unc", "ERROR")
+        cdata["WAVELENGTH"] *= u.micron
         cdata["NPTS"] = np.full((len(cdata["FLUX"])), 1.0)
         cdata["NPTS"][cdata["FLUX"] == 0.0] = 0.0
         stable.append(cdata)
@@ -51,3 +56,60 @@
         outname = args.starname.lower()
     mrs_file = f"{outname}_miri_ifu.fits"
     rb_mrs.write(f"{args.outpath}/{mrs_file}", overwrite=True)
+
+    # plot the original and merged Spectra
+    fontsize = 14
+    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)
+
+    fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(10, 5.5))
+
+    for ctable in stable:
+        gvals = ctable["NPTS"] > 0
+        cfluxes = (
+            ctable["FLUX"]
+            .to(fluxunit, equivalencies=u.spectral_density(ctable["WAVELENGTH"]))
+            .value
+        )
+        ax.plot(
+            ctable["WAVELENGTH"][gvals],
+            cfluxes[gvals],
+            "k-",
+            alpha=0.5,
+            label="orig",
+        )
+    gvals = rb_mrs["NPTS"] > 0
+    cfluxes = (
+        rb_mrs["FLUX"]
+        .to(fluxunit, equivalencies=u.spectral_density(ctable["WAVELENGTH"]))
+        .value
+    )
+
+    ax.plot(
+        rb_mrs["WAVELENGTH"][gvals].to(u.micron),
+        rb_mrs["FLUX"][gvals],
+        "b-",
+        alpha=0.5,
+        label="merged",
+    )
+
+    ax.set_xlabel(r"$\lambda$ [$\AA$]")
+    ax.set_ylabel(r"F($\lambda$)")
+
+    ax.set_xscale("log")
+    ax.set_yscale("log")
+
+    ax.legend()
+    fig.tight_layout()
+
+    fname = mrs_file.replace(".fits", "")
+    if args.png:
+        fig.savefig(f"{fname}.png")
+    elif args.pdf:
+        fig.savefig(f"{fname}.pdf")
+    else:
+        plt.show()