function p= prot3(n, param, hacc)
% function p= prot3(n, param, hacc)
%
% returns Eq. 18 with param= [a b ka0 ka1] and hacc (= 1.0e-7)
% n can be a vector
% written by V Shahrezaei & PS Swain

if nargin == 2
    hacc= 1.0e-7;
end;

a= param(1);
b= param(2);
ka0= param(3);
ka1= param(4);

ga= ka0+ka1;
s= sqrt(-4*a*ka0+(a+ga)^2);
al= (a+ga+s)/2;
be= (a+ga-s)/2;

for i= 1:length(n)
    
    lp(i)= gammaln(n(i)+1)-gammaln(al+n(i))-gammaln(be+n(i))...
        -gammaln(ga)+gammaln(al)+gammaln(be)+gammaln(ga+n(i))...
        -n(i)*log(b)+(n(i)+al)*log(1+b)...
        -hypergeometric2f1ln(n(i)+al, ga-be, n(i)+ga, b/(1+b), hacc);
    
end;

p= exp(-lp);


%%%%%%%%%%%%%%%%%%%%%%%%%%

function lz= hypergeometric2f1ln(a, b, c, x, acc)
% function lz= hypergeometric2f1ln(a, b, c, x, acc)
%
% note if vectors, a, b, c, and x must all be the same size
% returns log(function)


if max(x) > 1
    disp(' x must be less than one!!!');
    return;
end;
n= 1.0e5;

% find maximum term in series
lx= log(x);
lxm= mean(lx);
am= mean(a);
bm= mean(b);
cm= mean(c);
if meqn(0, am, bm, cm, lxm)*meqn(n, am, bm, cm, lxm) > 0
    rmax= log(a.*b./c)+lx;
    mmax= 1;
else
    mmax=round(fzero(@meqn,[0 n], [], am, bm, cm, lxm));
    rmax= gammaln(a+mmax)+gammaln(b+mmax)+gammaln(c)-gammaln(a)...
        -gammaln(b)-gammaln(c+mmax)-gammaln(1+mmax)+mmax*lx;
end;

% sum a_n/a_max
el= exp(-rmax);
s= el + 1; % as a_0= 1
if mmax == 1
    % first coeffient in series largest
    delta= 1;
    for m= 2:n
        delta= delta.*x.*(a+m-1).*(b+m-1)./m./(c+m-1);
        s= s + delta;
        if isempty(find(abs(delta./s) > acc, 1))
            break;
        end;
    end;
else
    % middle coefficient largest
    delta= 1;
    for m= mmax:-1:2
        delta= delta./x./(a+m-1)./(b+m-1).*m.*(c+m-1);
        s= s + delta;
        if isempty(find(abs(delta./s) > acc, 1))
            break;
        end;
    end;
    delta= 1;
    for m= mmax+1:n
        delta= delta.*x.*(a+m-1).*(b+m-1)./m./(c+m-1);
        s= s + delta;
        if isempty(find(abs(delta./s) > acc, 1))
            break;
        end;
    end;
end;

lz= rmax + log(s);



%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y= meqn(m, a, b, c, lx)
y= lx-psi(m+1)+psi(a+m)+psi(b+m)-psi(c+m);


