clc 
clear all
close all
%%
%inputs
%CHANGES TO MAKE-values of field length are in voxels ,have to convert

vl=0.25;%voxel length

%firr_thresh=0.2;%0.2

tilted=1;%check if tilted

%calculate stats from which neuron to which neuron?

if tilted
    %load("Tilted_data_props_b0.9b_std2A.mat")
    load("Tilted_data_props_b0.8pi_std3d1.5_std2d1.5.mat")
    %load("field_elongstats_tilted.mat")
    %neuron_end = 50;
    neuron_start = 1;
    neurons =   [5,6,7,10,14,18,25,26,33,36,39,41,46];
    %neurons = [2  6  7 10 11 12 19 21 24 25 27 30 32 33 37 41 45 47 48 50]; % with 0.9 std threshold
    %neurons = [0, 1, 2, 3, 5, 6, 7, 8, 10, 11, 14, 19, 20, 23, 24, 25, 26, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 40, 41, 42, 43, 45, 46, 47, 48, 49]+1;  %With aligned phi case
    %neurons = [ 1  3  4  6  9 10 11 13 14 15 17 21 26 27 32 33 34 35 37 39 41 42 43 44 46 48 49];

else
    load("Aligned_data_props_b2.51_t1.7_2Dt1.5.mat")
    neuron_start = 1;
    %load("field_elongstats_aligned.mat")
    %neuron_end = 47;
    neurons = [2, 4, 6, 8, 11, 16, 18, 24, 26, 30, 33, 34, 41] + 1;
end

minv=50;%minimum voxels to be considered as place field

%CHOOSE PLOTS TO PLOT
identify_placef=1;%set this to 1 to identify place fields ,check min voxels
calc_orientation_elongation=1;
plot_orientation=1; %plot orientation on sphere
plot_elongation=1;
plot_centroids=1;
plot_binary=1;%binary morphology
fields_per_cell=0;%plot fields per cell
plot_pf=0;%plot place field firing rate map
calc_proportion=0; %calculate fields parallel to x,y,z
plot_flength=0;
plotfirr=0;
neuron_firr=1;%neuron for which to plot firr
rotm = [[ 0.57357644,  0.,          0.81915204],
        [ 0.57922797,  0.70710678, -0.40557979],
        [-0.57922797,  0.70710678,  0.40557979]];
% rotm = [[ 0.70710678, -0.57922797, -0.40557979],
%        [ 0.        ,  0.57357644, -0.81915204],
%        [ 0.70710678,  0.57922797,  0.40557979]];
XYZ = [1 0 0; 0 1 0; 0 0 1];
ABC = rotm*XYZ;
axes_cords = [XYZ; ABC];
% set(gca,'FontSize',30);

%% %%%%%%% IDENTIFY PLACE CELLS %%%%%%%%%%%%

x1=1.125:0.25:5.875;
y1=1.125:0.25:5.875;
z1=1.125:0.25:5.875;

[X,Y,Z]=meshgrid(x1, y1 ,z1);
voxel_coords=[Y(:) X(:) Z(:)];%voxel coords increase along x then y then z

orientation=double.empty;
orientationplot=double.empty;

elongation=double.empty;

yfl=double.empty;
xfl=double.empty;
zfl=double.empty;
centroids=double.empty;

%%identify place fields from all neurons and store in pf structure
if identify_placef
    pf=struct;
    pf_counter=1;
%     for neuron=setdiff(neuron_start:neuron_end,avoid_neurons)
    for neuron = neurons
        propstable=FIRR_3dprops(neuron).props; %get 3d properties table of the neuron
        firrtable=FIRR_3dprops(neuron).firr;
        for i=1:height(propstable) 
            
            vol=propstable{i,"Volume"};
            
            if  vol>minv  %if volume of cluster is greater than minvoxels
                pf(pf_counter).props=propstable(i,:);
                pf(pf_counter).neuron=neuron;
                pf(pf_counter).firr=firrtable;
                pf_counter=pf_counter+1;
                
            end
        end
    end
end

%% %%%%%%% PLOT CENTROIDS FIG 4 %%%%%%%%%%%%


