利用MATLAB语言进行光学衍射现象的仿真

利用MATLAB语言进行光学衍射现象的仿真

储林华（安庆师范学院物理与电气工程学院安徽安庆246011）

指导教师：张杰

摘要：光的衍射是光的波动性的一种重要表现，因此对光的衍射现象的研究，不仅具有重要的理论意义，而且在光学仪器研制和成像分析等诸多实际应用方面均有重要价值，但是其衍射光强的计算非常复杂，对实验条件的要求非常高，通常情况下很难得到满意的效果，严重影响了光学的教学。本文从衍射的相关理论知识出发，首先介绍了惠更斯--菲涅耳原理及其数学表示形式，然后重点讨论了单色光经各种对称光学衍射元件（单缝，双缝，光栅，圆孔）的夫琅和费衍射情况，并分别给出了它们在焦平面上的衍射光强计算公式，最后利用科学计算软件MA TLAB对光的衍射现象进行了仿真，所得到的图样细致逼真，使整个物理过程变得直观形象，且与实验所得到的衍射图样进行了比较，两者吻合得很好，从而为光学的理论分析和实验教学提供了一种新的途径。

关键词：光的衍射，光栅衍射，圆孔衍射，Matlab，计算机仿真

0 引言

光的衍射现象是光具有波动性的重要特征，因此对衍射现象的研究无论在理论上还是在实践中都有很重要的意义。对光的衍射现象的研究，始于17世纪，当时著名的荷兰科学家惠更斯提出了光是一种波的假说，并根据波动理论提出了光的传播理论——即惠更斯原理[1]，根据这一原理，他解释了光的反射定律和折射定律，给出了折射率的意义，光在两种介质中的速度比。到了19世纪，法国年轻的科学家菲涅耳，根据叠加原理把惠更斯原理进一步具体化，给出了光在传播过程中光强学计算公式，这就是著名的惠更斯-菲涅耳原理[2]。但由于在实际应用过程中，障碍物形状的不规则性，导致光强的计算公式几乎无解析解，只能进行一些数值计算。

针对衍射计算中出现的困难，近代的研究人员想到运用科学的计算软件MA TLAB，利用其较强的绘图和图象功能，编写计算程序，使得多种衍射元件（单缝，双缝，光栅，矩孔，圆孔）下的衍射现象得以在计算机中形象地被模拟仿真。这种做法，条件限制较少，对于衍射的实验教学是一种的补充，起到了很不错的效果；但值得指出的是，许多前人撰写的论文，总是在系统化，可视化，条理化方面不够理想，本文将在他们工作的基础上，将此课题进一步做得完美。

1．惠更斯—菲涅耳原理

早在十七世纪后期，荷兰科学家惠更斯就提出了光是一种波动的假说，并阐述了关于波面传播的一种理论，既惠更斯原理。该原理认为，传播中的波面上任何一点都可以认为是一个新的次波源，由于这些次波源发出的次波是球面波，这些次波的公共包络面就是下一时刻的波面，根据这一原理，他解释了光的反射定律和折射定律，并给出了折射率的意义，光在两种介质中的速度比。

菲涅耳根据叠加原理把惠更斯原理进一步具体化了，他假设各次波是球面波，但这些球面只是等相面而不是等幅面，球面上各点振幅与传播方向有关，这就避免了次波的后向的传播；同时，他认为，下一时刻空间任一点的振动由各次波到达该点的振动叠加所决定，此定理称为惠更斯—菲涅耳定理。如右图（1）所示。菲涅耳假设：Q 点所发的次波对P 点贡献dU(p)正比于Q 点附近一个小面元的面积ds ， 正比于Q 点的复振

幅U(Q)，正比于e jkr /r (假定次波是球面波)，以及正比于一个与传播方向有关的函数f(θ)(θ是r 与小面元法线即波在Q 点的传播方向上的夹角)，即：

()()()jkr

u Q f e du p c ds r

θ=

