Source code for pynpoint.util.limits

"""
Functions for calculating detection limits.
"""

import math

from typing import Tuple

import numpy as np
from scipy.stats import t

from photutils.aperture import aperture_photometry, CircularAperture
from typeguard import typechecked

from pynpoint.util.analysis import student_t, fake_planet,\
    compute_aperture_flux_elements
from pynpoint.util.image import polar_to_cartesian, center_subpixel
from pynpoint.util.psf import pca_psf_subtraction
from pynpoint.util.residuals import combine_residuals


[docs]@typechecked def contrast_limit(path_images: str, path_psf: str, noise: np.ndarray, mask: np.ndarray, parang: np.ndarray, psf_scaling: float, extra_rot: float, pca_number: int, threshold: Tuple[str, float], aperture: float, residuals: str, snr_inject: float, position: Tuple[float, float]) -> Tuple[float, float, float, float]: """ Function for calculating the contrast limit at a specified position for a given sigma level or false positive fraction, both corrected for small sample statistics. Parameters ---------- path_images : str System location of the stack of images (3D). path_psf : str System location of the PSF template for the fake planet (3D). Either a single image or a stack of images equal in size to science data. noise : numpy.ndarray Residuals of the PSF subtraction (3D) without injection of fake planets. Used to measure the noise level with a correction for small sample statistics. mask : numpy.ndarray Mask (2D). parang : numpy.ndarray Derotation angles (deg). psf_scaling : float Additional scaling factor of the planet flux (e.g., to correct for a neutral density filter). Should have a positive value. extra_rot : float Additional rotation angle of the images in clockwise direction (deg). pca_number : int Number of principal components used for the PSF subtraction. threshold : tuple(str, float) Detection threshold for the contrast curve, either in terms of 'sigma' or the false positive fraction (FPF). The value is a tuple, for example provided as ('sigma', 5.) or ('fpf', 1e-6). Note that when sigma is fixed, the false positive fraction will change with separation. Also, sigma only corresponds to the standard deviation of a normal distribution at large separations (i.e., large number of samples). aperture : float Aperture radius (pix) for the calculation of the false positive fraction. residuals : str Method used for combining the residuals ('mean', 'median', 'weighted', or 'clipped'). snr_inject : float Signal-to-noise ratio of the injected planet signal that is used to measure the amount of self-subtraction. position : tuple(float, float) The separation (pix) and position angle (deg) of the fake planet. Returns ------- float Separation (pix). float Position angle (deg). float Contrast (mag). float False positive fraction. """ images = np.load(path_images) psf = np.load(path_psf) # Cartesian coordinates of the fake planet yx_fake = polar_to_cartesian(images, position[0], position[1]-extra_rot) # Determine the noise level noise_apertures = compute_aperture_flux_elements(image=noise[0, ], x_pos=yx_fake[1], y_pos=yx_fake[0], size=aperture, ignore=False) t_noise = np.std(noise_apertures, ddof=1) * \ math.sqrt(1 + 1 / (noise_apertures.shape[0])) # get sigma from fpf or fpf from sigma # Note that the number of degrees of freedom is given by nu = n-1 with n the number of samples. # See Section 3 of Mawet et al. (2014) for more details on the Student's t distribution. if threshold[0] == 'sigma': sigma = threshold[1] # Calculate the FPF for a given sigma level fpf = t.sf(sigma, noise_apertures.shape[0] - 1, loc=0., scale=1.) elif threshold[0] == 'fpf': fpf = threshold[1] # Calculate the sigma level for a given FPF sigma = t.isf(fpf, noise_apertures.shape[0] - 1, loc=0., scale=1.) else: raise ValueError('Threshold type not recognized.') # Aperture properties im_center = center_subpixel(images) # Measure the flux of the star ap_phot = CircularAperture((im_center[1], im_center[0]), aperture) phot_table = aperture_photometry(psf_scaling*psf[0, ], ap_phot, method='exact') star = phot_table['aperture_sum'][0] # Magnitude of the injected planet flux_in = snr_inject*t_noise mag = -2.5*math.log10(flux_in/star) # Inject the fake planet fake = fake_planet(images=images, psf=psf, parang=parang, position=(position[0], position[1]), magnitude=mag, psf_scaling=psf_scaling) # Run the PSF subtraction _, im_res = pca_psf_subtraction(images=fake*mask, angles=-1.*parang+extra_rot, pca_number=pca_number) # Stack the residuals im_res = combine_residuals(method=residuals, res_rot=im_res) flux_out_frame = im_res[0, ] - noise[0, ] # Measure the flux of the fake planet after PCA # the first element is the planet flux_out = compute_aperture_flux_elements(image=flux_out_frame, x_pos=yx_fake[1], y_pos=yx_fake[0], size=aperture, ignore=False)[0] # Calculate the amount of self-subtraction attenuation = flux_out/flux_in # the throughput can not be negative. However, this can happen due to numerical inaccuracies if attenuation < 0: attenuation = 0 # Calculate the detection limit contrast = (sigma*t_noise + np.mean(noise_apertures))/(attenuation*star) # The flux_out can be negative, for example if the aperture includes self-subtraction regions if contrast > 0.: contrast = -2.5*math.log10(contrast) else: contrast = np.nan # Separation [pix], position angle [deg], contrast [mag], FPF return position[0], position[1], contrast, fpf