matlab波动光学的模拟

用matlab进行波动光学的计算机模拟

学生高志强指导教师常山

物理与电子信息系03级1班江西上饶334001

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

关键词:Matlab;波动光学;程序设计;计算机模拟

1引言

1984年由美国的MathWorks公司推出了MATLAB仿真软件，很快成为应用学科计算机辅助分析、设计、仿真和教学不可缺少的软件，应用到生物医学工程、信号分析、语音处理、图像识别、航天航海工程、统计分析、计算机技术、控制和数学等领域中[1]。

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

2MATLAB的使用[1][2][3]

MATLAB有两种运行方式：命令行方式（直接在命令窗口中输入命令行来实现计算或作图功能）和M文件方式（将编写的程序保存成文本文件，要执行时，在Debug菜单中选择Run选项，将在命令窗口中输出程序的执行结果）。

2.1命令行方式

直接在命令窗口输入命令行来实现计算或作图功能。例如，输入矩阵A和B的和，其中

1223A=，B=4534

首先打开MATLAB界面，然后在命令窗口输入下面命令行

A=[12;3

B=[24];5];3;4

作者简介:高志强(1984-),男,江西婺源人,本科生,喜爱光电技术和程序设计。

指导教师:常山(1964-),男,山东阳谷人,硕士,讲师,主要从事光电精密测量和光学信息处理。

C=A+B

其运算结果为

2.2M文件形式

当然和其它程序语言一样，Matlab程序也可以用一般的文字编程器来编写，如MS-WORD、PE等都可以，但必须要将编写的程序保存成文本文件。此外，在MATLAB的命令行上也可以输入EDIT命令来编辑文件，例如：edit或editdddd编辑之后，文件名会被自动保存为dddd.m。也可以用鼠标在MATLAB界面中点击File菜单中选择New中的M-file。

M文件有两个形式：一种称为命令文件（scriptfile）就好像dos下的批处理文件一样。这类程序包含了一连串的MATLAB命令，执行时依次执行。另一种称为函数文件（functionfile）。

它的第一句可执行语句是以function引导的定义语句，在函数文件中的变量都是局部变量。

2.2.1MATLAB的命令文件

命令文件实际上只包括两部分，即注解和指令。注解部分开头必须用百分比符号％注明。命令文件中的语句，以访问MATLAB工作间（workspace）中的所有数据。运行过程中的所有变量均是全局变量，这些变量一旦生成，就一直保存在内存空间中，除非用户用clear指令将它们清除。

运行一个命令文件，等价于从命令窗（commandwindow）中按顺序连续运行文件里的指令。由于命令文件只是一串指令的集合，因此，程序不需要预先定义，只需要按指令窗口中的指令输入顺序，依次地将指令编辑在命令文件中就可以了[3]。

例如：画出三维曲面法线图

运用MATLAB编辑器建立.M文件，然后输入指令：

[x,y,z]=cylinder(1:10);

Surfnorm(x,y,z)

图形结果如图

2-1所示。

图2-1三维曲面法线图

2.2.2MATLAB的函数文件

如果M文件的第一行包含function，此文件就是函数文件。每一个函数文件都定义一个函数。事实上，MATLAB提供的函数指令大部分都是由函数文件定义的。这足以说明函数文件的重要性。从使用角度看，函数是一个“黑箱”，把一些数据送进去，经加工处理，把结果送出来。从形式上看，函数文件区别于命令文件之处是：命令文件的变量在文件执行完后保留在内存中；而函数文件内定义的变量仅在函数文件内部起作用，当函数文件执行完之后，这些内部变量将被消除。

MATLAB函数文件实际上包含五个部分：

(1)函数定义行；

(2)函数主体；

(3)H1行；

(4)函数说明；

(5)注解。

其中，完整的函数一定要包括函数定义行与函数主体，其余三个部分都是为了辅助使用的。例：下面即为一个函数的例子，它接受一个参数，并且执行comp(x)=(x+10)^2*10运算，然后返回一个值。在函数定义行上，除了function为关键字外，返回值就直接定义成y=comp(x)。

先将函数文件存放于磁盘中，并起名为exp42.m。如果要执行这个函数，可以在命令行上直接输入:

exp42(10)

则会得到结果：

ans=

4000

或输入：

k=10;

exp42(k)

得到结果：

ans=

4000

3波动光学基本理论[4][5]