其中C 是一个与r ，Q ，θ无关的比例系数，f(θ)叫倾斜因子，它随θ增加而缓慢减少于是，按照叠加原理，有：

()()()()jkr

s

s

u Q f u p du p c

e ds r

θ==?? （1） 这就是惠更斯—菲涅耳原理的数学表达式，积分表示整个波面S 上各点所发的次波传播到P 点的作用的叠加。

从上面的表述，我们可以看到菲涅耳的思想比惠更斯有了更大的进步，他着眼于下一时刻空间各点的振动情况，而惠更斯只着眼于下一时刻波面的形状与位置，因此惠更斯只能定性地描述光的传播方向，而菲涅耳却能定量地描述衍射后的光强分布。

2. 夫朗禾费衍射现象研究

光的衍射现象根据光源到衍射屏以及观察屏的距离远近，可以分为近场衍射和远场衍射。如果光源到衍射屏以及观察屏的距离为有限值，则称为近场衍射（菲涅耳衍射）；如果光源到衍射屏以及观察屏的距离为无穷远，则称远场衍射，由于这种衍射最先由夫朗禾费在探索天体成像时作了系统的研究，故亦称为夫朗禾费衍射。一般情况下，利用（1）式进行光强计算时，菲涅耳衍射情况比较复杂，而夫朗禾费衍射情况比较简单,本文仅讨论后者。 2.1 圆孔的夫朗禾费衍射

在Oxy 平面上，有一以O 为圆心、R 为半径的圆孔，如图2所示。现用一束平行于Z 轴的光线照射,经圆孔衍射后到达焦平面上，于是按照（1）式，P 点的光强为：

()()()Q

jkr Q

Q

u Q f u p c

e ds r θ=

?

圆孔

200

()()Q

R

jkr Q

u Q f c

e d d r π

ρψ

θρρψ===

??

图1 波面上各点复振幅的传播

图2 圆孔衍射示意图

(a)

(b)

设从圆心O 到P 点的距离为r 0，则从圆孔中任一点Q 到P 点的距离r Q 为：

'00sin cos sin Q r r oo r θρ?θ=+=+

于是,由于程差0Q r r ?=-比0r 小得多，而与波长λ可以比拟，因此，可以把振幅中的Q r 以0r 代替而只考虑相位因子的变化，即：

02()

'

00

()Q R jk r r u p c

e

d d π

ρρψ-=?? （2）

其中：0

'

()()jkr u Q f c c

e r θ=为常数，将r Q 代入（2）式中得： 2'

cos sin '2100

2()()R jk J x u p c e d d c R x π

ρψθρρψπ??

==??????

其中：1J 称为一阶贝塞尔函数，2sin x R πθ

λ=，由此，可知P 点的光强可表达为：

2

2

'2

2

4

12()()()J x I p u p c R x π??

==????

（3） 根据(3)式可以画出焦平面的光强分布图样，如下图(3)所示：

(a) (b)

图3 圆孔衍射图样

其中，图（b ）是夫朗禾费衍射图样的照片。显然，这是一个圆对称图形，中心为主极大，在F 处，

0000,x I I θ====’R ，，，其他各点的光强可通过 'R 来表示，由图(2)可知：

''

'

sin 2xf R f R

λ

θπ== 当：

''

00.61 1.22,0,f R x I I R λπ===时， 为第一暗环；

''

0.8180.0175,f R R λπ==0时，x=1.635,I I 为第一亮环；

''

01.116 2.233,0,f R x I I R λπ===时， 为第二暗环；

''

01.619 3.238,0f R x I I R λπ===时，， 为第三暗环。

这就是夫朗禾费圆孔衍射的光强分布，它的中心永远是亮的，并且在中心亮斑处的光能占总光能的约84%，中心亮斑的半径也是确定的，为'f R λ0.61，或角半径为0.61R λ。 2.2 单缝的夫朗禾费衍射

图4 单缝夫朗禾费衍射示意图

