Unexpected result for short-distance Fresnel propagation

[josePhoenix]

User @maciekgroch reports unusual behavior when propagating a 10nm wave a short distance past a circular aperture using the following code:

import poppy
import astropy.units as u

npix = 1024
wf = poppy.FresnelWavefront(
    5 * u.um,
    wavelength=10 * u.nm,
    npix=npix,
    oversample=4
)
wf *= poppy.CircularAperture(radius=800 * u.nm)
plt.figure()
wf.display()

z = 12. * u.um
plt.figure()
wf.propagate_fresnel(z, display_intermed=True)

POPPY produces this intensity pattern after propagating 12 um:

The user's other propagation software produces this pattern:

@douglase, any ideas?

[josePhoenix]

@maciekgroch Are you sure the lengths you're using are the same in both tools? What tool produced the second image? I'm not familiar with Fresnel propagation at very short wavelengths.

[maciekgroch]

I am using a custom-written tool we developed at the university (it's not public unfortunately). I thought it should work for every wavelength.

[mperrin]

Have you looked at the debug log from this propagation? If you have webbpsf installed the easiest way is you can use webbpsf.setup_logging('debug') to turn on maximum verbosity. If you don't have webbpsf then you can manually use the Python logging API to add a log handler to save the log to a file or display on screen.

[josePhoenix]

I attempted to run with logging, but didn't see any output from POPPY when running this code. I'm not sure if that's a separate issue, or if there are no logging calls in the code used for this short example...

[josePhoenix]

Ah, it's level='DEBUG'. Here we go:

[  poppy] Padded WF array for oversampling by 4, to (4096, 4096).
[  poppy]   Multiplied WF by phasor for Optic: Circle, radius=800.0 nm
[  poppy] Beginning Fresnel Prop. Waist at z = 0.0 m
[  poppy]   Plane to Plane Regime, dz=12.0 um
[  poppy]   Constant Pixelscale: 9.765624999999999e-09 m / pix
[  poppy]    Using numpy FFT
[  poppy]    Using numpy FFT
[  poppy] ------ Propagated to plane of type PlaneType.intermediate at z = 1.20e-05 m ------

[mperrin]

That's weird, propagate_fresnel is full of debug logging calls.

(edit: this message posted before I got josePhoenix's second post immediately above)

[josePhoenix]

Here's the test case as a notebook: https://gist.github.com/josePhoenix/edfe9082b0230fb063fda5202d07d6b1

[douglase]

Interesting. propagate_direct gives an answer that appears to agree
with @maciekgroch's figure (see below). I will look more at ptp and see if something is the matter, I don't think it is exercised well by any of the test cases. the intent was always to write a test case to check them against each other.

import poppy
import astropy.units as u
import logging
logging.basicConfig(format='%(levelname)s:%(message)s', level=logging.DEBUG)
npix = 1024
wf = poppy.FresnelWavefront(
    5 * u.um,
    wavelength=10 * u.nm,
    npix=npix,
    oversample=4
)
wf *= poppy.CircularAperture(radius=800 * u.nm)

z = 12. * u.um

wf.propagate_direct(z)

plt.imshow(wf.intensity.value,cmap="gray")
plt.colorbar()
plt.savefig("propagate_direct.png")

[mperrin]

hmm, are we sure that propagate_fresnel is correct to invoke the propagate_ptp function for this case? This may be related to the questions from Joyce Fang in email, regarding how to set up the beam for a Fresnel propagation. In this case the constructor sets up a Fresnel wavefront with a beam waist of 5 microns. Then propagating it through the CircularAperture changes the intensity pattern but does not update the beam waist value. As a result the Rayleigh distance is calculated for 5 microns waist even though the beam is only 800 nm in radius. Would the _propagate_wts be a better choice in this case for the default code?

[maciekgroch]

I tried the one with propagate_direct. Results look good, however, I found another issue:

In: wf.pixelscale
Out: <Quantity 5000000.0 nm pix um / m>

I guess there is something wrong with units conversion🤕

[douglase]

@maciekgroch, good catch on the pixel scale, just put in pull request #195 which should fix the units.

[josePhoenix]

Thank you @douglase for the quick response! Also thank you @maciekgroch for exercising this corner of the code and taking the time to tell us what you found 😄

[douglase]

@mperrin, for a 800 nm waist at 10nm the Rayleigh length is 200 microns so this is well within the plane to plane regime. but whether the plane to plane code captures the physics is another question.

[mperrin]

yeah thanks for double checking that about the Rayleigh length. I did some more calculations after my post above and came to a similar conclusion, so that's probably not it.

But doesn't our test_Circular_Aperture_PTP test case cover the plane to plane code? That's a test against values from a plot in the Anderson & Enmark reference.

[douglase]

yeah, I forgot we finished test_Circular_Aperture_PTP. the mystery deepens.

[maciekgroch]

Are you sure that you fixed the units issue?
I pulled master and I have: <Quantity 3264000.0000000005 nm um / m>
Which is not what I would expect :). I looked into the code:
s = self.n * u.pix * self.pixelscale
If you multiply here by u.pix then:
self.pixelscale = self.pixelscale.to(u.meter) / u.pix
IMHO this is the correct way (conversion to meter makes it clean -- instead of nm*um/m, just meters are created).

What do you think?

[douglase]

I have been trying to avoid extra .to() calls unless there is a specific reason. they seem like clutter but that is a matter of taste. I did notice today that there are a few .value calls without to(), which could lead to unexpected behavior.

astropy simplifies any combination of units to SI units if you call wavefront.pixelscale.decompose()

[maciekgroch]

Alright then, I have never worked with Astropy before so I didn't know about decompose method. What about division by u.pix?

[douglase]

I don't want to redefine s because it now has the correct units for the reference cited but that's a good point, everywhere else plate scales are length/pixels, so I divided by u.pix on line 537 in a new PR, #196.

[douglase]

@maciekgroch, can you check the latest version of master addresses the problems you saw?

[maciekgroch]

@douglase, Sorry for the late replay, I was away again. I will work on this today.

[maciekgroch]

Yeah it works perfect! Thank you for help. I think the issue may be closed.
Executed the following code:

import poppy
import astropy.units as u

npix = 1024
wf = poppy.FresnelWavefront(
    5 * u.um,
    wavelength=10 * u.nm,
    npix=npix,
    oversample=4
)
wf *= poppy.CircularAperture(radius=800 * u.nm)
z = 12. * u.um
plt.figure()
wf.propagate_fresnel(z)
wf.display(imagecrop=0.000005)

[mperrin]

@maciekgroch Thanks for your help in identifying and fixing this problem! Much appreciated. :-)

