|
| 1 | +# Licensed under a 3-clause BSD style license - see LICENSE.rst |
| 2 | +from __future__ import (absolute_import, division, print_function, |
| 3 | + unicode_literals) |
| 4 | + |
| 5 | +# Third-party |
| 6 | +import numpy as np |
| 7 | +from astropy import units as u |
| 8 | + |
| 9 | + |
| 10 | +__all__ = ['exptime_from_ccd_snr'] |
| 11 | + |
| 12 | + |
| 13 | +@u.quantity_input(waveset=u.angstrom, flux=u.erg / u.s / u.cm ** 2 / u.cm, |
| 14 | + npix=u.pixel, n_background=u.pixel, |
| 15 | + background_rate=u.ct / u.pixel / u.s, |
| 16 | + darkcurrent_rate=u.ct / u.pixel / u.s) |
| 17 | +def exptime_from_ccd_snr(snr, waveset, flux, observer, telescope, |
| 18 | + npix=1 * u.pixel, |
| 19 | + n_background=np.inf * u.pixel, |
| 20 | + background_rate=0 * (u.ct / u.pixel / u.s), |
| 21 | + darkcurrent_rate=0 * (u.ct / u.pixel / u.s), |
| 22 | + force_overlap='taper'): |
| 23 | + """ |
| 24 | + Returns the exposure time needed in units of seconds to achieve |
| 25 | + the specified (idealized theoretical) signal to noise ratio |
| 26 | + (from pg 57-58 of [1]_). |
| 27 | +
|
| 28 | + Parameters |
| 29 | + ---------- |
| 30 | + snr : float, int, or `~astropy.units.Quantity` |
| 31 | + The signal to noise ratio of the given observation in dimensionless |
| 32 | + units. |
| 33 | + waveset : array-like |
| 34 | + The wavelengths associated with the target's flux. |
| 35 | + flux : array-like |
| 36 | + The flux of the target. |
| 37 | + observer : `~astroplan.observer.Observer` |
| 38 | + The Observer object. |
| 39 | + telescope : `~astroplan.telescope.Telescope` |
| 40 | + The Telescope object. |
| 41 | + npix : `~astropy.units.Quantity`, optional |
| 42 | + Number of pixels under consideration for the signal with units of |
| 43 | + pixels. Default is 1 * astropy.units.pixel. |
| 44 | + n_background : `~astropy.units.Quantity`, optional |
| 45 | + Number of pixels used in the background estimation with units of |
| 46 | + pixels. Default is set to np.inf * astropy.units.pixel |
| 47 | + such that there is no contribution of error due to background |
| 48 | + estimation. This assumes that n_background will be >> npix. |
| 49 | + background_rate : `~astropy.units.Quantity`, optional |
| 50 | + Photons per pixel per second due to the backround/sky with units of |
| 51 | + counts/second/pixel. |
| 52 | + Default is 0 * (astropy.units.ct / |
| 53 | + astropy.units.second / astropy.units.pixel) |
| 54 | + darkcurrent_rate : `~astropy.units.Quantity`, optional |
| 55 | + Counts per pixel per second due to the dark current with units |
| 56 | + of counts/second/pixel. |
| 57 | + Default is 0 * (astropy.units.ct / |
| 58 | + astropy.units.second / astropy.units.pixel) |
| 59 | + force_overlap : {'taper', 'extrap', 'none', None}, optional |
| 60 | + Force the spectral element x source spectrum convolution by the |
| 61 | + specified method even when they don't fully overlap in wavelength. |
| 62 | + 'taper' forces the incomplete spectral component to zero for its |
| 63 | + missing wavelengths, while 'extrap' attempts to extrapolate the |
| 64 | + incomplete part of the spectrum. |
| 65 | + Default is 'taper' |
| 66 | +
|
| 67 | + References |
| 68 | + ---------- |
| 69 | + .. [1] Howell, S. B. 2000, *Handbook of CCD Astronomy* (Cambridge, UK: |
| 70 | + Cambridge University Press) |
| 71 | +
|
| 72 | + Returns |
| 73 | + ------- |
| 74 | + t : `~astropy.units.Quantity` |
| 75 | + The exposure time needed (in seconds) to achieve the given signal |
| 76 | + to noise ratio. |
| 77 | + """ |
| 78 | + # import synphot, which isn't a required package of astroplan |
| 79 | + from synphot.models import Empirical1D |
| 80 | + from synphot.spectrum import SourceSpectrum, SpectralElement |
| 81 | + from synphot.observation import Observation |
| 82 | + |
| 83 | + # set the source spectrum with synphot |
| 84 | + source_spec = SourceSpectrum(Empirical1D, points=waveset, |
| 85 | + lookup_table=flux) |
| 86 | + |
| 87 | + # set the spectral elements if given (quantum efficiency,skymodel,bandpass) |
| 88 | + qe = telescope.ccd_response |
| 89 | + skymodel = observer.skymodel |
| 90 | + |
| 91 | + spec_elements = _get_spectral_element(telescope.bandpass, |
| 92 | + SpectralElement, Empirical1D) |
| 93 | + if skymodel is not False: |
| 94 | + spec_elements *= _get_spectral_element(skymodel, |
| 95 | + SpectralElement, Empirical1D) |
| 96 | + if qe is not False: |
| 97 | + spec_elements *= _get_spectral_element(qe, |
| 98 | + SpectralElement, Empirical1D) |
| 99 | + |
| 100 | + # get the synphot observation object |
| 101 | + synphot_obs = Observation(source_spec, spec_elements, force=force_overlap) |
| 102 | + |
| 103 | + # make sure the gain is in the correct units |
| 104 | + gain = telescope.gain |
| 105 | + if gain.unit in (u.electron / u.adu, u.photon / u.adu): |
| 106 | + gain = gain.value * (u.ct / u.adu) |
| 107 | + elif gain.unit != u.ct / u.adu: |
| 108 | + raise u.UnitsError('gain must have units of (either ' |
| 109 | + 'astropy.units.ct, astropy.units.electron, or ' |
| 110 | + 'astropy.units.photon) / astropy.units.adu') |
| 111 | + |
| 112 | + # get the countrate from the synphot observation object |
| 113 | + countrate = synphot_obs.countrate(area=telescope.area) / telescope.gain |
| 114 | + |
| 115 | + # define counts to be in ADU, which are not technically convertible in |
| 116 | + # astropy.units: |
| 117 | + countrate = countrate.value * (u.ct / u.s) |
| 118 | + |
| 119 | + # necessary for units to work in countrate calculation: |
| 120 | + if not hasattr(snr, 'unit'): |
| 121 | + snr = snr * np.sqrt(1 * u.ct) |
| 122 | + readnoise = _get_shotnoise(telescope.readnoise) |
| 123 | + gain_err = _get_shotnoise(telescope.gain * telescope.ad_err) |
| 124 | + |
| 125 | + # solve t with the quadratic equation (pg. 57 of Howell 2000) |
| 126 | + A = countrate ** 2 |
| 127 | + B = (-1) * snr ** 2 * (countrate + npix * (background_rate + |
| 128 | + darkcurrent_rate)) |
| 129 | + C = (-1) * snr ** 2 * npix * readnoise ** 2 |
| 130 | + |
| 131 | + t = (-B + np.sqrt(B ** 2 - 4 * A * C)) / (2 * A) |
| 132 | + |
| 133 | + if gain_err.value > 1 or np.isfinite(n_background.value): |
| 134 | + from scipy.optimize import fsolve |
| 135 | + # solve t numerically |
| 136 | + t = fsolve(_t_with_small_errs, t, args=(background_rate, |
| 137 | + darkcurrent_rate, |
| 138 | + gain_err, readnoise, countrate, |
| 139 | + npix, n_background)) |
| 140 | + t = float(t) * u.s |
| 141 | + |
| 142 | + return t |
| 143 | + |
| 144 | + |
| 145 | +def _get_spectral_element(spec_tuple, SpectralElement, Empirical1D): |
| 146 | + """Returns a synphot SpectralElement made from the given 2D array""" |
| 147 | + points, lookup_table = spec_tuple |
| 148 | + return SpectralElement(Empirical1D, points=points, |
| 149 | + lookup_table=lookup_table) |
| 150 | + |
| 151 | + |
| 152 | +def _get_shotnoise(detector_property): |
| 153 | + """ |
| 154 | + Returns the shot noise (i.e. non-Poissonion noise) in the correct |
| 155 | + units. ``detector_property`` must be a Quantity. |
| 156 | + """ |
| 157 | + # Ensure detector_property is in the correct units: |
| 158 | + detector_property = detector_property.to(u.ct / u.pixel) |
| 159 | + return detector_property.value * np.sqrt(1 * (u.ct / u.pixel)) |
| 160 | + |
| 161 | + |
| 162 | +def _t_with_small_errs(t, background_rate, darkcurrent_rate, gain_err, |
| 163 | + readnoise, countrate, npix, n_background): |
| 164 | + """ |
| 165 | + Returns the full expression for the exposure time including the |
| 166 | + contribution to the noise from the background and the gain. |
| 167 | + """ |
| 168 | + if not hasattr(t, 'unit'): |
| 169 | + t = t * u.s |
| 170 | + |
| 171 | + detector_noise = (background_rate * t + darkcurrent_rate * t + |
| 172 | + gain_err ** 2 + readnoise ** 2) |
| 173 | + radicand = countrate * t + (npix * (1 + npix / n_background) * |
| 174 | + detector_noise) |
| 175 | + |
| 176 | + return countrate * t / np.sqrt(radicand) |
0 commit comments