From c2e4cbcfb33f9a4801bdf5bb5532267d09344d4f Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Tue, 1 May 2018 17:26:21 -0400
Subject: [PATCH 01/89] First pass at XraySpectrum1D

---
 specutils/spectra/xrayspectrum1d.py | 47 +++++++++++++++++++++++++++++
 1 file changed, 47 insertions(+)
 create mode 100644 specutils/spectra/xrayspectrum1d.py

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
new file mode 100644
index 000000000..b7556ad77
--- /dev/null
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -0,0 +1,47 @@
+import logging
+
+import numpy as np
+from astropy import units as u
+from .spectrum1d import Spectrum1D
+
+__all__ = ['XraySpectrum1D']
+
+# For dealing with varied unit string choices
+EV   = ['eV', 'ev']
+KEV  = ['kev', 'keV']
+ANGS = ['angs', 'Angs', 'Angstrom', 'angstrom', 'Angstroms', 'angstroms', 'A', 'a']
+
+def _unit_parser(unit_string):
+    if unit_string in EV:
+        return u.eV
+    if unit_string in KEV:
+        return u.keV
+    if unit_string in ANGS:
+        return u.angstrom
+
+class XraySpectrum1D(Spectrum1D):
+    """
+    Spectrum container for holding X-ray spectrum
+    WIP by eblur
+
+    Parameters
+    ----------
+    bin_lo
+    bin_hi
+    bin_unit
+    counts
+    arf
+    rmf
+    """
+    def __init__(self, bin_lo, bin_hi, bin_unit, counts, arf=None, rmf=None):
+        try:
+            axis_unit = u.Unit(bin_unit)
+        else:
+            axis_unit = _unit_parser(bin_unit)
+
+        bin_mid = 0.5 * (bin_lo + bin_hi) * unit
+        Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts)
+
+        self.arf = arf
+        self.rmf = rmf
+        return

From 20c8b289b212ab8a00aad758db7149dc9bb34965 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Tue, 1 May 2018 17:40:28 -0400
Subject: [PATCH 02/89] Can load an X-ray spectrum

---
 specutils/spectra/__init__.py       | 1 +
 specutils/spectra/xrayspectrum1d.py | 4 ++--
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/__init__.py b/specutils/spectra/__init__.py
index 711a16154..95c052b1f 100644
--- a/specutils/spectra/__init__.py
+++ b/specutils/spectra/__init__.py
@@ -5,3 +5,4 @@
 from __future__ import absolute_import
 
 from .spectrum1d import *  # noqa
+from .xrayspectrum1d import *
diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index b7556ad77..d57146f97 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -36,10 +36,10 @@ class XraySpectrum1D(Spectrum1D):
     def __init__(self, bin_lo, bin_hi, bin_unit, counts, arf=None, rmf=None):
         try:
             axis_unit = u.Unit(bin_unit)
-        else:
+        except:
             axis_unit = _unit_parser(bin_unit)
 
-        bin_mid = 0.5 * (bin_lo + bin_hi) * unit
+        bin_mid = 0.5 * (bin_lo + bin_hi) * axis_unit
         Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts)
 
         self.arf = arf

From 6ac1c0841b81ec3efe4531423e9ab25cf1415b32 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Tue, 1 May 2018 17:44:41 -0400
Subject: [PATCH 03/89] Added  attribute

---
 specutils/spectra/xrayspectrum1d.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index d57146f97..b939b40d1 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -45,3 +45,8 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, arf=None, rmf=None):
         self.arf = arf
         self.rmf = rmf
         return
+
+    # Convenience function for Xray people
+    @property
+    def counts(self):
+        return self.flux

From 4571350ce927f1dd0934a02b155a9295bb4f138a Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 13:53:39 -0400
Subject: [PATCH 04/89] Added ARF and RMF objects

---
 specutils/spectra/xrayspectrum1d.py | 326 +++++++++++++++++++++++++++-
 1 file changed, 325 insertions(+), 1 deletion(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index b939b40d1..e36a500a6 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -4,7 +4,7 @@
 from astropy import units as u
 from .spectrum1d import Spectrum1D
 
-__all__ = ['XraySpectrum1D']
+__all__ = ['XraySpectrum1D', 'ARF', 'RMF']
 
 # For dealing with varied unit string choices
 EV   = ['eV', 'ev']
@@ -50,3 +50,327 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, arf=None, rmf=None):
     @property
     def counts(self):
         return self.flux
+
+## ----  Supporting response file objects
+
+class RMF(object):
+    def __init__(self, filename):
+        self._load_rmf(filename)
+        pass
+
+    def _load_rmf(self, filename):
+        """
+        Load an RMF from a FITS file.
+
+        Parameters
+        ----------
+        filename : str
+            The file name with the RMF file
+
+        Attributes
+        ----------
+        n_grp : numpy.ndarray
+            the Array with the number of channels in each
+            channel set
+
+        f_chan : numpy.ndarray
+            The starting channel for each channel group;
+            If an element i in n_grp > 1, then the resulting
+            row entry in f_chan will be a list of length n_grp[i];
+            otherwise it will be a single number
+
+        n_chan : numpy.ndarray
+            The number of channels in each channel group. The same
+            logic as for f_chan applies
+
+        matrix : numpy.ndarray
+            The redistribution matrix as a flattened 1D vector
+
+        energ_lo : numpy.ndarray
+            The lower edges of the energy bins
+
+        energ_hi : numpy.ndarray
+            The upper edges of the energy bins
+
+        detchans : int
+            The number of channels in the detector
+
+        """
+        # open the FITS file and extract the MATRIX extension
+        # which contains the redistribution matrix and
+        # anxillary information
+        hdulist = fits.open(filename)
+
+        # get all the extension names
+        extnames = np.array([h.name for h in hdulist])
+
+        # figure out the right extension to use
+        if "MATRIX" in extnames:
+            h = hdulist["MATRIX"]
+
+        elif "SPECRESP MATRIX" in extnames:
+            h = hdulist["SPECRESP MATRIX"]
+
+        data = h.data
+        hdr = h.header
+        hdulist.close()
+
+        # extract + store the attributes described in the docstring
+        n_grp = np.array(data.field("N_GRP"))
+        f_chan = np.array(data.field('F_CHAN'))
+        n_chan = np.array(data.field("N_CHAN"))
+        matrix = np.array(data.field("MATRIX"))
+
+        self.energ_lo = np.array(data.field("ENERG_LO"))
+        self.energ_hi = np.array(data.field("ENERG_HI"))
+        self.energ_unit = data.columns["ENERG_LO"].unit
+        self.detchans = hdr["DETCHANS"]
+        self.offset = self.__get_tlmin(h)
+
+        # flatten the variable-length arrays
+        self.n_grp, self.f_chan, self.n_chan, self.matrix = \
+                self._flatten_arrays(n_grp, f_chan, n_chan, matrix)
+
+        return
+
+    def __get_tlmin(self, h):
+        """
+        Get the tlmin keyword for `F_CHAN`.
+
+        Parameters
+        ----------
+        h : an astropy.io.fits.hdu.table.BinTableHDU object
+            The extension containing the `F_CHAN` column
+
+        Returns
+        -------
+        tlmin : int
+            The tlmin keyword
+        """
+        # get the header
+        hdr = h.header
+        # get the keys of all
+        keys = np.array(list(hdr.keys()))
+
+        # find the place where the tlmin keyword is defined
+        t = np.array(["TLMIN" in k for k in keys])
+
+        # get the index of the TLMIN keyword
+        tlmin_idx = np.hstack(np.where(t))[0]
+
+        # get the corresponding value
+        tlmin = np.int(list(hdr.items())[tlmin_idx][1])
+
+        return tlmin
+
+    def _flatten_arrays(self, n_grp, f_chan, n_chan, matrix):
+
+        if not len(n_grp) == len(f_chan) == len(n_chan) == len(matrix):
+            raise ValueError("Arrays must be of same length!")
+
+        # find all non-zero groups
+        nz_idx = (n_grp > 0)
+
+        # stack all non-zero rows in the matrix
+        matrix_flat = np.hstack(matrix[nz_idx])
+
+        # stack all nonzero rows in n_chan and f_chan
+        #n_chan_flat = np.hstack(n_chan[nz_idx])
+        #f_chan_flat = np.hstack(f_chan[nz_idx])
+
+        # some matrices actually have more elements
+        # than groups in `n_grp`, so we'll only pick out
+        # those values that have a correspondence in
+        # n_grp
+        f_chan_new = []
+        n_chan_new = []
+        for i,t in enumerate(nz_idx):
+            if t:
+                n = n_grp[i]
+                f = f_chan[i]
+                nc = n_chan[i]
+                if np.size(f) == 1:
+                    f_chan_new.append(f)
+                    n_chan_new.append(nc)
+                else:
+                    f_chan_new.append(f[:n])
+                    n_chan_new.append(nc[:n])
+
+        n_chan_flat = np.hstack(n_chan_new)
+        f_chan_flat = np.hstack(f_chan_new)
+
+        # if n_chan is zero, we'll remove those as well.
+        nz_idx2 = (n_chan_flat > 0)
+        n_chan_flat = n_chan_flat[nz_idx2]
+        f_chan_flat = f_chan_flat[nz_idx2]
+
+        return n_grp, f_chan_flat, n_chan_flat, matrix_flat
+
+    def apply_rmf(self, spec):
+        """
+        Fold the spectrum through the redistribution matrix.
+
+        The redistribution matrix is saved as a flattened 1-dimensional
+        vector to save space. In reality, for each entry in the flux
+        vector, there exists one or more sets of channels that this
+        flux is redistributed into. The additional arrays `n_grp`,
+        `f_chan` and `n_chan` store this information:
+            * `n_group` stores the number of channel groups for each
+              energy bin
+            * `f_chan` stores the *first channel* that each channel
+              for each channel set
+            * `n_chan` stores the number of channels in each channel
+              set
+
+        As a result, for a given energy bin i, we need to look up the
+        number of channel sets in `n_grp` for that energy bin. We
+        then need to loop over the number of channel sets. For each
+        channel set, we look up the first channel into which flux
+        will be distributed as well as the number of channels in the
+        group. We then need to also loop over the these channels and
+        actually use the corresponding elements in the redistribution
+        matrix to redistribute the photon flux into channels.
+
+        All of this is basically a big bookkeeping exercise in making
+        sure to get the indices right.
+
+        Parameters
+        ----------
+        spec : numpy.ndarray
+            The (model) spectrum to be folded
+
+        Returns
+        -------
+        counts : numpy.ndarray
+            The (model) spectrum after folding, in
+            counts/s/channel
+
+        """
+        # get the number of channels in the data
+        nchannels = spec.shape[0]
+
+        # an empty array for the output counts
+        counts = np.zeros(nchannels)
+
+        # index for n_chan and f_chan incrementation
+        k = 0
+
+        # index for the response matrix incrementation
+        resp_idx = 0
+
+        # loop over all channels
+        for i in range(nchannels):
+
+            # this is the current bin in the flux spectrum to
+            # be folded
+            source_bin_i = spec[i]
+
+            # get the current number of groups
+            current_num_groups = self.n_grp[i]
+
+            # loop over the current number of groups
+            for j in range(current_num_groups):
+
+                current_num_chans = int(self.n_chan[k])
+
+                if current_num_chans == 0:
+                    k += 1
+                    resp_idx += current_num_chans
+                    continue
+
+
+                else:
+                    # get the right index for the start of the counts array
+                    # to put the data into
+                    counts_idx = int(self.f_chan[k] - self.offset)
+                    # this is the current number of channels to use
+
+                    k += 1
+                    # add the flux to the subarray of the counts array that starts with
+                    # counts_idx and runs over current_num_chans channels
+                    counts[counts_idx:counts_idx +
+                                      current_num_chans] += self.matrix[resp_idx:resp_idx +
+                                                                                 current_num_chans] * \
+                                                                np.float(source_bin_i)
+                    # iterate the response index for next round
+                    resp_idx += current_num_chans
+
+
+        return counts[:self.detchans]
+
+
+class ARF(object):
+
+    def __init__(self, filename):
+
+        self._load_arf(filename)
+        pass
+
+    def _load_arf(self, filename):
+        """
+        Load an ARF from a FITS file.
+        Parameters
+        ----------
+        filename : str
+            The file name with the RMF file
+        Attributes
+        ----------
+        """
+        # open the FITS file and extract the MATRIX extension
+        # which contains the redistribution matrix and
+        # anxillary information
+        hdulist = fits.open(filename)
+        h = hdulist["SPECRESP"]
+        data = h.data
+        hdr = h.header
+        hdulist.close()
+
+        # extract + store the attributes described in the docstring
+
+        self.e_low  = np.array(data.field("ENERG_LO"))
+        self.e_high = np.array(data.field("ENERG_HI"))
+        self.e_unit = data.columns["ENERG_LO"].unit
+        self.specresp = np.array(data.field("SPECRESP"))
+
+        if "EXPOSURE" in list(hdr.keys()):
+            self.exposure = hdr["EXPOSURE"]
+        else:
+            self.exposure = 1.0
+
+        if "FRACEXPO" in data.columns.names:
+            self.fracexpo = data["FRACEXPO"]
+        else:
+            self.fracexpo = 1.0
+
+        return
+
+    def apply_arf(self, spec, exposure=None):
+        """
+        Fold the spectrum through the ARF.
+        The ARF is a single vector encoding the effective area information
+        about the detector. A such, applying the ARF is a simple
+        multiplication with the input spectrum.
+        Parameters
+        ----------
+        spec : numpy.ndarray
+            The (model) spectrum to be folded
+        exposure : float, default None
+            Value for the exposure time. By default, `apply_arf` will use the
+            exposure keyword from the ARF file. If this exposure time is not
+            correct (for example when simulated spectra use a different exposure
+            time and the ARF from a real observation), one can override the
+            default exposure by setting the `exposure` keyword to the correct
+            value.
+        Returns
+        -------
+        s_arf : numpy.ndarray
+            The (model) spectrum after folding, in
+            counts/s/channel
+        """
+        assert spec.shape[0] == self.specresp.shape[0], "The input spectrum must " \
+                                                      "be of same size as the " \
+                                                      "ARF array."
+        if exposure is None:
+            return np.array(spec) * self.specresp * self.exposure
+        else:
+            return np.array(spec) * self.specresp * exposure

From a92bd09751e7595466ed49bd407d0a94216b585a Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 14:07:50 -0400
Subject: [PATCH 05/89] Added assign_rmf and assign_arf methods

---
 specutils/spectra/xrayspectrum1d.py | 29 +++++++++++++++++++++++++++--
 1 file changed, 27 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index e36a500a6..a371b94a5 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -42,8 +42,8 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, arf=None, rmf=None):
         bin_mid = 0.5 * (bin_lo + bin_hi) * axis_unit
         Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts)
 
-        self.arf = arf
-        self.rmf = rmf
+        self.assign_rmf(rmf)
+        self.assign_arf(rmf)
         return
 
     # Convenience function for Xray people
@@ -51,6 +51,21 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, arf=None, rmf=None):
     def counts(self):
         return self.flux
 
+    def assign_arf(self, arf_inp):
+        if isinstance(arf_inp, str):
+            self.arf = ARF(arf_inp)
+        else:
+            self.arf = arf_inp
+        return
+
+    def assign_rmf(self, rmf_inp):
+        if isinstance(rmf_inp, str):
+            self.rmf = RMF(rmf_inp)
+        else:
+            self.rmf = rmf_inp
+        return
+
+
 ## ----  Supporting response file objects
 
 class RMF(object):
@@ -69,6 +84,9 @@ def _load_rmf(self, filename):
 
         Attributes
         ----------
+        filename : str
+            The file name that the RMF was drawn from
+
         n_grp : numpy.ndarray
             the Array with the number of channels in each
             channel set
@@ -100,6 +118,7 @@ def _load_rmf(self, filename):
         # which contains the redistribution matrix and
         # anxillary information
         hdulist = fits.open(filename)
+        self.filename = filename
 
         # get all the extension names
         extnames = np.array([h.name for h in hdulist])
@@ -309,17 +328,23 @@ def __init__(self, filename):
     def _load_arf(self, filename):
         """
         Load an ARF from a FITS file.
+
         Parameters
         ----------
         filename : str
             The file name with the RMF file
+
         Attributes
         ----------
+        filename : str
+            The file name that the ARF was drawn from
         """
         # open the FITS file and extract the MATRIX extension
         # which contains the redistribution matrix and
         # anxillary information
         hdulist = fits.open(filename)
+        self.filename = filename
+        
         h = hdulist["SPECRESP"]
         data = h.data
         hdr = h.header

From 7077dc7187ed413b6e39040b821d57aaa645243d Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 14:16:13 -0400
Subject: [PATCH 06/89] Added exposure input

---
 specutils/spectra/xrayspectrum1d.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index a371b94a5..51de9e726 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -33,7 +33,7 @@ class XraySpectrum1D(Spectrum1D):
     arf
     rmf
     """
-    def __init__(self, bin_lo, bin_hi, bin_unit, counts, arf=None, rmf=None):
+    def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure, arf=None, rmf=None):
         try:
             axis_unit = u.Unit(bin_unit)
         except:
@@ -42,6 +42,7 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, arf=None, rmf=None):
         bin_mid = 0.5 * (bin_lo + bin_hi) * axis_unit
         Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts)
 
+        self.exposure = exposure
         self.assign_rmf(rmf)
         self.assign_arf(rmf)
         return
@@ -344,7 +345,7 @@ def _load_arf(self, filename):
         # anxillary information
         hdulist = fits.open(filename)
         self.filename = filename
-        
+
         h = hdulist["SPECRESP"]
         data = h.data
         hdr = h.header

From 9c8159df3f7b5ce2583459056fca6b7eb77a6a66 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 14:18:22 -0400
Subject: [PATCH 07/89] Fixed bugs

---
 specutils/spectra/xrayspectrum1d.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 51de9e726..fade12358 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -2,6 +2,7 @@
 
 import numpy as np
 from astropy import units as u
+from astropy.io import fits
 from .spectrum1d import Spectrum1D
 
 __all__ = ['XraySpectrum1D', 'ARF', 'RMF']
@@ -44,7 +45,7 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure, arf=None, rmf=Non
 
         self.exposure = exposure
         self.assign_rmf(rmf)
-        self.assign_arf(rmf)
+        self.assign_arf(arf)
         return
 
     # Convenience function for Xray people

From fb5f1b6ed5c088bb3aa7ce0e1787dd25f156ad42 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 14:30:44 -0400
Subject: [PATCH 08/89] Updated docstrings

---
 specutils/spectra/xrayspectrum1d.py | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index fade12358..72c5bd11e 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -112,6 +112,9 @@ def _load_rmf(self, filename):
         energ_hi : numpy.ndarray
             The upper edges of the energy bins
 
+        energ_unit : str
+            Description of the energy units used
+
         detchans : int
             The number of channels in the detector
 
@@ -340,6 +343,28 @@ def _load_arf(self, filename):
         ----------
         filename : str
             The file name that the ARF was drawn from
+
+        e_low : numpy.ndarray
+            The lower edges of the energy bins
+
+        e_high : numpy.ndarray
+            The upper edges of the energy bins
+
+        e_unit : str
+            Description of the energy units used
+
+        specresp : numpy.ndarray
+            Description of the energy dependent telescope response area
+
+        fracexpo : float or numpy.ndarray
+            Fractional exposure time for the spectrum
+            (sometimes constant, sometimes dependent on spectral channel).
+            These values are stored for reference; generally, they are already
+            accounted for in the specresp array.
+
+        exposure :
+            Average exposure time for the dataset
+            (takes telescope dithering into account)
         """
         # open the FITS file and extract the MATRIX extension
         # which contains the redistribution matrix and

From a08d8ae03050f359b7e2dd9b647939133adc2146 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 14:36:08 -0400
Subject: [PATCH 09/89] Store bin and energy units as astropy.units.Unit

---
 specutils/spectra/xrayspectrum1d.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 72c5bd11e..2b98164aa 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -31,6 +31,7 @@ class XraySpectrum1D(Spectrum1D):
     bin_hi
     bin_unit
     counts
+    exposure
     arf
     rmf
     """
@@ -43,6 +44,7 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure, arf=None, rmf=Non
         bin_mid = 0.5 * (bin_lo + bin_hi) * axis_unit
         Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts)
 
+        self.bin_unit = axis_unit # keep for reference
         self.exposure = exposure
         self.assign_rmf(rmf)
         self.assign_arf(arf)
@@ -112,7 +114,7 @@ def _load_rmf(self, filename):
         energ_hi : numpy.ndarray
             The upper edges of the energy bins
 
-        energ_unit : str
+        energ_unit : astropy.units.Unit
             Description of the energy units used
 
         detchans : int
@@ -147,7 +149,7 @@ def _load_rmf(self, filename):
 
         self.energ_lo = np.array(data.field("ENERG_LO"))
         self.energ_hi = np.array(data.field("ENERG_HI"))
-        self.energ_unit = data.columns["ENERG_LO"].unit
+        self.energ_unit = _unit_parser(data.columns["ENERG_LO"].unit)
         self.detchans = hdr["DETCHANS"]
         self.offset = self.__get_tlmin(h)
 
@@ -350,7 +352,7 @@ def _load_arf(self, filename):
         e_high : numpy.ndarray
             The upper edges of the energy bins
 
-        e_unit : str
+        e_unit : astropy.units.Unit
             Description of the energy units used
 
         specresp : numpy.ndarray
@@ -381,7 +383,7 @@ def _load_arf(self, filename):
 
         self.e_low  = np.array(data.field("ENERG_LO"))
         self.e_high = np.array(data.field("ENERG_HI"))
-        self.e_unit = data.columns["ENERG_LO"].unit
+        self.e_unit = _unit_parser(data.columns["ENERG_LO"].unit)
         self.specresp = np.array(data.field("SPECRESP"))
 
         if "EXPOSURE" in list(hdr.keys()):

From 20ee61df2a27bd8ed6b0a1448393906d3d27130a Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 16:31:14 -0400
Subject: [PATCH 10/89] Added bin_lo and hi because wcs.bin_edges did not exist

---
 specutils/spectra/xrayspectrum1d.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 2b98164aa..d209c988b 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -35,15 +35,18 @@ class XraySpectrum1D(Spectrum1D):
     arf
     rmf
     """
