clear
clc
close all

n = 3;  % Total number of files
max = [];
format long  

for i = 3:n
    disp([num2str(i) '.csv'])
    
    % Load data
    filename = sprintf('%d.csv', i);
    data = readtable(filename);  % Use readtable to read CSV
    
    num1 = 1;   % Range of the data
    num2 = 9999;  % Range of the data
    
    ff = data{num1:num2, 1} / 10^8;
    zz = 50 * (1 - data{num1:num2, 2}.^2 - data{num1:num2, 3}.^2) ./ ((1 - data{num1:num2, 2}).^2 + data{num1:num2, 3}.^2);
    
    % Fitting
    ft = fittype('a*f^4/(f^2*fs^2/q^2+(f^2-fs^2)^2)+b+c*f', ...
        'dependent', {'z'}, 'independent', {'f'}, ...
        'coefficients', {'a', 'b', 'c', 'fs', 'q'});
    
    f = fit(ff, zz, ft, 'StartPoint', [0, 0.2, 0, 4.5, 30]);
    coeff = coeffvalues(f);
    
    % Find maximum
    [~, locs] = findpeaks(f(ff));
    
    max(i, 1) = i;
    max(i, 2) = ff(locs(1)) * 10^8;
    max(i, 3) = f(ff(locs(1)));
    max(i, 4) = coeff(5);
    
    % Plot results
    array = [1, 2, 3, 4, 5, 6, 7, 8, 22, 1000];
    if ismember(i, array)
        figure
        plot(f, ff, zz)
        hold on
        plot(max(i, 2) / 10^8, max(i, 3), 'o', 'Color', 'g', 'LineWidth', 1.5)
        legend('Data', 'Maximum Point', 'Location', 'Northwest')
        grid on
        hold off
    end
    
    % Save frequency and Re Z to Excel
    dataToSave = [data{num1:num2, 1} / 10^6, real(zz)]; % Combine frequency and Re Z for saving
    writematrix(dataToSave, ['PBS_Frequency_ReZ_' num2str(i) '.xlsx']);
end

% Save maximum data to CSV
csvwrite(['Maximum_1-' num2str(n) '.csv'], max);

% Save maximum data to ASCII tab-delimited file
save('Maximum.dat', 'max', '-ascii', '-tabs');