clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           Discrimination optimization
%
% Calculation of optimal estimates s1, s2, signal to noise ratio (S/N),
% and S/N gain percentage (G%) assuming s1 = s2 =1, according to Hellström's (1986, 1989) 
% discrimination optimization model. For detials see...
%
% Patching, G.R., Englund, M.P., & Hellström, Å. (2012). Time- and
% space-order effects in timed discrimination of brightness and size of paired
% visual stimuli. Journal of Experimental Psychology: Human Perception and Performance, in press.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%series 1 estimates
V1=1;    %stdev(Psi1)
V2=0.5;  %stdev(Psi2)
V3=0.75; %std(ReL1)
V4=0.75; %std(ReL2)

R4=0.5;  %corr(Psi1,Psi2);
R5=-0.3; %corr(Psi1,ReL1);
R6=-0.3; %corr(Psi2,ReL2);
R7=-0.2; %corr(Psi1,ReL2);
R8=-0.2; %corr(Psi2,ReL1);
R9=0.5;  %corr(ReL1,ReL2);

A=sqrt(1 + V3.^2 - 2 .* R5 .* V3);
K=sqrt(V2.^2 + V4.^2 - 2 .* R6 .* V2 .* V4) ./ A;
C=sqrt(V3.^2 + V4.^2 - 2 .* R9 .* V3 .* V4) ./ A;

R1=(R4 .* V2 - R7 .* V4 - R8 .* V2 .* V3 + R9 .* V3 .* V4) ./ (K .* A.^2);
R2=(R5 .* V3 - R7 .* V4 - V3.^2 + R9 .* V3 .* V4) ./ (C .* A.^2);
R3=(R8 .* V2 .* V3 - R6 .* V2 .* V4 - R9 .* V3 .* V4 + V4.^2) ./ (C .* K .* A.^2);

S1=C .* (-K .* R3.^2 + K + R1 - R2 .*R3)./(K .* R1 .* R3 - K .* R2 - R1 .* R2 + R3);
S2=C .* (K .* R1 -K .* R2 .* R3 - R2.^2 + 1)./(K .* (K .* R1 .* R3 - K .* R2 - R1 .* R2 + R3));

Squot= S1 ./ S2;

yden=sqrt(S1 .^2 + K.^2 .* S2.^2 - 2 .* K .* R1 .* S1 .* S2  + C.^2 + 2 .* C .* (R2 .* S1 - K .* R3 .* S2));
Y=(S1 + S2)./yden;

S1b=1;
S2b=1;

ydenb=sqrt(S1b .^2 + K.^2 .* S2b.^2 - 2 .* K .* R1 .* S1b .* S2b  + C.^2 + 2 .* C .* (R2 .* S1b - K .* R3 .* S2b));
Yb=(S1b + S2b)./ydenb;

SNROPT=Y ./ A;
SNRCOMP=Yb ./ A;

EFFPROC=100 .* SNROPT / SNRCOMP;
G=EFFPROC-100;

fprintf('Stdev(PSI2)  ReL    S1     S2     S1/S2   S/N     G \n');
fprintf('  %1.1f        %1.1f   %1.2f   %1.2f    %1.2f    %1.2f   %2.2f \n',V2,V3,S1,S2,S1/S2,Y,G);

CRIT123=R1.^2 + R2.^2 + R3.^2 - 2 .* R1 .* R2 .* R3;
CRIT458=R4.^2 + R5.^2 + R8.^2 - 2 .* R4 .* R4 .* R8;
CRIT467=R4.^2 + R6.^2 + R7.^2 - 2 .* R4 .* R6 .* R7;
CRIT579=R5.^2 + R7.^2 + R9.^2 - 2 .* R5 .* R7 .* R9;
CRIT689=R6.^2 + R8.^2 + R9.^2 - 2 .* R6 .* R8 .* R9;

T1=R7 .* V4 ./ V1;
T2=R8 .* V3 ./ V2;
S1PRIM = S1 - (1 - S2) .* T2;
S2PRIM = S2 - (1 - S1) .* T1;
SPRIMQ = S1PRIM ./ S2PRIM;
A1FACT = (1-S1) .* (1-T2);
A2FACT = (1-S2) .* (1-T1);

if (CRIT123 < 0) || (CRIT123>1) disp('ERROR')
elseif (CRIT458 < 0) || (CRIT458>1) disp('ERROR')
elseif (CRIT467 < 0) || (CRIT467>1) disp('ERROR')
elseif (CRIT579 < 0) || (CRIT579>1) disp('ERROR')
elseif (CRIT689 < 0) || (CRIT689>1) disp('ERROR')    
end