%% Probably better version that directly mixes radians

function [AccuMix, BiasMix]=squircleBehave2(s,A,CardPrez,behB,makefig)
  
%% Settings/prep
n=16;           % number of equidistant rotations
testd=6.43;        % +/- displacement of test stims in deg (ask Juan)?
% % developer settings:
% s=0.2;         % noise (in memory/decision-making)
% A=0.2;          % Key Parameter (squircle) 0: all circle; 1: all square 
% behB input with behavioural biases
% CardPrez=0.6;   % Reduce noise near cardinal by this factor (1: off)
% makefig=1;

%% Get ori vectors (where oriC is "Circle")
[oriC, oriP, oriN]=orivectors(n,testd); 

%% Also prepare categorical vector (where oriS is "Square") 
oriS1=[repmat(45,1,4) repmat(3*45,1,4) repmat(5*45,1,4) repmat(7*45,1,4)]; % cat
oriS1=deg2rad(oriS1);
oriS2=[0 0 repmat(90,1,4) repmat(2*90,1,4) repmat(3*90,1,4) 360 360];      % anti-cat
oriS2=deg2rad(oriS2);
 
% mix oris directly
if A>=0
    oriM=A*oriS1+(1-A)*oriC;
elseif A<0
    fA=-A; % need to flip since its negative
    oriM=fA*oriS2+(1-fA)*oriC;
end

%% Calculate angular distances
% between mixed ori and test ori:
dTp=angdiff(oriM,oriP);
dTn=angdiff(oriM,oriN);

%% (optional) better precision at cardinal   
sW=repmat([CardPrez 1 1 CardPrez],1,4);
sW=sW/mean(sW); % ensure that the mean of the vector is 1
s=s*sW;
% (this could be decomposed further into horizontal and vertical,
% but probably better to keep it simple for now)

%% Predict choices under either model
CP_P=sigmoid(dTp,s,0); % p of sensing "ccw" when test was "ccw" (correct)
CP_N=sigmoid(dTn,s,0); % p of sensing "ccw" when test was "cw"  (error)


%% Accuracy (by objective measures as used in analysis) should be
AccuMix=(CP_P+(1-CP_N))/2; % for circle this is symmetric

%% Bias should simply be the average of CPs (CW/CCW)
BiasMix=(CP_P+CP_N)/2;

%% Plots
% note: I'll put 0° at the right (standard circle)
if makefig
%     subplot(2,2,1); quickplot(oriC,AccuMix,AccuMix,'Accuracy');
%     subplot(2,2,2);
    quickplotcompB(oriC,BiasMix,behB,'Bias');

%     plotori=[oriM oriM(1)]; % to draw complete squircle
%     subplot(2,2,[3 4]); polarplot(plotori,ones(1,n+1),'o-');
    title(['Squircle A=' num2str(A)]); rlim([0 1.5]);   
end

end