3.1光的干涉

干涉是指因波的叠加而引起强度相长或相消的现象。有分波面、分振幅和分振动面三种干涉类型（方法）。

3.1.1双光束干涉

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

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

屏

x

图3-1双缝干涉

明、暗纹的位置由两束光的光程差Δ决定：

k0,1,2,明纹xkdsindr02k1k0,1,2,暗纹2

r屏上条纹间距：x0d

其干涉条纹的特点条纹形状为一组与狭缝平行、等间隔的直线。

屏上光强分布为II1I22I1I2cosdr0（3-1）（3-2）x（3-3）

注：而较为严格的表达式为后面（3-16）式。模拟干涉图样见图4-3和4-4。菲涅耳双棱镜、菲涅耳双面镜和劳埃镜等的干涉，都与此类似。

3.1.2牛顿环

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

在图3-2中，设透镜的曲率半径为R，对应于某干涉环的空气层厚度为h，该环的半径为r。在光线正入射时，由

2hn2n0sin2i2

2（3-4）

可得从空气层两表面反射的相干光线的光程差为2h，式中前为正号，是因为相位2

突变发生空气膜的下表面，由ΔDEC可得R2r2(Rh)2再化简得

r22hRh2h(2Rh)

由于2Rh，上式可以继续化简为（3-5）

r22hR

得到（3-6）

r2

（3-7）

当k时，由上式可得第K亮环的半径

rkR(k=1,2,3，……)（3-8）当k时，可得第k暗环的半径

rkR(k=1，2，3，……)（3-9）

当k=0时，r=0，与环中心对应，即牛顿环中心为一暗点，它来源于空气层下表面反射光的相位跃变。=568nm、曲率半径R=8.01m、空气层厚度h=42nm时的干涉图样如图

3-3所示。

图3-3牛顿环模拟图

图3-4牛顿环的照片图

对比图3-3和3-4，可见，利用MATLAB软件得出的模拟图跟照片的干涉图样是一致的，且比照片更为清晰。

3.2光的衍射

凡是不能用反射或折射予以解释的光偏离直线传播的现象，称为光的衍射。通常所研究的是具有各种形状的障碍物的衍射，将这些障碍物称作屏。入射波前在什么方向受到限制，衍射图样就沿该方向扩展，受限制愈严重，在该方向上的扩展也愈厉害。

按光源、衍射屏和接收屏三者之间的相对位置，可将衍射现象分为两种类型：（1）光源和接收屏或二者之一距离衍射屏为有限远时，所观察到的衍射称为菲涅耳衍射。菲涅耳衍射图样是带有衍射条纹的衍射孔的投影像。（2）光源和接收屏距离都在无穷远或相当于无穷远，在衍射孔上的入射波及其在各个方向的衍射波都可看成平面波，所观察到的衍射称为夫琅禾费衍射。夫琅禾费衍射的图样是带衍射条纹的光源的投影像，与衍射孔的形状很少相似或者完全不相似。

3.2.1基尔霍夫衍射积分公式

单色点光源L发出的球面波照射具有开孔S0的衍射屏后，衍射场中任一点P的光振幅可表示为

1~E（p）=ieikrcos0cos~E(Q)dsr2s0（3-10）式（3-10）即为菲涅耳——基尔霍夫衍射出积分公式。式中0和分别为L到Q的矢径和Q点到P点的矢径与Q点处面元ds的法线之间的夹角（如图3-5所示）。

图3-5点光源L对衍射场中P点的作用

3.2.2单狭缝的夫琅禾费衍射

在图3-5中，取坐标系OXYZ的原点在狭缝中心，Z轴沿系统的光轴，Y轴沿狭缝的长方向，X轴沿狭缝的宽方向。一般说来，当狭缝的长度b较宽度a大得多的情况下，衍射图样基本上只沿与狭缝垂直的方向扩展。因此，计算接收屏上的夫琅禾费衍射场时，只需考虑X轴方向的光振动的复振幅分布，即将它们作为一维问题处理。将入射波前经BB缝露出部分分为许多平行于缝长方向的等宽窄条带ds，各带发出的衍射角为方向的次波经L2会聚于接收屏上的P

点，因而L2像方焦面上的一点是与一个衍射方向相对应的。在缝内距原点为X的Q处，取一宽为dx的窄条带作为积分面元，即ds=bdx，它到P点的光程为r，按照菲涅耳——基尔霍夫公式，P点的光场为式（3-10）。