设波长为λ的平面波射向缝宽AB=b 的狭缝，衍射经透镜L 会聚在焦平面F 上，取XZ 坐标如图(4)所示，根据惠—菲原理，在焦平面上任一点P 的复振幅为：

()()()jkr

u Q f u p c

e ds r

θ=

?

狭缝

把狭缝细分为垂直于X 轴的许多小面元，面积为ds=Ldx ，L 是缝的长度，由于是平面波入射，故U （Q ）为常数，在角度不大的情况下，0()1,sin f r r x θθ≈=+，由于0x r <<，故在振幅中的r 可视为0r ，只有相位因子中的0sin r r x θ?=-=不可略去，故：

2

sin 0

()()

b jkr jkx b e u p cLu Q e dx r θ-=?

把积分求出，得：

'sin sin '

2

2

2()()sin(sin )sin sin 2

b b jk jk

c b c

u p e

e

k jk k θθ

θθ

θ-=-=

令： sin sin 2

b

b k

πθθβλ== 则：'sin ()u p c b ββ=

2'22222

0()()sin sin I u p u p c b I ββββ*=== (4)

其中，0I 为0β=时的光强，即：衍射斑中心点的光强，β

的物理意义是狭缝的边缘与中心的光

线在P 点产生的相位差，以β或sin θ作横坐标，以0I I 作纵坐标，对(4)式作图可得夫朗禾费单缝衍射的光强分布，如下图所示：

（a ）单缝衍射的光强分布

-20

-15

-10

-5

5

10

15

20

00.20.40.60.81

(b)单缝衍射图样

图5 单缝衍射光强分布及图样

计算表明，在 πβπ-<<之间的主极大集中了90%的能量，主极大的半角宽为sin θθλ≈=，把狭缝改为矩孔，即在x ，y 两个方向上考虑衍射效应，则光强表达式为：

22220sin sin I I ββαα= 其中，α是在y 方向上与β对应的量，矩孔衍射图样如图(6)所示，

图6 矩孔衍射图样

2.3 双缝的夫朗禾费衍射

以上考虑的是一个缝的衍射，如果有两个相邻的缝，由于衍射，通过两个缝的光在观察屏上会相遇。试验结果告诉我们，在两条缝的衍射光相互交叠的区域，不是简单的呈现光强的叠加，而是出现了由光强重新分布而产生的明暗相间的条纹，称为光的干涉。

显然，光的干涉与光的衍射一样，都是波的叠加原理所必然导致的结果。

现在，我们仍然用惠—菲原理来分析两条缝所产生的夫朗禾费衍射的结果，如图(7)所示：

图7 双缝衍射示意图

如图(7)，衍射屏上A ，B 处各有一条宽为b 的缝，缝间距为d ， 由于透镜的作用，这两条缝的衍射光在焦平面F 上的光强分布是完全一样的，但它们的相位分布不同，把坐标原点分别放在A 与B 的

中心，由式： 2

'

sin 2()b jkx b u p c

e dx θ-=?

可得：

2

()sin '

2

2

()A b d

jk x A A b u p c e

dx θ

+-=?

2

()sin '

2

2

()B b d

jk x B B b u p c

e

dx θ

--=?

令：sin 2

d

r k

θ=为单缝中心与双缝中心的光在P 点产生的相位差，则：

2

'

sin 2

()()()()b jkx jr jr A B b u p u p u P C

e dx e e θ--=+=+?

='

sin ()jr jr cb e e ββ--+

于是，总光强为：

220()()sin ()()jr jr jr jr I u p u p I e e e e ββ*--==++

=22204sin cos I r β

β （5）

这就是夫朗禾费双缝衍射的光强分布的表达式，设d=5b ， 即：r=5β，将式(5)作图，可得下图所示的光强分布图样：

1234

图8 双缝衍射光强分布图样

显然，这是由2

2sin

