%% =============================================================================================================
%% READ ME READ ME READ ME READ ME READ ME READ ME READ ME READ ME READ ME READ ME READ ME
%% =============================================================================================================
%  To use the following script, a MATLAB install is required. An account can be made at https://www.mathworks.com. 
%  Many universities pay for a license, in which case an account can be made with your university email. If this 
%  is not the case, you can use the online version which is free.
%  
%  1) First open a .cif (or other file format) using Olex2 and remove disorder components such that 
%  each Pd atom only has four coordinating N atoms. This can usually be done with the PARTs toggle or 
%  by manual deletion. If your file name is “Test1”, use Olex2 to generate a .xyz file by typing "file Test1.xyz" 
%  into the console, then drag your file into the MATLAB folder. 
%
%  2) Modify the settings below to use the name of your .xyz file and copy paste the entire
%  script (e.g., Ctrl + A, Ctrl + C, and Ctrl + V) into the MATLAB command line and press Enter to run the script.
%
%  3) The script will print the angle deconvolution of your cage into the command line, then display 
%  the data on three figures. You should use these figures to check the script has detected the coordination 
%  connectivity correctly. 
%
%  4) If the logic fails to construct connectivity correctly, manually input the connectivity below. 
%  Figure 2 displays the Pd and N labels for this purpose. Note: The logic typically fails when used on 
%  helicate lanterns (Pd2L4) with a large twist. So these require manual input.
%
%  To access the tabulated angles of individual ligands, navigate to the “Workspace” tab in the bottom left. 
%  Scroll down until you see “PdN_N_PdN_Table” and double-click to load. This table can be copy pasted into Excel.  

%% ==============================================================
%  Settings/Options
%% ==============================================================
filename = 'Pd9PART1.xyz';  % <-- Replace with your .xyz filename here

% Manual override for coordination connectivity
custom_connectivity = [];   % <-- If the script can't find your ligand then input it here like this: 
			    % ["Pd1N2N3Pd4", "Pd5N6N7Pd8", ...]. Otherwise must be left blank as []

graphical_output = 'on';    % <-- Toggle MATLAB output figures on or off, default on
fog = 'off';		    % <-- Toggle fog for depth perception on or off in the figures 
			    % Default is off, as this increases the time to run the script

%% SCRIPT STARTS BELOW:
%% ==============================================================
%  Step 1: Read XYZ file
%  - Loads atomic symbols and coordinates
%% ==============================================================
fprintf('Beginning calculation:\n');
fid = fopen(filename, 'r');
if fid == -1
    error('Failed to open file.');
end

numAtoms = str2double(fgetl(fid)); % Read number of atoms directly from the file
fgetl(fid); % skip the comment line

elements_raw = strings(numAtoms, 1);
coords_raw   = zeros(numAtoms, 3);

for i = 1:numAtoms
    line = strsplit(strtrim(fgetl(fid)));
    elements_raw(i) = string(line{1});  
    coords_raw(i, :) = str2double(line(2:4));
end
fclose(fid);

%% ==============================================================
%  Step 2: Identify Pd and N atoms and label them
%% ==============================================================
pd_idx   = find(elements_raw == "Pd");
n_idx    = find(elements_raw == "N");
pd_coords = coords_raw(pd_idx, :);
n_coords  = coords_raw(n_idx, :);

numPd = numel(pd_idx);
numN  = numel(n_idx);

pd_labels = "Pd" + (1:numPd)'; % Number Pd atoms based on order of appearance 
n_labels  = "N"  + (1:numN)';  % Number N atoms based on order of appearance 

fprintf('Number of Pd atoms: %d\n', numPd);
fprintf('Number of N atoms:  %d\n', numN);


%% ==============================================================
% Start override if loop
if isempty(custom_connectivity) % override for custom input
%% ==============================================================

%% ==============================================================
%  Step 3: Pd–N Coordination Table
%  - Assign each Pd its 4 closest unassigned N atoms
%  - Store a table of strings like "Pd1N3N5N12N22"
%% ==============================================================
if numN < 4 * numPd
    warning('Not enough N atoms to assign 4 per Pd.');
end

assigned_N = false(numN, 1);
pd_n_table = strings(numPd, 1);

for iPd = 1:numPd
    pd_pos   = pd_coords(iPd, :);
    distsN   = vecnorm(n_coords - pd_pos, 2, 2);
    distsN(assigned_N) = inf;                    % exclude already assigned N
    [~, nearest_idx]   = mink(distsN, 4);        % get 4 nearest unique N
    
    if any(isinf(distsN(nearest_idx)))
        error('Not enough unique N atoms left for %s.', pd_labels(iPd));
    end
    
    assigned_N(nearest_idx) = true;
    n_ids = n_labels(nearest_idx);
    
    % e.g. "Pd1" + "N2" + "N3" + "N4" + "N5"
    pd_n_table(iPd) = pd_labels(iPd) + strjoin(n_ids, '');
end

pd_n_coord_table = table(pd_n_table, 'VariableNames', {'PdN_Assignments'});

%disp('--- Pd–N Coordination Table ---');
%disp(pd_n_coord_table);

%% ==============================================================
%  Step 4: Evaluate Pd–Pd Proximity
%  - Setup the flag regular_cage = false if PdnL2n, n < 6 to indicate it does not have tetravalent framework topology
%  - Find 4 nearest Pd neighbours for each Pd atom
%  - Build table of unique Pd-Pd pairs like "Pd1Pd2"
%  - Also set flag false if any Pd-Pd distance is outside of 1.25 - 0.75 of the mean (indicating a irregular geometry)
%% ==============================================================

