AllFigScript.m

% UNIFIED MATLAB SCRIPT FOR MANUSCRIPT FIGURES 3, 4, & 5
%
% Date: Dec 9, 2025
% Author: Noah Muscat
%
% Description:
% This master script combines the logic from three separate scripts to
% generate all components for Figures 3, 4, and 5. It centralizes the
% configuration, loads and processes the data only once, and then
% generates each figure in a clearly defined section. All figures are
% saved as PNG files.
%
% Script Sections:
% 1. Configuration & Setup
% 2. Load and Preprocess Data
% 3. Generate Figure 3 Components
% 4. Generate Figure 4 Components
% 5. Generate Figure 5 Components
% 6. Helper Functions
%
% ASSUMPTIONS: violinplot.m is on the MATLAB path.
% ---------------------------------------------------------------------
%% --- 1. Configuration & Setup ---
clear; clc; close all; % Start fresh
% --- File and Path Parameters ---
filename = '/Users/noahmuscat/University of Michigan Dropbox/Noah Muscat/ActivityAnalysis/ActigraphyOnly/Unified_AO1to12_Combined_Data.csv';
savePath_Fig3 = '/Users/noahmuscat/Desktop/WatsonLab/AOPaper/Figure3Plots';
savePath_Fig4 = '/Users/noahmuscat/Desktop/WatsonLab/AOPaper/Figure4Plots';
savePath_Fig5 = '/Users/noahmuscat/Desktop/WatsonLab/AOPaper/Figure5Plots';
% --- Master list of Activity Metrics to Plot ---
activity_metrics_to_plot = {'SelectedPixelDifference', 'NormalizedActivity'};
metric_suffixes = {'PixelDiff', 'Normalized'};
metric_ylabels = {'Activity (Pixel Difference)', 'Normalized Activity'};
pixelDiffVar = 'SelectedPixelDifference';
normalizedActivityVar = 'NormalizedActivity';
% --- Key Data Column Names ---
animalVar = 'Animal';
conditionVar = 'Condition';
relativeDayVar = 'RelativeDay';
timeVarZT = 'DateZT';
ztHourVar = 'ZT_Time';
% --- Analysis Parameters ---
samplingIntervalMinutes = 5;
hoursPerDay = 24;
daysForPoolingFig3AB = 14;
targetConditionFig3AB = '300Lux';
ztLightsOnStart = 0; ztLightsOnEnd = 11;
ztLightsOffStart = 12; ztLightsOffEnd = 23;
% --- Condition Ordering ---
allConditions = {'300Lux', '1000Lux', 'FullDark', '300LuxEnd'};
% --- Animal Grouping ---
allAnimals = string({'AO1', 'AO2', 'AO3', 'AO4', 'AO5', 'AO6', 'AO7', 'AO8', 'AO9', 'AO10', 'AO11', 'AO12'});
maleAnimals = string({'AO1', 'AO2', 'AO3', 'AO7', 'AO9', 'AO12'});
femaleAnimals = string({'AO4', 'AO5', 'AO6', 'AO8', 'AO10', 'AO11'});
fourCondAnimals = string({'AO5', 'AO6', 'AO7', 'AO8', 'AO9','AO10','AO11','AO12'});
% --- Color Definitions ---
maleColor = [0.2 0.5470 0.8410];
femaleColor = [0.9500 0.4250 0.1980];
grayColor = [0.6 0.6 0.6];
% --- Setup Save Directories ---
all_paths = {savePath_Fig3, savePath_Fig4, savePath_Fig5};
for p_idx = 1:length(all_paths)
    current_path = all_paths{p_idx};
    disp(['Checking save directory: ', current_path]);
    if ~isfolder(current_path)
        try mkdir(current_path); disp('Save directory created.');
        catch ME, error('SaveDirErr:CouldNotCreate', 'Could not create save directory "%s". Err: %s', current_path, ME.message); end
    end
end
%% --- 2. Load and Preprocess Data ---
disp(['Loading data from: ', filename]);
try
    opts = detectImportOptions(filename);
    opts.VariableNamingRule = 'preserve';
    expectedVars = {pixelDiffVar, normalizedActivityVar, animalVar, conditionVar, relativeDayVar, timeVarZT, ztHourVar};
    if ~all(ismember(expectedVars, opts.VariableNames)), error('LoadErr:MissingColumn', 'One or more expected columns not found in %s.', filename); end
    
    opts = setvartype(opts, {pixelDiffVar, normalizedActivityVar, relativeDayVar, ztHourVar}, 'double');
    opts = setvartype(opts, {animalVar, conditionVar}, 'string');
    opts = setvartype(opts, timeVarZT, 'datetime');
    opts = setvaropts(opts, timeVarZT, 'InputFormat', 'dd-MMM-yyyy HH:mm:ss');
    
    dataTable = readtable(filename, opts);
    disp('Data loaded successfully.');
catch ME_load, error('LoadErr:FileProcessing', 'Failed to load/parse CSV "%s". Err: %s', filename, ME_load.message); end
disp('Starting data preprocessing...');
dataTable = sortrows(dataTable, timeVarZT);
dataTable.DayOfCondition = floor(dataTable.(relativeDayVar));
colsToClean = {pixelDiffVar, normalizedActivityVar, timeVarZT, relativeDayVar, ztHourVar, animalVar, conditionVar, 'DayOfCondition'};
nanRowsFilter = false(height(dataTable), 1);
for c_idx = 1:length(colsToClean)
    colData = dataTable.(colsToClean{c_idx});
    if isdatetime(colData), nanRowsFilter = nanRowsFilter | isnat(colData);
    elseif isnumeric(colData), nanRowsFilter = nanRowsFilter | isnan(colData);
    elseif isstring(colData), nanRowsFilter = nanRowsFilter | ismissing(colData);
    end
end
if any(nanRowsFilter), fprintf('%d rows removed due to missing values.\n', sum(nanRowsFilter)); dataTable(nanRowsFilter, :) = []; end
if isempty(dataTable), error('PreprocErr:NoData', 'No valid data after NaN removal.'); end
dataTable = dataTable(ismember(dataTable.(animalVar), allAnimals), :);
if isempty(dataTable), error('PreprocErr:NoAnimalData', 'No data found for the animals specified in `allAnimals`.'); end
dataTable.Condition = categorical(dataTable.Condition, allConditions, 'Ordinal', true);
dataTable.Sex = categorical(ismember(dataTable.(animalVar), maleAnimals), [0, 1], {'Female', 'Male'});
dataTable.LightPeriod = categorical((dataTable.(ztHourVar) >= ztLightsOnStart & dataTable.(ztHourVar) <= ztLightsOnEnd), [0, 1], {'Off', 'On'});
disp('Data preprocessing complete.');
%% --- 3. GENERATE FIGURE 3 COMPONENTS ---
disp('=====================================================');
disp('--- Starting Generation of Figure 3 Components ---');
disp('=====================================================');
initialBaselineData = [];
for i_an = 1:length(allAnimals)
    currentAnimalID = allAnimals(i_an);
    animalData = dataTable(dataTable.(animalVar) == currentAnimalID, :);
    animalConds = unique(animalData.Condition, 'stable');
    if ~isempty(animalConds) && animalConds(1) == string(targetConditionFig3AB)
        animalTargetCondData = animalData(animalData.Condition == string(targetConditionFig3AB), :);
        if ~isempty(animalTargetCondData)
            uniqueDays = unique(animalTargetCondData.DayOfCondition);
            if isempty(uniqueDays), continue; end
            daysToKeep = uniqueDays(max(1, end-daysForPoolingFig3AB+1):end);
            finalDataForPooling = animalTargetCondData(ismember(animalTargetCondData.DayOfCondition, daysToKeep), :);
            initialBaselineData = [initialBaselineData; finalDataForPooling];
        end
    end