ββ与24cos r 相乘而得的图样，在有双缝时，原来的单缝衍射图形不是平均

地增加到原来的两倍，而是在sin k d θλ=时，增加到原来的四倍，在sin (21)/2k d θλ=+处，降为零。

2.4 多缝夫朗禾费衍射

上面关于双缝衍射的讨论，可以进一步推广到更多的缝，如图(9)，有N 条等间距的缝，缝宽均为d ， 则相邻缝的对应程差为 sin d θ?=，相位差为：2sin /d σπθλ=，由(4)式知：

2000()()()()(1)()(1)jr jr jr jr j u p u p e e u p e e u p e σ---=+=+=+

图9 多缝夫朗禾费衍射示意图

如果把坐标原点放在第一个缝的中心，则0()jr u p e 就是它的单缝衍射振幅，而()0()j r u p e σ-则是另一个缝的衍射振幅，根据这个思想，立即可以把上式推广到多缝：

12()()()()N N u P u p u p u p =++

+

=2(1)1()(1)j j N j u p e

e e

σ

σ

σ

----++++=11()

1jN j e u p e σ

σ

-- 于是，光强分布为：

11

(1)(1)

()()()()()(1)(1)

jN jN N N N j j e e I P u p u p u p u p e e σσσσ

-*

*---==-- =22220sin sin (/2)sin (/2)I N β

βσσ （6）

设d=5b ， N=4， 按(5)式画出光强分布如下图所示：

图10 多缝衍射光强分布图

进而若N=100，则可得到的衍射图样为：

-20

-15

-10

-5

5

10

15

20

02000

40006000800010000

图11 光栅衍射图样

分析：多缝衍射与双缝衍射都是单缝衍射的结果。它们的区别在于多缝是干涉的结果，使极大值变细，峰值变高；而在两个极大值之间出现(N —1)个光强为零的β值和(N —2)个次峰。当N 很大时，各主极大变得十分尖锐，可以用来计量和分光，这种N 很大的能产生多缝衍射的光学元件称为光栅，d 叫做光栅常数。

3． 利用Matlab 软件进行衍射现象的仿真

3.1 Matlab 的简介

这里使用的软件MA TLAB6.5， 它是一个功能十分强大的应用软件，可以在很多学科中得到应用，它集数值计算，符号计算，数据可视化，系统动态仿真于一体，与其他计算机语言相比它更加灵活，更加接近科技人员的思维方式，因而编程效率更高。 3.2 仿真方法

下面我们应用Matlab 强大的图像处理与数值计算功能，仿真模拟光的夫朗禾费衍射（各种衍射屏）图样和强度分布曲线。仿真模拟首先是根据光的衍射光强分布的理论公式I(x ，y)及实验参数建立光强数据矩阵B(x ，y)，然后运用Matlab 的Image 命令绘制衍射图像，运用Plot 命令绘制光强分布曲线。 3.3 圆孔衍射的程序

% 圆孔衍射

% function ykys (r,f,lambda) r= 0.03; f=600;

lambda=600;

x= -20:0.05:20; y=-20:0.05:20;

[x1,y1]=meshgrid(x,y);

theta=atan(sqrt(x1.^2+y1.^2)./f);

x=1000000*(2*pi*r/lambda).*sin(theta+(theta==0)*eps);

I=(2*bessel(1,x)./x).^ 2;

subplot(2,1,1);

mesh(x1,y1,I);

axis([-20 20 -20 20 0 1])

subplot (2,1,2);

subimage (I*255)

axis off

圆孔衍射图样照片见图3（b）所示。

3.4 光栅衍射程序

% 光栅衍射

% function gsys (b，d，f，lambda，N)

b=0.05;

d=0.1;

f=600;

lambda=600;

N=1/2/100;

x=-20:0.1:20;

theta =atan (x./f);

beta=1000000*(pi*b/lambda).*sin(theta+(theta==0)*eps);

delta=beta.*(d/b);

