Category Archives: Uncategorized

Best friends from now on: shaped pupils and Lyot coronagraphs

How a shaped pupil solves the "classical" Lyot coronagraph problem of monochromatic field cancellation.
Figure from Zimmerman, et al., 2016: Illustration of how a shaped pupil solves the “classical” Lyot coronagraph problem of monochromatic field cancellation.

An article I authored with my colleagues from the Princeton high-contrast lab has been published by the Journal of Astronomical Telescopes, Instruments, and Systems: Shaped pupil Lyot coronagraphs: High-contrast solutions for restricted focal planes.

For 15+ years, Lyot coronagraphs and shaped pupils have evolved in parallel as competing tools for producing a zone of destructive interference next to a star’s image, with the goal of revealing orbiting exoplanets. After all this time, it turns out they are much better together than apart.

We learned this fact not a moment too soon. The instrument designs for WFIRST are getting locked in place rapidly, and that mission may be our last chance for many years to launch an ambitious space-based exoplanet imager. The new hybrid shaped pupil Lyot coronagraph offers a winning balance of performance, fabrication maturity, and error tolerance.

Generate an interactive gallery of sky survey cutouts from a list of targets

After filtering down some source catalog to a short list of interesting targets, it’s nice to have image cutouts handy, as either a sanity check, or as finder charts. Aladin Lite, hosted by CDS, has a javascript interface that allows you to display image cutouts from some major sky surveys. Here’s a short python script I wrote to generate a gallery of DSS and 2MASS cutouts, with hyperlinks to the full SIMBAD entry. As an input, it takes an ASCII file containing a list of 2MASS catalog IDs.

Here is an example gallery; the python source is below (hosted with gist).

Poster at 227th AAS meeting

AAS meeting, Jan 2016, Kissimmee, Florida
At 227th AAS meeting, Jan 5, 2016, Kissimmee, Florida. Photo by Heather D’Angelo.

Here I am at the American Astronomical Society meeting this week. I presented a poster about the first project I began after arriving at STScI, in collaboration with Marie Ygouf, Laurent Pueyo, Rémi Soummer, and others. I analyzed laboratory data from JPL’s High Contrast Imaging Testbed, where protoypes for the WFIRST coronagraph are being tested. The goal is to learn how well we can subtract residual speckle noise from the coronagraph images, to predict what kinds of planets we can detect.

Imaging the planet around GJ 504

Figure 2 from Skemer, Morley, Zimmerman, et al.: De-rotated co-add of all photometric quality GJ 504 images taken in the LNB7 (3.95 μm) filter. The planet, GJ 504 b, is the circled point source to the northwest of the star, at separation ∼2.′′5. Due to the large dynamic range of structure in the image, two regions are displayed with different flux scales: (1) within 2′′ of the star, the Airy rings are displayed with a logarithmic stretch; the outer region (2) that includes the planet is displayed with a linear stretch centered at zero. The corresponding grayscale bars on the right-hand side indicate the flux in units normalized to the peak of the unsaturated star’s PSF.

Andy Skemer, Caroline Morley, myself, and the rest of the LEECH Science Team posted a paper on the arXiv today: Characterization of the Coldest Directly Imaged Exoplanet, GJ 504 b, and Evidence for Super-Stellar Metallicity.

The article, soon to be published in ApJ, describes new observations of the planet candidate discovered to orbit GJ 504 a few years ago (Kuzuhara et al., 2013). By measuring the thermal emission from the companion in three custom filters in the range 3.6–4.1 microns, we’ve made new inferences about the temperature, density, and composition of the atmosphere. The observations were made possible with LMIRCam, an AO-fed infrared camera on the Large Binocular Telescope.

I was responsible for the entirety of the data reduction, taking steps familiar to other high-contrast observers in this wavelength regime:

  • Aligning and combing thousands of raw, short exposures, correcting for background sky noise and detector pattern noise
  • Angular differential imaging to subtract the star (here implemented with the KLIP algorithm)
  • Contrast curve calculation, PSF modeling, photometry, and error analysis

GJ 504 b is on the extreme faint end of all the planetary-mass companions imaged to date (as Figure 4 in the journal article illustrates). Fortunately, the combination of relatively high angular separation (2.5 arcsec from the host star) with the deep sensitivity of LMIRCam–as shown above in the co-add image–was high enough to secure a detection at each filter. The wavefront correction provided by the LBTI AO system was such that Airy rings spanning three orders of magnitude in intensity are evident, out beyond a 1-arcsec radius. Even without PSF subtraction, the planet is evident at a contrast of 1E-5, in the northwest corner of the co-add image.

