Diffraction From Thin Wire

[IvanOstr]

Hey there,
I'm trying to simulate diffraction from a thin wire.
Unfortunately, I'm not able to get the famous sinc from 1, 2.
It seems increasing the grid or spacing does not improve the results.
Slit diffraction worked perfectly.
Attaching pictures and code below.

from LightPipes import*
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams.update({'figure.max_open_warning': 0})
import LightPipes

def sinc(x, x0, d, l, wl):
    b = np.pi * d * (x -x0)/ (wl * l)
    return (np.sin(b) / b)**2

wavelength = 800*nm
k = 2*np.pi/wavelength
size = 30*mm
N = 5000
N2 = int(N/2)
spacing = size/N
n = (1.000293 + 0j)*np.ones((N,N))
center = round(N/2)
positions = np.linspace(0, size, N)

F = Begin(size, wavelength, N)
F = GaussAperture(F, 6*mm)
# F = PlaneWave(F, 20*mm)
F2 = LightPipes.Field.copy(F)
# F = RectAperture(F, 30*mm, 3*mm)
# F = CircAperture(F, 10*mm)

plt.figure()
plt.plot(np.array(Intensity(F, 0)[:, center]))
plt.title('Before Cut')

d = 0.1*mm
F = RectScreen(F, 25*mm, d)

plt.figure()
plt.plot(np.array(Intensity(F, 0)[:, center]))
plt.title('After Cut')

plt.figure()
plt.imshow(Intensity(F, 0))
plt.title('After Cut')

for i in np.linspace(50, 50, 1):
    F_prop = LightPipes.Field.copy(Fresnel(F, i * cm))
    F_no_wire = LightPipes.Field.copy(Fresnel(F2, i * cm))

    plt.figure(5)
    signal = Intensity(F_prop, 0)[:, center] - Intensity(F_no_wire, 0)[:, center]
    plt.plot(positions, signal/max(signal))
    plt.plot(positions, sinc(positions, size/2, d, i*cm, wavelength))
    plt.title(str(i) +'cm propagation')
    # plt.figure(6)
    # plt.plot(positions, Intensity(F_prop, 0)[center, :] / max(Intensity(F_prop, 0)[center, :]))
    # plt.figure()
    # plt.imshow(Intensity(F_prop, 0))
    # plt.title('After Cut')

Do you have any ideas why?

Update:
I've found more about the problem and updated the code - the problem is the oscillations in the sinc function. Notice I subtracted the gaussian without the wire to get the picture below.

[ldoyle]

Dear Ivan,
sorry that no-one got back to you quicker, I personally have been very busy in the past months and could not look after LightPipes much.
This is indeed an interesting problem and after I found no obvious errors and could reproduce the problem, I had to investigate longer. The worst part was I trusted all calculations to be correct (except for small grid effects etc), but I knew something must still be wrong since I have seen the experiment in real life (measuring the thickness of a human hair is a cool demo for a faculty open day).

To figure this out, I built a small script based on your code which compares the diffraction pattern of the screen and the complementary aperture automatically. According to Babinet's principle, they should be equal. Here are my conclusions:

• Not the intensity pattern as such must be equal, only the part created by the diffraction pattern related to the aperture/screen. This is mixed in with the intensity/field of the unperturbed beam. I'm sure it is simple enough to write this down/ understand without simulation, I have not tried it, though.
• I therefore tried to replicate as closely as possible the real parameters for the thin wire experiment. It turns out the beam diameter is the crucial part, setting a large beam diameter just "hides" the diffraction pattern behind much higher intensity on axis.
• The final beam diameter chosen is so small that the wire cuts ~10% of the beam intensity.
• Even then, the pattern is barely visible until you exaggerate it by "mimicking" an overexposed camera and applying some nonlinearity to exaggerate low intensity regions.
• At these low intensity levels, the grid does have an effect, e.g. some maxima vanish with the wrong combination of N/size (e.g. 1024, 30mm), they reappear with doubling the grid size (2048, 30mm). Just be aware of this while playing around.
• Beware that the wire is only a few grid points in size (in my sample script, 8 points), so "measuring" the thickness of the wire by measuring the distance of the diffraction maxima will have an error of ~1/8 of the actual thickness.

Here is an example result and the script:

