PSFs: circular and spatial methods#

Coronagraphic PSFs and the two ways the package images a target through them.

Circular method (cir_*): for rotationally symmetric optics and targets. A radial series of PSFs of sources along the +x axis (radial_psf_series, one per pixel of radius, out to the image’s edge; beyond it, out to the corners, the last one shifted) is weighted by the target’s radial profile, summed, and rotated about the star.

Spatial method (spatial_*): for any target. PSF kernels are sampled on a grid whose density follows how fast the PSF shape changes (spatial.equal_error_grid) and blended per pixel by spatial.spatial_psf. Kernels come either from direct propagation of every sample (‘direct’, any optics) or from x-axis PSFs interpolated in source position and rotated (‘rotate’, rotationally symmetric band-limited optics). spatial_error_report checks the model against exact PSFs.

Conventions (see toycoronagraph.optics): the star is at pixel coordinate c = (img_pixel-1)/2 on both axes; a source at pixel (row, col) is at x = (col-c)*pixel_scale, y = (row-c)*pixel_scale in lambda/D, with pixel_scale = 2*psf_range/img_pixel lambda/D per pixel. Intensities are in the package’s units (hcipy intensity / normal_factor), the same for frames, kernels and strips.

toycoronagraph.psf.radial_source_positions(img_pixel)[source]#

Source positions of the radial PSF series along the +x axis.

For odd img_pixel the star is on a pixel and the sources are at 0, 1, 2, … pixels. For even img_pixel the star lies between the four central pixels and source i is at sqrt(i^2+0.25) pixels (these radii lie between pixel centres).

