MATLAB光学模拟仿真实训报告模

学号：20163152616姓名：彭宏辉专业年级：光电16MATLAB光学模拟仿真实训目的1、熟悉matlab绘图和仿真功能2、复习物理光学等相关知识3、掌握运用matlab软件的仿真等功能与专业知识相结合实训内容MATLAB基础训练及光波在介质分界面的反射和折射相关原理当一个单色平面波射到两种不同介质的分界面时，将分为两个波，一个折射波一个反射波。从电磁场的边值关系可以求出它们的传播方向和入射波的振幅关系相位关系。实训任务已知界面两侧的折射率n2、n1和入射角，绘出在n1<n2（光由光疏介质射向光密介质）和n1>n2（光由光密介质射向光疏介质）两种情况下，反射系数、透射系数随入射角的变化曲线。程序如下：clear;%清空disp('请输入介质折射率n1和n2');%在显示括号内语句n1=input('n1=');%接受键盘任意输入合适的折射率n1n2=input('n2=');%接受键盘任意输入合适的折射率n2theta=0:0.1:90;%入射角范围范围0~90，步距0.1a=theta*pi/180;%角度化为弧度rp=(n2*cos(a)-n1*sqrt(1-(n1/n2*sin(a)).^2))./(n2*cos(a)+n1*sqrt(1-(n1/n2*sin(a)).^2));%p分量振幅反射率rs=(n1*cos(a)-n2*sqrt(1-(n1/n2*sin(a)).^2))./(n1*cos(a)+n2*sqrt(1-(n1/n2*sin(a)).^2));%s分量振幅反射率tp=2*n1*cos(a)./(n2*cos(a)+n1*sqrt(1-(n1/n2*sin(a)).^2));%p分量振幅透射率ts=2*n1*cos(a)./(n1*cos(a)+n2*sqrt(1-(n1/n2*sin(a)).^2));%s分量振幅透射率figure(1);%创建一个窗口subplot(1,2,1);%作图rp,rs,|rp|,|rs|随入射角的变化曲线plot(theta,rp,'-',theta,rs,'--',theta,abs(rp),':',theta,abs(rs),'-.','LineWidth',2);%用‘-’、‘--’，‘：’、‘-.’符号标注对应曲线legend('rp','rs','|rp|','|rs|');%标注曲线图例xlabel('入射角/theta_i');%命名x轴ylabel('振幅');%命名y轴title(['n_1=',num2str(n1),',n2=',num2str(n2),'时反射系数随入射角的变化曲线']);%命名图像axis([090-11]);%设定作图区间gridon;%添加网格subplot(1,2,2);%tp,ts,|tp|,|ts|随入射角的变化曲线plot(theta,tp,'-',theta,ts,'--',theta,abs(tp),':',theta,abs(ts),'-','LineWidth',2);%用‘-’、‘--’，‘：’、‘-.’符号标注对应曲线legend('tp','ts','|tp|','|ts|');%标注曲线图例xlabel('入射角/theta_i');%命名x轴ylabel('振幅');%命名y轴title(['n_1=',num2str(n1),',n2=',num2str(n2),'时透射系数随入射角的变化曲线']);%命名图像ifn1<n2%如果此时从光疏到光密axis([09001]);%设定作图区间else%否则axis([09003.5]);%设定作图区间end%结束gridon;%添加网格光波的叠加相关原理光波在空间某一区域相遇时，发生光波叠加现象。波的叠加服从叠加原理，波在相遇点产生的合振动是各个波单独产生的振动的矢量和。实训任务两列单色平面波的模拟。程序如下：w1=50;%波1频率w2=60;%波2频率k1=5;%波1波数k2=4;%波2波数t=0.1:0.2:1.3;%对时间进行等间隔取点a=1;%波动振幅x=0:0.001:5;%对传播方向x轴进行等间隔取点A2=a*cos(k2*x-w2*t(end));%A2波动函数A1=a*cos(k1*x-w1*t(end));%A1波动函数plot(x,A1,'-',x,A2,':')%绘制两个波形图set(gcf,'color',[111]);%添加颜色对比度set(gca,'YTick',[-1:0.5:1]);%标记Y轴坐标刻度set(gca,'XTicK',[0:1:5]);%标记X轴坐标刻度xlabel('变量X')%添加x轴说明ylabel('振幅')%添加y轴说明title('两列单色平面波的模拟')%添加标题legend('光波1','光波2')%添加图例光波1和光波2绘制合成波光强曲线。程序如下：w1=50,w2=60,k=0,k1=5,k2=4,a=1;%两列波的参数x=0:0.001:30;%传播方向X轴进行等间隔取点t=0.1:0.2:1.3%对时间进行等间隔取点k=k+1;%累加A=2*a*cos((k1-k2)/2*x-(w1-w2)/2*t(end));%波动函数I=A.*A;%合成光强plot(x,I);%绘制图形set(gca,'YTick',[0:1:4])%添加Y轴坐标刻度set(gca,'XTick',[0:5:30])%添加X轴做坐标刻度xlabel('变量X')%添加X轴说明ylabel('振幅变化')%添加Y轴说明title('合成波光强曲线')%添加标题两列单色平面波合成的动态仿真。程序如下：w1=50,w2=60,k1=5,k2=4,a=1;%两列波的参数x=0:0.001:30;%传播方向X轴进行等间隔取点k=0;%赋初值m2=moviein(length(0.1:0.2:1.3));%形成动态函数fort=0.1:0.2:1.3%对时间进行等间隔取点k=k+1;%累加A=2*a*cos((k1-k2)/2*x-(w1-w2)/2*t);%波动函数v=a*cos(k1*x-w1*t)+a*cos(k2*x-w2*t);%波动函数plot(x,v,'k-',x,A,'g:',x,-A,'b-.');%绘制图形axis([030-22]);%设置坐标轴set(gcf,'color',[111])%设置颜色对比度set(gca,'YTick',[-2:1:2])%添加Y轴坐标刻度set(gca,'XTick',[0:5:30])%添加X轴做坐标刻度xlabel('变量X')%添加X轴说明ylabel('振幅变化')%添加Y轴说明title('光学拍')%添加标题legend('合成波振幅','包络线1','包络线2')%添加图例m2(:,k)=getframe;%获取坐标为界的图像end%循环结束movie(m2,3)%绘制动态函数光波的干涉相关原理满足一定条件的两列相干波相遇叠加，在叠加区域某些点的振动始终加强，某些点的振动始终减弱，即在干涉区域内振动强度有稳定的空间分布。实训任务杨氏双缝干涉的模拟。程序如下：clearalllam=500e-9;%设置输入波长a=2e-3;%设置屏幕到双缝平面距离D=1;%设置双缝距离ym=5*lam*D/a;%设定光屏范围xs=ym;%赋值n=101;ys=linspace(-ym,ym,n);%把光屏Y方向分成101点fori=1:n%开始循环r1=sqrt((ys(i)-a/2).^2+D^2);%两个相干光源到到屏幕上任意点P的距离r2=sqrt((ys(i)+a/2)^2+D^2);%两个相干光源到到屏幕上任意点P的距离phi=2*pi*(r2-r1)/lam;%相位差B(i,:)=4*cos(phi/2).^2;%p点的光强end%循环结束N=255%确定灰度级为255级Br=(B/4.0)*N;%最大光强对应的灰度级subplot(1,2,1);%将一行二列的图从左到右从上到下的第一个位置image(xs,ys,Br);%画干涉条纹colormap(gray(N));%输出一个灰色的曲线图subplot(1,2,2);%将一行二列的图从左到右从上到下的第二个位置plot(B,ys);%画光强变化曲线模拟非定域干涉下点光源的干涉条纹。程序如下：xmax=40;ymax=40;%干涉条纹范围设定Lambad=632.8e-006;%光波长的选择，此处用氦氖激光器f=190;%扩束镜的焦距n=1.0%介质折射率，空气为1N=800;%精度或分辨率，决定x,y的数据量x=linspace(-xmax,xmax,N);%产生-xmax，xmax之间的N点行线性的矢量y=linspace(-ymax,ymax,N);%产生-xmax，xmax之间的N点行线性的矢量fork=0:15%模拟增加光程差次数d=0.39-0.00005*k;%模拟增加光程差fori=1:N%循环体forj=1:N%循环体r(i,j)=sqrt(x(i)*x(i)+y(j)*y(j));B(i,j)=cos(pi*(2*n*d*cos(asin(n*sin(atan(r(i,j)/f)))))/Lambad).^2;end%计算不同光程差条件下的条纹位置endfigure(gcf);%建立N1=255;%设置分辨率或精度Br=2.5*B*N1;%调整最终的条纹位置image(x,y,Br);%显示仿真结果colormap(gray(N1));%灰度图显示set(gca,'XTick',[]);%设置x坐标set(gca,'YTick',[]);%设置y坐标axisequal%坐标轴设置drawnow%立即绘图pause%保留图像end模拟定域干涉下的等倾干涉条纹。程序如下：d=2;%设置透镜厚度Lambad=632.8e-006;%光波长的选择，此处用氦氖激光器theta=0.15;%设置初始的入射倾角rmax=d*tan(theta/2)%设置标度N=800;%精度或分辨率，决定x,y的数据量x=linspace(-rmax,rmax,N);%设置x轴范围和标度y=linspace(-rmax,rmax,N);%设置y轴范围和标度fork=0:15%模拟增加倾角次数theta=theta-0.01*k%倾角减小rmax=d*tan(theta/2)%每条等倾条纹的新标度fori=1:N%循环体forj=1:N%循环体x(i)=(i-2)*2*rmax/(N-1)-rmax;%等倾条纹的x坐标y(j)=(j-2)*2*rmax/(N-1)-rmax;%等倾条纹的y坐标r(i,j)=sqrt(x(i)^2+y(j)^2);%等倾条纹的半径delta(i,j)=2*d/sqrt(1+r(i,j)^2/d^2);%中间参量角计算Phi(i,j)=2*pi*delta(i,j)/Lambad;%计算参量PhiB(i,j)=4*cos(Phi(i,j)/2)^2;%计算不同倾角条件下的条纹位置endendfigure(gcf);%创建图像显示窗口N1=255;%设置分辨率或精度Br=2.5*B*N1;%调整最终的条纹位置image(x,y,Br);%显示仿真结果colormap(gray(N1));%灰度图显示set(gca,'XTick',[]);%设置x坐标set(gca,'YTick',[]);%设置y坐标axisequal%坐标轴设置drawnow%立即绘图pause%保留图像end光的衍射相关原理光在传播路径中，遇到不透明或透明的障碍物或者小孔（窄缝），绕过障碍物，产生偏离直线传播的现象称为光的衍射。衍射时产生的明暗条纹或光环，叫衍射图样实训任务绘出单色光的衍射光强曲线。程序如下：clear;Lambda=input('输入光的波长(单位为nm):(取500)');%输入光波长Lambda=Lambda*(1e-9);%输入米*10^-9aWidth=input('输入缝宽(单位为mm):(取0.2)');%输入缝宽aWidth=aWidth*0.001;%输入米*0.0001Z=input('输入缝到屏的距离(单位为):(取1)');%输入缝到屏的距离ymax=3*Lambda*Z/aWidth;%屏幕范围沿沿y轴Ny=51;%屏幕上的点数ys=linspace(-ymax,ymax,Ny);%将屏幕分成Ny个点NPoints=51;%单缝的点数yPoint=linspace(-aWidth/2,aWidth/2,NPoints);%把单缝分成NP个光源forj=1:Ny%对屏幕上y方向各点作循环L=sqrt((ys(j)-yPoint).^2+Z^2);%光源到接收屏的距离Phi=2*pi.*(L-Z)./Lambda;%从距离差计算各波长的相位差SumCos=sum(cos(Phi));%余弦叠加SumSin=sum(sin(Phi));%正弦叠加B(j)=(SumCos^2+SumSin^2)/NPoints^2'%叠加各波长并影响计算光强end%循环结束plot(ys,B,'*',ys,B);grid;%画光强曲线axis([-ymax,ymax,0.0,1.0]);%y轴的区间控制在1内set(gcf,'color,','w');%将背景设置为白色并计算光强模拟单狭缝的衍射图样和光强分布曲线。程序如下：clearlam=500e-9;%设置输入波长a=1e-3;%设置屏幕到双缝平面距离D=1;%设置双缝距离ym=3*lam*D/a;%屏幕范围沿沿y轴ny=51;%屏幕上的点数ys=linspace(-ym,ym,ny);%将屏幕分成ny个点np=51;%单缝的点数yp=linspace(0,a,np);%把单缝分成NP个光源fori=1:ny%循环开始sinphi=ys(i)/D;%计算余弦值alpha=2*pi*yp*sinphi/lam;%计算角度值sumcos=sum(cos(alpha));%余弦叠加sumsin=sum(sin(alpha));%正弦叠加B(i,:)=(sumcos^2+sumsin^2)/np^2;%将光强赋值给B矩阵end%循环结束N=255;%确定灰度级为255级Br=(B/max(B))*N;%最大光强对应的灰度级subplot(1,2,1)%将一行二列的图从左到右从上到下的第一个位置image(ym,ys,Br);%画干涉条纹colormap(gray(N));%输出一个灰色的曲线图subplot(1,2,2)%将一行二列的图从左到右从上到下的第二个位置plot(B,ys);%画光强变化曲线模拟圆孔的夫琅禾费衍射。程序如下：clear;lam=632.8e-9;%波长R=5e-4;%孔径大小f=1;%透镜焦距ym=2.5e-3;%y轴的范围m=1000;%画的点数y=linspace(-ym,ym,m);%取值范围从-ym到ym，取值数量为mx=y%赋值fori=1:m%开始循环r=x(i).^2+y.^2;%圆孔半径s=sqrt(r./(r+f^2));%光源到接收屏距离x1=2*pi*R*s./lam;%计算参数x1I(:,i)=((2*besselj(1,x1)).^2./x1.^2).*5000;%将光强赋值给I矩阵end%循环结束image(x,y,I);%绘制图形title('圆孔夫琅禾费衍射分布');%添加标题n=100;%设置灰度级colormap(gray(n));%设置图像为灰色colorbar;%添加色标模拟光栅衍射图样。程序如下：clear;lam=500e-9;%波长N=2;%设置灰度级a=2e-4,D=5,d=5*a;%设置参数ym=2*lam*D/a;xs=ym;%Y轴的范围n=1001;%赋值ys=linspace(-ym,ym,n);%把光屏的y向分成1001个点fori=1:n%开始循环sinphi=ys(i)/D;%余弦值alpha=pi*a*sinphi/lam;%计算角度值beta=pi*d*sinphi/lam;%计算角度值B(i,:)=(sin(alpha)./alpha).^2.*(sin(N*beta)./sin(beta)).^2;%将光强赋值给B矩阵B1=B/max(B);%将最大光强设为1end%结束循环NC=255;%确定灰度级Br=(B/max(B))*NC;%最大光强对应的灰度级subplot(1,2,1)%选中1*2区中的1号区image(xs,ys,Br);%画出干涉条纹colormap(gray(NC));%设置图像为灰色subplot(1,2,2)%选中1*2区中的2号区plot(B1,ys);%画出沿y向的相对光强变化曲线光的偏振相关原理两个光源发出的单色光波频率相同，振动方向互相垂直，一个平行于x轴，一个平行于y轴，考察它们在z轴方向上的叠加。由于两叠加光波的角频率，很容易看出合矢量沿椭圆旋转的角频率不变。我们把光矢量周期性旋转，其末端的运动轨迹为一个椭圆的成为椭圆偏振光。实训任务绘制偏振光波轨迹的二维图形。程序如下：clearall:%复位c=3e+8;%设置光速lamd=5e-7;%波长T=lamd/c;t=linspace(0,T,1000);%设置t取值范围从0到1000，取值数量为T个z=linspace(0,5,1000);%设置z取值范围从0到1000，取值数量为5个w=2*pi/T%角频率k=2*pi/lamd;%波数Eox=10;Eoy=5;%两波的单位矢量Fx=0;%设置平行于X轴光波的位相i=1;%设置i参数forFy=0:pi/4:7*pi/4%进行for循环，设置平行于Y轴光波的位相Ex=Eox*cos(w*t+k*z);%平行于X轴的波的光振动表达式Ey=Eoy*cos(w*t+k*z+Fy);%平行于Y轴的波的光振动表达式subplot(2,4,i);%将2行4列的图绘制到一个图里i=i+1;%i在for循环中的叠加plot(Ex,Ey);%以Ex为横坐标，Ey为纵坐标绘制图像n=i-2%n=i-2xlabel('x');%命名横坐标名称ylabel('y');%命名纵坐标名称title(['Fy-Fx=',num2str(n),'*pi/4']);%命名图像名称gridon%显示网格end绘制偏振光波轨迹的三维图形。程序如下：c=3e+8,lamd=5e-7,T=lamd/c;%设置光波参数t=linspace(0,T,1000);%设置t取值范围从0到1000，取值数量为T个z=linspace(0,5,1000);%设置z取值范围从0到1000，取值数量为5个q=linspace(0,5,1000)%设置q取值范围从0到1000，取值数量为5个w=2*pi/T%角频率k=2*pi/lamd;%波数Eox=10;Eoy=5;%两波的单位矢量Fx=0;%设置平行于X轴光波的位相i=1;%设置i参数forFy=0:pi/4:7*pi/4%进行for循环，设置平行于Y轴光波的位相Ex=Eox*cos(w*t+k*z);%平行于X轴的波的光振动表达式Ey=Eoy*cos(w*t+k*z+Fy);%平行于Y轴的波的光振动表达式subplot(2,4,i);%将2行4列的图绘制到一个图里i=i+1;%i在for循环中的叠加plot3(Ex,Ey,q);%以Ex为x坐标，Ey为y坐标，Ez为坐标绘制图像n=i-2;%设置参数nxlabel('x');%命名x坐标名称ylabel('y');%命名y坐标名称zlabel('z');%命名z坐标名称title(['Fy-Fx=',num2str(n),'*pi/4']);%命名图像名称gridon%显示网格end傅里叶变换相关原理将满足一定条件的某个函数表示成三角函数（正弦和/或余弦函数）或者它们的积分的线性组合可分析信号的成分，也可用这些成分合成信号。许多波形可作为信号的成分，比如正弦波、方波、锯齿波等，傅立叶变换用正弦波作为信号的成分。实训任务绘制圆孔输入物。程序如下：clear;%清除变量N=1024;%定义NI64=zeros(N,N);%预定义平面I64的灰度值为0[m,n]=meshgrid(linspace(-N/2,N/2-1,N));%确定坐标系及坐标原点的位置r=input('请输入圆的半径r=');%输入圆的半径大小控制变量ra=input('请输入圆心位置的控制变量a=');%输入圆心位置控制变量ab=input('请输入远新位置的控制变量b=');%输入圆心位置控制变量bD=((m+a).^2+(n+b).^2).^(1/2);%圆函数关系式i=find(D<=r);%f返回满足条件D<=r的像素点的单下标值I64(i)=1;%像素点赋值imshow(I64);%显示图像三角孔和矩形孔的傅里叶频谱观察。程序如下：clearimga=imread('fft4.jpg');%读取图片gridonfigure(1);%建立第一幅图imshow(imga,[]);%将范围的灰度值变换到0-255区间显示imgray0=rgb2gray(imga);%将图像转换为灰度图imgray1=double(imgray0);%提高图像的精度值afft=fft2(imgray1);%将图像做二维变换afftI=fftshift(afft);%将afft零频率的分量移到频谱的中心RR=real(afftI);%取afftl的实数部分II=imag(afftI);%取afftl的虚数部分A=sqrt(RR.^2+II.^2);%取实部与虚部的算术平方根A=(A-min(min(A)))/(max(max(A))-min(min(A)))*255;%归一化figure(2);%建立第二幅图imshow(A);%绘A的图形[m,n]=size(A);%获取矩阵A的大小[x,y]=meshgrid(1:n,1:m);%获取三维图像figure(4)%建立第三幅图mesh(x,y,A)%绘制三维网格图卷积定理的仿真。程序如下：clear;A=zeros(200,200);%生成一个200*200的零矩阵Afori=1:10%开始循环A(:,20*i-9:20*i)=1;%将1填充到矩阵A的20*i-9到20*i列end%循环结束B=zeros(200,200);%生成一个200*200的零矩阵Bfori=1:10%开始循环B(20*i-9:20*i,:)=1;%将1填充到矩阵B的20*i-9到20*i行end%循环结束C=A.*B;%矩阵C等于矩阵A*Bc=abs(fftshift(fft2(C)));%二次傅里叶变换cc1=abs(fftshift(fft2(C)));%二次傅里叶变换aa=abs(fftshift(fft2(A)));%二次傅里叶变换bb=abs(fftshift(fft2(B)));%二次傅里叶变换cc=conv2(aa,bb);%求卷积figure(1);%建立图形subplot(1,2,1);imshow(B);%选中1*2区中的1号区并显示图像subplot(1,2,2);imshow(bb);%选中1*2区中的2号区并显示图像figure(2);%建立图形subplot(1,2,1);imshow(A);%选中1*2区中的1号区并显示图像subplot(1,2,2);imshow(aa);%选中1*2区中的2号区并显示图像figure(3);%建立图形imshow(c);%显示图像figure(4);%建