计算机模拟光学衍射

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

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

3、在方孔屏的后方紧贴放置 一焦距为F的凸透镜在凸透镜的焦平面上放置观测屏观测衍射花纹。 如图1所示，以方孔正中心为坐标原点0, X轴和Y轴分别平形于方孔的两组对边，Z轴 垂直于方孔并指向观测屏。观测屏上任一点P（x，y，F）在球坐标中可表示为：/k rsiFisiriip.Bur-clcin ( /2+y2 /z )r arc f-an (丫/工) 当 _rO 时ipvr+frc rin ( y/.T ) 当 rtV。时将方孔等分为NXN个小正方形，则每个小正方形的边长为a / N。只要N的取值足够大，则小正方形可以看成点波源发出 球状子波，每个子波中与OP平行的光线都将会聚到P点。OP与Y轴的

4、夹角的余弦为sine sin?,所以沿Y轴方向相邻的两 正方形发出子波到达P点的光程差为asine sin? / N。假定第1 行第1列的小正方形波源到达P点的初相位为0,则第i行第j 列的小正方形波源到达P点的初相位为2n (i l) asine cos? +(j 1) asine sin? ) / 入 N。假定每个小正方形波源到时达?点时电场强度的振幅为1, 则第市第1列的小正方形波源到达P点的电场强度可用谐振动 方程表示：厂EOS(皿+灯)P点总的电场强度仍然满足谐振动方程，且为所有小正方形波 源在P点引起的谐振动的合成：E=Acos(3 t+中) = EEEij分别为合成后P点总的电场

5、强度的振幅和初相位。在以上构成方孔的NXN个小正方形中画出任意形状孔，具体画 法是透光部分对应的小正方形用白色表示 不透光部分对应的小正方形用黑色表示用 一个NXN的矩阵R记录孔的形状即如果第i行第j列的小正方形为白色 则=1 如果第i行第j列的小正方形为黑色，则Rij = 0。在谐振动的合成 时只 让白色的小正方形波源在P点引起的电场强度谐振动方程参与叠加，则P点 总的电场强度E可用以下公式求得：NNE=Acos(3 t+中) = EERijEij理论上，根据以上公式，只要知道任意两个时刻的电场强度E 的大小，即可求得振幅A的大小。利用计算机计算出3 t = 0及3 t =n / 2时P点总

6、的电场强度E1和E2：舟g (1)2皿呼珀-1)如n财&物寺+前= 1杯帽m+2打W)心如gzi崩口典砂)2 I 此：2XN则龙M*+* ,而光强正比于，利用以上公式可以计算出观测屏上任一点的光强如果借助计算机编程计算出观测屏上足够多点的光强就可以模拟出任意形状孔的夫琅禾费衍射花纹方法：假定方孔的边长为0. 01mm，单色光波长为600nm，凸透镜焦距为0. 5m，方孔二值图像矩阵R的元素全取1,则方孔夫琅禾费衍射的矩阵解可由以下MATLAB程序求得R=ones (100, 100)；pi = 3. 1415926；F=0. 5；凸透镜焦距A=zeros (200, 200)；B=10 八一5

7、；for m= 1： 200 x=(100. 5-m)/ 1000for n= 1： 200y=(100. 5-n)/ 1000；c = 2*pi*B*sin (cita) / (100*600*10 9)；cita = atan (sqrt (x2 + y2)/ Z);if x0fai = atan (y / x)；elsefai=pi + atan (y / x)；endE1 = 0； E2 = 0；for i = 1： 100for j = 1： 100if R (i, j)E1 = E1 + cos (i 1) *cos (fai) + (j 1) *sin (fai) *c)E2 =

8、 E2 + cos (pi / 2+ (i 1) *cos ( fai) + ( j 1) *sin ( fai) *c)endendendA ( m, n)=sqrt ( E1 八 2 + E2 八 2)；endend矩阵A即为方孔夫琅禾费衍射的矩阵解.只要将上述方孔衍射程序中的矩阵R替换为任意形状孔的二值图像矩阵，即可求得任意形状孔夫琅禾费衍射的矩阵假定 a =0.01mm,入=600mm, F=0.5m, N = 200, 如 图 2 所示，由上节介绍方法在MATLAB中编程求得了方孔、圆孔、直角等 腰三角形孔、半圆孔的衍射花纹弓岭h|句方孔圆孔及三角孔等夫琅禾费衍射的计算机模拟结果与理

