function [t s] = TW_trace_ratio_threshold(p,n,beta,alpha)
%function [t s] = TW_trace_ratio_threshold(p,n,beta,alpha)
% Let L_j be the eigenvalues of the sample covariance matrix of multivariate
% Gaussian obserations with n samples in p dimensions, 
%
% INPUT:  p = dimension, n = number of samples, 
%         beta = 1/2 real/complex valued data, 
%         alpha = desired false alarm / right tail probability
%
% OUTPUT: t = An approximate value such that Pr[ L_1/(1/p sum(L_j) )> t ] ~  alpha
%         s = Pr[ (L_1/(1/p sum(L_j) - mu)/sigma > s] ~ alpha
%
%
% Threshold s is found by interpolation of precomumputed values of the
% Tracy-Widom distribution, and its second derivative. 
% These were generated using code by Prof. Folkmar Bornemann, based on his paper
% On the Numerical Evaluation of Distributions in Random Matrix Theory, 2010. 
%

if beta==1
    load TW_beta1.mat
else
    load TW_beta2.mat
end

% mu_np and sigma_np - centering and scaling
[mu_np sigma_np] = KN_mu_sigma(n,p,beta);   mu_np = mu_np ./ n; sigma_np = sigma_np ./ n; 

% Eq. 1.7 in paper [1]
U_cdf_complementary = 1-TW_s +1/beta/n/p * (mu_np/sigma_np)^2 * TW_s_tag_tag; 

% next find index where probability is close to alpha
[val idx] = min(abs(U_cdf_complementary-alpha) ); 

t = mu_np + sigma_np * x(idx); 
s = x(idx); 
