Merge e79a403 into afe790c

MANIFEST.in

@@ -1,15 +1,13 @@
 include WDmodel/TlustyGrids.hdf5
 include WDmodel/WDmodel_param_defaults.json
+include WDmodel/WDmodel_pb_obsmode_map.txt
+include WDmodel/reddening_custom.txt
+graft WDmodel/tests/
+graft WDmodel/passbands/
 include tesh.sh
 include test_WDmodel.py
 include LICENSE
-exclude docs/_build
-exclude docs/api
-exclude build/
-exclude cache/
-exclude products/
-exclude dist/
-exclude src/
+graft docs/
 global-include *.txt
 global-include *.rst
 global-include *.md
@@ -20,8 +18,13 @@ global-exclude *.prof
 global-exclude *.swp
 global-exclude *.png
 global-exclude out*/*
-graft docs/
-graft  WDmodel/tests/
+exclude docs/_build
+exclude docs/api
+exclude build/
+exclude cache/
+exclude products/
+exclude dist/
+exclude src/
 prune scripts/
 prune data/
 prune WDdata/

WDmodel/WDmodel_pb_obsmode_map.txt

@@ -1,7 +1,20 @@
-#pb obsmode magsys
-F275W WFC3,UVIS2,F275W abmag
-F336W WFC3,UVIS2,F336W abmag
-F475W WFC3,UVIS2,F475W abmag
-F625W WFC3,UVIS2,F625W abmag
-F775W WFC3,UVIS2,F775W abmag
-F160W WFC3,IR,F160W abmag
+#pb obsmode magsys magzero
+F275W WFC3,UVIS2,F275W abmag 0.
+F336W WFC3,UVIS2,F336W abmag 0.
+F475W WFC3,UVIS2,F475W abmag 0.
+F625W WFC3,UVIS2,F625W abmag 0.
+F775W WFC3,UVIS2,F775W abmag 0.
+F160W WFC3,IR,F160W abmag 0.
+ps1g passbands/PS1/g_filt_revised.txt abmag 0.
+ps1r passbands/PS1/r_filt_tonry.txt abmag 0.
+ps1i passbands/PS1/i_filt_tonry.txt abmag 0.
+ps1z passbands/PS1/z_filt_tonry.txt abmag 0.
+ps1y passbands/PS1/y_filt_tonry.txt abmag 0.
+sdssu passbands/SDSS/sdss_u.dat abmag 0.
+sdssg passbands/SDSS/sdss_g.dat abmag 0.
+sdssr passbands/SDSS/sdss_r.dat abmag 0.
+sdssi passbands/SDSS/sdss_i.dat abmag 0.
+sdssz passbands/SDSS/sdss_z.dat abmag 0.
+gaia_g passbands/Gaia/gaia_g.dat vegamag 0.
+gaia_bp passbands/Gaia/gaia_bp.dat vegamag 0.
+gaia_rp passbands/Gaia/gaia_rp.dat vegamag 0.

WDmodel/io.py

@@ -997,14 +997,48 @@ def get_pkgfile(infile):
         model grid file.
     """
 
-    pkgfile = pkg_resources.resource_filename('WDmodel',infile)
+    try:
+        pkgfile = pkg_resources.resource_filename('WDmodel',infile)
+    except KeyError as e:
+        message = 'Could not find package file {}'.format(infile)
+        raise IOError(message)
 
     if not os.path.exists(pkgfile):
         message = 'Could not find package file {}'.format(pkgfile)
         raise IOError(message)
     return pkgfile
 
 
+def get_filepath(infile):
+    """
+    Returns the full path to a file. If the path is relative, it is converted
+    to absolute. If this file does not exist, it is treated as a file within
+    the :py:mod:`WDmodel` package. If that file does not exist, an error is
+    raised.
+
+    Parameters
+    ----------
+    infile : str
+        The name of the file to set the full path for
+
+    Returns
+    -------
+    pkgfile : str
+        The path to the file
+
+    Raises
+    ------
+    IOError
+        If the ``infile`` could not be found at location or inside the
+        :py:mod:`WDmodel` package.
+    """
+    fullfile = os.path.abspath(infile)
+    if os.path.exists(fullfile):
+        return fullfile
+    else:
+        return get_pkgfile(infile)
+
+
 def write_fit_inputs(spec, phot, cont_model, linedata, continuumdata,\
         rvmodel, covtype, coveps, phot_dispersion, scale_factor, outfile):
     """

WDmodel/passband.py

@@ -6,7 +6,9 @@
 
 from __future__ import absolute_import
 from __future__ import unicode_literals
+import warnings
 import numpy as np
+from scipy.interpolate import interp1d
 import pysynphot as S
 from . import io
 from collections import OrderedDict
@@ -226,7 +228,7 @@ def chop_syn_spec_pb(spec, model_mag, pb, model):
     return outpb, outzp
 
 
-def get_pbmodel(pbnames, model, pbfile=None, mag_type=None):
+def get_pbmodel(pbnames, model, pbfile=None, mag_type=None, mag_zero=0.):
     """
     Converts passband names ``pbnames`` into passband models based on the
     mapping of name to ``pysynphot`` ``obsmode`` strings in ``pbfile``.
@@ -250,11 +252,17 @@ def get_pbmodel(pbnames, model, pbfile=None, mag_type=None):
         fullpath to a file that is readable by ``pysynphot``
     mag_type : str, optional
         One of ''vegamag'' or ''abmag''
-        Used to specify the standard that has 0 magnitude in the passband.
+        Used to specify the standard that has mag_zero magnitude in the passband.
         If ``magsys`` is specified in ``pbfile,`` that overrides this option.
         Must be the same for all passbands listed in ``pbname`` that do not
         have ``magsys`` specified in ``pbfile``
-        If ``pbnames`` require multiple ``mag_types``, concatentate the output.
+        If ``pbnames`` require multiple ``mag_type``, concatentate the output.
+    mag_zero : float, optional
+        Magnitude of the standard in the passband
+        If ``magzero`` is specified in ``pbfile,`` that overrides this option.
+        Must be the same for all passbands listed in ``pbname`` that do not
+        have ``magzero`` specified in ``pbfile``
+        If ``pbnames`` require multiple ``mag_zero``, concatentate the output.
 
     Returns
     -------
@@ -290,7 +298,7 @@ def get_pbmodel(pbnames, model, pbfile=None, mag_type=None):
 
         Trims the bandpass to entries with non-zero transmission and determines
         the ``VEGAMAG/ABMAG`` zeropoint for the passband - i.e. ``zp`` that
-        gives ``mag_Vega/AB=0.`` in all passbands.
+        gives ``mag_Vega/AB=mag_zero`` in all passbands.
 
     See Also
     --------
@@ -305,19 +313,27 @@ def get_pbmodel(pbnames, model, pbfile=None, mag_type=None):
     pbdata  = io.read_pbmap(pbfile)
     pbmap   = dict(list(zip(pbdata.pb, pbdata.obsmode)))
     sysmap  = dict(list(zip(pbdata.pb, pbdata.magsys)))
+    zeromap = dict(list(zip(pbdata.pb, pbdata.magzero)))
 
     # setup the photometric system by defining the standard and corresponding magnitude system
     if mag_type not in ('vegamag', 'abmag', None):
         message = 'Magnitude system must be one of abmag or vegamag'
         raise RuntimeError(message)
 
+    try:
+        mag_zero = float(mag_zero)
+    except ValueError as e:
+        message = 'Zero magnitude must be a floating point number'
+        raise RuntimeError(message)
+
     # define the standards
     vega = S.Vega
-    ab   = S.FlatSpectrum(3631, waveunits='angstrom', fluxunits='jy')
+    vega.convert('flam')
+    ab   = S.FlatSpectrum(0., waveunits='angstrom', fluxunits='abmag')
     ab.convert('flam')
 
     # defile the magnitude sysem
-    if mag_type is None or mag_type == 'vegamag':
+    if mag_type == 'vegamag':
         mag_type= 'vegamag'
     else:
         mag_type = 'abmag'
@@ -331,6 +347,7 @@ def get_pbmodel(pbnames, model, pbfile=None, mag_type=None):
         # load each passband
         obsmode = pbmap.get(pb, pb)
         magsys  = sysmap.get(pb, mag_type)
+        synphot_mag = zeromap.get(pb, mag_zero)
 
         if magsys == 'vegamag':
             standard = vega
@@ -340,24 +357,44 @@ def get_pbmodel(pbnames, model, pbfile=None, mag_type=None):
             message = 'Unknown standard system {} for passband {}'.format(magsys, pb)
             raise RuntimeError(message)
 
+        loadedpb = False
         # treat the passband as a obsmode string
         try:
             bp = S.ObsBandpass(obsmode)
+            loadedpb = True
         except ValueError:
-            # if that fails, try to load the passband interpreting obsmode as a file
             message = 'Could not load pb {} as an obsmode string {}'.format(pb, obsmode)
-            print(message)
+            warnings.warn(message, RuntimeWarning)
+            loadedpb = False
+
+        # if that fails, try to load the passband interpreting obsmode as a file
+        if not loadedpb:
             try:
-                bp = S.FileBandpass(obsmode)
-            except ValueError:
+                bandpassfile = io.get_filepath(obsmode)
+                bp = S.FileBandpass(bandpassfile)
+                loadedpb = True
+            except Exception as e:
                 message = 'Could not load passband {} from obsmode or file {}'.format(pb, obsmode)
-                raise RuntimeError(message)
+                warnings.warn(message, RuntimeWarning)
+                loadedpb = False
 
-        avgwave = bp.avgwave()
+        if not loadedpb:
+            message = 'Could not load passband {}. Giving up.'.format(pb)
+            raise RuntimeError(message)
 
-        # get the pysynphot standard magnitude (should be 0. on the standard magnitude system!)
-        ob = S.Observation(standard, bp)
-        synphot_mag = ob.effstim(magsys)
+        avgwave = bp.avgwave()
+        if standard.wave.min() > model._wave.min():
+            message = 'Standard does not extend past the blue edge of the model'
+            warnings.warn(message, RuntimeWarning)
+
+        if standard.wave.max() < model._wave.max():
+            message = 'Standard does not extend past the red edge of the model'
+            warnings.warn(message, RuntimeWarning)
+
+        # interpolate the standard onto the model wavelengths
+        sinterp = interp1d(standard.wave, standard.flux, fill_value='extrapolate')
+        standard_flux = sinterp(model._wave)
+        standard = np.rec.fromarrays([model._wave, standard_flux], names='wave,flux')
 
         # cut the passband to non-zero values and interpolate onto overlapping standard wavelengths
         outpb, outzp = chop_syn_spec_pb(standard, synphot_mag, bp, model)