[dgamiresult,a,x]=dgami(a,x);

function [dgamiresult,a,x]=dgami(a,x);
dgamiresult=[];
persistent factor ; 

;
%***BEGIN PROLOGUE  DGAMI
%***PURPOSE  Evaluate the incomplete Gamma function.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C7E
%***TYPE      doubleprecision (GAMI-S, DGAMI-D)
%***KEYWORDS  FNLIB, INCOMPLETE GAMMA FUNCTION, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% Evaluate the incomplete gamma function defined by
%
% DGAMI = integral from T = 0 to X of EXP(-T) * T**(A-1.0) .
%
% DGAMI is evaluated for positive values of A and non-negative values
% of X.  A slight deterioration of 2 or 3 digits accuracy will occur
% when DGAMI is very large or very small, because logarithmic variables
% are used.  The function and both arguments are doubleprecision.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  DGAMIT, DLNGAM, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   770701  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%***end PROLOGUE  DGAMI
if isempty(factor), factor=0; end;
%***FIRST EXECUTABLE STATEMENT  DGAMI
if( a<=0.0D0 )
xermsg('SLATEC','DGAMI','A MUST BE GT ZERO',1,2);
end;
if( x<0.0D0 )
xermsg('SLATEC','DGAMI','X MUST BE GE ZERO',2,2);
end;
%
dgamiresult = 0.0d0;
if( x==0.0d0 )
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',x); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',a); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
% THE ONLY ERROR POSSIBLE IN THE EXPRESSION BELOW IS A FATAL OVERFLOW.
factor = exp(dlngam(a)+a.*log(x));
%
dgamiresult = factor.*dgamit(a,x);
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',x); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',a); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK DGAMIT