%% 

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);