clear all;close all;
all_p_vals=[];
all_tuning_strength=[];

% --- Step 1: Load and Preprocess Data ---
vel_threshold=4;

cd 'D:\ACC Data\ACC\AllDays\635\day9'; 
heading=table2array(readtable('heading.csv'));
head_pos=heading(:,5:6);
body_pos=heading(:,8:9);

% Step 1.1: Calculate the heading angle for the entire session
delta_x = head_pos(:, 1) - body_pos(:, 1);
delta_y = head_pos(:, 2) - body_pos(:, 2);
heading_angle_rad = atan2(delta_x, delta_y);
heading_angle_deg = heading_angle_rad * (180 / pi);
heading_angle = heading_angle_deg;
heading_angle(heading_angle < 0) = heading_angle(heading_angle < 0) + 360;

% Step 1.2: Load spike data and extract relevant neurons
data=table2array(readtable('XYtotalspike2.csv'));
SCid=table2array(readtable('SC.csv'));
SCid=find(SCid(:,2));
spike_angle_ms1=cell(length(SCid),1);
bad_data=find(data(:,2)==1);
x_posn=data(:,2);
x_posn=smooth(x_posn,15);
y_posn=data(:,3);
y_posn=smooth(y_posn,15);
x_posn(bad_data)=[];
y_posn(bad_data)=[];
data(bad_data,:)=[];
vel_cmsec=abs(diff(x_posn))+abs(diff(y_posn))*(1/0.05);
vel_cmsec=smooth(vel_cmsec,30);
still_periods=find(vel_cmsec<vel_threshold);
data=data(1:end,6:end);
data=data(:,SCid);
notactive=find(sum(data)<5);
data(:,notactive)=[];
data(still_periods,:)=[];
heading_angle(still_periods)=[];
for cell_num=1:size(data,2)
    tempcell=data(:,cell_num);
    tempcell=find(tempcell);
    xy=heading_angle(tempcell);
    spike_angle_ms1{cell_num}=xy;
end

% --- Step 2: Directional Tuning Analysis and Shuffling Test ---

% Define plotting parameters
num_bins = 18; % Number of bins for a 360 degree circle (20 degrees per bin)
bin_edges_deg = linspace(0, 360, num_bins + 1);
bin_edges_rad = deg2rad(bin_edges_deg);

num_neurons = size(spike_angle_ms1,1);
tuning_strength = zeros(1, num_neurons);
p_values = zeros(1, num_neurons);
num_shuffles = 500;
num_time_bins_total = length(heading_angle); % Total number of time bins in the filtered session

% Loop through each neuron
for i = 1:num_neurons
    angles_deg = spike_angle_ms1{i};
    
    if isempty(angles_deg)
        disp(['Neuron ' num2str(i) ' has no spikes. Skipping plot and setting MVL to 0.']);
        tuning_strength(i) = 0;
        p_values(i) = 1;
        continue;
    end
    
    angles_rad = deg2rad(angles_deg);
    
    % Create a circular histogram plot for the neuron
    figure;
    p = polarhistogram(angles_rad, bin_edges_rad, 'Normalization', 'count');
    title(['Neuron ' num2str(i) ' - Head Directional Tuning']);
    set(gca, 'ThetaZeroLocation', 'top', 'ThetaDir', 'clockwise');
    
    % Quantify the tuning strength (MVL)
    mean_cos = sum(cos(angles_rad)) / length(angles_rad);
    mean_sin = sum(sin(angles_rad)) / length(angles_rad);
    tuning_strength(i) = sqrt(mean_cos^2 + mean_sin^2);
    
    % --- Step 3: Perform a shuffling test for significance (CORRECTED METHOD) ---
    % We will now compare the observed MVL to a distribution of MVLs from
    % shuffled data where the spike times are randomized.
    
    shuffled_mvls = zeros(1, num_shuffles);
    num_spikes = length(angles_deg);
    
    for k = 1:num_shuffles
        % Correct shuffling: Randomly sample angles from the full session
        % data, with the same number of spikes as the original.
        shuffled_indices = randi(num_time_bins_total, 1, num_spikes);
        shuffled_angles_rad = deg2rad(heading_angle(shuffled_indices));
        
        % Calculate the MVL of the shuffled data
        shuffled_mean_cos = sum(cos(shuffled_angles_rad)) / length(shuffled_angles_rad);
        shuffled_mean_sin = sum(sin(shuffled_angles_rad)) / length(shuffled_angles_rad);
        shuffled_mvls(k) = sqrt(shuffled_mean_cos^2 + shuffled_mean_sin^2);
    end
    
    % Calculate the p-value: the proportion of shuffled MVLs >= the observed MVL
    p_values(i) = sum(shuffled_mvls >= tuning_strength(i)) / num_shuffles;
    
    % Add the MVL and p-value to the plot title
    current_title = get(gca, 'Title');
    new_title = [current_title.String ' (MVL: ' num2str(tuning_strength(i), '%.3f') ', p = ' num2str(p_values(i), '%.4f') ')'];
    title(new_title);
end

disp(' ');
disp('Tuning Strength (Mean Resultant Vector Length) and P-Value for each neuron:');
disp([tuning_strength; p_values]);

[sorted_mvl, sorted_indices] = sort(tuning_strength, 'descend');

disp(' ');
disp('Neurons sorted by tuning strength:');
for j = 1:num_neurons
    disp(['Neuron ' num2str(sorted_indices(j)) ' (MVL: ' num2str(sorted_mvl(j), '%.3f') ', p = ' num2str(p_values(sorted_indices(j)), '%.4f') ')']);
