Speckle Simulation

Speckle simulation involves constructing a wrinkled wavefront and an aperture illumination function, then performing a Fourier transform using this amplitude/phase information to generate the far-field diffraction pattern arising from this input.

The most complicated part is generating a realistic wavefront. To do this, we use a Zernike basis, and a "power spectrum" of strength in each term. A Gaussian random variable for each Zernike term scales its influence in each instance. Each term is summed to produce a total wavefront. The series is terminated once a target precision is reached, where the amplitudes of the new additions become smaller than some threshold (default = 0.05 radians).

The coding was originally in IDL, recently ported to Python, with slight enhancements. First, I'll describe the common traits.

Common Parameters/Traits

One must first decide on a D/r₀ ratio to simulate. D is the aperture diameter, and r₀ is the characteristic scale (the Fried parameter) of wavefront distortions, so that the long-exposure "seeing" measure is λ/r₀ rather than λ/D.

The level to which the simulation is carried out will impact the fidelity (and speed) of the simulation. The default value is 0.05, meaning that higher-order terms contributing less than 0.05 radians to the sum will be truncated. This does not mean that the wavefront is solid to 0.05 radians, as smaller terms can add to something larger than this. But it probably does mean fidelity to about 3–4 times this level.

In principle, each simulation can be adjusted to the number of pixels covered (slightly easier in the Pyhton version); the default is 256.

The Zernike amplitude spectrum followed is as prescribed in Hardy (page 96) in both cases: the first 10 terms are from Table 3.3, and the remaining terms follow the power law set by Eq. 3.71.

Padding: In order to resolve the structure with adequate sampling, it is necessary to embed the aperture/phase function into a larger (blank) frame, perform the Fourier transform, then clip out the central region you care about. In the Python code, this is handled via the variable psf_scale. In IDL, it is somewhat more manually set (see code snippet below). It is easy to get confused about the pixel scale through these manipulations. If the scaling factor is α, then the final image spans an angle (in radians) of N_pixλ/αD. So for example, dealing with the 3.5 meter scope, 256 pixels, 532 nm, and a padding factor α=3 yields a frame 13 microradioans, or 2.67 arcseconds across.

Python Implementation

Two codes, wavefront.pyprog and speckle.pyprog (change names to *.py once on your machine) do all the work. The wavefront.py code contains utilities to: generate the Zernike basis functions; an aperture-generating bit that contains central obscuration and spider parameters; a plane-wave generator; a "seeing" wavefront generator (essentially the same as the IDL wavefront.pro described below); a generator of rhe far-field pattern; and an approximate calculator of the flux within an aperture. The speckle.py code is what you actualy run. It simply generates an example speckle pattern, given various parameters set within the program. It generates some gif images of the aperture, wavefront, the speckle pattern, and what the aperture would deliver for a perfect planar wavefront. It also makes fits images of the speckle and perfect diffraction pattern, as well as printing some diagnostic information.

IDL Implementation

Two codes, zernike.pro and wavefront.pro get you most of the way there. The wavefront code calls the zernike code. In order to turn the wavefront into a PSF, use code like the following (for hard-coded α=6):

npix = 256 x = findgen(npix)-(npix/2 - 1) y = findgen(npix)-(npix/2 - 1) xarr = x # replicate(1.0,npix) yarr = replicate(1.0,npix) # y rarr = sqrt(xarr^2 + yarr^2)/(npix/2) outside = where(rarr gt 1.0) ampl = fltarr(npix,npix)*0.0 + 1.0 ampl[outside]=0.0 input = complexarr(1536,1536) wf = wavefront(25.0,level=0.05) input[640:895,640:895] = complex(ampl*cos(wf),ampl*sin(wf)) im = abs((shift(fft(input),768,768))[640:895,640:895])^2 ;im = rebin(im[64:191,64:191],256,256)/4.0 ;; optional rebin

The first block makes a circular aperture (unobstructed). The second block defines a larger complex array (this is where the α=6 implicitly enters) and pastes the center with the aperture-limited complex phase, then takes the FFT and shifts and chops accordingly. The im output may then be written as a fits file, if desired.

Assignment

You are to use the code from either platform to understand something about wavefronts, seeing, etc. Whatever you do, you will need to describe your parameters like D, λ, r₀, etc. (and context, in some cases). Report the peak-to-peak wavefront error, in microns (will need to assume some wavelength).

Pick some topic to explore that involves configuring and running the code to learn something of interest to you. Possible topics might look like:

Pick a telescope/aperture and a seeing value; come up with the corresponding D/r₀ and evaluate the speckle pattern: you might look into:
- approximate number of speckles
- approximate size of blur
- characteristic offset of pattern from center
- long-term blur (adding many instances together) (see summing code)
- any impact of central obscuration/spider (more evident in diffraction-limited image or low D/r₀
Explore the appearance of speckle patterns as a function of wavelength if seeing scales like λ^−1/5
Truncate the series to only produce the tip/tilt Zernike terms so you can characterize how much of the seeing is caused by tip-tilt for a particular telescope/seeing.
other