regular_cage = true;  

if numPd > 5
    pd_pairs = strings(0, 1); % Store Pd-Pd pairs

    for i = 1:numPd
        distances = vecnorm(pd_coords - pd_coords(i, :), 2, 2);
        distances(i) = inf; % Ignore self
        [~, nearest_idx] = mink(distances, 4);

        for j = 1:4
            pair = sort([pd_labels(i), pd_labels(nearest_idx(j))]);
            pd_pairs(end+1) = pair(1) + pair(2); % Standard form e.g., Pd1Pd2
        end
    end

    unique_pd_pairs = unique(pd_pairs);
    pd_pd_table = table(unique_pd_pairs, 'VariableNames', {'PdPairs'});
    
    % Calculate mean Pd-Pd distance then check bounds
    all_pd_distances = [];
    for i = 1:numPd
        distances = vecnorm(pd_coords - pd_coords(i, :), 2, 2);
        distances(i) = inf; % Ignore self
        nearest_distances = mink(distances, 4);
        all_pd_distances = [all_pd_distances; nearest_distances];
    end
    
    mean_pd_distance = mean(all_pd_distances);
    lower_bound = 0.75 * mean_pd_distance;
    upper_bound = 1.25 * mean_pd_distance;
    
    % Check if any Pd-Pd distance falls outside the acceptable range
    if any(all_pd_distances < lower_bound) || any(all_pd_distances > upper_bound)
        regular_cage = false;
    end
else
    regular_cage = false;
end

if regular_cage == false
    fprintf('Irregular geometry detected - PATH B\n');
else
    fprintf('Assuming tetravalent framework topoology - PATH A.\n')
end

%% ==============================================================
%  Path A - Assuming tetravalent framework topoology
%% ==============================================================
if regular_cage == true

%% ==============================================================
%  Step 5A: Construct Pd–N–N–Pd Coordination
%  For each Pd-Pd pair:
%    - Find N atom in Pd1's set closest to Pd2
%    - Find N atom in Pd2's set closest to Pd1  
%    - Compose: Pd1 + N1 + N2 + Pd2, e.g., "Pd1N63N7Pd2"
%% ==============================================================

% Extract Pd-Pd pairs
pairs = pd_pd_table.PdPairs; % e.g., "Pd1Pd2"

% Map from Pd labels to their N labels 
pdToNs = containers.Map();
for i = 1:numPd
    pdEntry = pd_n_coord_table.PdN_Assignments(i);
    pdEntryStr = char(pdEntry);

    % Extract Pd label: 'Pd' followed by digits
    pdLabelTokens = regexp(pdEntryStr, 'Pd\d+', 'match');
    if isempty(pdLabelTokens)
        error('Could not extract Pd label from PdN_Assignment: %s', pdEntryStr);
    end
    pdLabel = strtrim(pdLabelTokens{1});

    % Extract all N atoms (e.g., 'N63', 'N33', etc.)
    nAtoms = regexp(pdEntryStr, 'N\d+', 'match');

    pdToNs(pdLabel) = nAtoms;
end

% Debug
% disp('Keys in pdToNs map:');
% disp(pdToNs.keys);

% Allocate one for each Pd-Pd pair
numPairs = length(pairs);
PdN_N_PdN = strings(numPairs, 1);

for k = 1:numPairs
    pairStr = pairs(k); % e.g., 'Pd1Pd2'

    % Parse Pd labels for extraction
    tokens = regexp(pairStr, 'Pd(\d+)Pd(\d+)', 'tokens');

    tokens = tokens{1};
    pd1 = "Pd" + tokens{1};
    pd2 = "Pd" + tokens{2};

    % Retrieve their N atom labels 
    nPd1 = pdToNs(pd1);
    nPd2 = pdToNs(pd2);

    % Convert N labels to indices to find their coords
    nPd1_idx = zeros(length(nPd1), 1);
    for iN = 1:length(nPd1)
        nPd1_idx(iN) = find(n_labels == nPd1{iN});
    end
    nPd2_idx = zeros(length(nPd2), 1);
    for iN = 1:length(nPd2)
        nPd2_idx(iN) = find(n_labels == nPd2{iN});
    end

    % Coordinates of N atoms for Pd1 and Pd2
    nPd1_coords = n_coords(nPd1_idx, :);
    nPd2_coords = n_coords(nPd2_idx, :);

    % Retrieve their Pd coordinates
    idx1 = sscanf(pd1, 'Pd%d');
    idx2 = sscanf(pd2, 'Pd%d');
    pd1_pos = pd_coords(idx1, :);
    pd2_pos = pd_coords(idx2, :);

    % Find N atom in Pd1's set closest to Pd2
    dists1 = vecnorm(nPd1_coords - pd2_pos, 2, 2);
    [~, minIdx1] = min(dists1);
    closestNtoPd2 = nPd1{minIdx1};

    % Find N atom in Pd2's set closest to Pd1
    dists2 = vecnorm(nPd2_coords - pd1_pos, 2, 2);
    [~, minIdx2] = min(dists2);
    closestNtoPd1 = nPd2{minIdx2};

    % Compose output string: e.g., "Pd1N63N7Pd2"
    PdN_N_PdN(k) = pd1 + closestNtoPd2 + closestNtoPd1 + pd2;