[mperrin]

Wait, upon closer examination I'm not 100% convinced everything is in agreement. It looks like the diffraction patterns at the start and end of this thread still don't quite match:

Note the differences in intensity of the various rings. Looking at the code it seems these are both for the same test case wavelength and propagation distance, etc? So that seems to show there is still some disagreement. Do we need to continue tracking this down more?

[maciekgroch]

Well I executed exactly the same params lambda = 10 nm, aperture of radius = 0.8 um, propagation on d = 12 um.
The results as following:
(To be honest I used a bit different wavelength (10.88 nm) in previous examples. That caused all the problems. Sorry for that.)

[douglase]

I don't have quick access a commercial package at the moment, but I ran PROPER in GDL with the prescription below and got a result that looks consistent with the current POPPY output. @maciekgroch, have you tested your pattern against other models or do you know of any 3rd party published figures we can try to recreate in the two packages?

PROPER:

gridsize=1024
sampling=4
wavelength=10e-9
print_it=1
PROP_VERBOSE= 1
use_fftw=0
prop_init_savestate
d_objective = 1600e-9           ;-- objective diameter in meters
beam_ratio = 0.4
prop_begin, wavefront, d_objective, wavelength, gridsize, beam_ratio
d_intermediate_image = 12e-6
prop_circular_aperture, wavefront, d_objective/2.0
prop_define_entrance, wavefront
prop_propagate, wavefront, d_intermediate_image, 'intermediate image'
prop_end, wavefront, sampling
writefits, "aperture_test_output.fits", wavefront```

[mperrin]

Aha! Great, @maciekgroch thanks for checking that. I had hoped it would be some small difference in calculation parameters like that which would explain the difference. Very glad to hear that my hunch was right! :-)

@douglase see above - @maciekgroch edited his comment after figuring out the solution. :-)

@douglase now that you've got a PROPER version of this, maybe we could read off a couple values like the peak intensity and add them as additional assertion checks in the new unit test for this? Would make it that much more robust of a test function.

[maciekgroch]

@douglase, I also experimented with Holopy. Can try it later with our setting.

[douglase]

@maciekgroch that's a relief. Thanks for reporting this and being so responsive with testing.

@mperrin what about including a FITS file for the test?

[josePhoenix]

Just checking, the functionality here has been sorted out and we're only waiting on tests to close the issue?

[maciekgroch]

I believe it is ready -- I experimented a bit and encountered no problems.

[mperrin]

This issue has been fixed since last October, we just never actually closed the Github issue then.