end
if isempty(initialBaselineData), error('DataErr:Fig3_NoBaseline', 'No baseline data found for Figure 3.'); end
for met_idx = 1:length(activity_metrics_to_plot)
    current_metric_var = activity_metrics_to_plot{met_idx};
    current_metric_suffix = metric_suffixes{met_idx};
    current_metric_ylabel = metric_ylabels{met_idx};
    fprintf('\n--- Generating Figure 3 components for metric: %s ---\n', current_metric_var);
    
    %% FIG 3A: Averaged Actogram
    disp('--- Starting Figure 3A: Averaged Actogram ---');
    try
        alignedAnimalData = NaN(daysForPoolingFig3AB, hoursPerDay, length(allAnimals));
        for i_an = 1:length(allAnimals)
            currentAnimalID = allAnimals(i_an);
            animalData = dataTable(dataTable.(animalVar) == currentAnimalID & dataTable.Condition == string(targetConditionFig3AB), :);
            if isempty(animalData), continue; end
            uniqueDays = unique(animalData.DayOfCondition);
            daysToKeep = uniqueDays(max(1, end-daysForPoolingFig3AB+1):end);
            for i_day_idx = 1:min(length(daysToKeep), daysForPoolingFig3AB)
                dayData = animalData(animalData.DayOfCondition == daysToKeep(i_day_idx), :);
                if ~isempty(dayData)
                    hourlyMeans = groupsummary(dayData, ztHourVar, 'mean', current_metric_var);
                    [lia, locb] = ismember(hourlyMeans.(ztHourVar), (0:23)');
                    alignedAnimalData(i_day_idx, locb(lia), i_an) = hourlyMeans.(['mean_', current_metric_var])(lia);
                end
            end
        end
        actogramMatrix = mean(alignedAnimalData, 3, 'omitnan');
        
        hFig3A_current = figure('Name', sprintf('Fig 3A: Averaged Actogram (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w');
        ax3A = axes('Parent', hFig3A_current);
        imagesc(ax3A, (0:23), (1:daysForPoolingFig3AB), actogramMatrix); colormap(ax3A, 'parula'); ax3A.YDir = 'reverse';
        xlabel(ax3A, 'ZT Hour'); ylabel(ax3A, sprintf('Day of %s (Aligned)', targetConditionFig3AB));
        ylim(ax3A,[0.5, daysForPoolingFig3AB + 0.5]); xticks(ax3A, [0,6,12,18,23]); yticks(ax3A, 2:2:daysForPoolingFig3AB);
        title(ax3A,sprintf('Averaged Actogram (%s)\nLast %d days of %s',strrep(current_metric_suffix,'_',' '),daysForPoolingFig3AB,targetConditionFig3AB));
        cb3A = colorbar(ax3A); ylabel(cb3A, current_metric_ylabel);
        hold(ax3A,'on'); patch(ax3A,[ztLightsOffStart-0.5, ztLightsOffEnd+0.5, ztLightsOffEnd+0.5, ztLightsOffStart-0.5],[0.5,0.5,daysForPoolingFig3AB+0.5,daysForPoolingFig3AB+0.5],[0.85 0.85 0.85],'FaceAlpha',0.25,'EdgeColor','none','HandleVisibility','off'); hold(ax3A,'off');
        set(ax3A,'Layer','top', 'Position', [0.13 0.11 0.7 0.815]); % Add margins for whitespace
        fig_filename = fullfile(savePath_Fig3, sprintf('Figure3A_Actogram_Averaged_%s.png', current_metric_suffix));
        exportgraphics(hFig3A_current, fig_filename); disp(['Saved: ', fig_filename]); close(hFig3A_current);
    catch ME_3A, warning('ErrGen:Fig3A','Error Fig 3A (%s): %s', current_metric_suffix, ME_3A.message); end
    
    %% FIG 3B: Pooled Periodogram
    disp('--- Starting Figure 3B: Pooled Periodogram ---');
    try
        hFig3B_current = figure('Name', sprintf('Fig 3B: Periodogram (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w');
        ax3B = axes('Parent', hFig3B_current);
        tempDataForPerio = initialBaselineData;
        tempDataForPerio.RoundedTime5Min = dateshift(tempDataForPerio.(timeVarZT),'start','minute',samplingIntervalMinutes);
        pooled5MinStats = groupsummary(tempDataForPerio,'RoundedTime5Min','mean',current_metric_var);
        activitySignal = pooled5MinStats.(['mean_',current_metric_var]);
        timeSignalHours = hours(pooled5MinStats.RoundedTime5Min - pooled5MinStats.RoundedTime5Min(1));
        validIdx = ~isnan(activitySignal); activitySignal = activitySignal(validIdx); timeSignalHours = timeSignalHours(validIdx);
        periodogramRangeHours = [0.5, 48];
        [pxx, f_hz] = plomb(activitySignal, timeSignalHours, 1/periodogramRangeHours(1), 4);
        periodHours = 1 ./ f_hz;
        validDataIdx = isfinite(periodHours) & isfinite(pxx);
        periodHours = periodHours(validDataIdx); pxx = pxx(validDataIdx);
        displayRangeIdx = periodHours >= periodogramRangeHours(1) & periodHours <= periodogramRangeHours(2);
        [periodHours, sortIdx] = sort(periodHours(displayRangeIdx)); pxx = pxx(displayRangeIdx); pxx = pxx(sortIdx);
        plot(ax3B,periodHours,pxx,'LineWidth',1.5,'Color','k');
        xlabel(ax3B,'Period (hours)'); ylabel(ax3B,'Power');
        title(ax3B,sprintf('Periodogram (%s)\nLast %dd %s',current_metric_suffix,daysForPoolingFig3AB,targetConditionFig3AB));
        grid(ax3B,'on'); xlim(ax3B,periodogramRangeHours);
        set(ax3B, 'Position', [0.13 0.11 0.775 0.75]); % Add margins for whitespace
        fig_filename = fullfile(savePath_Fig3, sprintf('Figure3B_Periodogram_Last14d_%s.png', current_metric_suffix));
        exportgraphics(hFig3B_current, fig_filename); disp(['Saved: ', fig_filename]); close(hFig3B_current);
    catch ME_3B, warning('ErrGen:Fig3B','Error Fig 3B (%s): %s',current_metric_suffix, ME_3B.message); end
    
%% FIG 3C: Pooled 24-Hour Activity Profile (WITH HOURLY STATS)
    disp('--- Starting Figure 3C: Pooled 24-Hour Profile ---');
    try
        hFig3C_current = figure('Name', sprintf('Fig 3C: 24h Profile (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w');
        ax3C = axes('Parent', hFig3C_current);
        
        % Calculate Plotting Data (Mean +/- SEM)
        hourlyStats_3C = groupsummary(initialBaselineData, ztHourVar, {'mean','std','nnz'}, current_metric_var);
        hourlyStats_3C = sortrows(hourlyStats_3C, ztHourVar);
        meanAct_3C = hourlyStats_3C.(['mean_',current_metric_var]);
        semAct_3C  = hourlyStats_3C.(['std_',current_metric_var]) ./ sqrt(hourlyStats_3C.(['nnz_',current_metric_var]));
        
        % Plot the Profile
        errorbar(ax3C, hourlyStats_3C.(ztHourVar), meanAct_3C, semAct_3C, 'o-', 'LineWidth', 1.5, 'MarkerSize', 4, 'CapSize', 3, 'Color', [0.2 0.2 0.2]);
        hold(ax3C,'on');
        
        % Stats: Hour X vs. Average of Other 23 Hours
        disp('--- STATS: Fig 3C (Hour vs Daily Average) ---');
        unique_animals = unique(initialBaselineData.Animal);
        p_values_3C = ones(24, 1);
        
        % Calculate animal-level means for the t-test to avoid pseudoreplication
        % Create a matrix: Rows = Animals, Cols = Hours (0-23)
        animal_hourly_matrix = NaN(length(unique_animals), 24);
        
        for i = 1:length(unique_animals)
            an_data = initialBaselineData(initialBaselineData.Animal == unique_animals(i), :);
            an_hourly = groupsummary(an_data, ztHourVar, 'mean', current_metric_var);
            % Map to 0-23 index
            [found, idx] = ismember((0:23)', an_hourly.(ztHourVar));
            animal_hourly_matrix(i, found) = an_hourly.(['mean_', current_metric_var])(idx(found));
        end
        
        % No bonferroni
        alpha_3C = 0.05;
        
        for h = 1:24 % Loop 1 to 24 (corresponding to ZT 0 to 23)
            col_target = animal_hourly_matrix(:, h);
            col_others = animal_hourly_matrix; 
            col_others(:, h) = []; % Remove target column
            mean_others = mean(col_others, 2, 'omitnan'); % Average of the other 23 hours per animal
            
            % Paired t-test: Is this hour significantly different from the animal's background average?
            if all(~isnan(col_target) & ~isnan(mean_others))
                [~, p] = ttest(col_target, mean_others);
                p_values_3C(h) = p;
            end
        end
        
        % 4. Plot Asterisks
        % Find significant indices
        sig_idx = find(p_values_3C < alpha_3C);
        
        if ~isempty(sig_idx)
            % Determine Y-position for asterisks (slightly above the error bar)
            y_points = meanAct_3C(sig_idx) + semAct_3C(sig_idx);
            % Add a buffer so it doesn't touch the bar (5% of range)
            y_range = max(meanAct_3C) - min(meanAct_3C);
            y_asterisks = y_points + (y_range * 0.05);
            
            plot(ax3C, sig_idx-1, y_asterisks, '*', 'Color', 'k', 'MarkerSize', 6);
            fprintf('Significant Hours (p < %.4f): ZT %s\n', alpha_3C, num2str(sig_idx' - 1));
        end

        % 5. Standard Formatting
        yl3C=ylim(ax3C); 
        % Adjust Y-lim if asterisks are cut off
        if ~isempty(sig_idx)
            ylim(ax3C, [yl3C(1), max(yl3C(2), max(y_asterisks) + y_range*0.05)]);
            yl3C=ylim(ax3C);
        end
        
        patch(ax3C,[ztLightsOffStart, ztLightsOffEnd+1, ztLightsOffEnd+1, ztLightsOffStart], [yl3C(1) yl3C(1) yl3C(2) yl3C(2)], [0.85 0.85 0.85],'FaceAlpha',0.25,'EdgeColor','none');
        if strcmp(current_metric_var, normalizedActivityVar), plot(ax3C,xlim(ax3C),[0 0],'k--'); end
        hold(ax3C,'off');
        xlabel(ax3C,'ZT Hour'); ylabel(ax3C,['Mean ',current_metric_ylabel]);
        title(ax3C,sprintf('Pooled 24h Activity Profile (%s - Baseline)',current_metric_suffix));
        xticks(ax3C,0:3:23); xlim(ax3C,[-0.5,23.5]); grid(ax3C,'on');
        set(ax3C, 'Position', [0.13 0.11 0.775 0.815]); 
        fig_filename = fullfile(savePath_Fig3, sprintf('Figure3C_Profile_Baseline_%s.png', current_metric_suffix));
        exportgraphics(hFig3C_current, fig_filename); disp(['Saved: ', fig_filename]); close(hFig3C_current);
    catch ME_3C, warning('ErrGen:Fig3C','Error Fig 3C (%s): %s',current_metric_suffix, ME_3C.message); end
    
    %% FIG 3D: Diurnality Quantification (Violin Plots) (WITH STATS)
    disp('--- Starting Figure 3D: Diurnality Violins ---');
    fprintf('\n--- STATS: Figure 3D (Last 14 days) for %s ---\n', current_metric_suffix);
    try
        if ~exist('violinplot.m','file'),error('ViolinPlotNotFound:Fig3D','violinplot.m not found.');end
        hFig3D_current = figure('Name', sprintf('Fig 3D: Diurnality Violins (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w');
        ax3D = axes('Parent', hFig3D_current); hold(ax3D, 'on');
        time_windows = {[0 23], [0 11], [3 9], [12 23], [15 21], [10 14], [22 2]};
        window_labels = {'Full Day', 'Bright (0-11)', 'Mid-Bright (3-9)', 'Dark (12-23)', 'Mid-Dark (15-21)', 'Bright-Dark (10-14)', 'Dark-Bright (22-2)'};
        animal_window_means_list = [];
        baselineAnimals = unique(initialBaselineData.Animal);
        for i_an = 1:length(baselineAnimals)
            T_animal = initialBaselineData(initialBaselineData.Animal == baselineAnimals(i_an), :);
            if isempty(T_animal), continue; end
            for i_win = 1:length(time_windows)
                window = time_windows{i_win};
                if window(1) <= window(2), idx_time = T_animal.(ztHourVar) >= window(1) & T_animal.(ztHourVar) <= window(2);
                else, idx_time = T_animal.(ztHourVar) >= window(1) | T_animal.(ztHourVar) <= window(2); end
                mean_act_this_win = mean(T_animal.(current_metric_var)(idx_time), 'omitnan');
                animal_window_means_list = [animal_window_means_list; {baselineAnimals(i_an), window_labels{i_win}, T_animal.Sex(1), mean_act_this_win}];
            end
        end
        plotTable_3D = cell2table(animal_window_means_list, 'VariableNames', {'Animal', 'Window', 'Sex', 'MeanActivity'});
        plotTable_3D.Window = categorical(plotTable_3D.Window, window_labels);
        violinplot(plotTable_3D.MeanActivity, plotTable_3D.Window, 'Parent', ax3D, 'GroupOrder', window_labels, 'ShowData', false, 'ViolinColor', grayColor);
        jitterAmount = 0.15;
        for i_win = 1:length(window_labels)
            window_data = plotTable_3D(plotTable_3D.Window == window_labels{i_win}, :);
            x_scatter = i_win + (rand(height(window_data), 1) - 0.5) * jitterAmount;
            male_idx = window_data.Sex == 'Male';
            scatter(ax3D, x_scatter(male_idx), window_data.MeanActivity(male_idx), 40, maleColor, 'filled', 'MarkerFaceAlpha', 0.9);
            scatter(ax3D, x_scatter(~male_idx), window_data.MeanActivity(~male_idx), 40, femaleColor, 'filled', 'MarkerFaceAlpha', 0.9);
        end
        
        num_tests = 3;
        corrected_alpha = 0.05; % not using Bonferroni
        fprintf('Using standard alpha (0.05)');
        
        y_max = max(plotTable_3D.MeanActivity, [], 'omitnan'); if isempty(y_max) || isnan(y_max), y_max=1; end
        current_ylim = ylim(ax3D);
        y_range = diff(current_ylim); if y_range==0 || isnan(y_range), y_range=y_max; end
        
        bright_data = plotTable_3D(plotTable_3D.Window == 'Bright (0-11)', :);
        dark_data = plotTable_3D(plotTable_3D.Window == 'Dark (12-23)', :);
        stats_table = outerjoin(bright_data(:,{'Animal', 'MeanActivity'}), dark_data(:,{'Animal', 'MeanActivity'}), 'Keys','Animal', 'MergeKeys', true);
        stats_table.Properties.VariableNames(2:3) = {'Bright', 'Dark'};
        [~, p_bvd] = ttest(stats_table.Bright, stats_table.Dark);
        fprintf('Paired T-Test (Bright vs Dark): p = %.4f\n', p_bvd);
        y_pos1 = current_ylim(2) + 0.05*y_range;
        plot_sig_bar(ax3D, [2, 4], p_bvd, y_pos1, corrected_alpha); 
        
        midbright_data = plotTable_3D(plotTable_3D.Window == 'Mid-Bright (3-9)', :);
        middark_data = plotTable_3D(plotTable_3D.Window == 'Mid-Dark (15-21)', :);
        stats_table_mid = outerjoin(midbright_data(:,{'Animal', 'MeanActivity'}), middark_data(:,{'Animal', 'MeanActivity'}), 'Keys','Animal', 'MergeKeys', true);
        stats_table_mid.Properties.VariableNames(2:3) = {'MidBright', 'MidDark'};
        [~, p_mbvmd] = ttest(stats_table_mid.MidBright, stats_table_mid.MidDark);
        fprintf('Paired T-Test (Mid-Bright vs Mid-Dark): p = %.4f\n', p_mbvmd);
        y_pos2 = y_pos1 + 0.15*y_range;
        plot_sig_bar(ax3D, [3, 5], p_mbvmd, y_pos2, corrected_alpha);
        
        transition_data = plotTable_3D(plotTable_3D.Window == 'Bright-Dark (10-14)', :);
        fullday_data = plotTable_3D(plotTable_3D.Window == 'Full Day', :);
        stats_table_trans = outerjoin(transition_data(:,{'Animal', 'MeanActivity'}), fullday_data(:,{'Animal', 'MeanActivity'}), 'Keys','Animal', 'MergeKeys', true);
        stats_table_trans.Properties.VariableNames(2:3) = {'Transition', 'FullDay'};
        [~, p_tvfd] = ttest(stats_table_trans.Transition, stats_table_trans.FullDay);
        fprintf('Paired T-Test (Transition vs Full Day): p = %.4f\n', p_tvfd);
        y_pos3 = y_pos2 + 0.15*y_range;
        plot_sig_bar(ax3D, [1, 6], p_tvfd, y_pos3, corrected_alpha);
        
        ylim(ax3D, [current_ylim(1), y_pos3 + 0.1*y_range]);
        
        xticklabels(ax3D, window_labels); xtickangle(ax3D, 45);
        ylabel(ax3D,['Mean ',current_metric_ylabel]);
        title(ax3D,sprintf('Diurnality (Last 14d Baseline - 300Lux)\n%s', strrep(current_metric_suffix,'_',' ')));
        if strcmp(current_metric_var, normalizedActivityVar), line(ax3D, xlim(ax3D), [0 0], 'Color', 'k', 'LineStyle', ':'); end
        set(ax3D, 'Position', [0.13 0.2 0.775 0.7]); 
        drawnow;
        fig_filename = fullfile(savePath_Fig3, sprintf('Figure3D_Violins_PooledSex_CorrectedStats_%s.png', current_metric_suffix));
        exportgraphics(hFig3D_current, fig_filename); disp(['Saved: ', fig_filename]); close(hFig3D_current);
    catch ME_3D, warning('ErrGen:Fig3D','Error Fig 3D (%s): %s', current_metric_suffix, ME_3D.message); end
    
    %% FIG 3D (Extended): Weekly Breakdown
    disp('--- Starting Figure 3D Weekly Breakdown ---');
    baselineDataOnly = dataTable(dataTable.Condition == targetConditionFig3AB, :);
    max_day = max(baselineDataOnly.DayOfCondition);
    num_weeks = floor((max_day+1)/7);
    for i_week = 1:num_weeks
        start_day = (i_week-1) * 7;
        end_day = i_week * 7 - 1;
        
        weekly_data = baselineDataOnly(baselineDataOnly.DayOfCondition >= start_day & baselineDataOnly.DayOfCondition <= end_day, :);
        
        if isempty(weekly_data)
            fprintf('Skipping Week %d: No data available.\n', i_week);
            continue;
        end
        
        fprintf('\n--- Generating Weekly Figure 3D for Week %d (Days %d-%d) ---\n', i_week, start_day, end_day);
        try
            if ~exist('violinplot.m','file'),error('ViolinPlotNotFound:Fig3D_Weekly','violinplot.m not found.');end
            hFig3D_weekly = figure('Name', sprintf('Fig 3D Weekly W%d (%s)', i_week, current_metric_suffix), 'Visible', 'off', 'Color', 'w');
            ax3D_w = axes('Parent', hFig3D_weekly); hold(ax3D_w, 'on');
            
            time_windows = {[0 23], [0 11], [3 9], [12 23], [15 21], [10 14], [22 2]};
            window_labels = {'Full Day', 'Bright (0-11)', 'Mid-Bright (3-9)', 'Dark (12-23)', 'Mid-Dark (15-21)', 'Bright-Dark (10-14)', 'Dark-Bright (22-2)'};
            
            animal_window_means_list_w = [];
            weekAnimals = unique(weekly_data.Animal);
            for i_an = 1:length(weekAnimals)
                T_animal_w = weekly_data(weekly_data.Animal == weekAnimals(i_an), :);
                if isempty(T_animal_w), continue; end
                for i_win = 1:length(time_windows)
                    window = time_windows{i_win};
                    if window(1) <= window(2), idx_time = T_animal_w.(ztHourVar) >= window(1) & T_animal_w.(ztHourVar) <= window(2);
                    else, idx_time = T_animal_w.(ztHourVar) >= window(1) | T_animal_w.(ztHourVar) <= window(2); end
                    mean_act_this_win = mean(T_animal_w.(current_metric_var)(idx_time), 'omitnan');
                    animal_window_means_list_w = [animal_window_means_list_w; {weekAnimals(i_an), window_labels{i_win}, T_animal_w.Sex(1), mean_act_this_win}];
                end
            end
            if isempty(animal_window_means_list_w), close(hFig3D_weekly); continue; end
            plotTable_3D_w = cell2table(animal_window_means_list_w, 'VariableNames', {'Animal', 'Window', 'Sex', 'MeanActivity'});
            plotTable_3D_w.Window = categorical(plotTable_3D_w.Window, window_labels);
            violinplot(plotTable_3D_w.MeanActivity, plotTable_3D_w.Window, 'Parent', ax3D_w, 'GroupOrder', window_labels, 'ShowData', false, 'ViolinColor', grayColor);
            jitterAmount = 0.15;
            for i_win = 1:length(window_labels)
                window_data = plotTable_3D_w(plotTable_3D_w.Window == window_labels{i_win}, :);
                x_scatter = i_win + (rand(height(window_data), 1) - 0.5) * jitterAmount;
                male_idx = window_data.Sex == 'Male';
                scatter(ax3D_w, x_scatter(male_idx), window_data.MeanActivity(male_idx), 40, maleColor, 'filled', 'MarkerFaceAlpha', 0.9);
                scatter(ax3D_w, x_scatter(~male_idx), window_data.MeanActivity(~male_idx), 40, femaleColor, 'filled', 'MarkerFaceAlpha', 0.9);
            end
            
            num_tests_weekly = 2;
            corrected_alpha_weekly = 0.05; % no bonferroni
            fprintf('--- STATS: Figure 3D WEEK %d for %s ---\n', i_week, current_metric_suffix);
            y_max_w = max(plotTable_3D_w.MeanActivity, [], 'omitnan'); if isempty(y_max_w) || isnan(y_max_w), y_max_w=1; end
            current_ylim_w = ylim(ax3D_w);
            y_range_w = diff(current_ylim_w); if y_range_w==0 || isnan(y_range_w), y_range_w=y_max_w; end
            
            bright_data_w = plotTable_3D_w(plotTable_3D_w.Window == 'Bright (0-11)', :);
            dark_data_w = plotTable_3D_w(plotTable_3D_w.Window == 'Dark (12-23)', :);
            stats_table_w1 = outerjoin(bright_data_w(:,{'Animal', 'MeanActivity'}), dark_data_w(:,{'Animal', 'MeanActivity'}), 'Keys','Animal', 'MergeKeys', true);
            stats_table_w1.Properties.VariableNames(2:3) = {'Bright', 'Dark'};
            [~, p_bvd_w] = ttest(stats_table_w1.Bright, stats_table_w1.Dark);
            fprintf('Paired T-Test (Bright vs Dark): p = %.4f\n', p_bvd_w);
            y_pos1_w = current_ylim_w(2) + 0.05 * y_range_w;
            plot_sig_bar(ax3D_w, [2, 4], p_bvd_w, y_pos1_w, corrected_alpha_weekly);
            
            midbright_data_w = plotTable_3D_w(plotTable_3D_w.Window == 'Mid-Bright (3-9)', :);
            middark_data_w = plotTable_3D_w(plotTable_3D_w.Window == 'Mid-Dark (15-21)', :);
            stats_table_w2 = outerjoin(midbright_data_w(:,{'Animal', 'MeanActivity'}), middark_data_w(:,{'Animal', 'MeanActivity'}), 'Keys','Animal', 'MergeKeys', true);
            stats_table_w2.Properties.VariableNames(2:3) = {'MidBright', 'MidDark'};
            [~, p_mbvmd_w] = ttest(stats_table_w2.MidBright, stats_table_w2.MidDark);
            fprintf('Paired T-Test (Mid-Bright vs Mid-Dark): p = %.4f\n', p_mbvmd_w);
            y_pos2_w = y_pos1_w + 0.15 * y_range_w;
            plot_sig_bar(ax3D_w, [3, 5], p_mbvmd_w, y_pos2_w, corrected_alpha_weekly);
            
            ylim(ax3D_w, [current_ylim_w(1), y_pos2_w + 0.1*y_range_w]);
            
            xticklabels(ax3D_w, window_labels); xtickangle(ax3D_w, 45);
            ylabel(ax3D_w,['Mean ',current_metric_ylabel]);
            title(ax3D_w,sprintf('Baseline Diurnality (300Lux) - Week %d (Days %d-%d)\n%s',i_week, start_day, end_day, strrep(current_metric_suffix,'_',' ')));
            if strcmp(current_metric_var, normalizedActivityVar), line(ax3D_w, xlim(ax3D_w), [0 0], 'Color', 'k', 'LineStyle', ':'); end
            set(ax3D_w, 'Position', [0.13 0.2 0.775 0.7]);
            drawnow;
            fig_filename = fullfile(savePath_Fig3, sprintf('Figure3D_WeeklyBaseline_W%d_CorrectedStats_%s.png', i_week, current_metric_suffix));
            exportgraphics(hFig3D_weekly, fig_filename); disp(['Saved: ', fig_filename]); close(hFig3D_weekly);
        catch ME_3D_weekly, warning('ErrGen:Fig3D_Weekly','Error Fig 3D Weekly (W%d, %s): %s', i_week, current_metric_suffix, ME_3D_weekly.message); end
    end
end
%% --- 4. GENERATE FIGURE 4 COMPONENTS ---
disp('=====================================================');
disp('--- Starting Generation of Figure 4 Components ---');
disp('=====================================================');
for met_idx = 1:length(activity_metrics_to_plot)
    current_metric_var = activity_metrics_to_plot{met_idx};
    current_metric_suffix = metric_suffixes{met_idx};
    current_metric_ylabel = metric_ylabels{met_idx};
    fprintf('\n--- Generating Figure 4 components for metric: %s ---\n', current_metric_var);
    
    %% FIG 4A: 48h Diurnality Line Plots by Condition (with Day/Night avg lines)
    disp('--- Starting Figure 4A: 48h Profiles ---');
    try
        lightingConditions_Fig4A = {'300Lux', '1000Lux', 'FullDark', '300LuxEnd'};
        hFig4A_Subplots = figure('Name', sprintf('Fig 4A: 48h Profiles by Condition (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w', 'Units', 'inches', 'Position', [1 1 7 9]);
        
        % Create a tiled layout
        t = tiledlayout(4, 1, 'TileSpacing', 'compact', 'Padding', 'compact');
        
        global_min_y = Inf; global_max_y = -Inf;
        all_profiles_4A = cell(length(lightingConditions_Fig4A), 1);
        % New storage for day/night averages per condition
        mean_levels_day = cell(length(lightingConditions_Fig4A), 1);
        mean_levels_night = cell(length(lightingConditions_Fig4A), 1);

        % --- Pre-computation Loop ---
        for i_cond = 1:length(lightingConditions_Fig4A)
            data_this_cond = dataTable(ismember(dataTable.Animal, fourCondAnimals) & dataTable.Condition == lightingConditions_Fig4A{i_cond}, :);
            if isempty(data_this_cond), continue; end
            
            % 1. Calculate Hourly Profile
            hourly_mean = groupsummary(data_this_cond, ztHourVar, 'mean', current_metric_var);
            mean_profile_24h = NaN(hoursPerDay, 1);
            [lia, locb] = ismember((0:23)', hourly_mean.(ztHourVar));
            mean_profile_24h(lia) = hourly_mean.(['mean_', current_metric_var])(locb(lia));
            all_profiles_4A{i_cond} = mean_profile_24h;
            
            % 2. Calculate Day (ZT 0-11) and Night (ZT 12-23) Grand Averages for this condition
            day_idx = data_this_cond.(ztHourVar) >= 0 & data_this_cond.(ztHourVar) <= 11;
            night_idx = data_this_cond.(ztHourVar) >= 12 & data_this_cond.(ztHourVar) <= 23;
            mean_levels_day{i_cond} = mean(data_this_cond.(current_metric_var)(day_idx), 'omitnan');
            mean_levels_night{i_cond} = mean(data_this_cond.(current_metric_var)(night_idx), 'omitnan');

            % 3. Update global limits for common Y-axis
            global_min_y = min(global_min_y, min(mean_profile_24h,[],'omitnan'));
            global_max_y = max(global_max_y, max(mean_profile_24h,[],'omitnan'));
        end
        
        % --- Plotting Loop ---
        for i_cond_4A = 1:length(lightingConditions_Fig4A)
            ax4A = nexttile; 
            hold(ax4A, 'on');
            
            % Plot the 48h Profile
            mean_profile_48h = [all_profiles_4A{i_cond_4A}; all_profiles_4A{i_cond_4A}];
            plot(ax4A, (0:47)', mean_profile_48h, 'k-', 'LineWidth', 1.5);
            
            % Plot the Day/Night Average horizontal dashed lines
            if isfinite(mean_levels_day{i_cond_4A})
                 yline(ax4A, mean_levels_day{i_cond_4A}, 'r--', 'LineWidth', 1, 'HandleVisibility', 'off');
            end
            if isfinite(mean_levels_night{i_cond_4A})
                 yline(ax4A, mean_levels_night{i_cond_4A}, 'r--', 'LineWidth', 1, 'HandleVisibility', 'off');
            end

            % Set a globally consistent Y-limit
            if isfinite(global_min_y) && isfinite(global_max_y)
                ylim(ax4A, [global_min_y - 0.1*abs(global_min_y), global_max_y + 0.1*abs(global_max_y)]);
            end
            
            % Add Shading and formatting
            current_ylim = ylim(ax4A);
            patch(ax4A, [ztLightsOffStart, ztLightsOffEnd+1, ztLightsOffEnd+1, ztLightsOffStart], [current_ylim(1) current_ylim(1) current_ylim(2) current_ylim(2)], [0.9 0.9 0.9], 'FaceAlpha', 0.3, 'EdgeColor', 'none');
            patch(ax4A, [ztLightsOffStart+24, ztLightsOffEnd+1+24, ztLightsOffEnd+1+24, ztLightsOffStart+24], [current_ylim(1) current_ylim(1) current_ylim(2) current_ylim(2)], [0.9 0.9 0.9], 'FaceAlpha', 0.3, 'EdgeColor', 'none');
            title(ax4A, lightingConditions_Fig4A{i_cond_4A});
            xticks(ax4A, 0:6:47); xlim(ax4A, [-0.5 47.5]); grid(ax4A, 'on');
        end
        
        % Add shared labels and a main title using the tiled layout handle
        xlabel(t, 'ZT Hour (Double Plotted)');
        ylabel(t, current_metric_ylabel);
        title(t, sprintf('48h Activity Profiles by Condition (%s)', current_metric_suffix));
        
        fig_filename = fullfile(savePath_Fig4, sprintf('Figure4A_SubplotProfiles_SameY_%s.png', current_metric_suffix));
        exportgraphics(hFig4A_Subplots, fig_filename); disp(['Saved: ', fig_filename]); close(hFig4A_Subplots);
    catch ME_4A, warning('ErrGen:Fig4A','Error Fig 4A (%s): %s', current_metric_suffix, ME_4A.message); end
    
    %% FIG 4B: Overlay Activity Profiles with Difference Line (N=12, SEM Shading)
    disp('--- Starting Figure 4B: Overlay Profiles (SEM) ---');
    try
        hFig4B_OverlayDiff = figure('Name', sprintf('Fig 4B: Overlay Profiles & Difference (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w');
        ax4B = axes('Parent', hFig4B_OverlayDiff); hold(ax4B, 'on');
        conditions_Fig4B = {'300Lux', '1000Lux', 'FullDark'};
        lineColors = {[0 0.4470 0.7410], [0.8500 0.3250 0.0980], [0.4660 0.6740 0.1880]};
        mean_profiles = table();
        
        % 1. Calculate and Plot Lines
        for i_cond = 1:length(conditions_Fig4B)
            animals_for_cond = allAnimals; 
            if strcmp(conditions_Fig4B{i_cond},'FullDark'), animals_for_cond=fourCondAnimals; end
            
            data_this_cond = dataTable(ismember(dataTable.Animal, animals_for_cond) & dataTable.Condition == conditions_Fig4B{i_cond}, :);
            hourly_stats = groupsummary(data_this_cond, ztHourVar, {'mean','std','nnz'}, current_metric_var);
            
            mean_prof = NaN(hoursPerDay,1); sem_prof = NaN(hoursPerDay,1);
            [lia,locb] = ismember((0:23)', hourly_stats.(ztHourVar));
            mean_prof(lia) = hourly_stats.(['mean_',current_metric_var])(locb(lia));
            
            % SEM Calculation
            sd_val = hourly_stats.(['std_',current_metric_var])(locb(lia));
            n_val  = hourly_stats.(['nnz_',current_metric_var])(locb(lia));
            sem_val = sd_val ./ sqrt(n_val); sem_val(n_val <= 1) = NaN; 
            sem_prof(lia) = sem_val;
            
            mean_profiles.(conditions_Fig4B{i_cond}) = mean_prof;
            
            plot(ax4B, (0:23)', mean_prof, 'Color', lineColors{i_cond}, 'LineWidth', 2, 'DisplayName', conditions_Fig4B{i_cond});
            fill(ax4B, [(0:23)'; flipud((0:23)')], [mean_prof-sem_prof; flipud(mean_prof+sem_prof)], ...
                lineColors{i_cond}, 'FaceAlpha', 0.2, 'EdgeColor', 'none', 'HandleVisibility', 'off');
        end
        
        % 2. Plot Difference Line
        diff_profile = mean_profiles.('1000Lux') - mean_profiles.('300Lux');
        plot(ax4B, (0:23)', diff_profile, 'k:', 'LineWidth', 1.5, 'DisplayName', 'Diff (1000L-300L)');
        
        % 3. STATS: Hourly T-Test (300 vs 1000) using ANIMAL MEANS (N=12)
        p_values = ones(24, 1);
        stats_animals = allAnimals; 
        
        fprintf('\n--- STATS: Figure 4B (300Lux vs 1000Lux) ---\n');
        fprintf('   Calculating statistics on N=%d animals (Paired T-Test)...\n', length(stats_animals));
        
        prof_300_N12 = NaN(length(stats_animals), 24);
        prof_1000_N12 = NaN(length(stats_animals), 24);
        
        for i = 1:length(stats_animals)
            d3 = dataTable(dataTable.Animal == stats_animals(i) & dataTable.Condition == '300Lux', :);
            if ~isempty(d3), h=groupsummary(d3,ztHourVar,'mean',current_metric_var); [~,l]=ismember(0:23,h.(ztHourVar)); prof_300_N12(i,l>0)=h.(['mean_',current_metric_var])(l(l>0)); end
            d1 = dataTable(dataTable.Animal == stats_animals(i) & dataTable.Condition == '1000Lux', :);
            if ~isempty(d1), h=groupsummary(d1,ztHourVar,'mean',current_metric_var); [~,l]=ismember(0:23,h.(ztHourVar)); prof_1000_N12(i,l>0)=h.(['mean_',current_metric_var])(l(l>0)); end
        end
        
        for zt = 1:24
            v3=prof_300_N12(:,zt); v1=prof_1000_N12(:,zt); idx=~isnan(v3)&~isnan(v1);
            if sum(idx)>=3, [~,p]=ttest(v3(idx),v1(idx)); p_values(zt)=p; end
        end
        
        % 4. PLOT ASTERISKS
        alpha_4B = 0.05; sig_indices = find(p_values < alpha_4B);
        current_ylim = ylim(ax4B); 
        
        if ~isempty(sig_indices)
            sig_zt = sig_indices - 1;
            fprintf('   Significant Differences (p<0.05) found at ZT: %s\n', num2str(sig_zt'));
            y_asterisks = NaN(size(sig_zt));
            for k = 1:length(sig_zt)
                idx = sig_indices(k);
                y_asterisks(k) = max(mean_profiles.('300Lux')(idx), mean_profiles.('1000Lux')(idx)) + (current_ylim(2)-current_ylim(1))*0.05;
            end
            plot(ax4B, sig_zt, y_asterisks, '*', 'Color', 'k', 'MarkerSize', 6, 'HandleVisibility', 'off');
            if max(y_asterisks) > current_ylim(2), ylim(ax4B, [current_ylim(1), max(y_asterisks)*1.05]); end
        else
            fprintf('   No significant hourly differences found.\n');
        end
        
        % 5. SHADING (Last)
        final_ylim = ylim(ax4B);
        patch(ax4B, [ztLightsOffStart, ztLightsOffEnd+1, ztLightsOffEnd+1, ztLightsOffStart], ...
            [final_ylim(1) final_ylim(1) final_ylim(2) final_ylim(2)], ...
            [0.9 0.9 0.9],'FaceAlpha',0.2,'EdgeColor','none', 'HandleVisibility', 'off');
        
        title(ax4B, sprintf('Fig 4B: Activity Profiles & Difference (%s)', current_metric_suffix));
        xlabel(ax4B, 'ZT Hour'); ylabel(ax4B, current_metric_ylabel);
        xticks(ax4B, 0:6:23); xlim(ax4B, [-0.5 23.5]); grid(ax4B, 'on'); legend(ax4B, 'Location', 'northeast');
        set(ax4B, 'Position', [0.13 0.11 0.775 0.815]);
        exportgraphics(hFig4B_OverlayDiff, fullfile(savePath_Fig4, sprintf('Figure4B_OverlayAndDiff_%s.png', current_metric_suffix))); close(hFig4B_OverlayDiff);
    catch ME_4B, warning('ErrGen:Fig4B','Error Fig 4B (%s): %s', current_metric_suffix, ME_4B.message); end
    
%% FIG 4C: Quantification of Diurnality (Violins)
    disp('--- Starting Figure 4C: Diurnality Violins ---');
    try
        if ~exist('violinplot.m','file'),error('ViolinPlotNotFound:Fig4C','violinplot.m not found.');end
        hFig4C_violins = figure('Name', sprintf('Fig 4C: Diurnality Violins by Cond (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w');
        ax4C = axes('Parent', hFig4C_violins); hold(ax4C, 'on');
        conditions_Fig4C = {'300Lux', '1000Lux', 'FullDark', '300LuxEnd'};
        plotData = []; groupData = []; groupOrder = {};
        
        % Data Containers for Cross-Condition Stats
        data_300_ON = []; data_1000_ON = [];
        data_300_OFF = []; data_1000_OFF = [];
        
        for i_cond = 1:length(conditions_Fig4C)
            for lp = ["On", "Off"]
                groupName = [conditions_Fig4C{i_cond}, '-', char(lp)];
                animals_for_cond = allAnimals; if ismember(conditions_Fig4C{i_cond},["FullDark","300LuxEnd"]), animals_for_cond=fourCondAnimals; end
                data_segment = dataTable(ismember(dataTable.Animal, animals_for_cond) & dataTable.Condition == conditions_Fig4C{i_cond} & dataTable.LightPeriod == lp, :);
                if isempty(data_segment), continue; end
                animal_means = groupsummary(data_segment, animalVar, 'mean', current_metric_var);
                
                % Store data for plotting
                plotData = [plotData; animal_means.(['mean_',current_metric_var])];
                groupData = [groupData; repmat(categorical({groupName}), height(animal_means),1)];
                groupOrder{end+1} = groupName;
                
                % Store data for Stats (Specific Comparisons)
                if strcmp(conditions_Fig4C{i_cond}, '300Lux') && lp == "On", data_300_ON = animal_means; end
                if strcmp(conditions_Fig4C{i_cond}, '1000Lux') && lp == "On", data_1000_ON = animal_means; end
                if strcmp(conditions_Fig4C{i_cond}, '300Lux') && lp == "Off", data_300_OFF = animal_means; end
                if strcmp(conditions_Fig4C{i_cond}, '1000Lux') && lp == "Off", data_1000_OFF = animal_means; end
            end
        end
        violinplot(plotData, groupData, 'Parent', ax4C, 'GroupOrder', groupOrder, 'ViolinColor', grayColor, 'ShowData', false);
        
        % Scatter overlay loop (Standard)
        jitterAmount = 0.15;
        for i_group = 1:length(groupOrder)
            group_label = groupOrder{i_group};
            parts = split(group_label, '-'); cond_str = string(parts{1}); lp_str = string(parts{2});
            animals_for_cond = allAnimals; if ismember(cond_str, ["FullDark","300LuxEnd"]), animals_for_cond=fourCondAnimals; end
            group_data = dataTable(ismember(dataTable.Animal, animals_for_cond) & dataTable.Condition == cond_str & dataTable.LightPeriod == lp_str, :);
            animal_means_this_group = groupsummary(group_data, {'Animal', 'Sex'}, 'mean', current_metric_var);
            if ~isempty(animal_means_this_group)
                x_scatter = i_group + (rand(height(animal_means_this_group), 1) - 0.5) * jitterAmount;
                male_idx = animal_means_this_group.Sex == 'Male';
                scatter(ax4C, x_scatter(male_idx), animal_means_this_group.(['mean_',current_metric_var])(male_idx), 40, maleColor, 'filled', 'MarkerFaceAlpha', 0.8, 'HandleVisibility', 'off');
                scatter(ax4C, x_scatter(~male_idx), animal_means_this_group.(['mean_',current_metric_var])(~male_idx), 40, femaleColor, 'filled', 'MarkerFaceAlpha', 0.8, 'HandleVisibility', 'off');
            end
        end

        % --- Cross-Condition (300 vs 1000) ---
        fprintf('\n--- STATS: Fig 4C Cross-Condition (Effect of Light Intensity) ---\n');
        
        % Align 300 and 1000 data (N=12)
        % Assuming animals are in same order or using join (safer)
        t_300_ON = data_300_ON(:, {animalVar, ['mean_', current_metric_var]}); t_300_ON.Properties.VariableNames{2} = 'Val_300';
        t_1000_ON = data_1000_ON(:, {animalVar, ['mean_', current_metric_var]}); t_1000_ON.Properties.VariableNames{2} = 'Val_1000';
        joined_ON = innerjoin(t_300_ON, t_1000_ON);
        
        [~, p_ON_cross] = ttest(joined_ON.Val_300, joined_ON.Val_1000);
        fprintf('Paired T-Test | 300Lux ON vs 1000Lux ON: p = %.4f (N=%d)\n', p_ON_cross, height(joined_ON));

        t_300_OFF = data_300_OFF(:, {animalVar, ['mean_', current_metric_var]}); t_300_OFF.Properties.VariableNames{2} = 'Val_300';
        t_1000_OFF = data_1000_OFF(:, {animalVar, ['mean_', current_metric_var]}); t_1000_OFF.Properties.VariableNames{2} = 'Val_1000';
        joined_OFF = innerjoin(t_300_OFF, t_1000_OFF);
        
        [~, p_OFF_cross] = ttest(joined_OFF.Val_300, joined_OFF.Val_1000);
        fprintf('Paired T-Test | 300Lux OFF vs 1000Lux OFF: p = %.4f (N=%d)\n', p_OFF_cross, height(joined_OFF));

        % --- Original Stats (On vs Off) ---
        disp('--- STATS: Fig 4C Diurnality (On vs Off) ---');
        
        % --- Plotting Significance Bars ---
        y_lims = get(ax4C, 'YLim'); y_range = diff(y_lims);
        y_top = y_lims(2);
        
        % 1. Cross-Condition ON (300 vs 1000)
        if p_ON_cross < 0.05
            x1 = find(strcmp(groupOrder, '300Lux-On'));
            x2 = find(strcmp(groupOrder, '1000Lux-On'));
            y_line = y_top + y_range * 0.15; % High arch
            plot(ax4C, [x1, x1, x2, x2], [y_top+y_range*0.02, y_line, y_line, y_top+y_range*0.02], 'k-', 'LineWidth', 1.2);
             if p_ON_cross < 0.001, str='***'; elseif p_ON_cross<0.01, str='**'; else, str='*'; end
            text(ax4C, mean([x1, x2]), y_line + y_range*0.02, str, 'HorizontalAlignment', 'center', 'FontSize', 10);
            y_top = y_line; % Bump up for next bars
        end

        xtickangle(ax4C, 45);
        ylabel(ax4C, ['Mean ', current_metric_ylabel]);
        title(ax4C, sprintf('Fig 4C: Diurnality Violins (%s)', current_metric_suffix));
        set(ax4C, 'Position', [0.13 0.2 0.775 0.7]); 
        fig_filename = fullfile(savePath_Fig4, sprintf('Figure4C_Violins_GrayPooled_%s.png', current_metric_suffix));
        exportgraphics(hFig4C_violins, fig_filename); disp(['Saved: ', fig_filename]); close(hFig4C_violins);
    catch ME_4C, warning('ErrGen:Fig4C','Error Fig 4C (%s): %s', current_metric_suffix, ME_4C.message); end

    %% FIG 4D: Shift in Activity Peaks During FullDark (CIRCULAR STATS + CoG PEAK)
    disp('--- Starting Figure 4D: Peak Shift (Center of Gravity) ---');
    if strcmp(current_metric_var, normalizedActivityVar)
        disp('Skipping Figure 4D for NormalizedActivity (Percentage change invalid for Z-scores).');
    else
        try
        hFig4D_peak = figure('Name', sprintf('Fig 4D: Peak Shift in FullDark (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w');
        ax4D = axes('Parent', hFig4D_peak);
        
        data_FD = dataTable(ismember(dataTable.Animal, fourCondAnimals) & dataTable.Condition == 'FullDark', :);
        animals_to_test = fourCondAnimals;
        unique_days_FD = unique(data_FD.DayOfCondition);
        
        % 1. Calculate Daily Peaks using Center of Gravity (CoG)
        animal_daily_peaks = NaN(length(unique_days_FD), length(animals_to_test));
        for i_animal = 1:length(animals_to_test)
            animal_id = animals_to_test(i_animal);
            data_this_animal = data_FD(data_FD.Animal == animal_id, :);
            for i_day = 1:length(unique_days_FD)
                day_val = unique_days_FD(i_day);
                data_this_day_animal = data_this_animal(data_this_animal.DayOfCondition == day_val, :);
                if isempty(data_this_day_animal), continue; end
                
                hourly_mean = groupsummary(data_this_day_animal, ztHourVar, 'mean', current_metric_var);
                
                % --- CENTER OF GRAVITY CALCULATION ---
                if ~isempty(hourly_mean)
                    hours = hourly_mean.(ztHourVar);
                    activity = hourly_mean.(['mean_', current_metric_var]);
                    
                    % Zero out negative values if using normalized data (cannot have negative mass)
                    activity(activity < 0) = 0; 
                    
                    if sum(activity) > 0
                        % Convert Time to Angles (0-24h -> 0-2pi)
                        angles = hours * (2*pi / 24);
                        
                        % Calculate Vector Components weighted by Activity
                        X = sum(activity .* cos(angles));
                        Y = sum(activity .* sin(angles));
                        
                        % Calculate Mean Phase (angle)
                        mean_angle = atan2(Y, X);
                        
                        % Convert back to Hours (0-24)
                        cog_zt = mod(mean_angle * (24 / (2*pi)), 24);
                        
                        animal_daily_peaks(i_day, i_animal) = cog_zt;
                    end
                end
            end
        end
        
        disp('--- STATS: Figure 4D (Circular Mean & Shift) ---');
        
        min_day_idx = 1; max_day_idx = length(unique_days_FD);
        if max_day_idx >= 7
            % Define Week Ranges
            first_week_indices = 1:7; 
            last_week_indices = (max_day_idx-6):max_day_idx;
            
            % Get Data (Rows=Days, Cols=Animals)
            first_week_data = animal_daily_peaks(first_week_indices, :);
            last_week_data = animal_daily_peaks(last_week_indices, :);
            
            % --- CIRCULAR MEAN FUNCTION ---
            % Converts ZT (0-24) to Radians, averages vectors, converts back
            calc_circ_mean = @(x) mod(atan2(mean(sin(x*2*pi/24),1,'omitnan'), mean(cos(x*2*pi/24),1,'omitnan')) * 24/(2*pi), 24);
            
            % Calculate Circular Mean for each ANIMAL (Columns)
            first_week_means = calc_circ_mean(first_week_data);
            last_week_means = calc_circ_mean(last_week_data);
            
            % --- CALCULATE SHIFT (Shortest Distance on Circle) ---
            raw_diff = last_week_means - first_week_means;
            circular_diff = mod(raw_diff + 12, 24) - 12; % Wraps result to [-12, +12] range
            
            % Calculate Grand Means for Display
            grand_mean_first = calc_circ_mean(first_week_means'); % Transpose to treat as column
            grand_mean_last = calc_circ_mean(last_week_means');
            
            fprintf('Comparing Mean Peak ZT (CoG) per animal in FullDark:\n');
            fprintf('  - First Week Circular Mean: %.2f\n', grand_mean_first);
            fprintf('  - Last Week Circular Mean:  %.2f\n', grand_mean_last);
            fprintf('  - Mean Shift: %.2f hours\n', mean(circular_diff));
            
            % Test if the SHIFT is significantly different from 0
            [h, p_peak, ci, stats_peak] = ttest(circular_diff, 0);
            fprintf('One-Sample T-Test on Phase Shift (diff vs 0): t=%.3f, p=%.4f\n', stats_peak.tstat, p_peak);
        end
        
        daily_peak_ZT_pooled = NaN(length(unique_days_FD),1);
        for i_day = 1:length(unique_days_FD)
            day_val = unique_days_FD(i_day);
            data_this_day_pooled = data_FD(data_FD.DayOfCondition == day_val, :);
            if ~isempty(data_this_day_pooled)
                hourly_mean = groupsummary(data_this_day_pooled, ztHourVar, 'mean', current_metric_var);
                if ~isempty(hourly_mean)
                    % --- POOLED CENTER OF GRAVITY ---
                    hours = hourly_mean.(ztHourVar);
                    activity = hourly_mean.(['mean_', current_metric_var]);
                    activity(activity < 0) = 0; % Remove negative weights
                    
                    if sum(activity) > 0
                        angles = hours * (2*pi / 24);
                        X = sum(activity .* cos(angles));
                        Y = sum(activity .* sin(angles));
                        mean_angle = atan2(Y, X);
                        daily_peak_ZT_pooled(i_day) = mod(mean_angle * (24 / (2*pi)), 24);
                    end
                end
            end
        end
        plot(ax4D, unique_days_FD(~isnan(daily_peak_ZT_pooled)), daily_peak_ZT_pooled(~isnan(daily_peak_ZT_pooled)), 'o-k', 'LineWidth', 1.5);
        title(ax4D, sprintf('Fig 4D: Peak Activity ZT (CoG) (%s)', current_metric_suffix));
        xlabel(ax4D, 'Day'); ylabel(ax4D, 'Peak ZT'); ylim(ax4D, [-1 24]);
        fig_filename = fullfile(savePath_Fig4, sprintf('Figure4D_PeakShift_CoG_%s.png', current_metric_suffix));
        exportgraphics(hFig4D_peak, fig_filename); disp(['Saved: ', fig_filename]); close(hFig4D_peak);
    catch ME_4D, warning('ErrGen:Fig4D','Error Fig 4D (%s): %s', current_metric_suffix, ME_4D.message); end
    end
%% FIG 4E: Percentage Change (ON vs OFF) Violin Plot
    disp('--- Starting Figure 4E: ON/OFF Pct Change Violins ---');
    if strcmp(current_metric_var, normalizedActivityVar)
        disp('Skipping Figure 4E for NormalizedActivity (Percentage change invalid for Z-scores).');
    else
        try
        % 1. Setup Figure
        if ~exist('violinplot.m','file'),error('ViolinPlotNotFound:Fig4E','violinplot.m not found.');end
        hFig4E_violins = figure('Name', sprintf('Fig 4E: ON/OFF Pct Change Violins (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w');
        ax4E = axes('Parent', hFig4E_violins); hold(ax4E, 'on');
        
        % 2. Define Conditions and Data Containers
        conditions_Fig4E = {'300Lux', '1000Lux', 'FullDark', '300LuxEnd'};
        plotData_4E = []; groupData_4E = []; groupOrder_4E = {};
        
        % Temporary storage to avoid re-calculating for the scatter loop
        grouped_pct_data = containers.Map; 
        
        % 3. Main Data Processing Loop
        for i_cond = 1:length(conditions_Fig4E)
            cond_name = conditions_Fig4E{i_cond};
            
            % Apply specific animal selection logic
            animals_for_cond = allAnimals; 
            if ismember(cond_name, ["FullDark","300LuxEnd"]), animals_for_cond=fourCondAnimals; end
            
            % Get raw data for this condition
            data_segment = dataTable(ismember(dataTable.Animal, animals_for_cond) & dataTable.Condition == cond_name, :);
            if isempty(data_segment), continue; end
            
            % Aggregate mean by Animal, Sex, and LightPeriod
            animal_means = groupsummary(data_segment, {'Animal', 'Sex', 'LightPeriod'}, 'mean', current_metric_var);
            
            % Pivot data so 'On' and 'Off' are columns for each animal
            try
                pivoted_data = unstack(animal_means, ['mean_' current_metric_var], 'LightPeriod', 'GroupingVariables', {'Animal', 'Sex'});
                
                % Check if both On and Off exist after unstacking
                if all(ismember({'On', 'Off'}, pivoted_data.Properties.VariableNames))
                    % Calculate INDIVIDUAL Percentage Change: ((ON - OFF) / OFF) * 100
                    % Add epsilon to denominator to prevent division by zero
                    epsilon = 1e-9; 
                    pivoted_data.PctChange = ((pivoted_data.On - pivoted_data.Off) ./ (pivoted_data.Off + epsilon)) * 100;
                    
                    % Remove NaNs or Infs if OFF was 0
                    valid_rows = isfinite(pivoted_data.PctChange);
                    pivoted_data = pivoted_data(valid_rows, :);
    
                    % Append to main plotting data
                    plotData_4E = [plotData_4E; pivoted_data.PctChange];
                    groupData_4E = [groupData_4E; repmat(categorical({cond_name}), height(pivoted_data), 1)];
                    groupOrder_4E{end+1} = cond_name;
                    
                    % Store for the scatter overlay loop
                    grouped_pct_data(cond_name) = pivoted_data;
                else
                    warning('Fig4E:IncompleteData', 'Missing ON or OFF data for condition: %s', cond_name);
                end
            catch ME_Pivot
                 warning('Fig4E:PivotError', 'Could not pivot data for %s: %s', cond_name, ME_Pivot.message);
            end
        end
        
        % 4. Create Base Gray Violin Plot
        violinplot(plotData_4E, groupData_4E, 'Parent', ax4E, 'GroupOrder', groupOrder_4E, 'ViolinColor', grayColor, 'ShowData', false);
        
        % 5. Reference line at 0%
        yline(ax4E, 0, 'k--', 'LineWidth', 1.2, 'HandleVisibility', 'off');
        
        % 6. Overlay Individual Data Points (Males vs Females)
        jitterAmount = 0.15;
        for i_group = 1:length(groupOrder_4E)
            group_label = groupOrder_4E{i_group};
            
            if isKey(grouped_pct_data, group_label)
                group_data = grouped_pct_data(group_label);
                
                % Jitter x-positions
                x_scatter = i_group + (rand(height(group_data), 1) - 0.5) * jitterAmount;
                
                % Split by Sex and Plot
                male_idx = group_data.Sex == 'Male';
                scatter(ax4E, x_scatter(male_idx), group_data.PctChange(male_idx), 40, maleColor, 'filled', 'MarkerFaceAlpha', 0.8, 'HandleVisibility', 'off');
                scatter(ax4E, x_scatter(~male_idx), group_data.PctChange(~male_idx), 40, femaleColor, 'filled', 'MarkerFaceAlpha', 0.8, 'HandleVisibility', 'off');
            end
        end
    
        try
            disp('--- STATS: Figure 4E (Corrected: Paired t-tests) ---');
            hold(ax4E, 'on');

            % --- 1. Create a single wide table to align all animals ---
            % Get data + Animal IDs from the map.
            temp_table_300 = grouped_pct_data('300Lux');
            table_300 = temp_table_300(:, {'Animal', 'PctChange'});
            table_300.Properties.VariableNames{'PctChange'} = 'Pct_300';
            
            temp_table_1000 = grouped_pct_data('1000Lux');
            table_1000 = temp_table_1000(:, {'Animal', 'PctChange'});
            table_1000.Properties.VariableNames{'PctChange'} = 'Pct_1000';
            
            temp_table_dark = grouped_pct_data('FullDark');
            table_dark = temp_table_dark(:, {'Animal', 'PctChange'});
            table_dark.Properties.VariableNames{'PctChange'} = 'Pct_Dark';
            
            temp_table_end = grouped_pct_data('300LuxEnd');
            table_end = temp_table_end(:, {'Animal', 'PctChange'});
            table_end.Properties.VariableNames{'PctChange'} = 'Pct_End';

            % Join them all.
            paired_table = outerjoin(table_300, table_1000, 'Keys', 'Animal', 'MergeKeys', true);
            paired_table = outerjoin(paired_table, table_dark, 'Keys', 'Animal', 'MergeKeys', true);
            paired_table = outerjoin(paired_table, table_end, 'Keys', 'Animal', 'MergeKeys', true);
            
            p_values = containers.Map('KeyType', 'char', 'ValueType', 'double');
            
            % --- 2. Run PAIRED t-tests (ttest) ---
            [~, p] = ttest(paired_table.Pct_300, paired_table.Pct_1000); 
            p_values('1000Lux_vs_300Lux') = p;
            fprintf('Paired T-Test (300Lux vs 1000Lux) [N=%d]: p = %.4f\n', sum(~isnan(paired_table.Pct_300) & ~isnan(paired_table.Pct_1000)), p);

            [~, p] = ttest(paired_table.Pct_300, paired_table.Pct_Dark); 
            p_values('FullDark_vs_300Lux') = p;
            fprintf('Paired T-Test (300Lux vs FullDark) [N=%d]: p = %.4f\n', sum(~isnan(paired_table.Pct_300) & ~isnan(paired_table.Pct_Dark)), p);

            [~, p] = ttest(paired_table.Pct_300, paired_table.Pct_End); 
            p_values('300LuxEnd_vs_300Lux') = p;
            fprintf('Paired T-Test (300Lux vs 300LuxEnd) [N=%d]: p = %.4f\n', sum(~isnan(paired_table.Pct_300) & ~isnan(paired_table.Pct_End)), p);

            [~, p] = ttest(paired_table.Pct_Dark, paired_table.Pct_1000);
            p_values('FullDark_vs_1000Lux') = p;
            fprintf('Paired T-Test (FullDark vs 1000Lux) [N=%d]: p = %.4f\n', sum(~isnan(paired_table.Pct_Dark) & ~isnan(paired_table.Pct_1000)), p);

            [~, p] = ttest(paired_table.Pct_End, paired_table.Pct_1000);
            p_values('300LuxEnd_vs_1000Lux') = p;
            fprintf('Paired T-Test (300LuxEnd vs 1000Lux) [N=%d]: p = %.4f\n', sum(~isnan(paired_table.Pct_End) & ~isnan(paired_table.Pct_1000)), p);

            [~, p] = ttest(paired_table.Pct_Dark, paired_table.Pct_End);
            p_values('FullDark_vs_300LuxEnd') = p;
            fprintf('Paired T-Test (FullDark vs 300LuxEnd) [N=%d]: p = %.4f\n', sum(~isnan(paired_table.Pct_Dark) & ~isnan(paired_table.Pct_End)), p);

            % --- 3. Plot Significance Bars (4C Style) ---
            x_control = find(strcmp(groupOrder_4E, '300Lux'));
            x_1000 = find(strcmp(groupOrder_4E, '1000Lux'));
            x_dark = find(strcmp(groupOrder_4E, 'FullDark'));
            x_end = find(strcmp(groupOrder_4E, '300LuxEnd'));
            
            y_lims_4E = get(ax4E, 'YLim');
            y_range_4E = diff(y_lims_4E);
            
            % Define 6 Y-levels
            y_level_1 = y_lims_4E(2) + y_range_4E * 0.05; y_text_1 = y_level_1 + y_range_4E * 0.01;
            y_level_2 = y_text_1 + y_range_4E * 0.05; y_text_2 = y_level_2 + y_range_4E * 0.01;
            y_level_3 = y_text_2 + y_range_4E * 0.05; y_text_3 = y_level_3 + y_range_4E * 0.01;
            y_level_4 = y_text_3 + y_range_4E * 0.05; y_text_4 = y_level_4 + y_range_4E * 0.01;
            y_level_5 = y_text_4 + y_range_4E * 0.05; y_text_5 = y_level_5 + y_range_4E * 0.01;
            y_level_6 = y_text_5 + y_range_4E * 0.05; y_text_6 = y_level_6 + y_range_4E * 0.01;
            
            new_y_max_4E = y_text_6 + y_range_4E * 0.03; 
            set(ax4E, 'YLim', [y_lims_4E(1), new_y_max_4E]);
                        
            % 1. 1000Lux vs 300Lux
            p = p_values('1000Lux_vs_300Lux');
            % Standard alpha 0.05 for this specific hypothesis
            if p < 0.05 
                if p < 0.001, str_sig = '***'; elseif p < 0.01, str_sig = '**'; else, str_sig = '*'; end
                plot(ax4E, [x_control, x_1000], [y_level_1, y_level_1], 'k-', 'LineWidth', 1.2);
                text(ax4E, mean([x_control, x_1000]), y_text_1, str_sig, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 10, 'FontWeight', 'bold');
            end
            
            % 2. FullDark vs 300Lux
            p = p_values('FullDark_vs_300Lux');
            if p < 0.001, str_sig = '***'; elseif p < 0.01, str_sig = '**'; elseif p < 0.05, str_sig = '*'; else, str_sig = 'n.s.'; end
            if ~strcmp(str_sig, 'n.s.')
                plot(ax4E, [x_control, x_dark], [y_level_2, y_level_2], 'k-', 'LineWidth', 1.2);
                text(ax4E, mean([x_control, x_dark]), y_text_2, str_sig, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 10, 'FontWeight', 'bold');
            end
            
            % 3. 300LuxEnd vs 300Lux
            p = p_values('300LuxEnd_vs_300Lux');
            if p < 0.001, str_sig = '***'; elseif p < 0.01, str_sig = '**'; elseif p < 0.05, str_sig = '*'; else, str_sig = 'n.s.'; end
            if ~strcmp(str_sig, 'n.s.')
                plot(ax4E, [x_control, x_end], [y_level_3, y_level_3], 'k-', 'LineWidth', 1.2);
                text(ax4E, mean([x_control, x_end]), y_text_3, str_sig, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 10, 'FontWeight', 'bold');
            end
            
            % 4. FullDark vs 1000Lux
            p = p_values('FullDark_vs_1000Lux');
            if p < 0.001, str_sig = '***'; elseif p < 0.01, str_sig = '**'; elseif p < 0.05, str_sig = '*'; else, str_sig = 'n.s.'; end
            if ~strcmp(str_sig, 'n.s.')
                plot(ax4E, [x_1000, x_dark], [y_level_4, y_level_4], 'k-', 'LineWidth', 1.2);
                text(ax4E, mean([x_1000, x_dark]), y_text_4, str_sig, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 10, 'FontWeight', 'bold');
            end

            % 5. 300LuxEnd vs 1000Lux
            p = p_values('300LuxEnd_vs_1000Lux');
            if p < 0.001, str_sig = '***'; elseif p < 0.01, str_sig = '**'; elseif p < 0.05, str_sig = '*'; else, str_sig = 'n.s.'; end
            if ~strcmp(str_sig, 'n.s.')
                plot(ax4E, [x_1000, x_end], [y_level_5, y_level_5], 'k-', 'LineWidth', 1.2);
                text(ax4E, mean([x_1000, x_end]), y_text_5, str_sig, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 10, 'FontWeight', 'bold');
            end

            % 6. FullDark vs 300LuxEnd
            p = p_values('FullDark_vs_300LuxEnd');
            if p < 0.001, str_sig = '***'; elseif p < 0.01, str_sig = '**'; elseif p < 0.05, str_sig = '*'; else, str_sig = 'n.s.'; end
            if ~strcmp(str_sig, 'n.s.')
                plot(ax4E, [x_dark, x_end], [y_level_6, y_level_6], 'k-', 'LineWidth', 1.2);
                text(ax4E, mean([x_dark, x_end]), y_text_6, str_sig, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 10, 'FontWeight', 'bold');
            end
            
        catch ME_Stats_4E
            warning('ErrGen:Fig4EStats', 'Could not plot significance for Fig 4E: %s', ME_Stats_4E.message);
        end

        % 7. Formatting and Saving
        ylabel(ax4E, {'% Change (ON vs OFF)', '((ON-OFF)/OFF \times 100)'});
        title(ax4E, sprintf('Fig 4E: %% Activity Change (ON vs OFF) (%s)', current_metric_suffix));
        set(ax4E, 'Position', [0.13 0.11 0.775 0.815]); 
    
        fig_filename = fullfile(savePath_Fig4, sprintf('Figure4E_PctChange_Violins_%s.png', current_metric_suffix));
        exportgraphics(hFig4E_violins, fig_filename); 
        disp(['Saved: ', fig_filename]); 
        close(hFig4E_violins);
        
    catch ME_4E
        warning('ErrGen:Fig4E','Error Fig 4E (%s): %s', current_metric_suffix, ME_4E.message); end
    end
    %% FIG 4E SUPPLEMENT: Descriptive Stats for Diurnality Index (Req 3)
    if strcmp(current_metric_var, normalizedActivityVar)
        disp('Skipping Figure 4E Supplement for NormalizedActivity (Percentage change invalid for Z-scores).');
    else
        try
        disp('--- STATS: Descriptive Means for Diurnality Index (Pct Change) ---');
        % 300 Lux
        % Note: This logic must be self-contained and recalculate PctChange
        data_300 = dataTable(dataTable.Condition == "300Lux" & ismember(dataTable.Animal, allAnimals), :);
        stats_300 = groupsummary(data_300, {'Animal','LightPeriod'}, 'mean', current_metric_var);
        piv_300 = unstack(stats_300, ['mean_', current_metric_var], 'LightPeriod', 'GroupingVariables', {'Animal'});
        % Ensure both On and Off columns exist
        if all(ismember({'On', 'Off'}, piv_300.Properties.VariableNames))
            pct_300 = ((piv_300.On - piv_300.Off) ./ piv_300.Off) * 100;
            pct_300_finite = pct_300(isfinite(pct_300));
            fprintf('300 Lux Diurnality Index: Mean = %.2f %%, SEM = %.2f %% (N=%d)\n', ...
                mean(pct_300_finite), std(pct_300_finite)/sqrt(length(pct_300_finite)), length(pct_300_finite));
        else
            disp('Could not calculate 300Lux Diurnality Index, On/Off data missing.');
        end
        
        % Full Dark
        data_FD = dataTable(dataTable.Condition == "FullDark" & ismember(dataTable.Animal, fourCondAnimals), :);
        stats_FD = groupsummary(data_FD, {'Animal','LightPeriod'}, 'mean', current_metric_var);
        piv_FD = unstack(stats_FD, ['mean_', current_metric_var], 'LightPeriod', 'GroupingVariables', {'Animal'});
        % Ensure both On and Off columns exist
        if all(ismember({'On', 'Off'}, piv_FD.Properties.VariableNames))
            pct_FD = ((piv_FD.On - piv_FD.Off) ./ piv_FD.Off) * 100;
            pct_FD_finite = pct_FD(isfinite(pct_FD));
            fprintf('FullDark Diurnality Index: Mean = %.2f %%, SEM = %.2f %% (N=%d)\n', ...
                mean(pct_FD_finite), std(pct_FD_finite)/sqrt(length(pct_FD_finite)), length(pct_FD_finite));
        else
            disp('Could not calculate FullDark Diurnality Index, On/Off data missing.');
        end
    catch ME_Desc
        warning('ErrGen:DescStats', 'Error calculating descriptive stats: %s', ME_Desc.message); 
    end
    end
    %% FIG 4F: FullDark Anticipatory Peak Analysis (Dusk & Dawn) - REVISED
    % Purpose: Test for anticipatory peaks at both Dusk (10-12) and Dawn (21-23) in FullDark.
    disp('--- Starting Figure 4F: FullDark Anticipatory Analysis (Dusk & Dawn) ---');
    try
        hFig4F_current = figure('Name', sprintf('Fig 4F: FullDark Anticipation (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w');
        ax4F = axes('Parent', hFig4F_current); hold(ax4F, 'on');
        
        % 1. Define Windows (Day, Dusk, Night, Dawn)
        window_labels_4F = {'Day Avg (0-11)', 'Dusk Anti (10-12)', 'Night Avg (12-23)', 'Dawn Anti (21-23)'};
        
        % 2. Extract Data (N=8 FullDark Animals)
        data_FD_4F = dataTable(ismember(dataTable.Animal, fourCondAnimals) & dataTable.Condition == 'FullDark', :);
        animal_means_4F = [];
        
        unique_animals_FD = unique(data_FD_4F.Animal);
        for i_an = 1:length(unique_animals_FD)
            current_an = unique_animals_FD(i_an);
            an_data = data_FD_4F(data_FD_4F.Animal == current_an, :);
            
            if ismember(current_an, maleAnimals), sex_val = 'Male'; else, sex_val = 'Female'; end
            
            % 1. Day Avg (ZT 0-11)
            idx_day = an_data.(ztHourVar) >= 0 & an_data.(ztHourVar) <= 11;
            mean_day = mean(an_data.(current_metric_var)(idx_day), 'omitnan');
            
            % 2. Dusk Anticipation (ZT 10-12)
            idx_dusk = an_data.(ztHourVar) >= 10 & an_data.(ztHourVar) <= 12;
            mean_dusk = mean(an_data.(current_metric_var)(idx_dusk), 'omitnan');
            
            % 3. Night Avg (ZT 12-23)
            idx_night = an_data.(ztHourVar) >= 12 & an_data.(ztHourVar) <= 23;
            mean_night = mean(an_data.(current_metric_var)(idx_night), 'omitnan');
            
            % 4. Dawn Anticipation (ZT 21-23)
            idx_dawn = an_data.(ztHourVar) >= 21 & an_data.(ztHourVar) <= 23;
            mean_dawn = mean(an_data.(current_metric_var)(idx_dawn), 'omitnan');
            
            % Store
            animal_means_4F = [animal_means_4F; {current_an, window_labels_4F{1}, sex_val, mean_day}];
            animal_means_4F = [animal_means_4F; {current_an, window_labels_4F{2}, sex_val, mean_dusk}];
            animal_means_4F = [animal_means_4F; {current_an, window_labels_4F{3}, sex_val, mean_night}];
            animal_means_4F = [animal_means_4F; {current_an, window_labels_4F{4}, sex_val, mean_dawn}];
        end
        
        % 3. Plotting
        plotTable_4F = cell2table(animal_means_4F, 'VariableNames', {'Animal', 'Window', 'Sex', 'MeanActivity'});
        plotTable_4F.Window = categorical(plotTable_4F.Window, window_labels_4F);
        
        violinplot(plotTable_4F.MeanActivity, plotTable_4F.Window, 'Parent', ax4F, 'GroupOrder', window_labels_4F, 'ShowData', false, 'ViolinColor', grayColor);
        
        % Scatter Overlay
        jitterAmount = 0.15;
        for i_win = 1:length(window_labels_4F)
            win_data = plotTable_4F(plotTable_4F.Window == window_labels_4F{i_win}, :);
            x_scatter = i_win + (rand(height(win_data), 1) - 0.5) * jitterAmount;
            
            male_idx = strcmp(win_data.Sex, 'Male');
            scatter(ax4F, x_scatter(male_idx), win_data.MeanActivity(male_idx), 40, maleColor, 'filled', 'MarkerFaceAlpha', 0.9);
            scatter(ax4F, x_scatter(~male_idx), win_data.MeanActivity(~male_idx), 40, femaleColor, 'filled', 'MarkerFaceAlpha', 0.9);
        end
        
        % 4. Stats (Paired T-Tests)
        disp('--- STATS: Figure 4F (Anticipatory Peaks) ---');
        piv_4F = unstack(plotTable_4F, 'MeanActivity', 'Window', 'GroupingVariables', 'Animal');
        
        % Dynamically find column names to be safe
        vars = piv_4F.Properties.VariableNames;
        col_day   = vars{contains(vars, 'Day')};
        col_dusk  = vars{contains(vars, 'Dusk')};
        col_night = vars{contains(vars, 'Night')};
        col_dawn  = vars{contains(vars, 'Dawn')};
        
        % Test 1: Day Avg vs Dusk Anti
        [~, p_dusk] = ttest(piv_4F.(col_day), piv_4F.(col_dusk));
        fprintf('Paired T-Test | Day Avg vs Dusk Anti:   p = %.4f\n', p_dusk);
        
        % Test 2: Night Avg vs Dawn Anti
        [~, p_dawn] = ttest(piv_4F.(col_night), piv_4F.(col_dawn));
        fprintf('Paired T-Test | Night Avg vs Dawn Anti: p = %.4f\n', p_dawn);
        
        % Plot Significance Bars (Alpha 0.05)
        y_max = max(plotTable_4F.MeanActivity);
        y_range = y_max - min(plotTable_4F.MeanActivity);
        y_pos = y_max + 0.05 * y_range;
        
        if p_dusk < 0.05
            plot_sig_bar(ax4F, [1, 2], p_dusk, y_pos, 0.05);
            y_pos = y_pos + 0.1 * y_range; % Stack bars if needed
        end
        
        if p_dawn < 0.05
            plot_sig_bar(ax4F, [3, 4], p_dawn, y_pos, 0.05);
        end
        
        ylabel(ax4F, ['Mean ', current_metric_ylabel]);
        title(ax4F, sprintf('FullDark Anticipation (%s)', strrep(current_metric_suffix,'_',' ')));
        set(ax4F, 'Position', [0.13 0.2 0.775 0.7]); 
        
        fig_filename = fullfile(savePath_Fig4, sprintf('Figure4F_FullDarkAnticipation_%s.png', current_metric_suffix));
        exportgraphics(hFig4F_current, fig_filename); disp(['Saved: ', fig_filename]); close(hFig4F_current);
        
    catch ME_4F
        warning('ErrGen:Fig4F','Error Fig 4F (%s): %s', current_metric_suffix, ME_4F.message); 
    end
end

%% --- 5. GENERATE FIGURE 5 COMPONENTS ---
disp('=====================================================');
disp('--- Starting Generation of Figure 5 Components ---');
disp('=====================================================');
for met_idx = 1:length(activity_metrics_to_plot)
    current_metric_var = activity_metrics_to_plot{met_idx};
    current_metric_suffix = metric_suffixes{met_idx};
    current_metric_ylabel = metric_ylabels{met_idx};
    fprintf('\n--- Generating Figure 5 components for metric: %s ---\n', current_metric_var);
    sex_groups = {maleAnimals, femaleAnimals};
    sex_labels = {'Male', 'Female'};
    
    %% FIG 5A: Sex differences in 24h profiles (INDIVIDUAL TRACES + SEM SHADING + LINKED AXES)
    disp('--- Starting Figure 5A: Sex Profiles with Individual Traces ---');
    conditionsFor5A = {'300Lux', '1000Lux', 'FullDark'};
    condColors_5A = {[0 0.4470 0.7410], [0.8500 0.3250 0.0980], [0.4660 0.6740 0.1880]};
    
    % --- PRE-CALCULATE GLOBAL Y-LIMITS ---
    global_y_max = -inf; global_y_min = inf;
    for i_s = 1:2
        for i_c = 1:3
            a_list = sex_groups{i_s}; 
            if strcmp(conditionsFor5A{i_c},'FullDark'), a_list = intersect(a_list, fourCondAnimals); end
            d_sub = dataTable(ismember(dataTable.Animal, a_list) & dataTable.Condition == conditionsFor5A{i_c}, :);
            if isempty(d_sub), continue; end
            h_st = groupsummary(d_sub, ztHourVar, {'mean','std','nnz'}, current_metric_var);
            mn = h_st.(['mean_',current_metric_var]);
            se = h_st.(['std_',current_metric_var]) ./ sqrt(h_st.(['nnz_',current_metric_var]));
            global_y_max = max(global_y_max, max(mn + se));
            global_y_min = min(global_y_min, min(mn - se));
        end
    end
    y_range_glob = global_y_max - global_y_min;
    global_ylim = [min(0, global_y_min - 0.05*y_range_glob), global_y_max + 0.15*y_range_glob];
    
    for sex_idx = 1:length(sex_groups)
        hFig5A_current = figure('Name', sprintf('Figure 5A: %s Profiles (%s)', sex_labels{sex_idx}, current_metric_suffix), 'Visible', 'off', 'Color', 'w');
        ax5A = axes('Parent', hFig5A_current); hold(ax5A, 'on');
        prof_300 = []; prof_1000 = []; 
        
        for i_cond = 1:length(conditionsFor5A)
            animals_for_cond = sex_groups{sex_idx};
            if strcmp(conditionsFor5A{i_cond}, 'FullDark'), animals_for_cond = intersect(animals_for_cond, fourCondAnimals); end
            sexCondData = dataTable(ismember(dataTable.Animal, animals_for_cond) & dataTable.Condition == conditionsFor5A{i_cond}, :);
            if isempty(sexCondData), continue; end
            
            % 1. Individual Traces
            unique_ans = unique(sexCondData.Animal);
            for ia = 1:length(unique_ans)
                an_d = sexCondData(sexCondData.Animal == unique_ans(ia), :);
                an_hourly = groupsummary(an_d, ztHourVar, 'mean', current_metric_var); an_hourly = sortrows(an_hourly, ztHourVar);
                plot(ax5A, an_hourly.(ztHourVar), an_hourly.(['mean_', current_metric_var]), 'Color', [condColors_5A{i_cond}, 0.15], 'LineWidth', 0.5, 'HandleVisibility', 'off');
            end
            % 2. Group Stats (SEM)
            hourlyStats = groupsummary(sexCondData, ztHourVar, {'mean','std','nnz'}, current_metric_var);
            meanProfile = NaN(hoursPerDay,1); semProfile = NaN(hoursPerDay,1);
            [lia,locb] = ismember((0:23)',hourlyStats.(ztHourVar));
            meanProfile(lia) = hourlyStats.(['mean_',current_metric_var])(locb(lia));
            semProfile(lia) = hourlyStats.(['std_',current_metric_var])(locb(lia)) ./ sqrt(hourlyStats.(['nnz_',current_metric_var])(locb(lia)));
            if strcmp(conditionsFor5A{i_cond}, '300Lux'), prof_300 = meanProfile; end
            if strcmp(conditionsFor5A{i_cond}, '1000Lux'), prof_1000 = meanProfile; end
            
            plot(ax5A, (0:23)', meanProfile, 'LineWidth', 2, 'DisplayName', conditionsFor5A{i_cond}, 'Color', condColors_5A{i_cond});
            fill(ax5A, [(0:23)'; flipud((0:23)')], [meanProfile-semProfile; flipud(meanProfile+semProfile)], condColors_5A{i_cond}, 'FaceAlpha', 0.2, 'EdgeColor', 'none', 'HandleVisibility', 'off');
        end
        
        if strcmp(current_metric_var, 'NormalizedActivity'), ylim(ax5A, [-1.5, 3]); else, ylim(ax5A, global_ylim); end
        current_ylim = ylim(ax5A);
        
        % --- STATS (N=6 Paired) ---
        fprintf('--- STATS: Fig 5A (%s) 300Lux vs 1000Lux (N=%d) ---\n', sex_labels{sex_idx}, length(sex_groups{sex_idx}));
        p_values = ones(24, 1);
        stats_animals = sex_groups{sex_idx}; 
        prof_300_sex = NaN(length(stats_animals), 24); prof_1000_sex = NaN(length(stats_animals), 24);
        
        for i = 1:length(stats_animals)
            d3=dataTable(dataTable.Animal==stats_animals(i)&dataTable.Condition=='300Lux',:); 
            if ~isempty(d3), h=groupsummary(d3,ztHourVar,'mean',current_metric_var); [~,l]=ismember(0:23,h.(ztHourVar)); prof_300_sex(i,l>0)=h.(['mean_',current_metric_var])(l(l>0)); end
            d1=dataTable(dataTable.Animal==stats_animals(i)&dataTable.Condition=='1000Lux',:); 
            if ~isempty(d1), h=groupsummary(d1,ztHourVar,'mean',current_metric_var); [~,l]=ismember(0:23,h.(ztHourVar)); prof_1000_sex(i,l>0)=h.(['mean_',current_metric_var])(l(l>0)); end
        end
        for zt = 1:24
            v3=prof_300_sex(:,zt); v1=prof_1000_sex(:,zt); idx=~isnan(v3)&~isnan(v1);
            if sum(idx)>=3, [~,p]=ttest(v3(idx),v1(idx)); p_values(zt)=p; end
        end
        
        alpha_5A = 0.05; sig_idx = find(p_values < alpha_5A);
        if ~isempty(sig_idx)
            fprintf('    Significant Hours (p < %.2f): ZT %s\n', alpha_5A, num2str(sig_idx' - 1));
            zt_sig = sig_idx - 1;
            y_pts = max([prof_300(sig_idx), prof_1000(sig_idx)], [], 2);
            y_buff = (current_ylim(2)-current_ylim(1)) * 0.05;
            plot(ax5A, zt_sig, y_pts + y_buff, '*', 'Color', 'k', 'MarkerSize', 5, 'HandleVisibility', 'off');
            if max(y_pts+y_buff) > current_ylim(2), ylim(ax5A, [current_ylim(1), max(y_pts+y_buff)+y_buff]); end
        else
            fprintf('    No significant hours found.\n');
        end
        
        final_ylim = ylim(ax5A);
        patch(ax5A, [ztLightsOffStart, ztLightsOffEnd+1, ztLightsOffEnd+1, ztLightsOffStart], ...
            [final_ylim(1) final_ylim(1) final_ylim(2) final_ylim(2)], ...
            [0.9 0.9 0.9],'FaceAlpha',0.2,'EdgeColor','none', 'HandleVisibility', 'off');
            
        title(ax5A,sprintf('Fig 5A: %s Profiles (%s) [Shading=SEM]', sex_labels{sex_idx}, current_metric_suffix));
        xlabel(ax5A,'ZT Hour'); ylabel(ax5A,['Mean ', current_metric_ylabel]);
        xticks(ax5A,0:6:23); xlim(ax5A,[-0.5,23.5]); grid(ax5A,'on'); legend(ax5A, 'Location','northwest');
        set(ax5A, 'Position', [0.13 0.11 0.775 0.815]);
        exportgraphics(hFig5A_current, fullfile(savePath_Fig5, sprintf('Figure5A_Profile_%s_%s_SEM.png', sex_labels{sex_idx}, current_metric_suffix))); close(hFig5A_current);
    end

%% FIG 5B: Diurnality Violins by Sex and Condition (12-Violin Plot) WITH DESCRIPTIVE STATS
    disp('--- Starting Figure 5B: 12-Violin Plot ---');
    try
        if ~exist('violinplot.m','file'),error('ViolinPlotNotFound:Fig5B','violinplot.m not found.');end
        hFig5B_current = figure('Name', sprintf('Fig 5B: Diurnality by Sex and Cond (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w', 'Units', 'inches', 'Position', [1 1 12 7]);
        ax5B = axes('Parent', hFig5B_current); hold(ax5B, 'on');
        
        plotTable_5B = table();
        conditions_5B = {'300Lux', '1000Lux', 'FullDark'};
        lightPeriods_5B = {'On', 'Off'};
        all_animal_means_5B = table(unique(dataTable.Animal), 'VariableNames', {'Animal'});
        
        % Data Collection Loop
        for i_cond = 1:length(conditions_5B)
            for i_lp = 1:length(lightPeriods_5B)
                cond_str = conditions_5B{i_cond}; lp_str = lightPeriods_5B{i_lp};
                animals_for_cond = allAnimals;
                if strcmp(cond_str, 'FullDark'), animals_for_cond = fourCondAnimals; end
                data_segment = dataTable(ismember(dataTable.Animal, animals_for_cond) & dataTable.Condition == cond_str & dataTable.LightPeriod == lp_str, :);
                if isempty(data_segment), continue; end
                animal_means = groupsummary(data_segment, {animalVar, 'Sex'}, 'mean', current_metric_var);
                animal_means.Group = repmat(categorical({sprintf('%s-%s', cond_str, lp_str)}), height(animal_means), 1);
                plotTable_5B = [plotTable_5B; animal_means];
                stat_col_name = sprintf('%s_%s', cond_str, lp_str);
                temp_stat_table = animal_means(:, {'Animal', ['mean_', current_metric_var]});
                temp_stat_table.Properties.VariableNames{2} = stat_col_name;
                all_animal_means_5B = outerjoin(all_animal_means_5B, temp_stat_table, 'Keys', 'Animal', 'MergeKeys', true);
            end
        end
        
        % Plotting
        groupOrder_5B = {'300Lux-On', '1000Lux-On', 'FullDark-On', '300Lux-Off', '1000Lux-Off', 'FullDark-Off'};
        finalGroupOrder = {};
        for i_grp = 1:length(groupOrder_5B)
            for i_sex = 1:length(sex_labels), finalGroupOrder{end+1} = sprintf('%s-%s', groupOrder_5B{i_grp}, sex_labels{i_sex}); end
        end
        plotTable_5B.PlotGroup = categorical(strcat(cellstr(plotTable_5B.Group), '-', cellstr(plotTable_5B.Sex)));
        violinplot(plotTable_5B.(['mean_', current_metric_var]), plotTable_5B.PlotGroup, 'GroupOrder', finalGroupOrder, 'ViolinColor', grayColor, 'ShowData', false, 'Parent', ax5B);
        
        % Scatter Overlay
        jitterAmount = 0.15;
        for i_plot_grp = 1:length(finalGroupOrder)
            group_label = finalGroupOrder{i_plot_grp};
            group_data = plotTable_5B(plotTable_5B.PlotGroup == group_label, :);
            if ~isempty(group_data)
                x_scatter = i_plot_grp + (rand(height(group_data), 1) - 0.5) * jitterAmount;
                if contains(group_label, 'Male'), scatter_color = maleColor; else, scatter_color = femaleColor; end
                scatter(ax5B, x_scatter, group_data.(['mean_', current_metric_var]), 40, scatter_color, 'filled', 'MarkerFaceAlpha', 0.9);
            end
        end
        ylabel(ax5B, ['Mean ' current_metric_ylabel]);
        title(ax5B, sprintf('Fig 5B: Diurnality by Condition and Sex (%s)', current_metric_suffix));
        xtickangle(ax5B, 45);
        set(ax5B, 'Position', [0.1 0.25 0.88 0.65]);
        
        % --- NEW: Descriptive Statistics Output ---
        fprintf('\n--- DESCRIPTIVE STATISTICS: Figure 5B (%s) ---\n', current_metric_suffix);
        fprintf('%-25s | %-10s | %-10s | %-10s\n', 'Group', 'Mean', 'Std Dev', 'CV (%)');
        fprintf('---------------------------------------------------------------\n');
        
        for i_grp = 1:length(finalGroupOrder)
            g_name = finalGroupOrder{i_grp};
            g_data = plotTable_5B.(['mean_', current_metric_var])(plotTable_5B.PlotGroup == g_name);
            if ~isempty(g_data)
                g_mean = mean(g_data, 'omitnan');
                g_std = std(g_data, 'omitnan');
                g_cv = (g_std / g_mean) * 100;
                fprintf('%-25s | %-10.2e | %-10.2e | %-10.1f%%\n', g_name, g_mean, g_std, g_cv);
            end
        end
        
        % --- Original Stats (Paired T-Tests) ---
        fprintf('\n--- STATS: Figure 5B (Paired On vs. Off by Sex/Condition) ---\n');
        tests_5B = { ...
            {'300Lux_On', '300Lux_Off', maleAnimals, 'Male 300Lux (On vs Off)'}, ...
            {'300Lux_On', '300Lux_Off', femaleAnimals, 'Female 300Lux (On vs Off)'}, ...
            {'1000Lux_On', '1000Lux_Off', maleAnimals, 'Male 1000Lux (On vs Off)'}, ...
            {'1000Lux_On', '1000Lux_Off', femaleAnimals, 'Female 1000Lux (On vs Off)'}, ...
            {'FullDark_On', 'FullDark_Off', intersect(maleAnimals, fourCondAnimals), 'Male FullDark (On vs Off)'}, ...
            {'FullDark_On', 'FullDark_Off', intersect(femaleAnimals, fourCondAnimals), 'Female FullDark (On vs Off)'} ...
        };
        num_tests_5B = length(tests_5B);
        alpha_5B = 0.05; % no bonferroni
        fprintf('Bonferroni Alpha: %.4f\n', alpha_5B);
        for i_test = 1:num_tests_5B
            test_info = tests_5B{i_test}; group_name = test_info{4};
            idx = ismember(all_animal_means_5B.Animal, test_info{3}); N = sum(idx);
            if N < 3, continue; end
            [~, p] = ttest(all_animal_means_5B.(test_info{1})(idx), all_animal_means_5B.(test_info{2})(idx));
            sig_str = 'n.s.'; if p < alpha_5B, sig_str = 'SIGNIFICANT'; end
            fprintf('%-30s | p = %.4f (N=%d) %s\n', group_name, p, N, sig_str);
        end

        stat_table_clean = all_animal_means_5B;
        stat_table_clean.Properties.VariableNames = matlab.lang.makeValidName(stat_table_clean.Properties.VariableNames);
        for i_sex = 1:length(sex_labels)
            sex_str = sex_labels{i_sex}; sex_animals = sex_groups{i_sex};
            t_wide_1 = stat_table_clean(ismember(string(stat_table_clean.Animal), string(sex_animals)), :);
            if height(t_wide_1) < 3, continue; end
            
            % 2-Way RM ANOVA
            withinDesign_1 = table({'x300Lux';'x300Lux';'x1000Lux';'x1000Lux'}, {'On';'Off';'On';'Off'}, 'VariableNames', {'Condition', 'LightPeriod'});
            rm_formula_1 = 'x300Lux_On-x1000Lux_Off ~ 1';
            rm_1 = fitrm(t_wide_1, rm_formula_1, 'WithinDesign', withinDesign_1);
            ranova_tbl_1 = ranova(rm_1, 'WithinModel', 'Condition*LightPeriod');
            
            % Extract Interaction P-Value
            p_idx_1 = contains(ranova_tbl_1.Properties.RowNames, 'Condition:LightPeriod');
            p_int = ranova_tbl_1.pValue(p_idx_1);
            fprintf('ANOVA (%s 300v1000): Condition*LightPeriod p = %.4f\n', sex_str, p_int(1));
        end

        fig_filename = fullfile(savePath_Fig5, sprintf('Figure5B_Violins_12group_%s.png', current_metric_suffix));
        exportgraphics(hFig5B_current, fig_filename); disp(['Saved: ', fig_filename]); close(hFig5B_current);
        
    catch ME_5B
        warning('ErrGen:Fig5B','Error Fig 5B (%s): %s', current_metric_suffix, ME_5B.message); 
    end

    %% FIG 5C: HOURLY DIFFERENCE COMPARISON (Male Delta vs Female Delta)
    disp('--- Starting Figure 5C: Hourly Difference Comparison (Male vs Female) ---');
    try
        hFig5C_Diff = figure('Name', sprintf('Fig 5C: Hourly Diff M vs F (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w', 'Position', [100, 100, 1200, 500]);
        ax5C = axes('Parent', hFig5C_Diff); hold(ax5C, 'on');
        
        % 1. Calculate Hourly Deltas
        sex_deltas = struct(); sex_deltas.Male = []; sex_deltas.Female = [];
        for i_sex = 1:2
            curr_group = sex_groups{i_sex}; sex_name = sex_labels{i_sex};
            group_delta_mat = NaN(length(curr_group), 24);
            for i_an = 1:length(curr_group)
                an_id = curr_group(i_an);
                d3=dataTable(dataTable.Animal==an_id&dataTable.Condition=='300Lux',:); h3=groupsummary(d3,ztHourVar,'mean',current_metric_var);
                d1=dataTable(dataTable.Animal==an_id&dataTable.Condition=='1000Lux',:); h1=groupsummary(d1,ztHourVar,'mean',current_metric_var);
                v3=NaN(1,24); v1=NaN(1,24);
                [~,l3]=ismember(0:23,h3.(ztHourVar)); v3(l3>0)=h3.(['mean_',current_metric_var])(l3(l3>0));
                [~,l1]=ismember(0:23,h1.(ztHourVar)); v1(l1>0)=h1.(['mean_',current_metric_var])(l1(l1>0));
                group_delta_mat(i_an, :) = v1 - v3;
            end
            sex_deltas.(sex_name) = group_delta_mat;
        end
        
        % 2. Means & SEMs
        male_mean = mean(sex_deltas.Male, 1, 'omitnan'); male_sem = std(sex_deltas.Male, 0, 1, 'omitnan') ./ sqrt(sum(~isnan(sex_deltas.Male)));
        fem_mean = mean(sex_deltas.Female, 1, 'omitnan'); fem_sem = std(sex_deltas.Female, 0, 1, 'omitnan') ./ sqrt(sum(~isnan(sex_deltas.Female)));
        
        % 3. Plot
        b = bar(ax5C, 0:23, [male_mean', fem_mean'], 'grouped');
        b(1).FaceColor = maleColor; b(1).EdgeColor = 'none'; b(2).FaceColor = femaleColor; b(2).EdgeColor = 'none';
        x_offset = [b(1).XEndPoints; b(2).XEndPoints];
        errorbar(ax5C, x_offset(1,:), male_mean, male_sem, 'k.', 'LineWidth', 1);
        errorbar(ax5C, x_offset(2,:), fem_mean, fem_sem, 'k.', 'LineWidth', 1);
        
        % 4. STATS with Verbose Output (KS Values Included)
        fprintf('\n--- STATS: Figure 5C (Male Delta vs Female Delta) ---\n');
        fprintf('   Comparing N=%d Males vs N=%d Females\n', size(sex_deltas.Male,1), size(sex_deltas.Female,1));
        p_vals_5C = ones(1, 24); sig_found = false;
        
        for zt = 1:24
            m_vals = sex_deltas.Male(:, zt); m_vals(isnan(m_vals))=[];
            f_vals = sex_deltas.Female(:, zt); f_vals(isnan(f_vals))=[];
            if length(m_vals)>=3 && length(f_vals)>=3
                % KS Normality Test (Normalize data first for valid KS test on small samples)
                if std(m_vals) > 0, [~, p_ks_m] = kstest((m_vals-mean(m_vals))/std(m_vals)); else, p_ks_m = 1; end
                if std(f_vals) > 0, [~, p_ks_f] = kstest((f_vals-mean(f_vals))/std(f_vals)); else, p_ks_f = 1; end
                
                % Print KS results for checking
                % fprintf('   ZT %d Normality: Male p=%.3f, Female p=%.3f -> ', zt-1, p_ks_m, p_ks_f);
                
                if p_ks_m > 0.05 && p_ks_f > 0.05
                    [~, p] = ttest2(m_vals, f_vals); test_type = 'T-Test';
                    % fprintf('Using T-Test\n');
                else
                    p = ranksum(m_vals, f_vals); test_type = 'RankSum';
                    % fprintf('Using RankSum\n');
                end
                p_vals_5C(zt) = p;
                
                % Report Significance
                if p < 0.05
                    fprintf('   ZT %d: SIGNIFICANT (%s) | p=%.4f (M:%.2f, F:%.2f) [KS: M=%.3f, F=%.3f]\n', ...
                        zt-1, test_type, p, mean(m_vals), mean(f_vals), p_ks_m, p_ks_f);
                    sig_found = true;
                end
            end
        end
        if ~sig_found, fprintf('   No significant hourly sex differences found.\n'); end
        
        % 5. Asterisks
        sig_idx = find(p_vals_5C < 0.05); y_limits = ylim(ax5C); y_max_plot = y_limits(2);
        for k = 1:length(sig_idx)
            zt_idx = sig_idx(k); h_star = max(male_mean(zt_idx)+male_sem(zt_idx), fem_mean(zt_idx)+fem_sem(zt_idx)) + 0.05*(y_limits(2)-y_limits(1));
            y_max_plot = max(y_max_plot, h_star);
            text(ax5C, zt_idx-1, h_star, '*', 'HorizontalAlignment', 'center', 'FontSize', 14, 'Color', 'k');
        end
        ylim(ax5C, [y_limits(1), y_max_plot + 0.1*(y_limits(2)-y_limits(1))]);
        
        % Shading
        yl_fin = ylim(ax5C);
        patch(ax5C, [ztLightsOffStart, ztLightsOffEnd+1, ztLightsOffEnd+1, ztLightsOffStart], ...
            [yl_fin(1) yl_fin(1) yl_fin(2) yl_fin(2)], [0.5 0.5 0.5], 'FaceAlpha', 0.1, 'EdgeColor', 'none'); 
            
        ylabel(ax5C, ['Mean Delta (1000L - 300L) ', current_metric_ylabel]); xlabel(ax5C, 'ZT Hour');
        title(ax5C, sprintf('Fig 5C: Hourly Response Difference (Male vs Female) (%s)', current_metric_suffix));
        legend(ax5C, {'Male \Delta', 'Female \Delta'}, 'Location', 'best'); xlim(ax5C, [-0.5 23.5]); xticks(ax5C, 0:23);
        exportgraphics(hFig5C_Diff, fullfile(savePath_Fig5, sprintf('Figure5C_HourlyDiff_Bar_%s.png', current_metric_suffix))); close(hFig5C_Diff);
    catch ME_5C, warning('ErrGen:Fig5C','Error Fig 5C (%s): %s', current_metric_suffix, ME_5C.message); end

    %% FIG 5D VARIATIONS: Smooth, Weekly, and Small Multiples (PixelDiff Only)
    % Purpose: Generate 3 variations of individual trajectory plots.
    
    if strcmp(current_metric_var, 'SelectedPixelDifference')
        disp('--- Starting Figure 5D Variations (Smooth, Weekly, Grid) ---');
        metric_to_plot = 'SelectedPixelDifference'; 
        ylabel_str = 'Daily Total Pixel Difference';
        
        % --- MANUAL Y-LIMIT SETTING ---
        manual_ylim = 4e8; 
        
        conditions_ordered = {'300Lux', '1000Lux', 'FullDark', '300LuxEnd'};
        sex_titles = {'Males', 'Females'};
        male_indiv_colors = parula(length(maleAnimals));
        female_indiv_colors = autumn(length(femaleAnimals));
        
        try
            % --- SHARED PRE-CALCULATIONS ---
            % 1. Condition Durations (for X-axis) - WITH 300LuxEnd TRUNCATION
            cond_durations = zeros(1, length(conditions_ordered));
            for c = 1:length(conditions_ordered)
                cond_idx = ismember(dataTable.Condition, conditions_ordered{c});
                cond_data = dataTable(cond_idx, :);
                if isempty(cond_data), continue; end
                
                if ismember('RelativeDay', cond_data.Properties.VariableNames)
                    cond_data.Day = floor(cond_data.RelativeDay);
                else, error('RelativeDay missing'); end
                
                % TRUNCATE 300LuxEnd to 7 days max for continuity logic
                if strcmp(conditions_ordered{c}, '300LuxEnd')
                    cond_data = cond_data(cond_data.Day <= 7, :); 
                end
                
                max_d = 0;
                u_ans = unique(cond_data.Animal);
                for ua = 1:length(u_ans)
                    d_list = cond_data.Day(cond_data.Animal == u_ans(ua));
                    days_pres = length(unique(d_list(~isnan(d_list))));
                    if days_pres > max_d, max_d = days_pres; end
                end
                cond_durations(c) = max_d;
            end
            x_end = cumsum(cond_durations);
            x_start = [1, x_end(1:end-1) + 1];
            x_mid = (x_start + x_end) / 2;
            
            
            %% VARIATION 1: SMOOTHED TRAJECTORIES (3-Day Moving Average)
            hFig_Smooth = figure('Name', 'Fig 5D_Smooth', 'Visible', 'off', 'Color', 'w', 'Position', [100, 100, 1000, 800]);
            sgtitle(hFig_Smooth, 'Figure 5D (Smooth): 3-Day Moving Average', 'FontSize', 14);
            
            for i_sex = 1:2
                ax = subplot(2, 1, i_sex, 'Parent', hFig_Smooth); hold(ax, 'on');
                curr_animals = sex_groups{i_sex};
                if i_sex == 1, c_map = male_indiv_colors; else, c_map = female_indiv_colors; end
                
                for i_an = 1:length(curr_animals)
                    an_id = curr_animals(i_an);
                    X_cat = []; Y_cat = [];
                    
                    for c = 1:length(conditions_ordered)
                        rows = (dataTable.Animal == an_id) & ismember(dataTable.Condition, conditions_ordered{c});
                        d = dataTable(rows, :);
                        d = d(~isnan(d.(metric_to_plot)), :);
                        if isempty(d), continue; end
                        d.Day = floor(d.RelativeDay);
                        
                        % APPLY FILTER: 300LuxEnd -> First 7 Days only
                        if strcmp(conditions_ordered{c}, '300LuxEnd')
                            d = d(d.Day <= 7, :);
                        end
                        
                        ds = groupsummary(d, 'Day', 'sum', metric_to_plot);
                        ds = sortrows(ds, 'Day');
                        
                        if height(ds)>0
                            x_pts = (x_start(c) : x_start(c) + height(ds) - 1)';
                            y_pts = ds.(['sum_', metric_to_plot]);
                            X_cat = [X_cat; x_pts]; Y_cat = [Y_cat; y_pts];
                        end
                    end
                    
                    if ~isempty(X_cat)
                        % Smoothing
                        Y_smooth = movmean(Y_cat, 3);
                        plot(ax, X_cat, Y_smooth, '-', 'Color', c_map(i_an,:), 'LineWidth', 2, 'DisplayName', char(an_id));
                    end
                end
                
                title(ax, sex_titles{i_sex}); ylabel(ax, ylabel_str);
                ylim(ax, [0, manual_ylim]); xlim(ax, [0, x_end(end)*1.05]);
                set(ax, 'XTick', x_mid, 'XTickLabel', conditions_ordered); grid(ax,'on');
                
                % Dividers (HandleVisibility off)
                yl=ylim(ax); 
                for c=1:length(conditions_ordered)-1
                    x_line=x_end(c)+0.5; 
                    line(ax,[x_line x_line],yl,'Color',[.4 .4 .4],'LineStyle','--', 'HandleVisibility', 'off'); 
                end
                
                % Legend
                legend(ax, 'Location', 'northeast', 'NumColumns', 1, 'FontSize', 8);
            end
            exportgraphics(hFig_Smooth, fullfile(savePath_Fig5, 'Figure5D_Smooth_PixelDiff.png')); close(hFig_Smooth);
    
            
            %% VARIATION 2: WEEKLY MEANS + SD (Visualizing Stability)
            hFig_Weekly = figure('Name', 'Fig 5D_Weekly', 'Visible', 'off', 'Color', 'w', 'Position', [100, 100, 1000, 800]);
            sgtitle(hFig_Weekly, 'Figure 5D (Weekly): Mean +/- SD per Week', 'FontSize', 14);
            
            for i_sex = 1:2
                ax = subplot(2, 1, i_sex, 'Parent', hFig_Weekly); hold(ax, 'on');
                curr_animals = sex_groups{i_sex};
                if i_sex == 1, c_map = male_indiv_colors; else, c_map = female_indiv_colors; end
                
                for i_an = 1:length(curr_animals)
                    an_id = curr_animals(i_an);
                    
                    for c = 1:length(conditions_ordered)
                        rows = (dataTable.Animal == an_id) & ismember(dataTable.Condition, conditions_ordered{c});
                        d = dataTable(rows, :);
                        d = d(~isnan(d.(metric_to_plot)), :);
                        if isempty(d), continue; end
                        d.Day = floor(d.RelativeDay);
                        
                        % APPLY FILTER: 300LuxEnd -> First 7 Days only (Results in 1 point)
                        if strcmp(conditions_ordered{c}, '300LuxEnd')
                            d = d(d.Day <= 7, :);
                        end
                        
                        ds = groupsummary(d, 'Day', 'sum', metric_to_plot);
                        
                        % Weekly Binning
                        num_days = height(ds);
                        num_weeks = ceil(num_days / 7);
                        
                        x_week = []; y_week = []; y_sd = [];
                        
                        for w = 1:num_weeks
                            start_d = (w-1)*7 + 1;
                            end_d = min(w*7, num_days);
                            idx = start_d:end_d;
                            w_data = ds.(['sum_', metric_to_plot])(idx);
                            
                            % X at center of week
                            x_center = x_start(c) + (start_d + end_d)/2 - 1;
                            x_week = [x_week; x_center];
                            y_week = [y_week; mean(w_data)];
                            y_sd   = [y_sd; std(w_data)];
                        end
                        
                        errorbar(ax, x_week, y_week, y_sd, '-o', 'Color', c_map(i_an,:), ...
                            'LineWidth', 1.5, 'MarkerSize', 5, 'MarkerFaceColor', c_map(i_an,:), ...
                            'DisplayName', char(an_id));
                    end
                end
                
                title(ax, sex_titles{i_sex}); ylabel(ax, ylabel_str);
                ylim(ax, [0, manual_ylim]); xlim(ax, [0, x_end(end)*1.05]);
                set(ax, 'XTick', x_mid, 'XTickLabel', conditions_ordered); grid(ax,'on');
                
                % Dividers (HandleVisibility off)
                yl=ylim(ax); 
                for c=1:length(conditions_ordered)-1
                    x_line=x_end(c)+0.5; 
                    line(ax,[x_line x_line],yl,'Color',[.4 .4 .4],'LineStyle','--', 'HandleVisibility', 'off'); 
                end
                
            end
            exportgraphics(hFig_Weekly, fullfile(savePath_Fig5, 'Figure5D_Weekly_PixelDiff.png')); close(hFig_Weekly);
    
            
            %% VARIATION 3: SMALL MULTIPLES (Grid Layout 6x2)
            hFig_Grid = figure('Name', 'Fig 5D_Grid', 'Visible', 'off', 'Color', 'w', 'Position', [100, 100, 1200, 1000]);
            sgtitle(hFig_Grid, 'Figure 5D (Individual): Trajectories', 'FontSize', 14);
            
            num_rows = max(length(maleAnimals), length(femaleAnimals));
            
            for r = 1:num_rows
                % -- PLOT MALE (Left Column) --
                if r <= length(maleAnimals)
                    ax = subplot(num_rows, 2, (r-1)*2 + 1, 'Parent', hFig_Grid); hold(ax,'on');
                    an_id = maleAnimals(r);
                    col = male_indiv_colors(r, :);
                    
                    X_cat=[]; Y_cat=[];
                    for c=1:length(conditions_ordered)
                        rows=(dataTable.Animal==an_id)&ismember(dataTable.Condition,conditions_ordered{c});
                        d=dataTable(rows,:); d=d(~isnan(d.(metric_to_plot)),:);
                        if ~isempty(d)
                            d.Day=floor(d.RelativeDay); 
                            % APPLY FILTER
                            if strcmp(conditions_ordered{c}, '300LuxEnd'), d=d(d.Day<=7, :); end
                            
                            ds=groupsummary(d,'Day','sum',metric_to_plot); ds=sortrows(ds,'Day');
                            if height(ds)>0
                                X_cat=[X_cat; (x_start(c):x_start(c)+height(ds)-1)'];
                                Y_cat=[Y_cat; ds.(['sum_',metric_to_plot])];
                            end
                        end
                    end
                    if ~isempty(X_cat)
                        plot(ax, X_cat, Y_cat, '-o', 'Color', col, 'LineWidth', 1.2, 'MarkerSize', 3, 'MarkerFaceColor', col);
                    end
                    
                    ylabel(ax, char(an_id), 'FontSize', 9, 'FontWeight','bold'); 
                    if r==1, title(ax, 'Males'); end
                    ylim(ax, [0, manual_ylim]); xlim(ax, [0, x_end(end)]);
                    
                    % Dividers
                    yl=ylim(ax); 
                    for c=1:length(conditions_ordered)-1
                        x_line=x_end(c)+0.5; 
                        line(ax,[x_line x_line],yl,'Color',[.8 .8 .8],'LineStyle','--', 'HandleVisibility', 'off'); 
                    end
                    
                    if r < num_rows, set(ax, 'XTick', []); else, set(ax, 'XTick', x_mid, 'XTickLabel', conditions_ordered); end
                end
                
                % -- PLOT FEMALE (Right Column) --
                if r <= length(femaleAnimals)
                    ax = subplot(num_rows, 2, (r-1)*2 + 2, 'Parent', hFig_Grid); hold(ax,'on');
                    an_id = femaleAnimals(r);
                    col = female_indiv_colors(r, :);
                    
                    X_cat=[]; Y_cat=[];
                    for c=1:length(conditions_ordered)
                        rows=(dataTable.Animal==an_id)&ismember(dataTable.Condition,conditions_ordered{c});
                        d=dataTable(rows,:); d=d(~isnan(d.(metric_to_plot)),:);
                        if ~isempty(d)
                            d.Day=floor(d.RelativeDay); 
                            % APPLY FILTER
                            if strcmp(conditions_ordered{c}, '300LuxEnd'), d=d(d.Day<=7, :); end
                            
                            ds=groupsummary(d,'Day','sum',metric_to_plot); ds=sortrows(ds,'Day');
                            if height(ds)>0
                                X_cat=[X_cat; (x_start(c):x_start(c)+height(ds)-1)'];
                                Y_cat=[Y_cat; ds.(['sum_',metric_to_plot])];
                            end
                        end
                    end
                    if ~isempty(X_cat)
                        plot(ax, X_cat, Y_cat, '-o', 'Color', col, 'LineWidth', 1.2, 'MarkerSize', 3, 'MarkerFaceColor', col);
                    end
                    
                    ylabel(ax, char(an_id), 'FontSize', 9, 'FontWeight','bold');
                    if r==1, title(ax, 'Females'); end
                    ylim(ax, [0, manual_ylim]); xlim(ax, [0, x_end(end)]);
                    
                    % Dividers
                    yl=ylim(ax); 
                    for c=1:length(conditions_ordered)-1
                        x_line=x_end(c)+0.5; 
                        line(ax,[x_line x_line],yl,'Color',[.8 .8 .8],'LineStyle','--', 'HandleVisibility', 'off'); 
                    end
                    
                    if r < num_rows, set(ax, 'XTick', []); else, set(ax, 'XTick', x_mid, 'XTickLabel', conditions_ordered); end
                end
            end
            
            exportgraphics(hFig_Grid, fullfile(savePath_Fig5, 'Figure5D_SmallMultiples_PixelDiff.png')); close(hFig_Grid);
            
        catch ME_5D_Var
            warning('ErrGen:Fig5D_Var','Error Fig 5D Variations: %s', ME_5D_Var.message);
        end
    else
        disp('Skipping Figure 5D Variations for Normalized Activity.');
    end

    %% FIG 5E: ON:OFF Activity Percent Change Plot (8-Violin Sex Comparison)
    disp('--- Starting Figure 5E: ON-OFF Pct Change by Sex ---');
    if strcmp(current_metric_var, normalizedActivityVar)
        disp('Skipping Figure 5E for NormalizedActivity (Percentage change invalid for Z-scores).');
    else
        try
        % 1. Setup Figure
        if ~exist('violinplot.m','file'),error('ViolinPlotNotFound:Fig5E','violinplot.m not found.');end
        hFig5E_current = figure('Name', sprintf('Figure 5E: ON-OFF Pct Change by Sex (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w');
        ax5E = axes('Parent', hFig5E_current); hold(ax5E, 'on');
        
        % 2. Data Preparation 
        pct_change_table = all_animal_means_5B; 
        pct_change_table.Sex = categorical(ismember(pct_change_table.Animal, maleAnimals), [0, 1], {'Female', 'Male'});
        
        pct_change_table.PctChg_300Lux = (pct_change_table.('300Lux_On') - pct_change_table.('300Lux_Off')) ./ pct_change_table.('300Lux_Off') * 100;
        pct_change_table.PctChg_1000Lux = (pct_change_table.('1000Lux_On') - pct_change_table.('1000Lux_Off')) ./ pct_change_table.('1000Lux_Off') * 100;
        pct_change_table.PctChg_FullDark = (pct_change_table.('FullDark_On') - pct_change_table.('FullDark_Off')) ./ pct_change_table.('FullDark_Off') * 100;
        
        if all(ismember({'300LuxEnd_On', '300LuxEnd_Off'}, pct_change_table.Properties.VariableNames))
            pct_change_table.PctChg_300LuxEnd = (pct_change_table.('300LuxEnd_On') - pct_change_table.('300LuxEnd_Off')) ./ pct_change_table.('300LuxEnd_Off') * 100;
            conditions_5E = {'300Lux', '1000Lux', 'FullDark', '300LuxEnd'};
        else
            warning('Fig5E:MissingData', '300LuxEnd data not found in table. Plotting 3 conditions only.');
            conditions_5E = {'300Lux', '1000Lux', 'FullDark'};
        end
        
        % --- (Data Prep for 8 Violins) ---
        plotData_5E = []; groupData_5E = []; groupOrder_5E = {};
        for i_cond = 1:length(conditions_5E)
            cond_name = conditions_5E{i_cond};
            cond_col_name = ['PctChg_', cond_name]; 
            animal_idx = true(height(pct_change_table), 1);
            if ismember(cond_name, ["FullDark", "300LuxEnd"])
                animal_idx = ismember(pct_change_table.Animal, fourCondAnimals);
            end
            for sex = ["Female", "Male"]
                group_name = [cond_name, '-', char(sex)];
                sex_idx = (pct_change_table.Sex == sex);
                final_idx = animal_idx & sex_idx & isfinite(pct_change_table.(cond_col_name));
                data_segment = pct_change_table(final_idx, :);
                if ~isempty(data_segment)
                    plotData_5E = [plotData_5E; data_segment.(cond_col_name)];
                    groupData_5E = [groupData_5E; repmat(categorical({group_name}), height(data_segment), 1)];
                    groupOrder_5E{end+1} = group_name;
                end
            end
        end
        
        % 3. Create Base Gray Violin Plot
        violinplot(plotData_5E, groupData_5E, 'Parent', ax5E, 'GroupOrder', groupOrder_5E, 'ViolinColor', grayColor, 'ShowData', false);
        
        % 4. Overlay Individual Data Points
        jitterAmount = 0.15;
        for i_group = 1:length(groupOrder_5E)
            group_label = groupOrder_5E{i_group};
            parts = split(group_label, '-'); cond_str = string(parts{1}); sex_str = string(parts{2});
            cond_col_name = ['PctChg_', char(cond_str)];
            animals_for_cond = allAnimals;
            if ismember(cond_str, ["FullDark", "300LuxEnd"]), animals_for_cond = fourCondAnimals; end
            data_idx = pct_change_table.Sex == sex_str & ...
                         ismember(pct_change_table.Animal, animals_for_cond) & ...
                         isfinite(pct_change_table.(cond_col_name));
            data_points = pct_change_table.(cond_col_name)(data_idx);
            if ~isempty(data_points)
                x_scatter = i_group + (rand(length(data_points), 1) - 0.5) * jitterAmount;
                scatterColor = (sex_str == "Male") * maleColor + (sex_str == "Female") * femaleColor;
                scatter(ax5E, x_scatter, data_points, 40, scatterColor, 'filled', 'MarkerFaceAlpha', 0.8, 'HandleVisibility', 'off');
            end
        end

        % 5. Reference line at 0%
        yline(ax5E, 0, 'k--');
        
        try
            disp('--- STATS: Figure 5E (M vs F Pct Change) ---');
            hold(ax5E, 'on');
            
            p_values_5E_MvF = containers.Map; % M vs F tests
            y_lims_5E = get(ax5E, 'YLim'); y_range_5E = diff(y_lims_5E);
            y_bar_level = y_lims_5E(2) + y_range_5E * 0.05; 
            y_text_level = y_bar_level + y_range_5E * 0.01; 
            new_y_max_5E = y_text_level + y_range_5E * 0.03;
            set(ax5E, 'YLim', [y_lims_5E(1), new_y_max_5E]); % Set YLim for first bar level

            for cond_name_str = conditions_5E
                cond_name = char(cond_name_str);
                cond_col_name = ['PctChg_', cond_name];
                animals_for_cond = allAnimals;
                if ismember(cond_name, ["FullDark", "300LuxEnd"]), animals_for_cond = fourCondAnimals; end
                
                male_data = pct_change_table.(cond_col_name)(pct_change_table.Sex=='Male' & ismember(pct_change_table.Animal, animals_for_cond) & isfinite(pct_change_table.(cond_col_name)));
                female_data = pct_change_table.(cond_col_name)(pct_change_table.Sex=='Female' & ismember(pct_change_table.Animal, animals_for_cond) & isfinite(pct_change_table.(cond_col_name)));
                
                if ~isempty(male_data) && ~isempty(female_data)
                    [~,p] = ttest2(male_data, female_data);
                    p_values_5E_MvF(cond_name) = p;
                    fprintf('Independent T-Test | %s PctChange (M vs F): p = %.4f\n', cond_name, p);
                else
                    fprintf('Skipping T-Test for %s: insufficient M or F data.\n', cond_name);
                end
            end

            % --- Plot M vs F bars (Level 1) ---
            for i_cond = 1:length(conditions_5E)
                cond_name = conditions_5E{i_cond};
                if isKey(p_values_5E_MvF, cond_name)
                    p = p_values_5E_MvF(cond_name);
                    x_female = find(strcmp(groupOrder_5E, [cond_name, '-Female']));
                    x_male = find(strcmp(groupOrder_5E, [cond_name, '-Male']));
                    if ~isempty(x_female) && ~isempty(x_male)
                        if p < 0.001, str_sig = '***'; elseif p < 0.01, str_sig = '**'; elseif p < 0.05, str_sig = '*'; else, str_sig = 'n.s.'; end
                        if ~strcmp(str_sig, 'n.s.')
                            plot(ax5E, [x_female, x_male], [y_bar_level, y_bar_level], 'k-', 'LineWidth', 1.2);
                            text(ax5E, mean([x_female, x_male]), y_text_level, str_sig, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 10, 'FontWeight', 'bold');
                        end
                    end
                end
            end
            
            if strcmp(current_metric_var, "SelectedPixelDifference")
                disp('--- STATS: Figure 5E (Cross-Condition Paired T-Tests) for PixelDiff ---');
                
                % Define new Y-levels based on the *current* max Y
                y_lims_2 = get(ax5E, 'YLim'); y_range_2 = diff(y_lims_2);
                y_bar_level_2 = y_lims_2(2) + y_range_2 * 0.05;
                y_text_level_2 = y_bar_level_2 + y_range_2 * 0.01;
                new_y_max_2 = y_text_level_2 + y_range_2 * 0.03;
                set(ax5E, 'YLim', [y_lims_2(1), new_y_max_2]); % Set new YLim for second bar level

                % 1. Female: 300Lux PctChg vs 1000Lux PctChg
                idx_f = ismember(pct_change_table.Animal, femaleAnimals);
                [~, p_f] = ttest(pct_change_table.PctChg_300Lux(idx_f), pct_change_table.PctChg_1000Lux(idx_f));
                fprintf('Paired T-Test | Female 300-PctChg vs 1000-PctChg: p = %.4f (N=%.0f)\n', p_f, sum(idx_f & ~isnan(pct_change_table.PctChg_300Lux) & ~isnan(pct_change_table.PctChg_1000Lux)));
                
                % 2. Male: 300Lux PctChg vs 1000Lux PctChg
                idx_m = ismember(pct_change_table.Animal, maleAnimals);
                [~, p_m] = ttest(pct_change_table.PctChg_300Lux(idx_m), pct_change_table.PctChg_1000Lux(idx_m));
                fprintf('Paired T-Test | Male 300-PctChg vs 1000-PctChg: p = %.4f (N=%.0f)\n', p_m, sum(idx_m & ~isnan(pct_change_table.PctChg_300Lux) & ~isnan(pct_change_table.PctChg_1000Lux)));

                % --- Plot Female Bar (Level 2) ---
                x_f_300 = find(strcmp(groupOrder_5E, '300Lux-Female'));
                x_f_1000 = find(strcmp(groupOrder_5E, '1000Lux-Female'));
                if ~isempty(x_f_300) && ~isempty(x_f_1000)
                    % Planned Comparison: Alpha = 0.05
                    if p_f < 0.05 
                        if p_f < 0.001, str = '***'; elseif p_f < 0.01, str = '**'; else, str = '*'; end
                        plot(ax5E, [x_f_300, x_f_1000], [y_bar_level_2, y_bar_level_2], 'k-', 'LineWidth', 1.2);
                        text(ax5E, mean([x_f_300, x_f_1000]), y_text_level_2, str, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 10, 'FontWeight', 'bold');
                    end
                end
                
                % --- Plot Male Bar (Level 2) ---
                x_m_300 = find(strcmp(groupOrder_5E, '300Lux-Male'));
                x_m_1000 = find(strcmp(groupOrder_5E, '1000Lux-Male'));
                if ~isempty(x_m_300) && ~isempty(x_m_1000)
                    % Planned Comparison: Alpha = 0.05
                    if p_m < 0.05 
                        if p_m < 0.001, str = '***'; elseif p_m < 0.01, str = '**'; else, str = '*'; end
                        plot(ax5E, [x_m_300, x_m_1000], [y_bar_level_2, y_bar_level_2], 'k-', 'LineWidth', 1.2);
                        text(ax5E, mean([x_m_300, x_m_1000]), y_text_level_2, str, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 10, 'FontWeight', 'bold');
                    end
                end
            end
            
        catch ME_Stats_5E
            warning('ErrGen:Fig5EStats', 'Could not plot significance for Fig 5E: %s', ME_Stats_5E.message);
        end

        % 8. Formatting and Saving
        title(ax5E, sprintf('Fig 5E: ON:OFF %% Change by Sex and Condition (%s)', current_metric_suffix));
        ylabel(ax5E, '% Change (ON vs OFF)');
        
        set(ax5E, 'XTick', 1:length(groupOrder_5E)); 
        set(ax5E, 'XTickLabel', groupOrder_5E);
        xtickangle(ax5E, 45); 
        
        grid(ax5E, 'on');
        
        fig_filename = fullfile(savePath_Fig5, sprintf('Figure5E_PctChangeViolins_BySex_%s.png', current_metric_suffix));
        exportgraphics(hFig5E_current, fig_filename); 
        disp(['Saved: ', fig_filename]); 
        close(hFig5E_current);
        
    catch ME_5E
        warning('ErrGen:Fig5E','Error Fig 5E (%s): %s', current_metric_suffix, ME_5E.message); end
    end
    %% FIG 5F: Peak Activity Shift (First vs Last Week) - CIRCULAR STATS
    disp('--- Starting Figure 5F: Peak Shift (First vs Last Week CoG) ---');
    if strcmp(current_metric_var, normalizedActivityVar)
        disp('Skipping Figure 5F for NormalizedActivity (Percentage change invalid for Z-scores).');
    else
        try
        hFig5F_peak = figure('Name', sprintf('Fig 5F: Peak Shift CoG (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w');
        ax5F = axes('Parent', hFig5F_peak); hold(ax5F, 'on');
        sex_colors = {maleColor, femaleColor}; 
        
        fprintf('\n--- STATS: Figure 5F (First Week vs Last Week Shift) ---\n');
        
        for i_sex = 1:length(sex_labels)
            sex_str = sex_labels{i_sex};
            sex_animals_all = sex_groups{i_sex};
            sex_color = sex_colors{i_sex};
            
            % Filter for animals that have FullDark data
            animals_to_test = intersect(sex_animals_all, fourCondAnimals);
            if isempty(animals_to_test), continue; end
            
            % Get FullDark Data
            data_FD_sex = dataTable(ismember(dataTable.Animal, animals_to_test) & dataTable.Condition == 'FullDark', :);
            unique_days_FD = unique(data_FD_sex.DayOfCondition);
            
            % Pre-allocate: Rows=Days, Cols=Animals
            animal_daily_peaks = NaN(length(unique_days_FD), length(animals_to_test));
            
            for i_an = 1:length(animals_to_test)
                an_id = animals_to_test(i_an);
                data_an = data_FD_sex(data_FD_sex.Animal == an_id, :);
                
                for i_day = 1:length(unique_days_FD)
                    d_val = unique_days_FD(i_day);
                    d_dat = data_an(data_an.DayOfCondition == d_val, :);
                    
                    if ~isempty(d_dat)
                        % CoG Calculation
                        h_mean = groupsummary(d_dat, ztHourVar, 'mean', current_metric_var);
                        if ~isempty(h_mean)
                            activity = h_mean.(['mean_', current_metric_var]); 
                            activity(activity < 0) = 0; % Remove negative weights
                            
                            if sum(activity) > 0
                                % Circular Center of Gravity
                                angles = h_mean.(ztHourVar) * (2*pi/24);
                                X = sum(activity .* cos(angles)); 
                                Y = sum(activity .* sin(angles));
                                mean_angle = atan2(Y, X);
                                peak_zt = mod(mean_angle * (24/(2*pi)), 24);
                                
                                animal_daily_peaks(i_day, i_an) = peak_zt;
                            end
                        end
                    end
                end
            end
            
            % --- PLOT DAILY TRAJECTORY (Visual Context) ---
            % Calculate daily circular means across animals for plotting
            daily_means_sex = NaN(length(unique_days_FD), 1);
            for d = 1:length(unique_days_FD)
                day_peaks = animal_daily_peaks(d, :);
                day_peaks = day_peaks(~isnan(day_peaks));
                if ~isempty(day_peaks)
                    rads = day_peaks * (2*pi/24);
                    mean_ang = atan2(mean(sin(rads)), mean(cos(rads)));
                    daily_means_sex(d) = mod(mean_ang * (24/(2*pi)), 24);
                end
            end
            plot(ax5F, unique_days_FD, daily_means_sex, 'o-', 'Color', sex_color, ...
                'LineWidth', 1.5, 'DisplayName', sex_str, 'MarkerFaceColor', sex_color);
            
            % --- STATS: COMPARE FIRST WEEK VS LAST WEEK ---
            % Define Weeks
            first_week_idx = 1:min(7, length(unique_days_FD));
            last_week_idx = max(1, length(unique_days_FD)-6):length(unique_days_FD);
            
            % Get Data subsets (CORRECTED VARIABLE NAMES)
            first_week_peaks = animal_daily_peaks(first_week_idx, :);
            last_week_peaks = animal_daily_peaks(last_week_idx, :);
            
            % Calculate Circular Mean PER ANIMAL for each block
            % Function to average angles:
            calc_circ_mean_cols = @(x) mod(atan2(mean(sin(x*2*pi/24),1,'omitnan'), mean(cos(x*2*pi/24),1,'omitnan')) * 24/(2*pi), 24);
            
            an_start_means = calc_circ_mean_cols(first_week_peaks); % 1 x N_animals
            an_end_means = calc_circ_mean_cols(last_week_peaks);    % 1 x N_animals
            
            % Calculate Shift (Shortest distance)
            % Shift = End - Start
            shifts = mod(an_end_means - an_start_means + 12, 24) - 12;
            
            % One-Sample T-Test vs 0 (Equivalent to Paired T-Test)
            [~, p_shift] = ttest(shifts, 0);
            
            % Grand Means for Reporting
            grand_start = mod(atan2(mean(sin(an_start_means*2*pi/24)), mean(cos(an_start_means*2*pi/24))) * 24/(2*pi), 24);
            grand_end = mod(atan2(mean(sin(an_end_means*2*pi/24)), mean(cos(an_end_means*2*pi/24))) * 24/(2*pi), 24);
            mean_shift_val = mean(shifts, 'omitnan');
            
            fprintf('%s Comparison (N=%d):\n', sex_str, length(shifts));
            fprintf('  First Week Mean ZT: %.2f\n', grand_start);
            fprintf('  Last Week Mean ZT:  %.2f\n', grand_end);
            fprintf('  Mean Shift:         %.2f hours\n', mean_shift_val);
            fprintf('  Paired T-Test:      p = %.4f\n', p_shift);
        end
        
        xlabel(ax5F, 'Day in Darkness'); 
        ylabel(ax5F, 'Peak Activity ZT (CoG)');
        title(ax5F, sprintf('Fig 5F: Daily Peak Activity (%s)', current_metric_suffix));
        legend(ax5F, 'Location', 'best');
        ylim(ax5F, [0 24]); yticks(ax5F, 0:4:24);
        grid(ax5F, 'on');
        
        fig_filename = fullfile(savePath_Fig5, sprintf('Figure5F_FirstVsLastWeek_%s.png', current_metric_suffix));
        exportgraphics(hFig5F_peak, fig_filename); disp(['Saved: ', fig_filename]); close(hFig5F_peak);
        
    catch ME_5F
        warning('ErrGen:Fig5F','Error Fig 5F (%s): %s', current_metric_suffix, ME_5F.message); 
    end
    end
    %% FIG 5G: Phase-Specific Percent Change (MID-PHASE ONLY)
    % Comparison: % Change in Mid-Light (ZT 3-9) vs Mid-Dark (ZT 15-21) 
    % phases for Males vs Females across conditions.
    disp('--- Starting Figure 5G: Phase-Specific Pct Change (Mid-Phase) ---');
    try
        hFig5G_current = figure('Name', sprintf('Fig 5G: Mid-Phase Pct Change (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w');
        ax5G = axes('Parent', hFig5G_current); hold(ax5G, 'on');
        
        % Comparison Definitions: {Name, BaseCond, TargCond, WindowType}
        comparisons = {
            '1000L-MidLight', '300Lux', '1000Lux', 'MidLight';
            'Dark-MidLight',  '300Lux', 'FullDark', 'MidLight';
            '1000L-MidDark',  '300Lux', '1000Lux', 'MidDark';
            'Dark-MidDark',   '300Lux', 'FullDark', 'MidDark';
        };
        
        plotTable_5G = table();
        
        % Calculate Data
        for i_comp = 1:size(comparisons, 1)
            comp_label = comparisons{i_comp, 1};
            base_cond  = comparisons{i_comp, 2};
            targ_cond  = comparisons{i_comp, 3};
            phase_type = comparisons{i_comp, 4};
            
            % If Dark is involved, use N=8 subset
            if contains(targ_cond, 'Dark'), animals_to_use = fourCondAnimals; else, animals_to_use = allAnimals; end
            
            for i_an = 1:length(animals_to_use)
                an_id = animals_to_use(i_an);
                
                % Get Baseline Data (ZT Filtered)
                base_data = dataTable(dataTable.Animal == an_id & dataTable.Condition == base_cond, :);
                if strcmp(phase_type, 'MidLight')
                    val_base = mean(base_data.(current_metric_var)(base_data.(ztHourVar) >= 3 & base_data.(ztHourVar) <= 9), 'omitnan');
                else
                    val_base = mean(base_data.(current_metric_var)(base_data.(ztHourVar) >= 15 & base_data.(ztHourVar) <= 21), 'omitnan');
                end
                
                % Get Target Data (ZT Filtered)
                targ_data = dataTable(dataTable.Animal == an_id & dataTable.Condition == targ_cond, :);
                if strcmp(phase_type, 'MidLight')
                    val_targ = mean(targ_data.(current_metric_var)(targ_data.(ztHourVar) >= 3 & targ_data.(ztHourVar) <= 9), 'omitnan');
                else
                    val_targ = mean(targ_data.(current_metric_var)(targ_data.(ztHourVar) >= 15 & targ_data.(ztHourVar) <= 21), 'omitnan');
                end
                
                % Calc % Change
                if ~isnan(val_base) && val_base ~= 0 && ~isnan(val_targ)
                    pct_chg = ((val_targ - val_base) / val_base) * 100;
                    
                    if ismember(an_id, maleAnimals), sex = 'Male'; else, sex = 'Female'; end
                    
                    % Store
                    new_row = table({an_id}, string(comp_label), string(sex), pct_chg, 'VariableNames', {'Animal', 'Comparison', 'Sex', 'PctChange'});
                    plotTable_5G = [plotTable_5G; new_row];
                end
            end
        end
        
        % Plot Violins
        plotTable_5G.Group = categorical(strcat(plotTable_5G.Comparison, '-', plotTable_5G.Sex));
        
        grp_order = {};
        for i = 1:size(comparisons, 1)
            grp_order{end+1} = [comparisons{i,1}, '-Female'];
            grp_order{end+1} = [comparisons{i,1}, '-Male'];
        end
        
        violinplot(plotTable_5G.PctChange, plotTable_5G.Group, 'Parent', ax5G, 'GroupOrder', grp_order, 'ShowData', false, 'ViolinColor', grayColor);
        yline(ax5G, 0, 'k:');
        
        % Scatter Overlay (Explicit Colors)
        jitterAmount = 0.15;
        for i_g = 1:length(grp_order)
            g_name = grp_order{i_g};
            g_data = plotTable_5G(plotTable_5G.Group == g_name, :);
            if ~isempty(g_data)
                x_scat = i_g + (rand(height(g_data),1)-0.5)*jitterAmount;
                if contains(g_name, 'Male'), col = maleColor; else, col = femaleColor; end
                scatter(ax5G, x_scat, g_data.PctChange, 40, col, 'filled', 'MarkerFaceAlpha', 0.9);
            end
        end
        
        % Stats (Male vs Female for each comp)
        disp('--- STATS: Figure 5G (Sex Differences in Mid-Phase Response) ---');
        y_lims = ylim(ax5G); y_range = diff(y_lims);
        y_pos = y_lims(2) + 0.05*y_range;
        
        for i_c = 1:size(comparisons, 1)
            comp_name = string(comparisons{i_c, 1});
            
            data_m = plotTable_5G.PctChange(plotTable_5G.Comparison == comp_name & plotTable_5G.Sex == "Male");
            data_f = plotTable_5G.PctChange(plotTable_5G.Comparison == comp_name & plotTable_5G.Sex == "Female");
            
            if ~isempty(data_m) && ~isempty(data_f)
                [~, p] = ttest2(data_m, data_f);
                fprintf('Independent T-Test | %s (M vs F): p = %.4f\n', comp_name, p);
                
                if p < 0.05
                    x_f = find(strcmp(grp_order, [char(comp_name), '-Female']));
                    x_m = find(strcmp(grp_order, [char(comp_name), '-Male']));
                    plot_sig_bar(ax5G, [x_f, x_m], p, y_pos, 0.05);
                end
            end
        end
        ylim(ax5G, [y_lims(1), y_pos + 0.1*y_range]);
        
        ylabel(ax5G, '% Change from 300Lux');
        title(ax5G, sprintf('Mid-Phase Response by Sex (%s)', strrep(current_metric_suffix,'_',' ')));
        xtickangle(ax5G, 45);
        set(ax5G, 'Position', [0.1 0.25 0.88 0.65]);
        
        fig_filename = fullfile(savePath_Fig5, sprintf('Figure5G_MidPhaseResponse_%s.png', current_metric_suffix));
        exportgraphics(hFig5G_current, fig_filename); disp(['Saved: ', fig_filename]); close(hFig5G_current);
        
    catch ME_5G
        warning('ErrGen:Fig5G','Error Fig 5G (%s): %s', current_metric_suffix, ME_5G.message); 
    end

    %% FIG 5G_New: ZT 0-5 Percent Change (Violin Plot)
    if strcmp(current_metric_var, 'SelectedPixelDifference')
        disp('--- Starting Figure 5G_New: ZT 0-5 Morning Masking % Change ---');
        try
            hFig5G_New = figure('Name', 'Fig 5G_New: ZT 0-5 % Change', 'Visible', 'off', 'Color', 'w');
            ax5G = axes('Parent', hFig5G_New); hold(ax5G, 'on');
            
            % Config
            sexes_5G = {'Male', 'Female'};
            morning_data = dataTable(dataTable.(ztHourVar) >= 0 & dataTable.(ztHourVar) <= 5, :);
            
            % Fallback for grayColor if not in workspace
            if ~exist('grayColor', 'var'), grayColor = [0.6 0.6 0.6]; end 
            
            plotTable_5G = table();
            
            % Data Collection
            for i_sex = 1:length(sexes_5G)
                sex_name = sexes_5G{i_sex};
                if strcmp(sex_name, 'Male'), current_animals = maleAnimals; else, current_animals = femaleAnimals; end
                
                for i_an = 1:length(current_animals)
                    an_id = current_animals(i_an);
                    
                    % 300Lux
                    d300 = morning_data(morning_data.Animal == an_id & morning_data.Condition == '300Lux', :);
                    val_300 = mean(d300.(current_metric_var), 'omitnan');
                    
                    % 1000Lux
                    d1000 = morning_data(morning_data.Animal == an_id & morning_data.Condition == '1000Lux', :);
                    val_1000 = mean(d1000.(current_metric_var), 'omitnan');
                    
                    if ~isnan(val_300) && ~isnan(val_1000) && val_300 > 0
                        pct_change = ((val_1000 - val_300) / val_300) * 100;
                        
                        % Append rows using STRING casting to prevent 'unique' error
                        row_pct = table(string(an_id), string(sex_name), pct_change, ...
                            'VariableNames', {'Animal', 'Sex', 'Value'});
                            
                        plotTable_5G = [plotTable_5G; row_pct];
                    end
                end
            end
            
            % Plotting
            grp_order = {'Male', 'Female'};
            plotTable_5G.Group = categorical(plotTable_5G.Sex, grp_order);
            
            violinplot(plotTable_5G.Value, plotTable_5G.Group, 'Parent', ax5G, 'GroupOrder', grp_order, 'ShowData', false, 'ViolinColor', grayColor);
            
            % Scatter Overlay (Matches 5H style)
            jitterAmount = 0.15;
            for i_g = 1:length(grp_order)
                g_name = grp_order{i_g};
                g_data = plotTable_5G(plotTable_5G.Group == g_name, :);
                if ~isempty(g_data)
                    x_scat = i_g + (rand(height(g_data),1)-0.5)*jitterAmount;
                    if contains(g_name, 'Male'), col = maleColor; else, col = femaleColor; end
                    scatter(ax5G, x_scat, g_data.Value, 40, col, 'filled', 'MarkerFaceAlpha', 0.9);
                end
            end
            
            % Stats: Independent Test (Male vs Female Pct Change)
            fprintf('\n--- STATS: Figure 5G_New (ZT 0-5 Pct Change) ---\n');
            m_data = plotTable_5G.Value(plotTable_5G.Group == 'Male');
            f_data = plotTable_5G.Value(plotTable_5G.Group == 'Female');
            
            y_lims = ylim(ax5G); y_range = diff(y_lims);
            y_pos = max(plotTable_5G.Value) + 0.10*y_range;
            
            if length(m_data) >= 3 && length(f_data) >= 3
                if std(m_data)>0, [~, pk_m] = kstest((m_data-mean(m_data))/std(m_data)); else, pk_m=1; end
                if std(f_data)>0, [~, pk_f] = kstest((f_data-mean(f_data))/std(f_data)); else, pk_f=1; end
                
                if pk_m > 0.05 && pk_f > 0.05
                    [~, p] = ttest2(m_data, f_data); test_name = 'T-Test';
                else
                    p = ranksum(m_data, f_data); test_name = 'RankSum';
                end
                
                fprintf('%s | Male vs Female Pct Change: p = %.4f\n', test_name, p);
                
                % Significance Markers (+ for trending, * for sig)
                if p < 0.10
                    plot(ax5G, [1, 2], [y_pos, y_pos], '-k', 'LineWidth', 1.5);
                    if p < 0.05
                        text(ax5G, 1.5, y_pos + 0.02*y_range, '*', 'FontSize', 20, 'HorizontalAlignment', 'center');
                    else
                        text(ax5G, 1.5, y_pos + 0.02*y_range, '+', 'FontSize', 16, 'HorizontalAlignment', 'center');
                    end
                end
            end
            
            ylim(ax5G, [y_lims(1), y_pos + 0.15*y_range]); 
            ylabel(ax5G, 'Percent Change (%) from 300 Lux (ZT 0-5)');
            title(ax5G, 'Morning Masking Reactivity (ZT 0-5)');
            set(ax5G, 'Position', [0.15 0.15 0.75 0.75]); % Slightly adjust padding
            
            fig_filename = fullfile(savePath_Fig5, 'Figure5G_ZT0to5_PctChange_Violin.png');
            exportgraphics(hFig5G_New, fig_filename); disp(['Saved: ', fig_filename]); close(hFig5G_New);
            
        catch ME_5G
            warning('ErrGen:Fig5G_New','Error Fig 5G_New: %s', ME_5G.message);
        end
    else
        disp('Skipping Figure 5G_New for Normalized Activity.');
    end

    %% FIG 5H: Mid-Phase Diurnality by Sex (FIXED TYPE ERROR)
    % Purpose: Test diurnality using only stable "Mid-Day" (ZT 3-9) vs "Mid-Night" (ZT 15-21) data.
    disp('--- Starting Figure 5H: Mid-Phase Diurnality ---');
    try
        hFig5H_current = figure('Name', sprintf('Fig 5H: Mid-Phase Diurnality (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w');
        ax5H = axes('Parent', hFig5H_current); hold(ax5H, 'on');
        
        % Config
        conds_5H = {'300Lux', '1000Lux', 'FullDark'};
        windows_5H = {'MidLight (3-9)', 'MidDark (15-21)'};
        sexes_5H = {'Male', 'Female'};
        
        plotTable_5H = table();
        
        % Data Collection
        for i_cond = 1:length(conds_5H)
            cond_name = conds_5H{i_cond};
            animals_for_cond = allAnimals;
            if strcmp(cond_name, 'FullDark'), animals_for_cond = fourCondAnimals; end
            
            for i_sex = 1:length(sexes_5H)
                sex_name = sexes_5H{i_sex};
                if strcmp(sex_name, 'Male'), current_animals = maleAnimals; else, current_animals = femaleAnimals; end
                
                % Intersect lists
                target_animals = intersect(animals_for_cond, current_animals);
                
                for i_an = 1:length(target_animals)
                    an_id = target_animals(i_an);
                    an_data = dataTable(dataTable.Animal == an_id & dataTable.Condition == cond_name, :);
                    
                    % Mid-Light (3-9)
                    idx_ml = an_data.(ztHourVar) >= 3 & an_data.(ztHourVar) <= 9;
                    val_ml = mean(an_data.(current_metric_var)(idx_ml), 'omitnan');
                    
                    % Mid-Dark (15-21)
                    idx_md = an_data.(ztHourVar) >= 15 & an_data.(ztHourVar) <= 21;
                    val_md = mean(an_data.(current_metric_var)(idx_md), 'omitnan');
                    
                    % Append rows using STRING casting to prevent 'unique' error
                    row_ml = table(string(an_id), string(sex_name), string(cond_name), string('MidLight'), val_ml, ...
                        'VariableNames', {'Animal', 'Sex', 'Condition', 'Phase', 'Value'});
                    row_md = table(string(an_id), string(sex_name), string(cond_name), string('MidDark'), val_md, ...
                        'VariableNames', {'Animal', 'Sex', 'Condition', 'Phase', 'Value'});
                        
                    plotTable_5H = [plotTable_5H; row_ml; row_md];
                end
            end
        end
        
        % Plotting
        grp_order = {};
        % Males First
        for i_c = 1:length(conds_5H)
            grp_order{end+1} = [conds_5H{i_c}, '-MidLight-Male'];
            grp_order{end+1} = [conds_5H{i_c}, '-MidDark-Male'];
        end
        % Females Second
        for i_c = 1:length(conds_5H)
            grp_order{end+1} = [conds_5H{i_c}, '-MidLight-Female'];
            grp_order{end+1} = [conds_5H{i_c}, '-MidDark-Female'];
        end
        
        plotTable_5H.Group = categorical(strcat(plotTable_5H.Condition, '-', plotTable_5H.Phase, '-', plotTable_5H.Sex));
        violinplot(plotTable_5H.Value, plotTable_5H.Group, 'Parent', ax5H, 'GroupOrder', grp_order, 'ShowData', false, 'ViolinColor', grayColor);
        
        % Scatter Overlay
        jitterAmount = 0.15;
        for i_g = 1:length(grp_order)
            g_name = grp_order{i_g};
            g_data = plotTable_5H(plotTable_5H.Group == g_name, :);
            if ~isempty(g_data)
                x_scat = i_g + (rand(height(g_data),1)-0.5)*jitterAmount;
                if contains(g_name, 'Male'), col = maleColor; else, col = femaleColor; end
                scatter(ax5H, x_scat, g_data.Value, 40, col, 'filled', 'MarkerFaceAlpha', 0.9);
            end
        end
        
        % Stats: Paired T-Test (MidLight vs MidDark) per Condition/Sex
        fprintf('\n--- STATS: Figure 5H (Mid-Phase Diurnality) ---\n');
        y_lims = ylim(ax5H); y_range = diff(y_lims);
        y_pos = y_lims(2) + 0.05*y_range;
        
        alpha_5H = 0.05;
        fprintf('Using Planned Comparison Alpha: %.2f\n', alpha_5H);
        
        piv_5H = unstack(plotTable_5H, 'Value', 'Phase', 'GroupingVariables', {'Animal', 'Sex', 'Condition'});
        
        for i_s = 1:length(sexes_5H)
            for i_c = 1:length(conds_5H)
                sex = sexes_5H{i_s}; cond = conds_5H{i_c};
                
                % Filter data
                d_sub = piv_5H(piv_5H.Sex == sex & piv_5H.Condition == cond, :);
                if height(d_sub) < 3, continue; end
                
                [~, p] = ttest(d_sub.MidLight, d_sub.MidDark);
                fprintf('Paired T-Test | %s %s (ML vs MD): p = %.4f\n', sex, cond, p);
                
                if p < alpha_5H
                    % Find X coords
                    x1 = find(strcmp(grp_order, [cond, '-MidLight-', sex]));
                    x2 = find(strcmp(grp_order, [cond, '-MidDark-', sex]));
                    if ~isempty(x1) && ~isempty(x2)
                        plot_sig_bar(ax5H, [x1, x2], p, y_pos, alpha_5H);
                    end
                end
            end
        end
        
        ylim(ax5H, [y_lims(1), y_pos + 0.15*y_range]); 
        ylabel(ax5H, ['Mean ', current_metric_ylabel]);
        title(ax5H, sprintf('Mid-Phase Diurnality (Stable 6h Windows) (%s)', strrep(current_metric_suffix,'_',' ')));
        xtickangle(ax5H, 45);
        set(ax5H, 'Position', [0.1 0.25 0.88 0.65]);
        
        fig_filename = fullfile(savePath_Fig5, sprintf('Figure5H_MidPhaseDiurnality_%s.png', current_metric_suffix));
        exportgraphics(hFig5H_current, fig_filename); disp(['Saved: ', fig_filename]); close(hFig5H_current);
        
    catch ME_5H
        warning('ErrGen:Fig5H','Error Fig 5H (%s): %s', current_metric_suffix, ME_5H.message); 
    end

    %% FIG 5I: Individual Trajectories & Variance Analysis (NEW)
    % Purpose: Visualize individual consistency (Slope) and test for differences in Variance (F-Test).
    disp('--- Starting Figure 5I: Individual Trajectories & Variance ---');
    try
        hFig5I_current = figure('Name', sprintf('Fig 5I: Individual Trajectories (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w');
        
        % Subplot 1: Females
        ax5I_f = subplot(1, 2, 1, 'Parent', hFig5I_current); hold(ax5I_f, 'on');
        title(ax5I_f, 'Females: 300L \rightarrow 1000L');
        xticks(ax5I_f, [1, 2]); xticklabels(ax5I_f, {'300 Lux', '1000 Lux'});
        xlim(ax5I_f, [0.5, 2.5]); ylabel(ax5I_f, current_metric_ylabel);
        
        % Subplot 2: Males
        ax5I_m = subplot(1, 2, 2, 'Parent', hFig5I_current); hold(ax5I_m, 'on');
        title(ax5I_m, 'Males: 300L \rightarrow 1000L');
        xticks(ax5I_m, [1, 2]); xticklabels(ax5I_m, {'300 Lux', '1000 Lux'});
        xlim(ax5I_m, [0.5, 2.5]);
        
        % Data Containers for Variance Test
        vals_f_1000 = [];
        vals_m_1000 = [];
        
        % --- Plotting Loop ---
        % FEMALES
        for i = 1:length(femaleAnimals)
            an = femaleAnimals(i);
            d300 = dataTable(dataTable.Animal == an & dataTable.Condition == '300Lux', :).(current_metric_var);
            d1000 = dataTable(dataTable.Animal == an & dataTable.Condition == '1000Lux', :).(current_metric_var);
            
            val300 = mean(d300, 'omitnan');
            val1000 = mean(d1000, 'omitnan');
            
            if ~isnan(val300) && ~isnan(val1000)
                plot(ax5I_f, [1, 2], [val300, val1000], 'o-', 'Color', femaleColor, 'LineWidth', 1.5, 'MarkerFaceColor', femaleColor);
                vals_f_1000 = [vals_f_1000; val1000];
            end
        end
        
        % MALES
        for i = 1:length(maleAnimals)
            an = maleAnimals(i);
            d300 = dataTable(dataTable.Animal == an & dataTable.Condition == '300Lux', :).(current_metric_var);
            d1000 = dataTable(dataTable.Animal == an & dataTable.Condition == '1000Lux', :).(current_metric_var);
            
            val300 = mean(d300, 'omitnan');
            val1000 = mean(d1000, 'omitnan');
            
            if ~isnan(val300) && ~isnan(val1000)
                plot(ax5I_m, [1, 2], [val300, val1000], 'o-', 'Color', maleColor, 'LineWidth', 1.5, 'MarkerFaceColor', maleColor);
                vals_m_1000 = [vals_m_1000; val1000];
            end
        end
        
        % Link Axes to same Y-scale for fair comparison
        linkaxes([ax5I_f, ax5I_m], 'y');
        
        % --- STATS: F-Test for Equality of Variances ---
        fprintf('\n--- STATS: Figure 5I (Variance Analysis) ---\n');
        
        % vartest2 performs an F-test (requires Normal distribution, but robust enough for this purpose)
        % H0: Variances are equal. H1: Variances are different.
        if length(vals_m_1000) > 2 && length(vals_f_1000) > 2
            [h_var, p_var, ci_var, stats_var] = vartest2(vals_m_1000, vals_f_1000);
            
            fprintf('F-Test for Equal Variances (Males vs Females in 1000 Lux):\n');
            fprintf('  Male Std Dev:   %.2e\n', std(vals_m_1000));
            fprintf('  Female Std Dev: %.2e\n', std(vals_f_1000));
            fprintf('  F-Statistic:    %.2f\n', stats_var.fstat);
            fprintf('  P-Value:        %.4f\n', p_var);
            
            if p_var < 0.05
                fprintf('  RESULT: Variances are SIGNIFICANTLY different.\n');
                subtitle(ax5I_m, sprintf('Male Variance Sig. Higher (p=%.3f)', p_var));
            else
                fprintf('  RESULT: Variances are not significantly different.\n');
            end
        else
            fprintf('Skipping Variance Test (N too low).\n');
        end
        
        fig_filename = fullfile(savePath_Fig5, sprintf('Figure5I_Trajectories_%s.png', current_metric_suffix));
        exportgraphics(hFig5I_current, fig_filename); disp(['Saved: ', fig_filename]); close(hFig5I_current);
        
    catch ME_5I
        warning('ErrGen:Fig5I','Error Fig 5I (%s): %s', current_metric_suffix, ME_5I.message); 
    end

    %% FIG 5J: Masking Effect in Mid-Phases by Sex (NEW VIOLINS)
    % Purpose: Direct paired comparison of 300Lux vs 1000Lux within stable Mid-Light (ML) 
    % and Mid-Dark (MD) windows, separated by sex.
    disp('--- Starting Figure 5J: Mid-Phase Masking Analysis ---');
    try
        hFig5J_current = figure('Name', sprintf('Fig 5J: Mid-Phase Masking (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w', 'Position', [100, 100, 1000, 600]);
        ax5J = axes('Parent', hFig5J_current); hold(ax5J, 'on');
        
        % Config
        conds_to_compare = {'300Lux', '1000Lux'}; % Only comparing these two
        windows_5J = {'MidLight', 'MidDark'}; % ZT 3-9 and 15-21
        sexes_5J = {'Male', 'Female'};
        
        plotTable_5J = table();
        
        % --- 1. Data Collection Loop ---
        for i_c = 1:length(conds_to_compare)
            cond_name = conds_to_compare{i_c};
            
            for i_s = 1:length(sexes_5J)
                sex_name = sexes_5J{i_s};
                if strcmp(sex_name, 'Male'), current_animals = maleAnimals; else, current_animals = femaleAnimals; end
                
                for i_an = 1:length(current_animals)
                    an_id = current_animals(i_an);
                    % Get data for this animal in this condition
                    an_data = dataTable(dataTable.Animal == an_id & dataTable.Condition == cond_name, :);
                    if isempty(an_data), continue; end
                    
                    % Calculate Mid-Light (ZT 3-9) Mean
                    idx_ml = an_data.(ztHourVar) >= 3 & an_data.(ztHourVar) <= 9;
                    val_ml = mean(an_data.(current_metric_var)(idx_ml), 'omitnan');
                    
                    % Calculate Mid-Dark (ZT 15-21) Mean
                    idx_md = an_data.(ztHourVar) >= 15 & an_data.(ztHourVar) <= 21;
                    val_md = mean(an_data.(current_metric_var)(idx_md), 'omitnan');
                    
                    % Append rows (using string casting for safety)
                    row_ml = table(string(an_id), string(sex_name), string(cond_name), string('MidLight'), val_ml, ...
                        'VariableNames', {'Animal', 'Sex', 'Condition', 'Phase', 'Value'});
                    row_md = table(string(an_id), string(sex_name), string(cond_name), string('MidDark'), val_md, ...
                        'VariableNames', {'Animal', 'Sex', 'Condition', 'Phase', 'Value'});
                        
                    plotTable_5J = [plotTable_5J; row_ml; row_md];
                end
            end
        end
        
        % --- 2. Plotting Setup ---
        % Define the specific order for paired comparisons
        grp_order_5J = { ...
            'Male-MidLight-300Lux', 'Male-MidLight-1000Lux', ...
            'Male-MidDark-300Lux',  'Male-MidDark-1000Lux', ...
            'Female-MidLight-300Lux', 'Female-MidLight-1000Lux', ...
            'Female-MidDark-300Lux',  'Female-MidDark-1000Lux'};
            
        % Create composite group variable
        plotTable_5J.Group = categorical(strcat(plotTable_5J.Sex, '-', plotTable_5J.Phase, '-', plotTable_5J.Condition));
        
        % Draw Violins
        violinplot(plotTable_5J.Value, plotTable_5J.Group, 'Parent', ax5J, 'GroupOrder', grp_order_5J, 'ShowData', false, 'ViolinColor', grayColor);
        
        % Scatter Overlay (colored by sex)
        jitterAmount = 0.15;
        for i_g = 1:length(grp_order_5J)
            g_name = grp_order_5J{i_g};
            g_data = plotTable_5J(plotTable_5J.Group == g_name, :);
            if ~isempty(g_data)
                x_scat = i_g + (rand(height(g_data),1)-0.5)*jitterAmount;
                if contains(g_name, 'Male'), col = maleColor; else, col = femaleColor; end
                scatter(ax5J, x_scat, g_data.Value, 40, col, 'filled', 'MarkerFaceAlpha', 0.9);
            end
        end
        
        % --- 3. Stats: Paired T-Tests (300 vs 1000 within Phase/Sex) ---
        fprintf('\n--- STATS: Figure 5J (Mid-Phase Masking: 300 vs 1000) ---\n');
        y_lims = ylim(ax5J); y_range = diff(y_lims);
        y_pos_bar = y_lims(2) + 0.05*y_range;
        
        % Planned comparisons: Alpha = 0.05
        alpha_5J = 0.05;
        
        % Pivot data wide to put 300 and 1000 side-by-side for paired tests
        piv_5J = unstack(plotTable_5J, 'Value', 'Condition', 'GroupingVariables', {'Animal', 'Sex', 'Phase'});
        
        % Loop through the 4 comparisons
        comparisons = {
            'Male', 'MidLight';
            'Male', 'MidDark';
            'Female', 'MidLight';
            'Female', 'MidDark'};
            
        for i = 1:size(comparisons, 1)
            curr_sex = comparisons{i,1}; curr_phase = comparisons{i,2};
            
            % Filter data for this pair
            d_sub = piv_5J(strcmp(piv_5J.Sex, curr_sex) & strcmp(piv_5J.Phase, curr_phase), :);
            
            % Ensure paired data (N>=3)
            if height(d_sub) >= 3 && all(ismember({'x300Lux', 'x1000Lux'}, d_sub.Properties.VariableNames))
                 % Paired T-Test
                 [~, p] = ttest(d_sub.x300Lux, d_sub.x1000Lux);
                 fprintf('Paired T-Test | %s %s (300 vs 1000): p = %.4f (N=%d)\n', curr_sex, curr_phase, p, height(d_sub));
                 
                 % Plot Significance Bar if needed
                 if p < alpha_5J
                    g1 = [curr_sex, '-', curr_phase, '-300Lux'];
                    g2 = [curr_sex, '-', curr_phase, '-1000Lux'];
                    x1 = find(strcmp(grp_order_5J, g1));
                    x2 = find(strcmp(grp_order_5J, g2));
                    plot_sig_bar(ax5J, [x1, x2], p, y_pos_bar, alpha_5J);
                 end
            end
        end
        
        % Formatting
        ylim(ax5J, [y_lims(1), y_pos_bar + 0.15*y_range]); 
        ylabel(ax5J, ['Mean ', current_metric_ylabel]);
        title(ax5J, sprintf('Masking Effect in Mid-Phases (300L vs 1000L) (%s)', strrep(current_metric_suffix,'_',' ')));
        xtickangle(ax5J, 45);
        set(ax5J, 'Position', [0.1 0.25 0.88 0.65]);
        grid(ax5J, 'on');
        
        fig_filename = fullfile(savePath_Fig5, sprintf('Figure5J_MidPhaseMasking_%s.png', current_metric_suffix));
        exportgraphics(hFig5J_current, fig_filename); disp(['Saved: ', fig_filename]); close(hFig5J_current);
        
    catch ME_5J
        warning('ErrGen:Fig5J','Error Fig 5J (%s): %s', current_metric_suffix, ME_5J.message); 
    end

    %% FIG 5K: Transition Masking Analysis (Dawn & Dusk) by Sex
    % Purpose: Direct paired comparison of 300Lux vs 1000Lux within 
    % Dawn (ZT 22-2) and Dusk (ZT 10-14) transition windows.
    disp('--- Starting Figure 5K: Transition Masking Analysis ---');
    try
        hFig5K_current = figure('Name', sprintf('Fig 5K: Transition Masking (%s)', current_metric_suffix), 'Visible', 'off', 'Color', 'w', 'Position', [100, 100, 1000, 600]);
        ax5K = axes('Parent', hFig5K_current); hold(ax5K, 'on');
        
        % Config
        conds_to_compare = {'300Lux', '1000Lux'};
        sexes_5K = {'Male', 'Female'};
        
        plotTable_5K = table();
        
        % --- 1. Data Collection Loop ---
        for i_c = 1:length(conds_to_compare)
            cond_name = conds_to_compare{i_c};
            
            for i_s = 1:length(sexes_5K)
                sex_name = sexes_5K{i_s};
                if strcmp(sex_name, 'Male'), current_animals = maleAnimals; else, current_animals = femaleAnimals; end
                
                % Use animals that have data for BOTH conditions (intersection with allAnimals)
                target_animals = intersect(current_animals, allAnimals);
                
                for i_an = 1:length(target_animals)
                    an_id = target_animals(i_an);
                    an_data = dataTable(dataTable.Animal == an_id & dataTable.Condition == cond_name, :);
                    if isempty(an_data), continue; end
                    
                    % --- Calculate Dusk Transition (ZT 10-14) ---
                    idx_dusk = an_data.(ztHourVar) >= 10 & an_data.(ztHourVar) <= 14;
                    val_dusk = mean(an_data.(current_metric_var)(idx_dusk), 'omitnan');
                    
                    % --- Calculate Dawn Transition (ZT 22-2) [WRAP AROUND] ---
                    idx_dawn = an_data.(ztHourVar) >= 22 | an_data.(ztHourVar) <= 2;
                    val_dawn = mean(an_data.(current_metric_var)(idx_dawn), 'omitnan');
                    
                    % Append rows (using string casting for safety)
                    row_dusk = table(string(an_id), string(sex_name), string(cond_name), string('DuskTrans'), val_dusk, ...
                        'VariableNames', {'Animal', 'Sex', 'Condition', 'Phase', 'Value'});
                    row_dawn = table(string(an_id), string(sex_name), string(cond_name), string('DawnTrans'), val_dawn, ...
                        'VariableNames', {'Animal', 'Sex', 'Condition', 'Phase', 'Value'});
                        
                    plotTable_5K = [plotTable_5K; row_dusk; row_dawn];
                end
            end
        end
        
        % --- 2. Plotting Setup ---
        % Define Order: Male Dusk -> Male Dawn -> Female Dusk -> Female Dawn
        grp_order_5K = { ...
            'Male-DuskTrans-300Lux', 'Male-DuskTrans-1000Lux', ...
            'Male-DawnTrans-300Lux', 'Male-DawnTrans-1000Lux', ...
            'Female-DuskTrans-300Lux', 'Female-DuskTrans-1000Lux', ...
            'Female-DawnTrans-300Lux', 'Female-DawnTrans-1000Lux'};
            
        % Create composite group variable
        plotTable_5K.Group = categorical(strcat(plotTable_5K.Sex, '-', plotTable_5K.Phase, '-', plotTable_5K.Condition));
        
        % Draw Violins
        violinplot(plotTable_5K.Value, plotTable_5K.Group, 'Parent', ax5K, 'GroupOrder', grp_order_5K, 'ShowData', false, 'ViolinColor', grayColor);
        
        % Scatter Overlay (Explicit Colors)
        jitterAmount = 0.15;
        for i_g = 1:length(grp_order_5K)
            g_name = grp_order_5K{i_g};
            g_data = plotTable_5K(plotTable_5K.Group == g_name, :);
            
            if ~isempty(g_data)
                x_scat = i_g + (rand(height(g_data),1)-0.5)*jitterAmount;
                
                % Assign Color
                if contains(g_name, 'Male'), col = maleColor; else, col = femaleColor; end
                
                scatter(ax5K, x_scat, g_data.Value, 40, col, 'filled', 'MarkerFaceAlpha', 0.9);
            end
        end
        
        % --- 3. Stats: Paired T-Tests (300 vs 1000 within Phase/Sex) ---
        fprintf('\n--- STATS: Figure 5K (Transition Masking: 300 vs 1000) ---\n');
        y_lims = ylim(ax5K); y_range = diff(y_lims);
        y_pos_bar = y_lims(2) + 0.05*y_range;
        
        alpha_5K = 0.05;
        
        % Pivot data wide
        piv_5K = unstack(plotTable_5K, 'Value', 'Condition', 'GroupingVariables', {'Animal', 'Sex', 'Phase'});
        
        % Loop through the 4 comparisons
        comparisons = {
            'Male', 'DuskTrans';
            'Male', 'DawnTrans';
            'Female', 'DuskTrans';
            'Female', 'DawnTrans'};
            
        for i = 1:size(comparisons, 1)
            curr_sex = comparisons{i,1}; curr_phase = comparisons{i,2};
            
            d_sub = piv_5K(strcmp(piv_5K.Sex, curr_sex) & strcmp(piv_5K.Phase, curr_phase), :);
            
            if height(d_sub) >= 3 && all(ismember({'x300Lux', 'x1000Lux'}, d_sub.Properties.VariableNames))
                 [~, p] = ttest(d_sub.x300Lux, d_sub.x1000Lux);
                 fprintf('Paired T-Test | %s %s (300 vs 1000): p = %.4f (N=%d)\n', curr_sex, curr_phase, p, height(d_sub));
                 
                 if p < alpha_5K
                    g1 = [curr_sex, '-', curr_phase, '-300Lux'];
                    g2 = [curr_sex, '-', curr_phase, '-1000Lux'];
                    x1 = find(strcmp(grp_order_5K, g1));
                    x2 = find(strcmp(grp_order_5K, g2));
                    plot_sig_bar(ax5K, [x1, x2], p, y_pos_bar, alpha_5K);
                 end
            end
        end
        
        % Formatting
        ylim(ax5K, [y_lims(1), y_pos_bar + 0.15*y_range]); 
        ylabel(ax5K, ['Mean ', current_metric_ylabel]);
        title(ax5K, sprintf('Transition Masking (Dusk 10-14, Dawn 22-2) (%s)', strrep(current_metric_suffix,'_',' ')));
        xtickangle(ax5K, 45);
        set(ax5K, 'Position', [0.1 0.25 0.88 0.65]);
        grid(ax5K, 'on');
        
        fig_filename = fullfile(savePath_Fig5, sprintf('Figure5K_TransitionMasking_%s.png', current_metric_suffix));
        exportgraphics(hFig5K_current, fig_filename); disp(['Saved: ', fig_filename]); close(hFig5K_current);
        
    catch ME_5K
        warning('ErrGen:Fig5K','Error Fig 5K (%s): %s', current_metric_suffix, ME_5K.message); 
    end

    %% FIG 5L: Morning Masking (ZT 0-5) Analysis
    if strcmp(current_metric_var, 'SelectedPixelDifference')
        disp('--- Starting Figure 5L: ZT 0-5 Morning Masking Analysis ---');
        try
            hFig5L = figure('Name', 'Fig 5L: ZT 0-5 Morning Masking', 'Visible', 'off', 'Color', 'w', 'Position', [100, 100, 1000, 500]);
            
            % Filter for morning hours only (ZT 0 through 5)
            morning_data = dataTable(dataTable.(ztHourVar) >= 0 & dataTable.(ztHourVar) <= 5, :);
            
            raw_1000 = struct('Male', [], 'Female', []);
            pct_changes = struct('Male', [], 'Female', []);
            
            for i_sex = 1:2
                curr_animals = sex_groups{i_sex};
                sex_str = sex_labels{i_sex};
                
                for i_an = 1:length(curr_animals)
                    an_id = curr_animals(i_an);
                    
                    % Get 300Lux Morning Mean
                    d300 = morning_data(morning_data.Animal == an_id & morning_data.Condition == '300Lux', :);
                    mean_300 = NaN;
                    if ~isempty(d300) && ~all(isnan(d300.(current_metric_var)))
                        mean_300 = mean(d300.(current_metric_var), 'omitnan');
                    end
                    
                    % Get 1000Lux Morning Mean
                    d1000 = morning_data(morning_data.Animal == an_id & morning_data.Condition == '1000Lux', :);
                    mean_1000 = NaN;
                    if ~isempty(d1000) && ~all(isnan(d1000.(current_metric_var)))
                        mean_1000 = mean(d1000.(current_metric_var), 'omitnan');
                        raw_1000.(sex_str)(end+1) = mean_1000; % Store raw 1000L data
                    end
                    
                    % Calculate Percent Change
                    if ~isnan(mean_300) && ~isnan(mean_1000) && mean_300 > 0
                        pct_change = ((mean_1000 - mean_300) / mean_300) * 100;
                        pct_changes.(sex_str)(end+1) = pct_change;
                    end
                end
            end
            
            % --- Subplot 1: Direct Compare 1000 Lux (ZT 0-5) ---
            ax1 = subplot(1, 2, 1, 'Parent', hFig5L); hold(ax1, 'on');
            m_raw = raw_1000.Male; f_raw = raw_1000.Female;
            
            bar(ax1, 1, mean(m_raw), 'FaceColor', maleColor, 'EdgeColor', 'none', 'BarWidth', 0.6);
            bar(ax1, 2, mean(f_raw), 'FaceColor', femaleColor, 'EdgeColor', 'none', 'BarWidth', 0.6);
            
            % Add individual data points
            plot(ax1, 1 + randn(size(m_raw))*0.05, m_raw, 'k.', 'MarkerSize', 15);
            plot(ax1, 2 + randn(size(f_raw))*0.05, f_raw, 'k.', 'MarkerSize', 15);
            
            ylabel(ax1, 'Mean Pixel Difference (ZT 0-5)');
            xticks(ax1, [1, 2]); xticklabels(ax1, {'Males', 'Females'});
            title(ax1, '1000 Lux Absolute Activity');
            
            % Stats for Subplot 1
            fprintf('\n--- STATS: ZT 0-5 Direct Compare 1000 Lux (M vs F) ---\n');
            if length(m_raw) >= 3 && length(f_raw) >= 3
                if std(m_raw)>0, [~, pk_m] = kstest((m_raw-mean(m_raw))/std(m_raw)); else, pk_m=1; end
                if std(f_raw)>0, [~, pk_f] = kstest((f_raw-mean(f_raw))/std(f_raw)); else, pk_f=1; end
                
                if pk_m > 0.05 && pk_f > 0.05, [~, p_raw] = ttest2(m_raw, f_raw); t_name = 'T-Test';
                else, p_raw = ranksum(m_raw, f_raw); t_name = 'RankSum'; end
                
                fprintf('   %s | Male Mean: %.1f vs Female Mean: %.1f | p = %.4f\n', t_name, mean(m_raw), mean(f_raw), p_raw);
                if p_raw < 0.05, text(ax1, 1.5, max([m_raw, f_raw])*1.05, '*', 'FontSize', 16, 'HorizontalAlignment', 'center'); end
            end
            
            % --- Subplot 2: Percent Change (300 to 1000 Lux) ---
            ax2 = subplot(1, 2, 2, 'Parent', hFig5L); hold(ax2, 'on');
            m_pct = pct_changes.Male; f_pct = pct_changes.Female;
            
            bar(ax2, 1, mean(m_pct), 'FaceColor', maleColor, 'EdgeColor', 'none', 'BarWidth', 0.6);
            bar(ax2, 2, mean(f_pct), 'FaceColor', femaleColor, 'EdgeColor', 'none', 'BarWidth', 0.6);
            
            plot(ax2, 1 + randn(size(m_pct))*0.05, m_pct, 'k.', 'MarkerSize', 15);
            plot(ax2, 2 + randn(size(f_pct))*0.05, f_pct, 'k.', 'MarkerSize', 15);
            
            ylabel(ax2, 'Percent Change (%)');
            xticks(ax2, [1, 2]); xticklabels(ax2, {'Males', 'Females'});
            title(ax2, 'Percent Change (300L to 1000L)');
            
            % Stats for Subplot 2
            fprintf('\n--- STATS: ZT 0-5 Percent Change (M vs F) ---\n');
            if length(m_pct) >= 3 && length(f_pct) >= 3
                if std(m_pct)>0, [~, pk_m] = kstest((m_pct-mean(m_pct))/std(m_pct)); else, pk_m=1; end
                if std(f_pct)>0, [~, pk_f] = kstest((f_pct-mean(f_pct))/std(f_pct)); else, pk_f=1; end
                
                if pk_m > 0.05 && pk_f > 0.05, [~, p_pct] = ttest2(m_pct, f_pct); t_name = 'T-Test';
                else, p_pct = ranksum(m_pct, f_pct); t_name = 'RankSum'; end
                
                fprintf('   %s | Male Mean: %.1f%% vs Female Mean: %.1f%% | p = %.4f\n', t_name, mean(m_pct), mean(f_pct), p_pct);
                if p_pct < 0.05, text(ax2, 1.5, max([m_pct, f_pct])*1.05, '*', 'FontSize', 16, 'HorizontalAlignment', 'center'); end
            end
            
            sgtitle(hFig5L, 'Figure 5L: Morning Masking (ZT 0-5) by Sex', 'FontSize', 14);
            fig_name_5L = fullfile(savePath_Fig5, 'Figure5L_MorningMasking_ZT0to5_PixelDiff.png');
            exportgraphics(hFig5L, fig_name_5L); disp(['Saved: ', fig_name_5L]); close(hFig5L);
            
        catch ME_5L
            warning('ErrGen:Fig5L','Error Fig 5L: %s', ME_5L.message);
        end
    else
        disp('Skipping Figure 5L for Normalized Activity.');
    end
end

%% FIG 6: Full Dark Free-Running Analysis
current_metric_var = 'SelectedPixelDifference';
if strcmp(current_metric_var, 'SelectedPixelDifference')
    disp('--- Starting Figure 6: Free-Running Period & Circadian Time Profiles ---');
    try
        % 1. Extract FullDark Data
        fd_data = dataTable(dataTable.Condition == 'FullDark', :);
        
        if isempty(fd_data)
            error('No FullDark data found for Free-Running analysis.');
        end
        
        % Ensure RelativeDay exists
        if ~ismember('RelativeDay', fd_data.Properties.VariableNames)
            error('Column "RelativeDay" is missing.');
        end
        
        % Initialize storage for tau
        tau_results = struct('Animal', {}, 'Sex', {}, 'Tau', {});
        ct_data_compiled = table(); % Will hold data re-aligned to CT
        
        % Search range for periodogram (22 to 26 hours)
        tau_search_range = linspace(22, 26, 2000); 
        f_vec = 1 ./ tau_search_range;
        
        % --- STEP 1: Calculate Tau for each animal ---
        fprintf('DEBUG: Calculating Free-Running Period (Tau) via Lomb-Scargle...\n');
        u_animals = unique(fd_data.Animal);
        
        for i_an = 1:length(u_animals)
            an_id = u_animals(i_an);
            an_data = fd_data(fd_data.Animal == an_id, :);
            
            % Determine Sex
            if ismember(an_id, maleAnimals), sex_str = 'Male';
            elseif ismember(an_id, femaleAnimals), sex_str = 'Female';
            else, continue; end
            
            % Create continuous elapsed hours vector
            % Subtract min to start at 0, multiply by 24 for hours
            t_hours = (an_data.RelativeDay - min(an_data.RelativeDay)) * 24;
            y_vals = an_data.(current_metric_var);
            
            % Remove NaNs for plomb
            valid_idx = ~isnan(y_vals) & ~isnan(t_hours);
            t_clean = t_hours(valid_idx);
            y_clean = y_vals(valid_idx);
            
            if length(y_clean) < 48 % Need at least 2 days of data
                fprintf('   Skipping Animal %s (Insufficient valid data)\n', char(an_id));
                continue;
            end
            
            % Mean-center data
            y_clean = y_clean - mean(y_clean);
            
            % Run Periodogram
            [pxx, f] = plomb(y_clean, t_clean, f_vec);
            [~, peak_idx] = max(pxx);
            best_tau = 1 / f(peak_idx);
            
            % Save Tau
            tau_results(end+1) = struct('Animal', an_id, 'Sex', sex_str, 'Tau', best_tau);
            fprintf('   Animal %s (%s): Tau = %.2f hrs\n', char(an_id), sex_str, best_tau);
            
            % --- STEP 2: Re-align to Circadian Time (CT) ---
            % Formula: Phase = mod(Elapsed, Tau) * (24/Tau)
            ct_float = mod(t_clean, best_tau) * (24 / best_tau);
            ct_hour = floor(ct_float); % Map to integer bins 0-23
            
            % Store in compiled table
            temp_t = table(repmat(an_id, length(y_clean), 1), repmat({sex_str}, length(y_clean), 1), ...
                           ct_hour, y_clean + mean(y_vals(valid_idx)), ... % add mean back for plotting magnitude
                           'VariableNames', {'Animal', 'Sex', 'CT_Hour', 'Metric'});
            ct_data_compiled = [ct_data_compiled; temp_t];
        end
        
        
        %% FIG 6A: Compare Tau by Sex ---
        hFig6A = figure('Name', 'Fig 6A: Free-Running Period (Tau)', 'Visible', 'off', 'Color', 'w', 'Position', [100, 100, 600, 500]);
        ax6A = axes('Parent', hFig6A); hold(ax6A, 'on');
        
        tau_males = [tau_results(strcmp({tau_results.Sex}, 'Male')).Tau];
        tau_females = [tau_results(strcmp({tau_results.Sex}, 'Female')).Tau];
        
        % Plot Bars/Points
        bar(ax6A, 1, mean(tau_males), 'FaceColor', maleColor, 'EdgeColor', 'none', 'BarWidth', 0.6);
        bar(ax6A, 2, mean(tau_females), 'FaceColor', femaleColor, 'EdgeColor', 'none', 'BarWidth', 0.6);
        
        plot(ax6A, 1 + randn(size(tau_males))*0.05, tau_males, 'k.', 'MarkerSize', 15);
        plot(ax6A, 2 + randn(size(tau_females))*0.05, tau_females, 'k.', 'MarkerSize', 15);
        
        % Formatting
        ylabel(ax6A, 'Free-Running Period \tau (Hours)');
        xticks(ax6A, [1, 2]); xticklabels(ax6A, {'Males', 'Females'});
        ylim(ax6A, [23, 25]); % Zoom in on circadian range
        title(ax6A, 'Fig 6A: Intrinsic Clock Speed in Full Dark');
        
        % Stats
        fprintf('\n--- STATS: Figure 6A (Tau Comparison) ---\n');
        if length(tau_males) >= 3 && length(tau_females) >= 3
            [~, p_ks_m] = kstest((tau_males-mean(tau_males))/std(tau_males));
            [~, p_ks_f] = kstest((tau_females-mean(tau_females))/std(tau_females));
            if p_ks_m > 0.05 && p_ks_f > 0.05
                [~, p_tau] = ttest2(tau_males, tau_females); test_name = 'T-Test';
            else
                p_tau = ranksum(tau_males, tau_females); test_name = 'RankSum';
            end
            fprintf('   %s | M: %.2f hrs vs F: %.2f hrs | p = %.4f\n', test_name, mean(tau_males), mean(tau_females), p_tau);
            
            % Add sig star if needed
            if p_tau < 0.05
                y_star = max([tau_males, tau_females]) + 0.2;
                plot(ax6A, [1, 2], [y_star, y_star], '-k', 'LineWidth', 1.5);
                text(ax6A, 1.5, y_star + 0.05, '*', 'FontSize', 16, 'HorizontalAlignment', 'center');
                ylim(ax6A, [23, y_star + 0.3]);
            end
        else
            fprintf('   Insufficient N for Tau stats.\n');
        end
        
        fig_name_6A = fullfile(savePath_Fig5, 'Figure6A_Tau_Comparison.png'); % Saving in Fig5 folder to keep together
        exportgraphics(hFig6A, fig_name_6A); disp(['Saved: ', fig_name_6A]); close(hFig6A);
        
        
        %% FIG 6B: Circadian Time (CT) Profiles ---
        hFig6B = figure('Name', 'Fig 6B: Circadian Time Profiles', 'Visible', 'off', 'Color', 'w', 'Position', [100, 100, 800, 500]);
        ax6B = axes('Parent', hFig6B); hold(ax6B, 'on');
        
        % Plot Males
        m_ct = ct_data_compiled(strcmp(ct_data_compiled.Sex, 'Male'), :);
        if ~isempty(m_ct)
            m_stats = groupsummary(m_ct, 'CT_Hour', {'mean', 'std', 'nnz'}, 'Metric');
            m_mean = m_stats.mean_Metric;
            m_sem = m_stats.std_Metric ./ sqrt(m_stats.nnz_Metric);
            
            plot(ax6B, m_stats.CT_Hour, m_mean, 'Color', maleColor, 'LineWidth', 2, 'DisplayName', 'Males');
            fill(ax6B, [m_stats.CT_Hour; flipud(m_stats.CT_Hour)], [m_mean-m_sem; flipud(m_mean+m_sem)], ...
                maleColor, 'FaceAlpha', 0.2, 'EdgeColor', 'none', 'HandleVisibility', 'off');
        end
        
        % Plot Females
        f_ct = ct_data_compiled(strcmp(ct_data_compiled.Sex, 'Female'), :);
        if ~isempty(f_ct)
            f_stats = groupsummary(f_ct, 'CT_Hour', {'mean', 'std', 'nnz'}, 'Metric');
            f_mean = f_stats.mean_Metric;
            f_sem = f_stats.std_Metric ./ sqrt(f_stats.nnz_Metric);
            
            plot(ax6B, f_stats.CT_Hour, f_mean, 'Color', femaleColor, 'LineWidth', 2, 'DisplayName', 'Females');
            fill(ax6B, [f_stats.CT_Hour; flipud(f_stats.CT_Hour)], [f_mean-f_sem; flipud(f_mean+f_sem)], ...
                femaleColor, 'FaceAlpha', 0.2, 'EdgeColor', 'none', 'HandleVisibility', 'off');
        end
        
        % Formatting
        title(ax6B, 'Fig 6B: Full Dark Activity Aligned to Circadian Time (CT)');
        xlabel(ax6B, 'Circadian Time (CT Hour)');
        ylabel(ax6B, 'Mean SelectedPixelDifference');
        xlim(ax6B, [-0.5, 23.5]); xticks(ax6B, 0:6:23);
        
        % Shading Subjective Night (Usually CT 12-24)
        yl = ylim(ax6B);
        patch(ax6B, [12, 24, 24, 12], [yl(1) yl(1) yl(2) yl(2)], [0.9 0.9 0.9], ...
            'FaceAlpha', 0.3, 'EdgeColor', 'none', 'HandleVisibility', 'off');
        
        legend(ax6B, 'Location', 'best'); grid(ax6B, 'on');
        
        fig_name_6B = fullfile(savePath_Fig5, 'Figure6B_CT_Profiles.png');
        exportgraphics(hFig6B, fig_name_6B); disp(['Saved: ', fig_name_6B]); close(hFig6B);
        
    catch ME_6
        warning('ErrGen:Fig6','Error Fig 6 (Tau/CT): %s', ME_6.message);
    end

    %% FIG 6C: DIAGNOSTIC - ZT vs CT Comparison for Full Dark (Week 4)
    % Purpose: Directly answer whether aligning to CT changes the statistical 
    % outcome of DD diurnality compared to standard ZT, highlighting AO5.
    
    if strcmp(current_metric_var, 'SelectedPixelDifference')
        disp('--- Starting Diagnostic: ZT vs CT in Last Week of Full Dark ---');
        try
            % 1. Isolate Last 7 Days of Full Dark
            fd_data = dataTable(dataTable.Condition == 'FullDark', :);
            
            % We need tau_results from Fig 6. If missing, throw error.
            if ~exist('tau_results', 'var')
                error('tau_results not found. Run Figure 6A/B first.');
            end
            
            m_zt_day = []; m_zt_night = []; m_ct_day = []; m_ct_night = [];
            f_zt_day = []; f_zt_night = []; f_ct_day = []; f_ct_night = [];
            
            ao5_stats = struct(); % Store AO5 specifically for Lily
            
            u_ans = unique(fd_data.Animal);
            for i = 1:length(u_ans)
                an_id = u_ans(i);
                d = fd_data(fd_data.Animal == an_id, :);
                d.Day = floor(d.RelativeDay);
                max_day = max(d.Day);
                d_last_week = d(d.Day > (max_day - 7), :); % Last 7 days
                
                if isempty(d_last_week), continue; end
                
                % --- FIX: Robust String Matching for Tau ---
                tau_animals = string({tau_results.Animal});
                tau_idx = (tau_animals == string(an_id));
                
                if ~any(tau_idx)
                    fprintf('   Warning: Could not find Tau match for Animal %s. Skipping...\n', string(an_id));
                    continue; 
                end
                an_tau = tau_results(tau_idx).Tau;
                
                % Calculate CT for last week
                % 1. Get all FullDark data for this animal to find the true start time
                an_fd_data = fd_data(fd_data.Animal == an_id, :);
                start_time = min(an_fd_data.RelativeDay);
                
                % 2. Find the exact ZT Hour when the FullDark condition started
                start_zt_array = an_fd_data.(ztHourVar)(an_fd_data.RelativeDay == start_time);
                start_zt = start_zt_array(1); % Extract the first value
                
                % 3. Calculate elapsed hours since the start of FullDark
                t_hours = (d_last_week.RelativeDay - start_time) * 24;
                
                % 4. Add the starting ZT phase to the elapsed time
                ct_float = mod(t_hours + start_zt, an_tau) * (24 / an_tau);
                d_last_week.CT_Hour = floor(ct_float);
                
                % --- ZT Averages (Mid-Phase 3-9 vs 15-21) ---
                zt_day_data = d_last_week.(current_metric_var)(d_last_week.(ztHourVar) >= 3 & d_last_week.(ztHourVar) <= 9);
                zt_night_data = d_last_week.(current_metric_var)(d_last_week.(ztHourVar) >= 15 & d_last_week.(ztHourVar) <= 21);
                zt_d_mean = mean(zt_day_data, 'omitnan');
                zt_n_mean = mean(zt_night_data, 'omitnan');
                
                % --- CT Averages (Mid-Phase 3-9 vs 15-21) ---
                ct_day_data = d_last_week.(current_metric_var)(d_last_week.CT_Hour >= 3 & d_last_week.CT_Hour <= 9);
                ct_night_data = d_last_week.(current_metric_var)(d_last_week.CT_Hour >= 15 & d_last_week.CT_Hour <= 21);
                ct_d_mean = mean(ct_day_data, 'omitnan');
                ct_n_mean = mean(ct_night_data, 'omitnan');
                
                % Store AO5
                if strcmp(string(an_id), "AO5")
                    ao5_stats.zt_day = zt_d_mean; ao5_stats.zt_night = zt_n_mean;
                    ao5_stats.ct_day = ct_d_mean; ao5_stats.ct_night = ct_n_mean;
                end
                
                % Store by Sex
                if ismember(an_id, maleAnimals)
                    m_zt_day(end+1) = zt_d_mean; m_zt_night(end+1) = zt_n_mean;
                    m_ct_day(end+1) = ct_d_mean; m_ct_night(end+1) = ct_n_mean;
                elseif ismember(an_id, femaleAnimals)
                    f_zt_day(end+1) = zt_d_mean; f_zt_night(end+1) = zt_n_mean;
                    f_ct_day(end+1) = ct_d_mean; f_ct_night(end+1) = ct_n_mean;
                end
            end
            
            % --- PRINT DIAGNOSTIC RESULTS ---
            fprintf('\n======================================================\n');
            fprintf('  DIAGNOSTIC: ZT vs CT in Last Week of Full Dark\n');
            fprintf('======================================================\n');
            
            fprintf('\n--- AO5 Specific Check (Tau = 23.81) ---\n');
            if isfield(ao5_stats, 'zt_day')
                fprintf('  Using Old ZT: Day Mean = %.0f, Night Mean = %.0f\n', ao5_stats.zt_day, ao5_stats.zt_night);
                fprintf('  Using New CT: Day Mean = %.0f, Night Mean = %.0f\n', ao5_stats.ct_day, ao5_stats.ct_night);
                pct_diff = abs(ao5_stats.ct_day - ao5_stats.zt_day) / ao5_stats.zt_day * 100;
                fprintf('  Difference in calculated "Day" activity: %.1f%%\n', pct_diff);
            else
                fprintf('  AO5 data not found in this segment.\n');
            end
            
            fprintf('\n--- Population Stats (Paired T-Test: Mid-Day vs Mid-Night) ---\n');
            
            if length(m_zt_day) >= 3
                [~, p_m_zt] = ttest(m_zt_day, m_zt_night);
                [~, p_m_ct] = ttest(m_ct_day, m_ct_night);
                fprintf('  MALES (N=%d):\n', length(m_zt_day));
                fprintf('    Old ZT Method p-value: %.4f\n', p_m_zt);
                fprintf('    New CT Method p-value: %.4f\n', p_m_ct);
            end
            
            if length(f_zt_day) >= 3
                [~, p_f_zt] = ttest(f_zt_day, f_zt_night);
                [~, p_f_ct] = ttest(f_ct_day, f_ct_night);
                fprintf('  FEMALES (N=%d):\n', length(f_zt_day));
                fprintf('    Old ZT Method p-value: %.4f\n', p_f_zt);
                fprintf('    New CT Method p-value: %.4f\n', p_f_ct);
            end
            fprintf('======================================================\n');
            
        catch ME_Diag
            warning('ErrGen:Diag','Error Diagnostic: %s', ME_Diag.message);
        end
    end
end

disp('=====================================================');
disp('--- All Figure Generation Complete ---');
disp('=====================================================');

%% --- 6. HELPER FUNCTIONS ---
function profile_data = get_last_2_weeks_profile(dataTable, animals, cond, ztHourVar, current_metric_var)
    all_animal_profiles = NaN(length(animals), 24);
    for an_idx = 1:length(animals)
        animal_data = dataTable(dataTable.Animal == animals(an_idx) & dataTable.Condition == cond, :);
        if isempty(animal_data), continue; end
        unique_days = unique(animal_data.DayOfCondition);
        if isempty(unique_days), continue; end
        last_14_days = unique_days(max(1, end-13):end);
        final_data = animal_data(ismember(animal_data.DayOfCondition, last_14_days), :);
        if isempty(final_data), continue; end
        hourly_means = groupsummary(final_data, ztHourVar, 'mean', current_metric_var);
        [lia, locb] = ismember(0:23, hourly_means.(ztHourVar));
        all_animal_profiles(an_idx, lia) = hourly_means.(['mean_', current_metric_var])(locb(lia));
    end
    profile_data = all_animal_profiles;
end
function plot_sig_bar(ax, x_coords, p_uncorrected, y_pos, corrected_alpha)
    % Assign default alpha if not provided by the call
    if nargin < 5
        corrected_alpha = 0.05;
    end
    
    % Only proceed to plot the bar if the result is significant after correction
    if p_uncorrected < corrected_alpha
        % Determine the number of asterisks based on the CONVENTIONAL (uncorrected) p-value
        if p_uncorrected < 0.001
            ast_str = '***';
        elseif p_uncorrected < 0.01
            ast_str = '**';
        else
            ast_str = '*';
        end
        
        hold(ax, 'on');
        plot(ax, x_coords, [y_pos, y_pos], '-k', 'LineWidth', 1.2);
        text(ax, mean(x_coords), y_pos, ast_str, ...
            'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 16);
        hold(ax, 'off');
    end
    % If p_uncorrected >= corrected_alpha, the function does nothing, and no bar is plotted.
end