%% MATLAB code for implementing transport of amplitude into phase based on Gerchberg-Saxton algorithm (TAP-GSA)


%% Step1-Define basic parameters for calibration 

N=1200; % Enter matrix size 
x=-N/2:N/2-1; % Define x-co-ordinate axis
y=-N/2:N/2-1; % Define y-co-ordinate axis
pixel=8*10^-6; % Enter pixel size 
lambda=0.63*10^-6; % Enter wavelength value
[X,Y]=meshgrid(x*pixel,y*pixel); % Create a calibrated workspace 
R=sqrt(X.*X+Y.*Y); % Define the radius
Aperture=zeros(N,N); % Create a matrix with zeros 
Aperture(R<N/2*pixel)=1; % Create a aperture
imagesc(Aperture) % Display the aperture

%% Step2-Design multiplexed lens and simulate  intensity distribution

z1=1000; % Enter the object distance value
f1=0.14; % Enter the focal length value of first lens
f2=0.25; % Enter the focal length value of second lens
Q1=exp(1i*(pi/(lambda*z1))*(X.*X+Y.*Y)); % Propagate to distance z1
Lens1=exp(-1i*(pi/(lambda*(f1)))*(X.*X+Y.*Y)); % Lens function 1
Lens2=exp(-1i*(pi/(lambda*(f2)))*(X.*X+Y.*Y)); % Lens function 2
LensF2=Lens1+Lens2.*exp(1i*(2)*2*pi/3); % Create a phase shift to Lens2 and add Lens1 and Lens2 ; % Enter the phase shift value in ( ) either 0 or 1 or 2. Here '2' is entered.
z2=0.178; % Enter the image distance value
Field1=Q1.*LensF2.*Aperture; % Interaction
Field_11=padarray(Field1, [N/2 N/2]); % Pads the array of Field1 with zeros
Q2=exp(1i*(pi/(lambda*z2))*(X1.*X1+Y1.*Y1)); % Propagate to distance z2 
Field41=ifftshift(ifft2(fft2(Field_11).*fft2(Q2))); % Take inverse fourier transform
I=Field41(N/2+1:3*N/2,N/2+1:3*N/2).*conj(Field41(N/2+1:3*N/2,N/2+1:3*N/2));% Calculate the intensity distribution
I3=I/max(max(I)); % Normalize the intensity distribution between 0 and 1
imagesc(I3) % Display the image of intensity distribution

%% Step3-(TAP-GSA)

Inputb=ones(N,N); % Create a matrix with ones
Matrix=angle(LensF2); % Take the Phase of LensF2 
I1=I/max(max(I)); % Normalize the intensity distribution between 0 and 1
% Matrix=rand(N,N)*2*pi; % Create a random matrix and multiply by 2pi
iter=100; % Enter iteration value 
Q1=exp(1i*(pi/(lambda*z1))*(X1.*X1+Y1.*Y1)); % Define Propagator to distance z1
Q2=exp(1i*(pi/(lambda*z2))*(X1.*X1+Y1.*Y1)); % Define Propagator to distance z2
Q3=exp(-1i*(pi/(lambda*z2))*(X1.*X1+Y1.*Y1)); % Define Propagator from distance z2
Q4=exp(-1i*(pi/(lambda*z1))*(X1.*X1+Y1.*Y1)); % Define Propagator from distance z1
for p=1:iter % Start for loop
    p % Show iteration number
    A1=Inputb.*exp(1i*Matrix).*Aperture; %  Create new complex field with phase of LensF2 and uniform intensity
    A11=padarray(A1,[N/2 N/2]).*Q1; % Pads the array of A1 with zeros and multiply with propagator Q1
    A2=ifftshift(ifft2(fft2(A11).*fft2(Q2))); % Take inverse fourier transform
    C1=angle(A2); % Take the phase of C2
    C2=angle(Field41); % Take the phase of Field41
    C1(601:1800,601:1800)=C2(601:1800,601:1800); % Define the DOF range 100% DOF different percentage can be set
    A2=padarray(sqrt(I1),[N/2 N/2]).*exp(1i*(C1)); % Pads the array of sqaure root of I1 with zeros
    Matrix1=angle(Q4.*ifftshift(ifft2(fft2(A2).*fft2(Q3)))); % Take the phase of calculated field
    Matrix=Matrix1(N/2+1:3*N/2,N/2+1:3*N/2); % Calculate matrix 
    imagesc(Matrix) % Display the image
    pause(0.1) % pause to see the evolution of the solution
end % end for loop
Phase=Matrix; % Assign Matrix to Phase
Lens33=exp(1i*Phase); % Final pure phase solution

%% Step4-Create mask

M2=(angle(Lens33)+pi)/(2*pi); % Take phase of of the obtained solution and normalize it
Mask2=padarray(M2,[0 360]); % Pads the array of M2 with zeros
imwrite(Mask2,'Enter file path and specify file format'); % Save as image