FUNCTION erf(x) RESULT(fn_val) !----------------------------------------------------------------------- ! EVALUATION OF THE REAL ERROR FUNCTION ! Based upon a Fortran 66 routine in the Naval Surface Warfare Center's ! Mathematics Library (1993 version). ! Adapted by Alan.Miller @ vic.cmis.csiro.au !----------------------------------------------------------------------- IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60) ! `Double precision' REAL (dp), INTENT(IN) :: x REAL (dp) :: fn_val ! Local variables REAL (dp), PARAMETER :: c = .564189583547756_dp, one = 1.0_dp, half = 0.5_dp, & zero = 0.0_dp REAL (dp), PARAMETER :: & a(5) = (/ .771058495001320D-04, -.133733772997339D-02, & .323076579225834D-01, .479137145607681D-01, & .128379167095513D+00 /), & b(3) = (/ .301048631703895D-02, .538971687740286D-01, & .375795757275549D+00 /), & p(8) = (/ -1.36864857382717D-07, 5.64195517478974D-01, & 7.21175825088309D+00, 4.31622272220567D+01, & 1.52989285046940D+02, 3.39320816734344D+02, & 4.51918953711873D+02, 3.00459261020162D+02 /), & q(8) = (/ 1.00000000000000D+00, 1.27827273196294D+01, & 7.70001529352295D+01, 2.77585444743988D+02, & 6.38980264465631D+02, 9.31354094850610D+02, & 7.90950925327898D+02, 3.00459260956983D+02 /), & r(5) = (/ 2.10144126479064D+00, 2.62370141675169D+01, & 2.13688200555087D+01, 4.65807828718470D+00, & 2.82094791773523D-01 /), & s(4) = (/ 9.41537750555460D+01, 1.87114811799590D+02, & 9.90191814623914D+01, 1.80124575948747D+01 /) REAL (dp) :: ax, bot, t, top, x2 !------------------------- ax = ABS(x) IF (ax <= half) THEN t = x*x top = ((((a(1)*t + a(2))*t + a(3))*t + a(4))*t + a(5)) + one bot = ((b(1)*t + b(2))*t + b(3))*t + one fn_val = x*(top/bot) RETURN END IF IF (ax <= 4.0_dp) THEN top = ((((((p(1)*ax + p(2))*ax + p(3))*ax + p(4))*ax + p(5))*ax & + p(6))*ax + p(7))*ax + p(8) bot = ((((((q(1)*ax + q(2))*ax + q(3))*ax + q(4))*ax + q(5))*ax & + q(6))*ax + q(7))*ax + q(8) fn_val = half + (half - EXP(-x*x)*top/bot) IF (x < zero) fn_val = -fn_val RETURN END IF IF (ax < 5.8_dp) THEN x2 = x*x t = one / x2 top = (((r(1)*t + r(2))*t + r(3))*t + r(4))*t + r(5) bot = (((s(1)*t + s(2))*t + s(3))*t + s(4))*t + one fn_val = (c - top/(x2*bot)) / ax fn_val = half + (half - EXP(-x2)*fn_val) IF (x < zero) fn_val = -fn_val RETURN END IF fn_val = SIGN(one, x) RETURN END FUNCTION erf