Wave propagation and interference phenomenom simulation with Matlab

Young's double-slit experiment


Animation, results, and observation

See and observe the following animation:



Figure is subdivised into 3 subplots (see using custom subplots),
  • left: whole simulation, wave propgation + double-slit
  • lower right: wave propagation after the slits
  • upper right: receiving screen (top of previous plots), showing interferences. Line curve is from theoritical computation.

Simulation & computation principles - Matlab code

Simulation and computationnal aspects rely on finite difference method, see for example (french) this page and also for example (also french) introduction to Matlab programming, and more precisely TP 9: Simulation bidimensionnelle de la propagation d'une onde, Diffraction par une simple ou double fente

Computation, simulation, and animation is performed using Matlab, from the following script:
% New beginning... clean...
clear all;close all
tic
c=10;

% Grid
Lx=100;
Ly=100;
T=15;

deltax=Lx/200;
deltay=Ly/200;
deltat=sqrt(deltax^2+deltay^2)/(2*c)
deltat=0.9*deltat

gax=(c*deltat)^2/(deltax^2);
gay=(c*deltat)^2/(deltay^2);

xx=[0:deltax:Lx];nx=length(xx);
yy=[0:deltay:Ly];ny=length(yy);
tt=[0:deltat:T];nt=length(tt);

Reflection_parois=0; % if 1, reflection on boundaries
                     % else absorbing waves on boundaries


Largeurx=1;Nlgx=round(Largeurx/Lx*nx);
posfx=50;nfx1=round(posfx/Lx*nx);nfx2=nfx1+Nlgx;

LargeurFy=2;NFy=round(LargeurFy/Ly*ny);
posF1y=Ly/2-LargeurFy/2;posF2y=Ly/2+LargeurFy/2;

cF1y=round(posF1y/Ly*ny);cF1y=96;
cF2y=round(posF2y/Ly*ny);cF2y=104;
nF1y1=cF1y-NFy/2;nF1y2=nF1y1+NFy/2;
nF2y1=cF2y-NFy/2;nF2y2=nF2y1+NFy/2;

Fentex=[nfx1:nfx2];
Fentey=[1:nF1y1,nF1y2:nF2y1,nF2y2:ny];

% Source location
sx=round(nx/4);sy=round((nF1y2+nF2y1)/2);

% Source frequency: 2 Hertz
sfrq=2;

Mat=zeros(nx,ny,nt);
% Mat Indexes, without boundaries:
i=[2:nx-1];j=[2:ny-1];

Mat(i,j,2)=Mat(i,j,1)+...
    gax/2*(Mat(i-1,j,1)+Mat(i+1,j,1)-2*Mat(i,j,1))+...
    gay/2*(Mat(i,j-1,1)+Mat(i,j+1,1)-2*Mat(i,j,1));


for k=3:nt-1
     Mat(i,j,k+1)=2*Mat(i,j,k)-Mat(i,j,k-1)+...
                gax*(Mat(i-1,j,k)+Mat(i+1,j,k)-2*Mat(i,j,k))+...
                gay*(Mat(i,j-1,k)+Mat(i,j+1,k)-2*Mat(i,j,k));
   
    Mat(sx,sy,k+1)=sin(2*pi*sfrq*k*deltat);
    
    if (Reflection_parois==0)
        Mat(1,:,k+1)=Mat(2,:,k);
        Mat(nx,:,k+1)=Mat(nx-1,:,k);
        Mat(:,1,k+1)=Mat(:,2,k);
        Mat(:,ny,k+1)=Mat(:,ny-1,k);
    end
    
    Mat(Fentex,Fentey,k+1)=0;
end
m=max(max(max(Mat)))
% Double-slit plot
Mat(Fentex,Fentey,:)=1;

% Plot command area: 
%  3 subplots: 
%    - left: whole simulation: wave propgation + double-slit
%    - lower right: wave after passing through the slits
%    - upper right: receiving screen, showing interferences
fig=figure(1);clf;whitebg('w');
set(fig,'position',[0 0 640 480]);

% Initialising Film variable
%  to contain all frames and to make finally a... film...
Film=[];

ecranx=200;
Interf=squeeze(Mat(ecranx,:,:));
MatEcran=Mat;
MatEcran(ecranx,:,:)=5;
MatEcranSup=MatEcran(101:end,:,:);

% Theoritical result, see below 
lbd=c/sfrq; % wave length
k=2*pi/lbd; % wave number
D=Lx/2;
D=(ecranx-max(Fentex))/200*Lx
d=(nF2y1-nF1y2)/200*Ly;
w=k*d/D;


