// ====================================================================
//  gammainc : Incomplete gamma function.
//    P = gammainc(X,A) evaluates the incomplete gamma function for
//    X and A.  X and A must be real matrix (or scalar). X must be non-negative. 
//    A must also be non-negative and non-equal to zero.
// 
//    The incomplete gamma function is defined as:
// 
//  P(x,a) = (1 / gamma(a)) * integral from 0 to x of exp(-t) t^(a-1) dt
//
//  Algorithm based on :
//          NUMERICAL RECIPES
//   The Art of Scientific Computing
//  Cambridge University Press - 1986
//
//
//  Mathieu OLIVIER
//  Météo-France CNRM/GMEI/LISA 2010
//
//  This software is governed by the CeCILL  license under French law and
//  abiding by the rules of distribution of free software.  You can  use,
//  modify and/ or redistribute the software under the terms of the CeCILL
//  license as circulated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.
// ====================================================================
function P = gammainc(X,A)

	if (argn(2)<2) then
		error('[gammainc] Not enough input arguments.');
	elseif (argn(2)>3) then
		error('[gammainc] Too many input arguments.');
	end
	if isempty(X) then
		error('[gammainc] Input argument ''X'' is empty.');
	elseif isempty(A) then
		error('[gammainc] Input argument ''A'' is empty.');
	elseif (type(X)~=1) then
		error(sprintf(...
		'[gammainc] Undefined function or method for input argument ''X'' of type ''%s''.'...
		,typeof(X)));
	elseif (type(A)~=1) then
		error(sprintf(...
		'[gammainc] Undefined function or method for input argument ''A'' of type ''%s''.'...
		,typeof(A)));
	elseif (~isreal(X)) then
		error('[gammainc] Input argument ''X'' must not be complex.');
	elseif (~isreal(A)) then
		error('[gammainc] Input argument ''A'' must not be complex.');
	elseif (~isempty(find(X<0))) then
		error('[gammainc] At least one element of X is less than zero.');
	elseif (~isempty(find(A<=0))) then
		error('[gammainc] At least one element of A is less or equal to zero.');
	elseif ((sum(length(A))~=1)&(length(X)~=length(A))) then
		error('[gammainc] Input arguments dimensions mismatch.');
	end
	if ((sum(length(A))~=1)&(size(X,1)~=size(A,1))) then
		A = A';
	end
	P = zeros(size(X,1),size(X,2));
	// -----------------------------------------------------------------------------
	// P with 0<X<A+1
	Xtmp = X((X<A+1)&(X>0));
	Psec = zeros(size(Xtmp,1),size(Xtmp,2));
	itmax = 100;
	if sum(length(A))==1 then
		Atmp = A.*ones(size(Xtmp,1),size(Xtmp,2));
	else
		Atmp = A((X<A+1)&(X>0));
	end
	Gln = gammaln(Atmp);
	AP = Atmp;
	SUM = 1../Atmp;
	DEL = 1../Atmp;
	for i = 1:length(Xtmp),
		n=1;
		while n<itmax,
			AP(i) = AP(i)+1;
			DEL(i) = DEL(i)*Xtmp(i)/AP(i);
			SUM(i)=SUM(i)+DEL(i);
			if abs(DEL(i))<abs(SUM(i))*epsilon(SUM(i)) then
				break;
			end
			n = n+1;
		end
		Psec(i) = SUM(i)*exp(-Xtmp(i)+Atmp(i)*log(Xtmp(i))-Gln(i));
	end
	P((X<A+1)&(X>0))=Psec;
	clear Psec; clear Ap; clear SUM; clear DEL;
	// -----------------------------------------------------------------------------
	// P with X>A+1
	Xtmp = X(X>=A+1);
	Psec = zeros(size(Xtmp,1),size(Xtmp,2));
	itmax = 100;
	if sum(length(A))==1 then
		Atmp = A.*ones(size(Xtmp,1),size(Xtmp,2));
	else
		Atmp = A(X>=A+1);
	end
	Gln = gammaln(Atmp);
	for i = 1:length(Xtmp),
		Gold = 0;
		A0 = 1;
		A1 = Xtmp(i);
		B0 = 0;
		B1 = 1;
		FAC = 1
		n=1;
		while n<itmax,
			ana = n-Atmp(i);
			A0 = (A1+A0*ana)*FAC;
			B0 = (B1+B0*ana)*FAC;
			ANF = n*FAC;
			A1 = Xtmp(i)*A0+ANF*A1;
			B1 = Xtmp(i)*B0+ANF*B1;
			if A1~=0 then
				FAC = 1/A1;
				G = B1*FAC;
				if abs((G-Gold)/G)<epsilon(G) then
					break;
				else
					Gold = G;
				end
			end
			n = n+1;
		end
		Psec(i) = exp(-Xtmp(i)+Atmp(i)*log(Xtmp(i))-Gln(i))*Gold;
	end
	P(X>=A+1)=1-Psec;

endfunction

// ====================================================================


function eps = epsilon(A)
	if length(A)~=1 then error(84,1);end
	eps = min([abs(A-nearfloat('pred',A));abs(A-nearfloat('succ',A))]);
endfunction