Unexpected result for short-distance Fresnel propagation

[mperrin]

Issue by josePhoenix
Tuesday Oct 04, 2016 at 20:55 GMT
Originally opened as mperrin/poppy#194

---

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?

[mperrin]

Comment by josePhoenix
Tuesday Oct 04, 2016 at 21:07 GMT

---

@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.

[mperrin]

Comment by maciekgroch
Tuesday Oct 04, 2016 at 21:09 GMT

---

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]

Comment by mperrin
Tuesday Oct 04, 2016 at 21:16 GMT

---

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.

[mperrin]

Comment by josePhoenix
Tuesday Oct 04, 2016 at 21:17 GMT

---

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...

[mperrin]

Comment by josePhoenix
Tuesday Oct 04, 2016 at 21:18 GMT

---

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]

Comment by mperrin
Tuesday Oct 04, 2016 at 21:19 GMT

---

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

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

[mperrin]

Comment by josePhoenix
Tuesday Oct 04, 2016 at 21:24 GMT

---

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

[mperrin]

Comment by douglase
Tuesday Oct 04, 2016 at 21:35 GMT

---

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]

Comment by mperrin
Tuesday Oct 04, 2016 at 21:44 GMT

---

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?

[mperrin]

Comment by maciekgroch
Tuesday Oct 04, 2016 at 21:47 GMT

---

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🤕

[mperrin]

Comment by douglase
Tuesday Oct 04, 2016 at 22:28 GMT

---

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

[mperrin]

Comment by josePhoenix
Tuesday Oct 04, 2016 at 22:30 GMT

---

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 😄

[mperrin]

Comment by douglase
Tuesday Oct 04, 2016 at 22:57 GMT

---

@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]

Comment by mperrin
Tuesday Oct 04, 2016 at 23:10 GMT

---

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.

[mperrin]

Comment by douglase
Wednesday Oct 05, 2016 at 01:01 GMT

---

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

[mperrin]

Comment by maciekgroch
Wednesday Oct 05, 2016 at 18:40 GMT

---

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?

[mperrin]

Comment by douglase
Wednesday Oct 05, 2016 at 18:48 GMT

---

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()

[mperrin]

Comment by maciekgroch
Wednesday Oct 05, 2016 at 18:51 GMT

---

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

[mperrin]

Comment by douglase
Wednesday Oct 05, 2016 at 19:49 GMT

---

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.

[mperrin]

Comment by douglase
Thursday Oct 13, 2016 at 15:15 GMT

---

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

[mperrin]

Comment by maciekgroch
Monday Oct 17, 2016 at 11:31 GMT

---

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

[mperrin]

Comment by maciekgroch
Tuesday Oct 18, 2016 at 08:23 GMT

---

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]

Comment by mperrin
Tuesday Oct 18, 2016 at 13:12 GMT

---

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

[mperrin]

Comment by mperrin
Tuesday Oct 18, 2016 at 13:21 GMT

---

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?

[mperrin]

Comment by maciekgroch
Tuesday Oct 18, 2016 at 18:34 GMT

---

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.)

[mperrin]

Comment by douglase
Tuesday Oct 18, 2016 at 19:07 GMT

---

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]

Comment by mperrin
Tuesday Oct 18, 2016 at 19:09 GMT

---

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.

[mperrin]

Comment by maciekgroch
Tuesday Oct 18, 2016 at 19:10 GMT

---

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

[mperrin]

Comment by douglase
Tuesday Oct 18, 2016 at 19:16 GMT

---

@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?

[mperrin]

Comment by josePhoenix
Friday Oct 28, 2016 at 16:30 GMT

---

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