Speckle Nulling

Speckle nulling utilizes two main functions, both located in gpilib/nulling/: nulling.pro and take_ifs_image.pro.


nulling.pro is called with an option input with four possible values:

  1. Create a new calibration
  2. Generate calibration from files on disk
  3. Run nulling.
  4. Simulate nulling from files on disk.


For options 0 and 3, you must set the imstart and imdat keywords to specify which files are used.

The only other required input is band, which specifies which coronagraph you are using and must be one of: Y, J, H, K1, or K2.

The optional keywords are split between calibration and nulling modes. There are only four common keywords between the two modes:

Keyword Description
intTime Integration time to use (seconds). Ignored for option = 1 and 3
refslice Cube slice to use. If unset, defaults to slice 9 (10th of 37). If set to -1, speckle aligns and sums all slices, keeping the speckles at their location in the default slice (9).
imdat Date string of image directory (i.e., ‘120618’). Ignored for option = 0 and 2
imstart Integer value of first image in set. Ignored for option = 0 and 2

Calibration Mode Keywords

There are three separate calibrations, described in detail in the next section. They are toggled by:

Keyword Description
/redolocations Redo the spatial frequency calibration
/redomatrix Redo the intensity calibration
/redo_scicam_values Redo image geometry calibration

Redoing the geometry calibration updates multiple internal values and requires that the other two calibrations be done. The most important of these values, available as optional outputs, are:

Optional Output Description
gpi_scicam_pix_to_ripple Number of pixels between PSF center to the MEMS ripple, along a major axis in the reference frame of the science image (i.e., horizontal, not diagonal distance)
gpi_scicam_offset Offset between PSF center and image center
gpi_scicam_rotation Rotation between image and PSF reference frames
gpi_scicam_satspot_cens 4x2 array of locations of dark hole satellite spot (in pixels in derotated, centered image)

If /redo_scicam_vals is set, an initial image will be taken with a waffle pattern applied to the MEMS. You will be asked to check the resulting value for gpi_scicam_pix_to_ripple, which should be around 130 pixels in H band. Running this should only be necessary if something has changed in the physical system, or you are switching the band or default slice to use.

All of these are calculated by take_ifs_image.pro when the /calcall flag is set. These values, and all of the calibrations described below, are saved in the GPI IDL private datapath as described below.

Nulling Mode Keywords

Keyword Description
phase_flat Previously found best phase. If unset, this value defaults to all zeros
fraction Reciprocal of fraction of possible modes to fix (i.e., fraction=2 will fix 50% of all controllable modes)
maxmodes Maximum number of modes to consider (overrides fraction keyword)
iterations Number of iterations
ignore_center Minimum radial spatial frequency to correct (defaults to 7)
masked_region Two element array of minimum and maximum radial spatial frequencies to correct (defaults to [10,20], overrides ignore_center)
spurious Inject fake speckle
gain Loop gain (defaults to 0.5)
integrator Loop integrator (defaults to 0.99)


  1. Generate a new calibration (will overwrite the old one - careful!):
    IDL> nulling, 0, 'H' , /redomatrix,/redolocations,/redo_scicam_vals,intTime=30
  2. Redo a calibration using files from disk:
    IDL> nulling, 1, 'H', imdat='120618',imstart=87,/redomatrix,/redolocations
  3. Do a run
    IDL> nulling, 2, 'H', constrained_region=[10,20],iterations=10,gain=0.5,intTime=30.
  4. Do a run reusing previous flat
    IDL> final_phase = readfits('/home/LAB/IDL/gpilib/private/NUL_HIST/130711/H_220714/NUL_results_iteration_0_final_phase.fits')
    IDL> nulling, 2, 'H', ignore_center=7,iterations=2,maxmodes=100,gain=0.5,intTime=60.,phase_flat=final_phase


The speckle nulling procedure requires three calibrations:

  • Image geometry
  • Pixel location to spatial frequency
  • Image intensity to amplitude


New calibrations overwrite the current defaults for the band. If you need to restore archived calibrations, see the Files section below.

Image Geometry

The nulling code assumes that all input images are centered and derotated such that the PSF core lands on the center pixel (140,140), and that the edges of the dark hole are aligned with the edges of the image. To do this, an open loop image is taken with a waffle pattern on the MEMS DM (control at the highest spatial frequency - neighboring actuators are poked in opposite directions), which produces waffle spots in the corners of the dark hole. These, along with the satellite spots, provide all of the geometric information needed.


