matlab仿真光束的传输特性

一、课程设计题目：

用matlab 仿真光束的传输特性。

二、任务和要求

用matlab 仿真光束通过光学元件的变换。

① 设透镜材料为k9玻璃，对1064nm 波长的折射率为1.5062，镜片中心厚度为3mm ，凸面曲率半径，设为100mm ，初始光线距离透镜平面20mm 。用matlab 仿真近轴光线（至少10条）经过平凸透镜的焦距，与理论焦距值进行对比，得出误差大小。

② 已知透镜的结构参数为101=r ，0.11=n ，51=d ，5163.121=='n n （K9玻璃），502-=r ，0.12='

n ，物点A 距第一面顶点的距离为100，由

A 点计算三条沿光轴夹角分别为10、20、30的光线的成像。试用Matlab 对以上三条光线光路和近轴光线光路进行仿真，并得出实际光线的球差大小。

③ 设半径为1mm 的平面波经凸面曲率半径为25mm ，中心厚度3mm 的平凸透镜。用matlab 仿真平面波在透镜几何焦平面上的聚焦光斑强度分布，计算光斑半径。并与理论光斑半径值进行对比，得出误差大小。（方法：采用波动理论，利用基尔霍夫—菲涅尔衍射积分公式。）

2、用MATLAB仿真平行光束的衍射强度分布图样。（夫朗和费矩形孔衍射、夫朗和费圆孔衍射、夫朗和费单缝和多缝衍射。）

3、用MATLAB仿真厄米—高斯光束在真空中的传输过程。（包括三维强度分布和平面的灰度图。）

4、（补充题）查找文献，掌握各类空心光束的表达式，采用费更斯-菲涅尔原理推导各类空心光束在真空中传输的光强表达式。用matlab对不同传输距离处的光强进行仿真。

三、理论推导部分

将坐标原点选在透镜中心处，θ1=arcsin(y1/r),由n1*sinθ1=n2*sin θ2可得出θ2=arcsin(n1/n2)*(y1/r)，由几何关系可得到θ=θ2－θ1，则出射光线的斜率k=tan(θ2－θ1),当入射直线y=y1时，x1=d－(r －)2^1

r )，并设出射直线为y=k*x+b;由直线经过（x1,y1）即

2^(y

可求出b值，从而就可以求出射直线。由单透镜焦点计算公式1/f=－(n-1)*(1/r1－1/r2)可求得f=193.6858。

利用近轴光学公式i1=(l1-r1)*u1/r1，i11=n1*i/n11，u11=u1+i1-i11

l11=r1+r1*i11/u11和转面公式u2=u11，l2=l11-d1可以求得u11、

u22、l22、h2等。

入射光线的夹角为u1,设入射光线为y1=k1*x1+b1其中的斜率k1=－u1又由于入射光线经过经过（-100，0）就可以求出b1。由h1=l1*u1即为y1，当y1为定值时就可以得到第一个横坐标x0，再利用最后的出射光线公式y3=k3*x3+b3，k3=－u22,又因为最终出射经过（d+l22,0）可求出b3，利用转面公式h2=h1－d*u11,即为y3可求出第二个横坐标x00。再求在透镜中的直线斜率k2=（（h2

－h1）/（x00－x0）），y2=k2*x2+b2经过（x0,h1）即可求得b2值，从而即可求得三条直线。 实际光束求法同理。

利用菲涅耳近似公式

11])2

