C @(#)erf.for	17.1.1.1 (ES0-DMD) 01/25/02 17:10:42
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
      DOUBLE PRECISION FUNCTION ERF(X)
C****************************************************************
C Error integral function
C
C ERF(x) = 1/sqrt(2*pi) * INTEGRAL [ exp(-0.5*u**2) du ]
C     from (-infinity) to x
C
C Reference: Ibbetson,D., 1963, Communications A.C.M. 6 (10), 616
C     1964, ibid. 7 (3), 148
C****************************************************************
C
C Authors: M. Rosa and O.-G. Richter, ESO Garching
C
      DOUBLE PRECISION Z,X,W,Y

      Z      = 0.0D0
      IF (X.EQ.0.0D0) GO TO 20
      Y      = 0.5D0*DABS(X)
      Z      = 1.0D0
      IF (Y.GE.3.0D0) GO TO 20
      IF (Y.GE.1.0D0) GO TO 10
      W      = Y*Y
      Z      = ((((((((0.000124818987*W-0.001075204047)*W+
     +         0.005198775019)*W-0.019198292004)*W+0.059054035642)*W-
     +         0.151968751364)*W+0.319152932694)*W-0.531923007300)*W+
     +         0.797884560593)*Y*2.0
      GO TO 20

   10 Y      = Y - 2.0D0
      Z      = (((((((((((((-0.000045255659*Y+0.000152529290)*Y-
     +         0.000019538132)*Y-0.000676904986)*Y+0.001390604284)*Y-
     +         0.000794620820)*Y-0.002034254874)*Y+0.006549791214)*Y-
     +         0.010557625006)*Y+0.011630447319)*Y-0.009279453341)*Y+
     +         0.005353579108)*Y-0.002141268741)*Y+0.000535310849)*Y +
     +         0.999936657524D0
   20 IF (X.GT.0.0D0) THEN
          ERF    = (Z+1.0D0)*0.5D0

      ELSE
          ERF    = (1.0D0-Z)*0.5D0
      END IF

      RETURN

      END