function pk=hypergeop(N,n,r,k)
% format function pk=hypergeop(N,n,r,k)
% input N=  total number of objects
%       n=  number of objects selected
%       r=  number of objects of type r
%       k=  value for which p(k) will be calculated
%   output: pk=p(sampling k objects of type r)
% See Larsen & Marx (2001) Theorem 3.3.1, page 128, or
% Larsen & Marx (2006, p. 139) Introduction to Mathematical Statistics,
% 4th edition
% Theorem 3.2.2,
% Written by E. Gallagher, EEOS, UMASS/Boston, revised 10/19/10
if length(k)==1
    if N<50
        pk=nchoosek(r,k).*nchoosek(N-r,n-k)./nchoosek(N,n);
        return
    else
        o=1;
        % will need to calculate with n!=exp(ln(gamma(n+1)) below
    end
elseif length(N(:))*length(n(:))*length(r(:))==1
    % This section will calculate all of the probabilities for vector k
    [rk,ck]=size(k);
    N=repmat(N,rk,ck);n=repmat(n,rk,ck);
    r=repmat(r,rk,ck);o=ones(rk,ck);
    % Use fact that n!=exp(gammaln(n+1))
else
    error('N, n, r must have size(k) or length 1')
end
lnnf=gammaln(n+o);lnNf=gammaln(N+o);
lnnmkf=gammaln(n-k+o);lnrf=gammaln(r+o);
lnwf=gammaln(N-r+o);lnrmkf=gammaln(r-k+o);
lnwmnmkf=gammaln(N-r-(n-k)+o);lnNmnf=gammaln(N-n+o);
lnkf=gammaln(k+o);
    % note that ln(a/b) is the same as ln(a)-ln(b)
     pk=exp(...
         (lnrf-lnrmkf-lnkf)+(lnwf-lnwmnmkf-lnnmkf)-...
         (lnNf-lnNmnf-lnnf)      );