Adding package

arcesetc/__init__.py

@@ -20,4 +20,5 @@ class UnsupportedPythonError(Exception):
 
 if not _ASTROPY_SETUP_:
     # For egg_info test builds to pass, put package imports here.
-    pass
+    from .utils import *
+    from .plots import *

arcesetc/data/sptype_dict.json

@@ -0,0 +1 @@
+{"G2IV": "51 Peg", "K4V": "HD 79555", "K7": "GJ 9781A", "K0": "HD 145742", "G0V": "HD39587", "G3VaHdel1": "HD86728", "B6IV": "HR 6943", "K0Ve": "HD87884", "K6VeFe-1": "HD 88230", "G8V": "HD62613", "G8IVk": "HD 222107", "K2": "GJ 705", "M2.0Ve": "GJ4099", "K5Ve": "EQ Vir", "G5V:": "HD134319", "K2III-IV": "HD 6497", "G7V": "HD 67767", "K3V": "HR 8832", "G1.5V": "HD 129333", "G0Vs": "KIC9652680", "G2V": "HD98230", "G3Va": "HD 10697", "K5V": "HD149957", "K7.5Ve": "HD151288", "K7V": "61CygB", "G5": "Kepler-63", "K0V": "sigma Draconis", "K2V": "HD 73667", "K2.5V": "HD200560", "G3V": "HD68017", "G5V": "Kepler-17", "G9V": "HD 42250", "K1V": "HD221639", "M0.5V": "HD 209290", "WN8h": "WR124", "K3.5V": "HD47752", "sdO2VIIIHe5": "BD+28 4211", "F7V": "HD 89744", "K6V": "HD148467", "G1.5IV-VFe-1": "HD34411", "B3V": "HR 3454", "B9.5V": "HR5501", "K2Ve": "HD45088"}

arcesetc/data/sptype_to_temp.json

@@ -0,0 +1 @@
+{"O5V": 54000, "O6V": 45000, "O7V": 43300, "O8V": 40600, "O9V": 37800, "B0V": 29200, "B1V": 23000, "B2V": 21000, "B3V": 17600, "B5V": 15200, "B6V": 14300, "B7V": 13500, "B8V": 12300, "B9V": 11400, "M7V": 2650, "M7.5V": 2600, "G0V": 5920, "M8V": 2500, "K0V": 5280, "G2V": 5770, "M2.5V": 3500, "F4V": 6640, "A4V": 8270, "M3.5V": 3250, "F1V": 7030, "F5V": 6510, "M4V": 3200, "A1V": 9200, "K4V": 4620, "K1V": 5170, "F3V": 6720, "M0V": 3870, "A5V": 8080, "F2V": 6810, "K6V": 4230, "G5V": 5660, "G1V": 5880, "M1V": 3700, "G9V": 5340, "G8V": 5490, "F6V": 6340, "M6.5V": 2710, "K3V": 4840, "A8V": 7500, "A3V": 8550, "M1.5V": 3650, "K9V": 3940, "A2V": 8840, "M0.5V": 3800, "F9V": 6040, "M6V": 2850, "K2V": 5040, "G3V": 5720, "K7V": 4070, "K5V": 4410, "M4.5V": 3100, "G6V": 5590, "M2V": 3550, "A6V": 8000, "M5.5V": 3000, "M5V": 3030, "F8V": 6170, "M3V": 3410, "F0V": 7220, "G7V": 5530, "A0V": 9700, "A7V": 7800, "A9V": 7440, "F7V": 6240, "G4V": 5680, "K8V": 4000, "WN2": 141000, "WN3": 85000, "WN4": 70000, "WN5": 60000, "WN5h": 50000, "WN6": 56000, "WN6h": 45000, "WN7": 50000, "WN7h" : 45000, "WN8h" : 40000, "WN9h" : 35000, "sdO2VIIIHe5": 82000}

arcesetc/plots.py