Figure: Slice of reconstructed open-loop cube with waffle pattern on MEMS. gpi_scicam_offset is defined as the difference between the PSF core and the center of the image, marked by the red and blue points, respectively.

The waffle and satellite spots are identified using the find_sat_spots code from the GPI pipeline. The dark hole satellite spots are identified first and masked out and the satellite spots outside the dark hole are masked with a Hanning window. The waffle spots are then identified using the masked image. The geometric values are then calculated as follows:

Value Calculation
PSF core location Mean of the waffle spot locations
gpi_scicam_pix_to_ripple Mean of the distances between each waffle spot and the PSF core, divided by the square root of 2 (to convert from diagonal to horizontal distance)
gpi_scicam_offset PSF core location minus (gpi_scicam_n - 1)/2, where gpi_scicam_n is the number of image pixels (281).
gpi_scicam_rotation Mean of the angle of the four vectors between the waffle spots and PSF core
gpi_mems_rotation Mean of the angle of the four vectors between the satellite spots and PSF core


The definition of the scicam offset is ‘neagtive’ so that positive values require shifts to the left and down. The scicam and mems rotation angles are defined such that a positive value is a clockwise rotation. These two angles should be very simillar, differing by only about 3 degrees for the H band apodizer.


Figure: Slice of reconstructed open-loop cube with waffle pattern on MEMS after derotation and centering. The length denoted by the arrow is equal to gpi_scicam_pix_to_ripple/2.

Spatial Frequency

As GPI’s MEMS has 48x48 actuators, we can potentially control up to 1152 spatial frequencies. We do not attempt to control any of the highest or lowest spatial frequencies, leaving 1102 controllable modes. These are enumerated using two coordinates, (k,l), and are mapped to a 2D grid in variable kl_valid:


Figure: Visualization of variable kl_valid. White areas are controllable. Coordinate (0,0) is in the bottom, right-hand corner.

When dealing with whole annular regions, we also use ‘radial spatial frequencies’, equal to \(\sqrt{k^2 + l^2}\).

Internally, these modes are also addressed by a single index (from 0 to 1101). To aid in converting between these two indexing systems, we carry two 1D arrays, backout_k and backout_l, whose entries map the 1D index to (k,l), and a 2D matrix forward_kl whose elements represent the 1D index of the controllable modes. Note that there is a slight ambiguity in the matrix, as all uncontrollable modes have value 0, but this does not cause any problems as mode (0,0) is set as uncontrollable.

In addition to the fundamentally uncontrollable modes, we may occasionally only wish to use a subset of controllable modes (for instance when focusing on a specific annulus or region of the focal plane). To aid in this, we define two arrays (of the same size as backout_k and backout_l), ignore_center_mask_long and constrained_region_long, which contain boolean values for which of the controllable modes may be used. The specific regions are set by inputs ignore_center (minimum radial spatial frequency) or constrained_region (minimum and maximum radial spatial frequencies).

In particular, we never want try to correct the satellite spots (which will be picked up by the code as particularly bright speckles), so the spatial frequencies corresponding to all of the locations within 9 pixels of the points in gpi_scicam_satspot_cens will be ignored:


Figure: Visualization of variable kl_valid with spatial frequencies corresponding to satellite spot locations masked.

Finally, to minimize cross-talk, we do not want to controll neighboring modes (i.e., if mode (4,4) is being controlled, we do not want to control any modes where k = 4 and l \(\in [3,5]\) and vice versa). To aid in this, we store two more arrays: modes_not_allowed_mask and modes_not_allowed_mask2. These have dimension 48x28x1102, and are 2D binary masks for each of the controllable modes. The second mask blocks neighbor-of-neighbor spatial frequencies as well as nearest neighbors.

In order to convert between pixel locations and spatial frequencies, we use arrays maskhandling_indices, maskhandling_values and maskhandling_lens. maskhandling_indices is an 1102xN array, where N is the largest number of pixels corresponding to any one spatial frequency. The columns of this array map the controllable frequencies to image pixels, represented by a 1D index into a gpi_scicam_n x gpi_scicam_n array. The maskhandling_values array is laid out in the same fashion, but contains weights on the individual pixels representing a 2D Gaussian spot at that spatial frequency on the image plane. maskhandling_lens is an 1102x1 array that stores the length of each column in the first two arrays to use (all of the columns are zero-padded to the size of the largest one).

Intensity Calibration

When generating new calibrations, the code will take a series of images with pure sinusoids applied at increasing spatial frequencies, such that k=l, producing a series of pairs of speckles moving outward radially from the center:


