This document is not intended as comprehensive documentation for the poppy
package. Please see the documentation for full details in the API, and the SPIE papers by Perrin et al. (2012, 2014) for further background information and discussion of algorithms. However this notebook should be enough to get you started.
Poppy implementsan object-oriented system for modeling physical optics propagation with diffraction, particularly for telescopic and coronagraphic imaging. Right now only image and pupil planes are supported using Fraunhofer (far-field) diffraction; Fresnel modeling of intermediate planes is a future goal.
You should be able to install poppy directly from PyPI in the usual manner:
pip install poppy --upgrade
import poppy
%pylab inline --no-import-all
matplotlib.rcParams['image.origin'] = 'lower'
Populating the interactive namespace from numpy and matplotlib
For all of the following examples, you will have more informative text output when running the code if you first enable Python’s logging mechanism to display log messages to screen. This can sometimes be more verbose than is desired, so feel free to turn this on or off as desired.
import logging
logging.getLogger('poppy').setLevel(logging.WARN)
#Can be logging.CRITICAL, logging.WARN, logging.INFO, logging.DEBUG for increasingly verbose output
OpticalElements can be instantiated from FITS files, or created by one of a large number of analytic function definitions implemented as AnalyticOpticalElement subclasses. Typically these classes take some number of arguments to set their properties. Once instantiated, any analytic function can be displayed on screen, sampled onto a numerical grid, and/or saved to disk.:
ap = poppy.CircularAperture(radius=2)
values = ap.sample(npix=512) # evaluate on 512 x 512 grid
plt.figure(figsize=(10,5))
ap.display(what='both') # display phase and amplitude transmission
Optics come in two basic types: pupil and image plane optics. The circular aperture above is a pupil plane optic, as are SquareAperture, HexagonalAperture, and many other classes. Pupil plane optics are dimensioned in meters.
Image plane optics include various types of coronagraph and spectrograph image plane masks. Image plane optics are dimensioned in arcseconds for convenience. Their physical scale is not modeled in detail, but is scaled proportionate to the implied angular scale of any images formed from the pupil plane optics.
bar = poppy.optics.BarOcculter()
bar.display()
my_slit = poppy.optics.RectangularFieldStop(width=2, height=0.2, name='Horizontal Slit')
my_slit.display()
Optics can have phase as well as amplitude. Depending on the particular class, phases may be specified using physical units for optical path difference in nanometers, or specified in terms of some number of waves for a given wavelength.
Note that wavelengths are always in meters.
poppy.ThinLens(nwaves=2, reference_wavelength=1e-6).display(what='phase')
In addition to image and pupil optics, there are Detector objects which define how to sample a PSF into an array (pixel scale, field of view), and there are Rotation objects, which represent a coordinate system rotation around the optical axis.
One or more OpticalElements can be arranged into an OpticalSystem, which functions in many ways like a python List, but distinguishes between pupil and image planes. OpticalSystems should start with at least one pupil plane optic, and must end with a Detector element that defines the desired sampling.
A basic circular pupil is simple to set up and to compute the PSF for, as it should be.
osys = poppy.OpticalSystem()
osys.addPupil( poppy.CircularAperture(radius=3)) # pupil radius in meters
osys.addDetector(pixelscale=0.010, fov_arcsec=2.0) # image plane coordinates in arcseconds
psf = osys.calcPSF(2e-6) # wavelength in microns
poppy.display_PSF(psf, title='The Airy Function')
Note that the returned PSF is in fact a FITS HDUlist object (from astropy.io.fits
), suitable for saving to disk for use with standard astronomical tools.
The FITS header contains information about the calculation that produced the PSF.
psf
[<astropy.io.fits.hdu.image.PrimaryHDU at 0x110e2cf10>]
psf[0].header
SIMPLE = T / conforms to FITS standard BITPIX = -64 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 400 NAXIS2 = 400 EXTEND = T PLANE1 = 'Wavefront Intensity' WAVELEN = 2E-06 / Weighted mean wavelength in meters DIFFLMT = 0.04583666666666666 / Diffraction limit lambda/D in arcsec OVERSAMP= 2 / Oversampling factor for FFTs in computation DET_SAMP= 2 / Oversampling factor for MFT to detector plane PIXELSCL= 0.005 / Scale in arcsec/pix (after oversampling) FOV = 2.0 / Field of view in arcsec (full array) NWAVES = 1 / Number of wavelengths used in calculation WAVE0 = 2E-06 / Wavelength 0 WGHT0 = 1.0 / Wavelength weight 0 FFTTYPE = 'pyFFTW ' / Algorithm for FFTs: numpy or fftw NORMALIZ= 'first ' / PSF normalization method HISTORY Created wavefront: wavelen=2e-06 m, diam=9.000000 m HISTORY using array size (1024, 1024) HISTORY Multiplied WF by phasor for Pupil plane: Circle, radius=3.00 m (Analyt HISTORY ic) HISTORY Propagating wavefront to Detector plane: Detector (200x200, 0.010000 a HISTORY rcsec/pixel). HISTORY Propagating w/ MFT: 0.0050"/pix fov=43.633 lam/D npix=400 HISTORY Calculation completed in 0.295 seconds
By combining multiple analytic optics together it is possible to create quite complex pupils:
ap = poppy.MultiHexagonAperture(rings=3, flattoflat=2) # 3 rings of 2 m segments yields 14.1 m circumscribed diameter
sec = poppy.SecondaryObscuration(secondary_radius=1.5, n_supports=4, support_width=0.1) # secondary with spiders
atlast = poppy.CompoundAnalyticOptic( opticslist=[ap, sec], name='Mock ATLAST') # combine into one optic
plt.figure(figsize=(12,8))
atlast.display(npix=1024, colorbar_orientation='vertical')
And here’s the PSF:
osys = poppy.OpticalSystem()
osys.addPupil(atlast)
osys.addDetector(pixelscale=0.010, fov_arcsec=2.0)
psf = osys.calcPSF(1e-6)
plt.figure(figsize=(12,8))
poppy.display_PSF(psf, title="Mock ATLAST PSF")
Many optical systems involve propagation between intermediate image and pupil planes. For instance, spectrographs will have an entrance slit at an intermediate image plane, while coronagraphs can have a dizzying variety of image and/or pupil plane stops, occulters, and phase masks to shape the diffraction pattern.
As an example of a more complicated calculation, here’s a NIRCam-style band limited coronagraph. In this example we intentionally create some blank planes in order to show the intensity both before and after the coronagraph stops.
oversample=2
pixelscale = 0.010 #arcsec/pixel
wavelength = 4.6e-6
osys = poppy.OpticalSystem("test", oversample=oversample)
osys.addPupil(poppy.CircularAperture(radius=3.25))
osys.addImage()
osys.addImage(poppy.BandLimitedCoron(kind='circular', sigma=5.0))
osys.addPupil()
osys.addPupil(poppy.CircularAperture(radius=3))
osys.addDetector(pixelscale=pixelscale, fov_arcsec=3.0)
plt.figure(figsize=(14,10))
corpsf = osys.calcPSF(wavelength=wavelength, display_intermediates=True)
print "The total intensity in the output PSF is only {0:.3f} of the input.".format(psf[0].data.sum())
The total intensity in the output PSF is only 0.985 of the input.
The OpticalSystem class has a few attributes that can be used to adjust the properties of the calculation, for instance the location of the source can be offset from the center of the field.
Let's repeat the same calculation, but with the source moved off-center.
osys.source_offset_theta = 45.
osys.source_offset_r = 0.1 # arcsec
plt.figure(figsize=(14,10))
corpsf2 = osys.calcPSF(wavelength=wavelength, display_intermediates=True)
poppy
includes several utility functions for evaluating and displaying PSFs. These all expect PSFs in the HDUlist format returned from poppy
calculations, and use the FITS header metadata when rendering the display.
We've already seen the display_PSF
function.
poppy.display_PSF(psf)
There is also a function for showing the difference between two PSFs.
poppy.display_PSF_difference(corpsf,corpsf2)
And one for showing radial profiles:
poppy.display_profiles(psf)
One can also make a variety of measurements, for instance the center.
In this case, we see that for an even-sized array, the PSF is by default centered exactly in the middle between 4 pixels symmetrically.
print "The PSF array shape is ", psf[0].data.shape
print "The centroid is located at ({0:.3f}, {1:.3f})".format(*poppy.measure_centroid(psf))
The PSF array shape is (400, 400) The centroid is located at (199.500, 199.500)
poppy.display_PSF(psf, imagecrop=0.1)
print "The FWHM is {0:.3f} arcsec at {1} microns ".format(poppy.measure_fwhm(psf), psf[0].header['WAVELEN']*1e6)
The FWHM is 0.017 arcsec at 1.0 microns
The package documentation lists the other available OpticalElement classes, and has some examples of more complex coronagraphic calculations.
You may find that you need an OpticalElement which is not represented among the various predefined classes. There are two ways to proceed in this case.
We encourage users to contribute their own models and simulations back into the poppy Examples page. Similarly if you create any generally useful OpticalElement classes, pull requests are very much welcomed.