There is a lot more we could learn from the infrared SED at wavelengths redward of 4 microns. For such a faint object, however, that part of the spectrum is not likely to be accessed from the ground. Therefore, in the future GJ 504 will make an excellent target for JWST NIRCam.

Sparse aperture mask wavefront sensor

An article I co-authored with Princeton grad student Hari Subedi has been published by the Journal of Astronomical Telescopes, Instruments, and Systems: Coronagraph-Integrated Wavefront Sensing with a Sparse Aperture Mask . A sparse aperture mask, similar to the kind now in common use for high-angular resolution astronomical imaging (a recent example being Hinkley et al. 2015), is also effective at converting phase aberrations into diverse intensity patterns. By intercepting, collimating, masking, and finally refocusing the starlight rejected from a coronagraph focal plane, we therefore have a simple approach to sensing dynamic low-order aberrations of the kind that are problematic for space telescope coronagraphs such as the one envisioned for WFIRST-AFTA.


A shaped pupil for a backyard telescope

Dr. Russ Genet at Cal Poly does speckle interferometry observations of binary stars. He came up with the idea to use a shaped pupil to improve the search zone around his primaries, and increase the kinds of targets accessible to his small telescope observing program. Our team agreed it would be an interesting experiment to see what shaped pupils can do for seeing-limited observations of binary stars.
I wrote some 2-D optimization programs (making yet more use out of the approach of Carlotti et al 2013) to come up with a design for Russ’ 11-inch Schmidt-Cassegrain pupil. There was one solution with a fairly simple layout. It improves the open aperture telescope PSF by suppressing the first sidelobe along one axis, as shown in the horizontal contrast cut plot below.C11_SP_12wa60_logc27_bowtie90degC11_contrast_cut_12wa60_logc27_bowtie90degIf Russ and his observing team, led by Cal Poly Master’s student Ed Foley, can freeze out the atmosphere and limit static aberrations, then they have a shot of effectively (almost) doubling the aperture of their Celestron setup.

naive simulated TPF coronagraph image

A student at Princeton is doing his senior thesis on simulated coronagraph data sets, to use as a testing ground for hypothesis testing of somewhat realistic (low SNR) planet signals.

Fomalhaut is the target star for the case study. For its distance of 7.61 pc, 1 AU and 5 AU projected separations correspond to 131.4 mas and 657.0 mas, respectively. For 8-meter telescope observing at 550 nm, \lambda/D = 14.2 mas, so maximum angular separations are 9.27 \lambda/D and 46.33 \lambda/D. Periastron for the inclined, circular 5 AU orbit is \cos(66^{\circ})\times 46.3 \lambda/D = 18.85 \lambda/D.
I handed off a Matlab Fourier propagation model of a concentric ring (spiderweb) coronagraph to my student. This design has a working angle range 6 - 60 \lambda/D to frame the system nicely. Also note that 1 AU projected separation at 1.34 pc distance of alpha Cen corresponds to 1/1.339 / 0.0142 = 52.6 \lambda/D.

This design surpasses 10^{-10} contrast for the ideal wavefront, so of course it’s fun to throw in the canonical Earth-like dot, shown below to the right of the star (whose PSF “saturates” the flux scale).

Floating point roundoff errors in Fourier optics calculations

Discrete Fourier transforms are essential for simulating light propagation in telescopes and coronagraphs. To cut down the computational burden, and maintain high spatial resolution, we can make use of well-established symmetry properties of Fourier transforms. For example, if we start with a real-valued electric field in the telescope pupil, with even symmetry about the y-axis, we should get an image plane electric field with even symmetry about the y-axis. Or so I thought, until a sanity check failed…

Here’s where my assumptions broke down. The above plot is the real part of the image formed by an ideal shaped telescope pupil. I computed this with a 2-D FFT in Matlab and shifted the result so that the zero-phase pivot is at the array center. Visually, it appears symmetric about the vertical meridian, so that \mathfrak{Re}\left\{ E(x, y) \right\} = \mathfrak{Re}\left\{ E(-x,y) \right\}. And that is indeed what we’d expect, transforming from a flat (zero-phase) pupil plane wavefront with the same even symmetry.

So when I take the right half of the image, and then flip and subtract it from the left half, I should get a wall of zeroes. Instead, disparities appear:


