#include "fpreproc/library_names.h"
#include "fpreproc/f77_dcomplx.h"
#include "fpreproc/fp_const.h"
*DECK DERF
      REAL*8 FUNCTION DERF (X)
C***BEGIN PROLOGUE  DERF
C***PURPOSE  Compute the error function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C8A, L5A1E
C***TYPE      REAL*8 (ERF-S, DERF-D)
C***KEYWORDS  ERF, ERROR FUNCTION, FNLIB, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C DERF(X) calculates the double precision error function for double
C precision argument X.
C
C Series for ERF        on the interval  0.          to  1.00000E+00
C                                        with weighted error   1.28E-32
C                                         log weighted error  31.89
C                               significant figures required  31.05
C                                    decimal places required  32.55
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, DCSEVL, DERFC, INITDS
C***REVISION HISTORY  (YYMMDD)
C   770701  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900727  Added EXTERNAL statement.  (WRB)
C   920618  Removed space from variable name.  (RWC, WRB)
C***END PROLOGUE  DERF
C============
C idecl:  explicitize implicit INTEGER declarations:
      IMPLICIT NONE
      INTEGER nterf,initds
C============
      REAL*8 X, ERFCS(21), SQEPS, SQRTPI, XBIG, Y, D1MACH,
     1  DCSEVL, DERFC
      LOGICAL FIRST
      EXTERNAL DERFC
      SAVE ERFCS, SQRTPI, NTERF, XBIG, SQEPS, FIRST
      DATA ERFCS(  1) / -CONST(.49046121234691808039984544033376,-1)/
      DATA ERFCS(  2) / -CONST(.14226120510371364237824741899631,+0)/
      DATA ERFCS(  3) / +CONST(.10035582187599795575754676712933,-1)/
      DATA ERFCS(  4) / -CONST(.57687646997674847650827025509167,-3)/
      DATA ERFCS(  5) / +CONST(.27419931252196061034422160791471,-4)/
      DATA ERFCS(  6) / -CONST(.11043175507344507604135381295905,-5)/
      DATA ERFCS(  7) / +CONST(.38488755420345036949961311498174,-7)/
      DATA ERFCS(  8) / -CONST(.11808582533875466969631751801581,-8)/
      DATA ERFCS(  9) / +CONST(.32334215826050909646402930953354,-10)/
      DATA ERFCS( 10) / -CONST(.79910159470045487581607374708595,-12)/
      DATA ERFCS( 11) / +CONST(.17990725113961455611967245486634,-13)/
      DATA ERFCS( 12) / -CONST(.37186354878186926382316828209493,-15)/
      DATA ERFCS( 13) / +CONST(.71035990037142529711689908394666,-17)/
      DATA ERFCS( 14) / -CONST(.12612455119155225832495424853333,-18)/
      DATA ERFCS( 15) / +CONST(.20916406941769294369170500266666,-20)/
      DATA ERFCS( 16) / -CONST(.32539731029314072982364160000000,-22)/
      DATA ERFCS( 17) / +CONST(.47668672097976748332373333333333,-24)/
      DATA ERFCS( 18) / -CONST(.65980120782851343155199999999999,-26)/
      DATA ERFCS( 19) / +CONST(.86550114699637626197333333333333,-28)/
      DATA ERFCS( 20) / -CONST(.10788925177498064213333333333333,-29)/
      DATA ERFCS( 21) / +CONST(.12811883993017002666666666666666,-31)/
      DATA SQRTPI / CONST(1.77245385090551602729816748334115,0)/
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DERF
      IF (FIRST) THEN
         NTERF = INITDS (ERFCS, 21, CONST(0.1,0)*AREAL(D1MACH(3)))
         XBIG = SQRT(-LOG(SQRTPI*D1MACH(3)))
         SQEPS = SQRT(CONST(2.0,0)*D1MACH(3))
      ENDIF
      FIRST = .FALSE.
C
      Y = ABS(X)
      IF (Y.GT.CONST(1.,0)) GO TO 20
C
C ERF(X) = 1.0 - ERFC(X)  FOR  -1.0 .LE. X .LE. 1.0
C
      IF (Y.LE.SQEPS) DERF = CONST(2.0,0)*X*X/SQRTPI
      IF (Y.GT.SQEPS) DERF = X*(CONST(1.0,0)+
     1   DCSEVL (CONST(2.,0)*X*X-CONST(1.,0),
     1   ERFCS, NTERF))
      RETURN
C
C ERF(X) = 1.0 - ERFC(X) FOR ABS(X) .GT. 1.0
C
 20   IF (Y.LE.XBIG) DERF = SIGN (CONST(1.0,0)-DERFC(Y), X)
      IF (Y.GT.XBIG) DERF = SIGN (CONST(1.0,0), X)
C
      RETURN
      END
! 22Jun2000 fgtok -s r8_precision.sub "r8con.csh conversion"
! 22Jun2000 fgtok
! 22Jun2000 fgtok
! 22Jun2000 fgtok