if plot_centroids

    %%%%%%%%%%%%% FIG 4 a1) Aligned or FIG 4 b1) Tilted %%%%%%%%%%%%%%
    set(0, 'DefaultFigureRenderer', 'painters');
    for i=1:size(pf,2)
        
        centroid=pf(i).props.Centroid;
        centroids=[centroids;centroid];
        centroids=centroids;
        
    end
%     centroids = rotm * centroids';
    centroids_scaled=(5/20)*(centroids)+1;
    if tilted == 1
    centroids_scaled1 = rotm * centroids_scaled';
    centroids_scaled1 = centroids_scaled1';
%     x_centroid = [10 0 0];
%     x_centroid_rot = x_centroid * rotm;
%     y_centroid = [0 10 0];
%     y_centroid_rot = y_centroid * rotm;
%     z_centroid = [0 0 10];
%     z_centroid_rot = z_centroid * rotm;
    end
    
    %%hold off
    
    %figure
    if tilted == 1
             V1 = [1 1 1];
             V2 = [6 1 1];
             V3 = [6 1 6];
             V4 = [1 1 6];
             V5 = [1 6 6];
             V6 = [1 6 1];
             V7 = [1 6 6];
             V8 = [6 6 6];
             V9 = [6 6 1];
             V10 = [1 6 1];
             V11 = [1 1 1];
             V12 = [1 1 6];
             V13 = [6 1 6];
             V14 = [6 6 6];
             V15 = [6 6 1];
             V16 = [6 1 1];
             C = [V1;V2;V3;V4;V5;V6;V7;V8;V9;V10;V11;V12;V13;V14;V15;V16];
             E = rotm*C';
             E = E';
             figure
             plot3(E(:,1),E(:,2),E(:,3),'black')
             hold on
    end
    if tilted
    scatter3(centroids_scaled1(:,1) ,centroids_scaled1(:,2),  centroids_scaled1(:,3),30, centroids_scaled1(:,3),"filled")
    title("Centroids of place fields- Tilted lattice")
    centroid = mean(E);

    % Define line directions in ABC axes
    line_directions = [ABC(:, 1)'; ABC(:, 2)'; ABC(:, 3)'];

    % Define line length
    line_length = 5;

    % Calculate line end points
    line_ends = centroid + line_length * line_directions;
    
    % Plot the lines
    hold on;
    line_colors = {'red', 'green', 'blue'};
    for i = 1:3
        line_coords = [centroid; line_ends(i, :)];
        plot3(line_coords(:, 1), line_coords(:, 2), line_coords(:, 3),  'Color', line_colors{i},'LineWidth', 2 );
    end
%     hold on
%     plot3(x_centroid_rot, y_centroid_rot, z_centroid_rot);
    daspect([1 1 1]);

        %hold on
    else
    scatter3(centroids_scaled(:,1) ,centroids_scaled(:,2),  centroids_scaled(:,3),30, centroids_scaled(:,3),"filled")
    title("Centroids of place fields- Aligned lattice")
    %hold on
        
    end
    xlabel('x','Fontsize',15)
    ylabel('y','Fontsize',15)
    zlabel('z','Fontsize',15)
    
    if tilted  
    medx=mean(centroids_scaled1(:,1));
    medy=mean(centroids_scaled1(:,2));
    medz=mean(centroids_scaled1(:,3));
    else
    medx=mean(centroids_scaled(:,1));
    
    medy=mean(centroids_scaled(:,2));
    medz=mean(centroids_scaled(:,3));
    end
    if tilted == 0
        hold on
%
      scatter3(medx,medy,medz,100,'filled','red')
      hold on
    end
    
    if tilted == 1
    hold on
      scatter3(medx,medy,medz,100,'filled','red')% scatter3(medx,3.5,3.5,100,'filled','green')

    end

    x2=3.5*ones(length(x1),1);% x2=3.5*ones(length(x1),1);
    y2=3.5*ones(length(y1),1);
    z2=3.5*ones(length(z1),1);
    xyz_rot = rotm * [x2 y2 z2]';

    if tilted ==1
%         plot3(x1, xyz_rot(2, :), xyz_rot(3, :), 'Color', 'red');
%         hold on
%         plot3(xyz_rot(1, :), y1, xyz_rot(3, :), 'Color', 'blue');
%         hold on
%         plot3(xyz_rot(1, :), xyz_rot(2, :), z1, 'Color', 'green');
%         hold on
%         plot3(x1,y2,z2,'Color','black');
%         hold on
%         plot3(x2,y1,z2,'Color','black');
%         hold on
%         plot3(x2,y2,z1,'Color','black');
%         hold on
    end
    
    hold off

    %%%%%%%%%%% Fig 4a2) Aligned or Fig 4b2) Tilted  %%%%%%%%%%%
    
    if tilted
        cc = inv(rotm)*centroids_scaled';
        cc = cc';
    end
    cc = centroids_scaled;
    center = [3.5, 3.5, 3.5];
    cc = cc - center;
    bw = 0.7;
    figure;
    xi = -5:0.1:5; 
    [~,~,~,f1] = histDENSITY(cc(:,1),xi,bw);
    [~,~,~,f2] = histDENSITY(cc(:,2),xi,bw);
    [~,~,~,f3] = histDENSITY(cc(:,3),xi,bw);
    leng_line = max([max(f1) max(f2) max(f3)]);
    p1 = plot(xi(:),f1,'LineWidth',1.5, 'color', 'red'); hold on;
    p2 = plot(xi(:),f2,'LineWidth',1.5, 'Color', 'green');    
    p3 = plot(xi(:),f3,'LineWidth',1.5, 'Color', 'blue'); 
    line([nanmedian(cc(:,1)),nanmedian(cc(:,1))], [0 leng_line], 'color', 'red');
    line([nanmedian(cc(:,2)),nanmedian(cc(:,2))], [0 leng_line], 'color', 'green');
    line([nanmedian(cc(:,3)),nanmedian(cc(:,3))], [0 leng_line], 'color', 'blue');
    legend('A', 'B', 'C');
    xlabel('Distance from center'); ylabel('Proportion of fields');


    %%%%%%% Fig 4a3) Aligned or Fig 4b3) Tilted   %%%%%%%%
    xi = -0.5:.001:0.5;
    med_rand_x =[]; med_rand_y =[]; med_rand_z=[];
    for i = 1:1000
        med_rand_x = [med_rand_x nanmedian((5*rand(1,length(cc)))-2.5)];
        med_rand_y = [med_rand_y nanmedian((5*rand(1,length(cc)))-2.5)];
        med_rand_z = [med_rand_z nanmedian((5*rand(1,length(cc)))-2.5)];
    end
    med_rand = [med_rand_x; med_rand_y; med_rand_z]';
    colors = ['r'; 'g'; 'b'];
    figure;
    for j =1:3
        subplot(3,1,j);
        temp = med_rand(:,j);
        [~,~,~,f] = histDENSITY(temp,xi,0.01);
        p1 = plot(xi(:),f,'LineWidth',1.5, 'color', colors(j)); %hold on;
        line([nanmedian(cc(:,j)),nanmedian(cc(:,j))], [0 max(f)], 'color', 'black', 'LineWidth', 1.5);
        line([prctile(temp, 97.5), prctile(temp, 97.5)], [0 max(f)], 'color', 'red', 'LineWidth', 1.5);
        line([prctile(temp, 2.5), prctile(temp, 2.5)], [0 max(f)], 'color', 'red','LineWidth', 1.5);
        xlabel("Median field position");
        ylabel("Proportion of fields");
    end