-    def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure, arf=None, rmf=None):
+    def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
+                 arf=None, rmf=None, **kwargs):
         try:
             axis_unit = u.Unit(bin_unit)
         except:
             axis_unit = _unit_parser(bin_unit)
 
         bin_mid = 0.5 * (bin_lo + bin_hi) * axis_unit
-        Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts)
+        Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts, **kwargs)
 
+        self.bin_lo = bin_lo
+        self.bin_hi = bin_hi
         self.bin_unit = axis_unit # keep for reference
         self.exposure = exposure
         self.assign_rmf(rmf)

From 6cd49f0ccb7290f97e6812f5a67cd1b1217050e2 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 16:34:35 -0400
Subject: [PATCH 11/89] Added apply_resp method

---
 specutils/spectra/xrayspectrum1d.py | 44 +++++++++++++++++++++++++++++
 1 file changed, 44 insertions(+)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index d209c988b..765e77241 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -72,6 +72,50 @@ def assign_rmf(self, rmf_inp):
             self.rmf = rmf_inp
         return
 
+    def apply_resp(self, mflux, exposure=None):
+        """
+        Given a model flux spectrum, apply the response. In cases where the
+        spectrum has both an ARF and an RMF, apply both. Otherwise, apply
+        whatever response is in RMF.
+
+        The model flux spectrum *must* be created using the same units and
+        bins as in the ARF (where the ARF exists)!
+
+        Parameters
+        ----------
+        mflux : iterable
+            A list or array with the model flux values in ergs/keV/s/cm^-2
+
+        exposure : float, default None
+            By default, the exposure stored in the ARF will be used to compute
+            the total counts per bin over the effective observation time.
+            In cases where this might be incorrect (e.g. for simulated spectra
+            where the pha file might have a different exposure value than the
+            ARF), this keyword provides the functionality to override the
+            default behaviour and manually set the exposure time to use.
+
+        Returns
+        -------
+        count_model : numpy.ndarray
+            The model spectrum in units of counts/bin
+
+        If no ARF file exists, it will return the model flux after applying the RMF
+        If no RMF file exists, it will return the model flux after applying the ARF (with a warning)
+        If no ARF and no RMF, it will return the model flux spectrum (with a warning)
+        """
+
+        if self.arf is not None:
+            mrate  = self.arf.apply_arf(mflux, exposure=exposure)
+        else:
+            mrate = mflux
+
+        if self.rmf is not None:
+            count_model = self.rmf.apply_rmf(mrate)
+            return count_model
+        else:
+            print("Caution: no response file specified")
+            return mrate
+
 
 ## ----  Supporting response file objects
 

From 9c08caaf324c4258e61b046b36d2b0f7ddea0eb4 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 16:52:47 -0400
Subject: [PATCH 12/89] Adjust apply_resp output if spectrum is stored by
 wavelength

---
 specutils/spectra/xrayspectrum1d.py | 13 ++++++++++---
 1 file changed, 10 insertions(+), 3 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 765e77241..b9f2df62d 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -110,11 +110,18 @@ def apply_resp(self, mflux, exposure=None):
             mrate = mflux
 
         if self.rmf is not None:
-            count_model = self.rmf.apply_rmf(mrate)
-            return count_model
+            result = self.rmf.apply_rmf(mrate)
         else:
             print("Caution: no response file specified")
-            return mrate
+            result = mrate
+
+        # The response is always stored according to energy
+        # If the spectrum is stored in wavelength, reverse the result so it
+        # matches in wavelength space
+        if self.bin_unit.physical_type == 'length':
+            return result[::-1]
+        else:
+            return result
 
 
 ## ----  Supporting response file objects

From 35d09321cf066a833adbec17e940dcf94c3b43f3 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 18:03:29 -0400
Subject: [PATCH 13/89] Don't need bin_unit

---
 specutils/spectra/xrayspectrum1d.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index b9f2df62d..bac3f6e3e 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -47,7 +47,6 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
 
         self.bin_lo = bin_lo
         self.bin_hi = bin_hi
-        self.bin_unit = axis_unit # keep for reference
         self.exposure = exposure
         self.assign_rmf(rmf)
         self.assign_arf(arf)
@@ -118,7 +117,8 @@ def apply_resp(self, mflux, exposure=None):
         # The response is always stored according to energy
         # If the spectrum is stored in wavelength, reverse the result so it
         # matches in wavelength space
-        if self.bin_unit.physical_type == 'length':
+        if self.spectral_axis.unit.physical_type == 'length':
+            print("reversing output")
             return result[::-1]
         else:
             return result

From da39e47faa109c2d18cf215955aacfa500a62c55 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 18:12:10 -0400
Subject: [PATCH 14/89] Because pha-file wavelength is in reverse order, the
 output from apply_resp is in correct order. No longer need to check

---
 specutils/spectra/xrayspectrum1d.py | 17 ++++++++---------
 1 file changed, 8 insertions(+), 9 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index bac3f6e3e..c054255df 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -34,6 +34,11 @@ class XraySpectrum1D(Spectrum1D):
     exposure
     arf
     rmf
+
+    Attributes
+    ----------
+    model : numpy.ndarray
+        Stores model values
     """
     def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
                  arf=None, rmf=None, **kwargs):
@@ -50,6 +55,7 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
         self.exposure = exposure
         self.assign_rmf(rmf)
         self.assign_arf(arf)
+        self.model_counts = np.zeros_like(bin_mid)
         return
 
     # Convenience function for Xray people
@@ -114,15 +120,8 @@ def apply_resp(self, mflux, exposure=None):
             print("Caution: no response file specified")
             result = mrate
 
-        # The response is always stored according to energy
-        # If the spectrum is stored in wavelength, reverse the result so it
-        # matches in wavelength space
-        if self.spectral_axis.unit.physical_type == 'length':
-            print("reversing output")
-            return result[::-1]
-        else:
-            return result
-
+        self.model_counts = result
+        return result
 
 ## ----  Supporting response file objects
 

From c5db843592a0d3d525f1bd2cdd587bf157317246 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 18:15:20 -0400
Subject: [PATCH 15/89] added save_model_counts option

---
 specutils/spectra/xrayspectrum1d.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index c054255df..cd78d2e49 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -77,7 +77,7 @@ def assign_rmf(self, rmf_inp):
             self.rmf = rmf_inp
         return
 
-    def apply_resp(self, mflux, exposure=None):
+    def apply_resp(self, mflux, exposure=None, store_model_counts=True):
         """
         Given a model flux spectrum, apply the response. In cases where the
         spectrum has both an ARF and an RMF, apply both. Otherwise, apply
@@ -99,9 +99,12 @@ def apply_resp(self, mflux, exposure=None):
             ARF), this keyword provides the functionality to override the
             default behaviour and manually set the exposure time to use.
 
+        store_model_counts : bool
+            If True, the output will also be stored in self.model_counts
+
         Returns
         -------
-        count_model : numpy.ndarray
+        model_counts : numpy.ndarray
             The model spectrum in units of counts/bin
 
         If no ARF file exists, it will return the model flux after applying the RMF
@@ -120,7 +123,8 @@ def apply_resp(self, mflux, exposure=None):
             print("Caution: no response file specified")
             result = mrate
 
-        self.model_counts = result
+        if store_model_counts:
+            self.model_counts = result
         return result
 
 ## ----  Supporting response file objects

From 291be38f142345a331d996655c8816ee4feea696 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Thu, 3 May 2018 10:06:41 -0400
Subject: [PATCH 16/89] Changed my mind on model_counts thing for now

---
 specutils/spectra/xrayspectrum1d.py | 10 +---------
 1 file changed, 1 insertion(+), 9 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index cd78d2e49..7b342cb96 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -34,11 +34,6 @@ class XraySpectrum1D(Spectrum1D):
     exposure
     arf
     rmf
-
-    Attributes
-    ----------
-    model : numpy.ndarray
-        Stores model values
     """
     def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
                  arf=None, rmf=None, **kwargs):
@@ -55,7 +50,6 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
         self.exposure = exposure
         self.assign_rmf(rmf)
         self.assign_arf(arf)
-        self.model_counts = np.zeros_like(bin_mid)
         return
 
     # Convenience function for Xray people
@@ -77,7 +71,7 @@ def assign_rmf(self, rmf_inp):
             self.rmf = rmf_inp
         return
 
-    def apply_resp(self, mflux, exposure=None, store_model_counts=True):
+    def apply_resp(self, mflux, exposure=None):
         """
         Given a model flux spectrum, apply the response. In cases where the
         spectrum has both an ARF and an RMF, apply both. Otherwise, apply
@@ -123,8 +117,6 @@ def apply_resp(self, mflux, exposure=None, store_model_counts=True):
             print("Caution: no response file specified")
             result = mrate
 
-        if store_model_counts:
-            self.model_counts = result
         return result
 
 ## ----  Supporting response file objects

From fcd9e25667931d27b1597ca1826429975481d354 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 23 May 2018 18:27:59 -0500
Subject: [PATCH 17/89] Hard coded in a rest_value

---
 specutils/spectra/xrayspectrum1d.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 7b342cb96..2c5f4d86e 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -36,14 +36,15 @@ class XraySpectrum1D(Spectrum1D):
     rmf
     """
     def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
-                 arf=None, rmf=None, **kwargs):
+                 arf=None, rmf=None, rest_value=0.0, **kwargs):
         try:
             axis_unit = u.Unit(bin_unit)
         except:
             axis_unit = _unit_parser(bin_unit)
 
         bin_mid = 0.5 * (bin_lo + bin_hi) * axis_unit
-        Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts, **kwargs)
+        Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts,
+                            rest_value=rest_value, **kwargs)
 
         self.bin_lo = bin_lo
         self.bin_hi = bin_hi

From fe8c1b5beec19de16bb86995a7a72af8842a6ef0 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 23 May 2018 18:29:31 -0500
Subject: [PATCH 18/89] Added unit to rest_value

---
 specutils/spectra/xrayspectrum1d.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 2c5f4d86e..d700837b7 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -36,7 +36,7 @@ class XraySpectrum1D(Spectrum1D):
     rmf
     """
     def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
-                 arf=None, rmf=None, rest_value=0.0, **kwargs):
+                 arf=None, rmf=None, rest_value=0.0 * u.angstrom, **kwargs):
         try:
             axis_unit = u.Unit(bin_unit)
         except:

From 1c7f27d1eb41ccdafe874565cfc724c17d8885a8 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 23 May 2018 18:47:34 -0500
Subject: [PATCH 19/89] Commented out distracting print statements

---
 specutils/spectra/xrayspectrum1d.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index d700837b7..ceb3f9fdf 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -110,12 +110,13 @@ def apply_resp(self, mflux, exposure=None):
         if self.arf is not None:
             mrate  = self.arf.apply_arf(mflux, exposure=exposure)
         else:
+            #print("Caution: no ARF file specified")
             mrate = mflux
 
         if self.rmf is not None:
             result = self.rmf.apply_rmf(mrate)
         else:
-            print("Caution: no response file specified")
+            #print("Caution: no RMF file specified")
             result = mrate
 
         return result

From c6f5aaa55f4cb784def23effb8f50a4f54ec6f69 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 23 May 2018 19:00:00 -0500
Subject: [PATCH 20/89] Wrote tests

---
 specutils/tests/test_xrayspectrum1d.py | 35 ++++++++++++++++++++++++++
 1 file changed, 35 insertions(+)
 create mode 100644 specutils/tests/test_xrayspectrum1d.py

diff --git a/specutils/tests/test_xrayspectrum1d.py b/specutils/tests/test_xrayspectrum1d.py
new file mode 100644
index 000000000..f8792a7ff
--- /dev/null
+++ b/specutils/tests/test_xrayspectrum1d.py
@@ -0,0 +1,35 @@
+import numpy as np
+from scipy import optimize
+import astropy.units as u
+from astropy.modeling.powerlaws import PowerLaw1D
+from specutils import XraySpectrum1D
+
+def test_create_from_arrays():
+    # Test that XraySpectrum1D can be initialized
+    amp0, alpha0 = 3.e-3, 2.0
+    powlaw0 = PowerLaw1D(amplitude=amp0, alpha=alpha0, x_0=1.e3)
+    energy = np.linspace(0.2, 10.0, 8000)
+    elo = energy[:-1]
+    ehi = energy[1:]
+    emid = 0.5 * (elo + ehi)
+    counts = np.random.poisson(lam=powlaw0(emid), size=len(emid))
+    test_spec = XraySpectrum1D(elo, ehi, u.keV, counts, exposure=1.0)
+    return test_spec
+
+def test_call_spectrum():
+    # Test that the XraySpectrum1D has attributes that make it unique
+    # to X-ray spectrum
+    test_spec = test_create_from_arrays()
+    assert len(test_spec.counts) == len(test_spec.spectral_axis)
+    assert len(test_spec.bin_lo) == len(test_spec.spectral_axis)
+    assert len(test_spec.bin_hi) == len(test_spec.spectral_axis)
+    return
+
+def test_apply_model():
+    # Test that one can evaluate XraySpectrum1D with a model
+    test_spec = test_create_from_arrays()
+    new_model = PowerLaw1D(amplitude=5.e-3, alpha=1.8, x_0=1.e3)
+    model_flux = new_model(test_spec.spectral_axis.value)
+    ymodel = test_spec.apply_resp(model_flux)
+    assert len(ymodel) == len(test_spec.spectral_axis)
+    return

From c1929350e17a1a62917e69834acc5cd042fcbacc Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 23 May 2018 19:17:48 -0500
Subject: [PATCH 21/89] Added docstrings

---
 specutils/spectra/xrayspectrum1d.py | 67 +++++++++++++++++++++++++++--
 1 file changed, 64 insertions(+), 3 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index ceb3f9fdf..0b72ccedf 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -13,6 +13,9 @@
 ANGS = ['angs', 'Angs', 'Angstrom', 'angstrom', 'Angstroms', 'angstroms', 'A', 'a']
 
 def _unit_parser(unit_string):
+    """
+    A function for matching unit strings to the correct Astropy unit
+    """
     if unit_string in EV:
         return u.eV
     if unit_string in KEV:
@@ -22,17 +25,51 @@ def _unit_parser(unit_string):
 
 class XraySpectrum1D(Spectrum1D):
     """
-    Spectrum container for holding X-ray spectrum
-    WIP by eblur
+    Spectrum properties specific to X-ray data
 
     Parameters
     ----------
+    bin_lo : numpy.ndarray
+        The left edges for bin values
+
+    bin_hi : numpy.ndarray
+        The right edges for bin values
+
+    bin_unit : string OR astropy.units.Unit
+        Unit for the bin values
+
+    counts : numpy.ndarray
+        Counts histogram for the X-ray spectrum
+
+    exposure : float
+        Exposure time for the dataset
+
+    arf : specutils.ARF or string, default None
+        Strings will be passed to ARF.__init__
+        All other input types are stored as arf attribute
+
+    rmf : specutils.RMF or string, default None
+        Strings will be passed to RMF.__init__
+        All other input types are stored as rmf attribute
+
+    rest_value : astropy.units.Quantity, default 0 Angstrom
+        See Spectrum1D rest_value input
+
+    Attributes
+    ----------
+    Inherits from Spectrum1D object.
+    The following inputs are stored as additional attributes.
+
     bin_lo
+
     bin_hi
-    bin_unit
+
     counts
+
     exposure
+
     arf
+
     rmf
     """
     def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
@@ -59,6 +96,18 @@ def counts(self):
         return self.flux
 
     def assign_arf(self, arf_inp):
+        """
+        Assign an ARF object to the XraySpectrum1D object
+
+        Input
+        -----
+        arf_inp : string
+            File name for the ARF (FITS file)
+
+        Returns
+        -------
+        Modifies the XraySpectrum1D.arf attribute
+        """
         if isinstance(arf_inp, str):
             self.arf = ARF(arf_inp)
         else:
@@ -66,6 +115,18 @@ def assign_arf(self, arf_inp):
         return
 
     def assign_rmf(self, rmf_inp):
+        """
+        Assign an RMF object to the XraySpectrum1D object
+
+        Input
+        -----
+        rmf_inp : string
+            File name for the RMF (FITS file)
+
+        Returns
+        -------
+        Modifies the XraySpectrum1D.rmf attribute
+        """
         if isinstance(rmf_inp, str):
             self.rmf = RMF(rmf_inp)
         else:

From 3314f7dfd7c30ea7389495e129f763091c319990 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 23 May 2018 19:19:32 -0500
Subject: [PATCH 22/89] Removed scipy dependence

---
 specutils/tests/test_xrayspectrum1d.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/specutils/tests/test_xrayspectrum1d.py b/specutils/tests/test_xrayspectrum1d.py
index f8792a7ff..53fbae5c0 100644
--- a/specutils/tests/test_xrayspectrum1d.py
+++ b/specutils/tests/test_xrayspectrum1d.py
@@ -1,5 +1,4 @@
 import numpy as np
-from scipy import optimize
 import astropy.units as u
 from astropy.modeling.powerlaws import PowerLaw1D
 from specutils import XraySpectrum1D

From 801e3b661213eb7ba273bf69d415204b2dbf5a8a Mon Sep 17 00:00:00 2001
From: eblur <lia@astro.wisc.edu>
Date: Fri, 1 Jun 2018 10:56:46 -0500
Subject: [PATCH 23/89] Spelled out acronyms, link to edu materials

---
 specutils/spectra/xrayspectrum1d.py | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 0b72ccedf..05e9cc892 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -12,6 +12,9 @@
 KEV  = ['kev', 'keV']
 ANGS = ['angs', 'Angs', 'Angstrom', 'angstrom', 'Angstroms', 'angstroms', 'A', 'a']
 
+## An introduction to X-ray analysis can be found here:
+## http://cxc.cfa.harvard.edu/xrayschool/talks/intro_xray_analysis.pdf
+
 def _unit_parser(unit_string):
     """
     A function for matching unit strings to the correct Astropy unit
@@ -97,12 +100,12 @@ def counts(self):
 
     def assign_arf(self, arf_inp):
         """
-        Assign an ARF object to the XraySpectrum1D object
+        Assign an area response file (ARF) object to the XraySpectrum1D object
 
         Input
         -----
         arf_inp : string
-            File name for the ARF (FITS file)
+            File name for the area response file (FITS file)
 
         Returns
         -------
@@ -116,12 +119,12 @@ def assign_arf(self, arf_inp):
 
     def assign_rmf(self, rmf_inp):
         """
-        Assign an RMF object to the XraySpectrum1D object
+        Assign a response matrix file (RMF) object to the XraySpectrum1D object
 
         Input
         -----
         rmf_inp : string
-            File name for the RMF (FITS file)
+            File name for the response matrix file (FITS file)
 
         Returns
         -------
@@ -191,7 +194,7 @@ def __init__(self, filename):
 
     def _load_rmf(self, filename):
         """
-        Load an RMF from a FITS file.
+        Load a response matrix file (RMF) from a FITS file.
 
         Parameters
         ----------
@@ -446,7 +449,7 @@ def __init__(self, filename):
 
     def _load_arf(self, filename):
         """
-        Load an ARF from a FITS file.
+        Load an area response file (ARF) from a FITS file.
 
         Parameters
         ----------
@@ -512,7 +515,7 @@ def _load_arf(self, filename):
 
     def apply_arf(self, spec, exposure=None):
         """
-        Fold the spectrum through the ARF.
+        Fold the spectrum through the area response file (ARF).
         The ARF is a single vector encoding the effective area information
         about the detector. A such, applying the ARF is a simple
         multiplication with the input spectrum.

From 4dff31686315a0342673ae173534a6d8f59a0ff7 Mon Sep 17 00:00:00 2001
From: eblur <lia@astro.wisc.edu>
Date: Fri, 1 Jun 2018 11:05:22 -0500
Subject: [PATCH 24/89] Added option to specify RMF FITS extension for response
 matrix

---
 specutils/spectra/xrayspectrum1d.py | 18 +++++++++++++++---
 1 file changed, 15 insertions(+), 3 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 05e9cc892..68460d8b6 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -136,7 +136,7 @@ def assign_rmf(self, rmf_inp):
             self.rmf = rmf_inp
         return
 
-    def apply_resp(self, mflux, exposure=None):
+    def apply_response(self, mflux, exposure=None):
         """
         Given a model flux spectrum, apply the response. In cases where the
         spectrum has both an ARF and an RMF, apply both. Otherwise, apply
@@ -192,7 +192,7 @@ def __init__(self, filename):
         self._load_rmf(filename)
         pass
 
-    def _load_rmf(self, filename):
+    def _load_rmf(self, filename, extension=None):
         """
         Load a response matrix file (RMF) from a FITS file.
 