end
all_p_vals=[all_p_vals;p_values'];
all_tuning_strength=[all_tuning_strength;tuning_strength'];



% --- Step 1: Load and Preprocess Data ---
vel_threshold=4;
cd 'D:\ACC Data\ACC\AllDays\681\day9'; 
heading=table2array(readtable('heading.csv'));
head_pos=heading(:,5:6);
body_pos=heading(:,8:9);

% Step 1.1: Calculate the heading angle for the entire session
delta_x = head_pos(:, 1) - body_pos(:, 1);
delta_y = head_pos(:, 2) - body_pos(:, 2);
heading_angle_rad = atan2(delta_x, delta_y);
heading_angle_deg = heading_angle_rad * (180 / pi);
heading_angle = heading_angle_deg;
heading_angle(heading_angle < 0) = heading_angle(heading_angle < 0) + 360;

% Step 1.2: Load spike data and extract relevant neurons
data=table2array(readtable('XYtotalspike2.csv'));
SCid=table2array(readtable('SC.csv'));
SCid=find(SCid(:,2));
spike_angle_ms2=cell(length(SCid),1);
bad_data=find(data(:,2)==1);
x_posn=data(:,2);
x_posn=smooth(x_posn,15);
y_posn=data(:,3);
y_posn=smooth(y_posn,15);
x_posn(bad_data)=[];
y_posn(bad_data)=[];
data(bad_data,:)=[];
vel_cmsec=abs(diff(x_posn))+abs(diff(y_posn))*(1/0.05);
vel_cmsec=smooth(vel_cmsec,30);
still_periods=find(vel_cmsec<vel_threshold);
data=data(1:end,6:end);
data=data(:,SCid);
notactive=find(sum(data)<5);
data(:,notactive)=[];
data(still_periods,:)=[];
heading_angle(still_periods)=[];
for cell_num=1:size(data,2)
    tempcell=data(:,cell_num);
    tempcell=find(tempcell);
    xy=heading_angle(tempcell);
    spike_angle_ms2{cell_num}=xy;
end

% --- Step 2: Directional Tuning Analysis and Shuffling Test ---

% Define plotting parameters
num_bins = 18; % Number of bins for a 360 degree circle (20 degrees per bin)
bin_edges_deg = linspace(0, 360, num_bins + 1);
bin_edges_rad = deg2rad(bin_edges_deg);

num_neurons = size(spike_angle_ms2,1);
tuning_strength = zeros(1, num_neurons);
p_values = zeros(1, num_neurons);
num_shuffles = 500;
num_time_bins_total = length(heading_angle); % Total number of time bins in the filtered session

% Loop through each neuron
for i = 1:num_neurons
    angles_deg = spike_angle_ms2{i};
    
    if isempty(angles_deg)
        disp(['Neuron ' num2str(i) ' has no spikes. Skipping plot and setting MVL to 0.']);
        tuning_strength(i) = 0;
        p_values(i) = 1;
        continue;
    end
    
    angles_rad = deg2rad(angles_deg);
    
    % Create a circular histogram plot for the neuron
    figure;
    p = polarhistogram(angles_rad, bin_edges_rad, 'Normalization', 'count');
    title(['Neuron ' num2str(i) ' - Head Directional Tuning']);
    set(gca, 'ThetaZeroLocation', 'top', 'ThetaDir', 'clockwise');
    
    % Quantify the tuning strength (MVL)
    mean_cos = sum(cos(angles_rad)) / length(angles_rad);
    mean_sin = sum(sin(angles_rad)) / length(angles_rad);
    tuning_strength(i) = sqrt(mean_cos^2 + mean_sin^2);
    
    % --- Step 3: Perform a shuffling test for significance (CORRECTED METHOD) ---
    % We will now compare the observed MVL to a distribution of MVLs from
    % shuffled data where the spike times are randomized.
    
    shuffled_mvls = zeros(1, num_shuffles);
    num_spikes = length(angles_deg);
    
    for k = 1:num_shuffles
        % Correct shuffling: Randomly sample angles from the full session
        % data, with the same number of spikes as the original.
        shuffled_indices = randi(num_time_bins_total, 1, num_spikes);
        shuffled_angles_rad = deg2rad(heading_angle(shuffled_indices));
        
        % Calculate the MVL of the shuffled data
        shuffled_mean_cos = sum(cos(shuffled_angles_rad)) / length(shuffled_angles_rad);
        shuffled_mean_sin = sum(sin(shuffled_angles_rad)) / length(shuffled_angles_rad);
        shuffled_mvls(k) = sqrt(shuffled_mean_cos^2 + shuffled_mean_sin^2);
    end
    
    % Calculate the p-value: the proportion of shuffled MVLs >= the observed MVL
    p_values(i) = sum(shuffled_mvls >= tuning_strength(i)) / num_shuffles;
    
    % Add the MVL and p-value to the plot title
    current_title = get(gca, 'Title');
    new_title = [current_title.String ' (MVL: ' num2str(tuning_strength(i), '%.3f') ', p = ' num2str(p_values(i), '%.4f') ')'];
    title(new_title);
end

disp(' ');
disp('Tuning Strength (Mean Resultant Vector Length) and P-Value for each neuron:');
disp([tuning_strength; p_values]);

[sorted_mvl, sorted_indices] = sort(tuning_strength, 'descend');

disp(' ');
disp('Neurons sorted by tuning strength:');
for j = 1:num_neurons
    disp(['Neuron ' num2str(sorted_indices(j)) ' (MVL: ' num2str(sorted_mvl(j), '%.3f') ', p = ' num2str(p_values(sorted_indices(j)), '%.4f') ')']);
end
all_p_vals=[all_p_vals;p_values'];
all_tuning_strength=[all_tuning_strength;tuning_strength'];



% --- Step 1: Load and Preprocess Data ---
vel_threshold=4;
cd 'D:\ACC Data\ACC\AllDays\691\day9'; 
heading=table2array(readtable('heading.csv'));
head_pos=heading(:,5:6);
body_pos=heading(:,8:9);

% Step 1.1: Calculate the heading angle for the entire session
delta_x = head_pos(:, 1) - body_pos(:, 1);
delta_y = head_pos(:, 2) - body_pos(:, 2);
heading_angle_rad = atan2(delta_x, delta_y);
heading_angle_deg = heading_angle_rad * (180 / pi);
heading_angle = heading_angle_deg;
heading_angle(heading_angle < 0) = heading_angle(heading_angle < 0) + 360;

% Step 1.2: Load spike data and extract relevant neurons
data=table2array(readtable('XYtotalspike2.csv'));
SCid=table2array(readtable('SC.csv'));
SCid=find(SCid(:,2));
spike_angle_ms3=cell(length(SCid),1);
bad_data=find(data(:,2)==1);
x_posn=data(:,2);
x_posn=smooth(x_posn,15);
y_posn=data(:,3);
y_posn=smooth(y_posn,15);
x_posn(bad_data)=[];
y_posn(bad_data)=[];
data(bad_data,:)=[];
vel_cmsec=abs(diff(x_posn))+abs(diff(y_posn))*(1/0.05);
vel_cmsec=smooth(vel_cmsec,30);
still_periods=find(vel_cmsec<vel_threshold);
data=data(1:end,6:end);
data=data(:,SCid);
notactive=find(sum(data)<5);
data(:,notactive)=[];
data(still_periods,:)=[];
heading_angle(still_periods)=[];
for cell_num=1:size(data,2)
    tempcell=data(:,cell_num);
    tempcell=find(tempcell);
    xy=heading_angle(tempcell);
    spike_angle_ms3{cell_num}=xy;
end

% --- Step 2: Directional Tuning Analysis and Shuffling Test ---

% Define plotting parameters
num_bins = 18; % Number of bins for a 360 degree circle (20 degrees per bin)
bin_edges_deg = linspace(0, 360, num_bins + 1);
bin_edges_rad = deg2rad(bin_edges_deg);

num_neurons = size(spike_angle_ms3,1);
tuning_strength = zeros(1, num_neurons);
p_values = zeros(1, num_neurons);
num_shuffles = 500;
num_time_bins_total = length(heading_angle); % Total number of time bins in the filtered session

% Loop through each neuron
for i = 1:num_neurons
    angles_deg = spike_angle_ms3{i};
    
    if isempty(angles_deg)
        disp(['Neuron ' num2str(i) ' has no spikes. Skipping plot and setting MVL to 0.']);
        tuning_strength(i) = 0;
        p_values(i) = 1;
        continue;
    end
    
    angles_rad = deg2rad(angles_deg);
    
    % Create a circular histogram plot for the neuron
    figure;
    p = polarhistogram(angles_rad, bin_edges_rad, 'Normalization', 'count');
    title(['Neuron ' num2str(i) ' - Head Directional Tuning']);
    set(gca, 'ThetaZeroLocation', 'top', 'ThetaDir', 'clockwise');
    
    % Quantify the tuning strength (MVL)
    mean_cos = sum(cos(angles_rad)) / length(angles_rad);
    mean_sin = sum(sin(angles_rad)) / length(angles_rad);
    tuning_strength(i) = sqrt(mean_cos^2 + mean_sin^2);
    
    % --- Step 3: Perform a shuffling test for significance (CORRECTED METHOD) ---
    % We will now compare the observed MVL to a distribution of MVLs from
    % shuffled data where the spike times are randomized.
    
    shuffled_mvls = zeros(1, num_shuffles);
    num_spikes = length(angles_deg);
    
    for k = 1:num_shuffles
        % Correct shuffling: Randomly sample angles from the full session
        % data, with the same number of spikes as the original.
        shuffled_indices = randi(num_time_bins_total, 1, num_spikes);
        shuffled_angles_rad = deg2rad(heading_angle(shuffled_indices));
        
        % Calculate the MVL of the shuffled data
        shuffled_mean_cos = sum(cos(shuffled_angles_rad)) / length(shuffled_angles_rad);
        shuffled_mean_sin = sum(sin(shuffled_angles_rad)) / length(shuffled_angles_rad);
        shuffled_mvls(k) = sqrt(shuffled_mean_cos^2 + shuffled_mean_sin^2);
    end
    
    % Calculate the p-value: the proportion of shuffled MVLs >= the observed MVL
    p_values(i) = sum(shuffled_mvls >= tuning_strength(i)) / num_shuffles;
    
    % Add the MVL and p-value to the plot title
    current_title = get(gca, 'Title');
    new_title = [current_title.String ' (MVL: ' num2str(tuning_strength(i), '%.3f') ', p = ' num2str(p_values(i), '%.4f') ')'];
    title(new_title);
end

disp(' ');
disp('Tuning Strength (Mean Resultant Vector Length) and P-Value for each neuron:');
disp([tuning_strength; p_values]);

[sorted_mvl, sorted_indices] = sort(tuning_strength, 'descend');

disp(' ');
disp('Neurons sorted by tuning strength:');
for j = 1:num_neurons
    disp(['Neuron ' num2str(sorted_indices(j)) ' (MVL: ' num2str(sorted_mvl(j), '%.3f') ', p = ' num2str(p_values(sorted_indices(j)), '%.4f') ')']);
end
all_p_vals=[all_p_vals;p_values'];
all_tuning_strength=[all_tuning_strength;tuning_strength'];



% --- Step 1: Load and Preprocess Data ---
vel_threshold=4;
cd 'D:\ACC Data\ACC\AllDays\747\day9'; 
heading=table2array(readtable('heading.csv'));
head_pos=heading(:,5:6);
body_pos=heading(:,8:9);

% Step 1.1: Calculate the heading angle for the entire session
delta_x = head_pos(:, 1) - body_pos(:, 1);
delta_y = head_pos(:, 2) - body_pos(:, 2);
heading_angle_rad = atan2(delta_x, delta_y);
heading_angle_deg = heading_angle_rad * (180 / pi);
heading_angle = heading_angle_deg;
heading_angle(heading_angle < 0) = heading_angle(heading_angle < 0) + 360;

% Step 1.2: Load spike data and extract relevant neurons
data=table2array(readtable('XYtotalspike2.csv'));
SCid=table2array(readtable('SC.csv'));
SCid=find(SCid(:,2));
spike_angle_ms4=cell(length(SCid),1);
bad_data=find(data(:,2)==1);
x_posn=data(:,2);
x_posn=smooth(x_posn,15);
y_posn=data(:,3);
y_posn=smooth(y_posn,15);
x_posn(bad_data)=[];
y_posn(bad_data)=[];
data(bad_data,:)=[];
vel_cmsec=abs(diff(x_posn))+abs(diff(y_posn))*(1/0.05);
vel_cmsec=smooth(vel_cmsec,30);
still_periods=find(vel_cmsec<vel_threshold);
data=data(1:end,6:end);
data=data(:,SCid);
notactive=find(sum(data)<5);
data(:,notactive)=[];
data(still_periods,:)=[];
heading_angle(still_periods)=[];
for cell_num=1:size(data,2)
    tempcell=data(:,cell_num);
    tempcell=find(tempcell);
    xy=heading_angle(tempcell);
    spike_angle_ms4{cell_num}=xy;
end

% --- Step 2: Directional Tuning Analysis and Shuffling Test ---

% Define plotting parameters
num_bins = 18; % Number of bins for a 360 degree circle (20 degrees per bin)
bin_edges_deg = linspace(0, 360, num_bins + 1);
bin_edges_rad = deg2rad(bin_edges_deg);

num_neurons = size(spike_angle_ms4,1);
tuning_strength = zeros(1, num_neurons);
p_values = zeros(1, num_neurons);
num_shuffles = 500;
num_time_bins_total = length(heading_angle); % Total number of time bins in the filtered session

% Loop through each neuron
for i = 1:num_neurons
    angles_deg = spike_angle_ms4{i};
    
    if isempty(angles_deg)
        disp(['Neuron ' num2str(i) ' has no spikes. Skipping plot and setting MVL to 0.']);
        tuning_strength(i) = 0;
        p_values(i) = 1;
        continue;
    end
    
    angles_rad = deg2rad(angles_deg);
    
    % Create a circular histogram plot for the neuron
    figure;
    p = polarhistogram(angles_rad, bin_edges_rad, 'Normalization', 'count');
    title(['Neuron ' num2str(i) ' - Head Directional Tuning']);
    set(gca, 'ThetaZeroLocation', 'top', 'ThetaDir', 'clockwise');
    
    % Quantify the tuning strength (MVL)
    mean_cos = sum(cos(angles_rad)) / length(angles_rad);
    mean_sin = sum(sin(angles_rad)) / length(angles_rad);
    tuning_strength(i) = sqrt(mean_cos^2 + mean_sin^2);
    
    % --- Step 3: Perform a shuffling test for significance (CORRECTED METHOD) ---
    % We will now compare the observed MVL to a distribution of MVLs from
    % shuffled data where the spike times are randomized.
    
    shuffled_mvls = zeros(1, num_shuffles);
    num_spikes = length(angles_deg);
    
    for k = 1:num_shuffles
        % Correct shuffling: Randomly sample angles from the full session
        % data, with the same number of spikes as the original.
        shuffled_indices = randi(num_time_bins_total, 1, num_spikes);
        shuffled_angles_rad = deg2rad(heading_angle(shuffled_indices));
        
        % Calculate the MVL of the shuffled data
        shuffled_mean_cos = sum(cos(shuffled_angles_rad)) / length(shuffled_angles_rad);
        shuffled_mean_sin = sum(sin(shuffled_angles_rad)) / length(shuffled_angles_rad);
        shuffled_mvls(k) = sqrt(shuffled_mean_cos^2 + shuffled_mean_sin^2);
    end
    
    % Calculate the p-value: the proportion of shuffled MVLs >= the observed MVL
    p_values(i) = sum(shuffled_mvls >= tuning_strength(i)) / num_shuffles;
    
    % Add the MVL and p-value to the plot title
    current_title = get(gca, 'Title');
    new_title = [current_title.String ' (MVL: ' num2str(tuning_strength(i), '%.3f') ', p = ' num2str(p_values(i), '%.4f') ')'];
    title(new_title);
end

disp(' ');
disp('Tuning Strength (Mean Resultant Vector Length) and P-Value for each neuron:');
disp([tuning_strength; p_values]);

[sorted_mvl, sorted_indices] = sort(tuning_strength, 'descend');

disp(' ');
disp('Neurons sorted by tuning strength:');
for j = 1:num_neurons
    disp(['Neuron ' num2str(sorted_indices(j)) ' (MVL: ' num2str(sorted_mvl(j), '%.3f') ', p = ' num2str(p_values(sorted_indices(j)), '%.4f') ')']);
end
all_p_vals=[all_p_vals;p_values'];
all_tuning_strength=[all_tuning_strength;tuning_strength'];


all_p_valsSC=all_p_vals;
all_tuning_strengthSC=all_tuning_strength;
all_p_vals=[];
all_tuning_strength=[];

% --- Step 1: Load and Preprocess Data ---
vel_threshold=4;

cd 'D:\ACC Data\ACC\AllDays\635\day9'; 
heading=table2array(readtable('heading.csv'));
head_pos=heading(:,5:6);
body_pos=heading(:,8:9);

% Step 1.1: Calculate the heading angle for the entire session
delta_x = head_pos(:, 1) - body_pos(:, 1);
delta_y = head_pos(:, 2) - body_pos(:, 2);
heading_angle_rad = atan2(delta_x, delta_y);
heading_angle_deg = heading_angle_rad * (180 / pi);
heading_angle = heading_angle_deg;
heading_angle(heading_angle < 0) = heading_angle(heading_angle < 0) + 360;

% Step 1.2: Load spike data and extract relevant neurons
data=table2array(readtable('XYtotalspike2.csv'));
SCid=table2array(readtable('pPC.csv'));
SCid=find(SCid(:,2));
spike_angle_ms1=cell(length(SCid),1);
bad_data=find(data(:,2)==1);
x_posn=data(:,2);
x_posn=smooth(x_posn,15);
y_posn=data(:,3);
y_posn=smooth(y_posn,15);
x_posn(bad_data)=[];
y_posn(bad_data)=[];
data(bad_data,:)=[];
vel_cmsec=abs(diff(x_posn))+abs(diff(y_posn))*(1/0.05);
vel_cmsec=smooth(vel_cmsec,30);
still_periods=find(vel_cmsec<vel_threshold);
data=data(1:end,6:end);
data=data(:,SCid);
notactive=find(sum(data)<5);
data(:,notactive)=[];
data(still_periods,:)=[];
heading_angle(still_periods)=[];
for cell_num=1:size(data,2)
    tempcell=data(:,cell_num);
    tempcell=find(tempcell);
    xy=heading_angle(tempcell);
    spike_angle_ms1{cell_num}=xy;
end

% --- Step 2: Directional Tuning Analysis and Shuffling Test ---

% Define plotting parameters
num_bins = 18; % Number of bins for a 360 degree circle (20 degrees per bin)
bin_edges_deg = linspace(0, 360, num_bins + 1);
bin_edges_rad = deg2rad(bin_edges_deg);

num_neurons = size(spike_angle_ms1,1);
tuning_strength = zeros(1, num_neurons);
p_values = zeros(1, num_neurons);
num_shuffles = 500;
num_time_bins_total = length(heading_angle); % Total number of time bins in the filtered session

% Loop through each neuron
for i = 1:num_neurons
    angles_deg = spike_angle_ms1{i};
    
    if isempty(angles_deg)
        disp(['Neuron ' num2str(i) ' has no spikes. Skipping plot and setting MVL to 0.']);
        tuning_strength(i) = 0;
        p_values(i) = 1;
        continue;
    end
    
    angles_rad = deg2rad(angles_deg);
    
    % Create a circular histogram plot for the neuron
    figure;
    p = polarhistogram(angles_rad, bin_edges_rad, 'Normalization', 'count');
    title(['Neuron ' num2str(i) ' - Head Directional Tuning']);
    set(gca, 'ThetaZeroLocation', 'top', 'ThetaDir', 'clockwise');
    
    % Quantify the tuning strength (MVL)
    mean_cos = sum(cos(angles_rad)) / length(angles_rad);
    mean_sin = sum(sin(angles_rad)) / length(angles_rad);
    tuning_strength(i) = sqrt(mean_cos^2 + mean_sin^2);
    
    % --- Step 3: Perform a shuffling test for significance (CORRECTED METHOD) ---
    % We will now compare the observed MVL to a distribution of MVLs from
    % shuffled data where the spike times are randomized.
    
    shuffled_mvls = zeros(1, num_shuffles);
    num_spikes = length(angles_deg);
    
    for k = 1:num_shuffles
        % Correct shuffling: Randomly sample angles from the full session
        % data, with the same number of spikes as the original.
        shuffled_indices = randi(num_time_bins_total, 1, num_spikes);
        shuffled_angles_rad = deg2rad(heading_angle(shuffled_indices));
        
        % Calculate the MVL of the shuffled data
        shuffled_mean_cos = sum(cos(shuffled_angles_rad)) / length(shuffled_angles_rad);
        shuffled_mean_sin = sum(sin(shuffled_angles_rad)) / length(shuffled_angles_rad);
        shuffled_mvls(k) = sqrt(shuffled_mean_cos^2 + shuffled_mean_sin^2);
    end
    
    % Calculate the p-value: the proportion of shuffled MVLs >= the observed MVL
    p_values(i) = sum(shuffled_mvls >= tuning_strength(i)) / num_shuffles;
    
    % Add the MVL and p-value to the plot title
    current_title = get(gca, 'Title');
    new_title = [current_title.String ' (MVL: ' num2str(tuning_strength(i), '%.3f') ', p = ' num2str(p_values(i), '%.4f') ')'];
    title(new_title);
end

disp(' ');
disp('Tuning Strength (Mean Resultant Vector Length) and P-Value for each neuron:');
disp([tuning_strength; p_values]);

[sorted_mvl, sorted_indices] = sort(tuning_strength, 'descend');

disp(' ');
disp('Neurons sorted by tuning strength:');
for j = 1:num_neurons
    disp(['Neuron ' num2str(sorted_indices(j)) ' (MVL: ' num2str(sorted_mvl(j), '%.3f') ', p = ' num2str(p_values(sorted_indices(j)), '%.4f') ')']);
end
all_p_vals=[all_p_vals;p_values'];
all_tuning_strength=[all_tuning_strength;tuning_strength'];



% --- Step 1: Load and Preprocess Data ---
vel_threshold=4;
cd 'D:\ACC Data\ACC\AllDays\681\day9'; 
heading=table2array(readtable('heading.csv'));
head_pos=heading(:,5:6);
body_pos=heading(:,8:9);

% Step 1.1: Calculate the heading angle for the entire session
delta_x = head_pos(:, 1) - body_pos(:, 1);
delta_y = head_pos(:, 2) - body_pos(:, 2);
heading_angle_rad = atan2(delta_x, delta_y);
heading_angle_deg = heading_angle_rad * (180 / pi);
heading_angle = heading_angle_deg;
heading_angle(heading_angle < 0) = heading_angle(heading_angle < 0) + 360;

% Step 1.2: Load spike data and extract relevant neurons
data=table2array(readtable('XYtotalspike2.csv'));
SCid=table2array(readtable('pPC.csv'));
SCid=find(SCid(:,2));
spike_angle_ms2=cell(length(SCid),1);
bad_data=find(data(:,2)==1);
x_posn=data(:,2);
x_posn=smooth(x_posn,15);
y_posn=data(:,3);
y_posn=smooth(y_posn,15);
x_posn(bad_data)=[];
y_posn(bad_data)=[];
data(bad_data,:)=[];
vel_cmsec=abs(diff(x_posn))+abs(diff(y_posn))*(1/0.05);
vel_cmsec=smooth(vel_cmsec,30);
still_periods=find(vel_cmsec<vel_threshold);
data=data(1:end,6:end);
data=data(:,SCid);
notactive=find(sum(data)<5);
data(:,notactive)=[];
data(still_periods,:)=[];
heading_angle(still_periods)=[];
for cell_num=1:size(data,2)
    tempcell=data(:,cell_num);
    tempcell=find(tempcell);
    xy=heading_angle(tempcell);
    spike_angle_ms2{cell_num}=xy;
end

% --- Step 2: Directional Tuning Analysis and Shuffling Test ---

% Define plotting parameters
num_bins = 18; % Number of bins for a 360 degree circle (20 degrees per bin)
bin_edges_deg = linspace(0, 360, num_bins + 1);
bin_edges_rad = deg2rad(bin_edges_deg);

num_neurons = size(spike_angle_ms2,1);
tuning_strength = zeros(1, num_neurons);
p_values = zeros(1, num_neurons);
num_shuffles = 500;
num_time_bins_total = length(heading_angle); % Total number of time bins in the filtered session

% Loop through each neuron
for i = 1:num_neurons
    angles_deg = spike_angle_ms2{i};
    
    if isempty(angles_deg)
        disp(['Neuron ' num2str(i) ' has no spikes. Skipping plot and setting MVL to 0.']);
        tuning_strength(i) = 0;
        p_values(i) = 1;
        continue;
    end
    
    angles_rad = deg2rad(angles_deg);
    
    % Create a circular histogram plot for the neuron
    figure;
    p = polarhistogram(angles_rad, bin_edges_rad, 'Normalization', 'count');
    title(['Neuron ' num2str(i) ' - Head Directional Tuning']);
    set(gca, 'ThetaZeroLocation', 'top', 'ThetaDir', 'clockwise');
    
    % Quantify the tuning strength (MVL)
    mean_cos = sum(cos(angles_rad)) / length(angles_rad);
    mean_sin = sum(sin(angles_rad)) / length(angles_rad);
    tuning_strength(i) = sqrt(mean_cos^2 + mean_sin^2);
    
    % --- Step 3: Perform a shuffling test for significance (CORRECTED METHOD) ---
    % We will now compare the observed MVL to a distribution of MVLs from
    % shuffled data where the spike times are randomized.
    
    shuffled_mvls = zeros(1, num_shuffles);
    num_spikes = length(angles_deg);
    
    for k = 1:num_shuffles
        % Correct shuffling: Randomly sample angles from the full session
        % data, with the same number of spikes as the original.
        shuffled_indices = randi(num_time_bins_total, 1, num_spikes);
        shuffled_angles_rad = deg2rad(heading_angle(shuffled_indices));
        
        % Calculate the MVL of the shuffled data
        shuffled_mean_cos = sum(cos(shuffled_angles_rad)) / length(shuffled_angles_rad);
        shuffled_mean_sin = sum(sin(shuffled_angles_rad)) / length(shuffled_angles_rad);
        shuffled_mvls(k) = sqrt(shuffled_mean_cos^2 + shuffled_mean_sin^2);
    end
    
    % Calculate the p-value: the proportion of shuffled MVLs >= the observed MVL
    p_values(i) = sum(shuffled_mvls >= tuning_strength(i)) / num_shuffles;
    
    % Add the MVL and p-value to the plot title
    current_title = get(gca, 'Title');
    new_title = [current_title.String ' (MVL: ' num2str(tuning_strength(i), '%.3f') ', p = ' num2str(p_values(i), '%.4f') ')'];
    title(new_title);
end

disp(' ');
disp('Tuning Strength (Mean Resultant Vector Length) and P-Value for each neuron:');
disp([tuning_strength; p_values]);

[sorted_mvl, sorted_indices] = sort(tuning_strength, 'descend');

disp(' ');
disp('Neurons sorted by tuning strength:');
for j = 1:num_neurons
    disp(['Neuron ' num2str(sorted_indices(j)) ' (MVL: ' num2str(sorted_mvl(j), '%.3f') ', p = ' num2str(p_values(sorted_indices(j)), '%.4f') ')']);
end
all_p_vals=[all_p_vals;p_values'];
all_tuning_strength=[all_tuning_strength;tuning_strength'];



% --- Step 1: Load and Preprocess Data ---
vel_threshold=4;
cd 'D:\ACC Data\ACC\AllDays\691\day9'; 
heading=table2array(readtable('heading.csv'));
head_pos=heading(:,5:6);
body_pos=heading(:,8:9);

% Step 1.1: Calculate the heading angle for the entire session
delta_x = head_pos(:, 1) - body_pos(:, 1);
delta_y = head_pos(:, 2) - body_pos(:, 2);
heading_angle_rad = atan2(delta_x, delta_y);
heading_angle_deg = heading_angle_rad * (180 / pi);
heading_angle = heading_angle_deg;
heading_angle(heading_angle < 0) = heading_angle(heading_angle < 0) + 360;

% Step 1.2: Load spike data and extract relevant neurons
data=table2array(readtable('XYtotalspike2.csv'));
SCid=table2array(readtable('pPC.csv'));
SCid=find(SCid(:,2));
spike_angle_ms3=cell(length(SCid),1);
bad_data=find(data(:,2)==1);
x_posn=data(:,2);
x_posn=smooth(x_posn,15);
y_posn=data(:,3);
y_posn=smooth(y_posn,15);
x_posn(bad_data)=[];
y_posn(bad_data)=[];
data(bad_data,:)=[];
vel_cmsec=abs(diff(x_posn))+abs(diff(y_posn))*(1/0.05);
vel_cmsec=smooth(vel_cmsec,30);
still_periods=find(vel_cmsec<vel_threshold);
data=data(1:end,6:end);
data=data(:,SCid);
notactive=find(sum(data)<5);
data(:,notactive)=[];
data(still_periods,:)=[];
heading_angle(still_periods)=[];
for cell_num=1:size(data,2)
    tempcell=data(:,cell_num);
    tempcell=find(tempcell);
    xy=heading_angle(tempcell);
    spike_angle_ms3{cell_num}=xy;
end

% --- Step 2: Directional Tuning Analysis and Shuffling Test ---

% Define plotting parameters
num_bins = 18; % Number of bins for a 360 degree circle (20 degrees per bin)
bin_edges_deg = linspace(0, 360, num_bins + 1);
bin_edges_rad = deg2rad(bin_edges_deg);

num_neurons = size(spike_angle_ms3,1);
tuning_strength = zeros(1, num_neurons);
p_values = zeros(1, num_neurons);
num_shuffles = 500;
num_time_bins_total = length(heading_angle); % Total number of time bins in the filtered session

% Loop through each neuron
for i = 1:num_neurons
    angles_deg = spike_angle_ms3{i};
    
    if isempty(angles_deg)
        disp(['Neuron ' num2str(i) ' has no spikes. Skipping plot and setting MVL to 0.']);
        tuning_strength(i) = 0;
        p_values(i) = 1;
        continue;
    end
    
    angles_rad = deg2rad(angles_deg);
    
    % Create a circular histogram plot for the neuron
    figure;
    p = polarhistogram(angles_rad, bin_edges_rad, 'Normalization', 'count');
    title(['Neuron ' num2str(i) ' - Head Directional Tuning']);
    set(gca, 'ThetaZeroLocation', 'top', 'ThetaDir', 'clockwise');
    
    % Quantify the tuning strength (MVL)
    mean_cos = sum(cos(angles_rad)) / length(angles_rad);
    mean_sin = sum(sin(angles_rad)) / length(angles_rad);
    tuning_strength(i) = sqrt(mean_cos^2 + mean_sin^2);
    
    % --- Step 3: Perform a shuffling test for significance (CORRECTED METHOD) ---
    % We will now compare the observed MVL to a distribution of MVLs from
    % shuffled data where the spike times are randomized.
    
    shuffled_mvls = zeros(1, num_shuffles);
    num_spikes = length(angles_deg);
    
    for k = 1:num_shuffles
        % Correct shuffling: Randomly sample angles from the full session
        % data, with the same number of spikes as the original.
        shuffled_indices = randi(num_time_bins_total, 1, num_spikes);
        shuffled_angles_rad = deg2rad(heading_angle(shuffled_indices));
        
        % Calculate the MVL of the shuffled data
        shuffled_mean_cos = sum(cos(shuffled_angles_rad)) / length(shuffled_angles_rad);
        shuffled_mean_sin = sum(sin(shuffled_angles_rad)) / length(shuffled_angles_rad);
        shuffled_mvls(k) = sqrt(shuffled_mean_cos^2 + shuffled_mean_sin^2);
    end
    
    % Calculate the p-value: the proportion of shuffled MVLs >= the observed MVL
    p_values(i) = sum(shuffled_mvls >= tuning_strength(i)) / num_shuffles;
    
    % Add the MVL and p-value to the plot title
    current_title = get(gca, 'Title');
    new_title = [current_title.String ' (MVL: ' num2str(tuning_strength(i), '%.3f') ', p = ' num2str(p_values(i), '%.4f') ')'];
    title(new_title);
end

disp(' ');
disp('Tuning Strength (Mean Resultant Vector Length) and P-Value for each neuron:');
disp([tuning_strength; p_values]);

[sorted_mvl, sorted_indices] = sort(tuning_strength, 'descend');

disp(' ');
disp('Neurons sorted by tuning strength:');
for j = 1:num_neurons
    disp(['Neuron ' num2str(sorted_indices(j)) ' (MVL: ' num2str(sorted_mvl(j), '%.3f') ', p = ' num2str(p_values(sorted_indices(j)), '%.4f') ')']);
end
all_p_vals=[all_p_vals;p_values'];
all_tuning_strength=[all_tuning_strength;tuning_strength'];



% --- Step 1: Load and Preprocess Data ---
vel_threshold=4;
cd 'D:\ACC Data\ACC\AllDays\747\day9'; 
heading=table2array(readtable('heading.csv'));
head_pos=heading(:,5:6);
body_pos=heading(:,8:9);

% Step 1.1: Calculate the heading angle for the entire session
delta_x = head_pos(:, 1) - body_pos(:, 1);
delta_y = head_pos(:, 2) - body_pos(:, 2);
heading_angle_rad = atan2(delta_x, delta_y);
heading_angle_deg = heading_angle_rad * (180 / pi);
heading_angle = heading_angle_deg;
heading_angle(heading_angle < 0) = heading_angle(heading_angle < 0) + 360;

% Step 1.2: Load spike data and extract relevant neurons
data=table2array(readtable('XYtotalspike2.csv'));
SCid=table2array(readtable('pPC.csv'));
SCid=find(SCid(:,2));
spike_angle_ms4=cell(length(SCid),1);
bad_data=find(data(:,2)==1);
x_posn=data(:,2);
x_posn=smooth(x_posn,15);
y_posn=data(:,3);
y_posn=smooth(y_posn,15);
x_posn(bad_data)=[];
y_posn(bad_data)=[];
data(bad_data,:)=[];
vel_cmsec=abs(diff(x_posn))+abs(diff(y_posn))*(1/0.05);
vel_cmsec=smooth(vel_cmsec,30);
still_periods=find(vel_cmsec<vel_threshold);
data=data(1:end,6:end);
data=data(:,SCid);
notactive=find(sum(data)<5);
data(:,notactive)=[];
data(still_periods,:)=[];
heading_angle(still_periods)=[];
for cell_num=1:size(data,2)
    tempcell=data(:,cell_num);
    tempcell=find(tempcell);
    xy=heading_angle(tempcell);
    spike_angle_ms4{cell_num}=xy;
end

% --- Step 2: Directional Tuning Analysis and Shuffling Test ---

% Define plotting parameters
num_bins = 18; % Number of bins for a 360 degree circle (20 degrees per bin)
bin_edges_deg = linspace(0, 360, num_bins + 1);
bin_edges_rad = deg2rad(bin_edges_deg);

num_neurons = size(spike_angle_ms4,1);
tuning_strength = zeros(1, num_neurons);
p_values = zeros(1, num_neurons);
num_shuffles = 500;
num_time_bins_total = length(heading_angle); % Total number of time bins in the filtered session

% Loop through each neuron
for i = 1:num_neurons
    angles_deg = spike_angle_ms4{i};
    
    if isempty(angles_deg)
        disp(['Neuron ' num2str(i) ' has no spikes. Skipping plot and setting MVL to 0.']);
        tuning_strength(i) = 0;
        p_values(i) = 1;
        continue;
    end
    
    angles_rad = deg2rad(angles_deg);
    
    % Create a circular histogram plot for the neuron
    figure;
    p = polarhistogram(angles_rad, bin_edges_rad, 'Normalization', 'count');
    title(['Neuron ' num2str(i) ' - Head Directional Tuning']);
    set(gca, 'ThetaZeroLocation', 'top', 'ThetaDir', 'clockwise');
    
    % Quantify the tuning strength (MVL)
    mean_cos = sum(cos(angles_rad)) / length(angles_rad);
    mean_sin = sum(sin(angles_rad)) / length(angles_rad);
    tuning_strength(i) = sqrt(mean_cos^2 + mean_sin^2);
    
    % --- Step 3: Perform a shuffling test for significance (CORRECTED METHOD) ---
    % We will now compare the observed MVL to a distribution of MVLs from
    % shuffled data where the spike times are randomized.
    
    shuffled_mvls = zeros(1, num_shuffles);
    num_spikes = length(angles_deg);
    
    for k = 1:num_shuffles
        % Correct shuffling: Randomly sample angles from the full session
        % data, with the same number of spikes as the original.
        shuffled_indices = randi(num_time_bins_total, 1, num_spikes);
        shuffled_angles_rad = deg2rad(heading_angle(shuffled_indices));
        
        % Calculate the MVL of the shuffled data
        shuffled_mean_cos = sum(cos(shuffled_angles_rad)) / length(shuffled_angles_rad);
        shuffled_mean_sin = sum(sin(shuffled_angles_rad)) / length(shuffled_angles_rad);
        shuffled_mvls(k) = sqrt(shuffled_mean_cos^2 + shuffled_mean_sin^2);
    end
    
    % Calculate the p-value: the proportion of shuffled MVLs >= the observed MVL
    p_values(i) = sum(shuffled_mvls >= tuning_strength(i)) / num_shuffles;
    
    % Add the MVL and p-value to the plot title
    current_title = get(gca, 'Title');
    new_title = [current_title.String ' (MVL: ' num2str(tuning_strength(i), '%.3f') ', p = ' num2str(p_values(i), '%.4f') ')'];
    title(new_title);
end

disp(' ');
disp('Tuning Strength (Mean Resultant Vector Length) and P-Value for each neuron:');
disp([tuning_strength; p_values]);

[sorted_mvl, sorted_indices] = sort(tuning_strength, 'descend');

disp(' ');
disp('Neurons sorted by tuning strength:');
for j = 1:num_neurons
    disp(['Neuron ' num2str(sorted_indices(j)) ' (MVL: ' num2str(sorted_mvl(j), '%.3f') ', p = ' num2str(p_values(sorted_indices(j)), '%.4f') ')']);
end
all_p_vals=[all_p_vals;p_values'];
all_tuning_strength=[all_tuning_strength;tuning_strength'];