计算机模拟光学衍射

1、 长江大学物理科学与技术学院题目： 计算机模拟光学衍射系 别 应用物理 专 业 光电子技术 班 级 应用物理 1001 姓 名 王飞 学 号 202106927 序 号 22 老 师 陈海燕 2021 年 11 月 28 日课题背景及研究意义衍射是光波的一种重要特征,对光波衍射现象的讨论是以惠更斯- 菲涅尔原理为根底的。在大多数光学问题中，可以采用光场的标量衍射理论来描述和计算。光场的频率非常高，只能测量在一个比光的周期大得多的时间间隔内的平均值，在计算衍射光波场的光振动和光强分布时那么利用菲涅尔- 基尔霍夫衍射积分公式。在物理光学的光的衍射教学过程中， 如不借助实验, 学生很难理解理论、实验

2、原理和过程。而光学实验的进行一般要稳定的环境、高精密的仪器以及光路的调试, 花费相当多的时间，并且很难调试到理想的状态，同时还会受到环境和温度的影响,导致实验结果将偏离理论预测。为了使物理光学课程教学过程形象生动，提高学生学习光学课程的兴趣，应该充分利用计算机软件功能为教学增添活力，为学生理解复杂理论和实验做好铺垫，可以用MATLAB软件为平台,对光学实验进行计算机仿真，具有两个方面意义：一方面，利用仿真结果指导实际实验。前期投资少，且可以减少贵重仪器的损伤等；另一方面，在教学上，将抽象难懂的概念、规律通过实验仿真生动、形象地表现出来，使学生更易于接受，具有明显的教学效果。夫琅禾费衍射：原理：

3、考虑一束波长为 的单色平行光照在边长为 a 的方孔上，在方孔屏的前方紧贴放置一焦距为 F 的凸透镜 在凸透镜的焦平 面 上 放 置 观 测 屏 观 测 衍 射 花 纹 。 如 图 1 所 示 ，以 方 孔 正 中 心为 坐 标 原 点 O，X 轴 和 Y 轴 分 别 平 形 于 方 孔 的 两 组 对 边 ，Z 轴垂直于方孔并指向观测屏。观测屏上任一点Px，y，F在球坐标中可表示为：将 方 孔 等 分 为 NN 个 小 正 方 形 ， 那么 每 个 小 正 方 形 的 边 长为 a N。 只 要 N 的 取 值 足 够 大 ，那么 小 正 方 形 可 以 看 成 点 波 源 发出 球 状 子

4、波 ， 每 个 子 波 中 与 OP 平 行 的 光 线 都 将 会 聚 到 P 点 。OP 与 Y 轴 的 夹 角 的 余 弦 为 sinsin， 所 以 沿 Y 轴 方 向 相 邻 的两 正 方 形 发 出 子 波 到 达 P 点 的 光 程 差 为 asinsin N。 假 定第 1行 第 1 列 的 小 正 方 形 波 源 到 达 P 点 的 初 相 位 为 0，那么 第 i 行 第 j列 的 小 正 方 形 波 源 到 达 P 点 的 初 相 位 为 2 i1asincosj1asinsin N。假 定 每 个 小 正 方 形 波 源 到 时 达P点 时 电 场 强 度 的 振 幅

5、为1，那么 第i行 第j列 的 小 正 方 形 波 源 到 达P点 的 电 场 强 度 可 用 谐 振 动方程表示 ：P 点 总 的 电 场 强 度 仍 然 满 足 谐 振 动 方 程 ，且 为 所 有 小 正 方形波源在 P 点引起的谐振动的合成 ：E=Acos(t+)=Eij分别为合成后 P 点总的电场强度的振幅和初相位 。在 以 上 构 成 方 孔 的 NN 个 小 正 方 形 中 画 出 任 意 形 状 孔 ，具体画法是透光局部对应的小正方形用白色表示 不透光局部对应的小正方形用黑色表示 用一个 NN 的矩阵 R 记录孔的形状 即如果第 i 行第 j 列的小正方形为白色 那么 1 如果

