function [derfresult,x]=derf(x);
derfresult=[];
persistent erfcs first firstCall nterf sqeps sqrtpi xbig y ; if isempty(firstCall),firstCall=1;end;
;
if isempty(nterf), nterf=0; end;
%***BEGIN PROLOGUE DERF
%***PURPOSE Compute the error function.
%***LIBRARY SLATEC (FNLIB)
%***CATEGORY C8A, L5A1E
%***TYPE doubleprecision (ERF-S, DERF-D)
%***KEYWORDS ERF, ERROR FUNCTION, FNLIB, SPECIAL FUNCTIONS
%***AUTHOR Fullerton, W., (LANL)
%***DESCRIPTION
%
% DERF(X) calculates the doubleprecision error function for double
% precision argument X.
%
% Series for ERF on the interval 0. to 1.00000E+00
% with weighted error 1.28E-32
% log weighted error 31.89
% significant figures required 31.05
% decimal places required 32.55
%
%***REFERENCES (NONE)
%***ROUTINES CALLED D1MACH, DCSEVL, DERFC, INITDS
%***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)
% 900727 Added EXTERNAL statement. (WRB)
% 920618 Removed space from variable name. (RWC, WRB)
%***end PROLOGUE DERF
if isempty(erfcs), erfcs=zeros(1,21); end;
if isempty(sqeps), sqeps=0; end;
if isempty(sqrtpi), sqrtpi=0; end;
if isempty(xbig), xbig=0; end;
if isempty(y), y=0; end;
if isempty(first), first=false; end;
if firstCall, erfcs(1)=[-.49046121234691808039984544033376d-1]; end;
if firstCall, erfcs(2)=[-.14226120510371364237824741899631d+0]; end;
if firstCall, erfcs(3)=[+.10035582187599795575754676712933d-1]; end;
if firstCall, erfcs(4)=[-.57687646997674847650827025509167d-3]; end;
if firstCall, erfcs(5)=[+.27419931252196061034422160791471d-4]; end;
if firstCall, erfcs(6)=[-.11043175507344507604135381295905d-5]; end;
if firstCall, erfcs(7)=[+.38488755420345036949961311498174d-7]; end;
if firstCall, erfcs(8)=[-.11808582533875466969631751801581d-8]; end;
if firstCall, erfcs(9)=[+.32334215826050909646402930953354d-10]; end;
if firstCall, erfcs(10)=[-.79910159470045487581607374708595d-12]; end;
if firstCall, erfcs(11)=[+.17990725113961455611967245486634d-13]; end;
if firstCall, erfcs(12)=[-.37186354878186926382316828209493d-15]; end;
if firstCall, erfcs(13)=[+.71035990037142529711689908394666d-17]; end;
if firstCall, erfcs(14)=[-.12612455119155225832495424853333d-18]; end;
if firstCall, erfcs(15)=[+.20916406941769294369170500266666d-20]; end;
if firstCall, erfcs(16)=[-.32539731029314072982364160000000d-22]; end;
if firstCall, erfcs(17)=[+.47668672097976748332373333333333d-24]; end;
if firstCall, erfcs(18)=[-.65980120782851343155199999999999d-26]; end;
if firstCall, erfcs(19)=[+.86550114699637626197333333333333d-28]; end;
if firstCall, erfcs(20)=[-.10788925177498064213333333333333d-29]; end;
if firstCall, erfcs(21)=[+.12811883993017002666666666666666d-31]; end;
if firstCall, sqrtpi=[1.77245385090551602729816748334115d0]; end;
if firstCall, first=[true]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT DERF
if( first )
[nterf ,erfcs]=initds(erfcs,21,0.1.*real(d1mach(3)));
xbig = sqrt(-log(sqrtpi.*d1mach(3)));
sqeps = sqrt(2.0d0.*d1mach(3));
end;
first = false;
%
y = abs(x);
if( y>1.0d0 )
%
% ERF(X) = 1.0 - ERFC(X) FOR ABS(X) .GT. 1.0
%
if( y<=xbig )
derfresult = (abs(1.0d0-derfc(y)).*sign(x));
end;
if( y>xbig )
derfresult = (abs(1.0d0).*sign(x));
end;
else;
%
% ERF(X) = 1.0 - ERFC(X) FOR -1.0 .LE. X .LE. 1.0
%
if( y<=sqeps )
derfresult = 2.0d0.*x.*x./sqrtpi;
end;
if( y>sqeps )
derfresult = x.*(1.0d0+dcsevl(2.0d0.*x.*x-1.0d0,erfcs,nterf));
end;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK DERKF