@@ -201,6 +201,10 @@ def _load_rmf(self, filename):
         filename : str
             The file name with the RMF file
 
+        extension : str (default None)
+            FITS file extinction keyword, if the response matrix is stored
+            under an extension other tham "MATRIX" or "SPECRESP MATRIX"
+
         Attributes
         ----------
         filename : str
@@ -246,12 +250,20 @@ def _load_rmf(self, filename):
         extnames = np.array([h.name for h in hdulist])
 
         # figure out the right extension to use
-        if "MATRIX" in extnames:
+        if extension is not None:
+            h = hdulist[extension]
+
+        elif "MATRIX" in extnames:
             h = hdulist["MATRIX"]
 
         elif "SPECRESP MATRIX" in extnames:
             h = hdulist["SPECRESP MATRIX"]
 
+        else:
+            print("Cannot find common FITS file extension for response matrix values")
+            print("Please set the `extension` keyword")
+            return
+
         data = h.data
         hdr = h.header
         hdulist.close()

From cdf337340b76068a7b90ba81add91c6efbbbdb60 Mon Sep 17 00:00:00 2001
From: eblur <lia@astro.wisc.edu>
Date: Fri, 1 Jun 2018 11:15:09 -0500
Subject: [PATCH 25/89] corrected typos

---
 specutils/spectra/xrayspectrum1d.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 68460d8b6..e2eb6ce9c 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -202,8 +202,8 @@ def _load_rmf(self, filename, extension=None):
             The file name with the RMF file
 
         extension : str (default None)
-            FITS file extinction keyword, if the response matrix is stored
-            under an extension other tham "MATRIX" or "SPECRESP MATRIX"
+            FITS file extension keyword, if the response matrix is stored
+            under an extension other than "MATRIX" or "SPECRESP MATRIX"
 
         Attributes
         ----------
@@ -466,7 +466,7 @@ def _load_arf(self, filename):
         Parameters
         ----------
         filename : str
-            The file name with the RMF file
+            The file name with the ARF file
 
         Attributes
         ----------

From 2e4724cc88652ee19ef19b450b7aa74152d4dbb3 Mon Sep 17 00:00:00 2001
From: eblur <lia@astro.wisc.edu>
Date: Fri, 1 Jun 2018 11:19:14 -0500
Subject: [PATCH 26/89] Removing superfluous returns and pass

---
 specutils/spectra/xrayspectrum1d.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index e2eb6ce9c..610dea245 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -91,7 +91,6 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
         self.exposure = exposure
         self.assign_rmf(rmf)
         self.assign_arf(arf)
-        return
 
     # Convenience function for Xray people
     @property
@@ -115,7 +114,6 @@ def assign_arf(self, arf_inp):
             self.arf = ARF(arf_inp)
         else:
             self.arf = arf_inp
-        return
 
     def assign_rmf(self, rmf_inp):
         """
@@ -190,7 +188,6 @@ def apply_response(self, mflux, exposure=None):
 class RMF(object):
     def __init__(self, filename):
         self._load_rmf(filename)
-        pass
 
     def _load_rmf(self, filename, extension=None):
         """

From a8a82450189dc3a4fb2863740e40797e58e0ddac Mon Sep 17 00:00:00 2001
From: eblur <lia@astro.wisc.edu>
Date: Fri, 1 Jun 2018 11:22:03 -0500
Subject: [PATCH 27/89] Catch ValueError

---
 specutils/spectra/xrayspectrum1d.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 610dea245..34b34e7d8 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -79,7 +79,7 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
                  arf=None, rmf=None, rest_value=0.0 * u.angstrom, **kwargs):
         try:
             axis_unit = u.Unit(bin_unit)
-        except:
+        except ValueError:
             axis_unit = _unit_parser(bin_unit)
 
         bin_mid = 0.5 * (bin_lo + bin_hi) * axis_unit

From 5d4d52c7453759d762544389f1239f582d0627d5 Mon Sep 17 00:00:00 2001
From: eblur <lia@astro.wisc.edu>
Date: Fri, 1 Jun 2018 11:27:25 -0500
Subject: [PATCH 28/89] Spelled out ARF and RMF class names

---
 specutils/spectra/xrayspectrum1d.py | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 34b34e7d8..3c452c7d3 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -5,7 +5,7 @@
 from astropy.io import fits
 from .spectrum1d import Spectrum1D
 
-__all__ = ['XraySpectrum1D', 'ARF', 'RMF']
+__all__ = ['XraySpectrum1D', 'AreaResponse', 'ResponseMatrix']
 
 # For dealing with varied unit string choices
 EV   = ['eV', 'ev']
@@ -47,12 +47,12 @@ class XraySpectrum1D(Spectrum1D):
     exposure : float
         Exposure time for the dataset
 
-    arf : specutils.ARF or string, default None
-        Strings will be passed to ARF.__init__
+    arf : specutils.AreaResponse or string, default None
+        Strings will be passed to AreaResponse.__init__
         All other input types are stored as arf attribute
 
-    rmf : specutils.RMF or string, default None
-        Strings will be passed to RMF.__init__
+    rmf : specutils.ResponseMatrix or string, default None
+        Strings will be passed to ResponseMatrix.__init__
         All other input types are stored as rmf attribute
 
     rest_value : astropy.units.Quantity, default 0 Angstrom
@@ -111,7 +111,7 @@ def assign_arf(self, arf_inp):
         Modifies the XraySpectrum1D.arf attribute
         """
         if isinstance(arf_inp, str):
-            self.arf = ARF(arf_inp)
+            self.arf = AreaResponse(arf_inp)
         else:
             self.arf = arf_inp
 
@@ -129,7 +129,7 @@ def assign_rmf(self, rmf_inp):
         Modifies the XraySpectrum1D.rmf attribute
         """
         if isinstance(rmf_inp, str):
-            self.rmf = RMF(rmf_inp)
+            self.rmf = ResponseMatrix(rmf_inp)
         else:
             self.rmf = rmf_inp
         return
@@ -185,7 +185,7 @@ def apply_response(self, mflux, exposure=None):
 
 ## ----  Supporting response file objects
 
-class RMF(object):
+class ResponseMatrix(object):
     def __init__(self, filename):
         self._load_rmf(filename)
 
@@ -449,7 +449,7 @@ def apply_rmf(self, spec):
         return counts[:self.detchans]
 
 
-class ARF(object):
+class AreaResponse(object):
 
     def __init__(self, filename):
 

From a297c452cdee7f22b15ebc52585a7712f95b29f9 Mon Sep 17 00:00:00 2001
From: eblur <lia@astro.wisc.edu>
Date: Fri, 1 Jun 2018 11:30:06 -0500
Subject: [PATCH 29/89] updated docstring

---
 specutils/spectra/xrayspectrum1d.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 3c452c7d3..cb99d21f8 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -137,8 +137,8 @@ def assign_rmf(self, rmf_inp):
     def apply_response(self, mflux, exposure=None):
         """
         Given a model flux spectrum, apply the response. In cases where the
-        spectrum has both an ARF and an RMF, apply both. Otherwise, apply
-        whatever response is in RMF.
+        spectrum has both an area response file (ARF) and a response matrix
+        file (RMF), apply both. Otherwise, apply whatever response is in RMF.
 
         The model flux spectrum *must* be created using the same units and
         bins as in the ARF (where the ARF exists)!

From 86f4a3b73772a05674f5fd7a2ed156aa018bab9b Mon Sep 17 00:00:00 2001
From: eblur <lia@astro.wisc.edu>
Date: Fri, 1 Jun 2018 11:34:47 -0500
Subject: [PATCH 30/89] Updated tests

---
 specutils/tests/test_xrayspectrum1d.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/specutils/tests/test_xrayspectrum1d.py b/specutils/tests/test_xrayspectrum1d.py
index 53fbae5c0..1e1f45bbe 100644
--- a/specutils/tests/test_xrayspectrum1d.py
+++ b/specutils/tests/test_xrayspectrum1d.py
@@ -29,6 +29,6 @@ def test_apply_model():
     test_spec = test_create_from_arrays()
     new_model = PowerLaw1D(amplitude=5.e-3, alpha=1.8, x_0=1.e3)
     model_flux = new_model(test_spec.spectral_axis.value)
-    ymodel = test_spec.apply_resp(model_flux)
+    ymodel = test_spec.apply_response(model_flux)
     assert len(ymodel) == len(test_spec.spectral_axis)
     return

From 769c923b6ba83719db4cf79f3443801e5b0384d9 Mon Sep 17 00:00:00 2001
From: Lia Corrales <lia@astro.wisc.edu>
Date: Tue, 12 Jun 2018 20:38:39 +0200
Subject: [PATCH 31/89] initialize attributes to None before loading from file

---
 specutils/spectra/xrayspectrum1d.py | 31 +++++++++++++++++++++--------
 1 file changed, 23 insertions(+), 8 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index cb99d21f8..5d98c79bc 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -187,9 +187,6 @@ def apply_response(self, mflux, exposure=None):
 
 class ResponseMatrix(object):
     def __init__(self, filename):
-        self._load_rmf(filename)
-
-    def _load_rmf(self, filename, extension=None):
         """
         Load a response matrix file (RMF) from a FITS file.
 
@@ -237,6 +234,19 @@ def _load_rmf(self, filename, extension=None):
             The number of channels in the detector
 
         """
+        self.filename = ""
+        self.offset = None
+        self.n_grp = None
+        self.f_chan = None
+        self.n_chan = None
+        self.matrix = None
+        self.energ_lo = None
+        self.energ_hi = None
+        self.energ_unit = None
+        self.detchans = None
+        self._load_rmf(filename)
+
+    def _load_rmf(self, filename, extension=None):
         # open the FITS file and extract the MATRIX extension
         # which contains the redistribution matrix and
         # anxillary information
@@ -452,11 +462,6 @@ def apply_rmf(self, spec):
 class AreaResponse(object):
 
     def __init__(self, filename):
-
-        self._load_arf(filename)
-        pass
-
-    def _load_arf(self, filename):
         """
         Load an area response file (ARF) from a FITS file.
 
@@ -492,6 +497,16 @@ def _load_arf(self, filename):
             Average exposure time for the dataset
             (takes telescope dithering into account)
         """
+        self.filename = ""
+        self.e_low = None
+        self.e_high = None
+        self.e_unit = None
+        self.specresp = None
+        self.fracexpo = None
+        self.exposure = None
+        self._load_arf(filename)
+
+    def _load_arf(self, filename):
         # open the FITS file and extract the MATRIX extension
         # which contains the redistribution matrix and
         # anxillary information

From a89af6bc769d05fabd86350a0a12ac3518adedc3 Mon Sep 17 00:00:00 2001
From: Lia Corrales <lia@astro.wisc.edu>
Date: Tue, 12 Jun 2018 20:39:55 +0200
Subject: [PATCH 32/89] save filename attribute for attempting load function

---
 specutils/spectra/xrayspectrum1d.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 5d98c79bc..58dcc3b73 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -234,7 +234,7 @@ def __init__(self, filename):
             The number of channels in the detector
 
         """
-        self.filename = ""
+        self.filename = filename
         self.offset = None
         self.n_grp = None
         self.f_chan = None
@@ -251,7 +251,6 @@ def _load_rmf(self, filename, extension=None):
         # which contains the redistribution matrix and
         # anxillary information
         hdulist = fits.open(filename)
-        self.filename = filename
 
         # get all the extension names
         extnames = np.array([h.name for h in hdulist])
@@ -497,7 +496,7 @@ def __init__(self, filename):
             Average exposure time for the dataset
             (takes telescope dithering into account)
         """
-        self.filename = ""
+        self.filename = filename
         self.e_low = None
         self.e_high = None
         self.e_unit = None
@@ -511,7 +510,6 @@ def _load_arf(self, filename):
         # which contains the redistribution matrix and
         # anxillary information
         hdulist = fits.open(filename)
-        self.filename = filename
 
         h = hdulist["SPECRESP"]
         data = h.data

From dfb9475112ea70dde9444eac8516fc5efaee1dea Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Tue, 1 May 2018 17:26:21 -0400
Subject: [PATCH 33/89] First pass at XraySpectrum1D

---
 specutils/spectra/xrayspectrum1d.py | 47 +++++++++++++++++++++++++++++
 1 file changed, 47 insertions(+)
 create mode 100644 specutils/spectra/xrayspectrum1d.py

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
new file mode 100644
index 000000000..b7556ad77
--- /dev/null
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -0,0 +1,47 @@
+import logging
+
+import numpy as np
+from astropy import units as u
+from .spectrum1d import Spectrum1D
+
+__all__ = ['XraySpectrum1D']
+
+# For dealing with varied unit string choices
+EV   = ['eV', 'ev']
+KEV  = ['kev', 'keV']
+ANGS = ['angs', 'Angs', 'Angstrom', 'angstrom', 'Angstroms', 'angstroms', 'A', 'a']
+
+def _unit_parser(unit_string):
+    if unit_string in EV:
+        return u.eV
+    if unit_string in KEV:
+        return u.keV
+    if unit_string in ANGS:
+        return u.angstrom
+
+class XraySpectrum1D(Spectrum1D):
+    """
+    Spectrum container for holding X-ray spectrum
+    WIP by eblur
+
+    Parameters
+    ----------
+    bin_lo
+    bin_hi
+    bin_unit
+    counts
+    arf
+    rmf
+    """
+    def __init__(self, bin_lo, bin_hi, bin_unit, counts, arf=None, rmf=None):
+        try:
+            axis_unit = u.Unit(bin_unit)
+        else:
+            axis_unit = _unit_parser(bin_unit)
+
+        bin_mid = 0.5 * (bin_lo + bin_hi) * unit
+        Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts)
+
+        self.arf = arf
+        self.rmf = rmf
+        return

From 2fa21425147157573e5c416de8d3a9678036be08 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Tue, 1 May 2018 17:40:28 -0400
Subject: [PATCH 34/89] Can load an X-ray spectrum

---
 specutils/spectra/__init__.py       | 4 ++++
 specutils/spectra/xrayspectrum1d.py | 4 ++--
 2 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/__init__.py b/specutils/spectra/__init__.py
index 8ce324ecb..05d5d14c4 100644
--- a/specutils/spectra/__init__.py
+++ b/specutils/spectra/__init__.py
@@ -3,6 +3,10 @@
 `~astropy.nddata.NDData`-inherited classes used for storing the spectrum data.
 """
 from .spectrum1d import *  # noqa
+<<<<<<< HEAD
 from .spectral_region import *  # noqa
 from .spectrum_collection import *  #noqa
 from .spectrum_list import * #noqa
+=======
+from .xrayspectrum1d import *
+>>>>>>> Can load an X-ray spectrum
diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index b7556ad77..d57146f97 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -36,10 +36,10 @@ class XraySpectrum1D(Spectrum1D):
     def __init__(self, bin_lo, bin_hi, bin_unit, counts, arf=None, rmf=None):
         try:
             axis_unit = u.Unit(bin_unit)
-        else:
+        except:
             axis_unit = _unit_parser(bin_unit)
 
-        bin_mid = 0.5 * (bin_lo + bin_hi) * unit
+        bin_mid = 0.5 * (bin_lo + bin_hi) * axis_unit
         Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts)
 
         self.arf = arf

From 0f37ccd0668398c6f31ae73f32fbd13b110a673a Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Tue, 1 May 2018 17:44:41 -0400
Subject: [PATCH 35/89] Added  attribute

---
 specutils/spectra/xrayspectrum1d.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index d57146f97..b939b40d1 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -45,3 +45,8 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, arf=None, rmf=None):
         self.arf = arf
         self.rmf = rmf
         return
+
+    # Convenience function for Xray people
+    @property
+    def counts(self):
+        return self.flux

From 56612ca4badb7127f5a1186aaa771c849d132797 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 13:53:39 -0400
Subject: [PATCH 36/89] Added ARF and RMF objects

---
 specutils/spectra/xrayspectrum1d.py | 326 +++++++++++++++++++++++++++-
 1 file changed, 325 insertions(+), 1 deletion(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index b939b40d1..e36a500a6 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -4,7 +4,7 @@
 from astropy import units as u
 from .spectrum1d import Spectrum1D
 
-__all__ = ['XraySpectrum1D']
+__all__ = ['XraySpectrum1D', 'ARF', 'RMF']
 
 # For dealing with varied unit string choices
 EV   = ['eV', 'ev']
@@ -50,3 +50,327 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, arf=None, rmf=None):
     @property
     def counts(self):
         return self.flux
+
+## ----  Supporting response file objects
+
+class RMF(object):
+    def __init__(self, filename):
+        self._load_rmf(filename)
+        pass
+
+    def _load_rmf(self, filename):
+        """
+        Load an RMF from a FITS file.
+
+        Parameters
+        ----------
+        filename : str
+            The file name with the RMF file
+
+        Attributes
+        ----------
+        n_grp : numpy.ndarray
+            the Array with the number of channels in each
+            channel set
+
+        f_chan : numpy.ndarray
+            The starting channel for each channel group;
+            If an element i in n_grp > 1, then the resulting
+            row entry in f_chan will be a list of length n_grp[i];
+            otherwise it will be a single number
+
+        n_chan : numpy.ndarray
+            The number of channels in each channel group. The same
+            logic as for f_chan applies
+
+        matrix : numpy.ndarray
+            The redistribution matrix as a flattened 1D vector
+
+        energ_lo : numpy.ndarray
+            The lower edges of the energy bins
+
+        energ_hi : numpy.ndarray
+            The upper edges of the energy bins
+
+        detchans : int
+            The number of channels in the detector
+
+        """
+        # open the FITS file and extract the MATRIX extension
+        # which contains the redistribution matrix and
+        # anxillary information
+        hdulist = fits.open(filename)
+
+        # get all the extension names
+        extnames = np.array([h.name for h in hdulist])
+
+        # figure out the right extension to use
+        if "MATRIX" in extnames:
+            h = hdulist["MATRIX"]
+
+        elif "SPECRESP MATRIX" in extnames:
+            h = hdulist["SPECRESP MATRIX"]
+
+        data = h.data
+        hdr = h.header
+        hdulist.close()
+
+        # extract + store the attributes described in the docstring
+        n_grp = np.array(data.field("N_GRP"))
+        f_chan = np.array(data.field('F_CHAN'))
+        n_chan = np.array(data.field("N_CHAN"))
+        matrix = np.array(data.field("MATRIX"))
+
+        self.energ_lo = np.array(data.field("ENERG_LO"))
+        self.energ_hi = np.array(data.field("ENERG_HI"))
+        self.energ_unit = data.columns["ENERG_LO"].unit
+        self.detchans = hdr["DETCHANS"]
+        self.offset = self.__get_tlmin(h)
+
+        # flatten the variable-length arrays
+        self.n_grp, self.f_chan, self.n_chan, self.matrix = \
+                self._flatten_arrays(n_grp, f_chan, n_chan, matrix)
+
+        return
+
+    def __get_tlmin(self, h):
+        """
+        Get the tlmin keyword for `F_CHAN`.
+
+        Parameters
+        ----------
+        h : an astropy.io.fits.hdu.table.BinTableHDU object
+            The extension containing the `F_CHAN` column
+
+        Returns
+        -------
+        tlmin : int
+            The tlmin keyword
+        """
+        # get the header
+        hdr = h.header
+        # get the keys of all
+        keys = np.array(list(hdr.keys()))
+
+        # find the place where the tlmin keyword is defined
+        t = np.array(["TLMIN" in k for k in keys])
+
+        # get the index of the TLMIN keyword
+        tlmin_idx = np.hstack(np.where(t))[0]
+
+        # get the corresponding value
+        tlmin = np.int(list(hdr.items())[tlmin_idx][1])
+
+        return tlmin
+
+    def _flatten_arrays(self, n_grp, f_chan, n_chan, matrix):
+
+        if not len(n_grp) == len(f_chan) == len(n_chan) == len(matrix):
+            raise ValueError("Arrays must be of same length!")
+
+        # find all non-zero groups
+        nz_idx = (n_grp > 0)
+
+        # stack all non-zero rows in the matrix
+        matrix_flat = np.hstack(matrix[nz_idx])
+
+        # stack all nonzero rows in n_chan and f_chan
+        #n_chan_flat = np.hstack(n_chan[nz_idx])
+        #f_chan_flat = np.hstack(f_chan[nz_idx])
+
+        # some matrices actually have more elements
+        # than groups in `n_grp`, so we'll only pick out
+        # those values that have a correspondence in
+        # n_grp
+        f_chan_new = []
+        n_chan_new = []
+        for i,t in enumerate(nz_idx):
+            if t:
+                n = n_grp[i]
+                f = f_chan[i]
+                nc = n_chan[i]
+                if np.size(f) == 1:
+                    f_chan_new.append(f)
+                    n_chan_new.append(nc)
+                else:
+                    f_chan_new.append(f[:n])
+                    n_chan_new.append(nc[:n])
+
+        n_chan_flat = np.hstack(n_chan_new)
+        f_chan_flat = np.hstack(f_chan_new)
+
+        # if n_chan is zero, we'll remove those as well.
+        nz_idx2 = (n_chan_flat > 0)
+        n_chan_flat = n_chan_flat[nz_idx2]
+        f_chan_flat = f_chan_flat[nz_idx2]
+
+        return n_grp, f_chan_flat, n_chan_flat, matrix_flat
+
+    def apply_rmf(self, spec):
+        """
+        Fold the spectrum through the redistribution matrix.
+
+        The redistribution matrix is saved as a flattened 1-dimensional
+        vector to save space. In reality, for each entry in the flux
+        vector, there exists one or more sets of channels that this
+        flux is redistributed into. The additional arrays `n_grp`,
+        `f_chan` and `n_chan` store this information:
+            * `n_group` stores the number of channel groups for each
+              energy bin
+            * `f_chan` stores the *first channel* that each channel
+              for each channel set
+            * `n_chan` stores the number of channels in each channel
+              set
+
+        As a result, for a given energy bin i, we need to look up the
+        number of channel sets in `n_grp` for that energy bin. We
+        then need to loop over the number of channel sets. For each
+        channel set, we look up the first channel into which flux
+        will be distributed as well as the number of channels in the
+        group. We then need to also loop over the these channels and
+        actually use the corresponding elements in the redistribution
+        matrix to redistribute the photon flux into channels.
+
+        All of this is basically a big bookkeeping exercise in making
+        sure to get the indices right.
+
+        Parameters
+        ----------
+        spec : numpy.ndarray
+            The (model) spectrum to be folded
+
+        Returns
+        -------
+        counts : numpy.ndarray
+            The (model) spectrum after folding, in
+            counts/s/channel
+
+        """
+        # get the number of channels in the data
+        nchannels = spec.shape[0]
+
+        # an empty array for the output counts
+        counts = np.zeros(nchannels)
+
+        # index for n_chan and f_chan incrementation
+        k = 0
+
+        # index for the response matrix incrementation
+        resp_idx = 0
+
+        # loop over all channels
+        for i in range(nchannels):
+
+            # this is the current bin in the flux spectrum to
+            # be folded
+            source_bin_i = spec[i]
+
+            # get the current number of groups
+            current_num_groups = self.n_grp[i]
+
+            # loop over the current number of groups
+            for j in range(current_num_groups):
+
+                current_num_chans = int(self.n_chan[k])
+
+                if current_num_chans == 0:
+                    k += 1
+                    resp_idx += current_num_chans
+                    continue
+
+
+                else:
+                    # get the right index for the start of the counts array
+                    # to put the data into
+                    counts_idx = int(self.f_chan[k] - self.offset)
+                    # this is the current number of channels to use
+
+                    k += 1
+                    # add the flux to the subarray of the counts array that starts with
+                    # counts_idx and runs over current_num_chans channels
+                    counts[counts_idx:counts_idx +
+                                      current_num_chans] += self.matrix[resp_idx:resp_idx +
+                                                                                 current_num_chans] * \
+                                                                np.float(source_bin_i)
+                    # iterate the response index for next round
+                    resp_idx += current_num_chans
+
+
+        return counts[:self.detchans]
+
+
+class ARF(object):
+
+    def __init__(self, filename):
+
+        self._load_arf(filename)
+        pass
+
+    def _load_arf(self, filename):
+        """
+        Load an ARF from a FITS file.
+        Parameters
+        ----------
+        filename : str
+            The file name with the RMF file
+        Attributes
+        ----------
+        """
+        # open the FITS file and extract the MATRIX extension
+        # which contains the redistribution matrix and
+        # anxillary information
+        hdulist = fits.open(filename)
+        h = hdulist["SPECRESP"]
+        data = h.data
+        hdr = h.header
+        hdulist.close()
+
+        # extract + store the attributes described in the docstring
+
+        self.e_low  = np.array(data.field("ENERG_LO"))
+        self.e_high = np.array(data.field("ENERG_HI"))
+        self.e_unit = data.columns["ENERG_LO"].unit
+        self.specresp = np.array(data.field("SPECRESP"))
+
+        if "EXPOSURE" in list(hdr.keys()):
+            self.exposure = hdr["EXPOSURE"]
+        else:
+            self.exposure = 1.0
+
+        if "FRACEXPO" in data.columns.names:
+            self.fracexpo = data["FRACEXPO"]
+        else:
+            self.fracexpo = 1.0
+
+        return
+
+    def apply_arf(self, spec, exposure=None):
+        """
+        Fold the spectrum through the ARF.
+        The ARF is a single vector encoding the effective area information
+        about the detector. A such, applying the ARF is a simple
+        multiplication with the input spectrum.
+        Parameters
+        ----------
+        spec : numpy.ndarray
+            The (model) spectrum to be folded
+        exposure : float, default None
+            Value for the exposure time. By default, `apply_arf` will use the
+            exposure keyword from the ARF file. If this exposure time is not
+            correct (for example when simulated spectra use a different exposure
+            time and the ARF from a real observation), one can override the
+            default exposure by setting the `exposure` keyword to the correct
+            value.
+        Returns
+        -------
+        s_arf : numpy.ndarray
+            The (model) spectrum after folding, in
+            counts/s/channel
+        """
+        assert spec.shape[0] == self.specresp.shape[0], "The input spectrum must " \
+                                                      "be of same size as the " \
+                                                      "ARF array."
+        if exposure is None:
+            return np.array(spec) * self.specresp * self.exposure
+        else:
+            return np.array(spec) * self.specresp * exposure