end
% Create the coordination table
PdN_N_PdN_Table = table(PdN_N_PdN, 'VariableNames', {'Pd_N_N_Pd_Coordination'});

% Debug
%disp('--- Pd-N-N-Pd Coordination Table ---');
%disp(PdN_N_PdN_Table);

%% ==============================================================
%  End Path A
%% ==============================================================
else
end

%% ==============================================================
%  Path B - Fallback logic for Non-tetravalent framework topology OR irregular geometry
%% ==============================================================
if regular_cage == false

%% ==============================================================
%  Step 5B Part 1: Build N–N Coordination 
%  For each Pd:
%    - For each of its 4 N atoms, find the nearest N among all other Pd’s N atoms.
%    - Use this logic to find the closest N-N pairs.
%  Store a table of unique pairs: "N3N4", "N5N24", ...
%% ==============================================================
% Array for N–N pairs
N_pairs = strings(0, 1);

for iPd = 1:numPd
    % N atoms attached to this Pd
    rowStr   = char(pd_n_coord_table.PdN_Assignments(iPd));
    nAtoms_i = regexp(rowStr, 'N\d+', 'match');   % Cell array like {'N3','N5'}
    
    % Indices & coords for this Pd’s N atoms
    nIdx_i = zeros(length(nAtoms_i), 1);
    for k = 1:length(nAtoms_i)
        nIdx_i(k) = find(n_labels == nAtoms_i{k});
    end
    nCoords_i = n_coords(nIdx_i, :);
    
    % Collect indices of N atoms attached to other Pd atoms
    otherIdx = [];
    for jPd = 1:numPd
        if jPd == iPd, continue; end
        rowStr_j   = char(pd_n_coord_table.PdN_Assignments(jPd));
        nAtoms_j   = regexp(rowStr_j, 'N\d+', 'match');
        for k = 1:length(nAtoms_j)
            otherIdx(end+1) = find(n_labels == nAtoms_j{k}); %#ok<AGROW>
        end
    end
    
    otherNCoords = n_coords(otherIdx, :);
    otherNAtoms  = n_labels(otherIdx); 
    
    % For each N of current Pd, find closest N from other Pd’s sets
    for a = 1:size(nCoords_i, 1)
        dists = vecnorm(otherNCoords - nCoords_i(a, :), 2, 2);
        [~, minIdx] = min(dists);
        
        N1 = string(nAtoms_i{a});
        N2 = otherNAtoms(minIdx);
        
        % Standardise pair order to enable dedupe later
        pair = sort([N1, N2]);
        N_pairs(end+1) = pair(1) + pair(2); %#ok<AGROW>
    end
end

% Dedupe N–N pairs
unique_N_pairs = unique(N_pairs);

N_N_Table = table(unique_N_pairs, 'VariableNames', {'NN_Pairs'});
%disp('--- N–N Coordination Table ---');
%disp(N_N_Table);

%% ==============================================================
%  Step 5B Part 2: Construct Pd–N–N–Pd from N–N + Pd–N map 
%  For each unique N-N pair "N1N2":
%    - Find Pd attached to N1 and Pd attached to N2 (from Pd–N table)
%    - Compose: Pd(N1) + N1 + N2 + Pd(N2), e.g., "Pd1N3N7Pd2"
%  Store Pd–N–N–Pd strings in coordination table
%% ==============================================================
% Build a mapping N -> Pd (each N is uniquely assigned to one Pd in Step 3)
N_to_Pd = containers.Map('KeyType','char','ValueType','char');
for iPd = 1:numPd
    rowStr   = char(pd_n_coord_table.PdN_Assignments(iPd));
    nAtoms   = regexp(rowStr, 'N\d+', 'match');
    for k = 1:length(nAtoms)
        N_to_Pd(nAtoms{k}) = char(pd_labels(iPd));
    end
end

% Compose output string: e.g., "Pd1N63Pd2N7"
PdN_N_PdN = strings(length(unique_N_pairs), 1);
for k = 1:length(unique_N_pairs)
    pairStr = char(unique_N_pairs(k));       % e.g., 'N3N4'
    Ns = regexp(pairStr, 'N\d+', 'match');   % {'N3','N4'}

    N1 = string(Ns{1});
    N2 = string(Ns{2});
    
    Pd1 = string(N_to_Pd(char(N1)));
    Pd2 = string(N_to_Pd(char(N2)));
    
    PdN_N_PdN(k) = Pd1 + N1 + N2 + Pd2;
end

% Create table
PdN_N_PdN_Table = table(PdN_N_PdN, 'VariableNames', {'Pd_N_N_Pd_Coordination'});

%disp('--- Pd-N-N-Pd Coordination Table ---');
%disp(PdN_N_PdN_Table);

%% ==============================================================
%  End Path B
%% ==============================================================
else
end % closes the if statement for Path B

%% ==============================================================
% End override if loop
%% ==============================================================
else
    fprintf('Using custom connectivity\n');
    
    if iscell(custom_connectivity)
        custom_connectivity = string(custom_connectivity);
    end
    
    % Ensure column vector
    custom_connectivity = custom_connectivity(:);
    
    % Create the coordination table
    PdN_N_PdN_Table = table(custom_connectivity, 'VariableNames', {'Pd_N_N_Pd_Coordination'});
    
    fprintf('Custom connectivity table created with %d entries:\n', height(PdN_N_PdN_Table));
    disp(PdN_N_PdN_Table);
