姓名 王锐 班级 12级电子信息工程二班 学号 1208421042 Matlab进行波动光学的计算机模拟

简述：在教学和科研中，用Matlab进行计算机模拟正越来越被重视。分析讨论了夫琅禾费衍射以及双光束,多缝,牛顿环和干涉平面光栅的衍射等理论，用Matlab编写出相应程序后，进行了计算机模拟，这有助于理解和研究衍射和干涉理论。

1功能介绍：光的波动性常常表现为干涉、衍射、偏振等现象。虽然关于这些现象的描述及其阐述有好多，但是未能配以精彩的直观形象图形。运用MATLAB获得模拟图形，将定性的语言描述和抽象的数学表示变为生动直观的表现，可以进一步分析和描述有关波动光学的现象和规律等理论。

光的衍射现象研究,不仅具有重要的理论意义,而且在光学仪器研制和成象分析等均有重要价值.针对衍射计算困难的问题,选择合适的计算公式,并基于科学计算软件MATLAB5.3编写计算程序,使计算得变得简洁,大大提高了实验的效率.用MATLAB编程得出的计算结果与用数码相机记录的衍射图样进行比较,两者吻合得很好。

2 功能实现：运用MATLAB获得模拟图形，将定性的语言描述和抽象的数学表示变为生动直观的表现，可以进一步分析和描述有关波动光学的现象和规律等理论。使用M文件形式。M文件有两个形式：一种称为命令文件（scriptfile）就好像dos下的批处理文件一样。这类程序包含了一连串的MATLAB命令，执行时依次执行。另一种称为函数文件（functionfile）。它的第一句可执行语句是以function引导的定义语句，在函数文件中的变量都是局部变量。

1

姓名 王锐 班级 12级电子信息工程二班 学号 1208421042

2.1 夫琅禾费衍射：

夫琅禾费衍射，是认为光源和观察屏离衍射屏（孔）处于无穷远处衍射现象。实验装置如图

S为单色点光源，放置在透镜L1的物方焦点处，所得平行光垂直入射到障碍物，借助于透镜L2将无穷远处的衍射图样移至L2的像方焦面上观察。

将衍射屏制作成输入图像，用imread()函数读入，然后利用傅里叶变换函数fft2()对其进行傅里叶变换，得到其傅里叶频谱。由函数fft2()实现的傅里叶变换频谱的直流分量位于图像的左上角，而由透镜实现的光学傅里叶变换的直流分量位于图像中心。因此，为了得到模拟的光学傅里叶变换，需调用函数fftshift()将零频移到频谱中心。

2

姓名 王锐 班级 12级电子信息工程二班 学号 1208421042

2.2双光束干涉:

双光束干涉，有分波面类型的杨氏干涉、劳埃镜干涉、菲涅耳双棱镜和双面镜干涉等，还有分振幅类型的麦克耳孙干涉等。

下图所示的是双缝干涉装置，是分波阵面的双光束干涉的典型实例。下面分析它的干涉原理与干涉条纹特点。

3

姓名 王锐 班级 12级电子信息工程二班 学号 1208421042

如图7.1.1所示，两狭缝间距为d，双缝所在平面与屏幕平行，两者之间的垂直距离为D，O为屏幕上的坐标原点且与两狭缝对称。两个狭缝光源满足振动方向相同、频率相同、相位差恒定的相干条件。故两列光在空间相遇，将产生干涉现象。屏幕上将出现明暗相间的干涉条纹。设OP=y，则由几何关系可知，两个相干光源到达屏幕上任意点P的距离分别为：

r1= D2+ y?d2^2 r2= D^2+ y+d2 ^2

这样两列相干光的光程差为?r=r2?r1，相位差为?φ=2π?r/λ。设两狭缝光源的光波单独到达屏幕P点处的振幅分别为E1和E2，光强分别为I1和I2。则两

列光波叠加后的振幅为E= E1^2+E2^2

+2E1*E2cos?φ：，而叠加后的光强为：

I=I1+I2+2 I1*I2cos?φ。设两列光波在屏幕上相遇点振幅相等，则P点光强为：I=4I0 cos

?φ2^2 ?φ=2kπ,k=0,1,2,3?时为干涉明条纹。 ?φ= 2k+1 π,k=0,1,2,3?时为干涉暗条纹。

4

姓名 王锐 班级 12级电子信息工程二班 学号 1208421042

