%Read data associated with tNORTs, perform stats analysis and generate Figure 4     
%  Ying Zheng, 04/03/2024. Version 1

clear all; close all; home;
%load data
tNORTdata(:,:,1) = xlsread('Data_tNORTs','sheet1','B3:E10'); % L/R/F/N, rough      
tNORTdata(:,:,2) = xlsread('Data_tNORTs','sheet1','H3:K10');  %dim: an x LRFN x paper        
tNORTdata(:,:,3) = xlsread('Data_tNORTs','sheet1','N3:Q10');     
tNORTdata(:,:,4) = xlsread('Data_tNORTs','sheet1','T3:W10');     
tNORTdata(:,:,5) = xlsread('Data_tNORTs','sheet1','Z3:AE10'); %smooth    
[Nrat,Nside,Npaper]=size(tNORTdata);

% calculate D index
for i=1:Npaper
    di_LR(:,i)=(tNORTdata(:,1,i)-tNORTdata(:,2,i))./(tNORTdata(:,1,i)+tNORTdata(:,2,i));
    di_FN(:,i)=(tNORTdata(:,4,i)-tNORTdata(:,3,i))./(tNORTdata(:,4,i)+tNORTdata(:,3,i));
end

%plot di
figure('pos',[600 200 500 100]);
m_di_LR=mean(di_LR); sem_di_LR=std(di_LR)/sqrt(Nrat);
bar(m_di_LR,0.5,'FaceColor', [1 1 1]);hold on; errorbar(m_di_LR,sem_di_LR,'.k','LineWidth',1);
axis([0 6 -1 1]);
for i=1:Npaper
    swarmchart(i+0.05*randn(8,1), di_LR(:,i),'.b');
end
set(gca,'XTicklabel',{''},'Xtick',0,...
    'YTicklabel',[{'-1'},{'-0.5'},{'0'},{'0.5'},{'1'}],'Ytick',[-1 -0.5 0 0.5 1],'FontName','Arial','FontSize',8);
box off;
set(gcf,'color','w');

figure('pos',[600 400 500 100]);
m_di_FN=mean(di_FN); sem_di_FN=std(di_FN)/sqrt(Nrat);
bar(m_di_FN,0.5,'FaceColor', [1 1 1]);hold on; errorbar(m_di_FN,sem_di_FN,'.k','LineWidth',1);
axis([0 6 -1 1]);
for i=1:Npaper
    swarmchart(i+0.05*randn(8,1), di_FN(:,i),'.b');
end
set(gca,'XTicklabel',{''},'Xtick',0,...
    'YTicklabel',[{'-1'},{'-0.5'},{'0'},{'0.5'},{'1'}],'Ytick',[-1 -0.5 0 0.5 1],'FontName','Arial','FontSize',8);
box off;
set(gcf,'color','w');

%screen for outliers
temp1=filloutliers(di_LR,"previous","mean");   max(max(abs(temp1-di_LR)))    %no outliers    
temp2=filloutliers(di_FN,"previous","mean");   max(max(abs(temp2-di_FN)))    %no outliers  
%normality test for di. All Gaussian
for i=1:Npaper
    [H(i,1)]=jbtest(di_LR(:,i));
    [H(i,2)]=jbtest(di_FN(:,i));
end

%t-test on di
for i=1:Npaper
    [h1(i),p1(i)]=ttest(di_LR(:,i));
    [h2(i),p2(i)]=ttest(di_FN(:,i));
end

%plot particle size
figure('pos',[600 600 500 100]);
sand=[201 125 82 58 35];
plot(1:5, sand,'ob');axis([0 6 0 205]);
set(gca,'XTicklabel',[{'P80'},{'P120'},{'P180'},{'P240'},{'P400'}],'Xtick',[1 2 3 4 5],...
    'YTicklabel',[{'0'},{'50'},{'100'},{'150'},{'200'}],'Ytick',[0 50 100 150 200],'FontName','Arial','FontSize',8);
box off;
set(gcf,'color','w');