用费涅尔波带片模拟透镜成像.m

用费涅尔波带片模拟透镜成像

%这个程序的入射屏面为V字母透射孔，用费涅尔波带片模拟透镜成像
close all

%*******设置参数********
r=632e-9; %入射光源波长
d0=200000;  %P1面到P2面距离
d1=d0/r*d0/256; %费涅尔波带片后焦距长度
k=2*pi/r;   
dx=1e-8;
du=r*d0/256/dx;  %抽样定理
dg=r*d1/256/du;  %抽样定理
%**********************

A=imread('E:\个人文件\毕业设计\衍射 matlab程序\V.bmp');
A=double(A);
subplot(2,2,1);imshow(A);title('P1物面的V字母入射光'); %P1物面的V字母入射函数
a=-128*dx:dx:127*dx; 
[x,y]=meshgrid(a);
b=-128*du:du:127*du; 
[u,v]=meshgrid(b);
B=1/(i*r*d0)*exp(i*k*d0)*exp(i*k/(2*d0)*(u.^2+v.^2)).*fft2(A.*exp(i*k*(x.^2+y.^2)/(2*d0)));  
B1=fftshift(abs(B))/max(max(abs(B)));
subplot(2,2,2);imshow(B1);title('P2面经费涅尔衍射后的光强分布');   %P2面经费涅尔衍射后的光强分布

r2=49;
a1=2.5*sqrt(r2);
b1=0;
%生成菲涅尔滤波片
L=zeros(256,256);
L(80:180,80:180)=1;
for k=1:56
    p=2.5*sqrt(k*r2);
for x=1:256
for y=1:256
    r1=sqrt((x-128)^2+(y-128)^2);
if r1>=a1&r1<=p
    L(x,y)=b1;
end
end
end
a1=p;
b1=abs(b1-1);
end
C=fftshift(B).*L;
C1=abs(C)/max(max(abs(C)));
subplot(2,2,3);imshow(C1); title('P3面经费涅尔波带片后的光强分布');    %P3面经费涅尔波带片后的光强分布

c=-128*dg:dg:127*dg; 
[g,h]=meshgrid(c);
D=1/(i*r*d1)*exp(i*k*d1)*exp(i*k/(2*d1)*(g.^2+h.^2)).*fft2(C.*exp(i*k*(u.^2+v.^2)/(2*d1)));
subplot(2,2,4);imshow(fftshift(abs(D))/max(max(abs(D))));title('P4像面经费涅尔衍射后的图像');    %P4像面经费涅尔衍射后的图像
figure,imshow(L);title('费涅尔波带片');
figure,imshow(fftshift(abs(D))*3/max(max(abs(D))));title('用费涅尔波带片模拟透镜成像—56环');