making pixelscale of direct have length/pix units

[mperrin]

Issue by douglase Wednesday Oct 05, 2016 at 19:45 GMT Originally opened as https://github.com/mperrin/poppy/pull/196

---

---

douglase included the following code: https://github.com/mperrin/poppy/pull/196/commits

[mperrin]

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

---

This PR also forces several .value calls for values with units to first convert to the base unit (meters in these cases) before taking the value.

[mperrin]

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

---

Hm, trying to make sure I understand the pitfalls of peeling the units off Astropy Quantity objects... should I go grepping through the source to find any .value that's not preceded by a .to(u.Thing)?

[mperrin]

Comment by josePhoenix Wednesday Oct 05, 2016 at 19:50 GMT

---

(Well, any user-supplied quantity...)

[mperrin]

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

---

try:

radius=6*u.imperial.furlong
radius.value

I just went through fresnel and didn't see any other cases except for in coordinates(), but there it is forced to specific units in a preceding line.

[mperrin]

Comment by josePhoenix Wednesday Oct 05, 2016 at 19:53 GMT

---

Okay, yep, that's what I thought would happen. Good to know about that particular footgun.

[mperrin]

Comment by coveralls Wednesday Oct 05, 2016 at 19:56 GMT

---

Coverage remained the same at 64.201% when pulling 95ef404e1e9394859688e734fbd2fdb129f06e7d on douglase:fix_direct_pixscl into 0715741a5b161d8f243ea8068777a1f2d633a03c on mperrin:master.

[mperrin]

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

---

alternatively, some of the .value calls could probably keep their units. I don't recall why the .value is being taken for the list of waists, but at least this way the behavior is consistent and won't break anything.

[mperrin]

Comment by mperrin Wednesday Oct 05, 2016 at 20:03 GMT

---

Yeah, for the list of waists is that actually used for anything? I thought that was mostly just kept as a sort of log of the optical propagation? It seems like that might be best to keep as a Quantity rather than trimming to just a value. I don't see the benefit in doing that, whereas keeping the units explicitly there seems like it could be useful.

[mperrin]

Comment by douglase Wednesday Oct 05, 2016 at 20:10 GMT

---

Yeah, it was for checking against Gaussian laser waist propagation codes, units could be helpful.

On Wed, Oct 5, 2016 at 4:03 PM, Marshall Perrin notifications@github.com wrote:

Yeah, for the list of waists is that actually used for anything? I thought that was mostly just kept as a sort of log of the optical propagation? It seems like that might be best to keep as a Quantity rather than trimming to just a value. I don't see the benefit in doing that, whereas keeping the units explicitly there seems like it could be useful.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/mperrin/poppy/pull/196#issuecomment-251783121, or mute the thread https://github.com/notifications/unsubscribe-auth/AA-nn_v6Cf3A4sV44IC6J6yYMp6pOjSMks5qxAKMgaJpZM4KPNYU .

Postdoctoral Research Associate, MIT AeroAstro.

[mperrin]

Comment by mperrin Thursday Oct 06, 2016 at 01:40 GMT

---

@douglase this revised propagate_direct doesn't seem to be working for me. I tried testing it, but the units don't come out properly in the end. I get a complex wavefront that's a Quantity with units of m/pix^2, and that's no good. Breaks all sorts of stuff since the code expects that to be dimensionless. Is this version of the propagate_direct working for you?

Here's a small test case which fails for me:

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

I think I've figured out the necessary fixes which I'll comment in the code review.

[mperrin]

Comment by douglase Sunday Oct 09, 2016 at 05:00 GMT

---

So, I found the problem with the short distance propagation, a scaling constant in the ptp function was wrong. Unfortunate that it gave such good agreement at long distances, shows that test line-by-line test coverage is not the same as testing physical extremes.

This latest update addresses #194, as demonstrated by an additional test within test_Circular_Aperture_PTP(). It also addresses the comments above about display() not working.

The Hubble propagation test by @mperrin is now failing, the Airy function disagrees with the Fresnel diffraction now by closer to 5%.