^122

)^1(2)^1(1[1ex p()1,1(1),(dy dx z y y x x ikz y x E z i y x E ??-+-+-

=λ 求衍射面上的光强要对孔径上的点求积分可以转换成对其x1,y1的微分求和，其中公式中的z1=f 。

2.（1）夫朗和费矩形孔衍射

若衍射孔为矩形则在透镜焦平面上得到的衍射图样如图，衍射图样的主要特征为衍射亮斑集中分布在两个相互垂直的方向上，并且x 轴上的亮斑宽度与y 轴亮斑宽度之比，恰与矩形孔在两个轴上的宽度相反。

其中的θ为θx,同样的β中的θ为θy,利用θx=x/f,θy=y/f进行求解。（2）夫朗和费圆形孔衍射

夫朗和费圆孔衍射的讨论方法和矩形孔衍射的讨论方法相同，只是由于圆孔的几何对称性，采用极坐标更为方便。

Ф=kaθ

（3）夫朗和费单缝衍射

对于前面讨论的夫朗和费矩形孔衍射，如果矩形的一个方向的尺寸比另一个方向大得多，则该矩形孔衍射就变成单缝衍射（如图），这时沿y方向的衍射效应不明显，只在x方向有亮暗变化的衍射图样。

实验中通过利用θ=x/f进行求解

（4）夫朗和费多缝衍射

夫朗和费多缝衍射装置如图，其每条狭缝均平行于y1方向，沿x1方向的缝宽为a，相邻狭缝的间距为d，在研究多缝衍射时，由于后透镜的存在使衍射屏上每个单缝的衍射条纹位置与位置无关。因此，

用平行光照射多缝时，其每一个单缝都要产生自己的衍射，形成各自一套衍射条纹。当每个单缝等宽时，各套衍射条纹在透镜焦平面上完全重叠，其总光强分布为它们的干涉叠加。

四、Matlab仿真部分

clear all

r=100;

n1=1.5163;

n2=1;%透镜的曲率半径为100mm，透镜的折射率n1=1.5，空气的折射率n2=1

d=3;

%x=77:0.1:320;

figure(1)

for n=-5:5

y1=0.1*n;

%hold on;

%plot(x1,y1);

a1=asin(y1/r);%入射角

a2=asin(n1/n2*(y1/r));%折射角

a=a2-a1;

k=tan(a);%出射光线的斜率

x1=sqrt(r^2-y1^2);

x2=x1-r+d;

b=y1+k*x2;

%出射光线经过(x2,y1)

x=-20:0.01:x2;%零坐标选在透镜中心，入射光线距透镜20mm，故x=-20

hold on

.

plot(x,y1);%平行光束

x3=x2:0.01:300;

y=-k*x3+b;%出射光线

h old on

plot(x3,y);

End

clear all

%透镜的结构参数

r1=10;

r2=-50;

l1=-100;

L1=-100;

n1=1.0;

d1=5;

n11=1.563

n2=1.563;

n22=1.0;

figure(1)

for n=-3:-1%沿光轴分别为1、2、3度的光线进行入射

%近轴光学成像公式第一个面

u1=n;

i1=(l1-r1)*u1/r1

i11=n1*i/n11;

u11=u1+i1-i11;

l11=r1+r1*i11/u11;

%转面公式

u2=u11;

l2=l11-d1;

%近轴光学成像公式第二个面

i2=(l2-r2)*u2/r2;

i22= n2*i2/n22;

u22=u2+i2-i22;

l22=r2+r2*i22/u22;

%入射光线与第一个透镜交点的纵坐标，坐标原点选在第一个透镜的顶点处

h1=l1*(u1*pi/180);

k1=-u1*pi/180;%入射光线的斜率

b1=100*k1;%因为入射光线经过（-100,0）点

x0=(h1-b1)/k1;%入射光线与第一个透镜交点的横坐标

x1=-100:0.01:x0;

y1=k1*x1+b1;

hold on

plot(x1,y1);%输出入射光线

k3=-u22*pi/180;%第二次折射后出射光线的斜率

b3=-k3*(d1+l22);%因为第二次折射后出射光线经过（d1+l22,0）点

h2=h1-d1*(u11*pi/180);%第一次折射后入射到第二个透镜的纵坐标

x00=(h2-b3)/k3;%第一次折射后入射到第二个透镜的横坐标

k2=(h2-h1)/(x00-x0);%第一次折射后光线的斜率

b2=h1-k2*x0;%因为第一次折射后光线经过（x0，h1）点

x2=x0:0.01:x00;

y2=k2*x2+b2;

hold on

plot(x2,y2);%输出第一次折射在两个透镜中的光线

x3=x00:0.01:30;%选在30是为了将输出图形看得更清晰些

y3=k3*x3+b3;

hold on

plot(x3,y3);%输出经过第二个透镜后的输出光线

%实际光路

U1=n*pi/180;

I1=asin((L1-r1)*sin(U1)/r1);

I11=asin(n1*sin(I1)/n11);

U11=U1+I1-I11;

L11=r1+r1*sin(I11)/sin(U11);

%转面公式

U2=U11;

L2=L11-d1;

%实际光学成像公式第二个面

I2=asin((L2-r2)*sin(U2)/r2);

I22=asin(n2*sin(I2)/n22);

U22=U2+I2-I22;

L22=r2+r2*sin(I22)/sin(U22);

%入射光线与第一个透镜交点的纵坐标，坐标原点选在第一个透镜的顶点处

h3=L1*tan(U1);

k4=-tan(U1);%入射光线的斜率

b4=100*k4;%因为入射光线经过（-100,0）点

x01=(h3-b4)/k4;%入射光线与第一个透镜交点的横坐标

x4=-100:0.01:x01;

y4=k4*x4+b4;

hold on

plot(x4,y4,'r');%输出入射光线

k6=-tan(U22);

b6=-k6*(d1+L22);%因为第二次折射后出射光线经过（d1+L22,0）点

h4=h3-d1*tan(U11);%第一次折射后入射到第二个透镜的纵坐标x02=(h4-b6)/k6;%第一次折射后入射到第二个透镜的横坐标

k5=(h4-h3)/(x02-x01);%第一次折射后光线的斜率

b5=h4-k5*x02;%因为第一次折射后光线经过（x02，h4）点

x5=x01:0.01:x02

y5=k5*x5+b5;

hold on

plot(x5,y5,'r');%输出第一次折射在两个透镜中的光线x6=x02:0.01:30;%选在30是为了将输出图形看得更清晰些

x6=x02:0.01:30;

y6=k6*x6+b6;

hold on

plot(x6,y6,'r');%输出经过第二个透镜后的输出光线

%球差

m=(L22+d1)-(l22+d1);

end

clear all

n=1.5062;%K9玻璃的折射率

d=3;%透镜的中心厚度

R=25;%透镜凸面曲率半径

f=R/(n-1);%透镜焦距

R0=1;%入射光束半径

lambda=1.064e-3;%波长

k=2*pi/lambda;

phy=lambda*0.61/R0;%角半径

w0=sqrt(f*lambda/pi);%实际光斑半径

data=w0-f*phy; %误差

z=f;

rmax=3*f*phy;%艾利斑半径

r=linspace(0,rmax,100);%产生从0到rmax之间的100点行矢量将衍射半径100等分

eta=linspace(0,2*pi,100);%将0到2*pi100等分

[rho,theta]=meshgrid(r,eta);%生成绘制3D图形所需的网格数据[x,y]=pol2cart(theta,rho);%衍射斑某点的坐标转换极坐标到直角坐标

r0=linspace(0,R0,100);%将入射光束半径100等分

eta0=linspace(0,2*pi,100);

[rho0,theta0]=meshgrid(r0,eta0);

[x0,y0]=pol2cart(theta0,rho0);

for dx=1:100%都是为了建立网格

for dy=1:100

Ep=-i/(lambda*z)*exp(i*k*z)*exp(i*k*((x-x0(dx,dy)).^2+(y-y0(d x,dy)).^2)/(2*z));

E2(dx,dy)=sum(Ep(:));%积分公式的求和表达end

end

Ie=conj(E2).*E2;%光强表达式

figure(1);

surf(x,y,Ie);

figure(2)

plot(x(50,:),Ie(50,:));

2.（1）夫朗和费矩形孔衍射

clear all;

lamda=500e-9;

a=1e-3;

b=1e-3;

f=1;

m=500;

ym=8000*lamda*f;

ys=linspace(-ym,ym,m)

xs=ys;

n=255;

sinth2=ys./sqrt(ys.^2+f^2);%相当于x/f

sinth1=xs(i)/sqrt(xs(i).^2+f^2);%xs(i）作用每给一个ys 值，要遍历到所有的x值

angleA=pi*a*sinth1/lamda;%相当于书上的alfa=kax/2f k=2*pi/lamda

angleB=pi*b*sinth2./lamda;

B(:,i)=(sin(angleA).^2.*sin(angleB).^2.*5000./(angleA.^2.*angl eB.^2));%光强度公式

end

subplot(1,2,1)

image(xs,ys,B)

colormap(gray(n))

subplot(1,2,2)

plot(B(m/2,:),ys)

（2）夫朗和费圆孔衍射

clear

lam=500e-9

a=1e-3

f=1

ym=5*0.61*lam*f/a;%取爱里光斑半径的5倍

ys=linspace(-ym,ym,m);

xs=ys;

n=200;

for i=1:m

r=xs(i)^2+ys.^2;%相当于r的平方

sinth=sqrt(r./(r+f^2));%角度

fai=2*pi*a*sinth./lam;%fai=k*a*sinth

hh=(2*BESSELJ(1,fai)).^2./fai.^2;%贝塞尔函数b(:,i)=hh.*5000;

end

subplot(1,2,1)

image(xs,ys,b)

colormap(gray(n))

subplot(1,2,2)

b(:,m/2)

plot(ys,b(:,m/2))

（3）夫朗和费单缝衍射

clear all

lam=500e-9;

xm=3*lam*f/a;

nx=50;

x=linspace(-xm,xm,nx);

ny=50;

y=linspace(0,a,ny);

for i=1:ny

sinphi=x/f;%角

af=(pi*a*sin(sinphi))/lam;

I(i,:)=5*(sin(af)./af).^2; end

N=255;%确定灰度等级

Br=(I/max(I(1,:)))*N; subplot(1,2,1)

image(x,y,Br);

colormap(gray(N));%颜色subplot(1,2,2)

plot(x,I(1,:));

（4）夫朗和费多缝衍射

clear all;

lamda=500e-9; %波长

N=2; %缝数，可以随意更改变换

a=2e-4;f=5;d=5*a;

ym=2*lamda*f/a;%选择坐标范围

xs=ym;

n=1001;

ys=linspace(-ym,ym,n);

for i=1:n

sinphi=ys(i)/f;

alpha=pi*a*sinphi/lamda;

fai=2*pi*d*sinphi/lamda;

I1=(sin(alpha)./alpha).^2;%单缝衍射因子

B(i,:)=I1*(sin(N*fai/2)./sin(fai/2)).^2;%多缝衍射光强的计算公式B1=B/max(B);%归一化光强

end

NC=256; %确定灰度的等级

Br=(B/max(B))*NC;

subplot(1,2,1)

image(xs,ys,Br);

colormap(gray(NC)); %色调处理

subplot(1,2,2)

plot(B1,ys,'k');

五、画出仿真图形①

②