spacetelescope / poppy

Physical Optics Propagation in Python
https://poppy-optics.readthedocs.io
BSD 3-Clause "New" or "Revised" License
220 stars 72 forks source link

Unexpected modification of the FresnelWavefront matrix size after applying rotation #572

Closed SteveOstile closed 1 year ago

SteveOstile commented 1 year ago

Dear Poppy team,

Thanks for nice work on Poppy !

I've found an unexpected modification of the FresnelWavefront matrix size after applying rotation using the following code :


import matplotlib.pyplot as plt
import numpy as np
import poppy
import astropy.units as u

w0 = 0.1 # in meter
lamb = 532e-9 # in meter

wf0 = poppy.FresnelWavefront(beam_radius=w0*u.m, wavelength=lamb)
ap = poppy.CircularAperture(radius=w0)
wf0 *= ap

#generate tilt to better see the following rotation
wf0.tilt(Xangle=10.0)

#Display before rotation
plt.figure(figsize=(12,5))
wf0.display('both',colorbar=True, showpadding=True)
plt.suptitle("Before Rotation", fontsize=18)

#Copy wavefront and apply 90° rotation
wf1 = wf0
print('Before rotation:',np.shape(wf1))
wf1.rotate(90.0)
print('After rotation:',np.shape(wf1))

#Display after rotation
plt.figure(figsize=(12,5))
wf1.display('both',colorbar=True, showpadding=True)
plt.suptitle("After rotation", fontsize=18)

The 2 print outputs are:

Before rotation: (2048, 2048) After rotation: (1024, 1024)

Displaying amplitude and phase before rotation: BeforeRotation

Displaying amplitude and phase after rotation: AfterRotation

As you can see, the matrix size goes from 2048 x 2048 to 1024x1024 and the width of the beam is doubled. Also this generates the following error while propagating the rotated wavefront #

Traceback (most recent call last):

  File "D:\Steve\ProgPython\Poppy\untitled1.py", line 53, in <module>
    wf1.propagate_fresnel(150)

  File "C:\ProgramData\Anaconda3\envs\spyder\lib\site-packages\poppy\utils.py", line 1436, in unit_check_wrapper
    return wrapped_function(*bound_args.args, **bound_args.kwargs)

  File "C:\ProgramData\Anaconda3\envs\spyder\lib\site-packages\poppy\fresnel.py", line 864, in propagate_fresnel
    self._propagate_ptp(delta_z)

  File "C:\ProgramData\Anaconda3\envs\spyder\lib\site-packages\poppy\utils.py", line 1436, in unit_check_wrapper
    return wrapped_function(*bound_args.args, **bound_args.kwargs)

  File "C:\ProgramData\Anaconda3\envs\spyder\lib\site-packages\poppy\fresnel.py", line 730, in _propagate_ptp
    self.wavefront *= exp_t  # eq. 6.68

ValueError: operands could not be broadcast together with shapes (1024,1024) (2048,2048) (1024,1024) 

Thanks for looking at this issue.

mperrin commented 1 year ago

Thanks for the report!

I agree this behavior is not obvious. That said it's at least partially intentional, though I agree there's definitely a bug still. Looking back at the code of the rotate function, it appears the logic was/is to avoid rotating large arrays of zero-padding, for performance reasons. The zero padding is supposed to be added back in later, if needed. See this part of code, from the definition of the rotate function:

        if self.ispadded:
            # pupil plane is padded - trim out the zeros since it's not needed in the rotation
            # If needed in later steps, the padding will be re-added automatically
            self.wavefront = utils.remove_padding(self.wavefront, self.oversample)
            self.ispadded = False

This all works as intended for the Fraunhofer code. But definitely something is not as intended for the Fresnel side of things. First, the apparent spatial extent of the illuminated beam certainly shouldn't change. Secondly, it's clearly not adding back in the zero padding before trying a subsequent Fresnel propagation.

The first problem appears to be related to extra internal index arrays the Fresnel code maintains, which don't get updated for the change in padding. The not adding back in the zero padding is simple: the propagate_fresnel code simply doesn't check for an unpadded array.

I think the most straightforward fix is likely to subclass the rotation function for FresnelWavefront and leave out the bit that removes the padding. Seems to me that any speed gains from not rotating the many zeros are not worth the added complexity in this case.

mperrin commented 1 year ago

Oh and @SteveOstile: welcome! Thanks for reporting this. And I see this looks like your first comment on a GitHub issue so let me extend an extra welcome in that case :-)

mperrin commented 1 year ago

Fix added in PR #573. @SteveOstile, let me know if you're comfortable enough with git/gitub to perhaps check out and test that code? For me, running with the modified version in the PR, your example code above runs with no surprises.

SteveOstile commented 1 year ago

Hi @mperrin,

Thanks for reacting so quickly to my comment! I understand the issue now. Well, I am so new on github that I am not very comfortable in reviewing your modification. I have to look how I can pip install your modification...

mperrin commented 1 year ago

Hi @SteveOstile we've merged the fix into the develop branch here. I believe you should be able to install from GitHub with something like this:

pip install git+https://github.com/spacetelescope/poppy.git@develop#egg=poppy
SteveOstile commented 1 year ago

@mperrin Excellent ! It works perfectly !