Fraunhofer_circ.rar_CIRC_Fraunhofer_夫琅禾费_夫琅禾费衍射

function result=Fraunhofer_circ(lam,N,w,d) %Computation of Fraunhofer diffraction of circular hole.lam is the incident %wavelength, N is the number of sampling, w is the radius of circular hole, %d is the diffraction distance. All the distance's units is mm. %meshgrid function is used. %Add the IFFT for recontruct object image. by Song Hao, 2018.03.31. k=2*pi/lam; L0=1;%the computation width n=1:N; x=-L0/2+L0/N.*(n-1);%define x axis y=x; [yy,xx]=meshgrid(y,x);%形成二维取样坐标 zz=sqrt(xx.^2+yy.^2);%define the r of polar axis Lz=k*w.*zz/d; J1=besselj(1,Lz);%compute the first order bessel function Lz(N/2+1,N/2+1)=1;%define the value of cental is 1 f=2*J1./Lz; f(N/2+1,N/2+1)=1;%令中心值为1，恢复准确理论值 C=pi*w*w/lam/d/1i; C2=(pi*w*w/lam/d)^2; U0=C*exp(1i*k*d)*f; U=U0.*exp(1i*k/2/d*zz.*zz); I=C2*f.*f;%the intensity of Fraunhofer diffraction Imax=max(max(I)); figstr=strcat('圆孔半径=',num2str(w),'mm,重建物平面宽度=',num2str(L0),'mm'); p=10; while p strp=num2str(p); figure(1); imshow(I,[0,Imax/p]); title(strcat('Imax/',strp,'限幅显示的衍射斑强度图像')); xlabel(figstr); p=input('Imax/p,限幅显示，p=''?(按Enter键结束）'); end D=1.22*lam*d/w;%define the diameter of the first circle strD=num2str(D); Ix=I(:,N/2);%取x axis value figure plot(x,Ix); title(strcat('第一环直径=',strD,'mm')); xlabel(figstr); %-----IFFT重建图像 Uf=ifft2(U0,N,N); Uf=ifftshift(Uf); figure imshow(abs(Uf),[]); colormap(gray); xlabel(strcat('初始平面宽度=',num2str(L0),'mm')); title('IFFT重建平面'); result=Imax;