end % closes the custom coordination if statement

%% ==============================================================
%  Compute Pd-N-N-Pd Dot Product (THETA)
%  This represents the angle between Pd-N coordination vectors;
%  Calculated using Equation 1.
%  
%  For each entry in "PdN_N_PdN_Table":
%	- Parse the coordination string, eg. "Pd1N3N7Pd2"
%	- Extract individual atom labels (Pd1, N3, N7, Pd2)
%	- Find their coordinates
%	- Find the two coordination vectors
%	- Compute the dot product between these two vectors and append it to the table
%% ==============================================================

% Count PdN_N_PdN_Table to determine the entries required for each angle
numEntries = height(PdN_N_PdN_Table);

angles_deg = zeros(numEntries, 1);

% Loop over all PdN_N_PdN_Table entries
for k = 1:numEntries
    entry = PdN_N_PdN_Table.Pd_N_N_Pd_Coordination(k);
    entryStr = char(entry); 

    % Extract labels
    % Pattern: Pd<number>N<number>N<number>Pd<number>
    tokens = regexp(entryStr, 'Pd(\d+)N(\d+)N(\d+)Pd(\d+)', 'tokens');
    tokens = tokens{1}; % cell array of strings: {'A','B','C','D'}

    % Construct full labels 
    PdA_label = "Pd" + tokens{1};
    NB_label = "N" + tokens{2};
    NC_label = "N" + tokens{3};
    PdD_label = "Pd" + tokens{4};

    % Find indices of these atoms in the global coordinate lists
    idx_PdA = find(pd_labels == PdA_label);
    idx_PdD = find(pd_labels == PdD_label);
    idx_NB = find(n_labels == NB_label);
    idx_NC = find(n_labels == NC_label);

    % Coordinates
    PdA_pos = pd_coords(idx_PdA, :);
    PdD_pos = pd_coords(idx_PdD, :);
    NB_pos = n_coords(idx_NB, :);
    NC_pos = n_coords(idx_NC, :);

    % Construct coordination vectors:
    v1 = NB_pos - PdA_pos; 
    v2 = NC_pos - PdD_pos; 

    cos_theta = dot(v1, v2) / (norm(v1) * norm(v2));

    % Clamp to [-1,1]
    cos_theta = max(min(cos_theta, 1), -1);

    theta_rad = acos(cos_theta);
    theta_deg = rad2deg(theta_rad);

    % Store result
    angles_deg(k) = theta_deg;
end

% Append to the table
PdN_N_PdN_Table.Angle_deg = angles_deg;

% Debug
%disp(PdN_N_PdN_Table);

%% ==============================================================
%  Compute Pd-N-N-Pd Dihedral angles (PHI)
%  The angle between PdA-NB-NC and NB-NC-PdD planes is the twist angle between Pd-N coordination vectors
%  Calculated using Equation 2.
%
%  For each entry in "PdN_N_PdN_Table":
%	- Use same coordinate extraction logic as before
%	- Calculate plane normals for PdA-NB-NC and NB-NC-PdD
%	- Calculate the angle between the two planes and append it to the table
%% ==============================================================

dihedral_deg = zeros(numEntries, 1);

% Loop over PdN_N_PdN_Table
for k = 1:numEntries
    entry = PdN_N_PdN_Table.Pd_N_N_Pd_Coordination(k);
    entryStr = char(entry);

    % Extract atom labels
    tokens = regexp(entryStr, 'Pd(\d+)N(\d+)N(\d+)Pd(\d+)', 'tokens');
    tokens = tokens{1};

    PdA_label = "Pd" + tokens{1};
    NB_label  = "N"  + tokens{2};
    NC_label  = "N"  + tokens{3};
    PdD_label = "Pd" + tokens{4};

    % Index lookups
    idx_PdA = find(pd_labels == PdA_label);
    idx_PdD = find(pd_labels == PdD_label);
    idx_NB  = find(n_labels == NB_label);
    idx_NC  = find(n_labels == NC_label);

    % Coordinates
    PdA_pos = pd_coords(idx_PdA, :);
    PdD_pos = pd_coords(idx_PdD, :);
    NB_pos  = n_coords(idx_NB, :);
    NC_pos  = n_coords(idx_NC, :);

    % Calculate vectors
    v1_plane1 = PdA_pos - NB_pos; 
    v2_plane1 = NC_pos - NB_pos;   
    v1_plane2 = NB_pos - NC_pos;   
    v2_plane2 = PdD_pos - NC_pos;  

    % Calculate plane normals
    n1 = cross(v1_plane1, v2_plane1);  % Normal to plane PdA-NB-NC
    n2 = cross(v1_plane2, v2_plane2);  % Normal to plane NB-NC-PdD
    
    % Normalise
    n1 = n1 / norm(n1);
    n2 = n2 / norm(n2);

    % Calculate angle between planes
    cos_angle = dot(n1, n2);
    
    % Clamp to [-1,1]
    cos_angle = max(min(cos_angle, 1), -1);
    
    angle_rad = acos(cos_angle);
    angle_deg = rad2deg(angle_rad);

    % Store the result
    dihedral_deg(k) = angle_deg;
    
end