From b7303e75966516657bb8ae88a5e1b5f4dd3cda2b Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 14:07:50 -0400
Subject: [PATCH 37/89] Added assign_rmf and assign_arf methods

---
 specutils/spectra/xrayspectrum1d.py | 29 +++++++++++++++++++++++++++--
 1 file changed, 27 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index e36a500a6..a371b94a5 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -42,8 +42,8 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, arf=None, rmf=None):
         bin_mid = 0.5 * (bin_lo + bin_hi) * axis_unit
         Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts)
 
-        self.arf = arf
-        self.rmf = rmf
+        self.assign_rmf(rmf)
+        self.assign_arf(rmf)
         return
 
     # Convenience function for Xray people
@@ -51,6 +51,21 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, arf=None, rmf=None):
     def counts(self):
         return self.flux
 
+    def assign_arf(self, arf_inp):
+        if isinstance(arf_inp, str):
+            self.arf = ARF(arf_inp)
+        else:
+            self.arf = arf_inp
+        return
+
+    def assign_rmf(self, rmf_inp):
+        if isinstance(rmf_inp, str):
+            self.rmf = RMF(rmf_inp)
+        else:
+            self.rmf = rmf_inp
+        return
+
+
 ## ----  Supporting response file objects
 
 class RMF(object):
@@ -69,6 +84,9 @@ def _load_rmf(self, filename):
 
         Attributes
         ----------
+        filename : str
+            The file name that the RMF was drawn from
+
         n_grp : numpy.ndarray
             the Array with the number of channels in each
             channel set
@@ -100,6 +118,7 @@ def _load_rmf(self, filename):
         # which contains the redistribution matrix and
         # anxillary information
         hdulist = fits.open(filename)
+        self.filename = filename
 
         # get all the extension names
         extnames = np.array([h.name for h in hdulist])
@@ -309,17 +328,23 @@ def __init__(self, filename):
     def _load_arf(self, filename):
         """
         Load an ARF from a FITS file.
+
         Parameters
         ----------
         filename : str
             The file name with the RMF file
+
         Attributes
         ----------
+        filename : str
+            The file name that the ARF was drawn from
         """
         # open the FITS file and extract the MATRIX extension
         # which contains the redistribution matrix and
         # anxillary information
         hdulist = fits.open(filename)
+        self.filename = filename
+        
         h = hdulist["SPECRESP"]
         data = h.data
         hdr = h.header

From 003b63696503f64d17f44671bc6995fc3b9b9224 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 14:16:13 -0400
Subject: [PATCH 38/89] Added exposure input

---
 specutils/spectra/xrayspectrum1d.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index a371b94a5..51de9e726 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -33,7 +33,7 @@ class XraySpectrum1D(Spectrum1D):
     arf
     rmf
     """
-    def __init__(self, bin_lo, bin_hi, bin_unit, counts, arf=None, rmf=None):
+    def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure, arf=None, rmf=None):
         try:
             axis_unit = u.Unit(bin_unit)
         except:
@@ -42,6 +42,7 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, arf=None, rmf=None):
         bin_mid = 0.5 * (bin_lo + bin_hi) * axis_unit
         Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts)
 
+        self.exposure = exposure
         self.assign_rmf(rmf)
         self.assign_arf(rmf)
         return
@@ -344,7 +345,7 @@ def _load_arf(self, filename):
         # anxillary information
         hdulist = fits.open(filename)
         self.filename = filename
-        
+
         h = hdulist["SPECRESP"]
         data = h.data
         hdr = h.header

From ce4349e4e76c20c69bbedc0fab8a154fbc165796 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 14:18:22 -0400
Subject: [PATCH 39/89] Fixed bugs

---
 specutils/spectra/xrayspectrum1d.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 51de9e726..fade12358 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -2,6 +2,7 @@
 
 import numpy as np
 from astropy import units as u
+from astropy.io import fits
 from .spectrum1d import Spectrum1D
 
 __all__ = ['XraySpectrum1D', 'ARF', 'RMF']
@@ -44,7 +45,7 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure, arf=None, rmf=Non
 
         self.exposure = exposure
         self.assign_rmf(rmf)
-        self.assign_arf(rmf)
+        self.assign_arf(arf)
         return
 
     # Convenience function for Xray people

From 753ab5eb08f712fa53ab2f78a41d941fcbd3da32 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 14:30:44 -0400
Subject: [PATCH 40/89] Updated docstrings

---
 specutils/spectra/xrayspectrum1d.py | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index fade12358..72c5bd11e 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -112,6 +112,9 @@ def _load_rmf(self, filename):
         energ_hi : numpy.ndarray
             The upper edges of the energy bins
 
+        energ_unit : str
+            Description of the energy units used
+
         detchans : int
             The number of channels in the detector
 
@@ -340,6 +343,28 @@ def _load_arf(self, filename):
         ----------
         filename : str
             The file name that the ARF was drawn from
+
+        e_low : numpy.ndarray
+            The lower edges of the energy bins
+
+        e_high : numpy.ndarray
+            The upper edges of the energy bins
+
+        e_unit : str
+            Description of the energy units used
+
+        specresp : numpy.ndarray
+            Description of the energy dependent telescope response area
+
+        fracexpo : float or numpy.ndarray
+            Fractional exposure time for the spectrum
+            (sometimes constant, sometimes dependent on spectral channel).
+            These values are stored for reference; generally, they are already
+            accounted for in the specresp array.
+
+        exposure :
+            Average exposure time for the dataset
+            (takes telescope dithering into account)
         """
         # open the FITS file and extract the MATRIX extension
         # which contains the redistribution matrix and

From 4f195e1be723c2b9e833156ac47de979fdd693f8 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 14:36:08 -0400
Subject: [PATCH 41/89] Store bin and energy units as astropy.units.Unit

---
 specutils/spectra/xrayspectrum1d.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 72c5bd11e..2b98164aa 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -31,6 +31,7 @@ class XraySpectrum1D(Spectrum1D):
     bin_hi
     bin_unit
     counts
+    exposure
     arf
     rmf
     """
@@ -43,6 +44,7 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure, arf=None, rmf=Non
         bin_mid = 0.5 * (bin_lo + bin_hi) * axis_unit
         Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts)
 
+        self.bin_unit = axis_unit # keep for reference
         self.exposure = exposure
         self.assign_rmf(rmf)
         self.assign_arf(arf)
@@ -112,7 +114,7 @@ def _load_rmf(self, filename):
         energ_hi : numpy.ndarray
             The upper edges of the energy bins
 
-        energ_unit : str
+        energ_unit : astropy.units.Unit
             Description of the energy units used
 
         detchans : int
@@ -147,7 +149,7 @@ def _load_rmf(self, filename):
 
         self.energ_lo = np.array(data.field("ENERG_LO"))
         self.energ_hi = np.array(data.field("ENERG_HI"))
-        self.energ_unit = data.columns["ENERG_LO"].unit
+        self.energ_unit = _unit_parser(data.columns["ENERG_LO"].unit)
         self.detchans = hdr["DETCHANS"]
         self.offset = self.__get_tlmin(h)
 
@@ -350,7 +352,7 @@ def _load_arf(self, filename):
         e_high : numpy.ndarray
             The upper edges of the energy bins
 
-        e_unit : str
+        e_unit : astropy.units.Unit
             Description of the energy units used
 
         specresp : numpy.ndarray
@@ -381,7 +383,7 @@ def _load_arf(self, filename):
 
         self.e_low  = np.array(data.field("ENERG_LO"))
         self.e_high = np.array(data.field("ENERG_HI"))
-        self.e_unit = data.columns["ENERG_LO"].unit
+        self.e_unit = _unit_parser(data.columns["ENERG_LO"].unit)
         self.specresp = np.array(data.field("SPECRESP"))
 
         if "EXPOSURE" in list(hdr.keys()):

From 67493c523e1f30624db50cca1e204319fc1fa942 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 16:31:14 -0400
Subject: [PATCH 42/89] Added bin_lo and hi because wcs.bin_edges did not exist

---
 specutils/spectra/xrayspectrum1d.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 2b98164aa..d209c988b 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -35,15 +35,18 @@ class XraySpectrum1D(Spectrum1D):
     arf
     rmf
     """
-    def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure, arf=None, rmf=None):
+    def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
+                 arf=None, rmf=None, **kwargs):
         try:
             axis_unit = u.Unit(bin_unit)
         except:
             axis_unit = _unit_parser(bin_unit)
 
         bin_mid = 0.5 * (bin_lo + bin_hi) * axis_unit
-        Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts)
+        Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts, **kwargs)
 
+        self.bin_lo = bin_lo
+        self.bin_hi = bin_hi
         self.bin_unit = axis_unit # keep for reference
         self.exposure = exposure
         self.assign_rmf(rmf)

From 7974939f7a92a44b437aad94a8c31b202e59e14f Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 16:34:35 -0400
Subject: [PATCH 43/89] Added apply_resp method

---
 specutils/spectra/xrayspectrum1d.py | 44 +++++++++++++++++++++++++++++
 1 file changed, 44 insertions(+)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index d209c988b..765e77241 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -72,6 +72,50 @@ def assign_rmf(self, rmf_inp):
             self.rmf = rmf_inp
         return
 
+    def apply_resp(self, mflux, exposure=None):
+        """
+        Given a model flux spectrum, apply the response. In cases where the
+        spectrum has both an ARF and an RMF, apply both. Otherwise, apply
+        whatever response is in RMF.
+
+        The model flux spectrum *must* be created using the same units and
+        bins as in the ARF (where the ARF exists)!
+
+        Parameters
+        ----------
+        mflux : iterable
+            A list or array with the model flux values in ergs/keV/s/cm^-2
+
+        exposure : float, default None
+            By default, the exposure stored in the ARF will be used to compute
+            the total counts per bin over the effective observation time.
+            In cases where this might be incorrect (e.g. for simulated spectra
+            where the pha file might have a different exposure value than the
+            ARF), this keyword provides the functionality to override the
+            default behaviour and manually set the exposure time to use.
+
+        Returns
+        -------
+        count_model : numpy.ndarray
+            The model spectrum in units of counts/bin
+
+        If no ARF file exists, it will return the model flux after applying the RMF
+        If no RMF file exists, it will return the model flux after applying the ARF (with a warning)
+        If no ARF and no RMF, it will return the model flux spectrum (with a warning)
+        """
+
+        if self.arf is not None:
+            mrate  = self.arf.apply_arf(mflux, exposure=exposure)
+        else:
+            mrate = mflux
+
+        if self.rmf is not None:
+            count_model = self.rmf.apply_rmf(mrate)
+            return count_model
+        else:
+            print("Caution: no response file specified")
+            return mrate
+
 
 ## ----  Supporting response file objects
 

From 7f54eaeaca2e3d132d76da83789986fce76abf7e Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 16:52:47 -0400
Subject: [PATCH 44/89] Adjust apply_resp output if spectrum is stored by
 wavelength

---
 specutils/spectra/xrayspectrum1d.py | 13 ++++++++++---
 1 file changed, 10 insertions(+), 3 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 765e77241..b9f2df62d 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -110,11 +110,18 @@ def apply_resp(self, mflux, exposure=None):
             mrate = mflux
 
         if self.rmf is not None:
-            count_model = self.rmf.apply_rmf(mrate)
-            return count_model
+            result = self.rmf.apply_rmf(mrate)
         else:
             print("Caution: no response file specified")
-            return mrate
+            result = mrate
+
+        # The response is always stored according to energy
+        # If the spectrum is stored in wavelength, reverse the result so it
+        # matches in wavelength space
+        if self.bin_unit.physical_type == 'length':
+            return result[::-1]
+        else:
+            return result
 
 
 ## ----  Supporting response file objects

From 5eded3d6b2716a25795144687a9cf35b1ce04939 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 18:03:29 -0400
Subject: [PATCH 45/89] Don't need bin_unit

---
 specutils/spectra/xrayspectrum1d.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index b9f2df62d..bac3f6e3e 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -47,7 +47,6 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
 
         self.bin_lo = bin_lo
         self.bin_hi = bin_hi
-        self.bin_unit = axis_unit # keep for reference
         self.exposure = exposure
         self.assign_rmf(rmf)
         self.assign_arf(arf)
@@ -118,7 +117,8 @@ def apply_resp(self, mflux, exposure=None):
         # The response is always stored according to energy
         # If the spectrum is stored in wavelength, reverse the result so it
         # matches in wavelength space
-        if self.bin_unit.physical_type == 'length':
+        if self.spectral_axis.unit.physical_type == 'length':
+            print("reversing output")
             return result[::-1]
         else:
             return result

From 6705068c79c934eec26af075c8a3aab9c45490df Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 18:12:10 -0400
Subject: [PATCH 46/89] Because pha-file wavelength is in reverse order, the
 output from apply_resp is in correct order. No longer need to check

---
 specutils/spectra/xrayspectrum1d.py | 17 ++++++++---------
 1 file changed, 8 insertions(+), 9 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index bac3f6e3e..c054255df 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -34,6 +34,11 @@ class XraySpectrum1D(Spectrum1D):
     exposure
     arf
     rmf
+
+    Attributes
+    ----------
+    model : numpy.ndarray
+        Stores model values
     """
     def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
                  arf=None, rmf=None, **kwargs):
@@ -50,6 +55,7 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
         self.exposure = exposure
         self.assign_rmf(rmf)
         self.assign_arf(arf)
+        self.model_counts = np.zeros_like(bin_mid)
         return
 
     # Convenience function for Xray people
@@ -114,15 +120,8 @@ def apply_resp(self, mflux, exposure=None):
             print("Caution: no response file specified")
             result = mrate
 
-        # The response is always stored according to energy
-        # If the spectrum is stored in wavelength, reverse the result so it
-        # matches in wavelength space
-        if self.spectral_axis.unit.physical_type == 'length':
-            print("reversing output")
-            return result[::-1]
-        else:
-            return result
-
+        self.model_counts = result
+        return result
 
 ## ----  Supporting response file objects
 

From 80d559d6c62b192fada371b9be12fa6905aed65e Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 2 May 2018 18:15:20 -0400
Subject: [PATCH 47/89] added save_model_counts option

---
 specutils/spectra/xrayspectrum1d.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index c054255df..cd78d2e49 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -77,7 +77,7 @@ def assign_rmf(self, rmf_inp):
             self.rmf = rmf_inp
         return
 
-    def apply_resp(self, mflux, exposure=None):
+    def apply_resp(self, mflux, exposure=None, store_model_counts=True):
         """
         Given a model flux spectrum, apply the response. In cases where the
         spectrum has both an ARF and an RMF, apply both. Otherwise, apply
@@ -99,9 +99,12 @@ def apply_resp(self, mflux, exposure=None):
             ARF), this keyword provides the functionality to override the
             default behaviour and manually set the exposure time to use.
 
+        store_model_counts : bool
+            If True, the output will also be stored in self.model_counts
+
         Returns
         -------
-        count_model : numpy.ndarray
+        model_counts : numpy.ndarray
             The model spectrum in units of counts/bin
 
         If no ARF file exists, it will return the model flux after applying the RMF
@@ -120,7 +123,8 @@ def apply_resp(self, mflux, exposure=None):
             print("Caution: no response file specified")
             result = mrate
 
-        self.model_counts = result
+        if store_model_counts:
+            self.model_counts = result
         return result
 
 ## ----  Supporting response file objects

From 73535d9344d037824362a7e9fb7cbe679661c089 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Thu, 3 May 2018 10:06:41 -0400
Subject: [PATCH 48/89] Changed my mind on model_counts thing for now

---
 specutils/spectra/xrayspectrum1d.py | 10 +---------
 1 file changed, 1 insertion(+), 9 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index cd78d2e49..7b342cb96 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -34,11 +34,6 @@ class XraySpectrum1D(Spectrum1D):
     exposure
     arf
     rmf
-
-    Attributes
-    ----------
-    model : numpy.ndarray
-        Stores model values
     """
     def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
                  arf=None, rmf=None, **kwargs):
@@ -55,7 +50,6 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
         self.exposure = exposure
         self.assign_rmf(rmf)
         self.assign_arf(arf)
-        self.model_counts = np.zeros_like(bin_mid)
         return
 
     # Convenience function for Xray people
@@ -77,7 +71,7 @@ def assign_rmf(self, rmf_inp):
             self.rmf = rmf_inp
         return
 
-    def apply_resp(self, mflux, exposure=None, store_model_counts=True):
+    def apply_resp(self, mflux, exposure=None):
         """
         Given a model flux spectrum, apply the response. In cases where the
         spectrum has both an ARF and an RMF, apply both. Otherwise, apply
