
close all;
clear all;
%rng(pi) % just to reproduce the random data I used
proj39=proj39collection;
t2A=readtable('T2AwithBiregression.csv');
t12A=readtable('T12A_20240207withBiregression.csv');
t12=readtable('T12_20240207CosSimilarity500.csv');
t1=readtable('T2AwithBiregression.csv');
x=categorical({'Within','Between'});
%x=categorical({'within'});
y1=t1.consistencyWithin;
y2=t1.consistencyBetween;
y=[y1, y2];
h=barFigure(x,y,'CImean',1);
ylabel('Consistency')
saveas(h, 'Consistency', 'fig');
%h2=barFigure(x,y,'stderr',1);
%legend('cos(x)','cos(2x)')

%% circular figure of heading errors


%t12A.E1headingError0=proj39.VectorCal.ChangeAngleTo_0to360(t12A.E1headingError0);
%t12A.E1headingError1=proj39.VectorCal.ChangeAngleTo_0to360(t12A.E1headingError1);

%change to 0-360 for oriana
oriana=table();
oriana.E1headingError0=proj39.VectorCal.ChangeAngleTo_0to360(t12A.E1headingError0);
oriana.E1headingError1=proj39.VectorCal.ChangeAngleTo_0to360(t12A.E1headingError1);
oriana.E2headingError0=proj39.VectorCal.ChangeAngleTo_0to360(t12A.E2headingError0);
oriana.E2headingError1=proj39.VectorCal.ChangeAngleTo_0to360(t12A.E2headingError1);
oriana.E3headingError0=proj39.VectorCal.ChangeAngleTo_0to360(t12A.E3headingError0);
oriana.E3headingError1=proj39.VectorCal.ChangeAngleTo_0to360(t12A.E3headingError1);
oriana.E4headingError0=proj39.VectorCal.ChangeAngleTo_0to360(t12A.E4headingError0);
oriana.E4headingError1=proj39.VectorCal.ChangeAngleTo_0to360(t12A.E4headingError1);
oriana.E5headingError0=proj39.VectorCal.ChangeAngleTo_0to360(t12A.E5headingError0);
oriana.E5headingError1=proj39.VectorCal.ChangeAngleTo_0to360(t12A.E5headingError1);
oriana.E6headingError0=proj39.VectorCal.ChangeAngleTo_0to360(t12A.E6headingError0);
oriana.E6headingError1=proj39.VectorCal.ChangeAngleTo_0to360(t12A.E6headingError1);
writetable(oriana,'orianaHeadingError.csv');

proj39.visualizeAnglePairs(t12A.E1headingError0,t12A.E1headingError1);
legendName={'Indoor','','','','Outdoor','','',''};
legend(legendName);
title('Bu');

proj39.visualizeAnglePairs(t12A.E2headingError0,t12A.E2headingError1);
legend(legendName);
title('CC');
proj39.visualizeAnglePairs(t12A.E3headingError0,t12A.E3headingError1);
legend(legendName);
title('Ca');
proj39.visualizeAnglePairs(t12A.E4headingError0,t12A.E4headingError1);
legend(legendName);
title('Hu');
proj39.visualizeAnglePairs(t12A.E5headingError0,t12A.E5headingError1);
legend(legendName);
title('Pe');
proj39.visualizeAnglePairs(t12A.E6headingError0,t12A.E6headingError1);
legend(legendName);
title('Ru');

%% for pointingerror

% y1=[t12A.E2abs_pointerror0 t12A.E4abs_pointerror0 t12A.E3abs_pointerror0 t12A.E6abs_pointerror0 t12A.E1abs_pointerror0 t12A.E5abs_pointerror0];
% y1=flip(y1,2);
% 
% y2=[t12A.E2abs_pointerror1 t12A.E4abs_pointerror1 t12A.E3abs_pointerror1 t12A.E6abs_pointerror1 t12A.E1abs_pointerror1 t12A.E5abs_pointerror1];
% y2=flip(y2,2)
% 
% y=[y1, y2]
% 
% x={'ccis4.46','hub4.43','cameron4.19','ruth3.978','butter3.63','pe1.52'}
% x=flip(x,2)

