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