unexpected ZernikeWFE.get_opd() behavior

[douglase]

Kudos to @tukyab for helping find this one.

setup of problem:

r=.0254
osys = poppy.OpticalSystem()
osys.add_pupil(poppy.CircularAperture(radius=r))    
osys.add_pupil(poppy.ZernikeWFE(radius=r,coefficients=[0,0,0.000001,0]))   

Expected behavior:

plt.imshow(osys.planes[1].get_opd(osys.input_wavefront(wavelength=600e-9)))
plt.colorbar()

gives:

unexpected behavior:

plt.imshow(osys.planes[1].get_opd(600e-9))
plt.colorbar()

This works when the aperture is a few meters across, but only there, it doesn't rescale at all (if you make r=8 then you can't see the circle anymore), the pixel_scale seems to be getting lost when you just pass a wavelength to get_opd().

The documentation says you can pass a wavelength as a float or a wavefront, but I don't actually see why it doesn't throw an error, since the self.circular_aperture.get_transmission(wave) function tests that wave is a wavefront. @josePhoenix ?

[josePhoenix]

Well, that's weird. I'll take a look.

[josePhoenix]

Something having to do with the physical dimensions of the input wavefront is getting set in osys.input_wavefront but not when creating a wavefront all on its own. Further investigation is required.

[josePhoenix]

By the way: thanks for the high-quality bug report, @douglase!

[josePhoenix]

The relevant section of OpticalSystem.input_wavefront is:

        # somewhat complicated logic here for historical reasons.
        # if we have a first optical plane, check and see if it specifies the entrance sampling.

        npix=None
        diam=None
        if len(self.planes) > 0:
            if self.planes[0].shape is not None:
                npix=self.planes[0].shape[0]
            if hasattr(self.planes[0], 'pupil_diam') and self.planes[0].pupil_diam is not None:
                diam = self.planes[0].pupil_diam

        # if sampling is still undefined, fall back to what is set for this optical system itself
        if npix is None:
            npix=self.npix if self.npix is not None else 1024
        if diam is None:
            diam = self.pupil_diameter if self.pupil_diameter is not None else 1

        # if the diameter was specified as an astropy.Quantity, cast it to just a scalar in meters
        if isinstance(diam, u.Quantity):
            diam = diam.to(u.m).value

What's happening is that in your first case, you're computing the Zernike disk on a 1024 px square array that corresponds to a diameter of 2 * 0.0254. Because the Zernike disk has the same diameter, it's well-sampled and looks like your first plot. From poppy.optics.CircularAperture.__init__:

        self.pupil_diam = pad_factor * 2 * self.radius

(Pad factor is 1.0 unless otherwise specified.)

In the second case, you're computing the Zernike disk on a Wavefront that's created "just in time" by a wrapper called _accept_wavefront_or_meters in wfe.py. It doesn't know the diameter of the pupil at initialization, so it defaults to 8 meters. This was probably a good choice for API ergonomics when POPPY was primarily used to simulate space telescopes with pupils several meters across 😄 Now, we need to decide how to proceed.

[josePhoenix]

Here's a workaround:

import poppy
zernike_optic = poppy.ZernikeWFE(radius=r,
                                 coefficients=[0, 0, 0.000001, 0])
r = 0.0254
osys = poppy.OpticalSystem()
osys.add_pupil(poppy.CircularAperture(radius=r))
osys.add_pupil(zernike_optic)
osys.display(what='both')
plt.imshow(zernike_optic.get_opd(poppy.Wavefront(600e-9, diam=2 * r)))

[douglase]

Thanks for clearing that up, I didn't see that _accept_wavefront_or_meters wrapper was creating a wavefront. It seems like that function should at least check for pixelscale and diameter before generating a wavefront? (or radius, or pupil_diam, we might want to make a list of all the ways we define the physical extent of a optic first).

[josePhoenix]

Agreed, this is not the least astonishing API. The decorator there should be used consistently as well; a lot of optics have some variant of that test copied and pasted in their various methods. It's on my to-do list, though far from the top...