buildingNames={'Bu','CC','Ca','Hu','Pe','Ru'};
ylimRange=[0 120];

PointingErrorArrays12=[t12A.E1abs_pointerror0 t12A.E2abs_pointerror0 t12A.E3abs_pointerror0 t12A.E4abs_pointerror0 t12A.E5abs_pointerror0 t12A.E6abs_pointerror0,t12A.E1abs_pointerror1 t12A.E2abs_pointerror1 t12A.E3abs_pointerror1 t12A.E4abs_pointerror1 t12A.E5abs_pointerror1 t12A.E6abs_pointerror1];
xlabelName={'Indoor','Outdoor'};
h=ErrorFigure (PointingErrorArrays12, buildingNames, 'mse',0,'Absolute pointing error','degree',xlabelName, ylimRange);
saveas(h, 'Absolute_pointing_error12A', 'fig');

y=[t2A.abs_pointerror0, t2A.abs_pointerror1];
x=categorical(xlabelName);
ylabelName='Absolute pointing error';
h=ErrorFigure (y, x, 'mse',1,ylabelName,'degree',{}, ylimRange);
saveas(h, 'Absolute_pointing_error2A', 'fig');

%%
% HotelingT=collectiveHotelingTsquared
% observedt=6.8
% n=46
% observedD=observedt*sqrt(2/n)
% withinOrBetween=1
% BF01=HotelingT.HotelingBF01BasedonD_cauchy(observedD, 1, n, withinOrBetween,1/sqrt(2)); 