@@ -0,0 +1,151 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import astropy.units as u
+
+from .utils import (_closest_target, available_sptypes, archive,
+                    _get_closest_order, _matrix_row_to_spectrum, _scale_flux,
+                    _sn_to_exp_time)
+
+__all__ = ['plot_order_counts', 'plot_order_sn', 'available_sptypes']
+
+
+@u.quantity_input(exp_time=u.s, wavelength=u.Angstrom)
+def plot_order_counts(sptype, wavelength, V, exp_time=None,
+                      signal_to_noise=None, **kwargs):
+    """
+    Plot the counts as a function of wavelength for the spectral
+    order nearest to ``wavelength`` for a star of spectral type ``sptype`` and
+    V magnitude ``V``.
+
+    Either ``exp_time`` or ``signal_to_noise`` should be supplied to the
+    function (but not both).
+
+    Parameters
+    ----------
+    sptype : str
+        Spectral type of the star. If
+    wavelength : `~astropy.units.Quantity`
+    V : float
+        V magnitude of the target
+    exp_time : None or float
+        If ``exp_time`` is a float, show the counts curve for that exposure
+        time. Otherwise, use ``signal_to_noise`` to compute the appropriate
+        exposure time.
+    signal_to_noise : None or float
+        If ``signal_to_noise`` is a float, compute the appropriate exposure time
+        to generate the counts curve that has S/N = ``signal_to_noise`` at
+        wavelength ``wavelength``. Otherwise, generate counts curve for
+        exposure time ``exp_time``.
+    kwargs : dict
+        All extra keyword arguments will be passed to the plot function
+
+    Returns
+    -------
+    fig : `~matplotlib.pyplot.Figure`
+        Matplotlib figure object.
+    ax : `~matplotlib.pyplot.Axes`
+        Matplotlib axes object
+    exp_time : `~astropy.units.Quantity`
+        Exposure time input, or computed to achieve S/N ratio
+        ``signal_to_noise`` at wavelength ``wavelength``
+    """
+    target, closest_spectral_type = _closest_target(sptype)
+
+    matrix = archive[target][:]
+
+    closest_order = _get_closest_order(matrix, wavelength)
+    wave, flux = _matrix_row_to_spectrum(matrix, closest_order)
+    flux *= _scale_flux(archive[target], V)
+
+    fig, ax = plt.subplots()
+
+    if exp_time is not None and signal_to_noise is None:
+        flux *= exp_time.to(u.s).value
+
+    elif exp_time is None and signal_to_noise is not None:
+        exp_time = _sn_to_exp_time(wave, flux, wavelength, signal_to_noise)
+        flux *= exp_time.value
+    else:
+        raise ValueError("Supply either the `exp_time` or the "
+                         "`signal_to_noise` keyword argument.")
+
+    ax.set_title('Sp. Type: {0}, Exposure time: {1:.1f}'
+                 .format(closest_spectral_type, exp_time.to(u.min)))
+    ax.plot(wave, flux, **kwargs)
+    ax.set_xlabel('Wavelength [Angstrom]')
+    ax.set_ylabel('Flux [DN]')
+    for s in ['right', 'top']:
+        ax.spines[s].set_visible(False)
+    ax.grid(ls=':', color='silver')
+    return fig, ax, exp_time
+
+
+@u.quantity_input(exp_time=u.s, wavelength=u.Angstrom)
+def plot_order_sn(sptype, wavelength, V, exp_time=None, signal_to_noise=None,
+                  **kwargs):
+    """
+    Plot the signal-to-noise ratio as a function of wavelength for the spectral
+    order nearest to ``wavelength`` for a star of spectral type ``sptype`` and
+    V magnitude ``V``.
+
+    Either ``exp_time`` or ``signal_to_noise`` should be supplied to the
+    function (but not both).
+
+    Parameters
+    ----------
+    sptype : str
+        Spectral type of the star. If
+    wavelength : `~astropy.units.Quantity`
+    V : float
+        V magnitude of the target
+    exp_time : None or float
+        If ``exp_time`` is a float, show the S/N curve for that exposure time.
+        Otherwise, use ``signal_to_noise`` to compute the appropriate exposure
+        time.
+    signal_to_noise : None or float
+        If ``signal_to_noise`` is a float, compute the appropriate exposure time
+        to generate the S/N curve that has S/N = ``signal_to_noise`` at
+        wavelength ``wavelength``. Otherwise, generate S/N curve for
+        exposure time ``exp_time``.
+    kwargs : dict
+        All extra keyword arguments will be passed to the plot function
+
+    Returns
+    -------
+    fig : `~matplotlib.pyplot.Figure`
+        Matplotlib figure object.
+    ax : `~matplotlib.pyplot.Axes`
+        Matplotlib axes object
+    exp_time : `~astropy.units.Quantity`
+        Exposure time input, or computed to achieve S/N ratio
+        ``signal_to_noise`` at wavelength ``wavelength``
+    """
+    target, closest_spectral_type = _closest_target(sptype)
+
+    matrix = archive[target][:]
+
+    closest_order = _get_closest_order(matrix, wavelength)
+    wave, flux = _matrix_row_to_spectrum(matrix, closest_order)
+    flux *= _scale_flux(archive[target], V)
+
+    fig, ax = plt.subplots()
+
+    if exp_time is not None:
+        flux *= exp_time.to(u.s).value
+    elif exp_time is None and signal_to_noise is not None:
+        exp_time = _sn_to_exp_time(wave, flux, wavelength, signal_to_noise)
+        flux *= exp_time.value
+    else:
+        raise ValueError("Supply either the `exp_time` or the "
+                         "`signal_to_noise` keyword argument.")
+    sn = flux / np.sqrt(flux)
+    ax.set_title('Sp. Type: {0}, Exposure time: {1:.1f}'
+                 .format(closest_spectral_type, exp_time.to(u.min)))
+    ax.plot(wave, sn, **kwargs)
+    ax.set_xlabel('Wavelength [Angstrom]')
+    ax.set_ylabel('Signal/Noise')
+    for s in ['right', 'top']:
+        ax.spines[s].set_visible(False)
+    ax.grid(ls=':', color='silver')
+    return fig, ax, exp_time
+

