Source code for dataval.noise_model

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Noise model as a function of magnitude and position

.. codeauthor:: Rasmus Handberg <rasmush@phys.au.dk>
"""

import numpy as np
from astropy import units as u
from astropy.coordinates import SkyCoord
from scipy.interpolate import InterpolatedUnivariateSpline
from .utilities import mag2flux

#--------------------------------------------------------------------------------------------------
[docs] def ZLnoise(gal_lat): """ RMS noise from Zodiacal background. .. codeauthor:: Mikkel N. Lund <mikkelnl@phys.au.dk> """ rms = (16 - 10) * (gal_lat/90 - 1)**2 + 10 # e-1 / pix in 2sec integration return rms
#--------------------------------------------------------------------------------------------------
[docs] def Pixinaperture(Tmag, cad=1800): """ Number of pixels in aperture as a function of Tmag .. codeauthor:: Mikkel N. Lund <mikkelnl@phys.au.dk> """ # Approximate relation for pixels in aperture (based on plot in Sullivan et al.) #pixels = (30 + (((3-30)/(14-7)) * (Tmag-7)))*(Tmag<14) + 3*(Tmag>=14) if cad in (1800, 600, 200): masksize = np.array([ [2.05920002, 1484.5], [2.95159999, 715], [3.84399996, 447], [4.73639993, 282.5], [5.62879990, 185], [6.52119987, 126], [7.41359984, 98], [8.30599982, 76], [9.19839979, 61], [10.09079976, 49], [10.98319973, 38], [11.8755997, 28], [12.76799967, 20], [13.66039964, 14], [14.55279961, 8] ]) elif cad in (120, 20): masksize = np.array([ [2.48170001, 473], [3.56310005, 210], [4, 174], [5.72590014, 85], [6.80730019, 69], [7.88870023, 61], [8.97010028, 50], [10.05150032, 38], [11.13290037, 26], [12.5, 13], [15.0, 4] ]) else: raise NotImplementedError() # Interpolate linearly in log-space: pix_log = InterpolatedUnivariateSpline(masksize[:,0], np.log10(masksize[:,1]), k=1, ext=0) pix = lambda x: np.round(10**(pix_log(x)), decimals=13) # noqa: E731 pixels = pix(Tmag) # Ensure lower limit of 4 pixels: pixels = np.clip(pixels, 4, None) return np.asarray(pixels, dtype='int32')
#--------------------------------------------------------------------------------------------------
[docs] def phot_noise(Tmag, timescale=3600, coord=None, sysnoise=60, Teff=5775, cadpix=1800): """ Photometric noise model in ppm/timescale. Parameters: Tmag (ndarray or float): TESS magnitudes to calculate noise model for. timescale (float): Timescale of noise model in seconds. Defaults to 3600 (1 hour). coord (SkyCoord or dict, optional): Sky coordinate used for zodiacal noise. Defaults to RA=0, DEC=0. sysnoise (float, optional): Systematic noise contribution in ppm/hr. Defaults to 60. Teff (float, optional): Effective temperture. Defaults to 5775. cadpix (TYPE, optional): Cadence where images are taken. Choices are 1800 and 120. Defaults to 1800. Returns: tuple: - ndarray: Total noise model as a function of TESS magnitide. - ndarray: Individual noise components in 4 columns: Shot, Zodiacal, Read and Systematic. Raises: ValueError: On invalid coord input. .. codeauthor:: Mikkel N. Lund <mikkelnl@phys.au.dk> .. codeauthor:: Rasmus Handberg <rasmush@phys.au.dk> """ # Make sure Tmag is a numpy array: Tmag = np.asarray(Tmag) # Create SkyCoord object from the coordinates supplied: if coord is None: # TODO: Nothing is given, so use complete dummy coordinates gc = SkyCoord(0*u.degree, 0*u.degree, frame='icrs') elif isinstance(coord, dict): if 'RA' in coord and 'DEC' in coord: gc = SkyCoord(coord['RA']*u.degree, coord['DEC']*u.degree, frame='icrs') elif 'ELON' in coord and 'ELAT' in coord: gc = SkyCoord(lon=coord['ELON']*u.degree, lat=coord['ELAT']*u.degree, frame='barycentrictrueecliptic') else: raise ValueError("Invalid coord in dict format") elif not isinstance(coord, SkyCoord): raise ValueError("Invalid coord") # Calculate galactic latitude for Zodiacal noise gc_gal = gc.transform_to('galactic') gal_lat0 = gc_gal.b.deg gal_lat = np.arcsin(np.abs(np.sin(gal_lat0*np.pi/180)))*180/np.pi # Number of 2 sec integrations in cadence integrations = timescale/2 # Number of pixels in aperture given Tmag pixels = Pixinaperture(Tmag, cadpix) # noise values are in rms, so square-root should be used when factoring up Flux_factor = np.sqrt(integrations * pixels) # Mean flux level in electrons per cadence mean_level_ppm = mag2flux(Tmag) * timescale # electrons (based on measurement) #, Teff # Shot noise shot_noise = 1e6/np.sqrt(mean_level_ppm) # Read noise read_noise = 10 * Flux_factor * 1e6/mean_level_ppm # ppm # Zodiacal noise zodiacal_noise = ZLnoise(gal_lat) * Flux_factor * 1e6/mean_level_ppm # ppm # Systematic noise in ppm systematic_noise = np.full_like(Tmag, sysnoise / np.sqrt(timescale/3600)) # ppm / sqrt(hr) # Put individual components together in single table: noise_vals = np.column_stack((shot_noise, zodiacal_noise, read_noise, systematic_noise)) noise_vals = np.clip(noise_vals, 0, None) # Calculate the total noise model by adding up the individual contributions: total_noise = np.sqrt(np.sum(noise_vals**2, axis=1)) total_noise = np.clip(total_noise, 0, None) return total_noise, noise_vals # ppm per cadence