function [IRmax, vector] = thresholdIR(input)

global Po
global Pc
Pf = 1 - Po - Pc;

vector = zeros(1,input);


for X=1:input
    pmf = 0;
    
    for k = 0:X
        for t = 0:(input-X)
            coef = factorial(input) / (factorial(k) * factorial(t) * factorial(input-t-k)) ;
            prob = (Pc^k) * (Pf^(input-t-k)) * (Po^t) ;
            tempsum = coef*prob;
            pmf = pmf + tempsum;
        end
    end
    vector(1,X) = pmf;
end

    
vector = floor(1000*vector);
 temp = vector(1);
 
 for k=2:numel(vector)          
     if temp >=vector(k)
         IRmax = k-1;
         break
     end    
     temp = vector(k);
 end


% [~, IRmax] = max(vector);

end