opticspy / lightpipes

LightPipes for Python, "Pure Python version"
https://opticspy.github.io/lightpipes/
BSD 3-Clause "New" or "Revised" License
228 stars 53 forks source link

switching the order of RectAperture and lens produces different results? #34

Open yellow-force opened 4 years ago

yellow-force commented 4 years ago

Hi all,

I have a simple system with only a lens and a rectangular aperture. I notice that when switching the order of them, I get different phase results. For example:

I suppose for the two scenarios the results should be same. But it seems that when the lens is behind the aperture the phase is not cut-off?

from LightPipes import * import matplotlib.pyplot as plt

size=40mm wavelength=1um N=256 w=20mm f=8m z=4*m

### lens first then aperture Field = Begin(size, wavelength, N) Field = Lens(f,0,0,Field) Field = RectAperture(w, w, 0, 0, 0, Field)

and

### aperture first then lens Field = Begin(size, wavelength, N) Field = RectAperture(w, w, 0, 0, 0, Field) Field = Lens(f,0,0,Field)

2 1

ldoyle commented 4 years ago

Thanks, this is an interesting observation, I gave it a try and found the following: The code evidently distinguishes between 0 and -0, resulting in a different complex phase angle. Even though the code is Cpp based, you can recreate the same behaviour with numpy:

import numpy as np
a = np.zeros((5,),dtype=complex)
a *= np.exp(1j*np.pi/2*np.array([0,1,2,3,4])) #pi/2 * N = phase angle rotation in 90° steps
phi = np.angle(a)
print(f'a = {a}') #f-strings work in python >= 3.6, else use '{}'.format(a)
print(f'phi = {phi}')
#output:
# a = [ 0.+0.j  0.+0.j -0.+0.j  0.-0.j  0.+0.j]
# phi = [ 0.          0.          3.14159265 -0.          0.        ]

As you can see, the resulting a "rotates" from (0+j0) via (-0+j0), (0-j0) back to (0,0) just as it would for (1+j0), (0+j1), (-1+j0), (0-j1), (1+j0) with magnitude 1 instead of 0. Indeed for some operations like 1/np.inf and 1/-np.inf it may be important to keep the sign, but of course here it is a side-effect.

Nevertheless, the code works fine, since at 0 intensity the phase is meaningless (and in this case just jumps between 0, -0 and 2pi, which are really the same). I checked and the intensity will definitely be 0.0 after RectAperture(), no matter the order of functions. In fact, also in your hologram examples in the other issue posted, a lot of data points will have seemingly random phases, but at sufficiently low intensity this is meaningless and for plotting often phase-blanking is convenient (i.e. setting phase=0 or NaN when Intensity<threshold).