mperrin / poppy

Physical Optics Propagation in Python
BSD 3-Clause "New" or "Revised" License
175 stars 41 forks source link

Unexpected result for short-distance Fresnel propagation #194

Closed josePhoenix closed 7 years ago

josePhoenix commented 8 years ago

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:

after_poppy

The user's other propagation software produces this pattern:

maciekgroch_pattern

@douglase, any ideas?

josePhoenix commented 8 years ago

@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 commented 8 years ago

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 commented 8 years ago

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 commented 8 years ago

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 commented 8 years ago

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 commented 8 years ago

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 commented 8 years ago

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

douglase commented 8 years ago

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

propagate_direct

mperrin commented 8 years ago

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 commented 8 years ago

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 commented 8 years ago

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

josePhoenix commented 8 years ago

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 commented 8 years ago

@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 commented 8 years ago

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 commented 8 years ago

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

maciekgroch commented 8 years ago

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 commented 8 years ago

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 commented 8 years ago

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

douglase commented 8 years ago

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 commented 8 years ago

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

maciekgroch commented 8 years ago

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

maciekgroch commented 8 years ago

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)

nice

mperrin commented 8 years ago

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

mperrin commented 8 years ago

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:

screen shot 2016-10-18 at 9 17 54 am

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 commented 8 years ago

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

douglase commented 8 years ago

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 commented 8 years ago

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 commented 8 years ago

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

douglase commented 8 years ago

@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 commented 8 years ago

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

maciekgroch commented 8 years ago

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

mperrin commented 7 years ago

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