d=8e-6;
a=1.55e-6;
b=-20:20;
c=b;[B,C]=meshgrid(b,c);
L=0.03;
i=sqrt(-1);
A=-((pi.*d)/(a.*L)).^2;
M=(b.^2+c.^2);
N=((-i).*((2*pi)/(d.*a)));
D=exp(N.*(b.*0+c.*0));
D1=exp(N.*(b.*0+c.*1));
D2=exp(N.*(b.*1+c.*0));
D3=exp(N.*(b.*1+c.*1));
I=abs(exp(A.*M).*(D+D1+D2+D3)).^2 ;
plot3(b,c,I)