function [ATS] = computeATS(P, L1, n1, alpha)
    A = 1 - (1 - P)^L1;
    ATS = (n1/P)*(1 - alpha*(1 - alpha)*A^2)...
        /(A^2*(1 + alpha*(1 - alpha)*(A - 2)));
endfunction

//------------------------------------------------------------------------------
function [k, L, n, ATS] = optGR(delta1, tau)
//------------------------------------------------------------------------------

m1 = tau
n1 = 1
k1 = 0.01
L1 = 1

while n1 >= m1
    k1 = 0.01
	while k1 < 3
        P0 = 1 - cdfnor(k1) + cdfnor(-k1);
        P1 = -cdfnor(k1 - delta1*sqrt(n1)) + cdfnor(-k1 - delta1*sqrt(n1))
        alpha0 = (1 - cdfnor(k1))/P0 ;
        alpha1 = (1 - cdfnor(k1 - delta1*sqrt(n1)))/P1
        L1 = 1
        while L1 < 20000
            ATS0 = computeATS(P0, L1, n1, alpha0)
            ATS1 = computeATS(P1, L1, n1, alpha1)
            if ATS1 > m1
                L1 = L1 + 1
            elseif ATS0 < tau
                k1 = k1 + 0.01
                L1 = 30000
            else
                m1 = ATS1
                m0 = ATS0
                L1 = L1 + 1
			end
		else
            k1 = k1 + 0.01
		end
	else
        n1 = n1 + 1
    end
end
k = k1
L = L1
n = n1
ATS = m1

endfunction