6、第 i行 第 j 列 的 小 正 方 形 为 黑 色 ，那么Rij0。 在 谐 振 动 的 合 成 时 只 让白色的小正方形波源在 P 点引起的电场强度谐振动方程参与叠 加 ，那么 P 点 总 的 电 场 强 度 E 可 用 以 下 公 式 求 得 ：NNE=Acos(t+)=RijEij理 论 上 ，根 据 以 上 公 式 ，只 要 知 道 任 意 两 个 时 刻 的 电 场 强度 E 的 大 小 ，即 可 求 得 振 幅 A 的 大 小 。利 用 计 算 机 计 算 出t0 及t 2 时 P 点 总 的 电 场 强 度 E1 和 E2：利用以上公式可以计算出观测屏上任一点的光强 如果借助计

7、算机编程计算出观测屏上足够多点的光强 就可以模拟出任意形状孔的夫琅禾费衍射花纹方法：假 定 方 孔 的 边 长 为 001mm，单 色 光 波 长 为 600nm，凸 透 镜焦 距 为 05m，方 孔 二 值 图 像 矩 阵 R 的 元 素 全 取 1，那么 方 孔 夫 琅禾 费 衍 射 的 矩 阵 解 可 由 以 下 MATLAB 程 序 求 得Rones100，100；pi31415926；F05； 凸透镜焦距Azeros200，200；B105；for m1：200 x1005m 1000for n1：200y1005n 1000；c2*pi*B*sincita 100*600*109；

8、citaatansqrtx2y2 Z; if x0faiatany x；elsefaipiatany x；endE10；E20；for i1：100for j1：100if Ri，jE1E1cosi1*cosfaij1*sinfai*cE2E2cospi 2i1*cosfaij1*sinfai*cendendendA m，nsqrtE12E22； endend矩阵 A 即为方孔夫琅禾费衍射的矩阵解.只要将上述方孔衍射程序中的矩阵 R 替换为任意形状孔的二值图像矩阵 ， 即可求得任意形状孔夫琅禾费衍射的矩阵假 定 0.01mm，=600mm，F0.5m，N200， 如 图 2 所 示 ，由 上

9、节 介 绍 方 法 在 MATLAB 中 编 程 求 得 了 方 孔 、圆 孔 、直 角 等腰 三 角 形 孔 、半 圆 孔 的 衍 射 花 纹方孔 圆孔及三角孔等夫琅禾费衍射的计算机模拟结果与理论计算的结果相符 说明这种计算机模拟方法是可行的 半圆孔夫琅禾费衍射的理论计算方法还未见报道 但通过这种方法可以得出模拟的衍射花纹菲涅尔衍射：如图 1 所示, 设 Q ( x , y , z ) 为单色光源, S 1 为衍射孔, d s 为 S 1 面上的一个面积元, P ( x 0, y 0, z )为前方的任一点. 根据惠更斯菲涅尔原理,P ( x 0, y 0, z ) 点的场是由开孔平面的无穷

10、多个虚设的次波源产生的, 而次波源的复振幅与入射波在该点的复振幅和面积元 d s 成正比, 与波长成反比由上述的假设, 单色光源 Q ( x , y , z ) 发出的光波照射到开孔面 S 1 之后, 任一点 P ( x 0, y 0, z )处产生的光振动的复振幅表示为基尔霍夫衍射公式其中, A 为常 数, 为光源的波长, k 为波矢, r , s分别为波面 S 1 到 Q ( x , y , z ) 和 P ( x 0, y 0, z ) 的距离, ( n , r) 和( n, s) 分别为开孔平面的法线与 r 和s 方向的夹角.基尔霍夫衍射公式, 它比拟复杂, 因而直接用它来计算非常困难

11、. 在实际进行运算时通常是对这个普遍理论作某些近似, 用所得的近似公式计算一定范重内的衍射场分布, 按照近似条件的不同, 可以分为菲涅尔衍射和夫琅禾费衍射. 在引进一系列近似以前, 先研究一下基尔霍夫衍射公式在直角坐标系中的一般形式如上图所示, 以任一个形状的开孔面为坐标原点 O, 观察屏与衍射屏( 开孔面) 相距为 z , 衍射圆孔的透过率函数为 f ( x , y ) , 那么位于( x , y , z ) 处的单色光源 Q 在观察屏上的任一点 P 处所产生的衍射光场为: 上式只是对式( 1. 1 ) 的 积分式中分母的 r 作了初步近似式, 其中 r 表示为衍射屏和观察面之间的距离. 对

