diff --git a/Data-Analyzer/computeBondOrderParameters.m b/Data-Analyzer/computeBondOrderParameters.m
new file mode 100644
index 0000000..fe3ecee
--- /dev/null
+++ b/Data-Analyzer/computeBondOrderParameters.m
@@ -0,0 +1,179 @@
+function [psi_global, order_ratio, Npoints] = computeBondOrderParameters(I)
+% Analyze bond-orientational order in a 2D image, restricted to a circular region
+%
+% Input:
+%   I - 2D grayscale image (double or uint8)
+%
+% Outputs:
+%   psi_global  - struct with fields psi2, psi4, psi6 (⟨|ψₙ|⟩ values)
+%   order_ratio - scalar = ⟨|ψ₆|⟩ / ⟨|ψ₂|⟩
+%   Npoints     - number of detected points used
+
+    %% Step 1: Create mask
+    [Ny, Nx] = size(I);
+    % Parameters for elliptical mask
+    a = Nx / 2.25;         % Semi-major axis (x-direction)
+    b = Ny / 2.25;         % Semi-minor axis (y-direction)
+    theta_deg = 0;     % Rotation angle in degrees
+    theta = deg2rad(theta_deg);  % Convert to radians
+    
+    % Meshgrid coordinates
+    [Ny, Nx] = size(I);
+    [X, Y] = meshgrid(1:Nx, 1:Ny);
+    cx = Nx / 2;
+    cy = Ny / 2;
+    
+    % Shifted coordinates
+    Xc = X - cx;
+    Yc = Y - cy;
+    
+    % Rotate coordinates
+    Xr =  cos(theta) * Xc + sin(theta) * Yc;
+    Yr = -sin(theta) * Xc + cos(theta) * Yc;
+    
+    % Elliptical mask equation
+    mask = (Xr / a).^2 + (Yr / b).^2 <= 1;
+    
+    % Apply mask to image
+    I = double(I) .* mask;
+
+    %% Step 2: Detect local maxima within mask
+    I_smooth = imgaussfilt(I, 2);
+    peaks = imregionalmax(I_smooth);
+    [y, x] = find(peaks);
+    Npoints = length(x);
+
+    if Npoints < 6
+        warning('Too few points (%d) to compute bond order meaningfully.', Npoints);
+        psi_global = struct('psi2', NaN, 'psi4', NaN, 'psi6', NaN);
+        order_ratio = NaN;
+        return;
+    end
+
+    %% Step 3: Try Delaunay triangulation
+    use_kNN = false;
+    try
+        tri = delaunay(x, y);
+    catch
+        warning('Delaunay triangulation failed - falling back to kNN neighbor search');
+        use_kNN = true;
+    end
+
+    % 3) Find neighbors
+    neighbors = cell(Npoints, 1);
+    if ~use_kNN
+        % Use Delaunay neighbors
+        for i = 1:size(tri,1)
+            for j = 1:3
+                idx = tri(i,j);
+                nbrs = tri(i,[1 2 3] ~= j);
+                neighbors{idx} = unique([neighbors{idx}, nbrs]);
+            end
+        end
+    else
+        % Use kNN neighbors (e.g. k=6)
+        k = min(6, Npoints-1);
+        pts = [x, y];
+        [idx, ~] = knnsearch(pts, pts, 'K', k+1); % +1 for self
+        for j = 1:Npoints
+            neighbors{j} = idx(j, 2:end);
+        end
+    end
+
+    % 4) Compute bond order parameters
+    psi2 = zeros(Npoints,1);
+    psi4 = zeros(Npoints,1);
+    psi6 = zeros(Npoints,1);
+
+    for j = 1:Npoints
+        nbrs = neighbors{j};
+        nbrs(nbrs == j) = [];
+        if isempty(nbrs)
+            continue;
+        end
+        angles = atan2(y(nbrs) - y(j), x(nbrs) - x(j));
+        psi2(j) = abs(mean(exp(1i * 2 * angles)));
+        psi4(j) = abs(mean(exp(1i * 4 * angles)));
+        psi6(j) = abs(mean(exp(1i * 6 * angles)));
+    end
+
+    psi_global.psi2 = mean(psi2);
+    psi_global.psi4 = mean(psi4);
+    psi_global.psi6 = mean(psi6);
+
+    order_ratio = psi_global.psi6 / psi_global.psi2;
+
+    % --- After neighbors are found and (x,y) extracted ---
+
+    % 1) Angular distribution of all bonds pooled from all points
+    all_angles = [];
+    for j = 1:Npoints
+        nbrs = neighbors{j};
+        nbrs(nbrs == j) = [];
+        if isempty(nbrs), continue; end
+        angles = atan2(y(nbrs)-y(j), x(nbrs)-x(j));
+        all_angles = [all_angles; angles(:)];
+    end
+    
+    % Histogram of bond angles over [-pi, pi]
+    nbins = 36; % 10 degree bins
+    edges = linspace(-pi, pi, nbins+1);
+    counts = histcounts(all_angles, edges, 'Normalization','probability');
+    
+    % Circular variance of bond angles
+    R = abs(mean(exp(1i*all_angles)));
+    circ_variance = 1 - R;  % 0 means angles clustered; 1 means uniform
+    
+    % 2) Positional correlation function g(r) estimate
+    
+    % Compute pairwise distances
+    coords = [x y];
+    D = pdist(coords);
+    max_r = min([Nx, Ny])/2; % max radius for g(r)
+    
+    % Compute histogram of distances
+    dr = 1; % bin width in pixels (adjust)
+    r_edges = 0:dr:max_r;
+    [counts_r, r_bins] = histcounts(D, r_edges);
+    
+    % Normalize g(r) by ideal gas distribution (circle annulus area)
+    r_centers = r_edges(1:end-1) + dr/2;
+    area_annulus = pi*((r_edges(2:end)).^2 - (r_edges(1:end-1)).^2);
+    density = Npoints / (Nx*Ny);
+    
+    g_r = counts_r ./ (density * area_annulus * Npoints);
+
+    figure(1);
+    clf
+    set(gcf,'Position',[500 100 1250 500])
+    t = tiledlayout(1, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); % 1x4 grid
+
+    nexttile
+    imshow(I, []);
+    hold on;
+    scatter(x, y, 30, 'w', 'filled');
+    hold off;
+    colormap(Helper.Colormaps.plasma()); % Example colormap for original image plot
+    cbar1 = colorbar;
+    cbar1.Label.Interpreter = 'latex';
+    xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14);
+    ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14);
+    title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14);
+    axis square;
+
+    % --- Plot angular distribution ---
+    nexttile
+    polarhistogram(all_angles, nbins, 'Normalization', 'probability');
+    title(sprintf('Bond angle distribution (circ. variance = %.3f)', circ_variance));
+    
+    % --- Plot g(r) ---
+    nexttile
+    plot(r_centers, g_r, 'LineWidth', 2);
+    xlabel('r (pixels)');
+    ylabel('g(r)');
+    title('Radial Pair Correlation Function g(r)');
+    grid on;
+    axis square
+
+
+end
\ No newline at end of file
diff --git a/Data-Analyzer/execution_scripts.m b/Data-Analyzer/execution_scripts.m
new file mode 100644
index 0000000..653cfda
--- /dev/null
+++ b/Data-Analyzer/execution_scripts.m
@@ -0,0 +1,63 @@
+%% 
+
+Data       = load('D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500/Run_000/psi_gs.mat','psi','Params','Transf','Observ');
+
+Params     = Data.Params;
+Transf     = Data.Transf; 
+Observ     = Data.Observ; 
+if isgpuarray(Data.psi)
+    psi = gather(Data.psi);
+else
+    psi = Data.psi;
+end
+if isgpuarray(Data.Observ.residual)
+    Observ.residual = gather(Data.Observ.residual); 
+else
+    Observ.residual = Data.Observ.residual; 
+end
+% Axes scaling and coordinates in micrometers
+x    = Transf.x * Params.l0 * 1e6; 
+y    = Transf.y * Params.l0 * 1e6;
+z    = Transf.z * Params.l0 * 1e6;
+dz   = z(2)-z(1);
+% Compute probability density |psi|^2
+n    = abs(psi).^2;
+nxy  = squeeze(trapz(n*dz,3));
+IMG  = nxy;
+
+extractHexaticOrder(IMG)
+
+%% 
+
+Data       = load('D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500/Run_000/psi_gs.mat','psi','Params','Transf','Observ');
+% Data       = load('D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta40/HighN/aS_9.562000e+01_theta_040_phi_000_N_508333/Run_000/psi_gs.mat','psi','Params','Transf','Observ');
+
+Params     = Data.Params;
+Transf     = Data.Transf; 
+Observ     = Data.Observ; 
+if isgpuarray(Data.psi)
+    psi = gather(Data.psi);
+else
+    psi = Data.psi;
+end
+if isgpuarray(Data.Observ.residual)
+    Observ.residual = gather(Data.Observ.residual); 
+else
+    Observ.residual = Data.Observ.residual; 
+end
+% Axes scaling and coordinates in micrometers
+x    = Transf.x * Params.l0 * 1e6; 
+y    = Transf.y * Params.l0 * 1e6;
+z    = Transf.z * Params.l0 * 1e6;
+dz   = z(2)-z(1);
+% Compute probability density |psi|^2
+n    = abs(psi).^2;
+nxy  = squeeze(trapz(n*dz,3));
+IMG  = nxy;
+
+%
+
+[psi, ratio, N] = computeBondOrderParameters(IMG);
+
+fprintf('Points: %d\n⟨|ψ₂|⟩ = %.3f, ⟨|ψ₄|⟩ = %.3f, ⟨|ψ₆|⟩ = %.3f\n', N, psi.psi2, psi.psi4, psi.psi6);
+fprintf('(⟨|ψ₆|⟩ / ⟨|ψ₂|⟩) = %.3f\n', ratio);
diff --git a/Data-Analyzer/extractHexaticOrder.m b/Data-Analyzer/extractHexaticOrder.m
new file mode 100644
index 0000000..6c81b0c
--- /dev/null
+++ b/Data-Analyzer/extractHexaticOrder.m
@@ -0,0 +1,98 @@
+function extractHexaticOrder(I)
+    % Input:
+    %   I = 2D image matrix (grayscale, double)
+    % Output:
+    %   Plots original image + hexatic order + histogram
+    
+    %% 1) Detect dots (local maxima)
+    I_smooth = imgaussfilt(I, 2);
+    bw = imregionalmax(I_smooth);
+    [y, x] = find(bw);
+    
+    if length(x) < 6
+        error('Too few detected dots (%d) for meaningful hexatic order.', length(x));
+    end
+    
+    %% 2) Compute hexatic order parameter psi6
+    
+    % Delaunay triangulation to find neighbors
+    tri = delaunay(x,y);
+    N = length(x);
+    neighbors = cell(N,1);
+    for i = 1:size(tri,1)
+        neighbors{tri(i,1)} = unique([neighbors{tri(i,1)} tri(i,2) tri(i,3)]);
+        neighbors{tri(i,2)} = unique([neighbors{tri(i,2)} tri(i,1) tri(i,3)]);
+        neighbors{tri(i,3)} = unique([neighbors{tri(i,3)} tri(i,1) tri(i,2)]);
+    end
+    
+    psi6_vals = zeros(N,1);
+    for j = 1:N
+        nbrs = neighbors{j};
+        nbrs(nbrs == j) = [];
+        if isempty(nbrs)
+            psi6_vals(j) = 0;
+            continue;
+        end
+        angles = atan2(y(nbrs)-y(j), x(nbrs)-x(j));
+        psi6_vals(j) = abs(mean(exp(1i*6*angles)));
+    end
+    psi6_global = mean(psi6_vals);
+    
+    %% 3) Plotting
+    
+    figure('Name','Hexatic Order Analysis','NumberTitle','off', 'Units','normalized', 'Position',[0.2 0.2 0.6 0.7]);
+    t = tiledlayout(2,2, 'Padding','none', 'TileSpacing','compact');
+    
+    % Top left: Original image with detected dots (square)
+    ax1 = nexttile(1);
+    imshow(I, []);
+    hold on;
+    scatter(x, y, 30, 'w', 'filled');
+    hold off;
+    colormap(ax1, Helper.Colormaps.plasma()); % Example colormap for original image plot
+    cbar1 = colorbar(ax1);
+    cbar1.Label.Interpreter = 'latex';
+    xlabel(ax1, '$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14);
+    ylabel(ax1, '$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14);
+    title(ax1, '$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14);
+    axis(ax1, 'square');
+    
+    % Top right: Histogram of local |psi6| (square)
+    ax2 = nexttile(2);
+    histogram(psi6_vals, 20, 'FaceColor', 'b');
+    xlabel(ax2, '$|\psi_6|$', 'Interpreter', 'latex', 'FontSize', 14);
+    ylabel(ax2, 'Count');
+    title(ax2, 'Histogram of Local Hexatic Order');
+    xlim(ax2, [0 1]);
+    grid(ax2, 'on');
+    axis(ax2, 'square');
+    
+    % Bottom: Hexatic order magnitude on dots with bond orientation arrows (span 2 tiles)
+    ax3 = nexttile([1 2]);
+    scatter_handle = scatter(x, y, 50, psi6_vals, 'filled');
+    colormap(ax3, 'jet');  % Different colormap for hexatic order magnitude
+    cbar3 = colorbar(ax3);
+    cbar3.Label.String = '|$\psi_6$|';
+    cbar3.Label.Interpreter = 'latex';
+    caxis(ax3, [0 1]);
+    axis(ax3, 'equal');
+    axis(ax3, 'tight');
+    title(ax3, sprintf('Local Hexatic Order |\\psi_6|, Global = %.3f', psi6_global));
+    hold(ax3, 'on');
+    
+    for j = 1:N
+        nbrs = neighbors{j};
+        nbrs(nbrs == j) = [];
+        if isempty(nbrs)
+            continue;
+        end
+        angles = atan2(y(nbrs)-y(j), x(nbrs)-x(j));
+        psi6_complex = mean(exp(1i*6*angles));
+        local_angle = angle(psi6_complex)/6;
+        arrow_length = 5;
+        quiver(ax3, x(j), y(j), arrow_length*cos(local_angle), arrow_length*sin(local_angle), ...
+            'k', 'MaxHeadSize', 2, 'AutoScale', 'off');
+    end
+    hold(ax3, 'off');
+    
+end