diff --git a/.travis.yml b/.travis.yml index ed0ef553..5f7a00e5 100644 --- a/.travis.yml +++ b/.travis.yml @@ -27,7 +27,7 @@ env: - ASTROPY_VERSION=stable - MAIN_CMD='python setup.py' - CONDA_DEPENDENCIES='pytz qt pyqt six' - - PIP_DEPENDENCIES='pytest-astropy' + - PIP_DEPENDENCIES='pytest-astropy synphot' - SETUP_CMD='test -V' - CONDA_CHANNELS='astropy' diff --git a/astroplan/__init__.py b/astroplan/__init__.py index ad305f82..1abf34c7 100644 --- a/astroplan/__init__.py +++ b/astroplan/__init__.py @@ -28,5 +28,8 @@ from .constraints import * from .scheduling import * from .periodic import * + from .telescope import * + from .exptime import * + from .skycalc import * get_IERS_A_or_workaround() diff --git a/astroplan/conftest.py b/astroplan/conftest.py index 2b1c9d2f..9524e275 100644 --- a/astroplan/conftest.py +++ b/astroplan/conftest.py @@ -44,6 +44,7 @@ PYTEST_HEADER_MODULES['pyephem'] = 'ephem' PYTEST_HEADER_MODULES['matplotlib'] = 'matplotlib' PYTEST_HEADER_MODULES['pytest-mpl'] = 'pytest_mpl' + PYTEST_HEADER_MODULES['synphot'] = 'synphot' del PYTEST_HEADER_MODULES['h5py'] except KeyError: pass diff --git a/astroplan/exceptions.py b/astroplan/exceptions.py index ab4037a6..0c7b898b 100644 --- a/astroplan/exceptions.py +++ b/astroplan/exceptions.py @@ -6,7 +6,8 @@ __all__ = ["TargetAlwaysUpWarning", "TargetNeverUpWarning", "OldEarthOrientationDataWarning", "PlotWarning", - "PlotBelowHorizonWarning", "AstroplanWarning"] + "PlotBelowHorizonWarning", "AstroplanWarning", + "SkyCalcError", "UserInputError"] class AstroplanWarning(AstropyWarning): @@ -34,5 +35,23 @@ class PlotWarning(AstroplanWarning): class PlotBelowHorizonWarning(PlotWarning): - """Warning for when something is hidden on a plot because it's below the horizon""" + """ + Warning for when something is hidden on a plot because + it's below the horizon + """ + pass + + +class SkyCalcError(AstroplanWarning): + """ + Raises an error when one of the parameters isn't accepted by + the SkyCalc Sky Model Calculator. + """ + pass + + +class UserInputError(AstroplanWarning): + """ + Raises an error when the user gives an array of the wrong shape + """ pass diff --git a/astroplan/exptime.py b/astroplan/exptime.py new file mode 100644 index 00000000..0890437a --- /dev/null +++ b/astroplan/exptime.py @@ -0,0 +1,176 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +from __future__ import (absolute_import, division, print_function, + unicode_literals) + +# Third-party +import numpy as np +from astropy import units as u + + +__all__ = ['exptime_from_ccd_snr'] + + +@u.quantity_input(waveset=u.angstrom, flux=u.erg / u.s / u.cm ** 2 / u.cm, + npix=u.pixel, n_background=u.pixel, + background_rate=u.ct / u.pixel / u.s, + darkcurrent_rate=u.ct / u.pixel / u.s) +def exptime_from_ccd_snr(snr, waveset, flux, observer, telescope, + npix=1 * u.pixel, + n_background=np.inf * u.pixel, + background_rate=0 * (u.ct / u.pixel / u.s), + darkcurrent_rate=0 * (u.ct / u.pixel / u.s), + force_overlap='taper'): + """ + Returns the exposure time needed in units of seconds to achieve + the specified (idealized theoretical) signal to noise ratio + (from pg 57-58 of [1]_). + + Parameters + ---------- + snr : float, int, or `~astropy.units.Quantity` + The signal to noise ratio of the given observation in dimensionless + units. + waveset : array-like + The wavelengths associated with the target's flux. + flux : array-like + The flux of the target. + observer : `~astroplan.observer.Observer` + The Observer object. + telescope : `~astroplan.telescope.Telescope` + The Telescope object. + npix : `~astropy.units.Quantity`, optional + Number of pixels under consideration for the signal with units of + pixels. Default is 1 * astropy.units.pixel. + n_background : `~astropy.units.Quantity`, optional + Number of pixels used in the background estimation with units of + pixels. Default is set to np.inf * astropy.units.pixel + such that there is no contribution of error due to background + estimation. This assumes that n_background will be >> npix. + background_rate : `~astropy.units.Quantity`, optional + Photons per pixel per second due to the backround/sky with units of + counts/second/pixel. + Default is 0 * (astropy.units.ct / + astropy.units.second / astropy.units.pixel) + darkcurrent_rate : `~astropy.units.Quantity`, optional + Counts per pixel per second due to the dark current with units + of counts/second/pixel. + Default is 0 * (astropy.units.ct / + astropy.units.second / astropy.units.pixel) + force_overlap : {'taper', 'extrap', 'none', None}, optional + Force the spectral element x source spectrum convolution by the + specified method even when they don't fully overlap in wavelength. + 'taper' forces the incomplete spectral component to zero for its + missing wavelengths, while 'extrap' attempts to extrapolate the + incomplete part of the spectrum. + Default is 'taper' + + References + ---------- + .. [1] Howell, S. B. 2000, *Handbook of CCD Astronomy* (Cambridge, UK: + Cambridge University Press) + + Returns + ------- + t : `~astropy.units.Quantity` + The exposure time needed (in seconds) to achieve the given signal + to noise ratio. + """ + # import synphot, which isn't a required package of astroplan + from synphot.models import Empirical1D + from synphot.spectrum import SourceSpectrum, SpectralElement + from synphot.observation import Observation + + # set the source spectrum with synphot + source_spec = SourceSpectrum(Empirical1D, points=waveset, + lookup_table=flux) + + # set the spectral elements if given (quantum efficiency,skymodel,bandpass) + qe = telescope.ccd_response + skymodel = observer.skymodel + + spec_elements = _get_spectral_element(telescope.bandpass, + SpectralElement, Empirical1D) + if skymodel is not False: + spec_elements *= _get_spectral_element(skymodel, + SpectralElement, Empirical1D) + if qe is not False: + spec_elements *= _get_spectral_element(qe, + SpectralElement, Empirical1D) + + # get the synphot observation object + synphot_obs = Observation(source_spec, spec_elements, force=force_overlap) + + # make sure the gain is in the correct units + gain = telescope.gain + if gain.unit in (u.electron / u.adu, u.photon / u.adu): + gain = gain.value * (u.ct / u.adu) + elif gain.unit != u.ct / u.adu: + raise u.UnitsError('gain must have units of (either ' + 'astropy.units.ct, astropy.units.electron, or ' + 'astropy.units.photon) / astropy.units.adu') + + # get the countrate from the synphot observation object + countrate = synphot_obs.countrate(area=telescope.area) / telescope.gain + + # define counts to be in ADU, which are not technically convertible in + # astropy.units: + countrate = countrate.value * (u.ct / u.s) + + # necessary for units to work in countrate calculation: + if not hasattr(snr, 'unit'): + snr = snr * np.sqrt(1 * u.ct) + readnoise = _get_shotnoise(telescope.readnoise) + gain_err = _get_shotnoise(telescope.gain * telescope.ad_err) + + # solve t with the quadratic equation (pg. 57 of Howell 2000) + A = countrate ** 2 + B = (-1) * snr ** 2 * (countrate + npix * (background_rate + + darkcurrent_rate)) + C = (-1) * snr ** 2 * npix * readnoise ** 2 + + t = (-B + np.sqrt(B ** 2 - 4 * A * C)) / (2 * A) + + if gain_err.value > 1 or np.isfinite(n_background.value): + from scipy.optimize import fsolve + # solve t numerically + t = fsolve(_t_with_small_errs, t, args=(background_rate, + darkcurrent_rate, + gain_err, readnoise, countrate, + npix, n_background)) + t = float(t) * u.s + + return t + + +def _get_spectral_element(spec_tuple, SpectralElement, Empirical1D): + """Returns a synphot SpectralElement made from the given 2D array""" + points, lookup_table = spec_tuple + return SpectralElement(Empirical1D, points=points, + lookup_table=lookup_table) + + +def _get_shotnoise(detector_property): + """ + Returns the shot noise (i.e. non-Poissonion noise) in the correct + units. ``detector_property`` must be a Quantity. + """ + # Ensure detector_property is in the correct units: + detector_property = detector_property.to(u.ct / u.pixel) + return detector_property.value * np.sqrt(1 * (u.ct / u.pixel)) + + +def _t_with_small_errs(t, background_rate, darkcurrent_rate, gain_err, + readnoise, countrate, npix, n_background): + """ + Returns the full expression for the exposure time including the + contribution to the noise from the background and the gain. + """ + if not hasattr(t, 'unit'): + t = t * u.s + + detector_noise = (background_rate * t + darkcurrent_rate * t + + gain_err ** 2 + readnoise ** 2) + radicand = countrate * t + (npix * (1 + npix / n_background) * + detector_noise) + + return countrate * t / np.sqrt(radicand) diff --git a/astroplan/observer.py b/astroplan/observer.py index 7a269121..40d6c172 100644 --- a/astroplan/observer.py +++ b/astroplan/observer.py @@ -114,7 +114,8 @@ class Observer(object): @u.quantity_input(elevation=u.m) def __init__(self, location=None, timezone='UTC', name=None, latitude=None, longitude=None, elevation=0*u.m, pressure=None, - relative_humidity=None, temperature=None, description=None): + relative_humidity=None, temperature=None, description=None, + skymodel=False): """ Parameters ---------- @@ -152,12 +153,20 @@ def __init__(self, location=None, timezone='UTC', name=None, latitude=None, description : str (optional) A short description of the telescope, observatory or observing location. + + skymodel : 2D-array (optional). default = False + A model for atmospheric transmission. Used when calling + `~astroplan.exptime.exptime_from_ccd_snr'. + Note: the zeroth index of the given array must be the + waveset of the response function, while the other + must contain the values of the function. """ self.name = name self.pressure = pressure self.temperature = temperature self.relative_humidity = relative_humidity + self.skymodel = skymodel # If lat/long given instead of EarthLocation, convert them # to EarthLocation diff --git a/astroplan/skycalc.py b/astroplan/skycalc.py new file mode 100644 index 00000000..20de2428 --- /dev/null +++ b/astroplan/skycalc.py @@ -0,0 +1,188 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +from __future__ import (absolute_import, division, print_function, + unicode_literals) + +# Standard library +import json +import os + +# Third-party +from astropy.io import fits +import astropy.units as u + +# Local +from .exceptions import SkyCalcError + + +__all__ = ['skymodel'] + + +def skymodel(airmass=1.0, pwv_mode='pwv', season=0, + time=0, pwv=3.5, msolflux=130.0, + incl_moon='Y', moon_sun_sep=90.0, + moon_target_sep=45.0, moon_alt=45.0, + moon_earth_dist=1.0, incl_starlight='Y', + incl_zodiacal='Y', + ecl_lon=135.0, ecl_lat=90.0, + incl_loweratm='Y', incl_upperatm='Y', + incl_airglow='Y', incl_therm='N', + therm_t1=0.0, therm_e1=0.0, + therm_t2=0.0, therm_e2=0.0, therm_t3=0.0, + therm_e3=0.0, vacair='vac', wmin=300.0, + wmax=2000.0, + wgrid_mode='fixed_wavelength_step', + wdelta=0.1, wres=20000, lsf_type='none', + lsf_gauss_fwhm=5.0, lsf_boxcar_fwhm=5.0, + observatory='paranal'): + """ + Returns the model atmospheric transmittance curve queried from the SkyCalc + Sky Model Calculator. The default parameters used here are the default + parameters provided by SkyCalc: + https://www.eso.org/observing/etc/doc/skycalc/helpskycalccli.html + + Parameters + ---------- + airmass : float within range [1.0, 3.0] + Airmass. Alt and airmass are coupled through the plane parallel + approximation airmass=sec(z), z being the zenith distance + z=90°−Alt + pwv_mode : {'pwv', 'season'}. default = 'pwv' + season : {0,1,2,3,4,5,6}. default = 0 + Time of year if not in pwv mode. + (0 = all year, 1 = dec/jan, 2 = feb/mar...) + time : {0, 1, 2, 3}. default = 0 + 0 = all year, 1, 2, 3 = third of night + pwv : {-1.0, 0.5, 1.0, 1.5, 2.5, 3.5, 5.0, 7.5, 10.0, 20.0}. default = 3.5 + Precipitable Water Vapor. + msolflux : float + Monthly Averaged Solar Flux, s.f.u float > 0 (default is 130.0) + incl_moon : {'Y', 'N'}. default = 'Y' + Flag for inclusion of scattered moonlight. + Moon coordinate constraints: |z – zmoon| ≤ ρ ≤ |z + zmoon| where + ρ=moon/target separation, z=90°−target altitude and + zmoon=90°−moon altitude. + moon_sun_sep : float within range [0.0, 360.0]. default = 90.0 + Degrees of separation between Sun and Moon as seen from Earth + (i.e. the "moon phase"). + moon_target_sep : float in range [0.0, 180.0]. defualt = 45.0 + Degrees of separation between the moon and the target. + moon_alt : float in range [-90.0, 90.0] defualt = 45.0 + Moon Altitude over Horizon in degrees. + moon_earth_dist : float in range [0.91, 1.08]. default = 1.0 + Moon-Earth Distance (mean=1). + incl_starlight : {'Y', 'N'}. default = 'Y' + Flag for inclusion of scattered starlight. + incl_zodiacal : {'Y', 'N'}. default = 'Y' + Flag for inclusion of zodiacal light. + ecl_lon : float in range [-180.0, 180.0]. default = 135.0 + Heliocentric ecliptic in degrees. + ecl_lat : float in range [-90.0, 90.0]. default = 90.0 + Ecliptic latitude in degrees. + incl_loweratm : {'Y', 'N'}. default = 'Y' + Flag for inclusion of molecular emission of lower atmosphere. + incl_upperatm : {'Y', 'N'}. default = 'Y' + Flag for inclusion of molecular emission of upper atmosphere. + incl_airglow : {'Y', 'N'}. default = 'Y' + Flag for inclusion of airglow continuum (residual continuum) + incl_therm : {'Y', 'N'}. default = 'N' + Flag for inclusion of instrumental thermal radiation. + Note: This radiance component represents an instrumental effect. + The emission is provided relative to the other model components. + To obtain the correct absolute flux, an instrumental response curve + must be applied to the resulting model spectrum. + See section 6.2.4 in the SkyCalc documentation + `here `_. + therm_t1, therm_t2, therm_t3 : float. default = 0.0 + Temperature in K + therm_e1, therm_e2, therm_e3 : float in range [0.0, 1.0]. default = 0.0 + vacair : {'vac', 'air'}. default = 'vac' + In regards to the wavelength grid. + wmin : float in range [300.0, 30000.0]. default = 300.0 + Minimum wavelength (nm) in the wavelength grid. + Must be < wmax. + wmax : float in range [300.0, 30000.0]. default = 2000.0 + Maximum wavelength (nm) in the wavelength grid. + Must be > wmin. + wgrid_mode : {'fixed_spectral_resolution', 'fixed_wavelength_step', + 'user'}. + default = 'fixed_wavelength_step' + Mode of the wavelength grid. + wdelta : float in range [0, 30000.0]. default = 0.1 + Wavelength sampling step dlam (i.e. nm/step) + wres : int in range [0,1.0e6]. default = 20000 + lam/dlam where dlam is wavelength step. + wgrid_user : list of floats + default = [500.0, 510.0, 520.0, 530.0, 540.0, 550.0] + lsf_type : {'none', 'Gaussian', 'Boxcar'}. default = 'none' + Line spread function type for convolution. + lsf_gauss_fwhm : float. default = 5.0 + Gaussian full-width half-max for line spread function + wavelength bins. + lsf_boxcar_fwhm : float. default = 5.0 + Boxcar full-width half-max for line spread function + wavelength bins. + observatory : {'paranal', 'lasilla', 'armazones'}. + default = 'paranal' + Observatory where observation takes place. + + Returns + ------- + trans_waves, transmission: tuple of arrays of floats + 'trans_waves' is an array of wavelengths in angstroms (float), + 'transmission' is an array of fractional atmospheric + transmittance (float). + """ + params = locals() + + if params['observatory'] == 'lasilla': + params['observatory'] = '2400' + elif params['observatory'] == 'paranal': + params['observatory'] = '2640' + elif (params['observatory'] == '3060m' or + params['observatory'] == 'armazones'): + params['observatory'] = '3060' + else: + raise ValueError('Wrong Observatory name, please refer to the ' + 'skycalc_cli documentation.') + + # import packages not required by astroplan + import requests + + # Use the bit from skycalc_cli which queries from the SkyCalc Sky Model + server = 'http://etimecalret-001.eso.org' + url = server + '/observing/etc/api/skycalc' + + response = requests.post(url, data=json.dumps(params)) + results = json.loads(response.text) + + status = results['status'] + tmpdir = results['tmpdir'] + tmpurl = server + '/observing/etc/tmp/' + tmpdir + '/skytable.fits' + + if status == 'success': + try: + response = requests.get(tmpurl, stream=True) + data = response.content + except requests.exceptions.RequestException as e: + print(e, 'could not retrieve FITS data from server') + else: + raise SkyCalcError('HTML request failed. A custom parameter you ' + + 'set is most likely not accepted by the ' + + 'SkyCalc Calculator. json.loads() returns: ' + + '{}'.format(results)) + + # Create a temporary file to write the binary results to + tmp_data_file = './tmp_skycalc_data.fits' + + with open(tmp_data_file, 'wb') as f: + f.write(data) + + hdu = fits.open(tmp_data_file) + trans_waves = hdu[1].data["LAM"] * u.um # wavelengths + transmission = hdu[1].data["TRANS"] * u.Unit('') + + # Delete the file after reading from it + os.remove(tmp_data_file) + + return trans_waves.to(u.angstrom), transmission diff --git a/astroplan/telescope.py b/astroplan/telescope.py new file mode 100644 index 00000000..a4dfd0c8 --- /dev/null +++ b/astroplan/telescope.py @@ -0,0 +1,163 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +from __future__ import (absolute_import, division, print_function, + unicode_literals) + +# Third-party +import numpy as np +import astropy.units as u +from astroplan.exceptions import UserInputError + + +__all__ = ["Telescope"] + +AD_ERR_DEFAULT = np.sqrt(0.289) * (u.adu / u.pixel) + + +def _default_ccd_response(): + """ + A generic response function of a silicon detector from the + table on the following page: + https://www.apo.nmsu.edu/arc35m/Instruments/ARCTIC/#3p5 + """ + wl = [100., 309.45600215, 319.53901519, 329.99173871, + 339.70504127, 349.78805431, 379.53294279, 390.00376402, + 399.57293121, 409.3114413, 419.5961146, 429.14136695, + 439.52687038, 449.9795939, 459.69289647, 469.60785929, + 479.85892255, 489.5638226, 499.59835962, 509.99161922, + 519.5607864, 529.30446727, 539.85285015, 549.59976275, + 559.51472558, 569.94676599, 579.68075166, 589.59571448, + 599.9631202, 609.56008031, 619.65269622, 630.09581687, + 639.67467926, 649.50561697, 659.84070534, 669.7685951, + 679.50258077, 689.9637068, 699.78494931, 709.53555533, + 719.83463294, 729.72858948, 739.88431656, 749.9673296, + 759.66253445, 769.59042422, 780.37854306, 789.59647942, + 799.60677842, 810.07759966, 819.65646205, 829.52341053, + 839.82248813, 849.68943661, 859.77244965, 870.07152726, + 879.71760973, 889.81096433, 899.94245339, 909.73137855, + 919.69435573, 930.06545485, 939.64431724, 949.72733028, + 959.95862293, 969.6772918, 979.76030484, 990.05938245, + 1100.] + response = [0., 70.30747425, 73.11442327, 75.26541144, 76.58538173, + 74.70200949, 77.10351031, 81.79762371, 84.528972, + 87.18136276, 89.0405892, 91.16038906, 92.49142617, + 93.66048522, 94.87582372, 95.64295585, 96.56602958, + 97.20551595, 98.01709918, 98.32588679, 98.26189462, + 98.79719419, 99.10133838, 99.09379281, 99.35788748, + 99.30100555, 99.14661175, 98.81209184, 98.74843825, + 98.30855134, 98.04495971, 97.98459522, 97.24513015, + 96.85779131, 96.40002722, 96.03435769, 95.87183789, + 95.09275863, 94.89671913, 94.49854563, 94.12126754, + 93.52139537, 92.57268607, 91.77633908, 90.30410094, + 88.68846298, 86.74856762, 84.85909033, 82.85400237, + 80.76562301, 78.18504085, 75.90628116, 72.63150731, + 68.93418199, 65.34249453, 61.79608045, 58.13799205, + 54.49429826, 49.87975185, 45.18326838, 41.57397462, + 36.82027063, 32.80603172, 28.7917928, 25.09446748, + 21.39714216, 18.4392819, 14.89286782, 0.] + return np.array(wl) * u.nm, np.array(response) * u.Unit('') / 100 + + +def _user_input_test(user_input, user_input_str): + """ + Raise appropriate errors if the user doesn't input the + correct array len or type. + """ + + if len(user_input) != 2: + raise UserInputError(user_input_str + + " ccd_response must be " + "a numpy array of shape = (2, n)") + + elif (type(user_input[0]) is not u.Quantity + or type(user_input[1]) is not u.Quantity): + raise TypeError(user_input_str + + " must be a 2D-array of" + "`~astropy.units.Quantity` objects") + + +class Telescope(object): + """A class to store telescope specifications""" + @u.quantity_input(diameter=u.m, gain=u.ct / u.adu, + readnoise=u.ct / u.pixel, ad_err=u.adu / u.pixel) + def __init__(self, diameter, bandpass, gain, ccd_response=False, + readnoise=0 * (u.ct / u.pixel), ad_err=AD_ERR_DEFAULT): + """ + Parameters + ---------- + diameter : `~astropy.units.Quantity` + Diameter of the telescope aperture. + bandpass : `synphot.spectrum.SpectralElement` or 2D-array + The bandpass of the instrument/telescope. Can either be + specified by a `~synphot.spectrum.SpectralElement` object + or an array with dimensions 2 x n, where one row is the + array of wavelengths associated with the bandpass, and the other + row contains the fractional values of the bandpass. + gain : `~astropy.units.Quantity`, optional + Gain of the CCD with units of counts/ADU. Default is + 1 * astropy.units.ct / astropy.units.adu such that the + contribution to the error due to the gain is assumed to be small. + ccd_response : 'default' or 2D-array of `~astropy.units.Quantity` + objects + Note: the zeroth index of the given array must be the + waveset of the response function, while the other + must contain the values of the function. + readnoise : `~astropy.units.Quantity`, optional + Counts per pixel from the read noise with units of counts/pixel. + Default is 0 * astropy.units.ct / astropy.units.pixel. + ad_err : `~astropy.units.Quantity`, optional + An estimate of the 1 sigma error within the A/D converter with + units of adu/pixel. Default is set to + sqrt(0.289) * astropy.units.adu / astropy.units.pixel, where + sqrt(0.289) comes from the assumption that "for a charge + level that is half way between two output ADU steps, there + is an equal chance that it will be assigned to the lower + or to the higher ADU value when converted to a digital number" + (text from footnote on page 56 of [1]_, see also [2]_). + + References + ---------- + .. [1] Howell, S. B. 2000, *Handbook of CCD Astronomy* (Cambridge, + UK: Cambridge University Press) + .. [2] Merline, W. & Howell, S. B. *A Realistic Model for + Point-sources Imaged on Array Detectors: The Model and + Initial Results*. ExA, 6:163 (1995) + """ + self.diameter = diameter + self.gain = gain + self.ad_err = ad_err + self.readnoise = readnoise + self._bandpass = bandpass + self._ccd_response = ccd_response + + @property + def area(self): + """Area of the telescope""" + return np.pi * (self.diameter / 2) ** 2 + + @property + def bandpass(self): + """Bandpass of the instrument""" + from synphot.spectrum import SpectralElement + + if isinstance(self._bandpass, SpectralElement): + waves = self._bandpass.waveset + return waves, self._bandpass(waves) + + else: + _user_input_test(self._bandpass, "bandpass") + + return self._bandpass + + @property + def ccd_response(self): + """Response function of the CCD, i.e. the quantum efficiency""" + if self._ccd_response == 'default': + return _default_ccd_response() + + elif self._ccd_response is False: + return self._ccd_response + + else: + _user_input_test(self._ccd_response, "ccd_response") + + return self._ccd_response[0], self._ccd_response[1] diff --git a/astroplan/tests/test_exptime.py b/astroplan/tests/test_exptime.py new file mode 100644 index 00000000..eded1296 --- /dev/null +++ b/astroplan/tests/test_exptime.py @@ -0,0 +1,190 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +from __future__ import (absolute_import, division, print_function, + unicode_literals) + +# Third-party +import astropy.units as u +from astropy.tests.helper import assert_quantity_allclose +import numpy as np +from astroplan.exceptions import UserInputError +from astropy.tests.helper import raises +import astropy.units as u +import numpy as np +from astroplan import FixedTarget, Observer, Telescope +from astroquery.gaia import Gaia +from astropy.utils.data import download_file + +# Local +from astroplan import exptime_from_ccd_snr + +@raises(TypeError) +def _raise_type_error_test(testcase): + return testcase + +class TestExptime(object): + """Tests for exptime.py""" + @pytest.mark.skipif(reason="synphot not installed") + def setup_class(self): + + from synphot.spectrum import SpectralElement + from synphot.models import ConstFlux1D + + def test_synphot_specmodel(self): + """ + Tests that exptime returns the expected result given a synphot + spectral model, for both numeric and analytic solutions. + This test is built off of the exptime tutorial in the + documentation. + """ + # numeric, i.e. when n_background != inf or gain_err > 1 + # note that gain_err > 1 by default, since gain_err = gain * ad_err, where + # ad_err is >1 by default + snr = 1 + + waveset = np.array([1, 2]) * u.angstrom + flux = SpectralElement(ConstFlux1D, 1 * (u.photon / u.s / u.cm ** 2 / u.angstrom)) + + area = 1 * u.cm ** 2 + diameter = np.sqrt(area / np.pi) * 2 + gain = 1 * u.photon / u.adu + throughput = np.ones(len(waveset)) * u.Unit('') + bandpass = SpectralElement(Empirical1D, points=waveset, lookup_table=throughput) + ad_err = 0 * u.adu / u.pixel + telescope = Telescope(diameter, bandpass, gain, ad_err=ad_err) + + result = exptime_from_ccd_snr(snr, waveset, flux, observer, telescope) + assert_quantity_allclose(result, 1 * u.s, rtol=1e-3) + + # numeric, i.e. when n_background != inf or gain_err > 1 + telescope = Telescope(diameter, bandpass, gain) + # gain error doesn't have that big of an effect so the result shouldn't + # be too different + assert_quantity_allclose(result, 1 * u.s, rtol=1e-1) + + def test_user_specmodel(self): + """ + Tests that exptime returns the expected result given a user + inputed spectral model, for both numeric and analytic solutions + """ + # numeric, i.e. when n_background != inf or gain_err > 1 + # note that gain_err > 1 by default, since gain_err = gain * ad_err, where + # ad_err is >1 by default + T_eff = 4800.0 + distance_scale = 6.701456436998426e-20 + + # setting up the synphot spectrum + flux_url = ('ftp://phoenix.astro.physik.uni-goettingen.de/' + 'v2.0/HiResFITS/PHOENIX-ACES-AGSS-COND-2011/Z-0.0/' + 'lte{T_eff:05d}-{log_g:1.2f}-0.0.PHOENIX-ACES-AGSS-' + 'COND-2011-HiRes.fits').format(T_eff=int(T_eff), + log_g=4.5) + wavelength_url = ('ftp://phoenix.astro.physik.uni-goettingen.de/' + 'v2.0/HiResFITS/WAVE_PHOENIX-ACES-AGSS-COND-' + '2011.fits') + flux_of_target = fits.getdata(flux_url) * distance_scale + waveset_of_target = fits.getdata(wavelength_url) + + # setting up the telescope + diameter = 3.5 * u.m + svo_link = ('http://svo2.cab.inta-csic.es/theory/fps3/fps.php?ID=') + filt_path = download_file(svo_link + 'SLOAN/SDSS.rprime_filter') + bandpass = SpectralElement.from_file(filt_path) + gain = 1.9 * (u.ct / u.adu) + + snr = 100 + flux = flux_of_target * (u.erg / u.s / u.cm ** 2 / u.cm) + waveset = waveset_of_target * u.Angstrom + observer = Observer.at_site('apo') + telescope = Telescope(diameter, bandpass, gain, ccd_response='default') + result = exptime_from_ccd_snr(snr, waveset, flux, observer, telescope) + assert_quantity_allclose(result, 0.0022556432 * u.s, rtol=1e-3) + + # analytic + waveset = np.array([1, 2]) * u.angstrom + throughput = np.ones(len(waveset)) * u.Unit('') + bandpass = SpectralElement(Empirical1D, points=waveset, lookup_table=throughput) + flux = np.ones(len(waveset)) * (u.photon / u.s / u.cm ** 2 / u.angstrom) + ad_err = 0 * u.adu / u.pixel + + @raises(TypeError) + def test_user_specmodel_typeerror(self): + """ + Tests that exptime raises an error when the user gives a spectral + model of the wrong type (ie not an array of + `~astropy.units.quantity.Quantity` objects) + """ + + @raises(UserInputError) + def test_user_specmodel_lenerror(self): + """ + Tests that exptime raises an error when the user gives a spectral + model of the wrong length + """ + + def test_optional_param_type_errors_raised(self): + """ + Tests that exptime raises type errors whenever any optional + parameters is given in the incorrect unit + """ + _raise_type_error_test(testcase1) + _raise_type_error_test(testcase2) + _raise_type_error_test(testcase3) + + +def test_t_exp_numeric(): + """ + A test to check that the numerical method in + observation.exptime_from_ccd_snr (i.e. when the error in background + noise or gain is non-negligible) is done correctly. Based on the worked + example on pg. 56-57 of [1]_. + + References + ---------- + .. [1] Howell, S. B. 2000, *Handbook of CCD Astronomy* (Cambridge, UK: + Cambridge University Press) + """ + t = 300 * u.s + snr = 342 + gain = 5 * (u.ct / u.adu) + countrate = 24013 * u.adu * gain / t + npix = 1 * u.pixel + n_background = 200 * u.pixel + background_rate = 620 * (u.adu / u.pixel) * gain / t + darkcurrent_rate = 22 * (u.ct / u.pixel / u.hr) + readnoise = 5 * (u.ct / u.pixel) + + result = exptime_from_ccd_snr(snr, countrate, npix=npix, + n_background=n_background, + background_rate=background_rate, + darkcurrent_rate=darkcurrent_rate, + readnoise=readnoise, gain=gain) + + assert_quantity_allclose(result, t, atol=1 * u.s) + + +def test_t_exp_analytic(): + """ + A test to check that the analytic method in + observation.exptime_from_ccd_snr is done correctly. + """ + snr_set = 50 + countrate = 1000 * (u.ct / u.s) + npix = 1 * u.pixel + background_rate = 100 * (u.ct / u.pixel / u.s) + darkcurrent_rate = 5 * (u.ct / u.pixel / u.s) + readnoise = 1 * (u.ct / u.pixel) + + t = exptime_from_ccd_snr(snr_set, countrate, npix=npix, + background_rate=background_rate, + darkcurrent_rate=darkcurrent_rate, + readnoise=readnoise) + + # if t is correct, ccd_snr() should return snr_set: + snr_calc = ccd_snr(countrate * t, + npix=npix, + background=background_rate * t, + darkcurrent=darkcurrent_rate * t, + readnoise=readnoise) + + assert_quantity_allclose(snr_calc, snr_set, + atol=0.5 * u.Unit("")) \ No newline at end of file diff --git a/astroplan/tests/test_skycalc.py b/astroplan/tests/test_skycalc.py new file mode 100644 index 00000000..900e19e0 --- /dev/null +++ b/astroplan/tests/test_skycalc.py @@ -0,0 +1,25 @@ +from astropy.tests.helper import raises +from astroplan import skymodel +from astroplan.exceptions import SkyCalcError + + +@raises(ValueError) +def test_null_observatory(): + """Tests that skymodel raises an error if an observatory + which SkyCalc doesn't support is given.""" + return skymodel(observatory='apo') + + +@raises(SkyCalcError) +def test_null_parameter(): + """Tests that skymodel raises an error if given + a parameter which SkyCalc doesn't support.""" + return skymodel(airmass='test') + + +def test_skymodel_units(): + """Tests that skymodel returns the atmospheric + transmission function as a `~astropy.units.quantity.Quantity` + """ + sm = skymodel() + assert hasattr(sm[0], "unit") and hasattr(sm[1], "unit") diff --git a/astroplan/tests/test_telescope.py b/astroplan/tests/test_telescope.py new file mode 100644 index 00000000..2c7b90bc --- /dev/null +++ b/astroplan/tests/test_telescope.py @@ -0,0 +1,197 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +from __future__ import (absolute_import, division, print_function, + unicode_literals) +import astropy.units as u +from ..telescope import Telescope +import pytest +import numpy as np +from astroplan.exceptions import UserInputError +from astropy.tests.helper import raises + + +class TestTelescope(object): + """Tests for telescope.py""" + @pytest.mark.skipif(reason="synphot not installed") + def setup_class(self): + + from synphot.spectrum import SpectralElement + from astropy.utils.data import download_file + + svo_link = ('http://svo2.cab.inta-csic.es/theory/fps3/fps.php?ID=') + filt_path = download_file(svo_link + 'SLOAN/SDSS.rprime_filter') + + self.diameter = 3.5 * u.m + self.bandpass = SpectralElement.from_file(filt_path) + self.gain = 1.9 * (u.ct / u.adu) + self.default_readnoise = 0 * (u.ct / u.pixel) + self.default_ad_err = np.sqrt(0.289) * (u.adu / u.pixel) + self.telescope = Telescope(self.diameter, self.bandpass, self.gain) + + def test_diameter(self): + """ + Tests that Telescope returns the correct diameter + """ + assert(self.telescope.diameter == self.diameter) + + def test_bandpass(self): + """ + Tests that Telescope returns the correct bandpass + """ + waves = self.bandpass.waveset + assert(all(self.telescope.bandpass[0] == waves) and + all(self.telescope.bandpass[1] == self.bandpass(waves))) + + def test_gain(self): + """ + Tests that Telescope returns the correct gain + """ + assert(self.telescope.gain == self.gain) + + def test_ccd_response(self): + """ + Tests that Telescope.ccd_response is False if a ccd response + function isn't given. + """ + assert(self.telescope.ccd_response is False) + + def test_readnoise(self): + """Tests that Telescope returns the default readnoise""" + assert(self.telescope.readnoise == self.default_readnoise) + + def test_ad_err(self): + """Tests that Telescope returns the default ad_err""" + assert(self.telescope.ad_err == self.default_ad_err) + + @raises(TypeError) + def test_diameter_type(self): + """ + Tests that Telescope raises an error if the user doesn't + gives the bandpass as a unit quantity + """ + return Telescope(1, self.bandpass, self.gain).diameter + + @raises(TypeError) + def test_bandpass_type(self): + """ + Tests that Telescope raises an error if the user doesn't + gives the bandpass as a unit quantity + """ + return Telescope(self.diameter, np.array([[], []]), self.gain).bandpass + + @raises(UserInputError) + def test_bandpass_len(self): + """ + Tests that Telescope raises an error if the user doesn't + give a bandpass of the correct shape + """ + return Telescope(self.diameter, np.array([]), self.gain).bandpass + + @raises(TypeError) + def test_gain_type(self): + """ + Tests that Telescope raises an error if gain isn't given + as a `~astropy.units.quantity.Quantity`. + """ + return Telescope(self.diameter, self.bandpass, 1).gain + + def test_default_ccd_response(self): + """ + Tests that Telescope.ccd_response returns the following array + if set to 'default'. + """ + wl = np.array([100., 309.45600215, 319.53901519, 329.99173871, + 339.70504127, 349.78805431, 379.53294279, 390.00376402, + 399.57293121, 409.3114413, 419.5961146, 429.14136695, + 439.52687038, 449.9795939, 459.69289647, 469.60785929, + 479.85892255, 489.5638226, 499.59835962, 509.99161922, + 519.5607864, 529.30446727, 539.85285015, 549.59976275, + 559.51472558, 569.94676599, 579.68075166, 589.59571448, + 599.9631202, 609.56008031, 619.65269622, 630.09581687, + 639.67467926, 649.50561697, 659.84070534, 669.7685951, + 679.50258077, 689.9637068, 699.78494931, 709.53555533, + 719.83463294, 729.72858948, 739.88431656, 749.9673296, + 759.66253445, 769.59042422, 780.37854306, 789.59647942, + 799.60677842, 810.07759966, 819.65646205, 829.52341053, + 839.82248813, 849.68943661, 859.77244965, 870.07152726, + 879.71760973, 889.81096433, 899.94245339, 909.73137855, + 919.69435573, 930.06545485, 939.64431724, 949.72733028, + 959.95862293, 969.6772918, 979.76030484, 990.05938245, + 1100.]) * u.nm + response = np.array([0., 70.30747425, 73.11442327, 75.26541144, + 76.58538173, 74.70200949, 77.10351031, + 81.79762371, 84.528972, 87.18136276, + 89.0405892, 91.16038906, 92.49142617, + 93.66048522, 94.87582372, 95.64295585, + 96.56602958, 97.20551595, 98.01709918, + 98.32588679, 98.26189462, 98.79719419, + 99.10133838, 99.09379281, 99.35788748, + 99.30100555, 99.14661175, 98.81209184, + 98.74843825, 98.30855134, 98.04495971, + 97.98459522, 97.24513015, 96.85779131, + 96.40002722, 96.03435769, 95.87183789, + 95.09275863, 94.89671913, 94.49854563, + 94.12126754, 93.52139537, 92.57268607, + 91.77633908, 90.30410094, 88.68846298, + 86.74856762, 84.85909033, 82.85400237, + 80.76562301, 78.18504085, 75.90628116, + 72.63150731, 68.93418199, 65.34249453, + 61.79608045, 58.13799205, 54.49429826, + 49.87975185, 45.18326838, 41.57397462, + 36.82027063, 32.80603172, 28.7917928, + 25.09446748, 21.39714216, 18.4392819, + 14.89286782, 0.]) * u.Unit('') / 100 + tele = Telescope(self.diameter, self.bandpass, self.gain, + ccd_response='default') + assert(all(wl == tele.ccd_response[0]) and + all(response == tele.ccd_response[1])) + + def test_custom_ccd_response(self): + """ + Tests that Telescope.ccd_response gives back the user specified + ccd response function. + """ + wl = np.arange(10) * u.angstrom + response = np.arange(10) * u.Unit('') + custom_ccd_response = (wl, response) + + tele = Telescope(self.diameter, self.bandpass, self.gain, + ccd_response=custom_ccd_response) + assert(tele.ccd_response == custom_ccd_response) + + @raises(TypeError) + def test_custom_ccd_response_type(self): + """ + Tests that telescope raises an error if ccd_response isn't passed in + as a `~astropy.units.quantity.Quantity`. + """ + wl = np.arange(10) + response = np.arange(10) + custom_ccd_response = (wl, response) + return Telescope(self.diameter, self.bandpass, self.gain, + ccd_response=custom_ccd_response).ccd_response + + @raises(UserInputError) + def test_custom_ccd_response_len(self): + """ + Tests that telescope raises an error if the given ccd_response isn't + the right length + """ + return Telescope(self.diameter, self.bandpass, self.gain, + ccd_response=np.array([])).ccd_response + + @raises(TypeError) + def test_readnoise_type(self): + """Tests that Telescope raises an error if readnoise isn't passed in + as a `~astropy.units.quantity.Quantity`. + """ + return Telescope(self.diameter, self.bandpass, self.gain, + readnoise=0) + + @raises(TypeError) + def test_ad_err_type(self): + """ + Tests that telescope raises an error if ad_err isn't given + as a `~astropy.units.quantity.Quantity`. + """ + return Telescope(self.diameter, self.bandpass, self.gain, + ad_err=1).ad_err diff --git a/docs/index.rst b/docs/index.rst index c0e4fea9..abc3542b 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -24,8 +24,9 @@ Features: observatories anywhere on Earth * Built-in plotting convenience functions for standard observation planning plots (airmass, parallactic angle, sky maps). -* Determining observability of sets of targets given an arbitrary set of +* Determine observability of sets of targets given an arbitrary set of constraints (i.e., altitude, airmass, moon separation/illumination, etc.). +* Calculate the exposure time needed to obtain a given signal to noise ratio * `Astropy`_ powered! Links @@ -71,3 +72,4 @@ Contributors * Erik Tollerud * Brigitta Sipocz * Karl Vyhmeister +* Tiffany Jansen diff --git a/docs/tutorials/exptime.rst b/docs/tutorials/exptime.rst new file mode 100644 index 00000000..d6d42c1c --- /dev/null +++ b/docs/tutorials/exptime.rst @@ -0,0 +1,498 @@ +.. _exptime_tutorial: + +.. doctest-skip-all + +**************************** +Predicting an Exposure Time +**************************** + +This module of astroplan requires that `scipy `_ +and `synphot `_ be installed. + +Contents +======== + + In this tutorial we will show how to predict the exposure time needed + to obtain a given signal to noise ratio for a target star. + +* :ref:`exptime-querying_target_properties` +* :ref:`exptime-obtaining_a_model_spectrum` +* :ref:`exptime-setting_the_location_of_the_observer` +* :ref:`exptime-creating_a_telescope_model` +* :ref:`exptime-getting_the_exposure_time` +* :ref:`exptime-adding_an_atmospheric_transmission_model` + +.. code-block:: python + + >>> import astropy.units as u + + >>> import numpy as np + >>> import matplotlib.pyplot as plt + +.. _exptime-querying_target_properties: + +Querying Target Properties +========================== + +In this example we will use HAT-P-11 as the target. We will model the +spectrum of HAT-P-11 with a PHOENIX stellar atmosphere model, which we +can obtain by specifying a stellar temperature and querying from `the +Gottingen Spectral Library. `_ + +To get the stellar temperature, we can use :meth:`astroplan.FixedTarget.from_name` +to get the coordinates of HAT-P-11, then give these coordinates to +`astroquery.gaia`: + +.. code-block:: python + + >>> from astroplan import FixedTarget + >>> from astroquery.gaia import Gaia + + >>> hatp11 = FixedTarget.from_name('HAT-P-11') + + >>> # width / height of search: + >>> width = u.Quantity(1, u.arcmin) + >>> height = u.Quantity(1, u.arcmin) + >>> search_results = Gaia.query_object_async(coordinate=hatp11.coord, + width=width, height=height) + + >>> # the queried star should be the one nearest to the given coordinates + >>> search_results.add_index('dist', unique=True) + >>> hatp11_info = search_results.loc['dist', min(search_results['dist'])] + + >>> T_eff = round(hatp11_info['teff_val'], -2) # round to nearest 100 K + +Since the spectrum is defined at the source, we have to scale for the +distance to the target to get the flux at the telescope: + +.. code-block:: python + + >>> stellar_radius = hatp11_info['radius_val'] * u.R_sun + >>> parallax = hatp11_info['parallax'] * u.mas + >>> # parallax is given by Gaia in milliarcseconds, so convert to parsec: + >>> distance = parallax.to(u.parsec, equivalencies=u.parallax()) + + >>> distance_scale = float(stellar_radius / distance) ** 2 / np.pi + +.. _exptime-obtaining_a_model_spectrum: + +Obtaining a Model Spectrum +========================== + +If you don't already have a model spectrum, you may query one (at least for +a star) from the Gottingen Spectral Library using `astropy.io.fits.getdata`: + +.. code-block:: python + + >>> from astropy.io import fits + + >>> flux_url = ('ftp://phoenix.astro.physik.uni-goettingen.de/' + 'v2.0/HiResFITS/PHOENIX-ACES-AGSS-COND-2011/Z-0.0/' + 'lte{T_eff:05d}-{log_g:1.2f}-0.0.PHOENIX-'ACES-AGSS-' + 'COND-2011-HiRes.fits').format(T_eff=int(T_eff), + log_g=4.5) + >>> wavelength_url = ('ftp://phoenix.astro.physik.uni-goettingen.de/' + 'v2.0/HiResFITS/WAVE_PHOENIX-ACES-AGSS-COND-' + '2011.fits') + + >>> flux_of_target = fits.getdata(flux_url) * distance_scale + >>> waveset_of_target = fits.getdata(wavelength_url) + +The units we use below are specified by the PHOENIX models. + +.. note:: + It is good practice (and for some packages, essential) to attach units to + quantities so that you can be sure the final results are scaled correctly. + +.. code-block:: python + + >>> flux_of_target = flux_of_target * (u.erg / u.s / u.cm ** 2 / u.cm) + >>> waveset_of_target = waveset_of_target * u.Angstrom + + >>> plt.plot(waveset_of_target, flux_of_target) + +.. plot:: + + import astropy.units as u + + import numpy as np + import matplotlib.pyplot as plt + + from astroplan import FixedTarget + from astroquery.gaia import Gaia + + hatp11 = FixedTarget.from_name('Hat-p-11') + + # width / height of search: + width = u.Quantity(1, u.arcmin) + height = u.Quantity(1, u.arcmin) + search_results = Gaia.query_object_async(coordinate=hatp11.coord, width=width, height=height) + + # the queried star should be the one nearest to the given coordinates + search_results.add_index('dist', unique=True) + hatp11_info = search_results.loc['dist', min(search_results['dist'])] + + T_eff = round(hatp11_info['teff_val'], -2) # round to nearest 100 K + + stellar_radius = hatp11_info['radius_val'] * u.R_sun + parallax = hatp11_info['parallax'] * u.mas + # parallax is given by Gaia in milliarcseconds, so convert to parsec: + distance = parallax.to(u.parsec, equivalencies=u.parallax()) + + distance_scale = float(stellar_radius / distance) ** 2 / np.pi + + from astropy.io import fits + + flux_url = ('ftp://phoenix.astro.physik.uni-goettingen.de/v2.0/HiResFITS/' + 'PHOENIX-ACES-AGSS-COND-2011/Z-0.0/lte{T_eff:05d}-{log_g:1.2f}-0.0.PHOENIX-' + 'ACES-AGSS-COND-2011-HiRes.fits').format(T_eff=int(T_eff), log_g=4.5) + wavelength_url = ('ftp://phoenix.astro.physik.uni-goettingen.de/v2.0/HiResFITS/' + 'WAVE_PHOENIX-ACES-AGSS-COND-2011.fits') + + flux_of_target = fits.getdata(flux_url) * distance_scale + waveset_of_target = fits.getdata(wavelength_url) + + flux_of_target = flux_of_target * (u.erg / u.s / u.cm ** 2 / u.cm) + waveset_of_target = waveset_of_target * u.Angstrom + + plt.plot(waveset_of_target, flux_of_target) + plt.show() + +.. _exptime-setting_the_location_of_the_observer: + +Setting the Location of the Observer +==================================== + +Create an :class:`astroplan.observer.Observer` object via the +:meth:`~astroplan.observer.Observer.at_site` method to set the +location of the observing site: + +.. code-block:: python + + >>> from astroplan import Observer + + >>> observer = Observer.at_site('apo') + +.. _exptime-creating_a_telescope_model: + +Creating a Telescope Model +========================== + +Create a :class:`astroplan.telescope.Telescope` object, which models the telescope +you wish to use (we model ours off of `the ARCTIC instrument `_ on APO's 3.5m telescope), and input: + +- the diameter of the telescope aperture +- a model of your bandpass (we use `synphot `_ + to generate a bandpass downloaded from `the Spanish Virtual Observatory `_) +- the gain of the instrument, which is usually found on the observatory's website + +.. code-block:: python + + >>> from synphot.spectrum import SpectralElement + >>> from astropy.utils.data import download_file + + >>> diameter = 3.5 * u.m + >>> gain = 1.9 * (u.ct / u.adu) + + >>> # model the bandpass: + >>> svo_link = ('http://svo2.cab.inta-csic.es/' + + 'theory/fps3/fps.php?ID=') + >>> filt_path = download_file(svo_link + 'SLOAN/SDSS.rprime_filter') + >>> bandpass = SpectralElement.from_file(filt_path) + >>> bandpass.plot() + +.. plot:: + + import astropy.units as u + + import numpy as np + import matplotlib.pyplot as plt + + from astroplan import FixedTarget + from astroquery.gaia import Gaia + + hatp11 = FixedTarget.from_name('Hat-p-11') + + # width / height of search: + width = u.Quantity(1, u.arcmin) + height = u.Quantity(1, u.arcmin) + search_results = Gaia.query_object_async(coordinate=hatp11.coord, width=width, height=height) + + # the queried star should be the one nearest to the given coordinates + search_results.add_index('dist', unique=True) + hatp11_info = search_results.loc['dist', min(search_results['dist'])] + + T_eff = round(hatp11_info['teff_val'], -2) # round to nearest 100 K + + stellar_radius = hatp11_info['radius_val'] * u.R_sun + parallax = hatp11_info['parallax'] * u.mas + # parallax is given by Gaia in milliarcseconds, so convert to parsec: + distance = parallax.to(u.parsec, equivalencies=u.parallax()) + + distance_scale = float(stellar_radius / distance) ** 2 / np.pi + + from astropy.io import fits + + flux_url = ('ftp://phoenix.astro.physik.uni-goettingen.de/v2.0/HiResFITS/' + 'PHOENIX-ACES-AGSS-COND-2011/Z-0.0/lte{T_eff:05d}-{log_g:1.2f}-0.0.PHOENIX-' + 'ACES-AGSS-COND-2011-HiRes.fits').format(T_eff=int(T_eff), log_g=4.5) + wavelength_url = ('ftp://phoenix.astro.physik.uni-goettingen.de/v2.0/HiResFITS/' + 'WAVE_PHOENIX-ACES-AGSS-COND-2011.fits') + + flux_of_target = fits.getdata(flux_url) * distance_scale + waveset_of_target = fits.getdata(wavelength_url) + + flux_of_target = flux_of_target * (u.erg / u.s / u.cm ** 2 / u.cm) + waveset_of_target = waveset_of_target * u.Angstrom + + from astroplan import Observer + + observer = Observer.at_site('apo') + + from synphot.spectrum import SpectralElement + from astropy.utils.data import download_file + from astroplan import Telescope + + diameter = 3.5 * u.m + gain = 1.9 * (u.ct / u.adu) + + # model the bandpass: + svo_link = ('http://svo2.cab.inta-csic.es/' + + 'theory/fps3/fps.php?ID=') + filt_path = download_file(svo_link + 'SLOAN/SDSS.rprime_filter') + bandpass = SpectralElement.from_file(filt_path) + + bandpass.plot() + +You can also give it a CCD response function, either from your own model or +the default model given by :class:`~astroplan.telescope.Telescope` as we've done here: + +.. code-block:: python + + >>> from astroplan import Telescope + + >>> telescope = Telescope(diameter, bandpass, gain, ccd_response='default') + +.. code-block:: python + + >>> ccd_response_wls, ccd_response = telescope.ccd_response + >>> plt.plot(ccd_response_wls, ccd_response) + >>> plt.xlabel(ccd_response_wls.unit) + +.. plot:: + + import astropy.units as u + + import numpy as np + import matplotlib.pyplot as plt + + from astroplan import FixedTarget + from astroquery.gaia import Gaia + + hatp11 = FixedTarget.from_name('HAT-P-11') + + # width / height of search: + width = u.Quantity(1, u.arcmin) + height = u.Quantity(1, u.arcmin) + search_results = Gaia.query_object_async(coordinate=hatp11.coord, width=width, height=height) + + # the queried star should be the one nearest to the given coordinates + search_results.add_index('dist', unique=True) + hatp11_info = search_results.loc['dist', min(search_results['dist'])] + + T_eff = round(hatp11_info['teff_val'], -2) # round to nearest 100 K + + stellar_radius = hatp11_info['radius_val'] * u.R_sun + parallax = hatp11_info['parallax'] * u.mas + # parallax is given by Gaia in milliarcseconds, so convert to parsec: + distance = parallax.to(u.parsec, equivalencies=u.parallax()) + + distance_scale = float(stellar_radius / distance) ** 2 / np.pi + + from astropy.io import fits + + flux_url = ('ftp://phoenix.astro.physik.uni-goettingen.de/v2.0/HiResFITS/' + 'PHOENIX-ACES-AGSS-COND-2011/Z-0.0/lte{T_eff:05d}-{log_g:1.2f}-0.0.PHOENIX-' + 'ACES-AGSS-COND-2011-HiRes.fits').format(T_eff=int(T_eff), log_g=4.5) + wavelength_url = ('ftp://phoenix.astro.physik.uni-goettingen.de/v2.0/HiResFITS/' + 'WAVE_PHOENIX-ACES-AGSS-COND-2011.fits') + + flux_of_target = fits.getdata(flux_url) * distance_scale + waveset_of_target = fits.getdata(wavelength_url) + + flux_of_target = flux_of_target * (u.erg / u.s / u.cm ** 2 / u.cm) + waveset_of_target = waveset_of_target * u.Angstrom + + from astroplan import Observer + + observer = Observer.at_site('apo') + + from synphot.spectrum import SpectralElement + from astropy.utils.data import download_file + from astroplan import Telescope + + diameter = 3.5 * u.m + gain = 1.9 * (u.ct / u.adu) + + # model the bandpass: + svo_link = ('http://svo2.cab.inta-csic.es/' + + 'theory/fps3/fps.php?ID=') + filt_path = download_file(svo_link + 'SLOAN/SDSS.rprime_filter') + bandpass = SpectralElement.from_file(filt_path) + + telescope = Telescope(diameter, bandpass, gain, ccd_response='default') + ccd_response_wls, ccd_response = telescope.ccd_response + plt.plot(ccd_response_wls, ccd_response) + plt.xlabel(ccd_response_wls.unit) + plt.show() + +.. _exptime-getting_the_exposure_time: + +Getting the Exposure Time +========================= + +Then specify the signal to noise ratio you wish to reach and +call the :func:`~astroplan.exptime_from_ccd_snr()` function: + +.. code-block:: python + + >>> from astroplan import exptime_from_ccd_snr + + >>> snr = 100 + >>> exptime_from_ccd_snr(snr, waveset_of_target, flux_of_target, observer, telescope) + 0.0022556432s + +.. note:: + SNR is derived from the equations on pages 57-58 of Howell (2000). See the + docstring of :func:`~astroplan.exptime_from_ccd_snr()` for more details on + the reference and optional input parameters. + +.. _exptime-adding_an_atmospheric_transmission_model: + +(optional) Adding an Atmospheric Transmission Model +=================================================== + +Say you want to add the effect of the atmosphere to make the exposure +time prediction more realistic. You can pass your own sky model as an +array into your observer object, or you can query from the `SKYCALC +Sky Model Calculator `_ with :func:`astroplan.skymodel()`: + +.. code-block:: python + + >>> from astroplan import skymodel + + >>> skymodel_1pt5airmass = skymodel(airmass=1.5) + +`skymodel` has many different sky parameters that can be set, including +the precipitable water vapor, whether or not to include the moon, and more. +For options and defaults, see :func:`astroplan.skymodel()`. + +The skymodel attribute of the `Observer` object must then be set: + +.. code-block:: python + + >>> observer.skymodel = skymodel_1pt5airmass + + >>> skymodel_waveset, skymodel_flux = observer.skymodel + >>> plt.plot(skymodel_waveset, skymodel_flux) + +.. note:: + + Alternatively, you could pass it in as an optional argument + when initiating the observer object: + + .. code-block:: python + + >>> observer = Observer.at_site('apo', skymodel=skymodel_1pt5airmass) + +.. plot:: + + import astropy.units as u + + import numpy as np + import matplotlib.pyplot as plt + + from astroplan import FixedTarget + from astroquery.gaia import Gaia + + hatp11 = FixedTarget.from_name('Hat-p-11') + + # width / height of search: + width = u.Quantity(1, u.arcmin) + height = u.Quantity(1, u.arcmin) + search_results = Gaia.query_object_async(coordinate=hatp11.coord, width=width, height=height) + + # the queried star should be the one nearest to the given coordinates + search_results.add_index('dist', unique=True) + hatp11_info = search_results.loc['dist', min(search_results['dist'])] + + T_eff = round(hatp11_info['teff_val'], -2) # round to nearest 100 K + + stellar_radius = hatp11_info['radius_val'] * u.R_sun + parallax = hatp11_info['parallax'] * u.mas + # parallax is given by Gaia in milliarcseconds, so convert to parsec: + distance = parallax.to(u.parsec, equivalencies=u.parallax()) + + distance_scale = float(stellar_radius / distance) ** 2 / np.pi + + from astropy.io import fits + + flux_url = ('ftp://phoenix.astro.physik.uni-goettingen.de/v2.0/HiResFITS/' + 'PHOENIX-ACES-AGSS-COND-2011/Z-0.0/lte{T_eff:05d}-{log_g:1.2f}-0.0.PHOENIX-' + 'ACES-AGSS-COND-2011-HiRes.fits').format(T_eff=int(T_eff), log_g=4.5) + wavelength_url = ('ftp://phoenix.astro.physik.uni-goettingen.de/v2.0/HiResFITS/' + 'WAVE_PHOENIX-ACES-AGSS-COND-2011.fits') + + flux_of_target = fits.getdata(flux_url) * distance_scale + waveset_of_target = fits.getdata(wavelength_url) + + flux_of_target = flux_of_target * (u.erg / u.s / u.cm ** 2 / u.cm) + waveset_of_target = waveset_of_target * u.Angstrom + + from astroplan import Observer + + observer = Observer.at_site('apo') + + from synphot.spectrum import SpectralElement + from astropy.utils.data import download_file + from astroplan import Telescope + + diameter = 3.5 * u.m + gain = 1.9 * (u.ct / u.adu) + + # model the bandpass: + svo_link = ('http://svo2.cab.inta-csic.es/' + + 'theory/fps3/fps.php?ID=') + filt_path = download_file(svo_link + 'SLOAN/SDSS.rprime_filter') + bandpass = SpectralElement.from_file(filt_path) + + telescope = Telescope(diameter, bandpass, gain, ccd_response='default') + ccd_response_wls, ccd_response = telescope.ccd_response + + from astroplan import exptime_from_ccd_snr + + snr = 100 + exptime_from_ccd_snr(snr, waveset_of_target, flux_of_target, observer, telescope) + + from astroplan import skymodel + + skymodel_1pt5airmass = skymodel(airmass=1.5) + + observer.skymodel = skymodel_1pt5airmass + # alternatively: >>> observer = Observer.at_site('apo', skymodel=skymodel_1pt5airmass) + skymodel_waveset, skymodel_flux = observer.skymodel + plt.plot(skymodel_waveset, skymodel_flux) + plt.show() + +And get the exposure time: + +.. code-block:: python + + >>> exptime_from_ccd_snr(snr, waveset_of_target, flux_of_target, observer, telescope) + 0.0026073144s + +.. note:: + The exposure time is slightly longer when accounting for the atmosphere than it was without an atmosphere due to slight atmospheric extinction in r' the band. diff --git a/docs/tutorials/index.rst b/docs/tutorials/index.rst index 7023a188..24798940 100644 --- a/docs/tutorials/index.rst +++ b/docs/tutorials/index.rst @@ -23,4 +23,5 @@ We currently have the following tutorials: plots periodic constraints - scheduling \ No newline at end of file + scheduling + exptime \ No newline at end of file diff --git a/setup.py b/setup.py index d61c2a51..e710b18a 100755 --- a/setup.py +++ b/setup.py @@ -108,7 +108,7 @@ docs=['sphinx_rtd_theme'], test=['pytest-astropy'], ), - tests_require=['pytest-astropy'], + tests_require=['pytest-astropy', 'synphot'], provides=[PACKAGENAME], author=AUTHOR, author_email=AUTHOR_EMAIL,