1!      program main
2!      integer i,ndata
3!      double precision x,erf
4!      ndata=20
5!      x=0.0d0
6!      do i=1,ndata
7!       call merf(x,erf)
8!       write(*,1) i,x,erf
9!1      format('(x,erf)(',i3,')=(',f15.8,',',f15.8,')')
10!       x=x+0.10d0
11!      enddo
12!      end
13       subroutine merf (y, erf)
14
15!  ======================================================================
16!  purpose:  evaluation of the error function.
17!            two related entries are provided: the complementary
18!            error function and the normal distribution function.
19!            (entries merfc and mdnor)
20!
21!  input  :  y      - the argument of the function
22!  output :  erf    - the function value
23!                     merf : ranging from -1 to +1 for increasing y
24!                     merfc  : ranging from +2 to 0
25!                     mdnor  : max. for y=0, falling off to zero for
26!                              increasing abs(y)
27!
28!  *=====================================================================
29!  ======================================================================
30      double precision c,eight,half,one,r0p47,
31     &                 r5p5,small,sqr1d2,tenth,thirdm,
32     &                 two,zero,a,p,q,x,y,erf,sn,xx
33      integer kret,md,isw
34!
35      parameter (c = 1.1283791670955d0, eight = 8.0d0,
36     &           half = 0.5d0, one = 1.0d0, r0p47 = 0.47d0,
37     &           r5p5 = 5.5d0, small = 1.0d-15,
38     &           sqr1d2 = 0.70710678118655d0,
39     &           tenth = 0.1d0, thirdm =
40     &          -one/3.0d0, ttt = -one/1320.0d0, two = 2.0d0,
41     &           zero = 0.0d0)
42!
43      dimension a(10), p(8), q(9)
44!
45      data a /one, thirdm, tenth, -0.02380923809238d0,
46     &        0.4629629629630d-2,
47     &       ttt, 1.0683760683761d-4,
48     &       -1.3227513227513d-5, 1.4589169000934d-6,
49     &       -1.45038552223150d-7/
50!
51      data p /883.4789426085d0, 1549.6793124037d0,
52     &        1347.1941340976d0, 723.04000277753d0,
53     &        255.50049469496d0, 59.240010112914d0,
54     &        8.3765310814197d0, 0.56418955944261d0/
55!
56      data q /883.4789426085d0, 2546.5785458098d0,
57     &        3337.2213699893d0, 2606.7120152651d0,
58     &        1333.5699756800d0, 460.28512369160d0,
59     &        105.50025439769d0, 14.847012237523d0, one/
60
61             kret = 2
62          !
63             x  = y
64             md = 0
65          !
66             if (x<zero) then
67                isw = 1
68                x   = - x
69             else
70                isw = 2
71             end if
72          !
73             if (kret/=2) goto 120
74             if (x<=r5p5) goto 40
75          !
76             erf = one
77             if (isw==1) erf = - one
78             goto 150
79
80          !  -----------------------------------------------------------
81          !  abs(x) less than .47 compute erf by taylor series expansion
82          !  -----------------------------------------------------------
83
84   40        if (x<=r0p47) goto 90
85             kret = 1
86
87          !  -----------------------------------------------------
88          !  abs(x) between .47 and 10. compute complemented error
89          !  function by a rational function of x
90          !  -----------------------------------------------------
91
92   50        sn = p(8)
93            do i = 7, 1, -1
94               sn = sn*x + p(i)
95            end do
96          !
97             sd = q(9)
98             do i = 8, 1, -1
99                sd = sd*x + q(i)
100             end do
101          !
102             erf = (sn/sd)*exp(-x*x)
103             if (kret/=1) goto 80
104
105          !  -----------------------------------
106          !  compute complemented error function
107          !  -----------------------------------
108
109             erf = one - erf
110             if (isw/=2) erf = - erf
111             goto 150
112          !
113   80         if (isw/=2) erf = two - erf
114             if (md/=0) erf = half*erf
115             goto 150
116          !
117   90        xx = x*x
118          !
119             erf = a(10)
120             do i = 9, 1, -1
121                erf = erf*xx + a(i)
122             end do
123          !
124             erf = c*erf*x
125          !
126             if (kret==2) then
127                if (isw/=2) erf = - erf
128                goto 150
129             else
130                erf = one - erf
131                goto 80
132             end if
133          !
134  120        if (x>eight) then
135                erf = zero
136                goto 80
137          !
138             else if (x>r0p47) then
139                kret = 2
140                goto 50
141          !
142             else if (x>small) then
143                goto 90
144          !
145             else
146                erf = one
147                goto 80
148             end if
149          !
150  150     return
151          end
152c $Id$
153