其实绝对的单色光只是理论上存在，实际的“单色光”都是有一定谱线宽度的准单色光，这种准单色光入射到干涉装置后，其中的每一种波长成分都将产生自己的干涉条纹，由于波长不同，除零级明条纹外，其他级次的条纹将彼此错开，并发生不同级条纹的重叠。在重叠处总的光强为各种波长的条纹的光强的非相干相加。随着级次的增大，干涉条纹的明暗对比减小，当级次增大到某一值后，干涉条纹就消失了。对于谱线宽度为?λ的准单色光，干涉条纹消失的位置应是波长 λ+?λ2的成分的k级明条纹与波长为λ??λ2的成分的k+1级明纹重合的位置。由于两成分在此位置上有同一光程差，根据光程差与明纹级次的关系可知，条纹消失时，应满足：(λ+?λ2)k= λ??λ2 k+1 。由于?λ?λ，于是可得：k=λ

?λ由此可知，?λ愈大，即光的单色性愈差，能观察到得干涉条纹级次就愈小。

clearall

a=-4*pi:0.01*pi:4*pi;%通过改变PI的倍数而改变条纹数

P=1-sin(2*a).^2./sin(a).^2;

%当要求P的曲线分布图时P=sin(2*a).^2./sin(a).^2plot(a,P)

lgray=zeros(256,3);fori=0:255

lgray(i+1,:)=(255-i)/255;endimagesc(P)colormap(lgray)

5

姓名 王锐 班级 12级电子信息工程二班 学号 1208421042

2.3牛顿环:

将一个球面曲率很大的平凸透镜A的面紧贴在一块平板玻璃B上，球面和平面相切于O点，在它们之间形成一空气薄膜，从切点到边缘膜的厚度逐渐增加。在以切点为中心的圆周上，空气层的厚度相等，将它放在观察等厚干涉的实验装置中，让单色光正投射于其上，从反射光和透射光中都可以观察到以O为中心的一系列明暗相间的同心圆环。牛顿环属于分振幅干涉类型。

6

姓名 王锐 班级 12级电子信息工程二班 学号 1208421042

使用imshow建立灰度图像的使用格式：imshow(I,N) 。其中参数N

为正整数，指定灰度的层次，当缺省该参数时，系统默认为256级的

灰度级；参数I为数值矩阵，imshow的作用就是将数值矩阵I的元素

值用N个灰度级的黑白图像可视化显示出来。实际上是在数值矩阵I

和N个灰色调之间建立了一种颜色映射关系：I当中元素值最大者映

射为白色（将该元素值作为白色显示），元素值最小者映射为黑色（将

该元素值作为黑色显示），元素值介于最大和最小之间的则按照某种

约定的规则映射到其它的灰度级（显示为不同灰度的灰色）。

为了方便叙述，假定再现的图像尺寸2mm×2mm,使用上述指令可

以很方便的将牛顿环干涉条纹在该区域内再现。

首先，利用(4)式获取干涉光强I在该区 （0.0010.001,0.0010.001 xy

假定观察屏是xy平面）的数值分 布 x=linspace(-0.001,0.001,200); y=linspace(-0.001,0.001,200);

[X,Y]=meshgrid(x,y); %将xy平面2mm×2mm的区域分割为200×200

的网

%格（像素），矩阵X、Y分别输出格点的x和y坐标

r2=X.^2+Y.^2 I=abs(sin(pi*r2/R/)).^2;

7

姓名 王锐 班级 12级电子信息工程二班 学号 1208421042

% 200×200的数值矩阵；计算格点上光强I,得到光强的数值分布 接下来使用imshow指令将数值矩阵I可视化

imshow(I)

2.4 平面光栅的衍射：

有大量等宽度，等间距的平行狭缝组成的光学系统称为衍射光栅。单缝宽度a和刻痕宽度b之和称为光栅常数d，d=a+b。观察到的实验现象中衍射图样的强度分布具有如下一些特征：

（1）与单缝衍射图样相比，多缝衍射的图样中出现一系列新的强度最大值和最小值，其中那些较强的亮线称为主最大，较弱的亮线叫做次最大。

（2）主最大的位置与缝数N无关，但它们的宽度随N的增大而减小，其强度正比于2N

（3）相邻主最大之间有N-1条暗纹和N-2个次最大。（注：具体情形如图4-6的a、b所示）（4）强度分布中保留了单缝衍射的因子，那 8

姓名 王锐 班级 12级电子信息工程二班 学号 1208421042