图3-6计算P点的复振幅用图

在傍轴条件下，cos0≈cos≈1,而且被积函数分母中的r可用缝中心至P点的光程r0代替，~~~E(Q)E0~t(x),E0为入射平面波在BB上有复振幅，它是一个复常数。~t(x)为单缝的复振幅透射率函数，它可以写为：

1xt(x)=0x

于是可得（3-11）

~bE0~E(p)ir0t(x)eikrdx（3-12）

又r0rsin，于是（3-12）式又可写为

~bE0eikr0~E(p)ir0_t(x)eikxsindx（3-13）

~表示，则得将(3-13)式中与x无关的量归并到一起用复常数c

sinkasin~~aE(p)c（3-14）ksinkasinasin若令，它表示狭缝边缘两点在P点所产生的振动相位差之半，则2

得：E(p)A0sin~~

~，P点的振幅为A(p)A，式中A0ac0~sin，P点的光强度则为

（3-15）2~~I(p)E(p)E(p)I0sin2

~式中I0A02~2表示屏中心处的光强，ac（3-15）式即单狭缝夫琅禾费衍射光强分布公式。

用Matlab编程进行模拟，得到的强度分布曲线和模拟现象如图3-7和3-8所示。

图3-7单缝衍射强度分布图

图3-8单缝衍射模拟图

3.2.3多缝干涉

单色光通过N个窄缝S1,S2SN射向屏幕，相当于位置不同的N个同频率同相位的光源向屏幕照射的叠加，由于到达屏幕各点的距离不同引起相位差，叠加的结果是：有的点加强、有的点抵消的现象而出现干涉条纹。表达式是

2II0sinn2（3-16）

多缝干涉的光强分布曲线图和干涉现象模拟图见第四章的图4-5。

3.2.4平面光栅的衍射

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

（1）与单缝衍射图样相比，多缝衍射的图样中出现一系列新的强度最大值和最小值，其中

那些较强的亮线称为主最大，较弱的亮线叫做次最大。

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

（3）相邻主最大之间有N-1条暗纹和N-2个次最大。（注：具体情形如图4-6的a、b所示）

（4）强度分布中保留了单缝衍射的因子，那就是曲线的包迹与单缝衍射强度曲线一样。光栅衍射条纹，理应是单缝衍射和缝间干涉的共同结果。

设光栅有N条狭缝，透镜焦距为f'，理论分析可以得到，夫琅禾费衍射的光屏上任意一点2P的光强为

Ipsin02sinNsin2（3-17）

式中dsin，上式的前一部分与单缝衍射的强度相同，它来源于单缝衍射，是整个衍射图样的轮廓，称为单缝衍射因子。后一部分分数可改写如下：因为dsin为相邻两缝对应点到达观察点的光程差（如图所示），这个光程差所引起的相位差为2dsin

（3-17）式的分数部分来源于缝间干涉，叫做缝间干涉因子，可写为

2sin(N)sinN2sinsin2()22。

可见，光栅衍射的光强是单缝衍射因子和缝间干涉因子的乘积，即是单缝衍射因子对干涉主最大起调制作用。

对于一定的波长来说，根据dsinj(j0,1,2,)，可知各级谱线之间的距离由光栅常量d决定，各级谱线的强度分布将随b与d的比值而改变。由

