diff --git a/Chapter_3_Eight_level_phase_grating.m b/Chapter_3_Eight_level_phase_grating.m
new file mode 100644
index 0000000..a6152c7
--- /dev/null
+++ b/Chapter_3_Eight_level_phase_grating.m
@@ -0,0 +1,28 @@
+%% 4-level 1D grating
+clear %Clear all memory
+
+% Defining grating parameters
+N=500; % Matrix size
+P=100; % Grating period
+A1=ones(P,N); % Size of fundamental building block of grating
+g=8; % Number of phase levels
+delphase= 2*pi/g; %Phase step size
+
+% Constructing one n-level section of the phase grating
+sub =round(P/g)-1;
+for count = 1:g;
+ A1((count -1)*sub+1:count*sub,:)=exp(1i*(count-1)*delphase);
+end
+
+%Constructing the full grating
+A2=repmat(A1,N/P,1);
+
+%Observation of the diffraction pattern
+E=fftshift(fft2(A2)); %fftshift to re-order the terms to the natural order
+IN=(abs(E)/(N*N)).*(abs(E)/(N*N)); %Normalize the intensity values
+figure(1)
+colormap(gray);
+imagesc(angle(A2))
+figure(2)
+colormap(gray);
+imagesc(IN);
diff --git a/Chapter_3_Four_level_axicon.m b/Chapter_3_Four_level_axicon.m
new file mode 100644
index 0000000..37be4e4
--- /dev/null
+++ b/Chapter_3_Four_level_axicon.m
@@ -0,0 +1,32 @@
+% 4-level axicon
+ clear; %Clear all memory
+
+%Defining axicon parameters
+N=480; % Matrix size
+A=zeros(N,N);
+P=40; % Grating period
+g=4; %Define the number of phase levels in each period
+w=P/g; %Define the width of each phase levels
+delphase=2*pi/g;
+%Constructing the 8-level axicon
+ x=1:N;
+ y=1:N;
+ [X,Y]=meshgrid(x,y);
+ r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+for n=1:g;
+ A(rem(r+(n-2)*w,P)N/2)=0;
+
+
+%Observation of the diffraction pattern
+E=fftshift(fft2(A)); %fftshift to re-order the terms to the natural order
+I=abs(E).*abs(E);
+figure(1)
+colormap(gray);
+imagesc(angle(A))
+figure(2)
+colormap(gray);
+imagesc(I);
+
diff --git a/Chapter_3_Sixteen_level_axicon.m b/Chapter_3_Sixteen_level_axicon.m
new file mode 100644
index 0000000..32b489d
--- /dev/null
+++ b/Chapter_3_Sixteen_level_axicon.m
@@ -0,0 +1,32 @@
+%% 4-level axicon
+ clear; %Clear all memory
+
+%Defining axicon parameters
+N=480; % Matrix size
+A=zeros(N,N);
+P=80; % Grating period
+g=16; %Define the number of phase levels in each period
+w=P/g; %Define the width of each phase levels
+delphase=2*pi/g;
+%Constructing the 8-level axicon
+ x=1:N;
+ y=1:N;
+ [X,Y]=meshgrid(x,y);
+ r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+for n=1:g;
+ A(rem(r+(n-2)*w,P)
N/2)=0;
+
+
+%Observation of the diffraction pattern
+E=fftshift(fft2(A)); %fftshift to re-order the terms to the natural order
+I=abs(E).*abs(E);
+figure(1)
+colormap(gray);
+imagesc(angle(A))
+figure(2)
+colormap(gray);
+imagesc(I);
+
diff --git a/Chapter_6_Binary_Axicon_Binary_2D_phase_grating.m b/Chapter_6_Binary_Axicon_Binary_2D_phase_grating.m
new file mode 100644
index 0000000..66d382d
--- /dev/null
+++ b/Chapter_6_Binary_Axicon_Binary_2D_phase_grating.m
@@ -0,0 +1,34 @@
+%Multifunctional DOE – Blazed FZP and 2-d grating
+clear;%Clear all memory
+
+%Define grating and FZP parameters
+N=500;%Define Matrix size
+f=10000;
+lambda=0.632;
+Pr=50;
+Px=25;
+Py=25;
+FFr=0.5;
+FFy=0.5;
+FFx=0.5;
+A1=zeros(N,N);%Define the matrices assigning ones to all pixels
+A2=zeros(N,N);
+A3=zeros(N,N);
+
+%Grating and Axicon construction and multiplexing
+ x=1:N;
+ y=1:N;
+ [X,Y]=meshgrid(x,y);
+ r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+ A1(rem(r,Pr)150)=0;
+
+%Observation of diffraction pattern
+C=fftshift(fft2(B));
+D=abs(C).*abs(C);
+colormap(gray)
+imagesc(D)
diff --git a/Chapter_6_Blazed_Axicon_Binary_2D_phase_grating.m b/Chapter_6_Blazed_Axicon_Binary_2D_phase_grating.m
new file mode 100644
index 0000000..1b412da
--- /dev/null
+++ b/Chapter_6_Blazed_Axicon_Binary_2D_phase_grating.m
@@ -0,0 +1,35 @@
+%Multifunctional DOE – Blazed FZP and 2-d grating
+clear;%Clear all memory
+
+%Define grating and FZP parameters
+N=500;%Define Matrix size
+f=10000;
+lambda=0.632;
+Pr=50;
+Px=25;
+Py=25;
+FFr=0.5;
+FFy=0.5;
+FFx=0.5;
+A1=zeros(N,N);%Define the matrices assigning ones to all pixels
+A2=zeros(N,N);
+A3=zeros(N,N);
+
+%Grating and Axicon construction and multiplexing
+ x=1:N;
+ y=1:N;
+ [X,Y]=meshgrid(x,y);
+ r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+ A1=(rem(r,Pr))*(2*pi)/Pr;
+ A2(rem(Y,Py)150)=0;
+
+%Observation of diffraction pattern
+C=fftshift(fft2(B));
+D=abs(C).*abs(C);
+colormap(gray)
+imagesc(D)
diff --git a/Chapter_6_Blazed_FZP_Binary_2D_phase_grating.m b/Chapter_6_Blazed_FZP_Binary_2D_phase_grating.m
new file mode 100644
index 0000000..d77eabc
--- /dev/null
+++ b/Chapter_6_Blazed_FZP_Binary_2D_phase_grating.m
@@ -0,0 +1,30 @@
+%Multifunctional DOE – Blazed FZP and 2-d grating
+clear;%Clear all memory
+
+%Define grating and FZP parameters
+N=500;%Define Matrix size
+f=10000;
+lambda=0.632;
+P=100;
+FFy=0.5;
+FFx=0.5;
+A1=zeros(N,N);%Define the matrices assigning ones to all pixels
+A2=zeros(N,N);
+A3=zeros(N,N);
+
+%Grating and FZP construction and multiplexing
+ x=1:N;
+ y=1:N;
+ [X,Y]=meshgrid(x,y);
+ r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+ A1=(f-sqrt(f*f+r.*r))*(2*pi)/(0.632);
+ A1=rem(A1,2*pi);
+ A2(rem(Y,P)100)=0;
+ E=fftshift(fft2(Grating_axicon));
+ I2=(abs(E).*abs(E));
+ colormap(gray);%%Display result
+ imagesc(I2);
+
diff --git a/Chapter_7_Helical_axicon.m b/Chapter_7_Helical_axicon.m
new file mode 100644
index 0000000..783b711
--- /dev/null
+++ b/Chapter_7_Helical_axicon.m
@@ -0,0 +1,27 @@
+%Forked grating
+%Define parameters
+N=500;%%Define size of the matrix
+Angle=1;%Define angle of the second plane wave
+V=0.5;%%Visibility controller
+lambda=0.632*1e-6;%Define wavelength
+L=5;
+
+%Create sampled space
+ del=1*1e-6;%sampling
+ x=-N/2:N/2-1;
+ y=-N/2:N/2-1;
+ [X,Y]=meshgrid(x*del,y*del);
+ r=sqrt(X.^2+Y.^2);
+ A=V*exp(1i*L*(atan2(Y,X)));
+ B=V*exp(1i*(r/lambda)*tand(Angle)*2*pi);
+ D=A+B;%Interference of the object and reference wave
+
+%%Intensity profile
+ I=abs(D).*abs(D);
+ I1=im2bw(I);
+ Heli_axicon=exp(1i*pi*I1);
+ Heli_axicon(r>100*1e-6)=0;
+ E=fftshift(fft2(Heli_axicon));
+ I2=(abs(E).*abs(E));
+ colormap(gray);%%Display result
+ imagesc(I2);
\ No newline at end of file
diff --git a/E_5_2.m b/E_5_2.m
new file mode 100644
index 0000000..69546d5
--- /dev/null
+++ b/E_5_2.m
@@ -0,0 +1,27 @@
+%%program to calculate back transverse focal aberration
+u=5000;
+v=30000;
+l=0.632;
+na=1;
+ng=1.5;
+for n=1:1000;
+ a(n)=n*n*l*l+2*n*l*(u+v)+2*u*v;
+ r(n)=sqrt(((a(n)*a(n))-4*u*u*v*v)/(4*(a(n)+u*u+v*v)));
+ if r(n)>=500 && r(n)<=700;
+ t=3000;
+ else
+ t=1000;
+ end
+ theta(n)=atan(r(n)/u);
+ rr(n)=r(n)*((u-t)/u);
+ r1(n)=rr(n)+t*tan(asin(na/ng*(sin(atan(r(n)/u)))));
+ u1(n)=u*(r1(n)/r(n));
+ a1(n)=n*n*l*l+2*n*l*(u1(n)+v)+2*u1(n)*v;
+ rr(n)=sqrt(((a1(n)*a1(n))-4*u1(n)*u1(n)*v*v)/(4*(a1(n)+u1(n)*u1(n)+v*v)));
+ b(n)=(4*n*n*l*l+8*n*l*u1(n)-4*rr(n)*rr(n));
+ c(n)=(4*n^3*l^3+12*n*n*l*l*u1(n)+8*n*l*u1(n)*u1(n)-8*rr(n)*rr(n)*n*l-8*rr(n)*rr(n)*u1(n));
+ d(n)=(n^4*l^4+4*n*n*l*l*u1(n)*u1(n)+4*n^3*l^3*u1(n)-4*rr(n)*rr(n)*u1(n)*u1(n)-4*rr(n)*rr(n)*n*n*l*l-8*rr(n)*rr(n)*n*l*u1(n));
+ v2(n)=(-c(n)+sqrt(c(n)*c(n)-4*b(n)*d(n)))/(2*b(n));
+ the1(n)=atan(rr(n)/v2(n));
+end
+plot(the1)
diff --git a/E_5_3.m b/E_5_3.m
new file mode 100644
index 0000000..dca59a5
--- /dev/null
+++ b/E_5_3.m
@@ -0,0 +1,24 @@
+%%program to calculate back transverse focal aberration
+u=5000;
+v=30000;
+l=0.632;
+na=1;
+ng=1.5;
+for n=1:1000;
+ a(n)=n*n*l*l+2*n*l*(u+v)+2*u*v;
+ r(n)=sqrt(((a(n)*a(n))-4*u*u*v*v)/(4*(a(n)+u*u+v*v)));
+ t(n)=1000+(r(n)/1000)*1000;
+ theta(n)=atan(r(n)/u);
+ rr(n)=r(n)*((u-t(n))/u);
+ r1(n)=rr(n)+t(n)*tan(asin(na/ng*(sin(atan(r(n)/u)))));
+ u1(n)=u*(r1(n)/r(n));
+ a1(n)=n*n*l*l+2*n*l*(u1(n)+v)+2*u1(n)*v;
+ rr(n)=sqrt(((a1(n)*a1(n))-4*u1(n)*u1(n)*v*v)/(4*(a1(n)+u1(n)*u1(n)+v*v)));
+ b(n)=(4*n*n*l*l+8*n*l*u1(n)-4*r(n)*r(n));
+ c(n)=(4*n^3*l^3+12*n*n*l*l*u1(n)+8*n*l*u1(n)*u1(n)-8*r(n)*r(n)*n*l-8*r(n)*r(n)*u1(n));
+ d(n)=(n^4*l^4+4*n*n*l*l*u1(n)*u1(n)+4*n^3*l^3*u1(n)-4*r(n)*r(n)*u1(n)*u1(n)-4*r(n)*r(n)*n*n*l*l-8*r(n)*r(n)*n*l*u1(n));
+ v2(n)=(-c(n)+sqrt(c(n)*c(n)-4*b(n)*d(n)))/(2*b(n));
+ the1(n)=atan(rr(n)/v2(n));
+end
+plot(rr)
+
diff --git a/E_6_3_Spiral_phase_plate_binary_axicon.m b/E_6_3_Spiral_phase_plate_binary_axicon.m
new file mode 100644
index 0000000..480aa53
--- /dev/null
+++ b/E_6_3_Spiral_phase_plate_binary_axicon.m
@@ -0,0 +1,37 @@
+%Spiral phase plate + axicon
+%Clear all memory
+clear;
+%Define Matrix size
+N=500;
+%Define Matrices by assigning 0 or 1 to all elements
+A1=zeros(N,N);
+A2=ones(N,N);
+r=zeros(N,N);
+%Define focal length and wavelength (in micrometers)
+P=10;
+M=(N/P)*0.5;
+L=10;
+r1=zeros(M,M);
+%Calculate the widths of the grating lines
+for n=1:M;
+ r1(n)=n*P;
+end
+%Scan element by element using two for loops
+%Define pixels within grating lines with pi
+for n=1:2:M;
+for p=1:N;
+ for q=1:N;
+ r(p,q)=sqrt((p-N/2)*(p-N/2)+(q-N/2)*(q-N/2));
+ if r(p,q)<100;
+ if r(p,q) > r1(n) && r(p,q) < r1(n+1);
+ A1(p,q)=exp(1i*L*(atan2((q-N/2),(p-N/2))));
+ end
+ end
+ end
+end
+end
+%Observing the diffraction pattern
+E=fftshift(fft2(A1));
+I=abs(E).*abs(E);
+colormap(gray)
+imagesc(I)
diff --git a/E_6_4_Ring_lens_2d_grating.m b/E_6_4_Ring_lens_2d_grating.m
new file mode 100644
index 0000000..6413ba3
--- /dev/null
+++ b/E_6_4_Ring_lens_2d_grating.m
@@ -0,0 +1,57 @@
+%2D Amplitude grating
+%Clear all memory
+clear;
+%Define Matrix size
+N=2000;
+M=50;
+%Define a Matrix by assigning 0 to all elements
+A1=zeros(N,N);
+A2=zeros(N,N);
+A3=zeros(N,N);
+r=zeros(N,N);
+%Define the period of the grating
+f=30000;
+Px=50;
+Py=50;
+FFx=0.5;
+FFy=0.5;
+R=100;
+lambda=0.633;
+%Define fill factors for x and y periodicity
+FFr=0.5;
+FFx=0.5;
+FFy=0.5;
+%Scan element by element using two for loops
+%Use reminder function 'rem' to construct the grating
+x=1:N;
+y=1:N;
+[X,Y]=meshgrid(x,y);
+r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+A1(rem(Y,Py) r(n)
+ r2(n)=R-r(n);
+ end
+end
+for n=1:2:M;
+for p=1:N;
+ for q=1:N;
+ r(p,q)=sqrt((p-N/2)*(p-N/2)+(q-N/2)*(q-N/2));
+ if r(p,q) r1(n) && r(p,q) < r1(n+1);
+ A3(p,q)=1;
+ end
+ end
+ end
+end
+end
+%XOR operation between A1 and A2
+A=exp(1i*pi*xor(A3,xor(A1,A2)));
+A(r>N/2)=0;
+%Observation of the DOE
+colormap(gray)
+imagesc(angle(A));
diff --git a/Norm_outputA.m b/Norm_outputA.m
new file mode 100644
index 0000000..1b6569e
--- /dev/null
+++ b/Norm_outputA.m
@@ -0,0 +1,13 @@
+function Norm_outputA(x, N)
+% This function generates the output intensity pattern with intensity normalized for amplitude DOEs
+% function arguments (required when calling the function) are
+% DOE amplitude (x) and matrix size (N)
+
+E=fftshift(fft2(x)); %fftshift is used to re-order the terms in their natural order
+IN=(abs(E)/(N*N)).*(abs(E)/(N*N)); % Calculating intensity
+figure (1)
+colormap(gray); %colormap(gray) is used to display greyscale image
+imagesc(x);% imagesc is used to display a high constrast image
+figure (2)
+colormap(gray);
+imagesc(IN);
diff --git a/Norm_outputP.m b/Norm_outputP.m
new file mode 100644
index 0000000..324d105
--- /dev/null
+++ b/Norm_outputP.m
@@ -0,0 +1,13 @@
+function Norm_outputP(x, N)
+% This function generates the output intensity pattern with intensity normalized for phase DOEs
+% function arguments (required when calling the function) are
+% DOE amplitude (x) and matrix size (N)
+
+E=fftshift(fft2(x)); %fftshift is used to re-order the terms in their natural order
+IN=(abs(E)/(N*N)).*(abs(E)/(N*N)); % Calculating intensity
+figure (1)
+colormap(gray); %colormap(gray) is used to display greyscale image
+imagesc(angle(x));% imagesc is used to display a high constrast image
+figure (2)
+colormap(gray);
+imagesc(IN);
diff --git a/Table_2_10_Circular_FZP.m b/Table_2_10_Circular_FZP.m
new file mode 100644
index 0000000..bf555d6
--- /dev/null
+++ b/Table_2_10_Circular_FZP.m
@@ -0,0 +1,35 @@
+%%Circular FZP%%
+clear; %Clear all memory
+
+% Defining Grating Parameters
+N=500; %Define Matrix size
+M=32;
+r1=zeros(M);
+A=zeros(N,N);
+r=zeros(N,N);
+f=3000;
+lambda=0.632;
+
+% Constructing the FZP
+for n=1:M; %Calculate the width of the grating lines
+ r1(n)=sqrt(n*f*lambda);
+end
+for n=1:2:M;
+ for p=1:N;
+ for q=1:N;
+ r(p,q)=sqrt((p-N/2)*(p-N/2)+(q-N/2)*(q-N/2));
+ if r(p,q) > r1(n) && r(p,q) < r1(n+1);
+ A(p,q)=exp(1i*pi);
+ end
+ end
+ end
+end
+%Observing the grating output in the far-field
+E=fftshift(fft2(A)); %fftshift is used to re-order the terms in their natural order
+IN=(abs(E)/(N*N)).*(abs(E)/(N*N)); % Calculating intensity
+figure(1)
+colormap(gray);
+imagesc(angle(A))
+figure(2)
+colormap(gray);
+imagesc(IN);
\ No newline at end of file
diff --git a/Table_2_1_Amplitude_1D_grating.m b/Table_2_1_Amplitude_1D_grating.m
new file mode 100644
index 0000000..90eb251
--- /dev/null
+++ b/Table_2_1_Amplitude_1D_grating.m
@@ -0,0 +1,32 @@
+%%1D Amplitude grating%%
+clear %Clear all memory
+
+% Defining Grating Parameters
+ N=500; %Define Matrix size
+ A=zeros(1,N); %Define a row Matrix by assigning 0 to all pixels
+ P=100; %Define the period of the grating
+ FF=0.5; %Define fill factor
+
+% Constructing the Grating
+for q=1:N;
+ if rem(q,P)N/2-2)=0;
+
+%Observing the grating output in the far-field
+E=fftshift(fft2(A));
+IN=(abs(E)/(N*N)).*(abs(E)/(N*N));
+figure(1)
+colormap(gray);
+imagesc(angle(A))
+figure(2)
+colormap(gray);
+imagesc(IN);
diff --git a/Table_2_9_1_D_FZP.m b/Table_2_9_1_D_FZP.m
new file mode 100644
index 0000000..831b669
--- /dev/null
+++ b/Table_2_9_1_D_FZP.m
@@ -0,0 +1,32 @@
+%%1-d FZP%%
+clear; %Clear all memory
+
+%Defining FZP parameters
+N=500; %Define Matrix sizes
+M=50;%Define the number of grating lines
+A=ones(N,N); %Define Matrices by assigning 1 to all pixels
+x=zeros(M,M);
+f=3000;%Define the focal length of FZP in micrometers
+lambda=0.633;%Define wavelength in micrometers
+
+% Constructing the FZP
+for n=1:M;
+ x(n)=sqrt(n*f*lambda);
+end
+for n=1:2:M;
+ for q=1:N;
+ if abs(q-N/2) > x(n) && abs(q-N/2) < x(n+1);
+ A(:,q)=exp(1i*pi);
+ end
+ end
+end
+
+%Observing the FZP and farfield patern
+E=fftshift(fft2(A)); %fftshift is used to re-order the terms in their natural order
+IN=(abs(E)/(N*N)).*(abs(E)/(N*N)); % Calculating intensity
+figure(1)
+colormap(gray);
+imagesc(angle(A))
+figure(2)
+colormap(gray);
+imagesc(IN);
\ No newline at end of file
diff --git a/Table_3_10_IFTA.m b/Table_3_10_IFTA.m
new file mode 100644
index 0000000..8e03f39
--- /dev/null
+++ b/Table_3_10_IFTA.m
@@ -0,0 +1,35 @@
+%% Iterative Fourier Transform Algorithm%
+clear; % Clear all memory
+
+% Loading the target image
+N=500;% Matrix size
+target=imread(‘File location’); % Read image from file
+m=size(target,1); % Size of the image
+scale=N/m; % Estimate the necessary scaling factor
+target=imresize(target,scale); % Resize image to the matrix size
+target=double(target); % Convert symbolic object to a numeric object
+target=target/(max(max(target))); % Normalize matrix
+
+% Defining DOE phase
+DOE=2*pi*rand(N,N); % Generate a N x N matrix of random phase between 0 and 2?
+s=5;
+
+% IFTA algorithm
+for t=1:s; %Iterate to calculate the phase value
+ DOEphase=exp(1i*DOE);
+ % Forward iteration
+ iterf=fft2(DOEphase);
+ intf=abs(iterf);
+ angf=angle(iterf);
+ A=target.*exp(1i*angf);
+ % Backward iteration
+ iterb=ifft2(A);
+ angb=angle(iterb);
+ DOE=angb;
+ error=target-intf/max(max(intf)); %Calculate error
+ E=sum(sum(abs(error)))/(N*N);
+ if E<0.05;
+ iteration=t;
+ break
+ end
+end
diff --git a/Table_3_11_Mesh_technique.m b/Table_3_11_Mesh_technique.m
new file mode 100644
index 0000000..1c0b5b1
--- /dev/null
+++ b/Table_3_11_Mesh_technique.m
@@ -0,0 +1,82 @@
+%program creates a mesh of zones of equal energy for a Gaussian beam
+%mesh is rotated away from axes to avoid zeros
+clear all
+close all
+% all units in mm
+
+global I J R xh yh sigma
+sigma=0.8; % Gaussian spot size
+R=1; % radius of beam
+%mesh co-ordinates
+I=20; %related to rows and therefore also to y coords
+J=20; %related to columns and therefore also to x coords
+
+%N=(I*J); %total number of zones
+TP=pi*sigma*sigma/2; %total power
+P=(1-exp(-2*R*R/(sigma*sigma)))*TP; %total power in radius R
+rings=I/2; %no of rings
+Pring=P/rings; %Power per ring
+
+%calculating radii of rings
+n=1:rings;
+r=sigma*sqrt(0.5*log(1./(1-n.*Pring)));
+
+%%add very small value for first ring (close to zero)
+rad(1)=1e-10;
+
+i=1:n(end);
+rad(i+1)=r(i);
+r=rad;
+
+% number of zones in a ring-decided by angles chosen
+angle_step=pi/I;
+angle_shift=pi/100;
+
+count1=0;
+for theta=pi-angle_shift:-angle_step:-angle_shift,
+count1=count1+1;
+ if theta==(pi-angle_shift)
+ xh=[r.*cos(theta)];
+ yh=[r.*sin(theta)];
+ else
+ xh=[xh;r.*cos(theta)];
+ yh=[yh;r.*sin(theta)];
+end
+end
+
+xh=flipud(xh');
+yh=flipud(yh');
+
+theta=(pi-angle_shift)+angle_step;
+for count2=1:count1,
+ if theta==((pi-angle_shift)+angle_step)
+ xh1=[r.*cos(theta)];
+ yh1=[r.*sin(theta)];
+ else
+ xh1=[xh1;r.*cos(theta)];
+ yh1=[yh1;r.*sin(theta)];
+ end
+ theta=theta+angle_step;
+end
+
+xh1=(xh1');
+yh1=(yh1');
+
+xh=[xh(1:end-1,:);xh1];
+yh=[yh(1:end-1,:);yh1];
+
+figure(1)
+plot(xh,yh,'k.')
+title('Locations of nodes of a Gaussian input')
+xlabel('X-direction (arb. units)')
+ylabel('Y-direction (arb. units)')
+
+figure(2)
+plot(xh',yh')
+title('Annular rings of a Gaussian input')
+xlabel('X-direction (arb. units)')
+ylabel('Y-direction (arb. units)')
+
+% Saves mesh details and variables in a file called GaussianMesh
+% for use by further programs
+save GaussianMesh I J R sigma xh yh
diff --git a/Table_3_12_Flat_top_intensity.m b/Table_3_12_Flat_top_intensity.m
new file mode 100644
index 0000000..0b420fd
--- /dev/null
+++ b/Table_3_12_Flat_top_intensity.m
@@ -0,0 +1,31 @@
+% I, J will be available since defined as global variables
+% if memory was cleared, their values must be entered
+
+%creation of image plane co-ordinates for a square
+S=2; %size of square at image plane is 2S
+stepx=2*S/J;
+stepy=2*S/I;
+X=-S:stepx:S;
+Y=-S:stepy:S;
+[xsq,ysq]=meshgrid(X,Y);
+xsq=flipud(xsq);
+ysq=flipud(ysq);
+figure(2)
+plot(xsq,ysq,'b.')
+xlabel('X-direction (arb. units)')
+ylabel('Y-direction (arb. units)')
+
+clear X Y stepx stepy
+
+bits=512;
+stepx=2*S/bits;
+
+indexx=-S:stepx:S; %regular spaced points
+indexx=indexx(1:end-1); %even number of points
+%indexy=-Sy:stepy:Sy; ; %regular spaced points
+%indexy=indexy(1:end-1); %even number of points
+[XS,YS] = meshgrid(indexx,indexx); %creates grid of points with limits given by index
+
+% Saves mesh details and variables in a file called square
+% for use by further programs
+save square xsq ysq XS YS
diff --git a/Table_3_13_Eikonal_DOE_phase.m b/Table_3_13_Eikonal_DOE_phase.m
new file mode 100644
index 0000000..0702020
--- /dev/null
+++ b/Table_3_13_Eikonal_DOE_phase.m
@@ -0,0 +1,108 @@
+%This program reads mesh data from the files GaussianMesh and square
+%Phase over the hologram is calculated using the Eikonal technique
+%Phase is calculated at a limited number of irregularly spaced data points
+close all
+clear all
+load GaussianMesh
+load square
+bits=512;
+x=xsq;
+y=ysq;
+
+%inputs
+f=300; %distance between input and output planes in mm
+lambda=0.000633; %wavelength in mm
+
+figure(1)
+plot(xh,yh,'r.')
+title('Mesh over hologram plane')
+
+%Eikonal retrieval
+%delta psi
+ del_psix=(x-xh)./f;
+ del_psiy=(y-yh)./f;
+
+%polynomial calculation to obtain phase of hologram
+D=5; %polynomial of degree D
+M=(D+2)*(D+1)/2;
+
+%calculates k and l values for each m value
+count=0;
+for k=0:D
+ for l=0:D-k
+ count=count+1;
+ m=k*((2*D-k+3)/2)+l+1;
+ indices(count,1)=m;
+ indices(count,2)=k;
+ indices(count,3)=l;
+ end
+end
+
+for n=2:M
+ for m=2:M
+ m1=m-1;
+ n1=n-1;
+ c(n1,m1)=0;
+ b(n1,1)=0;
+ for j1=1:J+1
+ for i1=1:I+1
+ k=indices(m,2);
+ k1=indices(n,2);
+ l=indices(m,3);
+ l1=indices(n,3);
+ c(n1,m1)=c(n1,m1)+(k*k1*(xh(j1,i1)^(k+k1-2))*(yh(j1,i1)^(l+l1)))+(l*l1*(xh(j1,i1)^(k+k1))*(yh(j1,i1)^(l+l1-2)));
+ b(n1,1)=b(n1,1)+del_psix(j1,i1)*k1*(xh(j1,i1)^(k1-1))*(yh(j1,i1)^l1)+del_psiy(j1,i1)*l1*(xh(j1,i1)^k1)*(yh(j1,i1)^(l1-1));
+ end
+ end
+ end
+end
+
+%coefficients
+a=c\b;
+clear del_psix del_psiy x y
+
+% eikonal -psi
+for j1=1:J+1
+ for i1=1:I+1
+ psi(j1,i1)=0;
+ for m=2:M
+ psi(j1,i1)=psi(j1,i1)+a(m-1)*(xh(j1,i1)^indices(m,2))*(yh(j1,i1)^indices(m,3));
+ end
+ end
+end
+clear a b c
+
+k=2*pi/lambda;
+psi=k*psi%phase of hologram
+colormap(gray)
+figure(3)
+surf(xh,yh,psi)
+title('Phase of hologram')
+xlabel('X-cordinate in units of length')
+ylabel('Y-cordinate in units of length')
+zlabel('Phase (radians)')
+%save phase psi xh yh
+step=R*2/bits;
+
+%------------------------------------------------------
+%converting arrays to vectors – required for fit process
+[m,n]=size(xh);
+xh1=xh(:,1);
+yh1=yh(:,1);
+psi1=psi(:,1);
+for i=2:n
+ xh1=[xh1;xh(:,i)];
+ yh1=[yh1;yh(:,i)];
+ psi1=[psi1;psi(:,i)];
+end
+
+%saves data in files that can be retrieved by software that will carry out regression
+% and generate coefficients to help generate phase at equidistant points
+output=[xh1 yh1 psi1];
+save save('testdata.txt', 'output', '-ascii')
+%-----------------------------------------------------
+index=-R:step:R; %regular spaced points
+index=index(1:end-1); %even number of points
+[XH1,YH1] = meshgrid(index,index); %creates grid of points with limits given by index
+ %store regularly spaced points and phase for later use
+save phase_circle XH1 YH1 sigma k f lambda xh1 yh1 psi1
diff --git a/Table_3_14_Grid_phase.m b/Table_3_14_Grid_phase.m
new file mode 100644
index 0000000..1db2be3
--- /dev/null
+++ b/Table_3_14_Grid_phase.m
@@ -0,0 +1,77 @@
+% This program calculates the phase over a grid
+%loads uniformly spaced x and y coordinates
+%calculates coefficients using Matlab fit
+
+clear all
+close all
+
+load phase_circle %contains regularly spaced grid points
+
+choice = input('Enter 1 to generate fit coefficients in matlab and 2 to upload coefficients from other programme: ');
+
+if choice == 1
+ ft = fittype('poly55');
+ dataFit = fit([xh1,yh1],psi1, ft);
+ coeff1=coeffvalues(dataFit);
+ coeffN=coeffnames(ft);
+ R=confint(dataFit);
+
+ c=size(coeff1);
+ psi1=0;
+
+ for count=1:c(2)
+
+ coeffI=coeffN(count);
+ coeffI=coeffI{:};
+ m=str2double(coeffI(2));
+ n=str2double(coeffI(3));
+ psi1=psi1+coeff1(count).*(XH1.^m).*(YH1.^n);
+
+ end
+
+elseif choice ==2
+
+ %coefficients obtained with R=1, sigma=0.8 and I=J=20, S=2 from a data handling software such as excel
+ %SQUARE
+ a0 = 0.106394007
+ a1 = -0.836150168
+ a2 = 0.020604165
+ a3 = 17.61722947
+ a4 = 15.83427723
+ a5 = -7.266911537
+ a6 = -4.752015829
+
+ % Calculating phase at regularly spaced points
+ psi1=a0+a1*XH1+a2*YH1+a3*XH1.^2+a4*YH1.^2+a5*XH1.^4+a6*YH1.^4;
+end
+
+% Output
+figure(1)
+surf(XH1,YH1,psi1)
+title('Phase over DOE')
+ %______________________________________________________%
+% to simulate output intensity
+i=sqrt(-1);
+g=exp(-((XH1).^2+(YH1).^2)/(sigma^2)); %Gaussian beam
+amp1=exp(i*psi1).*g; %Amplitude just after hologram
+sing=find(isnan(amp1)); %locates singularities
+amp1(sing)=zeros(size(sing)); %replaces singularities with zeros
+
+amp1=amp1.*exp(i*(k/(2*f))*(XH1.^2+YH1.^2));
+figure(2)
+colormap(gray)
+imagesc(amp1.*conj(amp1))
+contour(XH1,YH1,amp1.*conj(amp1))
+title('Intensity just after DOE')
+xlabel('X-direction')
+ylabel('Y-direction')
+
+amp2=fft2((amp1));
+amp2=fftshift(amp2);
+
+figure(3)
+colormap(gray)
+contour(abs(amp2)/max(max(abs(amp2))))
+title('Output intensity - Gaussian to square')
+xlabel('X-direction')
+ylabel('Y-direction')
diff --git a/Table_3_1_Four_level_phase_grating.m b/Table_3_1_Four_level_phase_grating.m
new file mode 100644
index 0000000..e00bce9
--- /dev/null
+++ b/Table_3_1_Four_level_phase_grating.m
@@ -0,0 +1,28 @@
+%% 4-level 1D grating
+clear %Clear all memory
+
+% Defining grating parameters
+N=500; % Matrix size
+P=100; % Grating period
+A1=ones(P,N); % Size of fundamental building block of grating
+g=4; % Number of phase levels
+delphase= 2*pi/g; %Phase step size
+
+% Constructing one n-level section of the phase grating
+sub =round(P/g)-1;
+for count = 1:g;
+ A1((count -1)*sub+1:count*sub,:)=exp(1i*(count-1)*delphase);
+end
+
+%Constructing the full grating
+A2=repmat(A1,N/P,1);
+
+%Observation of the diffraction pattern
+E=fftshift(fft2(A2)); %fftshift to re-order the terms to the natural order
+IN=(abs(E)/(N*N)).*(abs(E)/(N*N)); %Normalize the intensity values
+figure(1)
+colormap(gray);
+imagesc(angle(A2))
+figure(2)
+colormap(gray);
+imagesc(IN);
diff --git a/Table_3_2_Blazed_phase_grating.m b/Table_3_2_Blazed_phase_grating.m
new file mode 100644
index 0000000..5609ee4
--- /dev/null
+++ b/Table_3_2_Blazed_phase_grating.m
@@ -0,0 +1,26 @@
+%%Blazed 1D grating
+clear; %Clear all memory
+
+% Defining grating parameters
+N=500; % Matrix size
+P=100; % Grating period
+A1=ones(P,N); % Size of fundamental building block of grating
+
+% Constructing the fundamental block of the blazed grating
+for p=1:P;
+ A1(p,:)=exp(1i*(p/P)*2*pi);
+end
+
+%Constructing the full grating
+A2=repmat(A1,N/P,1);
+
+%Observation of the diffraction pattern
+E=fftshift(fft2(A2)); %fftshift to re-order the terms to the natural order
+IN=(abs(E)/(N*N)).*(abs(E)/(N*N)); %Normalize the intensity values
+figure(1)
+colormap(gray);
+imagesc(angle(A2))
+figure(2)
+colormap(gray);
+imagesc(IN);
+
diff --git a/Table_3_3_Eight_level_axicon.m b/Table_3_3_Eight_level_axicon.m
new file mode 100644
index 0000000..7392f40
--- /dev/null
+++ b/Table_3_3_Eight_level_axicon.m
@@ -0,0 +1,32 @@
+%% 8-level axicon
+ clear; %Clear all memory
+
+%Defining axicon parameters
+N=480; % Matrix size
+A=zeros(N,N);
+P=40; % Grating period
+g=8; %Define the number of phase levels in each period
+w=P/g; %Define the width of each phase levels
+delphase=2*pi/g;
+%Constructing the 8-level axicon
+x=1:N;
+y=1:N;
+[X,Y]=meshgrid(x,y);
+r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+for n=1:g;
+ A(rem(r+(n-2)*w,P)N/2)=0;
+
+
+%Observation of the diffraction pattern
+E=fftshift(fft2(A)); %fftshift to re-order the terms to the natural order
+I=abs(E).*abs(E);
+figure(1)
+colormap(gray);
+imagesc(angle(A))
+figure(2)
+colormap(gray);
+imagesc(I);
+
diff --git a/Table_3_4_Blazed_axicon.m b/Table_3_4_Blazed_axicon.m
new file mode 100644
index 0000000..f951d75
--- /dev/null
+++ b/Table_3_4_Blazed_axicon.m
@@ -0,0 +1,25 @@
+%%Blazed axicon
+clear;%Clear all memory
+%Define axicon parameters
+N=500; % Matrix size
+P=100; % Grating period
+A=ones(N,N); % Set up matrix, assigning 1’s to all pixels
+
+%Constructing the blazed axicon
+ x=1:N;
+ y=1:N;
+ [X,Y]=meshgrid(x,y);
+ r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+ A=exp(1i*(rem(r,P))*(2*pi)/P);
+ A(r>N/2)=0;
+
+%Observation of the diffraction pattern
+E=fftshift(fft2(A)); %fftshift to re-order the terms to the natural order
+I=abs(E).*abs(E);
+figure(1)
+colormap(gray);
+imagesc(angle(A))
+figure(2)
+colormap(gray);
+imagesc(I);
+
diff --git a/Table_3_5_Blazed_FZP.m b/Table_3_5_Blazed_FZP.m
new file mode 100644
index 0000000..974cc63
--- /dev/null
+++ b/Table_3_5_Blazed_FZP.m
@@ -0,0 +1,19 @@
+% Blazed FZP
+clear; %Clear all memory
+
+% Defining FZP parameters
+N=500; % Matrix size
+f=10000; % focal length in microns
+lambda=0.632; % Wavelength in microns
+
+% Constructing the blazed FZP
+ x=1:N;
+ y=1:N;
+ [X,Y]=meshgrid(x,y);
+ r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+ A=exp(1i*(f-sqrt(f*f+r.*r))*(2*pi)/(0.632));
+ A(r>N/2)=0;
+
+%Observing the phase profile
+colormap(gray)
+imagesc(angle(A))
\ No newline at end of file
diff --git a/Table_3_6_Four_level_FZP.m b/Table_3_6_Four_level_FZP.m
new file mode 100644
index 0000000..bd7f50b
--- /dev/null
+++ b/Table_3_6_Four_level_FZP.m
@@ -0,0 +1,33 @@
+%% 4-level FZP
+clear;%Clear all memory
+
+%Define FZP parameters
+N=500; % Matrix size
+f=10000; % Focal length in microns
+lambda=0.632;% Wavelength in microns
+g=4; % Number of phase levels
+delphase=(2*pi)/g;
+% Set up matrices, assigning 0’s or 1’s to all pixels
+C=ones(N,N);
+
+%Constructing the FZP
+x=1:N;
+y=1:N;
+[X,Y]=meshgrid(x,y);
+r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+A=(f-sqrt(f*f+r.*r))*(2*pi)/(0.632);
+B=rem(A,2*pi);
+B(r>N/2)=0;
+for p=1:N;%Construct the n-level FZP
+ for q=1:N;
+ for n1=1:g;
+ if B(p,q)> (-2*pi+(n1-1)*delphase) && B(p,q) <= (-3*pi/2+(n1-1)*delphase);
+ C(p,q)=exp(1i*(-2*pi+(n1-1)*pi/2));
+ end
+ end
+ end
+end
+
+%Observing the phase profile
+colormap(gray)
+imagesc(angle(C))
\ No newline at end of file
diff --git a/Table_3_7_Greyscale_SPP.m b/Table_3_7_Greyscale_SPP.m
new file mode 100644
index 0000000..1da58bc
--- /dev/null
+++ b/Table_3_7_Greyscale_SPP.m
@@ -0,0 +1,21 @@
+%Greyscale SPP
+clear; %Clear all memory
+
+%Defining SPP parameters
+N=500; % Matrix size
+L=1;% Topological charge number
+
+%Constructing the SPP
+ x=1:N;
+ y=1:N;
+ [X,Y]=meshgrid(x,y);
+ theta=atan2((X-N/2),(Y-N/2));
+ r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+ A1=exp(1i*L*(theta));
+ A1(r>30)=0;
+
+%Observation of the far-field pattern
+E=fftshift(fft2(A1)); %Calculate the Fourier transform and rearrange terms
+I=abs(E).*abs(E);
+colormap(gray)
+imagesc(I)
diff --git a/Table_3_8_Multilevel_SPP.m b/Table_3_8_Multilevel_SPP.m
new file mode 100644
index 0000000..a0cbad9
--- /dev/null
+++ b/Table_3_8_Multilevel_SPP.m
@@ -0,0 +1,36 @@
+%Multilevel SPP
+clear;%Clear all memory
+
+%Define SPP parameters
+N=500; % Matrix size
+L=1; % Topological charge number
+g=5; % Number of phase levels
+delphase=(L*2*pi)/g; % Phase increment
+
+%Constructing the SPP
+x=1:N;
+y=1:N;
+[X,Y]=meshgrid(x,y);
+theta=atan2((X-N/2),(Y-N/2));
+r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+A1=L*(theta)+L*pi;
+A1(r>30)=0;
+for n=1:g;
+ for p=1:N;%
+ for q=1:N;
+ if A1(p,q) > (delphase*(n-1)) && A1(p,q) <= (delphase*n)
+ A1(p,q) = exp(1i*delphase*(n-1));
+ end
+ end
+ end
+end
+
+%Observation of the diffraction pattern
+E=fftshift(fft2(A1)); %fftshift to re-order the terms to the natural order
+I=abs(E).*abs(E);
+figure(1)
+colormap(gray);
+imagesc(angle(A1))
+figure(2)
+colormap(gray);
+imagesc(I);
diff --git a/Table_3_9_Blazed_Axilens.m b/Table_3_9_Blazed_Axilens.m
new file mode 100644
index 0000000..4f7dfbb
--- /dev/null
+++ b/Table_3_9_Blazed_Axilens.m
@@ -0,0 +1,26 @@
+%% Blazed axilens
+clear; % Clear all memory
+
+% Define axilens parameters
+% All lengths in microns
+N=500; % Matrix size
+f0=5000;% focal length
+delf=1000;% focal depth
+R=1000;% radius of the element
+lambda=0.632;% wavelength
+
+%Constructing the axilens
+ x=1:N;
+ y=1:N;
+ [X,Y]=meshgrid(x,y);
+ r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+ f=f0+(delf/R)*r; % Focal length equation
+ A=(f-sqrt(f.*f+r.*r))*(2*pi)/(0.632);
+ B=exp(1i*rem(A,2*pi));
+ B(r>N/2)=0;
+
+
+%Observation of phase of axilens
+colormap(gray);
+imagesc(angle(B));
+
diff --git a/Table_4_1_Circular_aperture.m b/Table_4_1_Circular_aperture.m
new file mode 100644
index 0000000..2bf3eaf
--- /dev/null
+++ b/Table_4_1_Circular_aperture.m
@@ -0,0 +1,23 @@
+%Fresnel diffraction of a circular aperture%
+clear; %Clear all memory
+%Defining the parameters
+N=500;% Define the matrix size
+lambda=0.633*10^-6;%Define wavelength in meters
+z=0.02;%Propagation distance = 20 mm
+r=10^-3;%Radius of aperture = 1 mm
+wsamp=10*lambda;%sampling period or width
+%Creating sampled space
+ x=1:N;
+ y=1:N;
+ [X,Y]=meshgrid(x,y);%Sampling
+ Rsamp=sqrt((X-N/2).^2+(Y-N/2).^2).*wsamp;%Define sampled radius
+%Constructing the aperture
+ A=ones(N,N);%Define matrix by assigning zeros to all pixels
+ A(Rsamp>=r)=0;
+% Calculating the Fresnel diffraction pattern
+ PPF=exp(1i*pi/(lambda*z).*Rsamp.*Rsamp); %Calculate the parabolic phase factor
+ A1=A.*PPF;%Multiply the circular aperture function with the parabolic phase factor
+ E=abs(fftshift(fft2(fftshift(A1)))); %Calculate Fourier Transform
+%Observation of the diffraction pattern
+colormap(gray)%Display greyscale image
+imagesc(E)%Display scaled image
diff --git a/Table_4_2_Binary_axicon.m b/Table_4_2_Binary_axicon.m
new file mode 100644
index 0000000..ec2df4c
--- /dev/null
+++ b/Table_4_2_Binary_axicon.m
@@ -0,0 +1,26 @@
+%Fresnel diffraction of binary axicon%
+clear; %Clear all memory
+
+% Defining the parameters
+N=500;% Define the matrix size
+lambda=0.633*10^-6;%Define wavelength in meters
+z=0.025;%Propagation distance = 25mm and 50 m
+P=10^-4;%Radius of aperture = 0.1 mm
+wsamp=10*lambda;%sampling period or width
+%Creating sampled space
+ x=1:N;
+ y=1:N;
+ [X,Y]=meshgrid(x,y);%Sampling
+ Rsamp=sqrt((X-N/2).^2+(Y-N/2).^2).*wsamp;%Define sampled radius
+
+%Constructing the DOE
+ A=ones(N,N);%Define matrix by assigning ones to all pixels
+ A(rem(Rsamp,P) < P/2)=exp(1i*pi);
+
+%Calculating the Fresnel diffraction pattern
+ PPF=exp(1i*pi/(lambda*z).*Rsamp.*Rsamp); %Calculate the parabolic phase factor
+ A1=A.*PPF; %Multiply the circular aperture function with the parabolic phase factor
+ E=abs(fftshift(fft2(fftshift(A1)))); %Calculate Fourier Transform
+ colormap(gray)
+ imagesc(E)
+
diff --git a/Table_4_3_Axilens.m b/Table_4_3_Axilens.m
new file mode 100644
index 0000000..5db9d36
--- /dev/null
+++ b/Table_4_3_Axilens.m
@@ -0,0 +1,37 @@
+% Fresnel diffraction of an axilens
+clear; %Clear all memory
+
+%Defining all parameters and sampling
+N=500; %Matrix size
+lambda=0.633e-6; %Wavelength
+wsamp=1e-6;%sampling period
+ x=-N/2:N/2-1;
+ y=-N/2:N/2-1;
+ [X,Y]=meshgrid(x*wsamp,y*wsamp);
+ Rsamp=sqrt(X.^2+Y.^2);
+ f0=0.005;%Focal length
+ delf=0.003;%Focal depth
+ R=10^-3;%Radius of the axilens
+
+%Design of axilens
+ f=(f0+(delf/R)*sqrt(X.^2+Y.^2));%Focal length calculation
+ FZA=exp(-1i*(pi/(lambda))*((X.^2+Y.^2)./f));%Phase profile of Axilens
+
+%Calculation of Fresnel diffraction pattern
+ m=100;
+ n=1:m;
+ zs2=0.003+(0.005/m).*n; % Propagation distance
+ PPF=zeros(N,N,m);
+ A1=zeros(N,N,m);
+ E=zeros(N,N,m);
+ Field1=zeros(100,m);
+ for counter1=1:1:m;
+ PPF(:,:,counter1)=exp(1i*pi/(lambda*zs2(counter1)).*Rsamp.*Rsamp); % Parabolic phase factor
+ A1(:,:,counter1)=FZA.*PPF(:,:,counter1); %Multiply the axilens function with parabolic phase factor
+ E(:,:,counter1)=abs(fftshift(fft2(fftshift(A1(:,:,counter1))))); %Calculate Fourier Transform
+ % imagesc(E(201:300,201:300,counter1));
+ % pause(1.0)
+ Field1(:,counter1)=E(N/2+1,201:300,counter1); %Accumulate the intensity profile
+ end
+colormap(gray)
+imagesc(Field1)
\ No newline at end of file
diff --git a/Table_4_4_IFTA.m b/Table_4_4_IFTA.m
new file mode 100644
index 0000000..fd61aa8
--- /dev/null
+++ b/Table_4_4_IFTA.m
@@ -0,0 +1,29 @@
+% Fresnel diffraction of a DOE designed by Gerchberg-Saxton algorithm
+% Use table 3.10 to obtain the “DOE phase” and save it as DOE
+% Multiply the “DOE phase” with a lens function
+%Defining all parameters and sampling
+DOE1=exp(1i*DOE);
+wsamp=1e-6;%sampling period
+ x=-N/2:N/2-1;
+ y=-N/2:N/2-1;
+ [X,Y]=meshgrid(x*wsamp,y*wsamp);
+ Rsamp=sqrt(X.^2+Y.^2);
+ lambda=0.633e-6; %Wavelength
+ Q=exp(-1i*(pi/(lambda*d))*(X.^2+Y.^2));
+ DOE2=DOE1.*Q;
+%Calculation of Fresnel diffraction pattern
+ m=1;
+ n=1:m;
+ zs2=0.05; % Propagation distance
+ PPF=zeros(N,N,m);
+ A1=zeros(N,N,m);
+ E=zeros(N,N,m);
+ Field1=zeros(100,m);
+ for counter1=1:1:m;
+ PPF(:,:,counter1)=exp(1i*pi/(lambda*zs2(counter1)).*Rsamp.*Rsamp); % Parabolic phase factor
+ A1(:,:,counter1)=DOE2.*PPF(:,:,counter1); %Multiply the axilens function with the parabolic phase factor
+ E(:,:,counter1)=abs((fft2((A1(:,:,counter1))))); %Calculate Fourier Transform
+ % imagesc(E(:,:,counter1)); title (zs2)
+ % pause(1.0)
+ end
+
diff --git a/Table_4_5_Talbot_effect.m b/Table_4_5_Talbot_effect.m
new file mode 100644
index 0000000..27f8eee
--- /dev/null
+++ b/Table_4_5_Talbot_effect.m
@@ -0,0 +1,50 @@
+%program to generate Talbot images of amplitude grating
+%Output of this file is very dependent on values of r0, P, N, w
+clear all
+
+%constants
+lambda=600e-9; %in m
+k=2*pi/lambda;
+
+% Defining Grating Parameters
+N=2000; %Define Matrix size
+P=75; %Define the period of the grating in pixels
+FF=0.25; %Define fill factor
+period=P*1e-6; %period in distance units
+A=ones(1,N); %Define a Matrix by assigning 1 to all pixels
+
+% Constructing the Grating
+ q=1:N;
+ A(rem(q,P)150)=0;
+
+%Observing the grating output in the far-field
+E=fftshift(fft2(A3));
+I=abs(E).*abs(E);
+figure(1)
+colormap(gray)
+imagesc(angle(A3))
+figure(2)
+colormap(gray)
+imagesc(I)
\ No newline at end of file
diff --git a/Table_6_3_Blazed_FZP_Binary_1D_phase_grating.m b/Table_6_3_Blazed_FZP_Binary_1D_phase_grating.m
new file mode 100644
index 0000000..aa833a0
--- /dev/null
+++ b/Table_6_3_Blazed_FZP_Binary_1D_phase_grating.m
@@ -0,0 +1,28 @@
+%Multifunctional DOE – Blazed FZP and 1-d grating
+clear;%Clear all memory
+
+%Define grating and FZP parameters
+N=500;%Define Matrix size
+f=10000;
+lambda=0.632;
+P=200;
+FFy=0.5;
+A1=zeros(N,N);%Define the matrices assigning ones to all pixels
+A2=zeros(N,N);
+A3=zeros(N,N);
+
+%Grating and FZP construction
+ x=1:N;
+ y=1:N;
+ [X,Y]=meshgrid(x,y);
+ r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+ A1=(f-sqrt(f*f+r.*r))*(2*pi)/(0.632);
+ A2(rem(Y,P)150)=0;
+
+%Observation of diffraction pattern
+E=(fftshift(fft2(B)));
+I=abs(E).*abs(E);
+colormap(gray)
+imagesc(I)
diff --git a/Table_6_5_SPP_FZP.m b/Table_6_5_SPP_FZP.m
new file mode 100644
index 0000000..564b1a2
--- /dev/null
+++ b/Table_6_5_SPP_FZP.m
@@ -0,0 +1,29 @@
+%Multifunctional DOE – Blazed circular grating and blazed 1-d grating
+clear; %Clear all memory
+%Define grating parameters
+N=500; %Define Matrix size
+M=32; %Define number of zones
+A1=zeros(N,N); %Define Matrices by assigning 0 or 1 to all elements
+r1=zeros(M,M);
+r=zeros(N,N);
+f=3000; %Define focal length and wavelength (in micrometers)
+lambda=0.633;
+L=1;%Define topological charge
+% Construction of the spiral FZP
+for n=1:M; %Calculate the widths of the zones
+ r1(n)=sqrt(n*f*lambda);
+end
+for n=1:2:M; %Scan element by element using two for loops
+ for p=1:N;
+ for q=1:N;
+ r(p,q)=sqrt((p-N/2)*(p-N/2)+(q-N/2)*(q-N/2));
+ if r(p,q) > r1(n) && r(p,q) < r1(n+1);
+ A1(p,q)=exp(1i*L*(atan2((q-N/2),(p-N/2))));
+ end
+ end
+ end
+end
+
+%Observation of the DOE
+colormap(gray)
+imagesc(angle(A1))
\ No newline at end of file
diff --git a/Table_7_1_grating_1d.m b/Table_7_1_grating_1d.m
new file mode 100644
index 0000000..d6df693
--- /dev/null
+++ b/Table_7_1_grating_1d.m
@@ -0,0 +1,26 @@
+%%1-d grating
+%Define parameters
+N=500;%%Define size of the matrix
+Angle=1;%Define angle of the second plane wave
+V=0.5;%%Visibility controller
+lambda=0.632*1e-6;%Define wavelength
+
+%Create sampled space
+ del=1*1e-6;%sampling
+ x=-N/2:N/2-1;
+ y=-N/2:N/2-1;
+ [X,Y]=meshgrid(x*del,y*del);
+
+%Simulate interference
+ A=V*exp(1i*(Y/lambda)*tand(Angle)*2*pi);
+ B=V*exp(1i*2*pi);
+ D=A+B;%Interference of the object and reference wave
+ I=abs(D).*abs(D);
+
+%%Construct grating
+I1=im2bw(I);%Binarize the matrix
+Grating=exp(1i*pi*I1);%Generate the phase grating
+
+%Observation of the HOE
+colormap(gray)
+imagesc(angle(Grating))
\ No newline at end of file
diff --git a/Table_7_2_FZP.m b/Table_7_2_FZP.m
new file mode 100644
index 0000000..d9abfca
--- /dev/null
+++ b/Table_7_2_FZP.m
@@ -0,0 +1,30 @@
+%%Fresnel zone plate
+%Define parameters
+N=500;%%Define size of the matrix
+Angle=1;%Define angle of the second plane wave
+V=0.5;%%Visibility controller
+lambda=0.632*1e-6;%Define wavelength
+
+%Create sampled space
+ del=1*1e-6;%sampling
+ x=-N/2:N/2-1;
+ y=-N/2:N/2-1;
+ [X,Y]=meshgrid(x*del,y*del);
+ f=0.01; %Focal length of 1 cm
+
+%Simulate interference
+ A=V*exp(1i*((2*pi)/(lambda))*sqrt(X.^2+Y.^2+f*f));
+ B=V*exp(1i*2*pi);
+ D=A+B;%Interference of the object and reference wave
+ I=abs(D).*abs(D);
+ figure (2)
+ colormap gray
+ imagesc(I)
+
+%%Construct FZP
+ I1=im2bw(I);%Binarize the matrix
+ FZP=exp(1i*pi*I1);%Generate the FZP
+
+%Observation of the HOE
+colormap(gray)
+imagesc(angle(FZP))
diff --git a/Table_7_3_Axicon.m b/Table_7_3_Axicon.m
new file mode 100644
index 0000000..62253a6
--- /dev/null
+++ b/Table_7_3_Axicon.m
@@ -0,0 +1,29 @@
+%Diffractive axicon
+%Define parameters
+N=500;%%Define size of the matrix
+Angle=1;%Define angle of the second plane wave
+V=0.5;%%Visibility controller
+lambda=0.632*1e-6;%Define wavelength
+
+%Create sampled space
+ del=1*1e-6;%sampling
+ x=-N/2:N/2-1;
+ y=-N/2:N/2-1;
+ [X,Y]=meshgrid(x*del,y*del);
+
+%Simulate interference
+ A=V* exp(1i*(sqrt(X.^2+Y.^2)/lambda)*tand(Angle)*2*pi);
+ B=V*exp(1i*2*pi);
+ D=A+B;%Interference of the object and reference wave
+ I=abs(D).*abs(D);
+ figure (2)
+ colormap gray
+ imagesc(I)
+
+%%Construct FZP
+ I1=im2bw(I);%Binarize the matrix
+ Grating=exp(1i*pi*I1);%Generate the Axicon
+
+%Observation of the HOE
+colormap(gray)
+imagesc(angle(Grating))
\ No newline at end of file
diff --git a/Table_7_4_Forked_grating.m b/Table_7_4_Forked_grating.m
new file mode 100644
index 0000000..a8345a2
--- /dev/null
+++ b/Table_7_4_Forked_grating.m
@@ -0,0 +1,30 @@
+%Forked grating
+%Define parameters
+N=500;%%Define size of the matrix
+Angle=1;%Define angle of the second plane wave
+V=0.5;%%Visibility controller
+lambda=0.632*1e-6;%Define wavelength
+L=1;
+
+%Create sampled space
+ del=1*1e-6;%sampling
+ x=-N/2:N/2-1;
+ y=-N/2:N/2-1;
+ [X,Y]=meshgrid(x*del,y*del);
+ r=sqrt(X.^2+Y.^2);
+ A=V*exp(1i*L*(atan2(Y,X)));
+ B=V*exp(1i*(Y/lambda)*tand(Angle)*2*pi);
+ D=A+B;%Interference of the object and reference wave
+
+%%Intensity profile
+ I=abs(D).*abs(D);
+ I1=im2bw(I);
+ Forked_grating=exp(1i*pi*I1);
+ Forked_grating(r>100*1e-6)=0;
+ E=fftshift(fft2(Forked_grating));
+ I2=(abs(E).*abs(E));
+
+%%Display result
+colormap(gray);
+imagesc(I2);
+
diff --git a/Table_7_5_Off_axis_axilens.m b/Table_7_5_Off_axis_axilens.m
new file mode 100644
index 0000000..f42b9d9
--- /dev/null
+++ b/Table_7_5_Off_axis_axilens.m
@@ -0,0 +1,33 @@
+%Off axis axilens
+%Define parameters
+N=500;%%Define size of the matrix
+Angle=1;%Define angle of the second plane wave
+V=0.5;%%Visibility controller
+lambda=0.632*1e-6;%Define wavelength
+L=3;
+f0=0.01;
+delf=0.003;
+
+%Create sampled space
+ del=1*1e-6;%sampling
+ x=-N/2:N/2-1;
+ y=-N/2:N/2-1;
+ [X,Y]=meshgrid(x*del,y*del);
+ f=f0+delf*sqrt(X.^2+Y.^2);
+
+%Simulate interference
+ A=V*exp(1i*((2*pi)/lambda)*(f-sqrt(f.*f-(X.^2+Y.^2))));
+ B=V*exp(1i*(2*pi/lambda)*tand(Angle)*Y);
+ D=A+B;%Interference of the object and reference wave
+
+%Far field diffraction pattern
+ I=abs(D).*abs(D);
+ I1=im2bw(I);
+ Off_axilens=exp(1i*pi*I1);
+ Off_axilens(r>100*1e-6)=0;
+ E=fftshift(fft2(Off_axilens));
+ I2=(abs(E).*abs(E));
+
+%%Display result
+colormap(gray);
+imagesc(I2);
diff --git a/Table_7_6_Accelerating_beams.m b/Table_7_6_Accelerating_beams.m
new file mode 100644
index 0000000..9b7aa6e
--- /dev/null
+++ b/Table_7_6_Accelerating_beams.m
@@ -0,0 +1,29 @@
+%Accelerating beams
+%Define parameters
+N=500;%%Define size of the matrix
+Angle=0.1;%Define angle of the second plane wave
+V=0.5;%%Visibility controller
+lambda=0.632*1e-6;%Define wavelength
+
+%Create sampled space
+ del=1*1e-6;%sampling
+ x=-N/2:N/2-1;
+ y=-N/2:N/2-1;
+ [X,Y]=meshgrid(x*del,y*del);
+ r=sqrt(X.^2+Y.^2);
+ A=V*exp(1i*(r/lambda)*tand(Angle)*2*pi);
+ B=V*exp(1i*(10*sqrt(Y/del+N/2)+10*sqrt(X/del+N/2)));
+ D=A+B;%Interference of the object and reference wave
+
+%%Intensity profile
+ I=abs(D).*abs(D);
+ I1=im2bw(I);
+ AB=exp(1i*pi*I1);
+% AB(r>100*del)=0;
+ E=fftshift(fft2(AB));
+ I2=(abs(E).*abs(E));
+
+%%Display result
+colormap(gray);
+imagesc(I2);
+
diff --git a/Table_7_7_Multiple_beam_interference.m b/Table_7_7_Multiple_beam_interference.m
new file mode 100644
index 0000000..8c2312d
--- /dev/null
+++ b/Table_7_7_Multiple_beam_interference.m
@@ -0,0 +1,26 @@
+%%Multiple beam interference
+%Define parameters
+N=500;%%Define size of the matrix
+Angle=1;%Define angle of the second plane wave
+V=0.25;%%Visibility controller
+lambda=0.632*1e-6;%Define wavelength
+
+%Sampling
+ del=1*1e-6;%sampling
+ x=-N/2:N/2-1;
+ y=-N/2:N/2-1;
+ [X,Y]=meshgrid(x*del,y*del);
+
+ %Simulate interference
+ A=V*exp(1i*(Y/lambda)*tand(Angle)*2*pi);
+ B=V*exp(1i*(X/lambda)*tand(-Angle)*2*pi);
+ C=V*exp(1i*(X/lambda)*tand(Angle)*2*pi);
+ D=V*exp(1i*(Y/lambda)*tand(-Angle)*2*pi);
+ I=A+B+C+D;%Interference of the object and reference wave
+
+%%Intensity profile
+ I1=abs(I).*abs(I);
+
+%Observation of the HOE
+colormap(gray)
+imagesc(I1)
diff --git a/Table_7_8_Fourier_hologram.m b/Table_7_8_Fourier_hologram.m
new file mode 100644
index 0000000..fe8c25f
--- /dev/null
+++ b/Table_7_8_Fourier_hologram.m
@@ -0,0 +1,26 @@
+%Fourier Hologram
+% Load the object and create the object and reference matrices
+N=500;%Define the size of the matrix
+A=imread('Image location');%Loading the image
+A=double(A);
+A=A(1:N,1:N);
+A=imresize(A,[N,N]);
+B=zeros(N,N);%Dirac delta function
+B(250,250)=100;
+
+%Calculate the far field diffraction patterns of the object and reference
+ A1=fftshift(fft2(A));%Calculate the diffraction pattern for the object
+ I1=abs(A1).*abs(A1);
+ B1=fftshift(fft2(B));%Calculate the diffraction pattern for the reference
+ I2=abs(B1).*abs(B1);
+
+%Create the hologram by interfering the far field diffracted fields of object and %reference
+ D1=A1+B1;
+ I3=abs(D1).*abs(D1);
+
+%Filtering
+I4=(I3-I1);%Filtering of autocorrelation term
+
+%Hologram reconstruction
+D2=fftshift(fft2(I4));
+I5=abs(D2).*abs(D2);
diff --git a/Table_7_9_Fresnel_hologram.m b/Table_7_9_Fresnel_hologram.m
new file mode 100644
index 0000000..0765493
--- /dev/null
+++ b/Table_7_9_Fresnel_hologram.m
@@ -0,0 +1,39 @@
+%Fresnel Hologram
+%Load object and create the object and reference matrices
+N=500;% Define the matrix size
+B=ones(N,N);%Reference wave
+A=imread('Image location');%Loading the image
+A=double(A);%Convert symbolic object to numeric object
+A=A(1:N,1:N);
+A=A/max(max(A));%Normalizing the matrix
+
+%Sampling
+ del=1;%
+ x=-N/2:N/2-1;
+ y=-N/2:N/2-1;
+ [X,Y]=meshgrid(x*del,y*del);
+
+%Calculate the Fresnel diffraction pattern of the object
+ A1=A.*exp(-1i*(pi/N)*(X.^2+Y.^2));
+ A2=fftshift(fft2(fftshift(A1)));
+ A3=A2.*exp(1i*(pi/N)*(X.^2+Y.^2));
+
+%Create interference between the object and reference waves
+ H1=B+A3;%Interference pattern
+
+%Filtering
+ I1=abs(A3).*abs(A3);
+ I2=abs(H1).*abs(H1);
+ H2=I2-I1;%Filtering
+ H2=H2/max(max(H2));%Normalizing the hologram
+
+%Reconstruct the hologram
+ H3=H2.*exp(-1i*(pi/N)*(X.^2+Y.^2));
+ H4=fftshift(fft2(fftshift(H3)));
+ H5=H4.*exp(1i*(pi/N)*(X.^2+Y.^2));
+ H6=(abs(H5).*abs(H5));
+ H7=rot90(H6,2);
+
+%Display reconstructed image
+colormap(gray)
+imagesc(H7)
diff --git a/Table_8_1_Line_data_lithography.m b/Table_8_1_Line_data_lithography.m
new file mode 100644
index 0000000..fc087e1
--- /dev/null
+++ b/Table_8_1_Line_data_lithography.m
@@ -0,0 +1,23 @@
+%Generating line data for lithography
+N=2000; % Size of the matrix
+A=zeros(N,N); % Define matrices
+B=zeros(N,N);
+P=100; % Define period of the grating
+p1=0; % Initialize line number
+for q=1:N;
+ if rem(q,P)<=P/2;
+ p1=p1+1;
+ X1(p1)=1;
+ Y1(p1)=q;
+ X2(p1)=N;
+ Y2(p1)=q;
+ end
+end
+fid = fopen('File location','wt');%Open a text file
+for s=1:p1;
+ fprintf(fid,'Draw Line %d,%d',X1(s),Y1(s));
+ fprintf(fid,'\n');
+ fprintf(fid,' %d,%d',X2(s),Y2(s));
+ fprintf(fid,'\n');
+end
+fclose(fid);%Close the text file
diff --git a/Table_8_1_Polygon_data_lithography.m b/Table_8_1_Polygon_data_lithography.m
new file mode 100644
index 0000000..b4b8644
--- /dev/null
+++ b/Table_8_1_Polygon_data_lithography.m
@@ -0,0 +1,24 @@
+%Generating polygon data for lithography
+N=2000; % Size of the DOE
+A=zeros(N,N); % Define matrices
+B=zeros(N,N);
+M=1000; % Define the number of points in the circle
+theta1=zeros(M);
+P=100; % Define period of the grating
+p1=0; % Initialize the iteration variable
+for r=1:300; % Generate point data of every circle
+ for theta=1:M;
+ if rem(r,P) x(n) && abs(q-N/2) < x(n+1);
+ A1(p,q)=1;
+ end
+ if abs(p-N/2) > y(n) && abs(p-N/2) < y(n+1);
+ A2(p,q)=1;
+ end
+ end
+ end
+end
+A3=exp(1i*pi*xor(A1,A2));
+
+%Observing the grating output in the far-field
+E=fftshift(fft2(A3)); %fftshift is used to re-order the terms in their natural order
+IN=(abs(E)/(N*N)).*(abs(E)/(N*N)); % Calculating intensity
+figure(1)
+colormap(gray);
+imagesc(angle(A3))
+figure(2)
+colormap(gray);
+imagesc(IN);
diff --git a/Table_E_2_5_Elliptical_FZP.m b/Table_E_2_5_Elliptical_FZP.m
new file mode 100644
index 0000000..db8a7b7
--- /dev/null
+++ b/Table_E_2_5_Elliptical_FZP.m
@@ -0,0 +1,37 @@
+%%Elliptical FZP%%
+clear; %Clear all memory
+
+%Define parameters
+N=500; %Define Matrix sizes
+M=25;%Define the number of grating lines
+A=ones(N,N);%Define Matrices by assigning zeros to all pixels
+rx=zeros(M,M);
+r=zeros(N,N);
+fx=3000; %Define focal lengths (in micrometers)
+fy=4500;
+lambda=0.633; %Define wavelength (in micrometers)
+
+%Construct the FZP
+for n=1:M; %Calculate the width of the grating lines
+ rx(n)=sqrt((n-1)*fx*lambda);
+end
+for n=1:2:M;%Construct elliptical FZP
+ for p=1:N;
+ for q=1:N;
+ r(p,q)=sqrt(((p-N/2)*(p-N/2))*(fx/fy)+((q-N/2)*(q-N/2)));
+ if r(p,q) > rx(n) && r(p,q) < rx(n+1);
+ A(p,q)=exp(1i*pi);
+ end
+ end
+ end
+end
+
+%Observing the grating output in the far-field
+E=fftshift(fft2(A)); %fftshift is used to re-order the terms in their natural order
+IN=(abs(E)/(N*N)).*(abs(E)/(N*N)); % Calculating intensity
+figure(1)
+colormap(gray);
+imagesc(angle(A))
+figure(2)
+colormap(gray);
+imagesc(IN);
\ No newline at end of file
diff --git a/Table_E_2_6_Axicon_array.m b/Table_E_2_6_Axicon_array.m
new file mode 100644
index 0000000..03d8116
--- /dev/null
+++ b/Table_E_2_6_Axicon_array.m
@@ -0,0 +1,33 @@
+%Axicon array
+clear; %Clear all memory
+
+%Define parameters
+N1=100; %Define matrix sizes
+N2=500;
+ratio=N2/N1;
+A1=zeros(N1,N1); %Define a Matrix by assigning 0 to all elements
+r=zeros(N1,N1);
+P=10; %Define the period of the axicon
+
+%Construct the axicon array
+for p=1:N1; %Generate the fundamental building block – single axicon
+ for q=1:N1;
+ r(p,q)=sqrt((p-N1/2)*(p-N1/2)+(q-N1/2)*(q-N1/2));
+ if r(p,q) (n1-1)*P/g && p <= (n1)*P/g ;
+ AF(p,q)=exp(1i*n1*delphase);
+ end
+ end
+ end
+end
+
+%Construct the full grating using repmat
+A=repmat(AF,5,1);
+
+%Observation of diffraction pattern
+AD=fftshift(fft2(A));
+IN=(abs(AD)/(N*N)).*(abs(AD)/(N*N));
+figure
+colormap(gray);
+imagesc(angle(A))
+figure
+colormap(gray);
+imagesc(IN);
\ No newline at end of file
diff --git a/Table_E_3_2_Four_level_Negative_FZP.m b/Table_E_3_2_Four_level_Negative_FZP.m
new file mode 100644
index 0000000..e84c0bd
--- /dev/null
+++ b/Table_E_3_2_Four_level_Negative_FZP.m
@@ -0,0 +1,32 @@
+%%4-level 1D FZP
+clear;%Clear all memory
+
+%Define the FZP parameters
+N=500;%Define Matrix size
+f=10000;%Define the focal length in ?m
+lambda=0.632;%Define the wavelength in ?m
+C=ones(N,N);
+g=4;
+
+%Constructing the negative FZP
+x=1:N;
+y=1:N;
+[X,Y]=meshgrid(x,y);
+r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+A=(f+sqrt(f*f+r.*r))*(2*pi)/(0.632);
+B=rem(A,2*pi);
+B1=exp(1i*B);
+B(r>N/2)=0;
+for p=1:N;%Construct the 4-level FZP
+ for q=1:N;
+ for n=1:g;
+ if B(p,q)> (n-1)*pi/2 && B(p,q) <= n*pi/2;
+ C(p,q)=exp(1i*n*pi/2);
+ end
+ end
+ end
+end
+
+%Observation of phase profile
+colormap(gray)
+imagesc(angle(C))
\ No newline at end of file
diff --git a/Table_E_3_3_Binary_SPP.m b/Table_E_3_3_Binary_SPP.m
new file mode 100644
index 0000000..086f5b3
--- /dev/null
+++ b/Table_E_3_3_Binary_SPP.m
@@ -0,0 +1,32 @@
+%Binary SPP with L=5
+clear;
+
+%Defining the SPP parameters
+N=500;%Define the matrix size
+L=5;%Define the topological charge number
+A3=zeros(N,N);
+
+%Constructing the SPP
+ x=1:N;
+ y=1:N;
+ [X,Y]=meshgrid(x,y);
+ theta=atan2((X-N/2),(Y-N/2));
+ r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+ A1=L*(theta+pi);
+ A2=rem(A1,2*pi);
+for p=1:N;%Construct the SPP using atan2 function and binarize
+ for q=1:N;
+ if rem(A2(p,q),2*pi)<=pi;
+ A3(p,q)=exp(1i*pi);
+ else
+ A3(p,q)=exp(1i*0);
+ end
+ end
+end
+A3(r>30)=0;
+
+%Observation of far field pattern
+E=fftshift(fft2(A3));%Calculate the fourier transform and rearrange terms
+I=abs(E).*abs(E);
+colormap(gray)
+imagesc(I)
diff --git a/Table_E_3_4_IFTA_Three_level.m b/Table_E_3_4_IFTA_Three_level.m
new file mode 100644
index 0000000..3f944cf
--- /dev/null
+++ b/Table_E_3_4_IFTA_Three_level.m
@@ -0,0 +1,55 @@
+%% Iterative Fourier Transform Algorithm%
+clear; % Clear all memory
+
+% Loading the target image
+N=500;% Matrix size
+target=imread('File location'); % Read image from file
+m=size(target,1); % Size of the image
+scale=N/m; % Estimate the necessary scaling factor
+target=imresize(target,scale); % Resize image to the matrix size
+target=double(target); % Convert symbolic object to a numeric object
+target=target/(max(max(target))); % Normalize matrix
+
+% Defining DOE phase
+DOE=2*pi*rand(N,N); % Generate a N x N matrix of random phase between 0 and 2?
+ns=5;
+
+% IFTA algorithm
+for t=1:ns; %Iterate to calculate the phase value
+ DOEphase=exp(1i*DOE);
+
+ % Forward iteration
+ iterf=fft2(DOEphase);
+ intf=abs(iterf);
+ angf=angle(iterf);
+ A=target.*exp(1i*angf);
+
+ % Backward iteration
+ iterb=ifft2(A);
+ angb=angle(iterb);
+ DOE=angb;
+ error=target-intf/max(max(intf)); %Calculate error
+ E=sum(sum(abs(error)))/(N*N);
+ if E<0.05;
+ iteration=t;
+ break
+ end
+end
+%%
+DOE1=zeros(N,N);
+g=3;%Define the number of phase levels
+delphase=2*pi/3;%Define the phase increment
+for p=1:N;%Convert the greyscale phase profile into a 3-level phase profile
+ for q=1:N;
+ for n=1:g;
+ if DOE(p,q)>-pi+(n-1)*delphase && DOE(p,q)<=-pi+(n)*delphase;
+ DOE1(p,q)=(n-1)*delphase;
+ end
+ end
+ end
+end
+%Verification of result
+DOE2=exp(1i*DOE1);
+I=abs(fft2(DOE2)); %Calculate the Fourier Transform
+colormap(gray)
+imagesc(I)
\ No newline at end of file
diff --git a/Table_E_3_5_Ringlens.m b/Table_E_3_5_Ringlens.m
new file mode 100644
index 0000000..8d475d1
--- /dev/null
+++ b/Table_E_3_5_Ringlens.m
@@ -0,0 +1,21 @@
+%%Greyscale ring lens
+clear;%Clear all memory
+
+%Define lens parameters
+N=500;%Define Matrix size
+f=2000;%Define the focal length in um
+r0=100;%Define radius of ring in um
+lambda=0.632;%Define the wavelength in um
+
+%Constructing the lens
+ x=1:N;
+ y=1:N;
+ [X,Y]=meshgrid(x,y);
+ r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+ A=(f-sqrt(f*f+(r-r0).*(r-r0)))*(2*pi)/(0.632);
+ B=rem(A,2*pi);
+ B1=exp(1i*B);
+
+%Observation of phase profile
+colormap(gray)
+imagesc(B)
diff --git a/Table_E_4_1_Blazed_axicon.m b/Table_E_4_1_Blazed_axicon.m
new file mode 100644
index 0000000..3cd4637
--- /dev/null
+++ b/Table_E_4_1_Blazed_axicon.m
@@ -0,0 +1,28 @@
+%Fresnel diffraction of blazed axicon%
+clear; %Clear all memory
+
+% Defining the parameters
+N=500;% Define the matrix size
+lambda=0.633*10^-6;%Define wavelength in meters
+z=50;%Propagation distance = 50 m
+P=10^-4;%Period of axicon = 0.1 mm
+wsamp=10*lambda;%sampling period or width
+
+%Sampling the space
+ x=1:N;
+ y=1:N;
+ [X,Y]=meshgrid(x,y);%Sampling
+ Rsamp=sqrt((X-N/2).^2+(Y-N/2).^2).*wsamp;%Define sampled radius
+
+%Constructing the DOE
+ A=exp(1i*(rem(Rsamp,P))*(2*pi)/P);
+ A(Rsamp>N/2*wsamp)=0;
+
+% Calculating the Fresnel diffraction
+ PPF=exp(1i*pi/(lambda*z).*Rsamp.*Rsamp); %Calculate the parabolic phase factor
+ A1=A.*PPF; %Multiply the circular aperture function with the parabolic phase factor
+ E=abs(fftshift(fft2(fftshift(A1)))); %Calculate Fourier Transform
+
+%Observation of diffraction pattern
+ colormap(gray)
+ imagesc(E)
\ No newline at end of file
diff --git a/Table_E_4_2_DOE.m b/Table_E_4_2_DOE.m
new file mode 100644
index 0000000..ae43a35
--- /dev/null
+++ b/Table_E_4_2_DOE.m
@@ -0,0 +1,27 @@
+%Fresnel diffraction of DOE
+clear; %Clear all memory
+
+% Defining the parameters
+N=500;% Define the matrix size
+lambda=0.633*10^-6;%Define wavelength in meters
+z=10*1e-3;%Propagation distance
+del=1*1e-6;%sampling period or width
+
+%Sampling the space
+ x=-N/2:N/2-1;
+ y=-N/2:N/2-1;
+ [X,Y]=meshgrid(x*del,y*del);%Sampling
+ R=sqrt(X.^2+Y.^2);%Define sampled radius
+
+%Constructing the DOE
+ A=exp(1i*(X+250*1e-6).^3*2*pi*5*10^11);
+
+%Calculating the Fresnel diffraction
+ PPF=exp(1i*pi/(lambda*z)*R.*R); %Calculate the parabolic phase factor
+ A1=A.*PPF; %Multiply the circular aperture function with the parabolic phase factor
+ E=abs(fftshift(fft2(fftshift(A1)))); %Calculate Fourier Transform
+ E=E/max(max(E));
+ figure (1)
+ imagesc(E);
+ figure (2)
+ plot(E(N/2,:)/(max(max(E(N/2,:)))))
diff --git a/Table_E_6_1_binary_FZP_Axicon.m b/Table_E_6_1_binary_FZP_Axicon.m
new file mode 100644
index 0000000..fd78368
--- /dev/null
+++ b/Table_E_6_1_binary_FZP_Axicon.m
@@ -0,0 +1,39 @@
+%Multifunctional DOE – Circular grating and FZP – Modulo-2? phase addition method
+clear; %Clear all memory
+
+%Define grating parameters
+N=1000;%Define matrix size
+M=10;%Define number of half period zones of FZP
+f=30000; %Define focal length of FZP
+lambda=0.633;%Define wavelength
+A1=zeros(N,N); %Define a Matrix by assigning 0 to all elements
+A2=zeros(N,N);
+r=zeros(N,N);
+r1=zeros(M);
+Pr=38; %Define the period of the grating
+FFr=0.5; %Define fill factors for x and y periodicity
+
+%Construct the binary FZP and binary circular grating
+for n=1:M;
+ r1(n)=sqrt(n*f*lambda);
+end
+for n=1:2:M;
+ for p=1:N;
+ for q=1:N;
+ r(p,q)=sqrt((p-N/2)*(p-N/2)+(q-N/2)*(q-N/2));
+ if r(p,q) r1(n) && r(p,q) < r1(n+1);
+ A2(p,q)=1;
+ end
+ end
+ end
+ end
+end
+A3=exp(1i*pi*xor(A1,A2)); %XOR operation between A1 and A2
+
+%Observation of the DOE
+colormap(gray)
+imagesc(angle(A3))
diff --git a/Table_E_6_2_Blazed_Axicon_binary_grating.m b/Table_E_6_2_Blazed_Axicon_binary_grating.m
new file mode 100644
index 0000000..62cf2c3
--- /dev/null
+++ b/Table_E_6_2_Blazed_Axicon_binary_grating.m
@@ -0,0 +1,22 @@
+%%Multifunctional DOE- Blazed axicon and binary grating
+clear;%Clear all memory
+%Define grating parameters
+N=500;%Define Matrix size
+Pr=51;
+Px=7;
+FFx=0.5;
+A2=zeros(N,N);%Define the matrices assigning ones to all pixels
+
+%Construction of blazed axicon and binary grating
+ x=1:N;
+ y=1:N;
+ [X,Y]=meshgrid(x,y);
+ r=sqrt((X-N/2).*(X-N/2)+(Y-N/2).*(Y-N/2));
+ P1=rem(r,Pr);
+ A1=(P1/Pr)*2*pi;
+ A2(rem(X,Px)