12、式( 1. 2 ) 中的指数 r 作二项式展开得到将上式代入( 1. 2) 式, 衍射叠加积分式化为式( 1. 3) 成为菲涅耳近似. 这个近似式成立的区域成为菲涅耳衍射区, 由此近似式导出的公式( 1. 4 )那么成为菲涅耳衍射公式将菲涅耳衍射公式( 1. 4) 中的指数因子展开并稍加整理后有式中, f x = x 0/ z , f y = y 0/ z . 将式( 1. 5) 代入式( 1. 4) , 衍射公式写为为了把式( 1. 6) 写成和卷积相同的形式, 利用卷积的符号, 将式( 1. 6) 改写成uP ( x 0, y 0) = f ( x , y ) * hz ( x 0, y

13、0) , ( 1. 7)式中由于 f ( x , y ) 是 把点光 源放在 原点时, 式 子 exp( j kr / r ) 中的 r 进行二 项式展开的近 似式, 所以,称之为点光源传递函数. 采用式( 1. 6) 计算, 或采用式( 1 . 7) 计算, 所得的结果应该是 完全相等的. 可是, 在各种不同的场合, 选择何种计算方法, 计算的难易程度是不同的. 一般, 在必须对 u p ( x , y ) 进行傅里叶变换时, 往往式( 1. 7) 更为方便.显然, 菲涅耳衍射公式也可以看成是乘积的傅里叶变换在菲涅耳近似近似条件下, 如果进一步增大观察面与衍射屏之间的距离 z , 使式( 1

14、. 3 ) 中因子( x 2+ y 2) 对相位的影响可以忽略, 即满足这就是夫琅 和费近似条件. 由它规 定的 z 值范围称为夫琅和费衍射区, 在这个区域内产生的衍射,称为夫琅和费衍射.将菲涅耳衍射公式( 1. 4 ) 中的 e 指数展开并忽略去 项，就得到夫琅和费尔衍射公式显然, 夫琅和费衍射公式也可以看成是 f ( x , y ) 的傅里叶变换. 对于具有规那么形状的衍射物来说, 利用傅里叶变换式( 1 . 9) 来计算夫琅和费衍射把戏要简单得多矩孔菲涅耳衍射的模拟计算菲涅耳衍射是在菲涅耳近似条件成立的距离范围内所观察的衍射现象. 相对于夫琅和费衍射而言, 观察屏距衍射屏不太远. 在菲涅

15、尔衍射中, 输入变量和输出变量分别为衍射孔径平面的光场分布和观察平面的光场以及光强分布. 考虑到这三个量都是二维分布, 而且 M at lab 主要应用于矩阵数值运算, 故为了提高运算速度和精度以及便于程序编写, 本程序选择用二维矩阵来存储衍射孔径平面和观察平面的场分布, 并分别以矩阵的列数和行数来对应平面的直角坐标值( x , y ) 以及( x 0, y 0) .选取直角坐标系原点与矩孔的中心重合( 如图3 所示) , 矩形的边分别与 x 和 y 轴平行, 边长分别为 2 W x 和 2 W y , 矩孔的振幅透射率函数设有一单位振幅的单色平面波垂直照射上述孔径, 在孔径之后的复场分布为由

16、( 1. 4) 式, 观察平面上矩孔菲涅耳衍射的复振幅为( 常数 A 可以忽略)( 1. 12)式所表示的积分可以别离为两个单一变量积分的乘积, 形成如下式中因而得到( a) 矩孔边长 a= b= 0. 5 10- 3 m, 孔径距离 d= 0. 005 m( b) 矩孔边长 a= b= 0. 5 10 - 3 m, 孔径距离 d= 0. 05 m( c) 矩孔边长 a= b= 0. 5 10- 3 m, 孔径距离 d= 0. 155 m( d) 矩孔边长 a= b= 0. 5 10 m, 孔径距离 d= 2 m变量代换后的积分限由以下四式给出:为 计 算 ( 1. 17) 式中 的 积分,

