[dbetaresult,a,b]=dbeta(a,b);

function [dbetaresult,a,b]=dbeta(a,b);
dbetaresult=[];
persistent alnsml first firstCall xmax xmin ; if isempty(firstCall),firstCall=1;end; 

;
%***BEGIN PROLOGUE  DBETA
%***PURPOSE  Compute the complete Beta function.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C7B
%***TYPE      doubleprecision (BETA-S, DBETA-D, CBETA-C)
%***KEYWORDS  COMPLETE BETA FUNCTION, FNLIB, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% DBETA(A,B) calculates the doubleprecision complete beta function
% for doubleprecision arguments A and B.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  D1MACH, DGAMLM, DGAMMA, DLBETA, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   770601  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)
%   900727  Added EXTERNAL statement.  (WRB)
%***end PROLOGUE  DBETA
if isempty(alnsml), alnsml=0; end;
if isempty(xmax), xmax=0; end;
if isempty(xmin), xmin=0; end;
if isempty(first), first=false; end;
if firstCall,   first=[true];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  DBETA
if( first )
[xmin,xmax]=dgamlm(xmin,xmax);
alnsml = log(d1mach(1));
end;
first = false;
%
if( a<=0.0D0 || b<=0.0D0 )
xermsg('SLATEC','DBETA','BOTH ARGUMENTS MUST BE GT 0',2,2);
end;
%
if( a+b<xmax )
dbetaresult = dgamma(a).*dgamma(b)./dgamma(a+b);
end;
if( a+b<xmax )
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',b); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',a); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
[dbetaresult ,a,b]=dlbeta(a,b);
if( dbetaresult<alnsml )
%
dbetaresult = 0.0d0;
xermsg('SLATEC','DBETA','A AND/OR B SO BIG BETA UNDERFLOWS',1,1);
else;
dbetaresult = exp(dbetaresult);
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',b); 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',b); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',a); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK DBETAI