dsin)sinApA0sin(sin)sinN(

幅为（3-18）不难决定各主最大值的振幅分布，把由光栅方程决定的sin的值代入上式，得到j级谱线的振

AjA0NdsinbjbjA0Ndbsin(jjbd（3-19）

又因为单缝衍射最小值位置sink()(k1,2,)。由此可见，若k=1，2，…为衍射最小值的级数，j为主最大值的级数且jk，则当dj时，bksin(j)sink0，因而Aj0，故jkd时，级数为j的谱线消失，也就是缺

级。

关于上面论述的知识点，我们所学的教材[4]只给出了平面光栅衍射的强度分布曲线，仅由光栅衍射此强度函数曲线图，读者会还是会感觉抽象，难以对此现象形成更为直接的直观印象，特别是平面光栅衍射中涉及的缺级现象更是难以理解（课本中只给了以表格形式表示的光强分布公式）。本文下一章利用Matlab编写相应的程序，对此现象进行了模拟。由图4-7a和4-7b（d=4b，N=6）的情况可知：j=4k（k=1,2,…）时对应的亮条纹消失，也就是缺级。

3.2.5圆孔的夫琅和费衍射

将图3-6中的单狭缝换成圆孔，就可在接收屏上观察到圆孔的夫琅禾费衍射图样，它的中央是一个很亮的圆斑，外面分布着几圈很能淡的亮环。设圆孔的半径为R，与P点对应的衍射角为，令2Rsin

2J102表示圆孔边缘与中心的次边缘在P点所产生的振动的相位差，根据菲涅耳——基尔霍夫衍射积分公式可求得P点的光强为IpI（3-20）式中I0为图样中心P0点的光强，J1为一级贝塞耳函数。衍射图样如图4-8所示。4波动光学的Matlab模拟

将以上几种波动现象的模拟所需要的程序编写如下（其中P代表光强）。

4.1单缝衍射clc

clearall

a=-2*pi:0.0001*pi:2*pi;P=(1-sinc(a)).^2;%

当要求P的曲线分布图时P=sinc(a).^2

plot(a,P)

lgray=zeros(256,3);

fori=0:255

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

end

imagesc(P)

colormap(lgray)

4.2多缝干涉

4.2.1双缝干涉当式（3-16）的n=2时即我们所熟知的杨氏双缝干涉。

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);

fori=0:255

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

end

imagesc(P)

colormap(lgray)

4.2.2多缝干涉当式（3-16）的n>2时，即多缝干涉,现在取n=4

clearall

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

P=1-sin(4*a).^2./sin(a).^2;%当要求P

的曲线分布图时P=sin(4*a).^2./sin(a).^2

plot(a,P)

lgray=zeros(256,3);

fori=0:255

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

end

imagesc(P)

colormap(lgray)

4.3平面光栅衍射clc

clearall

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);

fori=0:99

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

end

imagesc(P)

colormap(lgray)

上饶师范学院

优秀本科毕业论文图4-6多缝衍射（N=6，d=15b）模拟图

b

在图a和b中可见，主最大的位置与缝数N无关，但它们的宽度随N的增大而减小，强度正比于N，同时也可见相邻主最大之间有N-1条暗纹和N-2个次最大。2

图4-7d=4b,N=6多缝衍射模拟图，由图可知j=4时是缺级

由于屏幕分辨率小限制，故将图4-7的灰度图中的次最大条纹的亮度调得较大和其宽度调得较宽。

4.4圆孔夫琅禾费衍射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;

form=1:N;

forn=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,[])

图4-8圆孔衍射模拟图

4.5牛顿环

牛顿环演示的MATLAB程序

closeall;

figure('Position',[[1**********]]);

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;

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',...

'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','12','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],...

'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，牛顿环的MATLAB模拟结果如图4-9所示：

图4-9牛顿环模拟图

在电子文档中，当点击按钮改变波长，曲率半径，空气厚度时，图案将出现相应的变换：①改变波长时，条纹的颜色会变化；

②改变曲率半径时，条纹的形状会发生变化；

③改变空气层厚度使其增大时有条纹陷入，而使空气层厚度减小时会有条纹涌出和实验现象一致且符合理论结果。

参考文献

[1]闻新.MATLAB科学图形构建基础应用(6.X)(第一版)[M].北京:科学出版社,2002。

[2]胡守信.基于MATLAB的数学实验(第一版)[M].北京:科学出版社,2004.

[3]陈扬.MATLAB6.X图形编程与图像处理(第一版)[M].西安:西安电子科技大学出版社,2002.

[4]姚启钧.光学教程(第三版)[M].北京:高等教育出版社,2002.

[5]郭永康.光学(第一版)[M].北京:高等教育出版社,2005.

SimulatingWaveOpticswithMatlabonComputerWriter:GaoZhi-qiang,Tutor:ChangShan

(ShangraoNormalCollege,ShangraoJiangxi334001,China)

Abstract:Inteachingandresearch,computersimulationisincreasinglybeingtakenseriously.AfterthetheoryofinterferenceofFraunhoferdiffraction,two-beam,multi-slotandNewtonringsisanalyzed,thenmatlabisusedtocompilethecorrespondingproceduresandtosimulatebycomputer.Thisisusefultounderstandandstudyofthediffractionandinterferencetheory.Keywords:matlab;waveoptics;programming;

