[gamitresult,a,x]=gamit(a,x);

function [gamitresult,a,x]=gamit(a,x);
gamitresult=[];
persistent aeps ainta algap1 alneps alng alx bot first firstCall h sga sgngam sqeps t ; if isempty(firstCall),firstCall=1;end; 

;
if isempty(aeps), aeps=0; end;
if isempty(ainta), ainta=0; end;
if isempty(algap1), algap1=0; end;
if isempty(alneps), alneps=0; end;
if isempty(alng), alng=0; end;
if isempty(alx), alx=0; end;
if isempty(bot), bot=0; end;
if isempty(h), h=0; end;
if isempty(sga), sga=0; end;
if isempty(sgngam), sgngam=0; end;
if isempty(sqeps), sqeps=0; end;
if isempty(t), t=0; end;
%***BEGIN PROLOGUE  GAMIT
%***PURPOSE  Calculate Tricomi's form of the incomplete Gamma function.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C7E
%***TYPE      SINGLE PRECISION (GAMIT-S, DGAMIT-D)
%***KEYWORDS  COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB,
%             SPECIAL FUNCTIONS, TRICOMI
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
%   Evaluate Tricomi's incomplete gamma function defined by
%
%   GAMIT = X**(-A)/GAMMA(A) * integral from 0 to X of EXP(-T) *
%             T**(A-1.)
%
%   for A .GT. 0.0 and by analytic continuation for A .LE. 0.0.
%   GAMMA(X) is the complete gamma function of X.
%
%   GAMIT is evaluated for arbitrary real values of A and for non-
%   negative values of X (even though GAMIT is defined for X .LT.
%   0.0), except that for X = 0 and A .LE. 0.0, GAMIT is infinite,
%   which is a fatal error.
%
%   The function and both arguments are REAL.
%
%   A slight deterioration of 2 or 3 digits accuracy will occur when
%   GAMIT is very large or very small in absolute value, because log-
%   arithmic variables are used.  Also, if the parameter  A  is very
%   close to a negative integer (but not a negative integer), there is
%   a loss of accuracy, which is reported if the result is less than
%   half machine precision.
%
%***REFERENCES  W. Gautschi, A computational procedure for incomplete
%                 gamma functions, ACM Transactions on Mathematical
%                 Software 5, 4 (December 1979), pp. 466-481.
%               W. Gautschi, Incomplete gamma functions, Algorithm 542,
%                 ACM Transactions on Mathematical Software 5, 4
%                 (December 1979), pp. 482-489.
%***ROUTINES CALLED  ALGAMS, ALNGAM, GAMR, R1MACH, R9GMIT, R9LGIC,
%                    R9LGIT, XERCLR, 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)
%   920528  DESCRIPTION and REFERENCES sections revised.  (WRB)
%***end PROLOGUE  GAMIT
if isempty(first), first=false; end;
if firstCall,   first=[true];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  GAMIT
if( first )
alneps = -log(r1mach(3));
sqeps = sqrt(r1mach(4));
bot = log(r1mach(1));
end;
first = false;
%
if( x<0.0 )
xermsg('SLATEC','GAMIT','X IS NEGATIVE',2,2);
end;
%
if( x~=0.0 )
alx = log(x);
end;
sga = 1.0;
if( a~=0.0 )
sga = (abs(1.0).*sign(a));
end;
ainta = fix(a+0.5.*sga);
aeps = a - ainta;
%
if( x<=0.0 )
gamitresult = 0.0;
if( ainta>0.0 || aeps~=0.0 )
[ gamitresult ]=gamr(a+1.0);
end;
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;
%
elseif( x<=1.0 ) ;
if( a>=(-0.5) || aeps~=0.0 )
[dumvar1,algap1,sgngam]=algams(a+1.0,algap1,sgngam);
end;
[gamitresult ,a,x,algap1,sgngam,alx]=r9gmit(a,x,algap1,sgngam,alx);
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;
%
elseif( a<x ) ;
%
[alng ,a,x,alx]=r9lgic(a,x,alx);
%
% EVALUATE GAMIT IN TERMS OF LOG(GAMIC(A,X))
%
h = 1.0;
if( aeps~=0.0 || ainta>0.0 )
[dumvar1,algap1,sgngam]=algams(a+1.0,algap1,sgngam);
t = log(abs(a)) + alng - algap1;
if( t>alneps )
%
t = t - a.*alx;
if( t<bot )
xerclr;
end;
gamitresult = -sga.*sgngam.*exp(t);
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;
else;
if( t>(-alneps) )
h = 1.0 - sga.*sgngam.*exp(t);
end;
if( abs(h)<=sqeps )
xerclr;
xermsg('SLATEC','GAMIT','RESULT LT HALF PRECISION',1,1);
end;
end;
end;
%
t = -a.*alx + log(abs(h));
if( t<bot )
xerclr;
end;
gamitresult = (abs(exp(t)).*sign(h));
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;
else;
[t ,a,x]=r9lgit(a,x,alngam(a+1.0));
if( t<bot )
xerclr;
end;
gamitresult = exp(t);
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;
%
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 GAMLIM