I=((sin(beta)./beta).^2).*((sin(N.*delta)./sin(delta)).^2);

subplot(2,1,1);

plot(x,I);

subplot(2, 1 ,2);

colormap(gray(10))

image(I*255);

axis off

运行后得到的图样分别见图5（b），图8，图11所示。

3.5 矩孔衍射的程序

% 矩孔衍射

% function jkys (a,b,f,lambda)

a=0.05;

b=0.05;

f=600;

lambda=600;

x= -20:0.05:20;

y= -20:0.05:20;

[X,Y]=meshgrid(x,y);

theta1=atan (X./f);

theta2=atan(Y./f);

beta=1000000*(pi*b/lambda).*sin(theta1+(theta1==0)*eps);

alpha=1000000*(pi*a/lambda).*sin(theta2+(theta2==0)*eps);

I=((sin(beta)./beta).^2).*((sin(alpha)./alpha).^2);

subplot(2,1,1);

mesh(X,Y,I);

axis([ -20 20 -20 20 0 1])

subplot(2,1,2);

subimage(I*255);

运行后得到矩孔衍射图样见图6所示。

4．结束语：

经过Matlab 软件的仿真，得到各种形状孔屏的夫朗禾费衍射图样，这与实验室做实验得到的结果符合得很好；而且通过仿真，也更进一步使夫朗禾费衍射的相关知识形象化，可视化，加深了对这些内容的理解。

参考文献：

[1] 章志鸣，光学，北京:高等教育出版社,2000.

[2] 姚启钧，光学教程，北京:高等教育出版社,2002.

[3] 谢嘉宁，陈伟成，matlab在光学信息处理仿真实验中的应用，物理实验，V ol. 24,No. 6,2004:23~26.

[4] 谢嘉宁，赵建林，夫朗禾费衍射的计算机仿真，大学物理，V ol. 23,No. 3,2004:51~54

[5]楼顺天，于卫，阎华粱，matlab程序设计语言，西安：西安电子科技大学出版社，2000

[6] 孙向阳，圆孔衍射与象分辨本领的matlab仿真分析，湖北工业大学学报，Vol. 19，No. 5，2004：47~49.

[7] 王炳武，MATLAB5.3实际教程，北京：中国水利水电出版社，2000.

[8] 赵凯华，钟锡华，光学，北京大学出版社，1984.

[9] Grant R. Fowles. Introduction to Modern Optics. New York: Dover Publications, Inc. ,1989.

[10] Portmann S, Luthi H P. Molekel: An Interractive Molecular Graphics Tool. CHIMIA, 2000.

Simulation of Light Diffraction by MATLAB

Linhua Chu

(School of Physics and Electrical Engineering of Anqing Normal College，Anqing 246011) Abstact: The light’s diffraction is a very important behavior of the light’s f luctuation character. So the study

of the light’s diffraction phenomenon has not only the important theory value but also the important value on many practical using，such as optics instrument’s development and imagery analysis and so on. However，the calculation of the diffraction light intensity is very complicated and has a very high standard on the experiment condition. Usually, it’s not easy to get the pleased results, which affects the optics teaching seriously. This thesis is based on the related reason knowledge of diffraction. Firstly，it introduces Huygens-Fresnel principle and the math expression form. Then it emphatically discusses the single frequency’s light of Fraunhofer diffraction cases through many kinds of the symmetric element of the ligh t’s diffraction (single crack, double crack, grating, rectangular hole, round hole),furthermore, the thesis gives obtain the diffraction light intensity calculation on the focal plane. Finally，it uses the scientific calculation software MATLAB to simulate the light’s diffraction phenomenon. The pattern we got is delicate and realistic，which makes the whole physical process intuitional and visual. And compared it with the diffraction pattern which is be gotten from experiment, they coincide with each other well，which supply a new way for the optics’ theory analysis and experiment teaching.

Key words:light’s diffraction，grating diffraction，round hole diffraction，Matlab ，computer simulation