就是曲线的包迹与单缝衍射强度曲线一样。光栅衍射条纹，理应是单缝衍射和缝间干涉的共同结果。关于上面论述的知识点，我们所学的教材[4]只给出了平面光栅衍射的强度分布曲线，仅由光栅衍射此强度函数曲线图，读者会还是会感觉抽象，难以对此现象形成更为直接的直观印象，特别是平面光栅衍射中涉及的缺级现象更是难以理解（课本中只给了以表格形式表示的光强分布公式）。本文下一章利用Matlab编写相应的程序，对此现象进行了模拟。由图4-7a和4-7b

（d=4b，N=6）的情况可知：j=4k（k=1,2,…）时对应的亮条纹消失，也就是缺级。

单缝衍射的Matlab编程

运用Matlab编制单缝衍射程序(张智星,2002),程序运行后在坐标区可以生成单缝夫琅和费衍射图样。

为了满足程序的普遍性和通用性,主程序中共设置4个输入参数,分别为光波波长、透镜焦距、单缝宽度、最大坐标范围，长度以mm为单位。主程序如下：

clear all

d=-4*pi:0.001*pi:4*pi;b=d./15;

P=1-(sinc(b).*sin(4*d)./sin(d)).^2;

%当要求P的曲线分布图时P=(sinc(b).*sin(4*d)./sin(d)).^2;4是N值可调plot(d,P);

lgray=zeros(100,3);for i=0:99

lgray(i+1,:)=(99-i)/99;

9

姓名 王锐 班级 12级电子信息工程二班 学号 1208421042 endimagesc(P)

colormap(lgray)

3 程序总结：

（1）本MATLAB程序具有方便的数据可视化功能，以将向量和距阵用图形表现出来。高层次的作图包括二维和三维的可视化、图象处理、动画和表达式作图。可以方便的看出光学模拟图形，有利于研究。

（2）程序简洁，逻辑性强，操作方便，提示信息较多，一般不会出现太大的错误。

（3）功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ,为用户提供了大量方便实用的处理工具。

（4）但是MATLAB也存在着一些不足之处，MATLAB的缺点是，它和其他高级程序相比，程序的执行速度较慢，由于MATLAB的程 10

姓名 王锐 班级 12级电子信息工程二班 学号 1208421042

序不用编译等预处理，也不生成可执行文件，程序为解释执行，所以速度相对较慢。由此可见，Matlab的语言执行效率是比较低的。

4课程总结：

通过建立实验的仿真平台，使那些过于抽象光学概念和不能用数学方法得到解析解复杂公式变得清晰简单；由于图形用户界面的应用，提高到了学习的可视性和可操作性，使学生可以自行灵活地修改光学类型和相关的参数，实现人机交互，使学习过程变得轻松自由,且能让学生对艰深难懂的光学概念和公式有一个更深的理解。在教学方面，它打破了课时的限制和光学实验条件的苛刻性，使得光学演示实验能在课堂中轻松进行，大大提高了教学效果，使教学变得更加主动和开放。

在光学领域里，计算机仿真技术也在发挥着愈来愈重要的作用。特别是在光学教学过程中，对光学现象的理解离不开光学实验。采用计算机仿真技术对其进行仿真，得到满足各种条件的光学实验结果图像，有助于学生在学习过程中建立清晰正确的观念。但是MATLAB也存在着一些不足之处，MATLAB的缺点是，它和其他高级程序相比，程序的执行速度较慢，由于MATLAB的程序不用编译等预处理，也不生成可执行文件，程序为解释执行，所以速度相对较慢。由此可见，Matlab的语言执行效率是比较低的。

大学教学课程中引入计算机模拟技术正日益受到重视，用MATLAB软件做光学实验的模拟，只需要用数学方式表达和描述， 11

姓名 王锐 班级 12级电子信息工程二班 学号 1208421042

省去了大量繁琐的编程过程。下面来介绍利用MATLAB进行光学模拟的方法。

在利用计算机来实现傅里叶变换的过程中，由于计算机中图像的存储使用是数字形式，连续的傅里叶变换不适用于计算机的处理，所以在计算机中的傅里叶变换一般都采用离散傅里叶变换（DiscreteFourierTransform--DFT）。在计算机中离散傅里叶变换主要是基于以下几个原因：

（1）离散傅里叶变换的输入输出都是离散值，方便计算机的运算操作。

（2）采用离散傅里叶变换，可以用快速傅里叶变换（FFT）算法来实现，来提高计算速度。

