function i = Interference(width, height, sizeh,sizev,z)
% parameters explained
% width: is the horizontal pixel pitch in microns
% height: is the vertical pixel pitch in microns
% size is the width=height of the CCD in number of pixels
% z is distance from source to image plane
A0 = 3; %# amplitude of reference wave
A1 = 1; %# amplitude of object wave A0>>A1: A0/A1>=3
lambda = 635 * 10^-9; % wavelength in nanometers
%the linspace was wrong
x=linspace(0,width*sizeh,sizeh); % vector from 0 to width of size 'size'
y=linspace(0,height*sizev,sizev); % vector from 0 to height of size 'size'
[X,Y]=meshgrid(x,y); % matrices of x and y values at each position
s=Super(A0, A1, X-((width*sizeh)/2), Y-((height*sizev)/2), z, lambda); % size-by-size (1024x1024)
r=rand(size); % 1024x1024 matrix of random values on [0 1]
%i=s;
im = zeros(size);
im(r<(s/(A0+A1))) = 1; %# do this all at once instead of pixel-by-pixel
i=im;
end % end of function Interference
% Super is now vectorized, so you can give it a matrix of values for x and y
function S = Super(refamp,objamp,x,y,a,lambda)
r1 = sqrt(a.*a+x.*x+y.*y); % dot notation: multiply element-wise
S = refamp+(objamp*cos(2*pi*r1/(lambda)));
end % end of function Super
usage of the function
width = 2.8 * 10^-6;
height= 2.8 * 10^-6; %pixel size
% sizeh = 16; %image size in pixels
% sizev = 16;
sizeh = 1600; %image size in pixels
sizev = 1200;
int_z = 100*10^-3; % z dist in m
% xes1 = 100;
%xes2 = ;
int_im = Interference(width,height,sizeh, sizev,int_z);
int_im = int_im/max(max(int_im)); % normalize
int_im = (int_im-0.5)*2; % enhance visualy
% display the image
figure
imshow(int_im,[])