end
        

%% %%%%%%% CALCULATE ORIENTATION AND ELONGATION %%%%%%%%%%%%
%orientation ,elongation

if calc_orientation_elongation
    orientation=double.empty;
    for i=1:length(pf)%calculate stats
%         i
        %calculate orientation
        b=pf(i).props.ConvexHull{1,1};
        a=pca(b);
        
        if tilted
            a = rotm * a;
            orientation=[orientation a(:,1)];
            orientationplot=[orientationplot a(:,1) -a(:,1)];
        else
            orientation=[orientation a(:,1)];
            orientationplot=[orientationplot a(:,1) -a(:,1)];
        end
        %calculate elongation
        t=pf(i).props;
        arr=t.PrincipalAxisLength;
        elong=2*arr(1)/(arr(2)+arr(3));
        elongation=[elongation ;elong];
        
    end
end

extent_field = [];

%for neuron=setdiff(neuron_start:neuron_end,avoid_neurons)
for neuron = neurons
    propstable=FIRR_3dprops(neuron).props;
    for i=1:height(propstable) 
            vol=propstable{i,"Volume"};
            if  vol>minv
             extent_field = [extent_field; propstable.BoundingBox(i,4:6)];   
             
            end
    end
end
figure
if tilted
    xi = 0:1:25;
    window_size = 10;
    [~,~,~,f1] = histDENSITY(extent_field(:,1),xi,1);
    [~,~,~,f2] = histDENSITY(extent_field(:,2),xi,1);
    [~,~,~,f3] = histDENSITY(extent_field(:,3),xi,1);
    smoothed_f1 = movmean(f1, window_size);
    smoothed_f2 = movmean(f2, window_size);
    smoothed_f3 = movmean(f3, window_size);
    p1 = plot(xi(:), smoothed_f1, 'color', 'r', 'LineWidth', 1.5);hold on;
    p2 = plot(xi(:), smoothed_f2, 'color', 'g', 'LineWidth', 1.5);
    p3 = plot(xi(:), smoothed_f3, 'color', 'b', 'LineWidth', 1.5);
    xlabel("Field extent");
    ylabel("Fields (%)");
    legend('A', 'B', 'C');