But I'm not sure this is a problem, the HST Fresnel diffraction pattern agrees well with the Fraunhofer pattern for the same aperture (the source code for this figure in the demo notebook in this branch, https://github.com/douglase/poppy/blob/e8c2a90aeec81fbe2f3b205850c64860e5c6bcfd/notebooks/Fresnel_Propagation_Demo.ipynb):

Is there some reason to believe the diffractive HST core is actually as close to an Airy pattern as the test found previously, or was that just a coincidence?

[mperrin]

Comment by mperrin Sunday Oct 09, 2016 at 14:03 GMT

---

Actually it looks like when I wrote that test originally, I forgot to set the obscuration parameter for the size of the secondary mirror in the call to airy_2d. I expect if you were to fix that the agreement would be much better. The presence of the obscuration redistributes more and less light in the different rings of the Airy function and I'm sure that's the reason for more intensity in the first Airy ring.

Glad to hear you found the bug in the scaling factor! Yeah, I agree it's a bit surprising the code was working as well as it did before with that wrong. You make a very good point about testing at extreme cases motivated by the physics.

[mperrin]

Comment by mperrin Monday Oct 10, 2016 at 17:01 GMT

---

I took a closer look at this again today and my above comment isn't the answer. There's a difference between the notebook and the unit test versions of the HST model. The version in the notebook does have the SM obscuration and SM supports. The version in the unit test has an unobscured aperture (despite having an SM in the optical path - so the simulation is physically inconsistent: the SM doesn't block any light on the way in but does reflect it after it bounces off the primary. Nonphysical but thinking back I'm sure I wrote the test intentionally that way so the comparison with the unobscured Airy function was more exact.

To answer your above question "is there some reason to believe the diffractive HST core is actually as close to an Airy pattern as the test found previously?" in this case yes. Unobscured aperture at the focus so the far-field Fresnel should asymptotically approach the Fraunhofer result, and the analytic Airy function as well.

I wonder if the issues is not fine enough pixel sampling? Hmm, I just tried a calculation with beam_ratio=0.125 so in effect twice as fine sampling in the Fresnel calculation. But the test still doesn't pass. So I'm still unsure what's going on.

beam ratio = 0.25; 4x oversampling:

beam ratio = 0.125; 8x oversamplingL

[mperrin]

Comment by douglase Monday Oct 10, 2016 at 17:38 GMT

---

could aliasing from the apertures be filling in the dark rings?

[mperrin]

Comment by mperrin Wednesday Oct 12, 2016 at 16:52 GMT

---

could aliasing from the apertures be filling in the dark rings?

Maybe? I just worked a bit on the long-postponed subpixel_geometry branch, to bring it up to date with master. But some tests are still failing on that, including the Fresnel PTP one now. I'm out of time to work on this today, need to switch to other tasks, but if you have time you can experiment with what's on that branch.

[mperrin]

Comment by mperrin Thursday Oct 13, 2016 at 03:14 GMT

---

Ha, would you believe it's just pure defocus error? The value for d_sec_to_focus leaves the final focal plane off by 330 microns from the beam waist.

With d_sec_to_focus = 6.3916645 * u.m instead, so the focus is precisely at the beam waist, the discrepancy vanishes.

Here's the result of a calculation done with that value and beam ratio=0.125 (8x oversampling). Great agreement. Note, this is for a case where I removed the SM obscuration and struts from the notebook calc for simplicity.

I've just checked and with the d_sec_to_focus value changed to the above, the test passes just fine. No other changes needed. Storing the list of beam waists has just paid off! That was essential for debugging this.

One other issue, but a way more minor effect, is the difference in PSF alignment between the Fresnel and Fraunhofer codes. (PSF center is at e.g. 2047.5 pixels in Fraunhofer, but at 2048 pixels in Fresnel). The radial_profile function isn't smart enough to figure that out on its own so you have to pass it in as a parameter. But that's a far less significant effect than the defocus. It matters a tiny bit for making the radial profile plots look nice, but doesn't affect the unit test since that comparison is done on the array pixel values directly instead of making profiles.

[mperrin]

Comment by douglase Thursday Oct 13, 2016 at 03:50 GMT

---

nice!

[mperrin]

Comment by mperrin Thursday Oct 13, 2016 at 14:52 GMT

---

I'm going to go ahead and merge this then fix the propagation distance.