@@ -123,8 +117,6 @@ def apply_resp(self, mflux, exposure=None, store_model_counts=True):
             print("Caution: no response file specified")
             result = mrate
 
-        if store_model_counts:
-            self.model_counts = result
         return result
 
 ## ----  Supporting response file objects

From e9be1099fea9034f52ba8ba8e6e30e1b2610b682 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 23 May 2018 18:27:59 -0500
Subject: [PATCH 49/89] Hard coded in a rest_value

---
 specutils/spectra/xrayspectrum1d.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 7b342cb96..2c5f4d86e 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -36,14 +36,15 @@ class XraySpectrum1D(Spectrum1D):
     rmf
     """
     def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
-                 arf=None, rmf=None, **kwargs):
+                 arf=None, rmf=None, rest_value=0.0, **kwargs):
         try:
             axis_unit = u.Unit(bin_unit)
         except:
             axis_unit = _unit_parser(bin_unit)
 
         bin_mid = 0.5 * (bin_lo + bin_hi) * axis_unit
-        Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts, **kwargs)
+        Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts,
+                            rest_value=rest_value, **kwargs)
 
         self.bin_lo = bin_lo
         self.bin_hi = bin_hi

From baa613c99e45bd8b78412e4c0b999f5143962a07 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 23 May 2018 18:29:31 -0500
Subject: [PATCH 50/89] Added unit to rest_value

---
 specutils/spectra/xrayspectrum1d.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 2c5f4d86e..d700837b7 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -36,7 +36,7 @@ class XraySpectrum1D(Spectrum1D):
     rmf
     """
     def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
-                 arf=None, rmf=None, rest_value=0.0, **kwargs):
+                 arf=None, rmf=None, rest_value=0.0 * u.angstrom, **kwargs):
         try:
             axis_unit = u.Unit(bin_unit)
         except:

From 56b4c1afdd8f667365ac005d804a30ca27ae10d1 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 23 May 2018 18:47:34 -0500
Subject: [PATCH 51/89] Commented out distracting print statements

---
 specutils/spectra/xrayspectrum1d.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index d700837b7..ceb3f9fdf 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -110,12 +110,13 @@ def apply_resp(self, mflux, exposure=None):
         if self.arf is not None:
             mrate  = self.arf.apply_arf(mflux, exposure=exposure)
         else:
+            #print("Caution: no ARF file specified")
             mrate = mflux
 
         if self.rmf is not None:
             result = self.rmf.apply_rmf(mrate)
         else:
-            print("Caution: no response file specified")
+            #print("Caution: no RMF file specified")
             result = mrate
 
         return result

From afb6333c6337a5104eaf4684cd2be1e7aa53105f Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 23 May 2018 19:00:00 -0500
Subject: [PATCH 52/89] Wrote tests

---
 specutils/tests/test_xrayspectrum1d.py | 35 ++++++++++++++++++++++++++
 1 file changed, 35 insertions(+)
 create mode 100644 specutils/tests/test_xrayspectrum1d.py

diff --git a/specutils/tests/test_xrayspectrum1d.py b/specutils/tests/test_xrayspectrum1d.py
new file mode 100644
index 000000000..f8792a7ff
--- /dev/null
+++ b/specutils/tests/test_xrayspectrum1d.py
@@ -0,0 +1,35 @@
+import numpy as np
+from scipy import optimize
+import astropy.units as u
+from astropy.modeling.powerlaws import PowerLaw1D
+from specutils import XraySpectrum1D
+
+def test_create_from_arrays():
+    # Test that XraySpectrum1D can be initialized
+    amp0, alpha0 = 3.e-3, 2.0
+    powlaw0 = PowerLaw1D(amplitude=amp0, alpha=alpha0, x_0=1.e3)
+    energy = np.linspace(0.2, 10.0, 8000)
+    elo = energy[:-1]
+    ehi = energy[1:]
+    emid = 0.5 * (elo + ehi)
+    counts = np.random.poisson(lam=powlaw0(emid), size=len(emid))
+    test_spec = XraySpectrum1D(elo, ehi, u.keV, counts, exposure=1.0)
+    return test_spec
+
+def test_call_spectrum():
+    # Test that the XraySpectrum1D has attributes that make it unique
+    # to X-ray spectrum
+    test_spec = test_create_from_arrays()
+    assert len(test_spec.counts) == len(test_spec.spectral_axis)
+    assert len(test_spec.bin_lo) == len(test_spec.spectral_axis)
+    assert len(test_spec.bin_hi) == len(test_spec.spectral_axis)
+    return
+
+def test_apply_model():
+    # Test that one can evaluate XraySpectrum1D with a model
+    test_spec = test_create_from_arrays()
+    new_model = PowerLaw1D(amplitude=5.e-3, alpha=1.8, x_0=1.e3)
+    model_flux = new_model(test_spec.spectral_axis.value)
+    ymodel = test_spec.apply_resp(model_flux)
+    assert len(ymodel) == len(test_spec.spectral_axis)
+    return

From a15b101ec3b083b34924569bc6cd704542572ddf Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 23 May 2018 19:17:48 -0500
Subject: [PATCH 53/89] Added docstrings

---
 specutils/spectra/xrayspectrum1d.py | 67 +++++++++++++++++++++++++++--
 1 file changed, 64 insertions(+), 3 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index ceb3f9fdf..0b72ccedf 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -13,6 +13,9 @@
 ANGS = ['angs', 'Angs', 'Angstrom', 'angstrom', 'Angstroms', 'angstroms', 'A', 'a']
 
 def _unit_parser(unit_string):
+    """
+    A function for matching unit strings to the correct Astropy unit
+    """
     if unit_string in EV:
         return u.eV
     if unit_string in KEV:
@@ -22,17 +25,51 @@ def _unit_parser(unit_string):
 
 class XraySpectrum1D(Spectrum1D):
     """
-    Spectrum container for holding X-ray spectrum
-    WIP by eblur
+    Spectrum properties specific to X-ray data
 
     Parameters
     ----------
+    bin_lo : numpy.ndarray
+        The left edges for bin values
+
+    bin_hi : numpy.ndarray
+        The right edges for bin values
+
+    bin_unit : string OR astropy.units.Unit
+        Unit for the bin values
+
+    counts : numpy.ndarray
+        Counts histogram for the X-ray spectrum
+
+    exposure : float
+        Exposure time for the dataset
+
+    arf : specutils.ARF or string, default None
+        Strings will be passed to ARF.__init__
+        All other input types are stored as arf attribute
+
+    rmf : specutils.RMF or string, default None
+        Strings will be passed to RMF.__init__
+        All other input types are stored as rmf attribute
+
+    rest_value : astropy.units.Quantity, default 0 Angstrom
+        See Spectrum1D rest_value input
+
+    Attributes
+    ----------
+    Inherits from Spectrum1D object.
+    The following inputs are stored as additional attributes.
+
     bin_lo
+
     bin_hi
-    bin_unit
+
     counts
+
     exposure
+
     arf
+
     rmf
     """
     def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
@@ -59,6 +96,18 @@ def counts(self):
         return self.flux
 
     def assign_arf(self, arf_inp):
+        """
+        Assign an ARF object to the XraySpectrum1D object
+
+        Input
+        -----
+        arf_inp : string
+            File name for the ARF (FITS file)
+
+        Returns
+        -------
+        Modifies the XraySpectrum1D.arf attribute
+        """
         if isinstance(arf_inp, str):
             self.arf = ARF(arf_inp)
         else:
@@ -66,6 +115,18 @@ def assign_arf(self, arf_inp):
         return
 
     def assign_rmf(self, rmf_inp):
+        """
+        Assign an RMF object to the XraySpectrum1D object
+
+        Input
+        -----
+        rmf_inp : string
+            File name for the RMF (FITS file)
+
+        Returns
+        -------
+        Modifies the XraySpectrum1D.rmf attribute
+        """
         if isinstance(rmf_inp, str):
             self.rmf = RMF(rmf_inp)
         else:

From 4ba7b5829c400d13d2bc976d241b6fb80981881d Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 23 May 2018 19:19:32 -0500
Subject: [PATCH 54/89] Removed scipy dependence

---
 specutils/tests/test_xrayspectrum1d.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/specutils/tests/test_xrayspectrum1d.py b/specutils/tests/test_xrayspectrum1d.py
index f8792a7ff..53fbae5c0 100644
--- a/specutils/tests/test_xrayspectrum1d.py
+++ b/specutils/tests/test_xrayspectrum1d.py
@@ -1,5 +1,4 @@
 import numpy as np
-from scipy import optimize
 import astropy.units as u
 from astropy.modeling.powerlaws import PowerLaw1D
 from specutils import XraySpectrum1D

From f5a996a0175e8c36d1a12c91c6109ebb27fe1495 Mon Sep 17 00:00:00 2001
From: eblur <lia@astro.wisc.edu>
Date: Fri, 1 Jun 2018 10:56:46 -0500
Subject: [PATCH 55/89] Spelled out acronyms, link to edu materials

---
 specutils/spectra/xrayspectrum1d.py | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 0b72ccedf..05e9cc892 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -12,6 +12,9 @@
 KEV  = ['kev', 'keV']
 ANGS = ['angs', 'Angs', 'Angstrom', 'angstrom', 'Angstroms', 'angstroms', 'A', 'a']
 
+## An introduction to X-ray analysis can be found here:
+## http://cxc.cfa.harvard.edu/xrayschool/talks/intro_xray_analysis.pdf
+
 def _unit_parser(unit_string):
     """
     A function for matching unit strings to the correct Astropy unit
@@ -97,12 +100,12 @@ def counts(self):
 
     def assign_arf(self, arf_inp):
         """
-        Assign an ARF object to the XraySpectrum1D object
+        Assign an area response file (ARF) object to the XraySpectrum1D object
 
         Input
         -----
         arf_inp : string
-            File name for the ARF (FITS file)
+            File name for the area response file (FITS file)
 
         Returns
         -------
@@ -116,12 +119,12 @@ def assign_arf(self, arf_inp):
 
     def assign_rmf(self, rmf_inp):
         """
-        Assign an RMF object to the XraySpectrum1D object
+        Assign a response matrix file (RMF) object to the XraySpectrum1D object
 
         Input
         -----
         rmf_inp : string
-            File name for the RMF (FITS file)
+            File name for the response matrix file (FITS file)
 
         Returns
         -------
@@ -191,7 +194,7 @@ def __init__(self, filename):
 
     def _load_rmf(self, filename):
         """
-        Load an RMF from a FITS file.
+        Load a response matrix file (RMF) from a FITS file.
 
         Parameters
         ----------
@@ -446,7 +449,7 @@ def __init__(self, filename):
 
     def _load_arf(self, filename):
         """
-        Load an ARF from a FITS file.
+        Load an area response file (ARF) from a FITS file.
 
         Parameters
         ----------
@@ -512,7 +515,7 @@ def _load_arf(self, filename):
 
     def apply_arf(self, spec, exposure=None):
         """
-        Fold the spectrum through the ARF.
+        Fold the spectrum through the area response file (ARF).
         The ARF is a single vector encoding the effective area information
         about the detector. A such, applying the ARF is a simple
         multiplication with the input spectrum.

From eed588cdb5c6c785d059ba77d6d120ed8c6208f6 Mon Sep 17 00:00:00 2001
From: eblur <lia@astro.wisc.edu>
Date: Fri, 1 Jun 2018 11:05:22 -0500
Subject: [PATCH 56/89] Added option to specify RMF FITS extension for response
 matrix

---
 specutils/spectra/xrayspectrum1d.py | 18 +++++++++++++++---
 1 file changed, 15 insertions(+), 3 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 05e9cc892..68460d8b6 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -136,7 +136,7 @@ def assign_rmf(self, rmf_inp):
             self.rmf = rmf_inp
         return
 
-    def apply_resp(self, mflux, exposure=None):
+    def apply_response(self, mflux, exposure=None):
         """
         Given a model flux spectrum, apply the response. In cases where the
         spectrum has both an ARF and an RMF, apply both. Otherwise, apply
@@ -192,7 +192,7 @@ def __init__(self, filename):
         self._load_rmf(filename)
         pass
 
-    def _load_rmf(self, filename):
+    def _load_rmf(self, filename, extension=None):
         """
         Load a response matrix file (RMF) from a FITS file.
 
@@ -201,6 +201,10 @@ def _load_rmf(self, filename):
         filename : str
             The file name with the RMF file
 
+        extension : str (default None)
+            FITS file extinction keyword, if the response matrix is stored
+            under an extension other tham "MATRIX" or "SPECRESP MATRIX"
+
         Attributes
         ----------
         filename : str
@@ -246,12 +250,20 @@ def _load_rmf(self, filename):
         extnames = np.array([h.name for h in hdulist])
 
         # figure out the right extension to use
-        if "MATRIX" in extnames:
+        if extension is not None:
+            h = hdulist[extension]
+
+        elif "MATRIX" in extnames:
             h = hdulist["MATRIX"]
 
         elif "SPECRESP MATRIX" in extnames:
             h = hdulist["SPECRESP MATRIX"]
 
+        else:
+            print("Cannot find common FITS file extension for response matrix values")
+            print("Please set the `extension` keyword")
+            return
+
         data = h.data
         hdr = h.header
         hdulist.close()

From 0fa3718a048314f120c9d2bb63a88af574db2496 Mon Sep 17 00:00:00 2001
From: eblur <lia@astro.wisc.edu>
Date: Fri, 1 Jun 2018 11:15:09 -0500
Subject: [PATCH 57/89] corrected typos

---
 specutils/spectra/xrayspectrum1d.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 68460d8b6..e2eb6ce9c 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -202,8 +202,8 @@ def _load_rmf(self, filename, extension=None):
             The file name with the RMF file
 
         extension : str (default None)
-            FITS file extinction keyword, if the response matrix is stored
-            under an extension other tham "MATRIX" or "SPECRESP MATRIX"
+            FITS file extension keyword, if the response matrix is stored
+            under an extension other than "MATRIX" or "SPECRESP MATRIX"
 
         Attributes
         ----------
@@ -466,7 +466,7 @@ def _load_arf(self, filename):
         Parameters
         ----------
         filename : str
-            The file name with the RMF file
+            The file name with the ARF file
 
         Attributes
         ----------

From aa458114d39d123c8fb0e93a00f84a7e492bbf94 Mon Sep 17 00:00:00 2001
From: eblur <lia@astro.wisc.edu>
Date: Fri, 1 Jun 2018 11:19:14 -0500
Subject: [PATCH 58/89] Removing superfluous returns and pass

---
 specutils/spectra/xrayspectrum1d.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index e2eb6ce9c..610dea245 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -91,7 +91,6 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
         self.exposure = exposure
         self.assign_rmf(rmf)
         self.assign_arf(arf)
-        return
 
     # Convenience function for Xray people
     @property
@@ -115,7 +114,6 @@ def assign_arf(self, arf_inp):
             self.arf = ARF(arf_inp)
         else:
             self.arf = arf_inp
-        return
 
     def assign_rmf(self, rmf_inp):
         """
@@ -190,7 +188,6 @@ def apply_response(self, mflux, exposure=None):
 class RMF(object):
     def __init__(self, filename):
         self._load_rmf(filename)
-        pass
 
     def _load_rmf(self, filename, extension=None):
         """

From 674ae640e7627528b241b780513ac15d9f55a882 Mon Sep 17 00:00:00 2001
From: eblur <lia@astro.wisc.edu>
Date: Fri, 1 Jun 2018 11:22:03 -0500
Subject: [PATCH 59/89] Catch ValueError

---
 specutils/spectra/xrayspectrum1d.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 610dea245..34b34e7d8 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -79,7 +79,7 @@ def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
                  arf=None, rmf=None, rest_value=0.0 * u.angstrom, **kwargs):
         try:
             axis_unit = u.Unit(bin_unit)
-        except:
+        except ValueError:
             axis_unit = _unit_parser(bin_unit)
 
         bin_mid = 0.5 * (bin_lo + bin_hi) * axis_unit

From 029f8abee7243cdec11506162860e90fec4e5e99 Mon Sep 17 00:00:00 2001
From: eblur <lia@astro.wisc.edu>
Date: Fri, 1 Jun 2018 11:27:25 -0500
Subject: [PATCH 60/89] Spelled out ARF and RMF class names

---
 specutils/spectra/xrayspectrum1d.py | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 34b34e7d8..3c452c7d3 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -5,7 +5,7 @@
 from astropy.io import fits
 from .spectrum1d import Spectrum1D
 
-__all__ = ['XraySpectrum1D', 'ARF', 'RMF']
+__all__ = ['XraySpectrum1D', 'AreaResponse', 'ResponseMatrix']
 
 # For dealing with varied unit string choices
 EV   = ['eV', 'ev']
@@ -47,12 +47,12 @@ class XraySpectrum1D(Spectrum1D):
     exposure : float
         Exposure time for the dataset
 
-    arf : specutils.ARF or string, default None
-        Strings will be passed to ARF.__init__
+    arf : specutils.AreaResponse or string, default None
+        Strings will be passed to AreaResponse.__init__
         All other input types are stored as arf attribute
 
-    rmf : specutils.RMF or string, default None
-        Strings will be passed to RMF.__init__
+    rmf : specutils.ResponseMatrix or string, default None
+        Strings will be passed to ResponseMatrix.__init__
         All other input types are stored as rmf attribute
 
     rest_value : astropy.units.Quantity, default 0 Angstrom
@@ -111,7 +111,7 @@ def assign_arf(self, arf_inp):
         Modifies the XraySpectrum1D.arf attribute
         """
         if isinstance(arf_inp, str):
-            self.arf = ARF(arf_inp)
+            self.arf = AreaResponse(arf_inp)
         else:
             self.arf = arf_inp
 
@@ -129,7 +129,7 @@ def assign_rmf(self, rmf_inp):
         Modifies the XraySpectrum1D.rmf attribute
         """
         if isinstance(rmf_inp, str):
-            self.rmf = RMF(rmf_inp)
+            self.rmf = ResponseMatrix(rmf_inp)
         else:
             self.rmf = rmf_inp
         return
@@ -185,7 +185,7 @@ def apply_response(self, mflux, exposure=None):
 
 ## ----  Supporting response file objects
 
-class RMF(object):
+class ResponseMatrix(object):
     def __init__(self, filename):
         self._load_rmf(filename)
 
@@ -449,7 +449,7 @@ def apply_rmf(self, spec):
         return counts[:self.detchans]
 
 
-class ARF(object):
+class AreaResponse(object):
 
     def __init__(self, filename):
 

From 2590b7a05e957ca5cf77213916a35a7e9982fdfc Mon Sep 17 00:00:00 2001
From: eblur <lia@astro.wisc.edu>
Date: Fri, 1 Jun 2018 11:30:06 -0500
Subject: [PATCH 61/89] updated docstring

---
 specutils/spectra/xrayspectrum1d.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 3c452c7d3..cb99d21f8 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -137,8 +137,8 @@ def assign_rmf(self, rmf_inp):
     def apply_response(self, mflux, exposure=None):
         """
         Given a model flux spectrum, apply the response. In cases where the
-        spectrum has both an ARF and an RMF, apply both. Otherwise, apply
-        whatever response is in RMF.
+        spectrum has both an area response file (ARF) and a response matrix
+        file (RMF), apply both. Otherwise, apply whatever response is in RMF.
 
         The model flux spectrum *must* be created using the same units and
         bins as in the ARF (where the ARF exists)!

From f6b66f52196879cdce5f63854b19f563c6172944 Mon Sep 17 00:00:00 2001
From: eblur <lia@astro.wisc.edu>
Date: Fri, 1 Jun 2018 11:34:47 -0500
Subject: [PATCH 62/89] Updated tests

---
 specutils/tests/test_xrayspectrum1d.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/specutils/tests/test_xrayspectrum1d.py b/specutils/tests/test_xrayspectrum1d.py
index 53fbae5c0..1e1f45bbe 100644
--- a/specutils/tests/test_xrayspectrum1d.py
+++ b/specutils/tests/test_xrayspectrum1d.py
@@ -29,6 +29,6 @@ def test_apply_model():
     test_spec = test_create_from_arrays()
     new_model = PowerLaw1D(amplitude=5.e-3, alpha=1.8, x_0=1.e3)
     model_flux = new_model(test_spec.spectral_axis.value)
-    ymodel = test_spec.apply_resp(model_flux)
+    ymodel = test_spec.apply_response(model_flux)
     assert len(ymodel) == len(test_spec.spectral_axis)
     return

From 4111295c4a990873bb7f75ac0c7d221e212e3fc2 Mon Sep 17 00:00:00 2001
From: Lia Corrales <lia@astro.wisc.edu>
Date: Tue, 12 Jun 2018 20:38:39 +0200
Subject: [PATCH 63/89] initialize attributes to None before loading from file

---
 specutils/spectra/xrayspectrum1d.py | 31 +++++++++++++++++++++--------
 1 file changed, 23 insertions(+), 8 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index cb99d21f8..5d98c79bc 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -187,9 +187,6 @@ def apply_response(self, mflux, exposure=None):
 
 class ResponseMatrix(object):
     def __init__(self, filename):
-        self._load_rmf(filename)
-
-    def _load_rmf(self, filename, extension=None):
         """
         Load a response matrix file (RMF) from a FITS file.
 