% Append to the table
PdN_N_PdN_Table.Dihedral_deg = dihedral_deg;

% Debug
%disp(PdN_N_PdN_Table);

%% ============================================================== 
% Compute in-plane angle (GAMMA) from THETA (dot-product) and PHI (dihedral) and append it to the table
% This angle represents the in-plane divergence and is calculated using Equation 3 from the manuscript.
%% ==============================================================

gamma_deg = zeros(numEntries, 1);

for k = 1:numEntries
    theta_rad = deg2rad(PdN_N_PdN_Table.Angle_deg(k));     
    phi_rad   = deg2rad(PdN_N_PdN_Table.Dihedral_deg(k));   

    % Compute gamma using Equation 3
    numerator   = cos(theta_rad) + 1;
    denominator = cos(phi_rad / 2)^2;
    arg         = numerator / denominator - 1;

    % Clamp to [-1,1]
    arg = max(min(arg, 1), -1);

    gamma_rad = acos(arg);
    gamma_deg(k) = rad2deg(gamma_rad);  
end

% Append to table
PdN_N_PdN_Table.InPlaneAngle_deg = gamma_deg;


%% ============================================================== 
%  Calculate and print statistics: 
%% ============================================================== 

% Display table with nice headings then restore the original
original_table = PdN_N_PdN_Table;
PdN_N_PdN_Table.Properties.VariableNames = {'Pd-N-N-Pd Coordination', 'Coordination Angle (θ) /°', 'Dihedral Angle (𝜑) /°', 'In-plane Angle (γ) /°'};
disp(PdN_N_PdN_Table);
PdN_N_PdN_Table = original_table;

% Calculate mean and population standard deviation values
mean_angle = mean(angles_deg);    
std_angle = std(angles_deg, 1);  

mean_dihedral = mean(dihedral_deg);
std_dihedral = std(dihedral_deg, 1);

mean_gamma = mean(gamma_deg);
std_gamma  = std(gamma_deg, 1);

fprintf('The following values are given as mean ± population standard deviation:\n')
fprintf('  θ = %.4f° ± %.4f°\n', mean_angle, std_angle);
fprintf('  𝜑 = %.4f° ± %.4f°\n', mean_dihedral, std_dihedral);
fprintf('  γ = %.4f° ± %.4f°\n', mean_gamma, std_gamma);

%% ============================================================== 
%  Debug Check: Ensure no N atom appears more than once in Pd-N-N-Pd Table
%% ============================================================== 
used_N_atoms = strings(0,1); 

for k = 1:height(PdN_N_PdN_Table)
    entry = PdN_N_PdN_Table.Pd_N_N_Pd_Coordination(k);
    tokens = regexp(entry, 'N\d+', 'match'); % Extract all N labels
    used_N_atoms = [used_N_atoms; string(tokens(:))]; % Append
end

% Check for duplicates
[unique_N, ~, ic] = unique(used_N_atoms);
counts = accumarray(ic, 1);
if any(counts > 1)
    fprintf('Valency violation. There is likely a mistake in your dataset\n');
end

%% ==============================================================
%% PLOTTING SECTION BELOW:
%% ==============================================================
if strcmpi(graphical_output, 'off')
    % No plotting, just exit
else
    % Time Figure plotting
    plottingTimer = tic;
    fprintf('Graphical output...\n'); % and the below section runs
%% ==============================================================
%% Common setup for all figures
%% ==============================================================

%% Extract filename for titles
[~, nameWithoutExt, ~] = fileparts(filename);

% Collect unique Pd and N atoms from Pd-N-N-Pd table
allPdNLabels = strings(0);
for k = 1:height(PdN_N_PdN_Table)
    entryStr = PdN_N_PdN_Table.Pd_N_N_Pd_Coordination(k);
    entryStr = char(entryStr);
    pdLabels = regexp(entryStr, 'Pd\d+', 'match');
    nLabels  = regexp(entryStr, 'N\d+', 'match');
    allPdNLabels = [allPdNLabels; pdLabels(:); nLabels(:)];
end
allPdNLabels = unique(allPdNLabels);

% Combine Pd and N coords and labels into one list
allCoords = [];
allAtomTypes = {};
for i = 1:length(allPdNLabels)
    label = allPdNLabels(i);
    if startsWith(label, 'Pd')
        idx = find(pd_labels == label);
        allCoords(end+1, :) = pd_coords(idx, :);
        allAtomTypes{end+1} = 'Pd';
    elseif startsWith(label, 'N')
        idx = find(n_labels == label);
        allCoords(end+1, :) = n_coords(idx, :);
        allAtomTypes{end+1} = 'N';
    end
end
allAtomTypes = allAtomTypes';

[sx, sy, sz] = sphere(7);  % Resolution of spheres
baseRadius = 0.3;
radius_pd = 2.5 * baseRadius;  % For Figures 1 & 2
radius_n = baseRadius;
radius_pd_full = 3 * baseRadius;  % For Figure 3 

colors = containers.Map({'Pd', 'N'}, {[0.000 0.412 0.522], [0.188, 0.314, 0.973]});

% View settings
margin = 4.5;
minCoords = min(allCoords);
maxCoords = max(allCoords);
ranges = maxCoords - minCoords;
maxRange = max(ranges);
center = (maxCoords + minCoords) / 2;
view_limits = [center(1) - maxRange/2 - margin, center(1) + maxRange/2 + margin;
                center(2) - maxRange/2 - margin, center(2) + maxRange/2 + margin;
                center(3) - maxRange/2 - margin, center(3) + maxRange/2 + margin];
