doublettgratingxyganyan.m

模拟了双周期光栅的衍射情况

clear;
Xc=37e-6;
a1=2.6e-6;b1=1.1e-6;d1=(a1+b1);D=0.1;
a2=2.3e-6;b2=1e-6;d2=(a2+b2);D=0.1;  N=(round((Xc/d1)*(Xc/d2)));
%为光栅各个参数赋值
nx=667;
nx1=200;
ny=420;
nz=560;
lamda1=12e-6; %光的波长
lamda2=10.8e-6; %光的波长
n=80;
x=-0.1:0.0003:0.1;y=-1:0.003:1; %屏上观察平面的大小
Cy=zeros(nx,nx);
E=zeros(nx,nx);
E1=zeros(nx1,nx1);
F=zeros(nx1,nx1);
G=zeros(nx1,nx1);
a=imread('F:\t.bmp');
grid on;
figure(1);
imshow(a,[]);
afft=fft2(a);
aabs=abs(afft);
aabss=fftshift(aabs);
max1=0;
t1=0;
t2=0;
for i=1:ny
   for j=1:nz
       if aabss(i,j)>=max1
           max1=aabss(i,j);
           t1=i;
           t2=j;
       end;
   end;
end;
max1
t1
t2
for i=1:nx1
   for j=1:nx1
       F(i,j)=aabss(i+(t1-nx1/2-1),j+(t2-nx1/2-1));
   end;
end;
figure(2);
imshow(aabss,[]);
colormap(gray);
colorbar;
figure(3);
plot(aabss);
colormap(gray);
for i=1:n
    for ii=1:n
        QB(i,ii)=aabss(210-n/2+i,280-n/2+ii);
    end;
end;
figure(4);
%meshc(QB)
%maxx1=max(max(QB));
surfl(QB);
colormap(hot);

[X,Y]=meshgrid(x,y); %组成数据网格
sinsita=X./(D^2+X.^2).^0.5; %计算sin(θ)
beita1=pi*d1*sinsita./lamda1; %计算β
beita2=pi*d2*sinsita./lamda2; %计算β
I1=(sin(N.*beita1).^2./(sin(beita1).^2)); %相对光强分布
I2=(sin(N.*beita2).^2./(sin(beita2).^2)); %相对光强分布
for i=1:nx
   for j=1:nx
       Cy(j,i)=I2(i,j);
   end;
end;
for i=1:nx
   for j=1:nx
      E(i,j)=I1(i,j)+Cy(i,j);
   end;
end;
xs=linspace(-0.1,0.1,100);
ys=linspace(-0.5,0.5,100);
N=255;
Br=E/2*N;
figure(5);
image(xs,ys,Br); %绘制干涉图样
colormap(gray(N));
%EE=fftshift(E);
max2=0;
t11=0;
t22=0;
for i=1:nx
   for j=1:nx
       if E(i,j)>=max2
           max2=E(i,j);
           t11=i;
           t22=j;
       end;
   end;
end;
max2
t11
t22
for i=1:nx1
   for j=1:nx1
       E1(i,j)=E(i+(t11-nx1/2-1),j+(t22-nx1/2-1));
   end;
end;
  N=255;
Br1=E/2*N;
figure(6);
mesh(Br1); %绘制干涉图样
colormap(winter);
for i=1:nx1
   for j=1:nx1
      G(i,j)=F(i,j)+E1(i,j);
   end;
end;
 N=255;
Br2=G/2*N;
add=0;
for i=1:nx1
   for j=1:nx1
      add=add+Br2(i,j);
   end;
end;
add
max3=0;
t111=0;
t222=0;
for i=1:nx1
   for j=1:nx1
       if Br2(i,j)>=max3
           max3=Br2(i,j);
           t111=i;
           t222=j;
       end;
   end;
end;
max3;
t111;
t222;
S=size(Br2);
Br2(t111,t222)=add/10;
figure(7);
surfl(Br2); %绘制干衍图样
colormap(hot);
shading interp;
title('周期=2.8,孔宽=2.0，波长＝9')