計算機模擬光學衍射

1、長江大學物理科學與技術學院題目：計算機模擬光學衍射系別應用物理專業(yè)光電子技術班級應用物理1001姓名王飛學號201006927序號22老師陳海燕2012年11月28日課題背景及研究意義衍射是光波的一種重要特征，對光波衍射現(xiàn)象的討論是以惠更斯-菲涅爾原理為基 礎的。在大多數(shù)光學問題中，可以采用光場的標量衍射理論來描述和計算。光場的頻率 非常高，只能測量在一個比光的周期大得多的時間間隔內(nèi)的平均值，在計算衍射光波場 的光振動和光強分布時則利用菲涅爾-基爾霍夫衍射積分公式。在物理光學的光的衍射 教學過程中，如不借助實驗，學生很難理解理論、實驗原理和過程。而光學實驗的進 行一般要穩(wěn)定的環(huán)境、高精密的儀器

2、以及光路的調(diào)試，花費相當多的時間，并且很難調(diào) 試到理想的狀態(tài)，同時還會受到環(huán)境和溫度的影響，導致實驗結果將偏離理論預測。為了使物理光學課程教學過程形象生動，提高學生學習光學課程的興趣，應該充分 利用計算機軟件功能為教學增添活力，為學生理解復雜理論和實驗做好鋪墊，可以用 MATLAB軟件為平臺，對光學實驗進行計算機仿真，具有兩個方面意義：一方面，利用仿 真結果指導實際實驗。前期投資少，且可以減少貴重儀器的損傷等；另一方面，在教學 上，將抽象難懂的概念、規(guī)律通過實驗仿真生動、形象地表現(xiàn)出來，使學生更易于接受， 具有明顯的教學效果。夫瑯禾費衍射：原理：考慮一束波長為的單色平行光照在邊長為a的方孔上，

3、在方孔屏的后方緊貼放置 一焦距為F的凸透鏡在凸透鏡的焦平面上放置觀測屏觀測衍射花紋。 如圖1所示，以方孔正中心為坐標原點0, X軸和Y軸分別平形于方孔的兩組對邊，Z軸 垂直于方孔并指向觀測屏。觀測屏上任一點P（x，y，F(xiàn)）在球坐標中可表示為：/k rsiFisiriip.Bur-clcin ( /2+y2 /z )r arc f-an (丫/工) 當 _rO 時ipvr+frc rin ( y/.T ) 當 rtV。時將方孔等分為NXN個小正方形，則每個小正方形的邊長為a / N。只要N的取值足夠大，則小正方形可以看成點波源發(fā)出 球狀子波，每個子波中與OP平行的光線都將會聚到P點。OP與Y軸的

4、夾角的余弦為sine sin?,所以沿Y軸方向相鄰的兩 正方形發(fā)出子波到達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理論上，根據(jù)以上公式，只要知道任意兩個時刻的電場強度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 所示，由上節(jié)介紹方法在MATLAB中編程求得了方孔、圓孔、直角等 腰三角形孔、半圓孔的衍射花紋弓嶺h|句方孔圓孔及三角孔等夫瑯禾費衍射的計算機模擬結果與理

9、論計算的結果相符說明這 種計算機模擬方法是可行的 半圓孔夫瑯禾費衍射的理論計算方法還未見報道 但通過 這種方法可以得出模擬的衍射花紋菲涅爾衍射：如圖1所示，設Q ( x , y , z )為單色光源，S 1為衍射孔，d s為S 1面上 的一個面積元，P ( x 0, y 0, z為前方的任一點.根據(jù)惠更斯菲涅爾原理,P ( x 0, y 0, z )點的場是由開孔平面的無窮多個虛設的次波源產(chǎn)生的，而次波源的復 振幅與入射波在該點的復振幅和面積元d s成正比，與波長成反比由上述的假設，單色光源Q ( x , y , z )發(fā)出的光波照射到開孔面S 1之后，任一 點P ( x 0, y 0, z

10、)處產(chǎn)生的光振動的復振幅表示為“M邛刃=2pJJs 1 I12=(C(句C(+S g) Sf aJ72； (C ( 5 - C j /2 + S P2 ； - S( 1- 19；對于給定的參數(shù)和坐標X0和y 0,先由(1. 18 )式求出積分限1, 2, 1, 2之值，然 后由菲涅耳積分的數(shù)值表查出C ( ) , C ( ) , S ( ) , S ()的值，最后將這些值代 入(1. 19 )式就得到I ( x 0, y 0)的數(shù)值解.但這是一種非常繁雜的計算過程，為 避免上述麻煩，對(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);菲涅爾衍射：根據(jù)惠更斯菲涅爾原理的次波假設和次波相十疊加原理，當次級波源處于同一波真 面上時，由它們發(fā)出的子波必然是彼此相干的，在波傳播的后面空間中任何一點處的光 振動則是這些次級波源產(chǎn)生的子波疊加的結果。菲涅爾原理數(shù)學表達式為：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