Parameters:: img_pixel (int) – Number of pixels on one edge of the image.
Returns:: ((img_pixel+1)//2,) float64, distance of each radial PSF’s source from the star, in pixels.
Return type:: np.ndarray

toycoronagraph.psf.extended_frame_size(img_pixel)[source]#

Size of radial PSF frames that, rotated about the star, still cover the image corners.

Parameters:: img_pixel (int) – Number of pixels on one edge of the image.
Returns:: ceil(sqrt(2)*img_pixel)+8 pixels, plus 1 if needed so it has the parity of img_pixel (the frame then shares the image’s pixel lattice and centre).
Return type:: int

toycoronagraph.psf.radial_psf_series(optics, num_cores, frame_size=None)[source]#

PSFs of sources along the +x axis (at radial_source_positions).

Parameters:

optics (Optics) – The optical system and image sampling.
num_cores (int) – Number of CPU cores to use for parallel calculation (joblib n_jobs).
frame_size (int) – None for frames of the image size (float64); otherwise the size (same parity as the image) of larger float32 frames centred on the star, e.g. extended_frame_size(img_pixel) so rotated PSFs cover the image corners.

Returns:

((img_pixel+1)//2, size, size) array of PSFs (size = img_pixel or frame_size), in the package’s intensity units.

Return type:

np.ndarray

Raises:

ValueError – If frame_size and img_pixel differ in parity.

toycoronagraph.psf.coronagraph_throughput(optics, num_cores, radii)[source]#

Throughput of point sources on the +x axis: their total light relative to the same sources without the starlight-suppression element (for the vortex: the power through the Lyot stop with and without the vortex; by Parseval this is the whole focal plane’s light, independent of the image frame).

Parameters:

optics (Optics) – The optical system and image sampling.
num_cores (int) – Number of CPU cores to use for parallel calculation.
radii (array-like) – Source distances from the star, in pixels.

Returns:

(len(radii),) float64, throughput (0-1) at each radius.

Return type:

np.ndarray

toycoronagraph.psf.tol_for_error_target(error_target, blend_share=0.2)[source]#

Grid tolerance whose blending error fits blend_share of the image error target.

Uses the measured law (largest blended-PSF shape error ~ 1.5 tol^1.6, vortex at 16 px per lambda/D); the error report checks the result for any optics. The kernel crop gets the other 80% of the target (see choose_kernel_size).

Parameters:

error_target (float) – Target relative error of the final image.
blend_share (float) – Share of the target given to blending.

Returns:

tol (L1 shape change of unit-sum kernels between neighbouring samples), clipped to [0.003, 0.05].

Return type:

float

toycoronagraph.psf.axis_products(optics, num_cores, kernel_size, cache_path=None)[source]#

Runs for sources on the +x axis, shared by every spatial PSF grid of rotationally symmetric, band-limited optics.

Sources are spaced at half the Nyquist spacing of the intensity’s band limit (so band-limited interpolation between them is exact), from the star to beyond the frame corners. Each run gives the throughput; runs beyond the reach of the radial frames also give a “far strip”: the intensity on the block of the image pixel lattice (possibly outside the image) that interpolation windows around nearby radii need. These depend only on the optics, the image sampling and the kernel size, not on the grid tolerance, so they can be cached.

Parameters:

optics (Optics) – The optical system and image sampling.
num_cores (int) – Number of CPU cores to use for parallel calculation.
kernel_size (int) – Kernel width in pixels (odd).
cache_path (str) – If given, load from / save to this .npz (reused when the settings match, otherwise recomputed and overwritten).

Returns:

axis_x ((n,) float64, source positions in pixels from the star), throughput ((n,) float64), far_x ((m,) float64, the positions with a far strip), far_strips ((m, 2wh+3, width) float64, intensity on rows far_row0..far_row0+2wh+2 and columns far_col0[k]..far_col0[k]+width-1 of the image pixel lattice), far_col0 ((m,) int), far_row0 (int), and the settings (optics.settings(), kernel_size, cache_version); wh is from _rotate_geometry.

Return type:

dict

toycoronagraph.psf.shape_change_rates(window, optics, kernel_size, num_cores=-1)[source]#

How fast the PSF shape changes with source radius and angle, for rotationally symmetric optics (from x-axis PSFs).

d_r(r) is the L1 difference between unit-sum kernels at neighbouring radii per pixel of radius, and d_phi(r) the L1 difference between the kernel at (r, 0) and at (r, 2 deg) per radian. Both are smoothed with a running maximum (conservative). Radii are 0.5 px apart out to 3 lambda/D (or the reach of the radial frames) and at the axis node spacing beyond.

Parameters:

window (callable) – From _window_builder.
optics (Optics) – The optical system and image sampling.
kernel_size (int) – Kernel width in pixels (odd).
num_cores (int) – Number of threads.

Returns:

r_metric ((n,) radii in pixels, out to past the image corners), d_r ((n,) L1 change per pixel of radius), d_phi ((n,) L1 change per radian); float64.

Return type:

tuple

toycoronagraph.psf.shape_change_rates_direct(optics, kernel_size, num_cores=-1, n_dirs=4)[source]#

How fast the PSF shape changes with source radius and angle, for any optics (from PSFs computed directly along several directions).

Same quantities as shape_change_rates, measured along n_dirs directions (off the pixel axes) and taking the largest, so the grid is dense enough in every direction. Radii are 1 px apart out to 3 lambda/D and 0.25 lambda/D beyond.

Parameters:

optics (Optics) – The optical system and image sampling.
kernel_size (int) – Kernel width in pixels (odd).
num_cores (int) – Number of CPU cores for the propagation.
n_dirs (int) – Number of directions.

Returns:

r_metric ((n,) radii in pixels), d_r ((n,) L1 change per pixel of radius) and d_phi ((n,) L1 change per radian), the maxima over the directions; float64.

Return type:

tuple

toycoronagraph.psf.choose_kernel_size(optics, pre_img=None, error_target=0.02, num_cores=-1, n_probe=48, crop_share=0.8, seed=0, verbose=True, radial_psfs=None)[source]#

The smallest kernel size whose crop keeps the image error within a target.

A kernel of K x K pixels misses each PSF’s light that lands on the image outside the kernel window. That loss is measured, not modelled: PSFs of the actual optics are computed on the whole image at n_probe pixels drawn where pre_img is bright, and for every odd K the light on the image outside the K x K window around each probe is summed. Weighted by each probe’s light on the image, this is the crop error of pre_img’s final image. The smallest K keeping it within crop_share * error_target is returned; the rest of the target is left for the other error sources (blending, transmission map), which the error report checks afterwards.

With radial_psfs (circularly symmetric optics), a probe’s PSF is the radial PSF rotated to it, as the circular method places planets (_cir_planet): no new propagation. Against propagated probes the kernel size was the same in 207 of 216 cases (4-16 px per lambda/D, charges 2 and 6, rings, arcs, 1/r^3 profiles) and one step (2 px) apart otherwise, where the true crop error sits at the threshold (within 1% of it). Probes beyond the last radial PSF (the image corners) are still propagated.

Parameters:

optics (Optics) – The optical system and image sampling.
pre_img (np.ndarray) – (img_pixel, img_pixel) target image (any unit; its absolute value weighs the probes); None, or an all-zero image, weighs every pixel equally.
error_target (float) – Target relative error of the final image (e.g. 0.02 for 2%).
num_cores (int) – Number of CPU cores for the propagation.
n_probe (int) – Number of probe pixels drawn (with replacement; repeats count several times).
crop_share (float) – Share of the error target given to the crop.
seed (int) – Random seed of the probe pixels.
verbose (bool) – Print the choice.
radial_psfs (np.ndarray or None) – ((img_pixel+1)//2, S, S) radial PSF series of circularly symmetric optics, on frames of S >= extended_frame_size(img_pixel) (smaller frames would lose the wings rotated into the image corners: then every probe is propagated); None: propagate every probe.

Returns:

the kernel size (int, odd, pixels; 2*img_pixel-1 if no size meets the target), and a dict with kernel_sizes ((img_pixel,) int: 1, 3, …, 2*img_pixel-1) and crop (their estimated crop error, a fraction of the image light).

Return type:

tuple

toycoronagraph.psf.spatial_psf_calculation(optics, num_cores=-1, radial_psfs=None, kernel_size=None, error_target=0.02, pre_img=None, tol=None, mode=None, axis_cache=None, save_path=None, validate=True, n_validate=6, max_kernel_gb=4.0)[source]#

Calculates the spatially-sampled PSF kernels, sample grid and blend weights for the spatial convolution method.

The sample grid (equal_error_grid) is built from how fast the PSF shape actually changes, so neighbouring samples differ by about tol; inside ~1-3 lambda/D every pixel is its own sample. Each kernel is centred on its source and divided by the source’s throughput, so it carries the PSF shape only; the throughput goes into the transmission map tmap: each pixel’s own throughput for symmetric, band-limited optics, otherwise the samples’ throughput blended onto the pixels like the kernels. Their product reproduces the PSF at every sample.

Parameters:

optics (Optics) – The optical system and image sampling (see toycoronagraph.optics).
num_cores (int) – Number of CPU cores to use for parallel calculation.
radial_psfs (np.ndarray) – ((img_pixel+1)//2, S, S) radial PSF series ( radial_psf_series; S >= img_pixel, larger frames are cropped); computed if needed and not given. Used for rotationally symmetric, band-limited optics, to measure the grid’s rates and by ‘rotate’.
kernel_size (int) – Kernel width in pixels (odd); None chooses the smallest one meeting error_target for pre_img (see choose_kernel_size).
error_target (float) – Target relative error of the final image (e.g. 0.02 for 2%).
pre_img (np.ndarray) – (img_pixel, img_pixel) target image, used only to choose the kernel size.
tol (float) – Target PSF-shape change (L1 of unit-sum kernels) between neighbouring samples; None derives it from error_target (see tol_for_error_target).
mode (str) – ‘rotate’ interpolates the x-axis PSFs in source position (band-limited weights) and rotates them (quintic spline): fast, only for optics that are rotationally symmetric with a known band limit. ‘direct’ propagates every sample (any optics, much slower). None picks ‘rotate’ when the optics allow it.
axis_cache (str) – .npz path to cache the x-axis runs (see axis_products); None for no cache.
save_path (str) – If given, save the result there with np.savez.
validate (bool) – Check the result against exact PSFs and print the error analysis per stage (see spatial_error_report); stored in the result.
n_validate (int) – Samples and pixels checked per radius band.
max_kernel_gb (float) – Refuse to build a kernel cube larger than this (GB).

Returns:

psfs ((n, K, K) float32 kernels centred on their sources, divided by their throughput; 0 where the throughput is <= 1e-12), throughput ((n,) float64), tmap ((img_pixel, img_pixel) float64 transmission map), the grid of equal_error_grid (positions ((n, 2) (row, col) pixels), rings, ring_radii, ring_theta0, ring_ntheta, r_pix, and the blend weights inds, wts ((img_pixel, img_pixel, 4))), the shape-change rates r_metric, d_r, d_phi, the settings (optics.settings(), mode, kernel_size, tol, error_target) and, if validated, the err_* entries of spatial_error_report.

Return type:

dict

Raises:

ValueError – If mode is unknown, ‘rotate’ is not allowed for the optics, or kernel_size is even.
MemoryError – If the kernel cube would exceed max_kernel_gb.
RuntimeError – If ‘direct’ kernels are not centred on their sources (tilt convention check).

toycoronagraph.psf.load_spatial_psfs(path)[source]#

Load spatial PSF data saved by spatial_psf_calculation (or axis_products).

Parameters:: path (str) – Path of the .npz file.
Returns:: The saved arrays, with 0-d entries (scalar settings) converted back to Python values.
Return type:: dict

toycoronagraph.psf.spatial_error_report(model, num_cores=-1, n_per_band=6, seed=0, verbose=True)[source]#

Error analysis of the spatial PSF model, stage by stage, against exact PSFs (propagated directly through the model’s optics).

Checks a few random samples in each radius band (per-pixel zone, r_pix-3, 3-8 and >8 lambda/D; empty bands are dropped), and pixels spread evenly in radius over each band, so the image estimate can follow the error with radius. Errors are L1 differences of the PSF on the image, relative to the PSF’s light that lands on the image (light falling outside the image is not part of it, so it is not an error):

1. kernel generation: a sample's kernel x throughput vs the exact PSF at
   the sample, over the kernel window;
2. blending between samples: the blended kernel (with the exact throughput)
   vs the exact PSF of a source at the pixel;
3. transmission map: its error, times the light in the kernel;
4. kernel crop: the PSF's light on the image but outside the kernel window;
5. cumulated: the pixel's effective PSF (transmission map x blended kernel,
   nothing outside the window) vs the exact PSF, end to end.

Stage 2 includes the kernel generation error of the neighbouring samples, and the cumulated error is at most about the sum of stages 2-4.

Parameters:

model (dict) – From spatial_psf_calculation (kernels, grid, tmap and settings).
num_cores (int) – Number of CPU cores for the hcipy runs.
n_per_band (int) – Samples and pixels checked per radius band.
seed (int) – Random seed of the choice of samples and pixels.
verbose (bool) – Print the table (largest error per stage and band).

Returns:

err_edges ((n_band+1,) band edges, pixels), err_stages ((5,) stage names), err_max and err_median ((5, n_band) relative errors; NaN where a band has no check), err_n ((n_band,) pixels checked per band), err_pix_r ((n_pix,) the checked pixels’ radii, pixels) and err_pix ((5, n_pix) their errors per stage; NaN for the kernel generation, which is checked at samples).

Return type:

dict

toycoronagraph.psf.spatial_psf_planets(planets_pos, planet_brightness, psf_scale, model)[source]#

Add one or more planets to an image with exact coronagraphic PSFs.

Each planet’s PSF is propagated through the model’s optics at its exact position (no interpolation), on the full image frame.

Parameters:

planets_pos (list) – Planet positions (radius in arcsec, angle in radian from +x towards +y).
planet_brightness (list) – Planet brightnesses, in Jy/arcsec^2 (as Target stores them).
psf_scale (float) – The scale of the PSF, in arcsec per pixel.
model (dict) – From spatial_psf_calculation (gives the optics settings).

Returns:

(img_pixel, img_pixel) float64 image of the planets (Jy/arcsec^2); planets outside the image are skipped with a message.

Return type:

np.ndarray

toycoronagraph.psf.exact_point_sources(optics, pixels, values)[source]#

Point sources through the optics, each propagated at its exact position on the full frame.

Parameters:

optics (Optics) – The optical system and image grid.
pixels (list) – (row, column) of each source on the image, in pixels (need not be integers).
values (list) – Each source’s flux per pixel area (e.g. Jy/arcsec^2: Jy / pixel_scale^2).

Returns:

(img_pixel, img_pixel) float64 image of the sources, in the unit of values. Sources more than half a pixel outside the image are skipped with a message.

Return type:

np.ndarray

toycoronagraph.psf.spatial_psf_contrast(pre_img, planet_pos, planet_brightness, psf_scale, model, use_gpu=False, ncores=-2, report=True)[source]#

Find the planet brightness and background brightness in the final image, with spatially-varying PSFs.

Same measurement as cir_psf_contrast: the planet’s core (pixels above half its peak) is the aperture; the planet and the dust (with and without the target inside the IWA) are summed inside it.

Parameters:

pre_img (np.ndarray) – (img_pixel, img_pixel) input image, in Jy/arcsec^2.
planet_pos (array-like) – Planet position (radius in arcsec, angle in radian).
planet_brightness (float) – Planet brightness in Jy/arcsec^2.
psf_scale (float) – The scale of the PSF, in arcsec per pixel.
model (dict) – From spatial_psf_calculation.
use_gpu (bool) – Use CuPy for the FFT convolutions.
ncores (int) – Number of threads.
report (bool) – Print the estimated error of the disk image.

Returns:

planet brightness, background brightness, background brightness (ignored dust inside IWA) in Jy (floats: sums over the planet’s core aperture times psf_scale^2).

Return type:

tuple

toycoronagraph.psf.cir_psf_contrast(pre_img, r, planet_angle, planet_brightness, psf_scale, img_pixel, rot_number, psfs_name)[source]#

Find the planet brightness and background brightness in the final image.

The planet’s core (pixels above half its peak) is the aperture; the planet and the dust (with and without the target inside the IWA, see cir_iwa) are summed inside it.

Parameters:

pre_img (np.ndarray) – (img_pixel, img_pixel) circularly symmetric input image, in Jy/arcsec^2.
r (float) – the distance between the planet and the origin, in pixel unit.
planet_angle (float) – The angle between the x-axis and the line of the planet and the origin (radian).
planet_brightness (float) – The planet brightness, in Jy/arcsec^2.
psf_scale (float) – Scale of the PSF, in arcsec per pixel.
img_pixel (int) – Number of pixels on one edge of the image.
rot_number (int) – Number of rotations for generating the circular PSF.
psfs_name (str or np.ndarray) – The .npy file of the radial PSFs (frames of the image size, or larger), or the PSFs themselves (e.g. memory-mapped once: Coronagraph.radial_psfs).

Returns:

planet brightness, background brightness, background brightness (ignored dust inside IWA) in Jy.

Return type:

tuple

toycoronagraph.psf.cir_psf(pre_img, planets_pos, planet_brightness, psf_scale, iwa_ignore, add_planet, img_pixel, rot_number, psfs_name)[source]#

Calculates the final image of a circular symmetric pre-image through circular symmetric PSF with added planets.

Each radial PSF is weighted by the light it stands for (see _cir_chunks: the profile times the annulus area, and near the star the pixels’ own light), and the PSFs are rotated with cubic interpolation (not clipped, so the disk image is linear in pre_img). With radial PSF frames larger than the image (radial_psf_series with frames of extended_frame_size) the rotated PSFs also fill the image corners; with frames of the image size, the corners (beyond the inscribed circle) miss part of their light.

Parameters:

pre_img (np.ndarray) – (img_pixel, img_pixel) circularly symmetric input image (any unit, e.g. Jy/arcsec^2).
planets_pos (list) – Planet positions (radius in arcsec, angle in radian).
planet_brightness (list) – Planet brightnesses, in the unit of pre_img.
psf_scale (float) – Scale of the PSF, in arcsec per pixel.
iwa_ignore (bool) – Whether to ignore the target inside the IWA (see cir_iwa).
add_planet (bool) – Whether to add planets to the PSF.
img_pixel (int) – Number of pixels on one edge of the image.
rot_number (int) – Number of rotations for generating the circular PSF.
psfs_name (str or np.ndarray) – The .npy file of the radial PSFs (frames of the image size, or larger), or the PSFs themselves (e.g. memory-mapped once: Coronagraph.radial_psfs).

Returns:

(img_pixel, img_pixel) float64 final image, in the unit of pre_img. Planets beyond the last radial PSF are skipped with a message.

Return type:

np.ndarray

toycoronagraph.psf.cir_psf_planets(planets_pos, planet_brightness, psf_scale, img_pixel, psfs_name)[source]#

Add one or more planets to an image through a series of circular symmetric PSFs.

Parameters:

planets_pos (list) – Planet positions (radius in arcsec, angle in radian).
planet_brightness (list) – Planet brightnesses, e.g. in Jy/arcsec^2 (the image has the same unit).
psf_scale (float) – The scale of the PSF, in arcsec per pixel.
img_pixel (int) – The number of pixels on one edge of the image.
psfs_name (str or np.ndarray) – The .npy file of the radial PSFs (frames of the image size, or larger), or the PSFs themselves (e.g. memory-mapped once: Coronagraph.radial_psfs).

Returns:

(img_pixel, img_pixel) float64 image of the planets. Planets beyond the last radial PSF are silently skipped.

Return type:

np.ndarray

toycoronagraph.psf.core_throughput_fun(psfs)[source]#

Calculate the core throughput of a series of PSFs.

Parameters:: psfs (numpy.ndarray) – (n, S, S) normalized PSFs (e.g. the radial PSF series).
Returns:: (n,) float64 core throughput of the psfs: the sum of each PSF’s pixels above half its peak.
Return type:: numpy.ndarray

toycoronagraph.psf.cir_iwa(psfs)[source]#

Calculate the inner working angle (IWA) of a series of circular symmetric PSFs.

The IWA is the index of the first radial PSF whose core throughput reaches half of the largest one (np.searchsorted, assuming the core throughput increases with radius).

Parameters:: psfs (numpy.ndarray) – (n, S, S) radial PSFs (see radial_psf_series).
Returns:: The IWA of the PSFs, as a radial PSF index (~ distance from the star in pixels).
Return type:: int

toycoronagraph.psf.cir_core_throughput_plot(psfs, plot_dpi=300, pixel_scale=None)[source]#

Plot the core throughput of a series of circular symmetric PSFs, with the IWA marked.

Parameters:

psfs (numpy.ndarray) – (n, S, S) radial PSFs (see radial_psf_series).
plot_dpi (int) – Dots per inch (DPI) for the plot.
pixel_scale (float) – Pixel size in lambda/D (e.g. Optics.pixel_scale): the distance axis in lambda/D; None: in pixels (the index of each radial PSF).

Returns:

None. Saves Core_throughput.png in the current directory and shows the plot.