function [tau_opt, R_tau, lambda] = optimal_tau_cal(k, gamma)
% This function optimizes the value of tau to obtained the maximum common
% throughput, according to the algorithm in Table I in the paper

R_min = 0; R_max = 1e2;
lambda = 1 ./ log10(1+gamma); % initialize the value of lambda

% calculate the optimal values of z and tau based on Proposition 4.1
[R_tau, tau_opt] = tau_opt_cal(k, lambda, gamma);
    
delta = 1e-3;
tau = 0.01;
while R_max - R_min > delta    
    flag_R_aver = 0;
    while flag_R_aver == 0
        R_aver = mean([R_min, R_max]);
        g_lambda = - sum(lambda .* R_tau) + sum(lambda) * R_aver; % compute g(lambda)
        
        if g_lambda > 0 % R_aver is infeasible, decrease R_aver
            R_max = R_aver;
            if R_max - R_min <= delta
                flag_R_aver =1;
            end
        elseif g_lambda <=0
            % update the value of lambda by using the gradient descent algorithm
            [tau, lambda, R_tau, tau_opt] = gradient_descent(tau, lambda, k, R_tau, R_aver, gamma);
            g_lambda = - sum(lambda .* R_tau) + sum(lambda) * R_aver; % compute g(lambda)
            if g_lambda <= 0
                R_min = R_aver;
            end  
            flag_R_aver = 1;
        end
        
    end
end