17、可用 式定义的菲涅耳积分和积分限( 1. 18) 式, 通过进行上述计算, 最后可得到观察屏幕上的光强分布对于给定的参数和坐标 x 0 和 y 0, 先由( 1. 18 )式求出积分限 1, 2, 1, 2 之值, 然后由菲涅耳积分的数值表查出 C ( ) , C ( ) , S ( ) , S ( ) 的值,最后将这些值代入( 1. 19 ) 式就得到 I ( x 0, y 0) 的数值解. 但这是一种非常繁杂的计算过程, 为防止上述麻烦, 对( 1. 19 ) 式用 M at lab 进行模拟很容易得到结果.在程序中波长取= 589. 3*10- 9 m, 矩孔边长 W x = W y =

18、 a = b= 0. 5*10 - 3 m, 观察屏幕到孔径距离 z = d = 0 . 005 m 2 m 取不同值, 输出的衍射图如图 4 所示.由( 1. 19) 式直接得到的结果, 与林继成 8 讨论结果完全一致. 程序如下:a= 0. 5E- 3;b= 0. 5E- 3; % apert ure lengt h and w idt hlambda= 589. 3E- 9; % incident w aveleng thd= 0. 005; % mask dist ancepoints= 500;x= linspace( - a* 2, a* 2 , points) ;y= linsp

19、ace( - b* 2 , b* 2 , points) ;u1= - sqrt ( 2/ ( lambda* d ) ) . * ( a+ x0) ;u2= sqrt ( 2/ ( lam bda* d) ) . * ( a- x0) ;Cu1= mf un( FresnelC , u1) ;Cu2= mf un( FresnelC , u2) ;Su1= mf un( FresnelC , u1) ;Su2= mf un( FresnelC , u2) ;Ix = ( 1/ 2) * ( ( ( Cu2- Cu1) . 2) + ( ( Su2-Su1 ) . 2 ) ) ;v 1= -

20、sqrt ( 2/ ( lambda* d ) ) . * ( b+ y0) ;v 2= sqrt ( 2/ ( lambda* d) ) . * ( b- y0) ;Cv1 = mf un( F resnelC , v 1) ;Cv2 = mfun( FresnelC , v2) ;Sv1 = mf un( F resnelC , v 1) ;Sv2 = mfun( FresnelC , v2) ;Iy = ( 1/ 2) * ( ( ( Cv2- Cv1) . 2) + ( ( Sv2- Sv1 ) . 2 ) ) ;I= Ix * Iy; subplot ( 121)imagesc( I

21、/ max ( I ( : ) ) , 0, 0. 9 ) ; colormap( hot ) ;subplot( 122)plot( x 0* 1E6 , I( end/ 2, : ) ) ;圆孔菲涅耳衍射的计算模拟对卷积符号表示的菲涅耳衍射公式( 1. 7) , 用Matlab 软件计算模拟, 可以得图 5、图 6 所示的相对光强度曲 线和衍射 图样. 在程序 中波长 取 =632. 8e- 6 m, 圆孔半径 r= 15 m, 观察屏幕到孔径距离 z= 1000000 m, 输出的衍射图如图 5 和图6 所示. 在 模拟过程中可改变光束传输距离 z 、圆孔半径 r 和光的波长 得到各种效