% for pointing error mixed linear model as testing envionrment has effect
%[abs_pointerror.continuous_unfam_fam.mdLME, abs_pointerror.continuous_unfam_fam.mdLME1,abs_pointerror.continuous_unfam_fam.mdl_angle]=ScatterFigure(t12,'continuous_unfam_fam','abs_pointerror','abs pointing error','continuous unfam fam','degree',proj39);
[h,abs_pointerror.continuous_unfam_fam_location.mdLME, abs_pointerror.continuous_unfam_fam_location.mdLME1,abs_pointerror.continuous_unfam_fam_location.mdl_angle]=ScatterFigure(t12,'continuous_unfam_fam_location','abs_pointerror', ylabelName,'Independent viewing building familarity','degree',proj39);
saveas(h, 'abs_pointerror_Independent_viewing_building_familarity', 'fig');
[h,abs_pointerror.envrionmentFamilarity.mdLME,abs_pointerror.envrionmentFamilarity.mdLME1,abs_pointerror.envrionmentFamilarity.mdl_angle]=ScatterFigure(t12,'envrionmentFamilarity','abs_pointerror', ylabelName,'Self viewing building familarity','degree',proj39);
saveas(h, 'abs_pointerror_Self_viewing_building_familarity', 'fig');
[h,abs_pointerror.envrionmentFamilarityMean.mdLME, abs_pointerror.envrionmentFamilarityMean.mdLME1,abs_pointerror.envrionmentFamilarityMean.mdl_angle]=ScatterFigure(t12,'envrionmentFamilarityMean','abs_pointerror', ylabelName,'Mean self viewing building familarity','degree',proj39);
saveas(h, 'abs_pointerror_Mean_Self_viewing_building_familarity', 'fig');
abs_pointerror.continuous_unfam_fam_location.mdLME1
abs_pointerror.envrionmentFamilarity.mdLME1
abs_pointerror.envrionmentFamilarityMean.mdLME1
%% beta and se for beta
x=t12.continuous_unfam_fam_location ;
X=[ones(length(x(:,1)),1), x];
Y=t12.abs_pointerror;
Beta=inv(X'*X)*X'*Y;
SSE=(Y-X*Beta)'*(Y-X*Beta);
p=2;
n=length(x(:,1));
MSE=SSE/(n-p-1);
errorvarBeta=MSE*inv(X'*X);
seBeta=sqrt([errorvarBeta(1,1);errorvarBeta(2,2)]);
t=Beta./seBeta;

%% for RT E2RT0

ylabelName='Response latency';
yRTlimRange=[0,5];
RTArrays12=[t12A.E1RT0 t12A.E2RT0 t12A.E3RT0 t12A.E4RT0 t12A.E5RT0 t12A.E6RT0,t12A.E1RT1 t12A.E2RT1 t12A.E3RT1 t12A.E4RT1 t12A.E5RT1 t12A.E6RT1];
h=ErrorFigure (RTArrays12, buildingNames, 'mse',0,ylabelName,'second',xlabelName, yRTlimRange);
saveas(h, 'RT12A', 'fig');
y=[t2A.RT0, t2A.RT1];
x=categorical(xlabelName);
h=ErrorFigure (y, x, 'mse',1,ylabelName,'second',{}, yRTlimRange);
saveas(h, 'RT2A', 'fig');
% for mixed linear model as testing envionrment has effect
[h,RT.continuous_unfam_fam_location.mdLME, RT.continuous_unfam_fam_location.mdLME1,RT.continuous_unfam_fam_location.mdl_angle]=ScatterFigure(t12,'continuous_unfam_fam_location','RT', ylabelName,'Independent viewing building familarity','second',proj39);
saveas(h, 'RT_Independent_viewing_building_familarity', 'fig');

[h,RT.envrionmentFamilarity.mdLME, RT.envrionmentFamilarity.mdLME1,RT.envrionmentFamilarity.mdl_angle]=ScatterFigure(t12,'envrionmentFamilarity','RT', ylabelName,'Self viewing building familarity','second',proj39);
saveas(h, 'RT_Self_viewing_building_familarity', 'fig');

[h,RT.envrionmentFamilarityMean.mdLME, RT.envrionmentFamilarityMean.mdLME1,RT.envrionmentFamilarityMean.mdl_angle]=ScatterFigure(t12,'envrionmentFamilarityMean','RT', ylabelName,'Mean self viewing building familarity','second',proj39);
saveas(h, 'RT_Mean_Self_viewing_building_familarity', 'fig');

RT.continuous_unfam_fam_location.mdLME1
RT.envrionmentFamilarity.mdLME1
RT.envrionmentFamilarityMean.mdLME1
%% for abs heading E2absheadingError0
%t12A.E1absheadingError0
ylabelName='Absolute heading error'
HeadingErrorArrays12=[t12A.E1absheadingError0 t12A.E2absheadingError0 t12A.E3absheadingError0 t12A.E4absheadingError0 t12A.E5absheadingError0 t12A.E6absheadingError0,t12A.E1absheadingError1 t12A.E2absheadingError1 t12A.E3absheadingError1 t12A.E4absheadingError1 t12A.E5absheadingError1 t12A.E6absheadingError1];
h=ErrorFigure (HeadingErrorArrays12, buildingNames, 'mse',0,ylabelName,'degree',xlabelName, ylimRange);
saveas(h, 'absheadingError12A', 'fig');
y=[t2A.absheadingError0, t2A.absheadingError1];
%x=categorical(xlabelName);
h=ErrorFigure (y, x, 'mse',1,ylabelName,'degree',{}, ylimRange);
saveas(h, 'absheadingError2A', 'fig');
% for mixed linear model as testing envionrment has effect
[h,absheadingError.continuous_unfam_fam_location.mdLME, absheadingError.continuous_unfam_fam_location.mdLME1,absheadingError.continuous_unfam_fam_location.mdl_angle]=ScatterFigure(t12,'continuous_unfam_fam_location','absheadingError', ylabelName,'Independent viewing building familarity','degree',proj39);
saveas(h, 'absheadingError_Independent_viewing_building_familarity', 'fig');

[h,absheadingError.envrionmentFamilarity.mdLME,absheadingError.envrionmentFamilarity.mdLME1,absheadingError.envrionmentFamilarity.mdl_angle]=ScatterFigure(t12,'envrionmentFamilarity','absheadingError', ylabelName,'Self viewing building familarity','degree',proj39);
saveas(h, 'absheadingError_Self_viewing_building_familarity', 'fig');

[h,absheadingError.envrionmentFamilarityMean.mdLME, absheadingError.envrionmentFamilarityMean.mdLME1,absheadingError.envrionmentFamilarityMean.mdl_angle]=ScatterFigure(t12,'envrionmentFamilarityMean','absheadingError', ylabelName,'Mean self viewing building familarity','degree',proj39);
saveas(h, 'absheadingError_Mean_Self_viewing_building_familarity', 'fig');

absheadingError.continuous_unfam_fam_location.mdLME1
absheadingError.envrionmentFamilarity.mdLME1
absheadingError.envrionmentFamilarityMean.mdLME1
%absheadingError.envrionmentFamilarityMean.mdLME;

%% for positionError E2abspositionAngularError1
ylabelName='Absolute position angular error'
HeadingErrorArrays12=[t12A.E1abspositionAngularError0 t12A.E2abspositionAngularError0 t12A.E3abspositionAngularError0 t12A.E4abspositionAngularError0 t12A.E5abspositionAngularError0 t12A.E6abspositionAngularError0,t12A.E1abspositionAngularError1 t12A.E2abspositionAngularError1 t12A.E3abspositionAngularError1 t12A.E4abspositionAngularError1 t12A.E5abspositionAngularError1 t12A.E6abspositionAngularError1];
h=ErrorFigure (HeadingErrorArrays12, buildingNames, 'mse',0,ylabelName,'degree',xlabelName, ylimRange);
saveas(h, 'abspositionAngularError12A', 'fig');

y=[t2A.abspositionAngularError0, t2A.abspositionAngularError1];
%x=categorical({'indoor','outdoor'});
h=ErrorFigure (y, x, 'mse',1,ylabelName,'degree',{}, ylimRange);
saveas(h, 'abspositionAngularError2A', 'fig');

% for mixed linear model as testing envionrment has effect
[h,abspositionAngularError.continuous_unfam_fam_location.mdLME, abspositionAngularError.continuous_unfam_fam_location.mdLME1,abspositionAngularError.continuous_unfam_fam_location.mdl_angle]=ScatterFigure(t12,'continuous_unfam_fam_location','abspositionAngularError', ylabelName,'Independent viewing building familarity','degree',proj39);
saveas(h, 'abspositionAngularError_Independent_viewing_building_familarity', 'fig');

[h,abspositionAngularError.envrionmentFamilarity.mdLME,abspositionAngularError.envrionmentFamilarity.mdLME1,abspositionAngularError.envrionmentFamilarity.mdl_angle]=ScatterFigure(t12,'envrionmentFamilarity','abspositionAngularError',ylabelName,'Self viewing building familarity','degree',proj39);
saveas(h, 'abspositionAngularError_Self_viewing_building_familarity', 'fig');

[h,abspositionAngularError.envrionmentFamilarityMean.mdLME, abspositionAngularError.envrionmentFamilarityMean.mdLME1,abspositionAngularError.envrionmentFamilarityMean.mdl_angle]=ScatterFigure(t12,'envrionmentFamilarityMean','abspositionAngularError', ylabelName,'Mean self viewing building familarity','degree',proj39);
saveas(h, 'abspositionAngularError_Mean_Self_viewing_building_familarity', 'fig');

abspositionAngularError.continuous_unfam_fam_location.mdLME1
abspositionAngularError.envrionmentFamilarity.mdLME1
abspositionAngularError.envrionmentFamilarityMean.mdLME1


%% for E2absPointingErrorRemovingPositionHeadingErrors

%% for  position distance error E2distanceError0

%% present targets and testing locations
t3=readtable('T6020231201CosSimilarity500.csv');
[M,I] = min(t3.leastErrSquared)
index=t3.participant_number==t3.participant_number(I)&strcmp(t3.envrionment_building,t3.envrionment_building(I));
targetNames=t3.target_building(index);
targetsLocations=[t3.TargetX(index),t3.TargetY(index)];

index1=t3.participant_number==t3.participant_number(I)&strcmp(t3.target_building,'CAB');
%environmentNames=t3.envrionment_building(index1);


environmentNames=t3.Enviornmentname(index1);
indoor_outdoorValue=t3.indoor_outdoorValue(index1);

for i=1:length(environmentNames)
     environmentNames{i}=[environmentNames{i}(1:2),'(',indoor_outdoorValue{i},')'];
end

environmentLocations=[t3.positionX(index1),t3.positionY(index1)];

heading=t3.facing_direction(index);
h=heading(1);



figure
arrowLength=80;
arrowLengthLong=100;
ArrowsizeS=1;

for i=1:length(environmentNames)  
    P=environmentLocations(i,:);
    PX=P(1);
    PY=P(2);
    plot (PX,PY,'k.')
    text (PX, PY,environmentNames{i},'Color',  'k')
    hold on
    drawArrowGivenPandAngle(PX,PY,h,arrowLength,ArrowsizeS,'h','k');
    hold on
end

for i=1:length(targetNames)  
    P=targetsLocations(i,:);
    PX=P(1);
    PY=P(2);
    plot (PX,PY,'k*');
    text (PX, PY,targetNames{i},'Color',  'k');
    hold on
end

drawArrowGivenPandAngle(-500,300,h,arrowLength,ArrowsizeS,'h','r');
text (-500, 300,'location building','Color',  'r');
hold on
plot (-500,300-arrowLengthLong,'r*');
text (-500+3, 300-arrowLengthLong,'target building','Color',  'r');

%% get indoor-outdoor for pointingerror headingerror positionerror motionerror and then do regression
t12A.abs_pointingerrorD=t2A.abs_pointerror0-t2A.abs_pointerror1;
t12A.absheadingErrorD=t2A.absheadingError0-t2A.absheadingError1;
t12A.abspositionAngularErrorD=t2A.abspositionAngularError0-t2A.abspositionAngularError1;
t12A.absPointingErrorRemovingPositionHeadingErrorsD=t2A.absPointingErrorRemovingPositionHeadingErrors0-t2A.absPointingErrorRemovingPositionHeadingErrors1;
t12A.rsqprime2cD=t12A.rsqinprime2inc-t12A.rsqoutprime2outc;

 r=corr([t12A.abs_pointingerrorD,t12A.absheadingErrorD,t12A.abspositionAngularErrorD])
corr([t12A.absheadingErrorD,t12A.abspositionAngularErrorD,t12A.absPointingErrorRemovingPositionHeadingErrorsD,t12A.rsqprime2cD]);

%X=[t12A.absheadingErrorD,t12A.abspositionAngularErrorD,t12A.absPointingErrorRemovingPositionHeadingErrorsD]
X=[t12A.absheadingErrorD,t12A.abspositionAngularErrorD];
y=t12A.abs_pointingerrorD;
mdl = fitlm(X,y);
T=array2table([t12A.abs_pointingerrorD,t12A.absheadingErrorD,t12A.abspositionAngularErrorD,t12A.absPointingErrorRemovingPositionHeadingErrorsD],"VariableNames",{'abs_pointingerrorD','absheadingErrorD','abspositionAngularErrorD','absPointingErrorRemovingPositionHeadingErrorsD'})
[mdregression1,tb1]=proj39.myregression (T, {'abs_pointingerrorD','absheadingErrorD','abspositionAngularErrorD'})
%[mdregression1,tb1]=proj39.myregression (T, {'abs_pointingerrorD','absheadingErrorD','abspositionAngularErrorD','absPointingErrorRemovingPositionHeadingErrorsD'})
%writetable(t12A,'T12A_20240207withBiregression.csv');

(mean(t2A.abs_pointerror0))^2;
(mean(t2A.absheadingError0))^2+(mean(t2A.abspositionAngularError0))^2+(mean(t2A.absPointingErrorRemovingPositionHeadingErrors0))^2;

mean(t2A.abs_pointerror0.^2);
mean(t2A.absheadingError0.^2)+mean(t2A.abspositionAngularError0.^2)+(mean(t2A.absPointingErrorRemovingPositionHeadingErrors0.^2));

Foutput=anova(mdregression1);
stdY=std(y);
stdX=std(X);
beta=mdregression1.Coefficients.Estimate(1:end)';
standardBeta=beta(2:end).*stdX./stdY;
X1=[ones(length(t12A.absheadingErrorD),1),X];
mse=Foutput{'Error','MeanSq'};
BetaSE= sqrt(mse*diag(inv(X1'*X1)));
Tvalues=beta'./BetaSE(1:end);
corrxy=corr([y,X]);
Rsq=standardBeta*corrxy(1,2:end)';

%% present an example of represented position and heading from five pointing directions
t3=readtable('T6020231201CosSimilarity500.csv');
[M,I] = min(t3.leastErrSquared);


%index=t3.participant_number==t3.participant_number(I)&strcmp(t3.envrionment_building,t3.envrionment_building(I))
index=t3.participant_number==9&strcmp(t3.envrionment_building,"Cameron Library(outside)");
targetNames=t3(index,'target_building');
targetsLocations=[t3(index,"TargetX") t3(index,"TargetY")];
responseAngle=t3(index,"pointing_angle");
respresentedLocation=[t3(index,'positionprimeX'),t3(index,'positionprimeY')];
representedHeading=t3(index,"headingprime");
location=[t3(index,"positionX"),t3(index,"positionY")];
heading=t3(index,"facing_direction");

figure
arrowLength=120;
arrowLengthLong=150;
ArrowsizeS=1;

P=[location(1,:).positionX, location(1,:).positionY];
PX=P(1);
PY=P(2);
h=heading(1,:).facing_direction;


plot (targetsLocations.TargetX, targetsLocations.TargetY,'k*');
hold on

text (PX, PY,'P','Color',  'k');
drawArrowGivenPandAngle(PX,PY,h,arrowLengthLong,ArrowsizeS,['h'],'k');

PprimeX=respresentedLocation(1,:).positionprimeX;
PprimeY=respresentedLocation(1,:).positionprimeY;
hprime=representedHeading(1,:).headingprime;
% plot (PprimeX, PprimeY,'b.')
% text(PprimeX+0.1, PprimeY,['P',char(39)], 'Color', 'b');
% %text(PprimeX+arrowLength*sind(hprime),PprimeY,+arrowLength*cosd(hprime),['h',char(39),'1'], 'Color', 'r');
% hold on
% 
% quiver(PprimeX, PprimeY,arrowLength*sind(hprime),arrowLength*cosd(hprime),  0, 'LineWidth',1, 'MaxHeadSize',ArrowsizeS,'Color', 'b');
% text(PprimeX+arrowLength*sind(hprime),PprimeY+arrowLength*cosd(hprime)+0.2,['h',char(39)], 'Color', 'b');
% hold on
text(PprimeX,PprimeY,['P',char(39)],'Color', 'r');
drawArrowGivenPandAngle(PprimeX,PprimeY,hprime,arrowLengthLong,ArrowsizeS,['h',char(39)],'r');

pointAngle=responseAngle.pointing_angle;

for i=1:length(pointAngle)
    drawArrowGivenPandAngle(PX,PY,h+pointAngle(i),arrowLength,ArrowsizeS,[targetNames.target_building{i}(1:2),char(39)],'r');
end

for i=1:length(pointAngle)
    targetLocation=[targetsLocations.TargetX(i),targetsLocations.TargetY(i)];
    targetName=targetNames.target_building{i};
    text(targetLocation(1),targetLocation(2),targetName);

    bearing = getbearingA ([PprimeX,PprimeY], targetLocation,hprime);
    drawArrowGivenPandAngle(PprimeX,PprimeY,hprime+bearing,arrowLength,ArrowsizeS,[targetNames.target_building{i}(1:2)],'k');
end


%%
function drawArrowGivenPandAngle(PprimeX,PprimeY,hprime,arrowLength,ArrowsizeS,textName,colorName)
    %PprimeX=respresentedLocation(1,:).positionprimeX
    %PprimeY=respresentedLocation(1,:).positionprimeY;
    %hprime=representedHeading(1,:).headingprime;
    %plot (PprimeX, PprimeY,[colorName,'.'])
    %text(PprimeX+0.1, PprimeY,['P',char(39)], 'Color', colorName);
    %text(PprimeX+arrowLength*sind(hprime),PprimeY,+arrowLength*cosd(hprime),['h',char(39),'1'], 'Color', 'r');
    hold on
    
    quiver(PprimeX, PprimeY,arrowLength*sind(hprime),arrowLength*cosd(hprime),  0, 'LineWidth',1, 'MaxHeadSize',ArrowsizeS,'Color', colorName);
    text(PprimeX+arrowLength*sind(hprime)+0.2,PprimeY+arrowLength*cosd(hprime)+0.2,textName, 'Color', colorName);
    hold on
end



function h=barFigure(x,y,errbarType,scatterFlag)%errbarType='stderr','CImean', ='mse', ='msej'
    %y=rand(10,2)*6;
    figure
    meany=mean(y);
    rowY=length(y(1,:))/length(x);
    columnX=length(x);
    %if length(y)>length(x)
    meany1=reshape(meany,[length(x),rowY]);
    meany1=meany1';
    %end
    h=bar(x, meany1');
    hold on
    
    %h.FaceColor='flat';
    % h.CData(1,:)=1
    % h.CData(2,:)=2
    % h.CData(1,:)=[0.5,0,0.5]
    % h.CData(2,:)=[0,1,0]
    % h(1).FaceColor='r';
    % h(2).FaceColor='y';
    n=length(y(:,1));
    df=n-1;
    switch errbarType
        case 'CImean'
            errbar=std(y)/sqrt(n)*tinv(0.975,df);
        case 'stderr'    
            errbar=std(y)/sqrt(n);
        case 'mse'
            [MSe, MSej]=produceErrorsWithoutSubjectVariance(y);
            errbar=repmat(sqrt(MSe/n), 1, length(y(1,:)));
        case 'msej'
            [MSe, MSej]=produceErrorsWithoutSubjectVariance(y);
            errbar=sqrt(MSej/n);
        otherwise
          %statements
    end
    % if strcmp(errbarType,'196stderr')
    %     n=length(y(:,1));
    %     df=n-1;
    %     errbar=std(y)/sqrt(n)*tinv(0.975,df);
    % end

    errbar1=reshape(errbar,[length(x),rowY]);
    errbar1=errbar1';
    for j=1:rowY
        errorbar(h(j).XEndPoints,meany1(j,:),errbar1(j,:),'LineStyle','none','Color','k','LineWidth',2);
    end
    %errorbar(h(2).XEndPoints,mean(y),std(y),'LineStyle','none','Color','k','LineWidth',2)
    % simple version
    if scatterFlag==1
        for i=1:columnX
            for j=1:rowY
                scatter(repmat(h(j).XEndPoints(i)-0.1, length(y(:,1)), 1),y(:,i+(j-1)*columnX),'MarkerFaceColor','r','MarkerEdgeColor','k','LineWidth',1);
            end
            
        end
    end
    
end

function  [MSe, MSej]=produceErrorsWithoutSubjectVariance(valueMat)
    %we follow two equations, first MSE=mean(MSEj) second MSS=MSE+c*varSub
    %so varSub=(MSS-MSE)/c 
    %then MSEj=MSWj-varSub; MSE is variance of noise only; MSW  below is variance
    %of noise + subVar
    
    %from the cov matrix varXY=[varX covX; covX varX] varX=varS+varE, covX=VarS
    %mse=mean(varX)-mean(covX)
    %p.264-265 Winer et al. (1991 Statistical principles in experimental design)
    
    varXY=cov(valueMat);
    c=length(varXY(1,:));
    meanvarX=mean(diag(varXY));
    meancovX=(sum(varXY,'all')-sum(diag(varXY)))/(c*(c-1));
    
    MSe=meanvarX-meancovX;
    varsub=meancovX;
    

    MSwj= diag(varXY);%var(valueMat); 
    
    MSej=MSwj-varsub;

end 

function bearing = getbearingA (x,A,i)

         P1A = [A(1) - x(1),0, A(2) - x(2)];
         P1Zplus = [0,0,1]; 
        
         bearingA_raw = getbearingbetween(P1A,P1Zplus);
         bearing = convertangleto360(bearingA_raw - i);
end

function h=ErrorFigure (y, x, errorBarName,scatter,titleName, ylableName,legendName, ylimRange)
    h=barFigure(x,y,errorBarName,scatter)
    %h2=barFigure2(x,y,'stderr',0)
    %title(titleName);
    ylabel([titleName,' (',ylableName, ')']);
    if length(legendName)>0
        legend(legendName);
    end
    ylim(ylimRange);
end

function [h,mdLME, mdLME1,mdl_angle]=ScatterFigure(t12,Xname,Yname, titleName,xlabelName,ylabelName,proj39)
    %mdLME=fitglme(t12, 'abs_pointerror~1+indoor_outdoor+continuous_unfam_fam+(1|participant_number)');

    command=[Yname,'~1+indoor_outdoor+',Xname,'+(1|participant_number)'];
    command1=[Yname,'~1+',Xname,'+(1|participant_number)'];
    mdLME=fitglme(t12, command);
    mdLME1=fitglme(t12, command1);

    %p=mdLME1.Coefficients(2,6);
    p=mdLME1.Coefficients.pValue(2);
    if p<0.001
        symbol='***';
    elseif p<0.01
        symbol='**';
    elseif p<0.05
         symbol='*' ;  
    else
        symbol='ns'  ; 
    end
    h=figure
    [equation, mdl_angle]=proj39.regressionLineWithoutFigure(t12.(Xname),t12.(Yname));
    %title(titleName);
    xlabel(xlabelName);
    ylabel([titleName,' (',ylabelName, ')']);
    indexInside=t12.indoor_outdoor==0;
    indexOutside=t12.indoor_outdoor==1;
    scatter(t12(indexInside,:).(Xname),t12(indexInside,:).(Yname));
    hold on
    scatter(t12(indexOutside,:).(Xname),t12(indexOutside,:).(Yname));
    
    legend('regession line','indoor','outdoor');

    x=2.5;
    y=mdl_angle.Coefficients.Estimate(2).*x+mdl_angle.Coefficients.Estimate(1);
    
    rsq_sl=mdLME1.Rsquared.Ordinary;% which has variance explained by both regressline and subject difference
    rsq_l=mdl_angle.Rsquared.Ordinary; % which has variance explained by regressline
    rsqremovingParticipants= rsq_l/(1-rsq_sl+rsq_l);
    
    r = sqrt(rsqremovingParticipants);
    index=find(equation(:,:)=='r');
    text(x,y, [equation(1:index),'=',num2str(r), symbol]);
    
    printEnvironmentName=1;
    if strcmp(Xname,'envrionmentFamilarity')
        printEnvironmentName=0;
    end
    if printEnvironmentName==1
        [meanM,buildingNames]=proj39.getMeanORP1(t12.(Xname),t12.('Enviornmentname'));
        for i=1:length(meanM)
            %index=strcmp(T12.Enviornmentname,buildingNames{i});
            %T12.envrionmentFamilarityMean(index)=meanM(i);
            text(meanM(i),0,buildingNames{i}(1:2));
        end
    end
end