We have disagreements as high as 2000, corresponding to almost 1% of the peak field amplitude in the image plane. Initially, I was very puzzled by what might cause such large errors. Was it an error in my pupil definition, or some array offset creeping in? Do I not understand the rules of Fourier transforms, after all? After some tests, I figured out that in each image plane pixel I was accumulating a large number of tiny floating-point roundoff errors.

\displaystyle X [ k, l ] = \frac{1}{\sqrt{MN}}\sum_{m=0}^{M-1} \sum_{n=0}^{N-1} x[m, n] e^{-i2\pi(mk/M + nl/N)}

Consider the above textbook formula for the 2-dimensional discrete Fourier transform. There are countless applications, but to an astronomer it describes the electric field diffraction pattern in the image plane of a telescope when the incoming starlight is monochromatic, unpolarized, and the optics are perfect. We can break this into inner and outer sums:

\displaystyle X [ k, l ] = \frac{1}{\sqrt{MN}} \sum_{m=0}^{M-1} \left( \sum_{n=0}^{N-1} x[m, n] e^{-i2\pi nl/N} \right) e^{-i2\pi mk/M} \\ = \frac{1}{\sqrt{MN}} \sum_{m=0}^{M-1} P[m, l] e^{-i2\pi mk/M}

In our case, the telescope pupil x[m,n] has even circular symmetry along the x-axis, so that x[m,n] = x[m,N-n]. When that is the case, and x[m,n] is real, it can be shown that the inner sum, P[m,l], has circular even symmetry: P[m,l] = P[m,N-l]. Furthermore, so does the image: X[k, l] = X[k, N-l].

The problem is, in practice there are several orders of magnitude of separation between the mirrored spatial frequency pair l and N-l. Before the complex exponential can be evaluated, its phase argument must be stored as a floating point value, which entails some roundoff error. To investigate my case, I chose some sample points to test how the inner complex exponential term behaves for my array (with size = 3000 x 3000).

With N = 3000, l = 10, I noticed a particularly high roundoff error at n = 1844. Since I only consider the real part for now, I simply compare cos(2 \pi nl/N) for this mirrored pair of spatial frequencies:

>> wrapToPi(2*pi*n*l/N)
ans =
>> wrapToPi(2*pi*n*(N-l)/N)
ans =
>> wrapToPi(2*pi*n*l/N) + wrapToPi(2*pi*n*(N-l)/N)
ans =
>> cos(2*pi*l*n/N) - cos(2*pi*(N-l)*n/N)
ans =

Because the two wrapped angles are symmetric about zero radians, or \left( nl/N~\%~1 \right) + \left( n(N-l)/N~\%~1 \right) = 1, the difference of cosines is mathematically zero. What we actually compute is the best approximation possible with 64-bit double precision floating point. The residual comes out exactly the same in both MATLAB and Python/Numpy. By contrast, if I take the difference of another pair of cosine evaluations, but now with the higher angle wrapping around only one period beyond nl/N, the error is several orders of magnitude smaller:

>> wrapToPi(2*pi*n*l/N) - wrapToPi(2*pi*(n*l/N + 1))
ans =
>> cos(2*pi*n*l/N) - cos(2*pi*(n*l/N + 1))
ans =

Although that 1e-12 above seems tiny, since we are adding up N phasors in the inner sum, and then adding up M of those, roundoff errors are accumulated 9 million times for each image pixel. Just left of the image center (the l=0 pivot axis of the DFT), N-l is at its highest. Conversely, just right of the image center, the roundoff error fed into the trig function will always be very low.

Smaller array sizes would reduce this problem. But depending on what kind of measurement you’re making, that may not be an option. A better solution is to exploit the pupil’s symmetry and only keep the right hand side of the FFT result (spatial frequencies l = 0, 1, 2,\dotsc, N/2-1). That would maintain a self-consistent symmetry between the pupil and the image. In fact, because the full complex image field of a real pupil is Hermetian, in this particular case the full image image is determined by one quadrant; the imaginary part in the lower half reflection just needs to be negated.

That leads to the question: is their an FFT library for MATLAB or Python/Numpy that either automatically detects, or allows you to specify the symmetry properties of your input array? It seems like an obvious way to eliminate wasteful operations which are particularly prone to quantization noise.

Some useful references on DFTs around the web:

DIY sun projector


Thank you, everyone who took time to stop and look through the sun projector today at Communiversity. Mike and I were thrilled with all the enthusiasm we encountered. Our highest hope was to inspire a few of you to make your own solar telescope at home. To help you along, I promise to post detailed instructions here this week. So please stay tuned, and in the mean time, feel free to pose questions through the Contact page.