@@ -237,6 +234,19 @@ def _load_rmf(self, filename, extension=None):
             The number of channels in the detector
 
         """
+        self.filename = ""
+        self.offset = None
+        self.n_grp = None
+        self.f_chan = None
+        self.n_chan = None
+        self.matrix = None
+        self.energ_lo = None
+        self.energ_hi = None
+        self.energ_unit = None
+        self.detchans = None
+        self._load_rmf(filename)
+
+    def _load_rmf(self, filename, extension=None):
         # open the FITS file and extract the MATRIX extension
         # which contains the redistribution matrix and
         # anxillary information
@@ -452,11 +462,6 @@ def apply_rmf(self, spec):
 class AreaResponse(object):
 
     def __init__(self, filename):
-
-        self._load_arf(filename)
-        pass
-
-    def _load_arf(self, filename):
         """
         Load an area response file (ARF) from a FITS file.
 
@@ -492,6 +497,16 @@ def _load_arf(self, filename):
             Average exposure time for the dataset
             (takes telescope dithering into account)
         """
+        self.filename = ""
+        self.e_low = None
+        self.e_high = None
+        self.e_unit = None
+        self.specresp = None
+        self.fracexpo = None
+        self.exposure = None
+        self._load_arf(filename)
+
+    def _load_arf(self, filename):
         # open the FITS file and extract the MATRIX extension
         # which contains the redistribution matrix and
         # anxillary information

From 8a1aa62787b3058068ca3e1c63fbaf9eaaf3f097 Mon Sep 17 00:00:00 2001
From: Lia Corrales <lia@astro.wisc.edu>
Date: Tue, 12 Jun 2018 20:39:55 +0200
Subject: [PATCH 64/89] save filename attribute for attempting load function

---
 specutils/spectra/xrayspectrum1d.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 5d98c79bc..58dcc3b73 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -234,7 +234,7 @@ def __init__(self, filename):
             The number of channels in the detector
 
         """
-        self.filename = ""
+        self.filename = filename
         self.offset = None
         self.n_grp = None
         self.f_chan = None
@@ -251,7 +251,6 @@ def _load_rmf(self, filename, extension=None):
         # which contains the redistribution matrix and
         # anxillary information
         hdulist = fits.open(filename)
-        self.filename = filename
 
         # get all the extension names
         extnames = np.array([h.name for h in hdulist])
@@ -497,7 +496,7 @@ def __init__(self, filename):
             Average exposure time for the dataset
             (takes telescope dithering into account)
         """
-        self.filename = ""
+        self.filename = filename
         self.e_low = None
         self.e_high = None
         self.e_unit = None
@@ -511,7 +510,6 @@ def _load_arf(self, filename):
         # which contains the redistribution matrix and
         # anxillary information
         hdulist = fits.open(filename)
-        self.filename = filename
 
         h = hdulist["SPECRESP"]
         data = h.data

From 4afde8f7e1d28f2146f8e50a57ac850a655a79c2 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Thu, 6 Dec 2018 13:30:05 -0700
Subject: [PATCH 65/89] [skip ci] provide xrayspectrum1d

---
 specutils/spectra/__init__.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/specutils/spectra/__init__.py b/specutils/spectra/__init__.py
index 05d5d14c4..d03dac32a 100644
--- a/specutils/spectra/__init__.py
+++ b/specutils/spectra/__init__.py
@@ -3,10 +3,8 @@
 `~astropy.nddata.NDData`-inherited classes used for storing the spectrum data.
 """
 from .spectrum1d import *  # noqa
-<<<<<<< HEAD
 from .spectral_region import *  # noqa
 from .spectrum_collection import *  #noqa
 from .spectrum_list import * #noqa
-=======
+
 from .xrayspectrum1d import *
->>>>>>> Can load an X-ray spectrum

From 90f8ad6871bde4a1b0353bc9a85f381682369458 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Thu, 6 Dec 2018 13:57:25 -0700
Subject: [PATCH 66/89] specutils flux input needs units

---
 specutils/tests/test_xrayspectrum1d.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/specutils/tests/test_xrayspectrum1d.py b/specutils/tests/test_xrayspectrum1d.py
index 1e1f45bbe..4317bc836 100644
--- a/specutils/tests/test_xrayspectrum1d.py
+++ b/specutils/tests/test_xrayspectrum1d.py
@@ -11,7 +11,7 @@ def test_create_from_arrays():
     elo = energy[:-1]
     ehi = energy[1:]
     emid = 0.5 * (elo + ehi)
-    counts = np.random.poisson(lam=powlaw0(emid), size=len(emid))
+    counts = np.random.poisson(lam=powlaw0(emid), size=len(emid)) * u.ct
     test_spec = XraySpectrum1D(elo, ehi, u.keV, counts, exposure=1.0)
     return test_spec
 

From fa80d2440a18430d62541d895da0681977d1650a Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Thu, 6 Dec 2018 14:32:03 -0700
Subject: [PATCH 67/89] attempting to fix docstrings

---
 specutils/spectra/xrayspectrum1d.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 58dcc3b73..c82d0cced 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -541,10 +541,12 @@ def apply_arf(self, spec, exposure=None):
         The ARF is a single vector encoding the effective area information
         about the detector. A such, applying the ARF is a simple
         multiplication with the input spectrum.
+
         Parameters
         ----------
         spec : numpy.ndarray
             The (model) spectrum to be folded
+
         exposure : float, default None
             Value for the exposure time. By default, `apply_arf` will use the
             exposure keyword from the ARF file. If this exposure time is not
@@ -552,6 +554,7 @@ def apply_arf(self, spec, exposure=None):
             time and the ARF from a real observation), one can override the
             default exposure by setting the `exposure` keyword to the correct
             value.
+
         Returns
         -------
         s_arf : numpy.ndarray

From f14e9afd4095372b2b0359d043517f62eb859021 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Thu, 6 Dec 2018 14:54:07 -0700
Subject: [PATCH 68/89] further docstring fixes for sphinx

---
 specutils/spectra/xrayspectrum1d.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index c82d0cced..ce0771053 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -41,7 +41,7 @@ class XraySpectrum1D(Spectrum1D):
     bin_unit : string OR astropy.units.Unit
         Unit for the bin values
 
-    counts : numpy.ndarray
+    counts : astropy.Quantity array
         Counts histogram for the X-ray spectrum
 
     exposure : float
@@ -374,10 +374,13 @@ def apply_rmf(self, spec):
         vector, there exists one or more sets of channels that this
         flux is redistributed into. The additional arrays `n_grp`,
         `f_chan` and `n_chan` store this information:
+        
             * `n_group` stores the number of channel groups for each
               energy bin
+
             * `f_chan` stores the *first channel* that each channel
               for each channel set
+
             * `n_chan` stores the number of channels in each channel
               set
 

From eecce756203afcafdaac91e80a05ce7769b56e25 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Thu, 6 Dec 2018 17:18:41 -0700
Subject: [PATCH 69/89] added extension keyword to _load_arf

---
 specutils/spectra/xrayspectrum1d.py | 25 +++++++++++++++++++------
 1 file changed, 19 insertions(+), 6 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index ce0771053..102bc16d4 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -374,7 +374,7 @@ def apply_rmf(self, spec):
         vector, there exists one or more sets of channels that this
         flux is redistributed into. The additional arrays `n_grp`,
         `f_chan` and `n_chan` store this information:
-        
+
             * `n_group` stores the number of channel groups for each
               energy bin
 
@@ -508,13 +508,26 @@ def __init__(self, filename):
         self.exposure = None
         self._load_arf(filename)
 
-    def _load_arf(self, filename):
-        # open the FITS file and extract the MATRIX extension
-        # which contains the redistribution matrix and
-        # anxillary information
+    def _load_arf(self, filename, extension=None):
+        # open the FITS file and extract the SPECRESP extension
+        # which contains the spectral response for the telescope
         hdulist = fits.open(filename)
 
-        h = hdulist["SPECRESP"]
+        # get all the extension names
+        extnames = np.array([h.name for h in hdulist])
+
+        # figure out the right extension to use
+        if extension is not None:
+            h = hdulist[extension]
+
+        elif 'SPECRESP' in extnames:
+            h = hdulist['SPECRESP']
+
+        else:
+            print("Cannot find common FITS file extension for spectral response")
+            print("Please set the `extension` keyword")
+            return
+
         data = h.data
         hdr = h.header
         hdulist.close()

From bdc5a418efbe49e659f015baffe7e85dbeef232d Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Thu, 6 Dec 2018 17:22:01 -0700
Subject: [PATCH 70/89] updated init and load calls with extension keyword

---
 specutils/spectra/xrayspectrum1d.py | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 102bc16d4..5ed8365e7 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -186,7 +186,7 @@ def apply_response(self, mflux, exposure=None):
 ## ----  Supporting response file objects
 
 class ResponseMatrix(object):
-    def __init__(self, filename):
+    def __init__(self, filename, extension=None):
         """
         Load a response matrix file (RMF) from a FITS file.
 
@@ -244,7 +244,7 @@ def __init__(self, filename):
         self.energ_hi = None
         self.energ_unit = None
         self.detchans = None
-        self._load_rmf(filename)
+        self._load_rmf(filename, extension=extension)
 
     def _load_rmf(self, filename, extension=None):
         # open the FITS file and extract the MATRIX extension
@@ -463,7 +463,7 @@ def apply_rmf(self, spec):
 
 class AreaResponse(object):
 
-    def __init__(self, filename):
+    def __init__(self, filename, extension=None):
         """
         Load an area response file (ARF) from a FITS file.
 
@@ -472,6 +472,10 @@ def __init__(self, filename):
         filename : str
             The file name with the ARF file
 
+        extension : str (default None)
+            FITS file extension keyword, if the spectral response is stored
+            under an extension other than "SPECRESP"
+
         Attributes
         ----------
         filename : str
@@ -506,7 +510,7 @@ def __init__(self, filename):
         self.specresp = None
         self.fracexpo = None
         self.exposure = None
-        self._load_arf(filename)
+        self._load_arf(filename, extension=extension)
 
     def _load_arf(self, filename, extension=None):
         # open the FITS file and extract the SPECRESP extension

From 3fef5b91c9fece08e5005887ab6b3b660510c3b2 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Fri, 7 Dec 2018 13:53:48 -0700
Subject: [PATCH 71/89] XraySpectrum1d uses astropy.Quantity for all inputs

---
 specutils/spectra/xrayspectrum1d.py | 17 +++++------------
 1 file changed, 5 insertions(+), 12 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 5ed8365e7..28bfb7742 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -32,19 +32,16 @@ class XraySpectrum1D(Spectrum1D):
 
     Parameters
     ----------
-    bin_lo : numpy.ndarray
+    bin_lo : astropy.Quantity
         The left edges for bin values
 
-    bin_hi : numpy.ndarray
+    bin_hi : astropy.Quantity
         The right edges for bin values
 
-    bin_unit : string OR astropy.units.Unit
-        Unit for the bin values
-
-    counts : astropy.Quantity array
+    counts : astropy.Quantity
         Counts histogram for the X-ray spectrum
 
-    exposure : float
+    exposure : astropy.Quantity
         Exposure time for the dataset
 
     arf : specutils.AreaResponse or string, default None
@@ -77,12 +74,8 @@ class XraySpectrum1D(Spectrum1D):
     """
     def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
                  arf=None, rmf=None, rest_value=0.0 * u.angstrom, **kwargs):
-        try:
-            axis_unit = u.Unit(bin_unit)
-        except ValueError:
-            axis_unit = _unit_parser(bin_unit)
 
-        bin_mid = 0.5 * (bin_lo + bin_hi) * axis_unit
+        bin_mid = 0.5 * (bin_lo + bin_hi)
         Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts,
                             rest_value=rest_value, **kwargs)
 

From 8b4b945e83c66c7edae91e3e6b30ca4d1e5c6207 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Fri, 7 Dec 2018 13:58:56 -0700
Subject: [PATCH 72/89] tested xrayspectrum1d with quantity inputs

---
 specutils/spectra/xrayspectrum1d.py    | 2 +-
 specutils/tests/test_xrayspectrum1d.py | 8 ++++----
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 28bfb7742..68f12da3f 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -72,7 +72,7 @@ class XraySpectrum1D(Spectrum1D):
 
     rmf
     """