simulatingwithcomputer

用matlab进行波动光学的计算机模拟

学生高志强指导教师常山

物理与电子信息系03级1班江西上饶334001

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

关键词:Matlab;波动光学;程序设计;计算机模拟

1引言

1984年由美国的MathWorks公司推出了MATLAB仿真软件，很快成为应用学科计算机辅助分析、设计、仿真和教学不可缺少的软件，应用到生物医学工程、信号分析、语音处理、图像识别、航天航海工程、统计分析、计算机技术、控制和数学等领域中[1]。

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

2MATLAB的使用[1][2][3]

MATLAB有两种运行方式：命令行方式（直接在命令窗口中输入命令行来实现计算或作图功能）和M文件方式（将编写的程序保存成文本文件，要执行时，在Debug菜单中选择Run选项，将在命令窗口中输出程序的执行结果）。

2.1命令行方式

直接在命令窗口输入命令行来实现计算或作图功能。例如，输入矩阵A和B的和，其中

1223A=，B=4534

首先打开MATLAB界面，然后在命令窗口输入下面命令行

A=[12;3

B=[24];5];3;4

作者简介:高志强(1984-),男,江西婺源人,本科生,喜爱光电技术和程序设计。

指导教师:常山(1964-),男,山东阳谷人,硕士,讲师,主要从事光电精密测量和光学信息处理。

C=A+B

其运算结果为

2.2M文件形式

当然和其它程序语言一样，Matlab程序也可以用一般的文字编程器来编写，如MS-WORD、PE等都可以，但必须要将编写的程序保存成文本文件。此外，在MATLAB的命令行上也可以输入EDIT命令来编辑文件，例如：edit或editdddd编辑之后，文件名会被自动保存为dddd.m。也可以用鼠标在MATLAB界面中点击File菜单中选择New中的M-file。

M文件有两个形式：一种称为命令文件（scriptfile）就好像dos下的批处理文件一样。这类程序包含了一连串的MATLAB命令，执行时依次执行。另一种称为函数文件（functionfile）。

它的第一句可执行语句是以function引导的定义语句，在函数文件中的变量都是局部变量。

2.2.1MATLAB的命令文件

命令文件实际上只包括两部分，即注解和指令。注解部分开头必须用百分比符号％注明。命令文件中的语句，以访问MATLAB工作间（workspace）中的所有数据。运行过程中的所有变量均是全局变量，这些变量一旦生成，就一直保存在内存空间中，除非用户用clear指令将它们清除。

运行一个命令文件，等价于从命令窗（commandwindow）中按顺序连续运行文件里的指令。由于命令文件只是一串指令的集合，因此，程序不需要预先定义，只需要按指令窗口中的指令输入顺序，依次地将指令编辑在命令文件中就可以了[3]。

例如：画出三维曲面法线图

运用MATLAB编辑器建立.M文件，然后输入指令：

[x,y,z]=cylinder(1:10);

Surfnorm(x,y,z)

图形结果如图

2-1所示。

图2-1三维曲面法线图

2.2.2MATLAB的函数文件

如果M文件的第一行包含function，此文件就是函数文件。每一个函数文件都定义一个函数。事实上，MATLAB提供的函数指令大部分都是由函数文件定义的。这足以说明函数文件的重要性。从使用角度看，函数是一个“黑箱”，把一些数据送进去，经加工处理，把结果送出来。从形式上看，函数文件区别于命令文件之处是：命令文件的变量在文件执行完后保留在内存中；而函数文件内定义的变量仅在函数文件内部起作用，当函数文件执行完之后，这些内部变量将被消除。

MATLAB函数文件实际上包含五个部分：

(1)函数定义行；

(2)函数主体；

(3)H1行；

(4)函数说明；

(5)注解。

其中，完整的函数一定要包括函数定义行与函数主体，其余三个部分都是为了辅助使用的。例：下面即为一个函数的例子，它接受一个参数，并且执行comp(x)=(x+10)^2*10运算，然后返回一个值。在函数定义行上，除了function为关键字外，返回值就直接定义成y=comp(x)。

先将函数文件存放于磁盘中，并起名为exp42.m。如果要执行这个函数，可以在命令行上直接输入:

exp42(10)

则会得到结果：

ans=

4000

或输入：

k=10;

exp42(k)

得到结果：

ans=

4000

3波动光学基本理论[4][5]

3.1光的干涉

