[sweet_catt]

Imax=length(x)%Imax=101
Jmax=length(y)%Jmax=101

a=d/2;

for I=1:Imax %initialize circular pupil
for J=1:Jmax
p(I,J)=(sqrt(x(I)^2+y(J)^2) <= a) ;
end
end