view_azimuth = 45;
view_elevation = 30;
%% ==============================================================
%% FIGURE 1: Skeleton with angle annotations
%% ==============================================================

fig1 = figure('Name', '\1: Plot', 'Color', [1 1 1]);
ax1 = axes('Parent', fig1);
hold(ax1, 'on');

% Batch plotting
pd_indices = find(strcmp(allAtomTypes, 'Pd'));
n_indices = find(strcmp(allAtomTypes, 'N'));

% Plot Pd atoms
for i = 1:length(pd_indices)
    idx = pd_indices(i);
    c = allCoords(idx, :);
    surf(radius_pd*sx + c(1), radius_pd*sy + c(2), radius_pd*sz + c(3), ...
        'FaceColor', colors('Pd'), 'EdgeColor', 'none');
end

% Plot N atoms
for i = 1:length(n_indices)
    idx = n_indices(i);
    c = allCoords(idx, :);
    surf(radius_n*sx + c(1), radius_n*sy + c(2), radius_n*sz + c(3), ...
        'FaceColor', colors('N'), 'EdgeColor', 'none');
end

% Bond plotting
bond_data = cell(height(PdN_N_PdN_Table), 1);
for k = 1:height(PdN_N_PdN_Table)
    entryStr = char(PdN_N_PdN_Table.Pd_N_N_Pd_Coordination(k));
    tokens = regexp(entryStr, '(Pd\d+|N\d+)', 'match');
    
    if length(tokens) == 4
        bond_data{k} = tokens;
        % Draw bonds: Pd-N, N-N, N-Pd
        pairs = [1 2; 2 3; 3 4];
        for p = 1:size(pairs,1)
            idx1 = find(allPdNLabels == tokens{pairs(p,1)});
            idx2 = find(allPdNLabels == tokens{pairs(p,2)});
            c1 = allCoords(idx1, :); c2 = allCoords(idx2, :);
            
            type1 = allAtomTypes{idx1}; type2 = allAtomTypes{idx2};
            
            if (strcmp(type1, 'Pd') && strcmp(type2, 'N')) || (strcmp(type1, 'N') && strcmp(type2, 'Pd'))
                plot3([c1(1), c2(1)], [c1(2), c2(2)], [c1(3), c2(3)], 'k-', 'LineWidth', 1.5);
            elseif strcmp(type1, 'N') && strcmp(type2, 'N')
                plot3([c1(1), c2(1)], [c1(2), c2(2)], [c1(3), c2(3)], 'r-', 'LineWidth', 1);
            end
        end
    end
end

% Label Pd atoms
for i = 1:length(pd_indices)
    idx = pd_indices(i);
    coord = allCoords(idx, :);
    text(coord(1), coord(2), coord(3) + 1, allPdNLabels(idx), ...
        'FontSize', 8, 'FontWeight', 'bold', 'Color', 'black', ...
        'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom');
end

% Angle annotations
for k = 1:height(PdN_N_PdN_Table)
    entryStr = char(PdN_N_PdN_Table.Pd_N_N_Pd_Coordination(k));
    tokens = regexp(entryStr, '(N\d+)', 'match');
    if numel(tokens) == 2
        idx1 = find(allPdNLabels == tokens{1});
        idx2 = find(allPdNLabels == tokens{2});
        midpoint = (allCoords(idx1, :) + allCoords(idx2, :)) / 2;
        
        annotationText = sprintf('θ=%.1f°\nφ=%.1f°\nγ=%.1f°', ...
            PdN_N_PdN_Table.Angle_deg(k), ...
            PdN_N_PdN_Table.Dihedral_deg(k), ...
            PdN_N_PdN_Table.InPlaneAngle_deg(k));
        
        text(midpoint(1), midpoint(2), midpoint(3), annotationText, ...
            'FontSize', 7, 'FontWeight', 'normal', ...
            'HorizontalAlignment', 'center', ...
            'VerticalAlignment', 'middle', ...
            'Color', [0.1 0.1 0.1], ...
            'BackgroundColor', [1 1 1 0.5], ...
            'Margin', 1);
    end
end

% Apply lighting and material settings
xlim(view_limits(1,:)); ylim(view_limits(2,:)); zlim(view_limits(3,:));
daspect([1 1 1]); 
camlight('right'); camlight('left'); camlight('headlight');
lighting gouraud;
material([0.4 0.6 0.5 40 1]);  % [ambient diffuse specular shininess]
set(gca, 'XGrid','off','YGrid','off','ZGrid','off','Box','on',...
         'XTickLabel',[],'YTickLabel',[],'ZTickLabel',[]);
xlabel('X'); ylabel('Y'); zlabel('Z');
title([nameWithoutExt ' Plot'], 'FontSize', 11, 'FontWeight', 'bold');
view(view_azimuth, view_elevation);

%% ==============================================================
%% FIGURE 2: Skeleton with all labels
%% ==============================================================

fig2 = figure('Name', '\2: Skeletal', 'Color', [1 1 1]);
ax2 = axes('Parent', fig2);
hold(ax2, 'on');

% Reuse plotting from Figure 1
for i = 1:length(pd_indices)
    idx = pd_indices(i); c = allCoords(idx, :);
    surf(radius_pd*sx + c(1), radius_pd*sy + c(2), radius_pd*sz + c(3), ...
        'FaceColor', colors('Pd'), 'EdgeColor', 'none');