干涉是指因波的叠加而引起强度相长或相消的现象。有分波面、分振幅和分振动面三种干涉类型（方法）。

3.1.1双光束干涉

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

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

屏

x

图3-1双缝干涉

明、暗纹的位置由两束光的光程差Δ决定：

k0,1,2,明纹xkdsindr02k1k0,1,2,暗纹2

r屏上条纹间距：x0d

其干涉条纹的特点条纹形状为一组与狭缝平行、等间隔的直线。

屏上光强分布为II1I22I1I2cosdr0（3-1）（3-2）x（3-3）

注：而较为严格的表达式为后面（3-16）式。模拟干涉图样见图4-3和4-4。菲涅耳双棱镜、菲涅耳双面镜和劳埃镜等的干涉，都与此类似。

3.1.2牛顿环

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

在图3-2中，设透镜的曲率半径为R，对应于某干涉环的空气层厚度为h，该环的半径为r。在光线正入射时，由

2hn2n0sin2i2

2（3-4）

可得从空气层两表面反射的相干光线的光程差为2h，式中前为正号，是因为相位2

突变发生空气膜的下表面，由ΔDEC可得R2r2(Rh)2再化简得

r22hRh2h(2Rh)

由于2Rh，上式可以继续化简为（3-5）

r22hR

得到（3-6）

r2

（3-7）

当k时，由上式可得第K亮环的半径

rkR(k=1,2,3，……)（3-8）当k时，可得第k暗环的半径

rkR(k=1，2，3，……)（3-9）

当k=0时，r=0，与环中心对应，即牛顿环中心为一暗点，它来源于空气层下表面反射光的相位跃变。=568nm、曲率半径R=8.01m、空气层厚度h=42nm时的干涉图样如图

3-3所示。

图3-3牛顿环模拟图

图3-4牛顿环的照片图

对比图3-3和3-4，可见，利用MATLAB软件得出的模拟图跟照片的干涉图样是一致的，且比照片更为清晰。

3.2光的衍射

凡是不能用反射或折射予以解释的光偏离直线传播的现象，称为光的衍射。通常所研究的是具有各种形状的障碍物的衍射，将这些障碍物称作屏。入射波前在什么方向受到限制，衍射图样就沿该方向扩展，受限制愈严重，在该方向上的扩展也愈厉害。

按光源、衍射屏和接收屏三者之间的相对位置，可将衍射现象分为两种类型：（1）光源和接收屏或二者之一距离衍射屏为有限远时，所观察到的衍射称为菲涅耳衍射。菲涅耳衍射图样是带有衍射条纹的衍射孔的投影像。（2）光源和接收屏距离都在无穷远或相当于无穷远，在衍射孔上的入射波及其在各个方向的衍射波都可看成平面波，所观察到的衍射称为夫琅禾费衍射。夫琅禾费衍射的图样是带衍射条纹的光源的投影像，与衍射孔的形状很少相似或者完全不相似。

3.2.1基尔霍夫衍射积分公式

单色点光源L发出的球面波照射具有开孔S0的衍射屏后，衍射场中任一点P的光振幅可表示为

1~E（p）=ieikrcos0cos~E(Q)dsr2s0（3-10）式（3-10）即为菲涅耳——基尔霍夫衍射出积分公式。式中0和分别为L到Q的矢径和Q点到P点的矢径与Q点处面元ds的法线之间的夹角（如图3-5所示）。

图3-5点光源L对衍射场中P点的作用

3.2.2单狭缝的夫琅禾费衍射

在图3-5中，取坐标系OXYZ的原点在狭缝中心，Z轴沿系统的光轴，Y轴沿狭缝的长方向，X轴沿狭缝的宽方向。一般说来，当狭缝的长度b较宽度a大得多的情况下，衍射图样基本上只沿与狭缝垂直的方向扩展。因此，计算接收屏上的夫琅禾费衍射场时，只需考虑X轴方向的光振动的复振幅分布，即将它们作为一维问题处理。将入射波前经BB缝露出部分分为许多平行于缝长方向的等宽窄条带ds，各带发出的衍射角为方向的次波经L2会聚于接收屏上的P

点，因而L2像方焦面上的一点是与一个衍射方向相对应的。在缝内距原点为X的Q处，取一宽为dx的窄条带作为积分面元，即ds=bdx，它到P点的光程为r，按照菲涅耳——基尔霍夫公式，P点的光场为式（3-10）。