22、果的衍射图样程序如下:N= 300;r= 15; a= 1; b= 1;I= zeros( N, N) ; m, n = meshg rid( linspace( - N/ 2, N/ 2- 1, N) ) ;D= ( ( m- a) . 2+ ( n- b) . 2) . ( 1/ 2) ;i= f ind( Dlmda=500e-9; r=1e-3; f=1; N=19; K=linspace(-0.1,0.1,N); lmda1=lmda*(1+K); xm=2000*lmda*f; xs=linspace(-xm,xm,2000); ys=xs; z0=zeros(2000); x,y

23、=meshgrid(xs); for i =1:19 s=2*pi*r*sqrt(x.2+y.2)./(lmda1(i); z=4*(besselj(1,s)./(s+eps).2; z0=z0+z; end z1=z0/19; subplot(1,2,1) imshow(z1*255); title(夫朗禾费圆孔衍射光强) xlabel(x) ylabel(y) subplot(1,2,2) mesh(x,y,z1) colormap(hot) xlabel(x) ylabel(y) zlabel(光强) = 2 * GB3 夫琅禾费矩孔衍射Matlab程序： f=1;a=0.0001;b=

24、0.0001;lmda=632.8e-9;xmax=0.05;ymax=xmax;def=0.0001;x=-xmax:def:xmax;y=-ymax:def:ymax;lenm=length(x);lenn=length(y);for m=1:lenm for n=1:lenn alpha=pi*x(m)*a/(lmda*f); beta=pi*y(n)*a/(lmda*f); I(m,n)=(sin(alpha)/(alpha)2*(sin(beta)/(beta)2; endend I=I/(max(max(I); X,Y=meshgrid(x,y); subplot(1,2,1) i

25、mshow(255*I); title(夫朗禾费矩孔衍射光强) xlabel(x); ylabel(y); subplot(1,2,2) mesh(X,Y,I); colormap(hot) xlabel(x); ylabel(y); zlabel(光强); = 3 * GB3 夫琅禾单缝缝衍射Matlab程序：lam=500e-9;a= 2e-4;z=5;d=5*a;xm=2*lam*z/a;y0=xm;n=1001;x0=linspace(-xm,xm,n);for i= 1: nsinphi=x0(i)/z;alpha=pi*a*sinphi/lam;beta=pi*d*sinphi/l

26、am;B(i,:)=(sin(alpha)./alpha).2.*(sin(beta)./sin(beta).2;B1=B/max(B);endNC=255;Br=(B/max(B)*NC;subplot(1,2,1)image(y0,x0,Br);colormap( hot );subplot(1,2,2)plot(Br); = 4 * GB3 夫琅禾费多缝衍射Matlab程序：lam=500e-9;N=2;a= 2e-4;z=5;d=5*a;xm=2*lam*z/a;y0=xm;n=1001;x0=linspace(-xm,xm,n);for i= 1: nsinphi=x0(i)/z;a

27、lpha=pi*a*sinphi/lam;beta=pi*d*sinphi/lam;B(i,:)=(sin(alpha)./alpha).2.*(sin(N*beta)./sin(beta).2;B1=B/max(B);endNC=255;Br=(B/max(B)*NC;subplot(1,2,1)image(y0,x0,Br);colormap( hot );subplot(1,2,2)plot(Br);菲涅尔衍射：根据惠更斯菲涅尔原理的次波假设和次波相干叠加原理，当次级波源处于同一波真面上时，由它们发出的子波必然是彼此相干的，在波传播的后面空间中任何一点处的光振动那么是这些次级波源产生的子

28、波叠加的结果。菲涅尔原理数学表达式为： 式中K(H)为倾斜因子 = 1 * GB3 菲涅耳圆孔衍射Matlab程序： N=300; r=15; a=1;b=1; I=zeros(N,N); m,n=meshgrid(linspace(-N/2,N/2-1,N); D=(m-a).2+(n-b).2).(1/2); i=find(D I(i)=1; subplot(2,2,1); imagesc(I) colormap(0 0 0;1 1 1) axis image title(衍射前图样) L=300; M=300; x,y=meshgrid(linspace(-L/2,L/2,M); lam

29、da=632.8e-6; k=2*pi/lamda; z=1000000; h=exp(j*k*z)*exp(j*k*(x.2+y.2)/(2*z)/(j*lamda*z); H=fftshift(fft2(h); B=fftshift(fft2(I); G=H.*B; U=fftshift(ifft2(G); Br=(U/max(U); subplot(2,2,2); imshow(abs(U); axis image; colormap(hot) title(衍射后图样); subplot(2,2,3); mesh(x,y,abs(U); subplot(2,2,4); plot(abs(Br) = 2 * GB3 菲涅耳单缝衍射Matlab程序：lam=500e-9;a= 1e-3; f=1;xm= 3*lam*f/a;nx= 51;xs=linspace(-xm,xm,nx);np=51;xp=linspace(0,a,np);for i=1:nxsinphi= xs(i)/f;alpha=2*pi*xp*sinphi/lam;sumcos=su