end
for i = 1:length(n_indices)
    idx = n_indices(i); c = allCoords(idx, :);
    surf(radius_n*sx + c(1), radius_n*sy + c(2), radius_n*sz + c(3), ...
        'FaceColor', colors('N'), 'EdgeColor', 'none');
end

for k = 1:height(PdN_N_PdN_Table)
    if ~isempty(bond_data{k})
        tokens = bond_data{k};
        pairs = [1 2; 2 3; 3 4];
        for p = 1:size(pairs,1)
            idx1 = find(allPdNLabels == tokens{pairs(p,1)});
            idx2 = find(allPdNLabels == tokens{pairs(p,2)});
            c1 = allCoords(idx1, :); c2 = allCoords(idx2, :);
            
            type1 = allAtomTypes{idx1}; type2 = allAtomTypes{idx2};
            
            if (strcmp(type1, 'Pd') && strcmp(type2, 'N')) || (strcmp(type1, 'N') && strcmp(type2, 'Pd'))
                plot3([c1(1), c2(1)], [c1(2), c2(2)], [c1(3), c2(3)], 'k-', 'LineWidth', 1.5);
            elseif strcmp(type1, 'N') && strcmp(type2, 'N')
                plot3([c1(1), c2(1)], [c1(2), c2(2)], [c1(3), c2(3)], 'r-', 'LineWidth', 1);
            end
        end
    end
end

% Label all atoms
for i = 1:size(allCoords, 1)
    coord = allCoords(i, :);
    verticalOffset = 1.2;  % Offset Pd
    if strcmp(allAtomTypes{i}, 'N')
        verticalOffset = 0.7;  % Offset N
    end
    text(coord(1), coord(2), coord(3) + verticalOffset, allPdNLabels(i), ...
        'FontSize', 8, 'FontWeight', 'bold', 'Color', 'black', ...
        'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom');
end

% Apply lighting and material settings
xlim(view_limits(1,:)); ylim(view_limits(2,:)); zlim(view_limits(3,:));
daspect([1 1 1]); 
camlight('right'); camlight('left'); camlight('headlight');
lighting gouraud;
material([0.4 0.6 0.5 40 1]);  
set(gca, 'XGrid','off','YGrid','off','ZGrid','off','Box','on',...
         'XTickLabel',[],'YTickLabel',[],'ZTickLabel',[]);
xlabel('X'); ylabel('Y'); zlabel('Z');
title([nameWithoutExt ' Skeletal'], 'FontSize', 11, 'FontWeight', 'bold');
view(view_azimuth, view_elevation);

%% ==============================================================
%% FIGURE 3: Full structure
%% ==============================================================

% Filter out Hydrogens
nonH_idx = elements_raw ~= "H";
filtered_elements = elements_raw(nonH_idx);
filtered_coords = coords_raw(nonH_idx, :);
atomTypes = cellstr(filtered_elements);

% Colour map
colors_full = containers.Map({'Ni', 'Pd', 'Pt', 'P', 'C', 'F', 'N', 'O', 'S', 'B', 'Cl', 'Br', 'I'}, ...
    {[0.000 0.412 0.522], [0.000 0.412 0.522], [0.000 0.412 0.522], ...
     [1 0.5 0], [0.3 0.3 0.3], [0.216, 0.784, 0.671], ...
     [0.188, 0.314, 0.973], [1, 0, 0], [1, 1, 0], ...
     [1, 1, 0], [0, 1, 0], [0.5, 0, 0], [0.58, 0, 0.82]});

fig3 = figure('Name', '\3: Actual', 'Color', [1 1 1]);
ax3 = axes('Parent', fig3);
hold(ax3, 'on');

% Batch plot atoms
unique_types = unique(atomTypes);
for t = 1:length(unique_types)
    type = unique_types{t};
    if isKey(colors_full, type)
        type_indices = find(strcmp(atomTypes, type));
        col = colors_full(type);
        
        radius = baseRadius;
        if any(strcmp(type, {'Ni', 'Pd', 'Pt'}))
            radius = radius_pd_full;  % Larger radius for metals in Figure 3
        end
        
        for i = 1:length(type_indices)
            idx = type_indices(i);
            c = filtered_coords(idx, :);
            surf(radius*sx + c(1), radius*sy + c(2), radius*sz + c(3), ...
                'FaceColor', col, 'EdgeColor', 'none');
        end
    end
end

% Bond calculation
n = size(filtered_coords, 1);
bond_pairs = [];

for i = 1:n
    type_i = atomTypes{i};
    if any(strcmp(type_i, {'Ni', 'Pd', 'Pt'}))   % For each atom of this type, connect it to its
        maxDist = 2.5; maxBonds = 4; 
    elseif ismember(type_i, {'P','C','N','O','F','S','Cl','Br','I'})
        maxDist = 2.0; maxBonds = 3; 
    else
        continue;
    end
    
    distances = sqrt(sum((filtered_coords - filtered_coords(i,:)).^2, 2));
    distances(i) = inf;
    
    potential_bonds = find(distances <= maxDist);
    [~, sorted_idx] = sort(distances(potential_bonds));
    bonds = potential_bonds(sorted_idx(1:min(maxBonds, length(sorted_idx))));
    
    for j = 1:length(bonds)
        bond_pairs = [bond_pairs; i, bonds(j)]; %#ok<AGROW>
    end
