Diffraction Grating using Fraunhofer diffraction approximation

Question

Minh Hoang 2024-1-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2068016-diffraction-grating-using-fraunhofer-diffraction-approximation

评论： Minh Hoang 2024-1-9

Greetings everyone, I am trying to make a diffraction grating pattern from an amplitude transmittance function given by: tA=1/2(1-cos(2*pi*x/P))*rect(x/D1)*rect(y/D1). I learned this method from "Computational Fourier Optics: a MATLAB Tutorial" by David G.Voelz (this book is available online: https://www.spiedigitallibrary.org/ebooks/TT/Computational-Fourier-Optics-A-MATLAB-Tutorial/eISBN-9780819482051/10.1117/3.858456?SSO=1#_=_ ). In this book it was mentioned that "When illuminated by a unit amplitude plane wave, the source field is U1(x1,y1) = tA(x1,y1)". I wonder in case the incident light is a Gaussian beam then can I consider this incident beam a plane wave as well? Or should I make any change to the source field U1?

Here is the code I learned from the book:

First I created a function for far-field approximation:

function[u2,L2]=propFF(u1,L1,lambda,z)
%propagation - Fraunhofer pattern
%assumes uniform sampling
%u1 - source plane field
%L1 - source plane side lenght
%lambda - wavelength
%z - propagation distance
%L2 - observation plane side length
%u2 - observation plane field
[M,N]=size(u1); %get input filed array size
dx1 = L1/M; %source sample interval
k=2*pi/lambda;
L2=lambda*z/dx1;    %obs sidelength
dx2=lambda*z/L1;    %obs sample interval
x2=-L2/2:dx2:L2/2-dx2;  %obs coords
[X2,Y2]=meshgrid(x2,x2);
c=1/(1i*lambda*z)*exp(1i*k/(2*z)*(X2.^2+Y2.^2));
u2=c.*ifftshift(fft2(fftshift(u1)))*dx1^2;
end

Then I use this approximation to plot diffraction grating pattern:

lambda= 0.5e-6; %wavelength
f=0.5;  %propagation distance
P=1e-4; %grating period
D1=1.02e-3; %grating side length
L1=1e-2;    %array side length
M=500;  %number of samples
dx1=L1/M;
x1=-L1/2:dx1:L1/2-dx1;  %source coordinates
[X1,Y1]=meshgrid(x1,x1);
% Grating field and irradiance
u1=1/2*(1-cos(2*pi*X1/P)).*rect(X1/D1).*rect(Y1/D1);
% Fraunhofer pattern
[u2,L2]=propFF(u1,L1,lambda,f);
dx2=L2/M;
x2=-L2/2:dx2:L2/2-dx2; 
y2=x2;
I2=abs(u2).^2;
figure(1)
plot(x2,I2(M/2+1,:));
xlabel('x(m)');
ylabel('Irradiance');