else
    xi = 0:1:25;
    [~,~,~,f1] = histDENSITY(extent_field(:,1),xi,2.5);
    [~,~,~,f2] = histDENSITY(extent_field(:,2),xi,2.5);
    [~,~,~,f3] = histDENSITY(extent_field(:,3),xi,2.5);
    p1 = plot(xi(:), f1, 'color', 'r', 'LineWidth', 1.5);hold on;
    p2 = plot(xi(:), f2, 'color', 'g', 'LineWidth', 1.5);
    p3 = plot(xi(:), f3, 'color', 'b', 'LineWidth', 1.5);
end  
%% %%%%%%% BINARY MORPHOLOGY FIG 7 %%%%%%%%%%%%
if plot_binary
    total=0;
    seq_len = 19;
    %create table like array that has SE length and correponding
    %proportions
    Px=double.empty;
    Py=double.empty;
    Pz=double.empty;
%     Px(1,:)=[3:2:19];
%     Py(1,:)=[3:2:19];
%     Pz(1,:)=[3:2:19];
    Px(1,:)=1:1:seq_len;
    Py(1,:)=1:1:seq_len;
    Pz(1,:)=1:1:seq_len;

    Px(2,:)=zeros(1,seq_len);
    Py(2,:)=zeros(1,seq_len);
    Pz(2,:)=zeros(1,seq_len);
    %seq_len_ratio = [];
    Seq_len_ratio = [];
    ind=1 ;   %index simply represents index of array l which contains structuring element length ;  
  for len=1:1:seq_len  %
    ex = strel('cuboid',[len 1 1]);    %elementx
    ey = strel('cuboid',[1 len 1]);
    ez = strel('cuboid',[1 1 len]);   
    total=0;
        
    %for neuron=setdiff(neuron_start:neuron_end,avoid_neurons) %binarize rate map of each neuron using 10% of maximum as thresh and make binary volumetric image
        
    %neuron
     for neuron = neurons
        bw=zeros(20,20,20);
        firr=FIRR_3dprops(neuron).firr;
        counter=1;
        threshold=0.1*max(firr);
        for k=1:20 % this loop is to create binary volumetric image for a neuron (with threshold)
            j=1;
            for j=1:20
                i=1;
                for i=1:20
                    %counter
                    if firr(counter)>threshold
                        bw(i,j,k)=1;
                    end
                    counter=counter+1;
                end
            end
        end
   
        Bx=imerode(bw,ex);
        By=imerode(bw,ey);
        Bz=imerode(bw,ez);
        
        total1=sum(Bx,"all")+sum(By,"all")+sum(Bz,"all");%total number of voxels with a particular connectivity for a particular neuron's firing rate
        total=total+total1; %total ,after complete iteration of all neurons,is the sum of all voxels remaining for that particular SE length 
        if neuron==neuron_start
        
            Px(2,ind)=sum(Bx,"all");
            Py(2,ind)=sum(By,"all");
            Pz(2,ind)=sum(Bz,"all");
            
        else
            Px(2,ind)=Px(2,ind)+sum(Bx,"all");%sum of remaining voxels of every iterated neuron along x axis
            Py(2,ind)=Py(2,ind)+sum(By,"all");
            Pz(2,ind)=Pz(2,ind)+sum(Bz,"all");
        end     
        neurons_ratio{neuron,1} = len;
        neurons_ratio{neuron,2} = sum(Bx,"all")/total1;
        neurons_ratio{neuron,3} = sum(By,"all")/total1;
        neurons_ratio{neuron,4} = sum(Bz,"all")/total1;
    end%neuron is inner loop
