%% Plot cavity signal
% Load the data from the CSV file
ScopeData    = readmatrix('.csv');

% Extract Time and CavitySignal with offsets
Time         = ScopeData(:, 1);
CavitySignal = ScopeData(:, 2);

% Calculate xoffset and yoffset
xoffset      = Time(1);
yoffset      = min(CavitySignal);

% Adjust CavitySignal and Time
CavitySignal = CavitySignal - yoffset;
Time         = Time - xoffset;

% Plot the data
figure;
scatter(Time, CavitySignal, 10, 'b', 'filled', 'MarkerFaceAlpha', 0.5);
grid on;

% Format the plot
xlabel('\bf Time (s)', 'FontSize', 16);
ylabel('\bf Voltage (V)', 'FontSize', 16);
title('\bf Cavity Signal', 'FontSize', 16);
set(gca, 'FontSize', 12);

%% Fit cavity signal

% Extract signal cavity mode
tstartIdx = 1;
tendIdx   = 50;
TruncatedScopeData = ScopeData(tstartIdx:tendIdx,:);

% Fit and plot Airy function, extracting characteristic parameters
fitAndplotCavityMode(TruncatedScopeData);


%% Extract distance between consecutive cavity modes and their amplitudes

% == Add two Airy functions to fit two consecutive cavity modes == %

function fitAndplotCavityMode(dataset)
    % Define the Airy function
    AiryFunc     = @(a, b, t) a ./ (1 + b * (sin(t / 2)).^2);

    % Perform non-linear fitting to find parameters a and b
    t            = dataset(:, 1);
    CavitySignal = dataset(:, 2);
    fitParams    = fit(t, CavitySignal, @(a, b, t) AiryFunc(a, b, t), ...
                    'StartPoint', [1, 1]);
    a            = fitParams.a;
    b            = fitParams.b; % Coefficient of finesse from fit

    % Calculate reflectivity (r)
    syms r;
    Reflectivity = solve(b == (4 * r) / (1 - r)^2, r);
    
    % Calculate finesse from reflectivity (r)
    Finesse      = (pi * sqrt(Reflectivity))/(1 - Reflectivity);

    % Generate fitted data
    fitData      = [t, AiryFunc(a, b, t)];

    % Find FWHM
    maxSignal    = max(fitData(:, 2));
    left         = find(fitData(:, 2) >= maxSignal / 2, 1, 'first');
    right        = find(fitData(:, 2) >= maxSignal / 2, 1, 'last');
    FWHM         = fitData(right, 1) - fitData(left, 1);
    
    % Calculate FSR from Finesse and FWHM
    FSR          = Finesse * FWHM;

    % Plot the original data and the fitted curve
    figure;
    plot(dataset(:, 1), dataset(:, 2), 'o', 'MarkerSize', 5, ...
         'MarkerEdgeColor', 'blue', 'MarkerFaceColor', 'blue', ...
         'DisplayName', 'Cavity Signal', 'MarkerFaceAlpha', 0.5);
    hold on;
    plot(fitData(:, 1), fitData(:, 2), '--', 'LineWidth', 1.5, ...
         'Color', [1, 0.5, 0], 'DisplayName', 'Best Fit');
    hold off;

    % Customize the plot
    grid on;
    xlabel('\bf Time (s)', 'FontSize', 16);
    ylabel('\bf Voltage (V)', 'FontSize', 16);
    title('\bf Airy Function Fit', 'FontSize', 16);
    legend('show', 'Location', 'Best');

    % Annotate the plot with fit results
    annotation('textbox', [0.6, 0.7, 0.3, 0.2], ...
               'String', sprintf(['Reflectivity = %.1f\n', ...
                                  'Finesse = %.1f\n', ...
                                  'FWHM = %.2f \n', ...
                                  'FSR = %.2f'], ...
                                 Reflectivity, Finesse, FWHM, FSR), ...
               'EdgeColor', 'none', ...
               'FontSize', 12);    
end