end

% Remove duplicates
bond_pairs = unique(sort(bond_pairs, 2), 'rows');

% Bond plotting
x_lines = []; y_lines = []; z_lines = [];
for i = 1:size(bond_pairs, 1)
    idx1 = bond_pairs(i, 1); idx2 = bond_pairs(i, 2);
    c1 = filtered_coords(idx1, :); c2 = filtered_coords(idx2, :);
    x_lines = [x_lines; c1(1); c2(1); NaN]; %#ok<AGROW>
    y_lines = [y_lines; c1(2); c2(2); NaN]; %#ok<AGROW>
    z_lines = [z_lines; c1(3); c2(3); NaN]; %#ok<AGROW>
end
plot3(x_lines, y_lines, z_lines, 'k-', 'LineWidth', 1.5);

% Apply lighting and material settings
xlim(view_limits(1,:)); ylim(view_limits(2,:)); zlim(view_limits(3,:));
daspect([1 1 1]); 
camlight('right'); camlight('headlight');
lighting gouraud;
material([0.4 0.6 0.5 40 1]);  
set(gca, 'XGrid','off','YGrid','off','ZGrid','off','Box','on',...
         'XTickLabel',[],'YTickLabel',[],'ZTickLabel',[]);
xlabel('X'); ylabel('Y'); zlabel('Z');
title([nameWithoutExt ' Actual'], 'FontSize', 11, 'FontWeight', 'bold');
view(view_azimuth, view_elevation);

hold(ax1, 'off');

%% ==============================================================
%% Add fake fog to all figures
%% ==============================================================

% Fog parameters
fogStartFraction = 0.25;  % Start position as fraction of bounding box depth
fogEndFraction = 1.0;     % End position as fraction of bounding box depth
nPlanes = 16;             % Number of fog planes
planeAlpha = 0.08;        % Opacity of each fog plane


% Only apply fog if enabled
if strcmpi(fog, 'on')		
    % Store figure and axes handles
    figHandles = [fig1, fig2, fig3];
    axHandles = [ax1, ax2, ax3];
    
    % Loop through all figures to add fog planes
    for figIdx = 1:3
        figure(figHandles(figIdx));
        ax = axHandles(figIdx);
        
        % Save current view
        [currentAz, currentEl] = view(ax);
        
        % Get axes limits
        xL = xlim(ax);
        yL = ylim(ax);
        zL = zlim(ax);
        
        % Camera information
        camPos = get(ax, 'CameraPosition');
        camTarget = get(ax, 'CameraTarget');
        camDir = camTarget - camPos;
        camDir = camDir / norm(camDir); 
        
        % Calculate plane normal (facing the camera)
        planeNormal = camDir;
        
        % Find a vector perpendicular to plane normal for plane orientation
        if abs(planeNormal(1)) > 0.5
            tempVec = [0, 1, 0];
        else
            tempVec = [1, 0, 0];
        end
        
        planeRight = cross(planeNormal, tempVec);
        planeRight = planeRight / norm(planeRight);
        planeUp = cross(planeRight, planeNormal);
        planeUp = planeUp / norm(planeUp);
        
        % Bounding box
        bboxCenter = [(xL(1)+xL(2))/2, (yL(1)+yL(2))/2, (zL(1)+zL(2))/2];
        bboxSize = [diff(xL), diff(yL), diff(zL)];
        maxDist = norm(bboxSize);
        
        % Plane must be big enough to cover view
        planeHalfSize = max(bboxSize) * 1.5;
        
        % Calculate start and end positions along camera direction
        startPos = bboxCenter + camDir * maxDist * (fogStartFraction - 0.5);
        endPos = bboxCenter + camDir * maxDist * (fogEndFraction - 0.5);
        
        % Create fog planes
        hold(ax, 'on');
        
        for i = 1:nPlanes
            % Compute plane position
            t = (i-1) / (nPlanes-1);  % 0 to 1
            planeCenter = startPos + t * (endPos - startPos);
            
            % Create plane vertices
            [u, v] = meshgrid([-planeHalfSize, planeHalfSize], [-planeHalfSize, planeHalfSize]);
            
            % Convert to 3D coordinates
            Xplane = planeCenter(1) + u * planeRight(1) + v * planeUp(1);
            Yplane = planeCenter(2) + u * planeRight(2) + v * planeUp(2);
            Zplane = planeCenter(3) + u * planeRight(3) + v * planeUp(3);
            
            % Draw plane
            surf(Xplane, Yplane, Zplane, ...
                'FaceColor', [1 1 1], ...
                'FaceAlpha', planeAlpha, ...
                'EdgeColor', 'none', ...
                'FaceLighting', 'none');
        end
        
        % Ensure planes render behind atoms/bonds
        fogPlanes = findobj(ax, 'Type', 'Surface', 'FaceColor', [1 1 1]);
        if ~isempty(fogPlanes)
            for i = 1:length(fogPlanes)
                uistack(fogPlanes(i), 'bottom');
            end
        end
        
        % Restore original view
        view(ax, currentAz, currentEl);
        hold(ax, 'off');
    end
end
%% ==============================================================
elapsedTime = toc(plottingTimer);
fprintf('...completed in %.2f seconds\n\n', elapsedTime);
end % end of plotting section