!****************************************************
!* Program to demonstrate ASYMERF *
!* ------------------------------------------------ *
!* Reference: BASIC Scientific Subroutines, Vol. II *
!* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].*
!* *
!* F90 Version by J.-P. Moreau, Paris. *
!* (www.jpmoreau.fr) *
!* ------------------------------------------------ *
!* SAMPLE RUN: *
!* *
!* Find the value of ERF(X)=2*Exp(-X*X)/SQRT(PI) *
!* *
!* Input X ? 3 *
!* ERF(X)= .9999779 with error estimate= -.00000000 *
!* Number of terms evaluated was 10 *
!* *
!* Input X ? 4 *
!* ERF(X)= 1.0000000 with error estimate= 0.0000000 *
!* Number of terms evaluated was 17 *
!* *
!****************************************************
PROGRAM ASYMERF
real*8 e, x, y
integer n
print *, ' '
write(*,'(" Input X: ")',advance='no')
read *, x
call Asym_erf(e,n,x,y)
write(*,50) y, e
write(*,55) n
stop
50 format(/' ERF(X)= ',F10.7,' with error estimate= ',F10.7/)
55 format(' Number of terms evaluated was ',i2//)
end
!***********************************************************
!* Asymptotic series expansion of the integral of *
!* 2 exp(-X*X)/(X*sqrt(PI)), the normalized error function *
!* (ASYMERF). This program determines the values of the *
!* above integrand using an asymptotic series which is *
!* evaluated to the level of maximum accuracy. *
!* The integral is from 0 to X. The input parameter, X *
!* must be > 0. The results are returned in Y and Y1, *
!* with the error measure in E. The number of terms used *
!* is returned in N. The error is roughly equal to first *
!* term neglected in the series summation. *
!* ------------------------------------------------------- *
!* Reference: A short table of integrals by B.O. Peirce, *
!* Ginn and Company, 1957. *
!***********************************************************
Subroutine Asym_erf(e,n,x,y)
! Labels used: 100, 200
real*8 e, x, y
real*8 c1, c2, pi, y1
integer n
pi = 4.d0*datan(1.d0)
n = 1; y = 1.d0
c2 = 1.d0 / (2.d0 * x * x)
100 y = y - c2
n = n + 2; c1 = c2
c2 = -c1 * n / (2.d0 * x * x)
! Test for divergence - The break point is roughly N=X*X
if (dabs(c2) > dabs(c1)) goto 200
! Continue summation
goto 100
200 n = (n + 1) / 2
e = dexp(-x * x) / (x * dsqrt(pi))
y1 = y * e
y = 1.d0 - y1
e = e * c2
return
end
! End of file asymerf.f90