-    def __init__(self, bin_lo, bin_hi, bin_unit, counts, exposure,
+    def __init__(self, bin_lo, bin_hi, counts, exposure,
                  arf=None, rmf=None, rest_value=0.0 * u.angstrom, **kwargs):
 
         bin_mid = 0.5 * (bin_lo + bin_hi)
diff --git a/specutils/tests/test_xrayspectrum1d.py b/specutils/tests/test_xrayspectrum1d.py
index 4317bc836..6a561307f 100644
--- a/specutils/tests/test_xrayspectrum1d.py
+++ b/specutils/tests/test_xrayspectrum1d.py
@@ -8,11 +8,11 @@ def test_create_from_arrays():
     amp0, alpha0 = 3.e-3, 2.0
     powlaw0 = PowerLaw1D(amplitude=amp0, alpha=alpha0, x_0=1.e3)
     energy = np.linspace(0.2, 10.0, 8000)
-    elo = energy[:-1]
-    ehi = energy[1:]
+    elo = energy[:-1] * u.keV
+    ehi = energy[1:] * u.keV
     emid = 0.5 * (elo + ehi)
-    counts = np.random.poisson(lam=powlaw0(emid), size=len(emid)) * u.ct
-    test_spec = XraySpectrum1D(elo, ehi, u.keV, counts, exposure=1.0)
+    counts = np.random.poisson(lam=powlaw0(emid.value), size=len(emid)) * u.ct
+    test_spec = XraySpectrum1D(elo, ehi, counts, exposure=1.0*u.second)
     return test_spec
 
 def test_call_spectrum():

From fd052cd5b1d6811477e1a9565e89d5083af2cfd2 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Fri, 7 Dec 2018 10:54:52 -0700
Subject: [PATCH 73/89] first pass at hetg loader

---
 specutils/io/default_loaders/chandra_hetg.py | 53 ++++++++++++++++++++
 1 file changed, 53 insertions(+)
 create mode 100644 specutils/io/default_loaders/chandra_hetg.py

diff --git a/specutils/io/default_loaders/chandra_hetg.py b/specutils/io/default_loaders/chandra_hetg.py
new file mode 100644
index 000000000..0c0136f84
--- /dev/null
+++ b/specutils/io/default_loaders/chandra_hetg.py
@@ -0,0 +1,53 @@
+import os
+
+from astropy.io import fits
+from astropy.units import Unit
+#from astropy.nddata import StdDevUncertainty
+
+from specutils.io.registers import data_loader
+from specutils import XraySpectrum1D
+
+__all__ = ['chandra_identify', 'hetg_spectrum_loader']
+
+def chandra_identify(origin, *args, **kwargs):
+    """Check whether given file contains Chandra HETG spectral data."""
+    with fits.open(args[0]) as hdu:
+        # Check that it is a Chandra spectrum file
+        if hdu[0].header['TELESCOP'] == 'CHANDRA' and 'SPECTRUM' in [h.name for h in hdu]:
+            # Check that it's an HETG file
+            if hdu['SPECTRUM'].header['GRATING'] == 'HETG':
+                return True
+    return False
+
+@data_loader(label="chandra_hetg", identifier=chandra_identify, extensions=['pha','pha2', 'FITS', 'FIT', 'fits', 'fit'])
+def hetg_spectrum_loader(file_name. arf=None, rmf=None):
+    """
+    Load Chandra HETG spectral data from a file into a spectrum object.
+
+    Parameters
+    ----------
+    file_name: str
+        The path to the FITS file
+
+    Returns
+    -------
+    data: Spectrum1D
+        The spectrum that is represented by the data in this table.
+    """
+
+    name = os.path.basename(file_name)
+
+    with fits.open(file_name) as hdu:
+        header = hdu[0].header
+        meta = {'header': header}
+
+        data     = hdu[1].data
+        bin_lo   = data['BIN_LO']
+        bin_hi   = data['BIN_HI']
+
+        bin_unit = data.columns['BIN_LO'].unit
+        counts   = data['COUNTS'] * u.ct
+        exposure = hdu[1].header['EXPOSURE']  # seconds
+
+    return XraySpectrum1D(bin_lo, bin_hi, bin_unit, counts,
+                          exposure=exposure, arf=arf, rmf=rmf)

From 7cc6af63109442edf9cda7d49e42b4ad5e5e1110 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Fri, 7 Dec 2018 10:58:03 -0700
Subject: [PATCH 74/89] fixed typo

---
 specutils/io/default_loaders/chandra_hetg.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/specutils/io/default_loaders/chandra_hetg.py b/specutils/io/default_loaders/chandra_hetg.py
index 0c0136f84..a7b486995 100644
--- a/specutils/io/default_loaders/chandra_hetg.py
+++ b/specutils/io/default_loaders/chandra_hetg.py
@@ -20,7 +20,7 @@ def chandra_identify(origin, *args, **kwargs):
     return False
 
 @data_loader(label="chandra_hetg", identifier=chandra_identify, extensions=['pha','pha2', 'FITS', 'FIT', 'fits', 'fit'])
-def hetg_spectrum_loader(file_name. arf=None, rmf=None):
+def hetg_spectrum_loader(file_name, arf=None, rmf=None):
     """
     Load Chandra HETG spectral data from a file into a spectrum object.
 

From ee33c6437a786efcb48b882a559ef748decd89c1 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Fri, 7 Dec 2018 10:59:22 -0700
Subject: [PATCH 75/89] it runs in my test notebook

---
 specutils/io/default_loaders/chandra_hetg.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/specutils/io/default_loaders/chandra_hetg.py b/specutils/io/default_loaders/chandra_hetg.py
index a7b486995..a4bcd7b4a 100644
--- a/specutils/io/default_loaders/chandra_hetg.py
+++ b/specutils/io/default_loaders/chandra_hetg.py
@@ -46,7 +46,7 @@ def hetg_spectrum_loader(file_name, arf=None, rmf=None):
         bin_hi   = data['BIN_HI']
 
         bin_unit = data.columns['BIN_LO'].unit
-        counts   = data['COUNTS'] * u.ct
+        counts   = data['COUNTS'] * Unit('ct')
         exposure = hdu[1].header['EXPOSURE']  # seconds
 
     return XraySpectrum1D(bin_lo, bin_hi, bin_unit, counts,

From 227520d5d0e83bf4bf34ecb29698c872660afb71 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Fri, 7 Dec 2018 11:07:53 -0700
Subject: [PATCH 76/89] updated docstrings

---
 specutils/io/default_loaders/chandra_hetg.py | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/specutils/io/default_loaders/chandra_hetg.py b/specutils/io/default_loaders/chandra_hetg.py
index a4bcd7b4a..3194973f2 100644
--- a/specutils/io/default_loaders/chandra_hetg.py
+++ b/specutils/io/default_loaders/chandra_hetg.py
@@ -29,9 +29,15 @@ def hetg_spectrum_loader(file_name, arf=None, rmf=None):
     file_name: str
         The path to the FITS file
 
+    arf : str OR AreaResponse
+        Filename for the area response file (ARF) or a pre-loaded AreaResponse object
+
+    rmf : str OR ResponseMatrix
+        Filename for the response matrix file (RMF) or a pre-loaded ResponseMatrix object
+
     Returns
     -------
-    data: Spectrum1D
+    data: XraySpectrum1D
         The spectrum that is represented by the data in this table.
     """
 
@@ -47,7 +53,7 @@ def hetg_spectrum_loader(file_name, arf=None, rmf=None):
 
         bin_unit = data.columns['BIN_LO'].unit
         counts   = data['COUNTS'] * Unit('ct')
-        exposure = hdu[1].header['EXPOSURE']  # seconds
+        exposure = hdu[1].header['EXPOSURE'] * Unit('second')
 
     return XraySpectrum1D(bin_lo, bin_hi, bin_unit, counts,
                           exposure=exposure, arf=arf, rmf=rmf)

From 66173d6f6311f1f8607241930cf5818c4af29eb7 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Fri, 7 Dec 2018 11:09:53 -0700
Subject: [PATCH 77/89] added file extensions to the list

---
 specutils/io/default_loaders/chandra_hetg.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/specutils/io/default_loaders/chandra_hetg.py b/specutils/io/default_loaders/chandra_hetg.py
index 3194973f2..80b0824af 100644
--- a/specutils/io/default_loaders/chandra_hetg.py
+++ b/specutils/io/default_loaders/chandra_hetg.py
@@ -19,7 +19,7 @@ def chandra_identify(origin, *args, **kwargs):
                 return True
     return False
 
-@data_loader(label="chandra_hetg", identifier=chandra_identify, extensions=['pha','pha2', 'FITS', 'FIT', 'fits', 'fit'])
+@data_loader(label="chandra_hetg", identifier=chandra_identify, extensions=['pha','pha2', 'PHA', 'PHA2', 'FITS', 'FIT', 'fits', 'fit'])
 def hetg_spectrum_loader(file_name, arf=None, rmf=None):
     """
     Load Chandra HETG spectral data from a file into a spectrum object.

From 0e814f38379364940f8ed7f7ecb99e0546f01f52 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Fri, 7 Dec 2018 14:02:00 -0700
Subject: [PATCH 78/89] rewrote to use astropy.Quantity at all times

---
 specutils/io/default_loaders/chandra_hetg.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/specutils/io/default_loaders/chandra_hetg.py b/specutils/io/default_loaders/chandra_hetg.py
index 80b0824af..d7366fdfd 100644
--- a/specutils/io/default_loaders/chandra_hetg.py
+++ b/specutils/io/default_loaders/chandra_hetg.py
@@ -45,15 +45,15 @@ def hetg_spectrum_loader(file_name, arf=None, rmf=None):
 
     with fits.open(file_name) as hdu:
         header = hdu[0].header
-        meta = {'header': header}
+        meta   = {'header': header}
+        data   = hdu[1].data
 
-        data     = hdu[1].data
-        bin_lo   = data['BIN_LO']
-        bin_hi   = data['BIN_HI']
+        bin_unit = Unit(data.columns['BIN_LO'].unit)
+        bin_lo   = data['BIN_LO'] * bin_unit
+        bin_hi   = data['BIN_HI'] * bin_unit
 
-        bin_unit = data.columns['BIN_LO'].unit
         counts   = data['COUNTS'] * Unit('ct')
         exposure = hdu[1].header['EXPOSURE'] * Unit('second')
 
-    return XraySpectrum1D(bin_lo, bin_hi, bin_unit, counts,
+    return XraySpectrum1D(bin_lo, bin_hi, counts,
                           exposure=exposure, arf=arf, rmf=rmf)

From 1edf177b5b871ff8c5a655b962efa2c106ee10db Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Wed, 26 Dec 2018 16:51:14 -0700
Subject: [PATCH 79/89] Moved arf loader to an outside function

---
 specutils/io/default_loaders/chandra_hetg.py |   7 +-
 specutils/spectra/xrayspectrum1d.py          | 193 +++++++++----------
 2 files changed, 98 insertions(+), 102 deletions(-)

diff --git a/specutils/io/default_loaders/chandra_hetg.py b/specutils/io/default_loaders/chandra_hetg.py
index d7366fdfd..306dab48c 100644
--- a/specutils/io/default_loaders/chandra_hetg.py
+++ b/specutils/io/default_loaders/chandra_hetg.py
@@ -2,14 +2,13 @@
 
 from astropy.io import fits
 from astropy.units import Unit
-#from astropy.nddata import StdDevUncertainty
 
 from specutils.io.registers import data_loader
 from specutils import XraySpectrum1D
 
-__all__ = ['chandra_identify', 'hetg_spectrum_loader']
+__all__ = ['chandra_hetg_identify', 'hetg_spectrum_loader']
 
-def chandra_identify(origin, *args, **kwargs):
+def chandra_hetg_identify(origin, *args, **kwargs):
     """Check whether given file contains Chandra HETG spectral data."""
     with fits.open(args[0]) as hdu:
         # Check that it is a Chandra spectrum file
@@ -19,7 +18,7 @@ def chandra_identify(origin, *args, **kwargs):
                 return True
     return False
 
-@data_loader(label="chandra_hetg", identifier=chandra_identify, extensions=['pha','pha2', 'PHA', 'PHA2', 'FITS', 'FIT', 'fits', 'fit'])
+@data_loader(label="chandra_hetg", identifier=chandra_hetg_identify, extensions=['pha','pha2', 'PHA', 'PHA2', 'FITS', 'FIT', 'fits', 'fit'])
 def hetg_spectrum_loader(file_name, arf=None, rmf=None):
     """
     Load Chandra HETG spectral data from a file into a spectrum object.
diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 68f12da3f..b70929caf 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -5,7 +5,7 @@
 from astropy.io import fits
 from .spectrum1d import Spectrum1D
 
-__all__ = ['XraySpectrum1D', 'AreaResponse', 'ResponseMatrix']
+__all__ = ['XraySpectrum1D', 'AreaResponse', 'ResponseMatrix', 'arf_loader']
 
 # For dealing with varied unit string choices
 EV   = ['eV', 'ev']
@@ -90,7 +90,7 @@ def __init__(self, bin_lo, bin_hi, counts, exposure,
     def counts(self):
         return self.flux
 
-    def assign_arf(self, arf_inp):
+    def assign_arf(self, arf_inp, **kwargs):
         """
         Assign an area response file (ARF) object to the XraySpectrum1D object
 
@@ -104,7 +104,7 @@ def assign_arf(self, arf_inp):
         Modifies the XraySpectrum1D.arf attribute
         """
         if isinstance(arf_inp, str):
-            self.arf = AreaResponse(arf_inp)
+            self.arf = arf_loader(arf_inp, **kwargs)
         else:
             self.arf = arf_inp
 
@@ -139,7 +139,7 @@ def apply_response(self, mflux, exposure=None):
         Parameters
         ----------
         mflux : iterable
-            A list or array with the model flux values in ergs/keV/s/cm^-2
+            A list or array with the model flux values in phot/s/cm^-2
 
         exposure : float, default None
             By default, the exposure stored in the ARF will be used to compute
@@ -163,15 +163,13 @@ def apply_response(self, mflux, exposure=None):
         """
 
         if self.arf is not None:
-            mrate  = self.arf.apply_arf(mflux, exposure=exposure)
+            mrate  = self.arf.apply_arf(mflux, exposure=exposure) # ct/bin with no RMF applied
         else:
-            #print("Caution: no ARF file specified")
-            mrate = mflux
+            mrate = mflux # assumes ct/bin is the input with no RMF applied
 
         if self.rmf is not None:
-            result = self.rmf.apply_rmf(mrate)
+            result = self.rmf.apply_rmf(mrate.value)
         else:
-            #print("Caution: no RMF file specified")
             result = mrate
 
         return result
@@ -358,27 +356,27 @@ def _flatten_arrays(self, n_grp, f_chan, n_chan, matrix):
 
         return n_grp, f_chan_flat, n_chan_flat, matrix_flat
 
-    def apply_rmf(self, spec):
+    def apply_rmf(self, model):
         """
         Fold the spectrum through the redistribution matrix.
 
         The redistribution matrix is saved as a flattened 1-dimensional
         vector to save space. In reality, for each entry in the flux
         vector, there exists one or more sets of channels that this
-        flux is redistributed into. The additional arrays `n_grp`,
-        `f_chan` and `n_chan` store this information:
+        flux is redistributed into. The additional arrays n_grp,
+        f_chan and n_chan store this information:
 
-            * `n_group` stores the number of channel groups for each
+            * n_group stores the number of channel groups for each
               energy bin
 
-            * `f_chan` stores the *first channel* that each channel
+            * f_chan stores the *first channel* that each channel
               for each channel set
 
-            * `n_chan` stores the number of channels in each channel
+            * n_chan stores the number of channels in each channel
               set
 
         As a result, for a given energy bin i, we need to look up the
-        number of channel sets in `n_grp` for that energy bin. We
+        number of channel sets in n_grp for that energy bin. We
         then need to loop over the number of channel sets. For each
         channel set, we look up the first channel into which flux
         will be distributed as well as the number of channels in the
@@ -391,7 +389,7 @@ def apply_rmf(self, spec):
 
         Parameters
         ----------
-        spec : numpy.ndarray
+        model : numpy.ndarray
             The (model) spectrum to be folded
 
         Returns
@@ -402,7 +400,7 @@ def apply_rmf(self, spec):
 
         """
         # get the number of channels in the data
-        nchannels = spec.shape[0]
+        nchannels = model.shape[0]
 
         # an empty array for the output counts
         counts = np.zeros(nchannels)
@@ -418,7 +416,7 @@ def apply_rmf(self, spec):
 
             # this is the current bin in the flux spectrum to
             # be folded
-            source_bin_i = spec[i]
+            source_bin_i = model[i]
 
             # get the current number of groups
             current_num_groups = self.n_grp[i]
@@ -455,100 +453,49 @@ def apply_rmf(self, spec):
 
 
 class AreaResponse(object):
-
-    def __init__(self, filename, extension=None):
+    def __init__(self, e_low, e_high, eff_area,
+                 exposure=None, fracexpo=1.0, filename=None):
         """
         Load an area response file (ARF) from a FITS file.
 
-        Parameters
-        ----------
-        filename : str
-            The file name with the ARF file
-
-        extension : str (default None)
-            FITS file extension keyword, if the spectral response is stored
-            under an extension other than "SPECRESP"
-
         Attributes
         ----------
-        filename : str
-            The file name that the ARF was drawn from
-
-        e_low : numpy.ndarray
+        e_low : astropy.units.Quantity
             The lower edges of the energy bins
 
-        e_high : numpy.ndarray
+        e_high : astropy.units.Quantity
             The upper edges of the energy bins
 
-        e_unit : astropy.units.Unit
-            Description of the energy units used
+        eff_area : astropy.units.Quantity
+            Description of the energy dependent telescope effective area
 
-        specresp : numpy.ndarray
-            Description of the energy dependent telescope response area
+        exposure :
+            Average exposure time for the dataset
+            (which takes telescope dithering into account).
 
         fracexpo : float or numpy.ndarray
             Fractional exposure time for the spectrum
             (sometimes constant, sometimes dependent on spectral channel).
             These values are stored for reference; generally, they are already
-            accounted for in the specresp array.
+            accounted for in the eff_area array.
 
-        exposure :
-            Average exposure time for the dataset
-            (takes telescope dithering into account)
+        filename : str
+            Stores the name of the file from which the ARF was loaded.
+
+        e_mid : Property that returns the middle of each energy bin
         """
         self.filename = filename
-        self.e_low = None
-        self.e_high = None
-        self.e_unit = None
-        self.specresp = None
-        self.fracexpo = None
-        self.exposure = None
-        self._load_arf(filename, extension=extension)
-
-    def _load_arf(self, filename, extension=None):
-        # open the FITS file and extract the SPECRESP extension
-        # which contains the spectral response for the telescope
-        hdulist = fits.open(filename)
-
-        # get all the extension names
-        extnames = np.array([h.name for h in hdulist])
-
-        # figure out the right extension to use
-        if extension is not None:
-            h = hdulist[extension]
-
-        elif 'SPECRESP' in extnames:
-            h = hdulist['SPECRESP']
-
-        else:
-            print("Cannot find common FITS file extension for spectral response")
-            print("Please set the `extension` keyword")
-            return
-
-        data = h.data
-        hdr = h.header
-        hdulist.close()
-
-        # extract + store the attributes described in the docstring
-
-        self.e_low  = np.array(data.field("ENERG_LO"))
-        self.e_high = np.array(data.field("ENERG_HI"))
-        self.e_unit = _unit_parser(data.columns["ENERG_LO"].unit)
-        self.specresp = np.array(data.field("SPECRESP"))
-
-        if "EXPOSURE" in list(hdr.keys()):
-            self.exposure = hdr["EXPOSURE"]
-        else:
-            self.exposure = 1.0
-
-        if "FRACEXPO" in data.columns.names:
-            self.fracexpo = data["FRACEXPO"]
-        else:
-            self.fracexpo = 1.0
+        self.e_low    = e_low
+        self.e_high   = e_high
+        self.eff_area = eff_area
+        self.fracexpo = fracexpo
+        self.exposure = exposure
 
-        return
+    @property
+    def e_mid(self):
+        return 0.5 * (self.e_low + self.e_high)
 
-    def apply_arf(self, spec, exposure=None):
+    def apply_arf(self, model, exposure=None):
         """
         Fold the spectrum through the area response file (ARF).
         The ARF is a single vector encoding the effective area information
@@ -557,8 +504,8 @@ def apply_arf(self, spec, exposure=None):
 
         Parameters
         ----------
-        spec : numpy.ndarray
-            The (model) spectrum to be folded
+        model : numpy.ndarray
+            The model spectrum to which the arf will be applied
 
         exposure : float, default None
             Value for the exposure time. By default, `apply_arf` will use the
@@ -574,10 +521,60 @@ def apply_arf(self, spec, exposure=None):
             The (model) spectrum after folding, in
             counts/s/channel
         """
-        assert spec.shape[0] == self.specresp.shape[0], "The input spectrum must " \
+        assert model.shape[0] == self.eff_area.shape[0], "The input spectrum must " \
                                                       "be of same size as the " \
                                                       "ARF array."
         if exposure is None:
-            return np.array(spec) * self.specresp * self.exposure
+            return model * self.eff_area * self.exposure
         else:
-            return np.array(spec) * self.specresp * exposure
+            return model * self.eff_area * exposure
+
+## ----- Convenience functions for loading ARFs and RMFs
+def arf_loader(filename, block='SPECRESP'):
+    """
+    Load a AreaResponse object from FITS file.
+
+    Parameters
+    ----------
+    filename : str
+        Path to the FITS file
+
+    block : str
+        FITS file block keyword, if the spectral response is stored
+        under an extension other than "SPECRESP"
+
+    Returns
+    -------
+    AreaResponse
+        The ARF that is represented by the FITS file
+    """
+    # open the FITS file and extract the SPECRESP block
+    # which contains the spectral response for the telescope
+    hdulist = fits.open(filename)
+
+    # Get the data from the appropriate extension
+    blocknames = np.array([h.name for h in hdulist])
+    assert block in blocknames, "Could not fine {} block in FITS file".format(block)
+    h = hdulist[block]
+
+    data = h.data
+    hdr = h.header
+    hdulist.close()
+
+    e_unit = u.Unit(data.columns["ENERG_LO"].unit)
+    e_low  = np.array(data.field("ENERG_LO")) * e_unit
+    e_high = np.array(data.field("ENERG_HI")) * e_unit
+
+    area_unit = u.Unit(data.columns['SPECRESP'].unit) * u.ct # usually ct cm^2 / phot
+    specresp  = np.array(data.field("SPECRESP")) * area_unit
+
+    exposure = None
+    if "EXPOSURE" in list(hdr.keys()):
+        exposure = hdr["EXPOSURE"]
+
+    fracexpo = 1.0
+    if "FRACEXPO" in data.columns.names:
+        fracexpo = data["FRACEXPO"]
+
+    return AreaResponse(e_low, e_high, specresp,
+                        exposure=exposure, fracexpo=fracexpo, filename=filename)

From fe052a08c137ab9f06af5e721cf9a10b3c04801f Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Thu, 27 Dec 2018 10:47:20 -0700
Subject: [PATCH 80/89] now possible to make custom RMF

---
 specutils/spectra/xrayspectrum1d.py | 28 ++++++++++++++++------------
 1 file changed, 16 insertions(+), 12 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index b70929caf..33a240e5f 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -177,7 +177,10 @@ def apply_response(self, mflux, exposure=None):
 ## ----  Supporting response file objects
 
 class ResponseMatrix(object):
-    def __init__(self, filename, extension=None):
+    def __init__(self, filename=None, extension=None,
+                 energ_lo=None, energ_hi=None, matrix=None, energ_unit=None,
+                 offset=0.0, n_grp=np.array([]), f_chan=np.array([]),
+                 n_chan=np.array([]), detchans=0):
         """
         Load a response matrix file (RMF) from a FITS file.
 
@@ -226,16 +229,17 @@ def __init__(self, filename, extension=None):
 
         """
         self.filename = filename
-        self.offset = None
-        self.n_grp = None
-        self.f_chan = None
-        self.n_chan = None
-        self.matrix = None
-        self.energ_lo = None
-        self.energ_hi = None
-        self.energ_unit = None
-        self.detchans = None
-        self._load_rmf(filename, extension=extension)
+        self.offset = offset
+        self.n_grp = n_grp
+        self.f_chan = f_chan
+        self.n_chan = n_chan
+        self.matrix = matrix
+        self.energ_lo = energ_lo
+        self.energ_hi = energ_hi
+        self.energ_unit = energ_unit
+        self.detchans = detchans
+        if filename is not None:
+            self._load_rmf(filename, extension=extension)
 
     def _load_rmf(self, filename, extension=None):
         # open the FITS file and extract the MATRIX extension
@@ -273,7 +277,7 @@ def _load_rmf(self, filename, extension=None):
 
         self.energ_lo = np.array(data.field("ENERG_LO"))
         self.energ_hi = np.array(data.field("ENERG_HI"))
-        self.energ_unit = _unit_parser(data.columns["ENERG_LO"].unit)
+        self.energ_unit = u.Unit(data.columns["ENERG_LO"].unit)
         self.detchans = hdr["DETCHANS"]
         self.offset = self.__get_tlmin(h)
 

From fd3455cff254275cb02c4ea88cfd9026fa2f01ad Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Thu, 27 Dec 2018 10:50:17 -0700
Subject: [PATCH 81/89] docstring reflects changes to RMF init

---
 specutils/spectra/xrayspectrum1d.py | 32 ++++++++++++-----------------
 1 file changed, 13 insertions(+), 19 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 33a240e5f..506b68462 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -184,8 +184,6 @@ def __init__(self, filename=None, extension=None,
         """
         Load a response matrix file (RMF) from a FITS file.
 
-        Parameters
-        ----------
         filename : str
             The file name with the RMF file
 
@@ -193,10 +191,19 @@ def __init__(self, filename=None, extension=None,
             FITS file extension keyword, if the response matrix is stored
             under an extension other than "MATRIX" or "SPECRESP MATRIX"
 
-        Attributes
-        ----------
-        filename : str
-            The file name that the RMF was drawn from
+        energ_lo : numpy.ndarray
+            The lower edges of the energy bins
+
+        energ_hi : numpy.ndarray
+            The upper edges of the energy bins
+
+        matrix : numpy.ndarray
+            The redistribution matrix as a flattened 1D vector
+
+        energ_unit : astropy.units.Unit
+            Description of the energy units used
+
+        offset : float
 
         n_grp : numpy.ndarray
             the Array with the number of channels in each
@@ -212,21 +219,8 @@ def __init__(self, filename=None, extension=None,
             The number of channels in each channel group. The same
             logic as for f_chan applies
 
-        matrix : numpy.ndarray
-            The redistribution matrix as a flattened 1D vector
-
-        energ_lo : numpy.ndarray
-            The lower edges of the energy bins
-
-        energ_hi : numpy.ndarray
-            The upper edges of the energy bins
-
-        energ_unit : astropy.units.Unit
-            Description of the energy units used
-
         detchans : int
             The number of channels in the detector
-
         """
         self.filename = filename
         self.offset = offset

From 0d9ce500eb91758784eff69fa294a792f922cd82 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Thu, 27 Dec 2018 15:06:04 -0700
Subject: [PATCH 82/89] Added unit to arf exposure time

---
 specutils/spectra/xrayspectrum1d.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 506b68462..0e801fc6b 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -568,7 +568,7 @@ def arf_loader(filename, block='SPECRESP'):
 
     exposure = None
     if "EXPOSURE" in list(hdr.keys()):
-        exposure = hdr["EXPOSURE"]
+        exposure = hdr["EXPOSURE"] * u.second
 
     fracexpo = 1.0
     if "FRACEXPO" in data.columns.names:

From 9d525623ea71d69462ceb564a18cf3b0920a9a3c Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Thu, 27 Dec 2018 16:23:24 -0700
Subject: [PATCH 83/89] updated docstring

---
 specutils/spectra/xrayspectrum1d.py | 16 +++++-----------
 1 file changed, 5 insertions(+), 11 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 0e801fc6b..736f252f9 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -138,10 +138,11 @@ def apply_response(self, mflux, exposure=None):
 
         Parameters
         ----------
-        mflux : iterable
-            A list or array with the model flux values in phot/s/cm^-2
+        mflux : astropy.units.Quantity
+            A list or array with the model flux values,
+            typically with units of phot/s/cm^-2
 
-        exposure : float, default None
+        exposure : astropy.units.Quantity (default: None)
             By default, the exposure stored in the ARF will be used to compute
             the total counts per bin over the effective observation time.
             In cases where this might be incorrect (e.g. for simulated spectra
@@ -149,9 +150,6 @@ def apply_response(self, mflux, exposure=None):
             ARF), this keyword provides the functionality to override the
             default behaviour and manually set the exposure time to use.
 
-        store_model_counts : bool
-            If True, the output will also be stored in self.model_counts
-
         Returns
         -------
         model_counts : numpy.ndarray
@@ -182,7 +180,7 @@ def __init__(self, filename=None, extension=None,
                  offset=0.0, n_grp=np.array([]), f_chan=np.array([]),
                  n_chan=np.array([]), detchans=0):
         """
-        Load a response matrix file (RMF) from a FITS file.
+        The response matrix file (RMF) for an X-ray observation.
 
         filename : str
             The file name with the RMF file
@@ -322,10 +320,6 @@ def _flatten_arrays(self, n_grp, f_chan, n_chan, matrix):
         # stack all non-zero rows in the matrix
         matrix_flat = np.hstack(matrix[nz_idx])
 
-        # stack all nonzero rows in n_chan and f_chan
-        #n_chan_flat = np.hstack(n_chan[nz_idx])
-        #f_chan_flat = np.hstack(f_chan[nz_idx])
-
         # some matrices actually have more elements
         # than groups in `n_grp`, so we'll only pick out
         # those values that have a correspondence in

From c760cf4b3f296eddc9c67b5ca1347f22e34f5730 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Thu, 27 Dec 2018 16:55:08 -0700
Subject: [PATCH 84/89] added basic tests for a fake arf

---
 specutils/tests/test_xrayspectrum1d.py | 25 +++++++++++++++++++++++--
 1 file changed, 23 insertions(+), 2 deletions(-)

diff --git a/specutils/tests/test_xrayspectrum1d.py b/specutils/tests/test_xrayspectrum1d.py
index 6a561307f..60bf8ff7f 100644
--- a/specutils/tests/test_xrayspectrum1d.py
+++ b/specutils/tests/test_xrayspectrum1d.py
@@ -1,7 +1,7 @@
 import numpy as np
 import astropy.units as u
 from astropy.modeling.powerlaws import PowerLaw1D
-from specutils import XraySpectrum1D
+from specutils import XraySpectrum1D, AreaResponse, ResponseMatrix
 
 def test_create_from_arrays():
     # Test that XraySpectrum1D can be initialized
@@ -27,8 +27,29 @@ def test_call_spectrum():
 def test_apply_model():
     # Test that one can evaluate XraySpectrum1D with a model
     test_spec = test_create_from_arrays()
-    new_model = PowerLaw1D(amplitude=5.e-3, alpha=1.8, x_0=1.e3)
+    new_model = PowerLaw1D(amplitude=5.e-3, alpha=1.8, x_0=1.0)
     model_flux = new_model(test_spec.spectral_axis.value)
     ymodel = test_spec.apply_response(model_flux)
     assert len(ymodel) == len(test_spec.spectral_axis)
+    # Since there's no arf or rmf, should get the same array out
+    assert np.all(ymodel == model_flux)
+
+    # Now add an arf
+    specresp = np.ones_like(test_spec.bin_lo.value) * 2.0
+    arf = AreaResponse(test_spec.bin_lo, test_spec.bin_hi,
+                       specresp, test_spec.exposure)
+    test_spec.assign_arf(arf)
+
+    # Test that the arf is applied properly
+    ymodel2 = test_spec.apply_response(model_flux)
+    # The arf also gets multiplied by the exposure time, so that needs
+    # to be included in this check
+    check_value = ymodel2 / (model_flux * test_spec.exposure)
+    assert np.all(check_value == 2.0)
+    return
+
+## Currently there is no RMF check because it requires a test data file
+## However, we can test that it can be initialized with empty values
+def test_init_rmf():
+    rmf = ResponseMatrix()
     return

From 2929b37bbb54af26007607ea4c29f6cd25035d3c Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Fri, 28 Dec 2018 09:47:13 -0700
Subject: [PATCH 85/89] made static read method for ARF and RMF

---
 specutils/spectra/xrayspectrum1d.py | 164 +++++++++++++---------------
 1 file changed, 76 insertions(+), 88 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 736f252f9..0858dc975 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -5,7 +5,7 @@
 from astropy.io import fits
 from .spectrum1d import Spectrum1D
 
-__all__ = ['XraySpectrum1D', 'AreaResponse', 'ResponseMatrix', 'arf_loader']
+__all__ = ['XraySpectrum1D', 'AreaResponse', 'ResponseMatrix']
 
 # For dealing with varied unit string choices
 EV   = ['eV', 'ev']
@@ -15,17 +15,6 @@
 ## An introduction to X-ray analysis can be found here:
 ## http://cxc.cfa.harvard.edu/xrayschool/talks/intro_xray_analysis.pdf
 
-def _unit_parser(unit_string):
-    """
-    A function for matching unit strings to the correct Astropy unit
-    """
-    if unit_string in EV:
-        return u.eV
-    if unit_string in KEV:
-        return u.keV
-    if unit_string in ANGS:
-        return u.angstrom
-
 class XraySpectrum1D(Spectrum1D):
     """
     Spectrum properties specific to X-ray data
@@ -104,7 +93,7 @@ def assign_arf(self, arf_inp, **kwargs):
         Modifies the XraySpectrum1D.arf attribute
         """
         if isinstance(arf_inp, str):
-            self.arf = arf_loader(arf_inp, **kwargs)
+            self.arf = AreaResponse.read(arf_inp, **kwargs)
         else:
             self.arf = arf_inp
 
@@ -122,7 +111,7 @@ def assign_rmf(self, rmf_inp):
         Modifies the XraySpectrum1D.rmf attribute
         """
         if isinstance(rmf_inp, str):
-            self.rmf = ResponseMatrix(rmf_inp)
+            self.rmf = ResponseMatrix.read(rmf_inp)
         else:
             self.rmf = rmf_inp
         return
@@ -175,19 +164,16 @@ def apply_response(self, mflux, exposure=None):
 ## ----  Supporting response file objects
 
 class ResponseMatrix(object):
-    def __init__(self, filename=None, extension=None,
-                 energ_lo=None, energ_hi=None, matrix=None, energ_unit=None,
+    def __init__(self, energ_lo=None, energ_hi=None, matrix=None, energ_unit=None,
                  offset=0.0, n_grp=np.array([]), f_chan=np.array([]),
                  n_chan=np.array([]), detchans=0):
         """
         The response matrix file (RMF) for an X-ray observation.
 
+        Attributes
+        ----------
         filename : str
-            The file name with the RMF file
-
-        extension : str (default None)
-            FITS file extension keyword, if the response matrix is stored
-            under an extension other than "MATRIX" or "SPECRESP MATRIX"
+            Stores the filename used with ResponseMatrix.read
 
         energ_lo : numpy.ndarray
             The lower edges of the energy bins
@@ -220,7 +206,7 @@ def __init__(self, filename=None, extension=None,
         detchans : int
             The number of channels in the detector
         """
-        self.filename = filename
+        self.filename = None
         self.offset = offset
         self.n_grp = n_grp
         self.f_chan = f_chan
@@ -230,10 +216,10 @@ def __init__(self, filename=None, extension=None,
         self.energ_hi = energ_hi
         self.energ_unit = energ_unit
         self.detchans = detchans
-        if filename is not None:
-            self._load_rmf(filename, extension=extension)
 
-    def _load_rmf(self, filename, extension=None):
+    @staticmethod
+    def read(filename, extension=None):
+        result = ResponseMatrix()
         # open the FITS file and extract the MATRIX extension
         # which contains the redistribution matrix and
         # anxillary information
@@ -267,17 +253,17 @@ def _load_rmf(self, filename, extension=None):
         n_chan = np.array(data.field("N_CHAN"))
         matrix = np.array(data.field("MATRIX"))
 
-        self.energ_lo = np.array(data.field("ENERG_LO"))
-        self.energ_hi = np.array(data.field("ENERG_HI"))
-        self.energ_unit = u.Unit(data.columns["ENERG_LO"].unit)
-        self.detchans = hdr["DETCHANS"]
-        self.offset = self.__get_tlmin(h)
+        result.energ_lo = np.array(data.field("ENERG_LO"))
+        result.energ_hi = np.array(data.field("ENERG_HI"))
+        result.energ_unit = u.Unit(data.columns["ENERG_LO"].unit)
+        result.detchans = hdr["DETCHANS"]
+        result.offset = result.__get_tlmin(h)
 
         # flatten the variable-length arrays
-        self.n_grp, self.f_chan, self.n_chan, self.matrix = \
-                self._flatten_arrays(n_grp, f_chan, n_chan, matrix)
+        result.n_grp, result.f_chan, result.n_chan, result.matrix = \
+                result._flatten_arrays(n_grp, f_chan, n_chan, matrix)
 
-        return
+        return result
 
     def __get_tlmin(self, h):
         """
@@ -446,12 +432,15 @@ def apply_rmf(self, model):
 
 class AreaResponse(object):
     def __init__(self, e_low, e_high, eff_area,
-                 exposure=None, fracexpo=1.0, filename=None):
+                 exposure=None, fracexpo=1.0):
         """
         Load an area response file (ARF) from a FITS file.
 
         Attributes
         ----------
+        filename : str
+            Stores the filename used with AreaResponse.read
+
         e_low : astropy.units.Quantity
             The lower edges of the energy bins
 
@@ -471,12 +460,9 @@ def __init__(self, e_low, e_high, eff_area,
             These values are stored for reference; generally, they are already
             accounted for in the eff_area array.
 
-        filename : str
-            Stores the name of the file from which the ARF was loaded.
-
         e_mid : Property that returns the middle of each energy bin
         """
-        self.filename = filename
+        self.filename = None
         self.e_low    = e_low
         self.e_high   = e_high
         self.eff_area = eff_area
@@ -487,6 +473,58 @@ def __init__(self, e_low, e_high, eff_area,
     def e_mid(self):
         return 0.5 * (self.e_low + self.e_high)
 
+    @staticmethod
+    def read(filename, block='SPECRESP'):
+        """
+        Load a AreaResponse object from FITS file.
+
+        Parameters
+        ----------
+        filename : str
+            Path to the FITS file
+
+        block : str
+            FITS file block keyword, if the spectral response is stored
+            under an extension other than "SPECRESP"
+
+        Returns
+        -------
+        AreaResponse
+            The ARF that is represented by the FITS file
+        """
+        # open the FITS file and extract the SPECRESP block
+        # which contains the spectral response for the telescope
+        hdulist = fits.open(filename)
+
+        # Get the data from the appropriate extension
+        blocknames = np.array([h.name for h in hdulist])
+        assert block in blocknames, "Could not fine {} block in FITS file".format(block)
+        h = hdulist[block]
+
+        data = h.data
+        hdr = h.header
+        hdulist.close()
+
+        e_unit = u.Unit(data.columns["ENERG_LO"].unit)
+        e_low  = np.array(data.field("ENERG_LO")) * e_unit
+        e_high = np.array(data.field("ENERG_HI")) * e_unit
+
+        area_unit = u.Unit(data.columns['SPECRESP'].unit) * u.ct # usually ct cm^2 / phot
+        specresp  = np.array(data.field("SPECRESP")) * area_unit
+
+        exposure = None
+        if "EXPOSURE" in list(hdr.keys()):
+            exposure = hdr["EXPOSURE"] * u.second
+
+        fracexpo = 1.0
+        if "FRACEXPO" in data.columns.names:
+            fracexpo = data["FRACEXPO"]
+
+        result = AreaResponse(e_low, e_high, specresp,
+                              exposure=exposure, fracexpo=fracexpo)
+        result.filename = filename
+        return result
+
     def apply_arf(self, model, exposure=None):
         """
         Fold the spectrum through the area response file (ARF).
@@ -520,53 +558,3 @@ def apply_arf(self, model, exposure=None):
             return model * self.eff_area * self.exposure
         else:
             return model * self.eff_area * exposure
-
-## ----- Convenience functions for loading ARFs and RMFs
-def arf_loader(filename, block='SPECRESP'):
-    """
-    Load a AreaResponse object from FITS file.
-
-    Parameters
-    ----------
-    filename : str
-        Path to the FITS file
-
-    block : str
-        FITS file block keyword, if the spectral response is stored
-        under an extension other than "SPECRESP"
-
-    Returns
-    -------
-    AreaResponse
-        The ARF that is represented by the FITS file
-    """
-    # open the FITS file and extract the SPECRESP block
-    # which contains the spectral response for the telescope
-    hdulist = fits.open(filename)
-
-    # Get the data from the appropriate extension
-    blocknames = np.array([h.name for h in hdulist])
-    assert block in blocknames, "Could not fine {} block in FITS file".format(block)
-    h = hdulist[block]
-
-    data = h.data
-    hdr = h.header
-    hdulist.close()
-
-    e_unit = u.Unit(data.columns["ENERG_LO"].unit)
-    e_low  = np.array(data.field("ENERG_LO")) * e_unit
-    e_high = np.array(data.field("ENERG_HI")) * e_unit
-
-    area_unit = u.Unit(data.columns['SPECRESP'].unit) * u.ct # usually ct cm^2 / phot
-    specresp  = np.array(data.field("SPECRESP")) * area_unit
-
-    exposure = None
-    if "EXPOSURE" in list(hdr.keys()):
-        exposure = hdr["EXPOSURE"] * u.second
-
-    fracexpo = 1.0
-    if "FRACEXPO" in data.columns.names:
-        fracexpo = data["FRACEXPO"]
-
-    return AreaResponse(e_low, e_high, specresp,
-                        exposure=exposure, fracexpo=fracexpo, filename=filename)

From 462c7fb285030d6745eb576b37a734695f63db79 Mon Sep 17 00:00:00 2001
From: Lia Corrales <liac@umich.edu>
Date: Fri, 28 Dec 2018 10:31:43 -0700
Subject: [PATCH 86/89] switch back to ARF and RMF convention

---
 specutils/spectra/xrayspectrum1d.py | 43 ++++++++++++++---------------
 1 file changed, 21 insertions(+), 22 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index 0858dc975..ee3c60afb 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -5,7 +5,7 @@
 from astropy.io import fits
 from .spectrum1d import Spectrum1D
 
-__all__ = ['XraySpectrum1D', 'AreaResponse', 'ResponseMatrix']
+__all__ = ['XraySpectrum1D', 'ARF', 'RMF']
 
 # For dealing with varied unit string choices
 EV   = ['eV', 'ev']
@@ -33,12 +33,12 @@ class XraySpectrum1D(Spectrum1D):
     exposure : astropy.Quantity
         Exposure time for the dataset
 
-    arf : specutils.AreaResponse or string, default None
-        Strings will be passed to AreaResponse.__init__
+    arf : specutils.ARF or string, default None
+        Strings will be passed to ARF.__init__
         All other input types are stored as arf attribute
 
-    rmf : specutils.ResponseMatrix or string, default None
-        Strings will be passed to ResponseMatrix.__init__
+    rmf : specutils.RMF or string, default None
+        Strings will be passed to RMF.__init__
         All other input types are stored as rmf attribute
 
     rest_value : astropy.units.Quantity, default 0 Angstrom
@@ -81,7 +81,7 @@ def counts(self):
 
     def assign_arf(self, arf_inp, **kwargs):
         """
-        Assign an area response file (ARF) object to the XraySpectrum1D object
+        Assign an auxiliary response file (ARF) object to the XraySpectrum1D object
 
         Input
         -----
@@ -93,13 +93,13 @@ def assign_arf(self, arf_inp, **kwargs):
         Modifies the XraySpectrum1D.arf attribute
         """
         if isinstance(arf_inp, str):
-            self.arf = AreaResponse.read(arf_inp, **kwargs)
+            self.arf = ARF.read(arf_inp, **kwargs)
         else:
             self.arf = arf_inp
 
     def assign_rmf(self, rmf_inp):
         """
-        Assign a response matrix file (RMF) object to the XraySpectrum1D object
+        Assign a redistribution matrix file (RMF) object to the XraySpectrum1D object
 
         Input
         -----
@@ -111,7 +111,7 @@ def assign_rmf(self, rmf_inp):
         Modifies the XraySpectrum1D.rmf attribute
         """
         if isinstance(rmf_inp, str):
-            self.rmf = ResponseMatrix.read(rmf_inp)
+            self.rmf = RMF.read(rmf_inp)
         else:
             self.rmf = rmf_inp
         return
@@ -119,7 +119,7 @@ def assign_rmf(self, rmf_inp):
     def apply_response(self, mflux, exposure=None):
         """
         Given a model flux spectrum, apply the response. In cases where the
-        spectrum has both an area response file (ARF) and a response matrix
+        spectrum has both an auxiliary response file (ARF) and a redistribution matrix
         file (RMF), apply both. Otherwise, apply whatever response is in RMF.
 
         The model flux spectrum *must* be created using the same units and
@@ -163,17 +163,17 @@ def apply_response(self, mflux, exposure=None):
 
 ## ----  Supporting response file objects
 
-class ResponseMatrix(object):
+class RMF(object):
     def __init__(self, energ_lo=None, energ_hi=None, matrix=None, energ_unit=None,
                  offset=0.0, n_grp=np.array([]), f_chan=np.array([]),
                  n_chan=np.array([]), detchans=0):
         """
-        The response matrix file (RMF) for an X-ray observation.
+        Redistribution Matrix File (RMF) for an X-ray spectrum.
 
         Attributes
         ----------
         filename : str
-            Stores the filename used with ResponseMatrix.read
+            Stores the filename used with RMF.read
 
         energ_lo : numpy.ndarray
             The lower edges of the energy bins
@@ -219,7 +219,7 @@ def __init__(self, energ_lo=None, energ_hi=None, matrix=None, energ_unit=None,
 
     @staticmethod
     def read(filename, extension=None):
-        result = ResponseMatrix()
+        result = RMF()
         # open the FITS file and extract the MATRIX extension
         # which contains the redistribution matrix and
         # anxillary information
@@ -430,16 +430,16 @@ def apply_rmf(self, model):
         return counts[:self.detchans]
 
 
-class AreaResponse(object):
+class ARF(object):
     def __init__(self, e_low, e_high, eff_area,
                  exposure=None, fracexpo=1.0):
         """
-        Load an area response file (ARF) from a FITS file.
+        Auxiliary Response File (ARF) for an X-ray spectrum.
 
         Attributes
         ----------
         filename : str
-            Stores the filename used with AreaResponse.read
+            Stores the filename used with ARF.read
 
         e_low : astropy.units.Quantity
             The lower edges of the energy bins
@@ -476,7 +476,7 @@ def e_mid(self):
     @staticmethod
     def read(filename, block='SPECRESP'):
         """
-        Load a AreaResponse object from FITS file.
+        Load an ARF object from FITS file.
 
         Parameters
         ----------
@@ -489,8 +489,7 @@ def read(filename, block='SPECRESP'):
 
         Returns
         -------
-        AreaResponse
-            The ARF that is represented by the FITS file
+        The ARF object that is represented by the FITS file
         """
         # open the FITS file and extract the SPECRESP block
         # which contains the spectral response for the telescope
@@ -520,14 +519,14 @@ def read(filename, block='SPECRESP'):
         if "FRACEXPO" in data.columns.names:
             fracexpo = data["FRACEXPO"]
 
-        result = AreaResponse(e_low, e_high, specresp,
+        result = ARF(e_low, e_high, specresp,
                               exposure=exposure, fracexpo=fracexpo)
         result.filename = filename
         return result
 
     def apply_arf(self, model, exposure=None):
         """
-        Fold the spectrum through the area response file (ARF).
+        Fold the spectrum through the auxiliary response file (ARF).
         The ARF is a single vector encoding the effective area information
         about the detector. A such, applying the ARF is a simple
         multiplication with the input spectrum.

From b00601c61394dc37d907d0cbf46610d6dc604b58 Mon Sep 17 00:00:00 2001
From: eblur <liac@umich.edu>
Date: Thu, 17 Jan 2019 16:50:10 -0500
Subject: [PATCH 87/89] dont need os module

---
 specutils/io/default_loaders/chandra_hetg.py | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/specutils/io/default_loaders/chandra_hetg.py b/specutils/io/default_loaders/chandra_hetg.py
index 306dab48c..532ada800 100644
--- a/specutils/io/default_loaders/chandra_hetg.py
+++ b/specutils/io/default_loaders/chandra_hetg.py
@@ -1,5 +1,3 @@
-import os
-
 from astropy.io import fits
 from astropy.units import Unit
 
@@ -40,8 +38,6 @@ def hetg_spectrum_loader(file_name, arf=None, rmf=None):
         The spectrum that is represented by the data in this table.
     """
 
-    name = os.path.basename(file_name)
-
     with fits.open(file_name) as hdu:
         header = hdu[0].header
         meta   = {'header': header}

From f137c02b35ad2204ef4deb75f819217f7ad266cd Mon Sep 17 00:00:00 2001
From: eblur <liac@umich.edu>
Date: Fri, 18 Jan 2019 14:45:12 -0500
Subject: [PATCH 88/89] fixed calls to ARF and RMF

---
 specutils/tests/test_xrayspectrum1d.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/specutils/tests/test_xrayspectrum1d.py b/specutils/tests/test_xrayspectrum1d.py
index 60bf8ff7f..fd583fb91 100644
--- a/specutils/tests/test_xrayspectrum1d.py
+++ b/specutils/tests/test_xrayspectrum1d.py
@@ -1,7 +1,7 @@
 import numpy as np
 import astropy.units as u
 from astropy.modeling.powerlaws import PowerLaw1D
-from specutils import XraySpectrum1D, AreaResponse, ResponseMatrix
+from specutils import XraySpectrum1D, ARF, RMF
 
 def test_create_from_arrays():
     # Test that XraySpectrum1D can be initialized
@@ -36,7 +36,7 @@ def test_apply_model():
 
     # Now add an arf
     specresp = np.ones_like(test_spec.bin_lo.value) * 2.0
-    arf = AreaResponse(test_spec.bin_lo, test_spec.bin_hi,
+    arf = ARF(test_spec.bin_lo, test_spec.bin_hi,
                        specresp, test_spec.exposure)
     test_spec.assign_arf(arf)
 
@@ -51,5 +51,5 @@ def test_apply_model():
 ## Currently there is no RMF check because it requires a test data file
 ## However, we can test that it can be initialized with empty values
 def test_init_rmf():
-    rmf = ResponseMatrix()
+    rmf = RMF()
     return

From 4225836daf5fe4551d3cae27fb8ce87ad726c81f Mon Sep 17 00:00:00 2001
From: eblur <liac@umich.edu>
Date: Fri, 18 Jan 2019 21:28:22 -0500
Subject: [PATCH 89/89] removed rest value keyword

---
 specutils/spectra/xrayspectrum1d.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/specutils/spectra/xrayspectrum1d.py b/specutils/spectra/xrayspectrum1d.py
index ee3c60afb..f4090ce5c 100644
--- a/specutils/spectra/xrayspectrum1d.py
+++ b/specutils/spectra/xrayspectrum1d.py
@@ -62,11 +62,10 @@ class XraySpectrum1D(Spectrum1D):
     rmf
     """
     def __init__(self, bin_lo, bin_hi, counts, exposure,
-                 arf=None, rmf=None, rest_value=0.0 * u.angstrom, **kwargs):
+                 arf=None, rmf=None, **kwargs):
 
         bin_mid = 0.5 * (bin_lo + bin_hi)
-        Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts,
-                            rest_value=rest_value, **kwargs)
+        Spectrum1D.__init__(self, spectral_axis=bin_mid, flux=counts, **kwargs)
 
         self.bin_lo = bin_lo
         self.bin_hi = bin_hi