图3-6计算P点的复振幅用图

在傍轴条件下，cos0≈cos≈1,而且被积函数分母中的r可用缝中心至P点的光程r0代替，~~~E(Q)E0~t(x),E0为入射平面波在BB上有复振幅，它是一个复常数。~t(x)为单缝的复振幅透射率函数，它可以写为：

1xt(x)=0x

于是可得（3-11）

~bE0~E(p)ir0t(x)eikrdx（3-12）

又r0rsin，于是（3-12）式又可写为

~bE0eikr0~E(p)ir0_t(x)eikxsindx（3-13）

~表示，则得将(3-13)式中与x无关的量归并到一起用复常数c

sinkasin~~aE(p)c（3-14）ksinkasinasin若令，它表示狭缝边缘两点在P点所产生的振动相位差之半，则2

得：E(p)A0sin~~

~，P点的振幅为A(p)A，式中A0ac0~sin，P点的光强度则为

（3-15）2~~I(p)E(p)E(p)I0sin2

~式中I0A02~2表示屏中心处的光强，ac（3-15）式即单狭缝夫琅禾费衍射光强分布公式。

用Matlab编程进行模拟，得到的强度分布曲线和模拟现象如图3-7和3-8所示。

图3-7单缝衍射强度分布图

图3-8单缝衍射模拟图

3.2.3多缝干涉

单色光通过N个窄缝S1,S2SN射向屏幕，相当于位置不同的N个同频率同相位的光源向屏幕照射的叠加，由于到达屏幕各点的距离不同引起相位差，叠加的结果是：有的点加强、有的点抵消的现象而出现干涉条纹。表达式是

2II0sinn2（3-16）

多缝干涉的光强分布曲线图和干涉现象模拟图见第四章的图4-5。

3.2.4平面光栅的衍射

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

（1）与单缝衍射图样相比，多缝衍射的图样中出现一系列新的强度最大值和最小值，其中

那些较强的亮线称为主最大，较弱的亮线叫做次最大。

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

（3）相邻主最大之间有N-1条暗纹和N-2个次最大。（注：具体情形如图4-6的a、b所示）

（4）强度分布中保留了单缝衍射的因子，那就是曲线的包迹与单缝衍射强度曲线一样。光栅衍射条纹，理应是单缝衍射和缝间干涉的共同结果。

设光栅有N条狭缝，透镜焦距为f'，理论分析可以得到，夫琅禾费衍射的光屏上任意一点2P的光强为

Ipsin02sinNsin2（3-17）

式中dsin，上式的前一部分与单缝衍射的强度相同，它来源于单缝衍射，是整个衍射图样的轮廓，称为单缝衍射因子。后一部分分数可改写如下：因为dsin为相邻两缝对应点到达观察点的光程差（如图所示），这个光程差所引起的相位差为2dsin

（3-17）式的分数部分来源于缝间干涉，叫做缝间干涉因子，可写为

2sin(N)sinN2sinsin2()22。

可见，光栅衍射的光强是单缝衍射因子和缝间干涉因子的乘积，即是单缝衍射因子对干涉主最大起调制作用。

对于一定的波长来说，根据dsinj(j0,1,2,)，可知各级谱线之间的距离由光栅常量d决定，各级谱线的强度分布将随b与d的比值而改变。由