arcesetc/utils.py

@@ -0,0 +1,194 @@
+
+import numpy as np
+from json import load
+import os
+import h5py
+import astropy.units as u
+
+__all__ = ['available_sptypes', 'signal_to_noise_to_exp_time']
+
+directory = os.path.dirname(__file__)
+
+sptypes = load(open(os.path.join(directory, 'data', 'sptype_dict.json'), 'r'))
+archive = h5py.File(os.path.join(directory, 'data', 'archive.hdf5'), 'r')
+
+sptype_to_temp = load(open(os.path.join(directory, 'data',
+                                        'sptype_to_temp.json'), 'r'))
+spectral_types = [key for key in sptype_to_temp.keys() if key in sptypes]
+temps = np.array([sptype_to_temp[key] for key in spectral_types])
+
+
+def _closest_sptype(sptype):
+    """
+    Return closest spectral type in the archive.
+
+    Parameters
+    ----------
+    sptype : str
+        Spectral type in the format: ``G2V``
+
+    Returns
+    -------
+    closest_spectral_type : str
+        Closest spectral type available in the archive
+    """
+    return spectral_types[np.argmin(np.abs(sptype_to_temp[sptype] - temps))]
+
+
+def _closest_target(sptype):
+    """
+    Return target with the closest spectral type in the archive.
+
+    Parameters
+    ----------
+    sptype : str
+        Spectral type in the format: ``G2V``
+
+    Returns
+    -------
+    target_name : str
+        Name of the target closest to spectral type ``sptype``
+    closest_spectral_type : str
+        Closest spectral type available in the archive
+    """
+    closest_spectral_type = _closest_sptype(sptype)
+    return sptypes[closest_spectral_type], closest_spectral_type
+
+
+def available_sptypes():
+    """
+    Return a list of available spectral types in the archive.
+
+    Returns
+    -------
+    sptypes : list
+        List of available spectral types
+    """
+    return sorted(spectral_types)
+
+
+def _get_closest_order(matrix, wavelength):
+    """
+    Return the spectral order index closest to wavelength ``wavelength``.
+
+    Parameters
+    ----------
+    matrix : `~np.ndarray`
+        Matrix of blaze function curves from the archive
+    wavelength : `~astropy.units.Quantity`
+        Wavelength of interest.
+
+    Returns
+    -------
+    closest_order : int
+        Closest spectral order to wavelength ``wavelength``
+    """
+    return np.argmin(np.abs(matrix[:, 0] - wavelength.to(u.Angstrom).value))
+
+
+def _matrix_row_to_spectrum(matrix, closest_order):
+    """
+    Given a ``matrix`` from the archive and a spectral order index
+    ``closest_order``, return the spectrum (wavelength and flux).
+
+    Parameters
+    ----------
+    matrix : `~np.ndarray`
+        Matrix of blaze function curves from the archive
+    closest_order : int
+        Closest spectral order to wavelength ``wavelength``
+
+    Returns
+    -------
+    wave : `~astropy.units.Quantity`
+        Wavelengths
+    flux : `~np.ndarray`
+        Fluxes in counts per second at each wavelength
+    """
+    lam_0, delta_lam, n_lam = matrix[closest_order][:3]
+    polynomial_coeffs = matrix[closest_order][3:]
+    wave = np.arange(lam_0 - n_lam*delta_lam/2, lam_0 + n_lam*delta_lam/2,
+                     delta_lam) * u.Angstrom
+    flux = np.polyval(polynomial_coeffs, wave-lam_0 * u.Angstrom)
+    return wave, flux
+
+
+def _scale_flux(dataset, V):
+    """
+    Parameters
+    ----------
+    dataset : `~h5py.File.dataset`
+        h5py dataset of the form ``archive[target]``.
+    V : float
+        V magnitude of the target of interest.
+    """
+    template_vmag = dataset.attrs['V'][0]
+    magnitude_scaling = 10**(0.4 * (template_vmag - V))
+    return magnitude_scaling
+
+
+def _sn_to_exp_time(wave, flux, wavelength, signal_to_noise):
+    """
+    Calculate the required exposure time to achieve signal-to-noise ratio
+    ``signal_to_noise`` given the count rates ``flux`` as a function of
+    wavelength ``wave``.
+
+    Parameters
+    ----------
+    wave : `~astropy.units.Quantity`
+        Wavelengths
+    flux : `~np.ndarray`
+        Flux in counts per second
+    wavelength : `~astropy.units.Quantity`
+        Wavelength of interest at which the S/N is ``signal_to_noise``
+    signal_to_noise : float
+        Desired signal-to-noise
+
+    Returns
+    -------
+    exp_time : `~astropy.units.Quantity`
+        Exposure time that will yield signal-to-noise ``signal_to_noise`` at
+        wavelength ``wavelength``.
+    """
+
+    # `flux` at the test wavelength
+    flux_0 = flux[np.argmin(np.abs(wave - wavelength))]
+    exp_time = signal_to_noise**2 / flux_0
+    return exp_time * u.s
+
+
+@u.quantity_input(exp_time=u.s, wavelength=u.Angstrom)
+def signal_to_noise_to_exp_time(sptype, wavelength, V, signal_to_noise):
+    """
+    Compute the exposure time required to collect signal-to-noise ratio
+    ``signal_to_noise`` at wavelength ``wavelength`` for a star of spectral type
+    ``sptype`` and V magnitude ``V``.
+
+    Parameters
+    ----------
+    sptype : str
+        Spectral type of the star. If
+    wavelength : `~astropy.units.Quantity`
+    V : float
+        V magnitude of the target
+    signal_to_noise : float
+        If ``signal_to_noise`` is a float, compute the appropriate exposure time
+        to generate the S/N curve that has S/N = ``signal_to_noise`` at
+        wavelength ``wavelength``. Otherwise, generate S/N curve for
+        exposure time ``exp_time``.
+    Returns
+    -------
+    exp_time : `~astropy.units.Quantity`
+        Exposure time input, or computed to achieve S/N ratio
+        ``signal_to_noise`` at wavelength ``wavelength``
+    """
+    target, closest_spectral_type = _closest_target(sptype)
+
+    matrix = archive[target][:]
+
+    closest_order = _get_closest_order(matrix, wavelength)
+    wave, flux = _matrix_row_to_spectrum(matrix, closest_order)
+    flux *= _scale_flux(archive[target], V)
+
+    exp_time = _sn_to_exp_time(wave, flux, wavelength, signal_to_noise)
+    return exp_time