!**************************************************** !* 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