function p= prot2(n, t, param, k)
%
% function p= prot2(n, t, param, k)
%
% returns Eq. 10 with param= [a b] and n a vector
% written by V Shahrezaei & PS Swain

p= nan(size(n));
r= [0:k];
v= exp(-r*t+gammaln(k+1)-gammaln(r+1)-gammaln(k-r+1)+(k-r)*log(1-exp(-t)));
for i= 1:length(n)
    if n(i) == 0
        p(i)= PproteinSt(0, t, param)*(1-exp(-t))^k;
    elseif n(i) > 0
        p(i)= sum(PproteinSt(n(i)-r, t, param).*v);
    else
        p(i)= 0;
    end;
end;



%%%%%%%%%%
function p= PproteinSt(n, t, param)
% function p= PproteinSt(n, t, param)
%
% returns prob at time t with param= [a b] and n a vector

a= param(1);
b= param(2);
te= exp(-t);
x= (1+b)*te/(1+b*te);
p= zeros(size(n));
fac= exp(-gammaln(n+1)+gammaln(n+a)-gammaln(a))...
    .*(b/(1+b)).^n*((1+b*te)/(1+b))^a;

for i= 1:length(n)
    if n(i) > 0
        
        % calculate 2F1(-n, -a, 1-a-n; x)
        z= 0;
        for m= 0:n(i)
            if m == 0
                delta= 1;
            else
                delta= delta*x*(-n(i)+m-1)*(-a+m-1)/m/(-a-n(i)+m);
            end;
            z= z + delta;
        end;
        p(i)= z*fac(i);
        
    elseif n(i) == 0
        
        p(i)= ((1+b*te)/(1+b))^a;
        
    else
        
        % negative n
        p(i)= 0;
        
    end;
end;