为了降低运算量，在计算离散傅里叶变换时，通常采用快速傅里叶变换来实现，其算法思想：

（1）首先将原图像进行转置。

（2）按行对转置后的图像矩阵做一维FFT，将变换后的中间矩阵再转置。

（3）对转置后的中间矩阵做一维FFT，最后得到的就是二维FFT。 物体图像的生成可以直接由矩阵运算生成，也可利用Windows下的画图工具，生成一副黑白图像，可调用命令函数imread()输入图像，输入的图像是一个巨大的二维矩阵，利用MATLAB的fft2（）命令对该矩阵进行二位离散傅里叶变换，得到图像的频谱，该频谱是一个复数矩阵，然后用取模函数abs（）对该复数矩阵取模，得到振幅谱 12

姓名 王锐 班级 12级电子信息工程二班 学号 1208421042

矩阵，利用函数fftshift（）对取模后的矩阵进行频谱位移，是直流分量移到频谱中心，从而使FFT频谱可视效果与实际图像相吻合。最后利用imshow（）函数将图像显示出来。

附录（部分程序）

夫琅禾费衍射：

clear;N=1;

tz=linspace(0,0.01*pi,1000*pi);

Pq=[];

K=6;%控制参数

[x,y]=meshgrid(linspace(0,N+1,800));z=x+i*y;

U=0;

For m=1:N;

For n=1:N;

zk=abs(z-[m+n*i])*K;

U=U+0.1*besselj(4,zk)./zk；

r=1-U;

A=1-abs(U).^2;end

Ip=imshow(A,[])

13

姓名 王锐 班级 12级电子信息工程二班 学号 1208421042

双光束干涉:

clearall

a=-4*pi:0.01*pi:4*pi;%通过改变PI的倍数而改变条纹数P=1-sin(2*a).^2./sin(a).^2;%当要求P的曲线分布图时

P=sin(2*a).^2./sin(a).^2

plot(a,P)；

lgray=zeros(256,3);

for i=0:255

lgray(i+1,:)=(255-i)/255;

endimagesc(P)

colormap(lgray)

牛顿环:

closeall;

figure('Position',[90164873483]);L=632.8;R=5.1;H=5;

a1=axes('Position',[0.4,0.16,0.4,0.7]);

[x,y]=meshgrid(linspace(-0.005,0.005,200));

r2=(x.^2+y.^2);

Di=[2*H+2*(R-sqrt(R^2-r2))*1e9]/L;

In=abs(cos(Di*pi*2));cr=abs(L-560)/200;cg=1-cr;

cb=abs(L-600)/240;Ik(:,:,1)=In*cr;Ik(:,:,2)=In*cg;Ik(:,:,3)=In*cb; 14

姓名 王锐 班级 12级电子信息工程二班 学号 1208421042 Pc=imshow(Ik,[]);

title('牛顿环模拟图','fontsize',18);

Lt=uicontrol(gcf,'style','text',...

'unit','normalized','position',[0.06,0.86,0.21,0.06],...

'BackgroundColor',0.7*[1,1,1],'ForegroundColor',[0.8,0.1,0.9],...'string','波

长:632.8nm','fontsize',16,'fontname','timesnewroman');s1=uicontrol(gcf,'style','slider',...'unit','normalized','position',[0.06,0.76,0.21,0.04],...

'BackgroundColor',0.7*[1,1,1],'ForegroundColor',[0.1,0.1,0.9],...'SliderStep',[0.01,0.01],'value',(632.8-360)/400,...'callback',['L=get(s1,''value'')*400+360;',...'set(Lt,''string'',[''波

长:'',num2str(L/1),''nm'']);',...'Di=[2*H+2*(R-sqrt(R^2-r2))*1e9]/L;',... 'In=abs(cos(Di*pi*2));

cr=abs(L-560)/200;

cg=1-cr;',...'cb=abs(L-600)/240;Ik(:,:,1)=In*cr;Ik(:,:,2)=In*cg;',...'Ik(:,:,3)=In*cb;set(Pc,''CData'',Ik);']);uicontrol(gcf,'style','text',...

'unit','normalized','position',[0.04,0.81,0.08,0.04],...

