Open IvanOstr opened 3 years ago
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:
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
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
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.
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.
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.
![image](https://user-images.githubusercontent.com/39938309/126974929-8f054d0d-4f35-4f66-a2b9-d2c39b31278c.png)
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.