The diffraction results show unexpected

[foxwb]

I simulated a tilted planewave, the results depended on the grid dimension. Here is the source:

import matplotlib.pyplot as plt
from LightPipes import *
f = Begin( 20*mm, 1*um, 1000 );  # grid size is 20*mm, grid dimension is N=1000, wavelength is 1um
#c = PlaneWave( f, 1*mm );
c = CircAperture( f, 0.5*mm );
c = Tilt( c, 25*mrad, 0*mrad );    # the planewave tilts 25mrad in x direction
#c = Fresnel( c, 100*mm );
c = Forvard( c, 100*mm );
I = Intensity( c );
plt.imshow(I);
plt.title( "size=20mm, lambda=1um, N=1000, tx=25mrad" );
plt.show()

In this way, the result shows the circle source splits into to circles:

However, if I change N from 1000 to 500, the result is totally different, as shown:

compare the two pictures, the first one (N=1000) is different from we have seen in experiment, the second one is still no unexpected, because the light is not tilted( if the direction of the light tilts, the position should change).

I programmed another program directly derived from numpy.fft , and got the same result. I don't know why this happens.

[ldoyle]

Hi foxwb, thanks for the report. It seems you ran into a problem with aliasing. The field is stored as a 2D array of complex numbers, the magnitude (abs of number) is the electric field strength and the angle is the phase. Due to the finite sampling, this can lead to ambiguous fields, since there is mathematically no way of differentiating between a phase value X and X+2pi, X+4pi, and so on. (Compare the image of an aliased/ambiguous sine wave on https://en.wikipedia.org/wiki/Aliasing )

In other words, your tilt is too steep for the grid spacing you have chosen. What makes things worse (or interesting from a learning perspective) is that your grid size, wavelength, spacing and tilt coincidentally match up so that each phase tilt on your discrete grid points amounts to a multiple of pi. For N=1000, this creates an alternating phase of 0 and pi, for N=500 every phase value amounts to 0 (or 2pi, but this is the same mathematically/numerically), which gives the same result as an untilted field!

I took your script and experimented with different tilt angles (with N=100 to show more clearly the phase jumps). See the results here for tx=0.25, 2.5, 5, 5.25mrad: As you can see, at tx=5mrad the tilt is indifferentiable from 0 tilt, and 5.25mrad looks just like 0.25. In other words, for the given grid spacing and wavelength, a 2pi periodicity is matched with a 5mrad periodicity. For odd combinations, e.g. N=1024 or lambda=633nm or 25.33mrad, it would seem that the tilt worked but it might still be too small, since a large part got aliased!

Experimenting further, I created a video scanning a range of tilts (this time for N=1000 as in your original question). Here you can see that around the critical values, interesting things happen.

Final thoughts:

• If you change from Forvard() to Fresnel(), the problem seems to get worse. In general it is always good to test both and see which seems better in the current situation, as this always depends on propagation length, size of zero-intensity borders and so on.
• I experienced the same thing in the past when simulating a lens with too fast focusing (short focal length compared to diameter).
• I am sure there is a formula to estimate when the calculation breaks down (for aliasing in general there is the Nyquist limit, not sure how to transfer this here correctly, though), but for me the easiest solution is to start with really small tilt/distortion, and gradually increase the magnitude but stop when the phase jump between two neighbouring field values approaches Pi.
• When you have reached this limit, make the field spacing finer by increasing N or making the field size smaller (if possible).
• If you hit a limit there where N is just too large, e.g. above 4000, try to reformulate the physical problem. For my short focal length lenses case, the LightPipes manual provides a better way of doing this with the LensForvard etc. For a tilt of a single beam, maybe the optical axis can be "redefined" by just leaving out the tilt and re-applied later by remembering/calculating the resulting offset (this of course depends on the concrete experiment/problem being simulated)
• Maybe we can add a section about this pitfall in the documentation.
• The scripts to generate the images and animation will be added to the Examples, for the moment you can find them in https://github.com/opticspy/lightpipes/tree/develop_LenDo/Examples/tests as tilt_test.py and tilt_test_anim.py

[FredvanGoor]

I had the same answer as Idoyle (but much shorter however!) The problem is indeed caused by the fact that the maximum spatial frequency of the system is larger than the maximum spatial frequency allowed. This maximum is determined by the distance between the grid points. It is a good idea to pay attention to this problem in the manual because it is often encountered by users. The difference between the Forvard and Fresnel propagation routines is roughly far- (Forvard) and near field (Fresnel). A useful criterium is the so-called Fresnel number: FN=a^2/(L*lambda). If it is larger than 1 it is the near field and lower than one the far field. a = the beam size (aperture radius), L is the distance to be propagated.

[foxwb]

THANKS ldoyle: After reading your answer, I think I have gotten the answer. The issue came from the coincidence between the tilt angle and grid. To show the process clearly, I programed to illuminate the relation between the position and the phase. Here is the source code:

import numpy as npy
import matplotlib.pyplot as plt;

size=20E-3    # grid size
lam = 1E-6    # wavelength
tx = 25E-3    # tilt angle
K = 2*npy.pi/lam

def foo( z, N ):  # function foo is used to calculate phase
    delta_z = tx * z;    # optical path difference
    return ( K * delta_z) % (2*npy.pi) 

for i,j in enumerate([1001,931,773,501]):  # grid N=1001, 931, 773, 501
    x = npy.linspace( -size/2, size/2, j )
    y = foo(x,j)
    plt.subplot( 2,2,i+1 )
    plt.plot( x,y,'.-')
    plt.xlabel( 'position')
    plt.ylabel( 'phase')
    plt.title( f'N={j}')

plt.show();

When the tilt angle is 25mrad, I tested different grid dimension(N=1001, 931, 773, 501), the result is shown as the figure below: From the figure, we can get the item of phase has only 3 values(exactly 2 values): 0, Pi, and 2Pi when N=1000, and 2 values(0, 2Pi) when N=500. But when N=931 and 773, more informations about phase could be gotten.

OK, We now know why this happens. Howevery, I think it is difficult to avoid the problem, especially when we simulate in sophisticated condition. Can we try to solve the problem?

[foxwb]

I agree the suggestion from Mr. FredvanGoor, the Fresnel number may help us to assess system and to choose the suitable grid dimension.