Assistance in generating MATLAB code for creating a phase mask that converts a Gaussian beam to a flat top using a Spatial Light Modulator (SLM)

I aim to generate a flat top beam with dimensions of 18*3 microns from an incident Gaussian beam using a Spatial Light Modulator (SLM). To achieve this, I have implemented the Gerchberg-Saxton algorithm. However, I am facing difficulties in obtaining the desired flat top shape, and I am unsure about the technical errors that might be causing this issue.

I have followed the general steps of the Gerchberg-Saxton algorithm, including initialization, Fourier transform, and applying the phase mask. However, the resulting beam does not exhibit the desired flat top shape. I would greatly appreciate it if someone could provide a detailed explanation of potential sources of error in my implementation

%% Parameters

wavelength = 782*10^-9;    
focal = 8*10^-3; %lens 
pixel_SLM = 8*10^-6; %pitch

rows = 2048*2;    % 4096 or 1080
columns = 2048*2;    % 4096 or 1080

in = 3*10^-6;    % Y-axis of flat-top
di = 18*10^-6;    % X-axis of flat-top

%% Target profile - flat top 

d = di / (wavelength * focal / (rows * pixel_SLM));
a = in / (wavelength * focal / (rows * pixel_SLM));

centerx = rows / 2;
centery = columns / 2;

Image_Target = zeros(rows, columns);

% Flat-top generation
flatTopWidth = round(2 *d);    % Width of the flat-top region

% Calculate the indices for the flat-top region
startRow = round(centery - flatTopWidth / 2);
endRow = round(centery + flatTopWidth / 2);
startCol = round(centerx - flatTopWidth / 2);
endCol = round(centerx + flatTopWidth / 2);

% Set the flat-top region to 1
Image_Target(startRow:endRow, startCol:endCol) = 1;

Amp_Target = sqrt(Image_Target);    

figure
imagesc(Image_Target)
colorMap = jet(256);
colormap(colorMap);

figure
imagesc(Amp_Target)
colorMap = jet(256);
colormap(colorMap);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Guassian beam

Image_SLM_0 = zeros(rows, columns);

SLM_windowSize=1080;
SLM_windowSize2=1920;

SLM_windowSize3=1080;

diameter_gaussian=(3.423*10^-3)/10;                                   % 1/e^2

diameter=diameter_gaussian/(8*10^-6);
sigma=diameter/4;

Gau = fspecial('gaussian', SLM_windowSize3, sigma);


normg = Gau - min(Gau(:));
normg = normg ./ max(normg(:));

Gx(1)=centerx;
Gy(1)=centery;

randomRow = round(Gy(1)-SLM_windowSize3/2); %positionX
randomCol = round(Gx(1)-SLM_windowSize3/2); %positionY

Image_SLM_0(randomRow:randomRow+SLM_windowSize3-1, randomCol:randomCol+SLM_windowSize3-1) =  ...
    Image_SLM_0(randomRow:randomRow+SLM_windowSize3-1, randomCol:randomCol+SLM_windowSize3-1) + normg;


Amp_SLM=sqrt(Image_SLM_0); 

figure
imagesc(Image_SLM_0)
colorMap = jet(256);
colormap(colorMap);

figure
imagesc(Amp_SLM)

colorMap = jet(256);

%% Initial guess mask 

% size of the phase mask
a=1;
b=1;
rows = 4096;
cols = 4096;

x = linspace(-diameter, diameter, rows);
y = linspace(-diameter, diameter, cols);
[X, Y] = meshgrid(x, y);
hyperbolic_mask= sqrt((X.^2 / a^2) - (Y.^2 / b^2));

% hyperbolic phase mask
figure
imagesc(real(hyperbolic_mask))
colormap(hsv)
colorbar

Image_SLM=Amp_SLM.*hyperbolic_mask;

%% Gerchberg Saxton Algorithm

fin=20;

for j=1:fin
    
    clc
    disp(j)
    disp(fin)
  
    Image_Lens=fftshift(fft2(ifftshift(Image_SLM)));
    Image_Lens_final=Image_Lens;
    Phase_Lens=angle(Image_Lens);
    Image_Lens=Amp_Target .* exp(1i .* Phase_Lens);
    Image_SLM=fftshift(ifft2(ifftshift(Image_Lens)));
    Phase_SLM=angle(Image_SLM);
    Image_SLM=Amp_SLM .* exp(1i .* Phase_SLM);
    
end

     1

    20

     2

    20

     3

    20

     4

    20

     5

    20

     6

    20

     7

    20

     8

    20

     9

    20

    10

    20

    11

    20

    12

    20

    13

    20

    14

    20

    15

    20

    16

    20

    17

    20

    18

    20

    19

    20

    20

    20

Submatrix_SLM = Phase_SLM(centerx-SLM_windowSize/2:centerx+SLM_windowSize/2-1, centerx-SLM_windowSize2/2:centerx+SLM_windowSize2/2-1);

figure
imagesc(Phase_SLM)
colorMap = gray(256);
colormap(colorMap);

ExportSLM=Submatrix_SLM+pi; 


figure
imagesc(ExportSLM)
colorMap = gray(256);
colormap(colorMap);