ErrorFunction.F

c***********************************************************************
#include "author.inc"
c*    $Id: ErrorFunction.F,v 1.3 1995/07/30 17:57:21 turner Exp $
c*
c*    Use rational approximation to compute error function.
c*
c*    erf x =
c*       1 - (a1*t + a2*t^2 + a3*t^3 + a4*t^4 + a5*t^5) e^(-x^2) + eps
c*
c*    where:
c*                  1
c*           t = ---------  ,  ABS(eps) <= 1.5e-7 ,
c*                1 + p*x
c*
c*    and p, a1, a2, a3, a4, and a5 are constants.
c*
c*    Reference:
c*      Handbook of Mathematical Functions
c*      M. Abramowitz and I. E. Stegun, ed.
c*      Tenth Printing, 12/72, p. 299
c*
c*    <PARAMETER LIST>
c*
c*     Input:
c*      x - argument
c*
c*     Output:
c*      JT_ErrorFunction - error function evaluated at x
c*
#include "copyright.inc"
c****************************************************************
      real function JT_ErrorFunction(x)
      implicit none
c
      real x, p, a1, a2, a3, a4, a5, t, t2, t3, t4, one
      parameter (one = 1.0d0,
     &           p  =  0.3275911d0,
     &           a1 =  0.254829592d0,
     &           a2 = -0.284496736d0,
     &           a3 =  1.421413741d0,
     &           a4 = -1.453152027d0,
     &           a5 =  1.061405429d0)
c
      t = one / (one + p*x)
      t2 = t*t
      t3 = t*t2
      t4 = t*t3
c
      JT_ErrorFunction = one -
     &   ( a1*t + a2*t2 + a3*t3  + a4*t4 + a5*t4*t )*EXP(-x*x)
c
      return
      end