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);