dsin)sinApA0sin(sin)sinN(

幅为（3-18）不难决定各主最大值的振幅分布，把由光栅方程决定的sin的值代入上式，得到j级谱线的振

AjA0NdsinbjbjA0Ndbsin(jjbd（3-19）

又因为单缝衍射最小值位置sink()(k1,2,)。由此可见，若k=1，2，…为衍射最小值的级数，j为主最大值的级数且jk，则当dj时，bksin(j)sink0，因而Aj0，故jkd时，级数为j的谱线消失，也就是缺

级。

关于上面论述的知识点，我们所学的教材[4]只给出了平面光栅衍射的强度分布曲线，仅由光栅衍射此强度函数曲线图，读者会还是会感觉抽象，难以对此现象形成更为直接的直观印象，特别是平面光栅衍射中涉及的缺级现象更是难以理解（课本中只给了以表格形式表示的光强分布公式）。本文下一章利用Matlab编写相应的程序，对此现象进行了模拟。由图4-7a和4-7b（d=4b，N=6）的情况可知：j=4k（k=1,2,…）时对应的亮条纹消失，也就是缺级。

3.2.5圆孔的夫琅和费衍射

将图3-6中的单狭缝换成圆孔，就可在接收屏上观察到圆孔的夫琅禾费衍射图样，它的中央是一个很亮的圆斑，外面分布着几圈很能淡的亮环。设圆孔的半径为R，与P点对应的衍射角为，令2Rsin

2J102表示圆孔边缘与中心的次边缘在P点所产生的振动的相位差，根据菲涅耳——基尔霍夫衍射积分公式可求得P点的光强为IpI（3-20）式中I0为图样中心P0点的光强，J1为一级贝塞耳函数。衍射图样如图4-8所示。4波动光学的Matlab模拟

将以上几种波动现象的模拟所需要的程序编写如下（其中P代表光强）。

4.1单缝衍射clc

clearall

a=-2*pi:0.0001*pi:2*pi;P=(1-sinc(a)).^2;%

当要求P的曲线分布图时P=sinc(a).^2

plot(a,P)

lgray=zeros(256,3);

fori=0:255

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

end

imagesc(P)

colormap(lgray)

4.2多缝干涉

4.2.1双缝干涉当式（3-16）的n=2时即我们所熟知的杨氏双缝干涉。

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);

fori=0:255

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

end

imagesc(P)

colormap(lgray)

4.2.2多缝干涉当式（3-16）的n>2时，即多缝干涉,现在取n=4

clearall

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

P=1-sin(4*a).^2./sin(a).^2;%当要求P

的曲线分布图时P=sin(4*a).^2./sin(a).^2

plot(a,P)

lgray=zeros(256,3);

fori=0:255

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

end

imagesc(P)

colormap(lgray)

4.3平面光栅衍射clc

clearall

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);

fori=0:99

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

end

imagesc(P)

colormap(lgray)

上饶师范学院

优秀本科毕业论文图4-6多缝衍射（N=6，d=15b）模拟图

b

在图a和b中可见，主最大的位置与缝数N无关，但它们的宽度随N的增大而减小，强度正比于N，同时也可见相邻主最大之间有N-1条暗纹和N-2个次最大。2

图4-7d=4b,N=6多缝衍射模拟图，由图可知j=4时是缺级

由于屏幕分辨率小限制，故将图4-7的灰度图中的次最大条纹的亮度调得较大和其宽度调得较宽。

4.4圆孔夫琅禾费衍射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;

form=1:N;

forn=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,[])

图4-8圆孔衍射模拟图

4.5牛顿环

牛顿环演示的MATLAB程序

closeall;

figure('Position',[[1**********]]);

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;

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',...

'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','12','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],...

'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，牛顿环的MATLAB模拟结果如图4-9所示：

图4-9牛顿环模拟图

在电子文档中，当点击按钮改变波长，曲率半径，空气厚度时，图案将出现相应的变换：①改变波长时，条纹的颜色会变化；

②改变曲率半径时，条纹的形状会发生变化；

③改变空气层厚度使其增大时有条纹陷入，而使空气层厚度减小时会有条纹涌出和实验现象一致且符合理论结果。

参考文献

[1]闻新.MATLAB科学图形构建基础应用(6.X)(第一版)[M].北京:科学出版社,2002。

[2]胡守信.基于MATLAB的数学实验(第一版)[M].北京:科学出版社,2004.

[3]陈扬.MATLAB6.X图形编程与图像处理(第一版)[M].西安:西安电子科技大学出版社,2002.

[4]姚启钧.光学教程(第三版)[M].北京:高等教育出版社,2002.

[5]郭永康.光学(第一版)[M].北京:高等教育出版社,2005.

SimulatingWaveOpticswithMatlabonComputerWriter:GaoZhi-qiang,Tutor:ChangShan

(ShangraoNormalCollege,ShangraoJiangxi334001,China)

Abstract:Inteachingandresearch,computersimulationisincreasinglybeingtakenseriously.AfterthetheoryofinterferenceofFraunhoferdiffraction,two-beam,multi-slotandNewtonringsisanalyzed,thenmatlabisusedtocompilethecorrespondingproceduresandtosimulatebycomputer.Thisisusefultounderstandandstudyofthediffractionandinterferencetheory.Keywords:matlab;waveoptics;programming;

simulatingwithcomputer