clear all; close all;

% -----------------------------------------
% Your current files + frequencies
% -----------------------------------------
files = {
    "0p1.csv"   % 0.1 Hz
    "0p22.csv"  % 0.2 Hz
    "0p3.csv"   % 1 Hz
    "0p4.csv"    % 10 Hz
    "0p5.csv"    % 20 Hz
};

freqs = [0.1 0.2 1 10 20];

Fs = 960;   % WinDaq sampling rate

numF = length(freqs);
phase_deg = zeros(1, numF);

% -----------------------------------------
% MAIN LOOP — PHASE ONLY
% -----------------------------------------
for k = 1:numF

    filename = files{k};
    f = freqs(k);

    fprintf("\nProcessing %.2f Hz from file %s\n", f, filename);

    % Load data
    data = readmatrix(filename);
    V = data(:,2);      % voltage only
    N = length(V);

    % Time vector
    t = (0:N-1)'/Fs;

    % Remove DC
    V_ac = V - mean(V);

    % Sine model
    w = 2*pi*f;
    sine_model = @(b,t) b(1)*sin(w*t + b(2));

    % Fit VOLTAGE
    beta0 = [range(V_ac)/2, 0];
    beta = nlinfit(t, V_ac, sine_model, beta0);

    A_out = beta(1);
    phi_out = beta(2);

    % Store phase
    phase_deg(k) = rad2deg(phi_out);

    % Plot voltage fit for this frequency
    t_fit = linspace(min(t), max(t), 2000);
    V_fit = A_out * sin(w*t_fit + phi_out);

    figure;
    plot(t, V_ac, 'b'); hold on;
    plot(t_fit, V_fit, 'r', 'LineWidth', 2);
    xlabel("Time (s)");
    ylabel("Voltage (AC)");
    title(sprintf("%.2f Hz | Phase Lag = %.2f°", f, phase_deg(k)));
    legend("Measured Voltage", "Voltage Fit");
    grid on;

end

% -----------------------------------------
% FINAL PHASE PLOT
% -----------------------------------------
figure;
semilogx(freqs, phase_deg, 'o-', 'LineWidth', 2);
xlabel("Frequency (Hz)");
ylabel("Phase (degrees)");
title("Phase vs Frequency (Voltage Only)");
grid on;