for kk=1:1:nt
    subplot('position',[0.05 0.2 0.45 0.6])
    pcolor(squeeze(MatEcran(:,:,kk)));
    shading interp,axis off
    caxis([-0.4 0.4])
    
    subplot('position',[0.55 0.05 0.4 0.4])
    pcolor(squeeze(MatEcranSup(:,:,kk)));
    shading interp
    caxis([-0.003 0.003]);
    axis off
    
    subplot('position',[0.55 0.55 0.4 0.4])
    axis off;hold on
    % The very theory:
    plot(xx,(.75e-6)*(cos(w*(xx-49))).^2,'-r','linewidth',2);

    plot(xx,Interf(:,kk).^2)
    axis([xx(30) xx(end-30) 0 1e-6])
    title('{\bf \''Ecran}','fontsize',16,'interpreter','Latex')
    
    %Film=[Film getframe(fig)];
    pause(0.05)
end


%movie2avi(Film,'Film.avi','compression','none','fps',15)

Details on theoritical and mathematical computations: interference phenomenom



\[\psset{unit=1cm,arrowsize=5pt} \begin{pspicture}(-.8,-3)(4.6,3) \psline[linewidth=1.5pt](0,-3)(0,-.4) \psline[linewidth=1.5pt](0,-.2)(0,.2) \psline[linewidth=1.5pt](0,.4)(0,3) \psline[linecolor=blue]{<->}(-.4,-.3)(-.4,.3) \rput[r](-.6,0){\blue$d$} \psline[linecolor=blue]{<->}(0,-2)(4,-2) \rput[r](2,-2.3){\blue$D$} \psline[linewidth=1.5pt](4,-3)(4,3) \psline(0,-.3)(4,1.5)(0,.3) \rput(2,1.15){$d_1$}\rput(2,.3){$d_2$} \psline[linestyle=dashed](-.25,0)(4.5,0) \psline[linecolor=blue]{->}(4.2,0)(4.2,1.5) \rput[l](4.3,.8){\blue$x$} \end{pspicture}\]


Amplitude $A(x)$ on a point of the screen is found by adding amplitudes of the two secondary coherent (same frequency) sources:
\[A(x)(x)=A_1(x)+A_2(x)=I_0e^{i(\omega t-kd_1)}+I_0e^{i(\omega t-kd_2)}\]

where $I_0$ is the intensity of light, $\omega$ is the pulsation, $t$ is the time of observation (with respect to an abitrary origin, out of interest here), $k$ is the wave number, and $d_1$ and $d_2$ distances as shown on previous illustration.
This sum can be written
\[A(x)=I_0e^{i\omega t}\left( e^{-ikd_1}+e^{-ikd_2}\rp\]

then, as is common with complex exponentials, factorising with the mean argument:
\[\begin{array}{ll} A(x) &=I_0e^{i\omega t}e^{-ik(d_1+d_2)/2}\Bigl[ e^{ik(d_2-d_1)/2}+e^{ik(d_2-d_1)/2}\Bigr]\\[.8em] &=I_0e^{i\omega t}e^{-ik(d_1+d_2)/2}\times 2\cos\left( k\dfrac{d_2-d_1}{2}\right) \enar\]


Then, the intensity of a light field is proportionnal to the square of the modulus of its amplitude:
\[I(x)=|A(x)|^2 =4I_0^2\cos^2\left( k\dfrac{d_2-d_1}{2}\right) \]


$\Delta=d_2-d_1$ is the so-called phase difference, and the previous expression shows thus that intensity just depend on the difference of distance travelled by each wave.
We have geometrically, with above notations,
\[\begin{array}{ll} d_1&=\sqrt{\left( x-\dfrac{d}{2}\rp^2+D^2} \\[1em] d_2&=\sqrt{\left( x+\dfrac{d}{2}\rp^2+D^2} \enar\]


Assuming that the distance $D$ from slits to the screen is large compared with the width sepraration $d$ of the slits, we have, using finite expansion (or finite Taylor series expansion), at first order,
\[\begin{array}{ll} d_1&=\sqrt{x^2-dx+\dfrac{d^2}{4}+D^2}\\[1em] &=D\sqrt{1+\dfrac{x^2}{D^2}-\dfrac{d}{D^2}x+\dfrac{d^2}{4D^2}}\\[1em] &\simeq D\left( 1+\dfrac{x^2}{2D^2}-\dfrac{d}{2D^2}x\rp\\ &\simeq D+\dfrac{x^2}{2D}-\dfrac{d}{2D}x \enar\]

and the same for $d_2$:
\[ d_2\simeq D+\dfrac{x^2}{2D}+\dfrac{d}{2D}x \]


Finally, at first order,
\[d_2-d_1=\dfrac{d}{D}x\]

so that, intensity on the screen is given by
\[I(x)=4I_0^2\cos^2\lp\dfrac{kd}{D}x\rp\]


Intensity is finally proportionnal to:
\[\psset{xunit=1cm,yunit=2cm,arrowsize=7pt} \begin{pspicture}(-5,-.3)(5.5,1.4) \psline[linewidth=1pt]{->}(-5,0)(5.5,0) \psline[linewidth=1pt]{->}(0,-.2)(0,1.4) \psplot[plotpoints=200,linecolor=red,linewidth=1.5pt]{-5}{5}{x 180 mul 3.14 div cos 2 exp} \end{pspicture}\]

which is coded in the Matlab simulation and can be observed on top of the page

See also:
Top Lien