"""

Test of Babinet's principle:
    "the diffraction pattern from an opaque body is identical
    to that from a hole of the same size and shape except for
    the overall forward beam intensity"[1]

As can be seen by playing around with the parameters here, the
effect may be well hidden when the overall beam effect is much larger.

It seems to me the best results may be found when the intensity that is
blocked is the same order of magnitude as the intensity that is passed by the 
complementary object, i.e. the object should obstruct 10-50% of the intensity.

[1] https://en.wikipedia.org/wiki/Babinet%27s_principle

Created on Mon Aug  2 12:31:33 2021

@author: Leonard.Doyle
"""

import numpy as np
import matplotlib.pyplot as plt

import LightPipes as lp
from LightPipes.units import * #mm, nm, PI, ...

wavelength = 532*nm
size = 25*mm
N = 2048
dx = size/N
center = round(N/2)
positions = np.linspace(0, size, N)

F0 = lp.Begin(size, wavelength, N)

#*** define beam ****
Fbeam = F0
Fbeam = lp.GaussAperture(Fbeam, .5*mm)
# Fbeam = lp.CircAperture(Fbeam, .3*mm) #also kinda works, Gauss looks better

#*** define aperture/screen/object ****
d = 0.1*mm
print('Number of grid points representing wire:',int(d/dx))
Fobj = lp.RectScreen(F0, d, 2*size)
Iobj = lp.Intensity(Fobj, 0)
Icomplementobj = 1-Iobj

#*** calculate complementary field ****
Fobj = lp.MultIntensity(Fbeam, Iobj)
Fcompl = lp.MultIntensity(Fbeam, Icomplementobj)

#*** propagate field ***
z = 0.5*m
Fprop_obj = lp.Forvard(Fobj, z) #Fresnel also works fine
Fprop_compl = lp.Forvard(Fcompl, z)

#**** Plot results ****
Iobj = lp.Intensity(Fobj)
Iprop_obj = lp.Intensity(Fprop_obj)
Icompl = lp.Intensity(Fcompl)
Iprop_compl = lp.Intensity(Fprop_compl)

Ibeam = lp.Intensity(Fbeam, 0)
print('Intensity blocked by object: {:.1f}%'.format(
    100*(1-Iprop_obj.sum()/Ibeam.sum())))
print('Intensity blocked by complementary aperture: {:.1f}%'.format(
    100*(1-Iprop_compl.sum()/Ibeam.sum())))

#visually exaggerate features to mimic overexposed camera
# most photographs of laser beams are actually heavile overexposed
# and possibly also show nonlinear response a little
enhance_powerlaw = 0
enhance_powerlaw = 1/3 #comment/uncomment to switch quickly
if enhance_powerlaw:
    Iprop_obj = Iprop_obj**enhance_powerlaw
    Iprop_compl = Iprop_compl**enhance_powerlaw

exaggerate_brightness = 0
exaggerate_brightness = 3 #comment/uncomment to switch quickly
if exaggerate_brightness:
    # clip image at some brightness to scale up brightness N times
    Iprop_obj[Iprop_obj>Iprop_obj.max()/exaggerate_brightness] = \
        Iprop_obj.max()/exaggerate_brightness
    Iprop_compl[Iprop_compl>Iprop_compl.max()/exaggerate_brightness] = \
        Iprop_compl.max()/exaggerate_brightness


fig, [(axtl, axtr), (axbl, axbr)] = plt.subplots(2,2, sharex=True, sharey=True)
fig.suptitle('Test of Babinet\'s principle')
axtl.set_title('Object'), axtl.imshow(Iobj)
axtr.set_title(f'Object, propagated to z={z:.1f}m')
axtr.imshow(Iprop_obj)
axbl.set_title('Complement object'), axbl.imshow(Icompl)
axbr.set_title('Complement, propagated same dist')
axbr.imshow(Iprop_compl)

All in all, I cannot claim I understand what is going on, but at least there are parameter sets with which the results make sense. Do let me know if you have any further insights.
Best, Lenny

[FredvanGoor]

Hi Lenny and Ivan,
Great answer Lenny! The borders of the grid can cause (numerical) reflections when the intensity at the border is not zero. Could that be the cause of the problems observed?
Once this demo of Babinet’s principle is ready I could put it in the LightPipes website. Do you agree?

Fred van Goor

[IvanOstr]

Fred - Unfortunately getting a bigger grid doesn't matter from what I tried.
Lenny - Great answer! I'm not sure I understand the non linearity of the camera part. I may check this in COMSOL and notify is if there is something interesting.

[IvanOstr]

After giving it a small thought what Fred says sounds in line with the fact that the best results are obtained with the smaller gaussian beam.