'BackgroundColor',0.8*[1,1,1],'ForegroundColor',[0.1,0.1,0.9],...'string','360','fontsize',16,'fontname','timesnewroman');uicontrol(gcf,'style','text',... 'unit','normalized','position',[0.22,0.81,0.08,0.04],...

'BackgroundColor',0.8*[1,1,1],'ForegroundColor',[0.1,0.1,0.9],...'string','760','fontsize',16,'fontname','timesnewroman');Rt=uicontrol(gcf,'style','text' 15

姓名 王锐 班级 12级电子信息工程二班 学号 1208421042 ,...

'unit','normalized','position',[0.06,0.66,0.23,0.06],...

'BackgroundColor',0.7*[1,1,1],'ForegroundColor',[0.8,0.1,0.9],...'string','曲率半径5m:','fontsize',16,'fontname','timesnewroman');

s2=uicontrol(gcf,'style','slider',...'unit','normalized','position',[0.06,0.56,0.21,0.04],...

'BackgroundColor',0.7*[1,1,1],'ForegroundColor',[0.1,0.1,0.9],...'SliderStep',[0.01,0.01],...

'callback',['R=get(s2,''value'')*7+5;',...

'set(Rt,''string'',[''曲率半

径:5m'',num2str(R),''m'']);',...'Di=[2*H+2*(R-sqrt(R^2-r2))*1e9]/L;',...'In=abs(cos(Di*pi*2));

cr=abs(L-560)/200;cg=1-cr;',...'cb=abs(L-600)/240;Ik(:,:,1)=In*cr; Ik(:,:,2)=In*cg;

',...'Ik(:,:,3)=In*cb;set(Pc,''CData'',Ik);']);

uicontrol(gcf,'style','text',...

'unit','normalized','position',[0.04,0.61,0.08,0.04],...

'BackgroundColor',0.8*[1,1,1],'ForegroundColor',[0.1,0.1,0.9],...'string','5','fontsize',16,'fontname','timesnewroman');

uicontrol(gcf,'style','text',...

'unit','normalized','position',[0.22,0.61,0.08,0.04],...

'BackgroundColor',0.8*[1,1,1],'ForegroundColor',[0.1,0.1,0.9],...'string','1 16

姓名 王锐 班级 12级电子信息工程二班 学号 1208421042 2','fontsize',16,'fontname','timesnewroman');

Ht=uicontrol(gcf,'style','text',...

'unit','normalized','position',[0.06,0.46,0.23,0.06],...

'BackgroundColor',0.7*[1,1,1],'ForegroundColor',[0.8,0.1,0.9],...'string','空气层厚度:5nm','fontsize',16,'fontname','timesnewroman');

s3=uicontrol(gcf,'style','slider',...

'unit','normalized','position',[0.06,0.36,0.21,0.04],...

'BackgroundColor',0.7*[1,1,1],'ForegroundColor',[0.1,0.1,0.9],...'SliderStep',[0.01,0.01],'value',0.05,...'callback',['H=get(s3,''value'')*100;',... 'set(Ht,''string'',[''空气层厚度:'',num2str(H),''nm'']);

',...'Di=[2*H+2*(R-sqrt(R^2-r2))*1e9]/L;',...

'In=abs(cos(Di*pi*2));cr=abs(L-560)/200;

cg=1-cr;',...'cb=abs(L-600)/240;

Ik(:,:,1)=In*cr;

Ik(:,:,2)=In*cg;',...'Ik(:,:,3)=In*cb;

set(Pc,''CData'',Ik);']);

uicontrol(gcf,'style','text',...

'unit','normalized','position',[0.04,0.41,0.08,0.04],...

'BackgroundColor',0.8*[1,1,1],'ForegroundColor',[0.1,0.1,0.9],...'string','0','fontsize',16,'fontname','timesnewroman');

uicontrol(gcf,'style','text',...

'unit','normalized','position',[0.22,0.41,0.08,0.04],...

17

姓名 王锐 班级 12级电子信息工程二班 学号 1208421042

'BackgroundColor',0.8*[1,1,1],'ForegroundColor',[0.1,0.1,0.9],...'string','100','fontsize',16,'fontname','timesnewroman')

在人机对话框中调节数据：=632.8nm,曲率半径R=8.64m，空气层厚度h=67nm。

平面光栅的衍射：

clc

clear all

d=-4*pi:0.001*pi:4*pi;

b=d./15;

p=1-(sinc(b).*sin(4*d)./sin(d)).^2;

%当要求P的曲线分布图时P=(sinc(b).*sin(4*d)./sin(d)).^2;

4是N值可调plot(d,P);

lgray=zeros(100,3);

for i=0:99

lgray(i+1,:)=(99-i)/99;

endimagesc(P)

colormap(lgray)

18