close all
clear all
home

% MATLAB code to obtain the phase element alpha given by Eq. 5

% SLM specifications
horizontal_resolution=1920;
vertical_resolution=1080;

% generate amplitude and phase patterns

% to load amplitude pattern (it must be a digital file "lenna.bmp" within "patterns" folder)
filename1=['.\patterns\lenna.bmp'];
amplitude=imread(filename1);
amplitude=double(amplitude);

% to load phase pattern (it must be a digital file "girl.bmp" within "patterns" folder)
filename2=['.\patterns\girl.bmp'];
phase=imread(filename2);
phase=double(phase);

% to check spatial resolution of patterns (Logical answer must be 1)
size(amplitude,1)==vertical_resolution
size(amplitude,2)==horizontal_resolution

size(phase,1)==vertical_resolution
size(phase,2)==horizontal_resolution

% to set/check minimum and maximun values of patterns

min(amplitude(:))       % should be >= 0
max(amplitude(:))       % should be <= 255

% to set phase between -pi/2 to pi/2
phase=phase-min(phase(:));
phase=phase/max(phase(:));
phase=(2*phase-1)*pi/2;

% to see the defined amplitude and phase patterns 
figure(1);subplot(1,2,1);imagesc(amplitude);title('amplitude');colorbar;axis off;
pbaspect([1 1 1])
figure(1);subplot(1,2,2);imagesc(phase);title('phase');colorbar;axis off;
pbaspect([1 1 1])
set(gcf, 'pos',[50 50 1200 600]);

% calculus of angles (here, named theta1 and theta2!) from Eqs 2 and 3
theta1=phase+acos(amplitude/max(amplitude(:)));
theta2=phase-acos(amplitude/max(amplitude(:)));

% to see the calculated angles theta1 and theta2 
figure(2);subplot(1,2,1);imagesc(theta1);title('theta1');colorbar;axis off;
pbaspect([1 1 1])
figure(2);subplot(1,2,2);imagesc(theta2);title('theta2');colorbar;axis off;
pbaspect([1 1 1])
set(gcf, 'pos',[50 50 1200 600]);

% to generate pairs checkerboard patterns M1 and M2 with different pixel cells
pixel_cell=2;  % typical values (for SLM resolution 1920x1080) 1,2,3,4,5,6,10,12
Mi=zeros(2);Mi(1,1)=1;Mi(2,2)=1;
M1=repmat(Mi,vertical_resolution/(2*pixel_cell),horizontal_resolution/(2*pixel_cell));
M1=imresize(M1, pixel_cell, 'nearest');
M2=not(M1);

% to see the calculated checkerboard patterns M1 and M2 
figure(3);subplot(1,2,1);imagesc(M1);title('M1');colorbar;axis off;
pbaspect([1 1 1])
figure(3);subplot(1,2,2);imagesc(M2);title('M2');colorbar;axis off;
pbaspect([1 1 1])
set(gcf, 'pos',[50 50 1200 600]);

% generation of the phase element alpha from Eq. 5
alpha=M1.*theta1+M2.*theta2;

% to see the calculated phase element alpha
figure(4);imagesc(alpha);title('phase element alpha in radians');colorbar;axis off;
set(gcf, 'pos',[50 50 1200 600]);