%    if total == 0
%        total = 0.00000001;
%    end
%      total = total + 0.0000001;
   Px(2,ind)=Px(2,ind)/total;
   Px
   Py(2,ind)=Py(2,ind)/total;
   Pz(2,ind)=Pz(2,ind)/total;
   ind=ind+1;
   %neurons_ratio = neurons_ratio(~isnan(neurons_ratio));
   Seq_len_ratio = vertcat(Seq_len_ratio, neurons_ratio);
   %seq_len_ratio{len} = neurons_ratio; 
  end %strel is outer loop
 
  
   Seq_len_ratio = cell2mat(Seq_len_ratio);
   %Seq_len_ratio = Seq_len_ratio(~isnan(Seq_len_ratio(:,2)), :);
    %Seq_len_ratio = Seq_len_ratio(~cellfun(@isempty, Seq_len_ratio), :);
    %Seq_len_ratio = Seq_len_ratio(~any(isnan(Seq_len_ratio), 2), :);     
    figure;
    plot(Px(1,:),Px(2,:), 'Color', 'r', 'LineWidth', 1.5)
    hold on 
    plot(Py(1,:),Py(2,:), 'Color', 'g', 'LineWidth', 1.5)
    hold on
    plot(Pz(1,:),Pz(2,:), 'Color', 'b', 'LineWidth', 1.5)
    ylabel('Proportion of voxels','FontSize',15)
    xlabel('Structuring element length','FontSize',15)
    if tilted
        legend('A','B','C')
    else
        legend('X','Y','Z')
    end
    xlim([0 20]); ylim([0 1]);
    title('Binary morphology','Fontsize',16)

end




%% %%%%%%% PLOT ELONGATION FIG 5 %%%%%%%%%%%%         
if plot_elongation
    if tilted
        load('field_elong_stats_TL.mat')
    else
        load('field_elong_stats_AL.mat')
    end

    %%%%%%%%%%% Fig 5a1) Aligned or Fig 5b1) Tilted  %%%%%%%%%%%
    hold off
     figure
    scatter(ones(size(elongation)),elongation,'filled','blue','jitter','on','jitterAmount',0.05);
    hold on 
    yline(1,'Color','black')
   
    hold on
    box1=boxplot(elongation','Whisker',5); hold on;
    set(box1,'Color','black');
    hold on
    ylim([0,5]); yticks([1, 2, 3, 4, 5]); xticks([1]);
    sph_fields = num2str(100 - round((sum(field_elong_stats(:,5))*100)/length(field_elong_stats),2)) + " %";
    hold on
    t = text(0.85,0.5,sph_fields); t.FontSize = 15;
    t=title("Elongation of place fields-Aligned",'FontSize',16);
    set(t);
    set(gca,'linewidth',2);
    ylabel("Field elongation",'FontSize',15);
    hold off
    %%%%%%%%%%% Fig 5a2) Aligned or Fig 5b2) Tilted  %%%%%%%%%%%

end

%% %%%%%%% PLOT ORIENTATION FIG 6 %%%%%%%%%%%%

