% =========================================================================
% 3D Impedance Surface - Series 2 Empirical Data
% =========================================================================
clear; clc; close all;
% --- 1. User Configuration ---
baseDir = 'C:\Users\anupa\Desktop\My files\PhD\experiments\pre_meas_data_2025\coating_exp_w_drying';
folderPath = fullfile(baseDir, 'Series_2');
% Series 2 Metadata
s2_names = {'C_Co_V1_t00'; 'C_Co_V1_t01'; 'C_Co_V1_t02'; 'C_Co_V1_t03'; 'C_Co_V1_t04'; 'C_Co_V1_t05'; 'C_Co_V1_t07'};
s2_moist = [0; 70; 66.2300885; 57.78761062; 40.09419152; 17.46214852; 9.858267717];
% Preallocate master arrays
all_F = []; all_M = []; all_Z = []; all_Phase = [];
% --- 2. Data Extraction ---
for i = 1:length(s2_names)
    fname = sprintf('%s.txt', s2_names{i});
    fullFile = fullfile(folderPath, fname);
    
    if isfile(fullFile)
        % Use the provided custom extraction function
        [f_m, mag, f_p, phase] = extractBodeData(fullFile);
        
        if ~isempty(f_m)
            % Clean duplicate frequencies to prevent interp1 crashes
            [f_p_unq, idx_unq] = unique(f_p);
            phase_unq = phase(idx_unq);
            
            % Align phase to magnitude frequencies
            phase_aligned = interp1(f_p_unq, phase_unq, f_m, 'linear', 'extrap');
            
            % Append to master arrays
            all_F = [all_F; f_m];
            all_Z = [all_Z; mag];
            all_Phase = [all_Phase; phase_aligned];
            all_M = [all_M; repmat(s2_moist(i), length(f_m), 1)];
        end
    else
        fprintf('File not found: %s\n', fullFile);
    end
end
if isempty(all_F)
    error('Master arrays are empty. Check file paths.');
end
% --- 3. FIXED: Structured Grid Interpolation ---
f_common = logspace(log10(min(all_F)), log10(max(all_F)), 200);
n_samples = length(s2_names);
Z_matrix     = NaN(n_samples, length(f_common));
Phase_matrix = NaN(n_samples, length(f_common));
% Sort samples by moisture for clean surface mesh
[m_sorted, sort_idx] = sort(s2_moist);
for k = 1:n_samples
    i = sort_idx(k);
    % Pull only the data rows belonging to this sample
    mask = abs(all_M - s2_moist(i)) < 1e-9;  
    f_s  = all_F(mask);
    z_s  = all_Z(mask);
    ph_s = all_Phase(mask);
    if length(f_s) < 2, continue; end
    [f_s, si] = sort(f_s);
    z_s  = z_s(si);
    ph_s = ph_s(si);
    % Interpolate each sample onto the common frequency grid
    Z_matrix(k,:)     = interp1(log10(f_s), log10(z_s), log10(f_common), ...
                                'linear', 'extrap');
    Phase_matrix(k,:) = interp1(log10(f_s), ph_s,       log10(f_common), ...
                                'linear', 'extrap');
end
[F_grid, M_grid] = meshgrid(f_common, m_sorted);

% --- 4. Render Plot ---
% Increased figure width from 900 to 1100 to accommodate huge fonts
figure('Color', 'w', 'Position', [100 100 1100 650]);
surf(F_grid, M_grid, Z_matrix, Phase_matrix);
shading interp;
colormap(jet);
clim([-90 0]);
view(135, 30);

% Manually squeeze the main plot to the left to make room for the colorbar
% [left bottom width height]
set(gca, 'Position', [0.12 0.15 0.60 0.75]); 

% formatting - Font size increased to 16 for better readability
set(gca, 'XScale', 'log', 'FontName', 'Times New Roman', 'FontSize', 16);
xlabel('Frequency (Hz)', 'FontSize', 18);
ylabel('Moisture (wt-%)', 'FontSize', 18);
zlabel('|Z| (\Omega)', 'FontSize', 18);
zlim([0 8]);
set(gca, 'XTick', [1e3 1e4 1e5 1e6], 'XTickLabel', {'10^3','10^4','10^5','10^6'});
set(gca, 'ZTick', 0:2:8, 'TickLabelInterpreter', 'tex');
set(gca, 'ZTickLabel', {'10^0','10^2','10^4','10^6','10^8'});

c = colorbar;
ylabel(c, 'Phase Angle (\circ)', 'FontSize', 18);
c.Ticks = [-90 -45 0];
c.TickLabels = {'-90 (Capacitive)', '-45 (Transition)', '0 (Resistive)'};
c.FontSize = 18; % Colorbar tick labels
grid on; box on;

% --- 5. Add Text Annotations ---
% Coordinates: text(Frequency, Moisture, log10(Impedance), 'String')

% Dry State: Placed on the high-impedance, low-moisture plane (back left)
% We bump Z up to 7.5 so it clears the surface mesh and faces the camera.
text(3e5, 0, 7, 'Dry State', 'Color', 'b', 'FontWeight', 'bold', ...
    'FontSize', 16, 'FontName', 'Times New Roman', ...
    'HorizontalAlignment', 'center');

% Wet State: Placed on the low-impedance, high-moisture plateau (front red area)
% We center it at 10^5 Hz and bump Z to 3.5 so it floats clearly above the plateau.
text(1e6, 70, 1.1, 'Wet State', 'Color', 'r', 'FontWeight', 'bold', ...
    'FontSize', 16, 'FontName', 'Times New Roman', ...
    'HorizontalAlignment', 'center');

% =========================================================================
% 6. HELPER FUNCTION
% =========================================================================
function [f_m, mag, f_p, phase] = extractBodeData(filepath)
    f_m = []; mag = []; f_p = []; phase = [];
    try
        opts = detectImportOptions(filepath);
        opts.CommentStyle = '%';
        opts.Delimiter = ';';
        opts.DataLines = [1 Inf];
        opts.VariableNamingRule = 'preserve';
        
        raw = readmatrix(filepath, opts);
        if isempty(raw), return; end
        
        freq = raw(:,1); 
        val = raw(:,2);
        
        % Slice the array precisely in half
        split_idx = floor(length(freq) / 2);
        
        f_m = freq(1:split_idx);
        mag = val(1:split_idx);
        
        f_p = freq(split_idx+1:end);
        phase = val(split_idx+1:end);
        
        % Sort frequencies ascending to prevent lines drawing backwards
        [f_m, idx_m] = sort(f_m);
        mag = mag(idx_m);
        
        [f_p, idx_p] = sort(f_p);
        phase = phase(idx_p);
    catch
        % Silent catch for clean output
    end
end