gusucode.com > 矩孔夫琅和费衍射仿真实现matlab源码 > Untitled.m
focallength=1000; lambda=500; a=1.0; b=1.0; resolution=64; center=(resolution)/2; A=zeros(resolution,resolution); for i=1:1:resolution for j=1:1:resolution if abs(i-center)<a*10/2 & abs(j-center)<b*10/2 A(j,i)=255; end end end E=ones(resolution,resolution); k=2*pi*10000/focallength/lambda; imag=sqrt(-1); for m=1:1:resolution x=m-center; for n=1:1:resolution y=n-center; C=ones(resolution,resolution); for i=1:1:resolution p=i-center; for j=1:1:resolution q=j-center; C(j,i)=A(j,i)*exp(-imag*k*(x*p+y*q)); end end E(n,m)=sum(C(:)); end end E=abs(E); I=E.^2; I=I.^(1/3); I=I.*255/max(max(I)); L=I; I=I+256; CM=[pink(255).^(2/3);gray(255)]; Colormap(CM); edge=(resolution-1)/20; [X,Y]=meshgrid([-edge:0.1:edge]); x=linspace(-edge,edge,resolution); y=linspace(-edge,edge,resolution); subplot(1,2,1); surf(x,y,L); axis([-edge,edge,-edge,edge,0,255]); caxis([0,511]); subplot(1,2,2); image(x,y,I); axis([-edge,edge,-edge,edge,0,511]); view(2);