[gamiresult,a,x]=gami(a,x);

function [gamiresult,a,x]=gami(a,x);
gamiresult=[];
persistent factor gami ; 

if isempty(factor), factor=0; end;
if isempty(gamiresult), gamiresult=0; end;
%***BEGIN PROLOGUE  GAMI
%***PURPOSE  Evaluate the incomplete Gamma function.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C7E
%***TYPE      SINGLE PRECISION (GAMI-S, DGAMI-D)
%***KEYWORDS  FNLIB, INCOMPLETE GAMMA FUNCTION, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% Evaluate the incomplete gamma function defined by
%
% GAMI = integral from T = 0 to X of EXP(-T) * T**(A-1.0) .
%
% GAMI 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 GAMI is very large or very small, because logarithmic variables
% are used.  GAMI, A, and X are single precision.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  ALNGAM, GAMIT, 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  GAMI
%***FIRST EXECUTABLE STATEMENT  GAMI
if( a<=0.0 )
xermsg('SLATEC','GAMI','A MUST BE GT ZERO',1,2);
end;
if( x<0.0 )
xermsg('SLATEC','GAMI','X MUST BE GE ZERO',2,2);
end;
%
gamiresult = 0.0;
if( x==0.0 )
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(alngam(a)+a.*log(x));
%
gamiresult = factor.*gamit(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 GAMIT