if plot_orientation    
    if tilted==0
        figure
        x1=-1.2:0.1:1.2;
        y1=-1.2:0.1:1.2;
        z1=-1.2:0.1:1.2;
        x2=zeros(25,1);
        y2=zeros(25,1);
        z2=zeros(25,1);

        plot3(x1,y2,z2,'Color','black'); hold on
        plot3(x2,y1,z2,'Color','black'); hold on
        plot3(x2,y2,z1,'Color','black'); hold on
       
        V=double.empty;
        V =orientationplot.';
        [X,Y,Z,F2] = projectEIGS(V,256,10);  
        mx = nanmax(F2(:));    
        F2m = (F2 ./ mx);  

        s = surf(X,Y,Z,F2m);
        s.EdgeColor = 'k';
        s.EdgeAlpha = 0.0;
        colormap(gca,'jet');
        axis on vis3d; 
        rotate3d on;
        title("Place field orientation",'FontSize',16)
        xlabel('Y'); ylabel('X'); zlabel('Z'); hold off
            
    elseif tilted==1
            
        scale=-1.5:1:1.5; %scale of reference axes (tilted)
            %axis A
            i=1;
            j=0;
            k=0;
       
            i=i*scale;
            j=j*scale;
            k=k*scale;
            x_axis = [i; j; k];
            x_axis = rotm * x_axis;
            i = x_axis(1,:);
            j = x_axis(2,:);
            k = x_axis(3,:);
            figure()
            plot3(i,j,k, 'r', 'LineWidth', 2);
            hold on
            
            
            
            %axis B
            i = 0;
            j = 1;
            k = 0;
            
            i=i*scale;
            j=j*scale;
            k=k*scale;
            y_axis = [i; j; k];
            y_axis = rotm * y_axis;
            i = y_axis(1,:);
            j = y_axis(2,:);
            k = y_axis(3,:);
            plot3(i,j,k, 'g', 'LineWidth', 2);
            hold on            
            
            %axis C
            i = 0;
            j = 0;
            k = 1;
            
            i=i*scale;
            j=j*scale;
            k=k*scale;
            z_axis = [i; j; k];
            z_axis = rotm * z_axis;
            i = z_axis(1,:);
            j = z_axis(2,:);
            k = z_axis(3,:);
            plot3(i,j,k, 'b', 'LineWidth', 2);
           
            
            
            V=double.empty;
            V =orientationplot';
            [X,Y,Z,F2] = projectEIGS(V,256,10);  
            mx = nanmax(F2(:));    
            F2m = (F2 ./ mx);  

            s = surf(X,Y,Z,real(F2m));
            s.EdgeColor = 'k';
            s.EdgeAlpha = 0.0;
            colormap(gca,'jet');
            axis on vis3d; 
            rotate3d on;
            title("Place field orientation")
            xlabel('Y')
            ylabel('X')
            zlabel('Z')
            
    end    


    %%Statistical analysis
    theta_bool = [];
    theta_bool1 = [];
    
    %collect all primary eigenvectors
    ev = [];
    %for neuron=setdiff(neuron_start:neuron_end,avoid_neurons)
    for neuron = neurons
        Vol = FIRR_3dprops(neuron).props.Volume;
        for i = 1:length(Vol) 
            if (Vol(i)>50)
                temp2 = FIRR_3dprops(neuron).props.EigenVectors(i);
                ev = [ev temp2{1,1}(:,1)];
            end
        end
    end
    
    %check if eigenvectors are in vicinity of axes    
    if tilted
        ev_tilted = rotm*ev;
        ev_tilted = ev_tilted';
        for i = 1:length(ev_tilted) 
            ev_temp = ev_tilted(i,:);
            theta = acosd(abs(ev_temp*XYZ'));
            temp = theta<30;
            theta_bool = [theta_bool; temp];         
        end
        ev= ev';
        for i = 1:length(ev) 
            ev_temp2 = ev(i,:);
            theta = acosd(abs(ev_temp2*XYZ'));
            temp = theta<30;
            theta_bool1 = [theta_bool1; temp];         
        end
    else
        ev= ev';
        for i = 1:length(ev) 
            ev_temp = ev(i,:);
            theta = acosd(abs(ev_temp*axes_cords'));
            temp = theta<30;
            theta_bool1 = [theta_bool1; temp];         
        end
    end

    theta_bool = [theta_bool theta_bool1];
    
    ratio_ev = sum(theta_bool,1)/length(theta_bool);
    xyz_abc = sum(ratio_ev(1:3))/sum(ratio_ev(4:6));
    
    %%error bars and confidence intervals
    % randomly choose pfs from uniform distribution and create error bars
    shuffled_ratio = [];
    if tilted
        for i=1:1000
            rindx = randi(length(ev), 1,length(ev));
            shuffl_pf_t = ev_tilted(rindx,:);
            theta_dat = []; theta_dat1 = [];
    
            for j = 1:length(shuffl_pf_t) 
                temp = shuffl_pf_t(j,:);
                theta = acosd(abs(temp*XYZ'));
                temp = theta<30;
                theta_dat = [theta_dat; temp];         
            end
            shuffl_pf = ev(rindx,:);
            for j = 1:length(shuffl_pf) 
                ev_temp2 = shuffl_pf(j,:);
                theta = acosd(abs(ev_temp2*XYZ'));
                temp = theta<30;
                theta_dat1 = [theta_dat1; temp];         
            end
            theta_dat2 = [theta_dat theta_dat1];
            shuffled_ratio = [shuffled_ratio; sum(theta_dat2,1)/length(theta_bool)];
        end
   else
        for i=1:1000
            rindx = randi(length(ev), 1,length(ev));
            shuffl_pf = ev(rindx,:);
            theta_dat2 = [];
            for j = 1:length(shuffl_pf) 
                temp = shuffl_pf(j,:);
                theta = acosd(abs(temp*axes_cords'));
                temp = theta<30;
                theta_dat2 = [theta_dat2; temp];         
            end
            shuffled_ratio = [shuffled_ratio; sum(theta_dat2,1)/length(theta_bool)];
        end
    end
    


    %%%%%%%%%%% Fig 6a1) Aligned or Fig 6b1) Tilted  %%%%%%%%%%%
    %%with random points
    % randomly generate points to know the by chance line and also the XYZ/ABC
    % distribution
    figure;
    total = [];
    for i = 1:1000
        TH = 2*pi*rand(1,1000);
        PH = asin(-1+2*rand(1,1000));
        [X,Y,Z] = sph2cart(TH,PH,1);
        %for j = 1:1000
        ev_rand = [X; Y; Z];
        ev_rand = ev_rand';
        theta = acosd(abs(ev_rand*axes_cords'));
        temp = theta<30;
        total = [total; sum(temp, 1)/length(temp)];   
    end
    
    rand_shuffl_ratio = sum(total(:,1:3), 2)./sum(total(:,4:6), 2);
    
    % plotting distribution
    xi = 0:0.01:2;
    dist = sum(total(:,1:3), 2)./sum(total(:,4:6), 2);
    llev = prctile(rand_shuffl_ratio, 1);
    hlev = prctile(rand_shuffl_ratio, 99);
    [~, ~, ~, fn] = histDENSITY(dist, xi, 0.09);
    % p = plot(xi(:),fn,'LineWidth',1.5, 'color', 'black');
    area(xi,fn,'FaceColor', [0.9 0.9 0.9],'LineWidth',1.5);
    line([xyz_abc xyz_abc], [min(fn), max(fn)], 'color', 'b', 'LineWidth', 1.5);
    line([llev llev], [min(fn), max(fn)], 'color', 'r', 'LineWidth', 1.5);
    line([hlev hlev], [min(fn), max(fn)], 'color', 'r', 'LineWidth', 1.5);
    xlabel('XYZ/ABC'); ylabel('Probability')
    



    %%%%%%%%%%% Fig 6a2) Aligned or Fig 6b1) Tilted  %%%%%%%%%%%
    %%plotting errorbar
    figure;
    prcr1 =  prctile(total, 2.5); % percentiles for chance line
    prcr2 =  prctile(total, 50);
    prcr3 =  prctile(total, 97.5);
    
    prcs1 = prctile(shuffled_ratio, 2.5); % percentile for error bar
    prcs2 = prctile(shuffled_ratio, 97.5);
    for i = 1:6       
        errorbar(i,ratio_ev(i),ratio_ev(i)-prcs1(i),ratio_ev(i)-prcs2(i), 'LineWidth', 1.0); hold on;
        line([0.5 7], [prcr1(i) prcr1(i)], 'color', 'red', 'LineWidth', 0.5);
        line([0.5 7], [prcr3(i) prcr3(i)], 'color', 'red', 'LineWidth', 0.5);
        line([0.5 7], [prcr2(i) prcr2(i)], 'color', 'red', 'LineWidth', 1.5);
        plot(i,ratio_ev(i), 'Marker','o', 'LineWidth', 1.5, 'MarkerSize', 15);
        xlim([0.5,7]);
        ylim([0,0.8]);
    end
    names = {'X'; 'Y'; 'Z'; 'A'; 'B'; 'C'};
    
    set(gca,'xtick',1:6,'xticklabel',names)
    xlabel('Maze Axis'); ylabel('Proportion of fields')
end