Merge pull request #47 from pllim/fix-ext-waveset

Ensure extinction waveset is not carried into composite spectrum

pysynphot/reddening.py

@@ -8,6 +8,29 @@
 from .spectrum import ArraySpectralElement
 from . import Cache
 from . import extinction #temporary(?) backwards compatibility
+from . import units
+
+
+# https://github.com/spacetelescope/pysynphot/issues/44
+class ExtinctionSpectralElement(ArraySpectralElement):
+    """Like :class:`~pysynphot.spectrum.ArraySpectralElement` but
+    with special ``waveset`` handling.
+    """
+    def GetWaveSet(self):
+        """Extinction curve ``waveset`` should not be merged."""
+        return None
+
+    def _getWaveProp(self):
+        """Return wavelength in user units."""
+        wave = units.Angstrom().Convert(self._wavetable, self.waveunits.name)
+        return wave
+
+    wave = property(_getWaveProp, doc="Wavelength property.")
+
+    def GetThroughput(self):
+        return self.__call__(self._wavetable)
+
+    throughput = property(GetThroughput, doc='Throughput property.')
 
 
 class CustomRedLaw(object):
@@ -78,11 +101,10 @@ def reddening(self,extval):
 
         """
         T = 10.0**(-0.4*extval*self.obscuration)
-        ans=ArraySpectralElement(wave=self.wave,
-                                 waveunits=self.waveunits,
-                                 throughput=T,
-                                 name='%s(Av=%g)'%(self.name,extval)
-                                 )
+        ans = ExtinctionSpectralElement(wave=self.wave,
+                                        waveunits=self.waveunits,
+                                        throughput=T,
+                                        name='%s(Av=%g)'%(self.name, extval))
         ans.citation = self.litref
         return ans
 

pysynphot/test/test_mergewavesets.py

@@ -1,10 +1,15 @@
 from __future__ import absolute_import, division, print_function
 
+import os
+
 import numpy as np
+from astropy.utils.data import get_pkg_data_filename
 
 from .utils import use_cdbs
+from ..obsbandpass import ObsBandpass
+from ..observation import Observation
 from ..reddening import Extinction
-from ..spectrum import BlackBody, MergeWaveSets, MERGETHRESH
+from ..spectrum import BlackBody, FileSourceSpectrum, MergeWaveSets, MERGETHRESH
 
 
 @use_cdbs
@@ -20,3 +25,21 @@ def test_merge_wave_sets():
     delta = new_wave[1:] - new_wave[:-1]
     assert np.all(delta > MERGETHRESH), \
         'Deltas should be < {}, min delta = {}'.format(MERGETHRESH, delta.min())  # noqa
+
+
+@use_cdbs
+def test_qso_countrate():
+    """
+    Extinction curve waveset should not be merged into composite spectrum
+    when applied. See https://github.com/spacetelescope/pysynphot/issues/44 .
+    """
+    bp = ObsBandpass('acs,hrc,f850lp')
+
+    fname = get_pkg_data_filename(os.path.join('data', 'qso_template.fits'))
+    qso = FileSourceSpectrum(fname)
+    sp_ext = qso * Extinction(1.0, 'mwavg')
+    sp = sp_ext.renorm(20, 'vegamag', ObsBandpass('johnson,v'), force=True)
+
+    obs = Observation(sp, bp, force='taper')
+    c = obs.countrate()
+    np.testing.assert_allclose(c, 2.3554364232173565e-05)