高斯光束的matlab仿真

1、題目：根據(jù)高斯光束數(shù)學(xué)模型，模擬仿真高斯光束在諧振腔中某一位置處的歸一化強(qiáng)度分布并給出其二維、三維強(qiáng)度分布仿真圖；用Matlab讀取實(shí)際激光光斑照片中所記錄的強(qiáng)度數(shù)據(jù)（讀取照片中光斑的一個(gè)直徑所記錄的強(qiáng)度數(shù)據(jù)即可，Matlab讀取照片數(shù)據(jù)命令為imread）,用該數(shù)據(jù)畫(huà)出圖片中激光光斑的強(qiáng)度二維分布圖，與之前數(shù)學(xué)模型仿真圖對(duì)比。（如同時(shí)考慮高斯光束光斑有效截面半徑和等相位面特點(diǎn)，仿真高斯光束光強(qiáng)、光斑有效截面半徑以及等相位面同時(shí)隨傳播距離z的變化并給出整體仿真圖可酌情加分。）原始光斑如圖1所示，用imread命令讀入matlab后直接用imshow命令讀取即可，圖1 CCD采集的高斯光束強(qiáng)度分

2、布讀入的數(shù)據(jù)是一個(gè)224 X 244的矩陣，矩陣中的數(shù)值代表光強(qiáng)分布。用讀入的數(shù)據(jù)取中間一行（122行）畫(huà)出強(qiáng)度分布如圖2所示。圖2 實(shí)驗(yàn)測(cè)量高斯曲線(xiàn)用理論上的高斯曲線(xiàn)公式畫(huà)出理論高斯曲線(xiàn)如圖3所示。圖3 理論高斯曲線(xiàn)M文件如下：A=imread('D:documents作業(yè)激光原理與應(yīng)用高斯.bmp');A1=A(:,122);x1=1:1:224;x2=-100:1:100;a2=exp(-x2.2/10);figureimshow(A);axis offtitle('fontsize12CCD采集的高斯光束光強(qiáng)分布');figureplot(x2,a2,&

3、#39;linewidth',1,'color','b');axis(-40 40 0 1.2)title('fontsize12實(shí)驗(yàn)測(cè)量高斯曲線(xiàn)')figureplot(x1,A1,'linewidth',1,'color','r')title('fontsize12理論高斯曲線(xiàn)')axis(50 200 0 180)畫(huà)三維強(qiáng)度分布。取圖片矩陣的中間層，用mesh命令畫(huà)出三維圖如圖4所示。圖4 三維強(qiáng)度分布由于讀入的圖片有一行白邊，需要手動(dòng)去除掉，否則三維圖會(huì)有一邊整體豎起

4、來(lái)，影響觀(guān)察。最終的M文件如下。A=imread('D:documents作業(yè)激光原理與應(yīng)用高斯.bmp');high, width, color = size(A);x=1:width;y=1:high-1;mesh(x', y', double(A(2:224,:,1); grid onxlabel('x'),ylabel('y'),zlabel('z');title('三維強(qiáng)度分布');再用matlab仿真理論上傳播過(guò)程中高斯光束的變化這次先給出M文件：%Gaussian_propagatio

5、n.m%Simulation of diffraction of Gaussian Beamclear;%Gaussian Beam%N:sampling numberN=input('Number of samples(enter from 100 to 500)=');L=10*10-3;Ld=input('wavelength of light in micrometers=');Ld=Ld*10-6;ko=(2*pi)/Ld;wo=input('Waist of Gaussian Beam in mm=');wo=wo*10-3;z_ra

6、y=(ko*wo2)/2*103;sprintf('Rayleigh range is %f mm',z_ray)z_ray=z_ray*10-3;z=input('Propagation length (z) in mm');z=z*10-3;%dx:step size dx=L/N;for n=1:N+1 for m=1:N+1 %Space axis x(m)=(m-1)*dx-L/2; y(n)=(n-1)*dx-L/2; %Gaussian Beam in space domain Gau(n,m)=exp(-(x(m)2+y(n)2)/(wo2);%

7、Frequency axis Kx(m)=(2*pi*(m-1)/(N*dx)-(2*pi*(N)/(N*dx)/2; Ky(n)=(2*pi*(n-1)/(N*dx)-(2*pi*(N)/(N*dx)/2; %Free space transfer function H(n,m)=exp(j/(2*ko)*z*(Kx(m)2+Ky(n)2); endend%Gaussian Beam in Frequency domain FGau=fft2(Gau);FGau=fftshift(FGau);%Propagated Gaussian beam in Frequency domain FGau

8、_pro=FGau.*H; %Peak amplitude of the initial Gaussian beam Peak_ini=max(max(abs(Gau);sprintf('Initial peak amplitude is %f mm',Peak_ini)%Propagated Gaussian beam in space domain Gau_pro=ifft2(FGau_pro);Gau_pro=Gau_pro;%Peak amplitude of the propagated Gaussian beam Peak_pro=max(max(abs(Gau_p

9、ro);sprintf('Propagated peak amplitude is %f mm',Peak_pro)%Calculated Beam Width N M=min(abs(x);Gau_pro1=Gau_pro(:,M);N1 M1=min(abs(abs(Gau_pro1)-abs(exp(-1)*Peak_pro);Bw=dx*abs(M1-M)*103;sprintf('Beam width(numerical) is %fmm',Bw)%Theoretical Beam Width W=(2*z_ray)/ko*(1+(z/z_ray)2)

10、;W=(W0.5)*103;sprintf('Beam width(theoretical) is %fmm',W)%axis in mm scale x=x*103;y=y*103;figure(1);mesh(x,y,abs(Gau)title('Initial Gaussian Beam')xlabel('x mm')ylabel('y mm')axis(min(x) max(x) min(y) max(y) 0 1)axis squarefigure(2);mesh(x,y,abs(Gau_pro)title('p

11、ropagated Gaussian Beam')xlabel('x mm')ylabel('y mm')axis(min(x) max(x) min(y) max(y) 0 1)axis square程序主要根據(jù)高斯光束的傳播規(guī)律計(jì)算傳播過(guò)程中任意z處的高斯光強(qiáng)分布。運(yùn)行結(jié)果：Number of samples(enter from 100 to 500)=500wavelength of light in micrometers=0.568Waist of Gaussian Beam in mm=1ans =Rayleigh range is 5530.972982 mmPropagation length (z) in mm100000ans =Initial peak amplitude is 1.000000 mmans =Propagated peak amplitude is 0.210252 mmans =Beam?width(numerical) is 1.940000mmans =Beam?width(theoretical) i