Figure: Slices of reconstructed closed-loop cubes with pure sinusoids of increasing spatial frequency applied. You can see the resulting bright speckles moving outwards from the center.

A reference image (with a flat DM shape) will be subtracted from each image, and the intensity of the added speckles recorded. The code will then attempt to fit the resulting intensities to a curve and present you with the results (first asking about the number of initial data points to drop - this should be around 3-4 for modes within the FPM). If you the fit looks bad, you can regenerate it, using the existing code for reference:

;; The radial spatial modes are stored in var klr2, intensity measurements in dr2 and 1 sigma errors in terrs2
IDL> baseval = median(dr2)
IDL> expr = 'P[0] + P[1]*exp(0.5*(X/P[2])^2.)'
IDL> res = MPFITEXPR(expr, klr2, dr2, terrs2 , [baseval,10,10])
IDL> new_response_kl = res[0] + res[1]*exp(0.5*(klr/res[2])^2.)

Once satisfied, a new calibration will be written:


Using this calibration two matrices will be generated: forward_matrix and inverse_matrix. As the names suggest, these will be used to map between image plane intensity to command amplitude.

The average flux of the dark hole satellite spots is recorded and stored along with these calibration matrices. When running the nulling mode, the average flux of the satellite spots in the baseline (pre-nulling) PSF is calculated, and the matrices are updated as:

forward_matrix *= (baseline_flux/cal_flux)
inverse_matrix /= (baseline_flux/cal_flux)


Nulling calibrations and pre-calculated values are stored in the gpi_idl_private_datapath. The code is set up to only store one set of calibrations per band at a time, but every calibration is archived and available for use.

Current calibrations are stored at the top level of gpi_idl_private_datapath, with file prefix NUL_private and file names:

Filenames Description
_BAND_gauss_maks_info_file?.fits The pixel to spatial frequency mapping arrays are stored in 9 files (3 each - indices, values, and lengths - for all the controllable frequencies, the nearest neighbor mask and the neighbor of neighbor mask)
_BAND_forward_matrix.fits Forward intensity mapping matrix
_BAND_inverse_matrix.fits Inverse intensity mapping matrix
_modes_not_allowed_mask?.fits The neighbor and neighbor of neighbor masks are stored in two files, indexed as 1 and 2, respectively
_kl_valid.sav IDL save file containing the kl_valid, backout_k, backout_l, and forward_kl arrays
_BAND_scicam_values.sav IDL save file containing all of the geometric information calculated by /redo_scicam_values
_BAND_final_phase.fits The best phase setting achieved by the latest nulling run


The values in NUL_private_BAND_scicam_values.sav are only used when redoing the spatial frequency and intensity calibrations without first redoing the geometric calibration. In normal nulling, these values are read from the header of the forward matrix fits file.

The spatial frequency mapping fits files include the following header values:

Keyword Calculation
RIPPLE gpi_scicam_pix_to_ripple - Horizontal distance from PSF core to MEMS ripple
OFFSET_X gpi_scicam_offset[0] - X offset between PSF core and image center
OFFSET_Y gpi_scicam_offset[1] - Y offset between PSF core and image center
AM_ROT gpi_scicam_rotation - Rotation of PSF in image
MEMS_ROT gpi_mems_rotation - Rotation of sat spots in image
REFSLICE refslice - Cube slice used, -1 for all
CEN_jk gpi_scicam_satspot_cens[j,k] - Location of spot in derotated cube

The inverse and forward matrix fits files contain the same information and a SATFLUX keyword giving the average sat spot flux in the calibration reference PSF. Their headers also contain HISTORY lines starting with ‘FILE:’ listing all of the files used in the calibration, starting with the reference PSF.

The archived calibrations and nulling results are stored in gpi_idl_private_datapath/NUL_HIST/. Each invocation of nulling.pro creates a new subdirectory of the form YYMMDD/BAND_HHMMSS (year,month,day/BAND_hour,minute,second). Thus, no results are ever overwritten, and even when you repeat a calibration (or run) using files already on disk (options 1 and 3), you will create a new results directory with the current time stamp. In these directories, the calibrations are written with the exact same filenames as at the top level, so to restore a previous calibration you can just copy everything to the top level:

> cp ./private/NUL_HIST/130101/121212/NUL_private* ./private/

Results are written with prefix NUL_results_. Every nulling run also generates a log file with all of the settings and a record of all of the IFS images generated. At the end of a run, the best phase setting achieved is copied to the top level _BAND_final_phase.fits file.