9、论计算的结果相符说明这 种计算机模拟方法是可行的 半圆孔夫琅禾费衍射的理论计算方法还未见报道 但通过 这种方法可以得出模拟的衍射花纹菲涅尔衍射：如图1所示，设Q ( x , y , z )为单色光源，S 1为衍射孔，d s为S 1面上 的一个面积元，P ( x 0, y 0, z为前方的任一点.根据惠更斯菲涅尔原理,P ( x 0, y 0, z )点的场是由开孔平面的无穷多个虚设的次波源产生的，而次波源的复 振幅与入射波在该点的复振幅和面积元d s成正比，与波长成反比由上述的假设，单色光源Q ( x , y , z )发出的光波照射到开孔面S 1之后，任一 点P ( x 0, y 0, z

10、)处产生的光振动的复振幅表示为“M邛刃=2pJJs 1 I12=(C(句C(+S g) Sf aJ72； (C ( 5 - C j /2 + S P2 ； - S( 1- 19；对于给定的参数和坐标X0和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,矩

11、孔边长W x = W y = a = b= 0. 5*10 - 3 m, 观察屏幕到孔径距离z = d = 0 . 005 m2 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= li

12、nspace( - 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= - sqrt (

13、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);ly = ( 1/ 2) * ( ( ( Cv2- Cv1) . 2) + ( ( Sv2- Sv1 ) . 2 );I= Ix * ly; subplot ( 121)imagesc( I/ max ( I ( : )

14、) , 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和光的波长得到各种效果的衍射 图样0 5-010015020025030001C0200300程序如下：N= 300;

15、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=meshgrid(xs); for i =1:19

16、s=2*pi*r*sqrt(x.A2+y.A2)./(lmda1(i);z=4*(besselj(1,s)./(s+eps).A2;z0=z0+z;endz1=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(光强)夫琅禾费矩孔衍射Mat lab程序： f=1;a=0.0001;b=0.0001;lmda=632.8e-9;xmax=0.05;ymax=xm

17、ax;def=0.0001;x=-xmax:def:xmax;y=-ymax:def:ymax;lenm=length(x);lenn=length(y);for m=1:lenmfor n=1:lennalpha=pi*x(m)*a/(lmda*f);beta=pi*y(n)*a/(lmda*f);I(m,n)=(sin(alpha)/(alpha)A2*(sin(beta)/(beta)A2;endend I=I/(max(max(I); X,Y=meshgrid(x,y); subplot(1,2,1) imshow(255*I); title(夫朗禾费矩孔衍射光强) xlabel(x)

18、; ylabel(y); subplot(1,2,2) mesh(X,YI); colormap(hot) xlabel(x); ylabel(y); zlabel(光强);天垃于查走导:元强夫琅禾单缝缝衍射Mat lab程序：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/lam;B(i,：)=(sin(alpha)./alpha).A2.*(sin(beta

19、)./sin(beta).A2;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);夫琅禾费多缝衍射Mat lab程序: 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;alpha=pi*a*sinphi/lam;beta=pi*d*sinphi/lam;B(i,

20、：)=(sin(alpha)./alpha).A2.*(sin(N*beta)./sin(beta).A2;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);菲涅尔衍射：根据惠更斯菲涅尔原理的次波假设和次波相十叠加原理，当次级波源处于同一波真 面上时，由它们发出的子波必然是彼此相干的，在波传播的后面空间中任何一点处的光 振动则是这些次级波源产生的子波叠加的结果。菲涅尔原理数学表达式为：U(P) = ckU(P0)k(H)dsa

21、 r式中K(H)为倾斜因子菲涅耳圆孔衍射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).A2+(n-b).A2).A(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); lamda=632.8e-6; k=2*pi/lamda;

22、 z=1000000; h=exp(j*k*z)*exp(j*k*(x.A2+y.A2)/(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)衍射前图杆EC心M匚ii i1C020J工

23、衍射后圈样菲涅耳单缝衍射Mat lab程序： 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=sum(cos(alpha);sumsin=sum(sin(alpha);B(i,:)=(sumcosA2+sumsinA2)/npA2;endN=255;Br=(B/max(B)*N;subplot(1,2,1)image(xm,xs,Br);colormap(hot);subplot(1,2,2) plot(B,xs);菲涅耳多缝衍射Mat lab程序：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;alpha=pi*a*sin