%% MASTER SCRIPT: All series merged — single calibration plot
clear; clc; close all;

% =========================================================================
% 1. USER CONFIGURATION
% =========================================================================
baseDir    = 'C:\Users\anupa\Desktop\My files\PhD\experiments\data_v_2025';
targetFreq = 100000; % 100 kHz

% =========================================================================
% 2. METADATA — all series merged into one
% =========================================================================

% Former Series 1 — only 1 clean point remains
names_all  = {'C_Co_V1_N1'};
thick_all  = [0.120];
cond_all   = [0];
folder_all = {'Series_1'};

% Former Series 2 — removed C_Ca_V1_N2 cond 1
names_all  = [names_all;  {'C_Co_V1_N0'; ...
                            'C_Ca_V1_N1';'C_Co_V1_N0'}];
thick_all  = [thick_all;   0.429; 0.445;0.429];
cond_all   = [cond_all;    1; 2;2];
folder_all = [folder_all; repmat({'Series_2'}, 3, 1)];

% Former Series 3 — valid points only (excludes 1.830mm delamination)
names_all  = [names_all;  {'C_Co_V1_N0_1';'C_Co_V1_N0_2';'C_Co_V1_N0_3';'C_Co_V1_N0_4'}];
thick_all  = [thick_all;   0.344;0.645;0.901;1.147];
cond_all   = [cond_all;    1;2;3;4];
folder_all = [folder_all; repmat({'Series_3'}, 4, 1)];

% Delamination point — plotted separately as excluded
delam_name   = 'C_Co_V1_N0_5';
delam_thick  = 1.830;
delam_cond   = 5;
delam_folder = 'Series_3';

% =========================================================================
% 3. EXTRACT ALL DATA
% =========================================================================
all_thick = []; all_mag = [];

for i = 1:length(names_all)
    if strcmp(folder_all{i}, 'Series_3')
        fname = sprintf('%s.txt', names_all{i});
    else
        fname = sprintf('%s_%d.txt', names_all{i}, cond_all(i));
    end
    fullFile = fullfile(baseDir, folder_all{i}, fname);
    if isfile(fullFile)
        m = extractData(fullFile, targetFreq);
        if ~isnan(m)
            all_thick(end+1,1) = thick_all(i);
            all_mag(end+1,1)   = m;
        end
    else
        fprintf('Not found: %s\n', fullFile);
    end
end

% Extract delamination point
delam_file = fullfile(baseDir, delam_folder, sprintf('%s.txt', delam_name));
delam_mag  = NaN;
if isfile(delam_file)
    delam_mag = extractData(delam_file, targetFreq);
end

if isempty(all_thick)
    error('No data loaded. Check file paths.');
end

% =========================================================================
% 4. S3 SLOPE — computed from S3 points only (thick >= 0.344mm)
% =========================================================================
s3_idx   = (all_thick >= 0.344);
p_s3     = polyfit(all_thick(s3_idx), all_mag(s3_idx), 1);
s3_slope = p_s3(1);

% =========================================================================
% 5. CENTROID LINE — S3 slope, intercept through centroid of all points
% =========================================================================
intercept = mean(all_mag) - s3_slope * mean(all_thick);

yfit   = s3_slope * all_thick + intercept;
ss_res = sum((all_mag - yfit).^2);
ss_tot = (length(all_mag)-1) * var(all_mag);
rsq    = 1 - ss_res/ss_tot;

fprintf('\n--- Fit Results ---\n');
fprintf('Slope     : %.4f MOhm/mm\n', s3_slope/1e6);
fprintf('Intercept : %.4f MOhm\n',    intercept/1e6);
fprintf('R²        : %.4f\n',         rsq);
fprintf('N points  : %d\n',           length(all_thick));

% =========================================================================
% 6. PLOT
% =========================================================================
figure('Color','w','Position',[100 100 820 600]);
hold on; grid on;

% All valid points
scatter(all_thick, all_mag/1e6, 80, [0.15 0.45 0.80], ...
    'filled', 'MarkerEdgeColor','k', 'LineWidth',0.8, ...
    'DisplayName','Valid Data (Linear Regime)');

% Delamination point
if ~isnan(delam_mag)
    scatter(delam_thick, delam_mag/1e6, 100, 'rx', ...
        'LineWidth', 2.5, 'DisplayName','Excluded Data (Delamination)');
end

% Linear fit line
x_line = linspace(0, 1.3, 200);
y_line = s3_slope * x_line + intercept;
plot(x_line, y_line/1e6, 'r--', 'LineWidth', 2.0, ...
    'DisplayName', sprintf('Linear Fit (R^2 = %.2f)', rsq));

set(gca, 'FontName','Times New Roman', 'FontSize',16);
xlabel('Coating Thickness (mm)',        'FontName','Times New Roman','FontSize',16);
ylabel('Impedance Magnitude (M\Omega)', 'FontName','Times New Roman','FontSize',16);
legend('Location','best','FontName','Times New Roman','FontSize',16);
xlim([0 2.0]);
ylim([1.25 1.35]);
hold off;

% =========================================================================
% 7. HELPER — must be at end of file
% =========================================================================
function magVal = extractData(filepath, freqTarget)
    magVal = NaN;
    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);
        split_idx = floor(length(freq)/2);
        f_m = freq(1:split_idx);
        mag = val(1:split_idx);

        [f_m, idx] = sort(f_m);
        mag = mag(idx);
        [f_m, u]   = unique(f_m);
        mag = mag(u);
        if length(f_m) < 2, return; end

        magVal = interp1(f_m, mag, freqTarget, 'linear', NaN);
    catch
    end
end