Error function - computed from the rational approximations of W. J. Cody, Math. Comp. 22 (1969), pages 631-637.

for abs(x) le 0.47 erf is calculated directly for abs(x) gt 0.47 erf is calculated via erf(x)=1-erfc(x)

Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|

real(kind=dp), | intent(in) | :: | x |

```
function qe_erf(x)
!---------------------------------------------------------------------
!
!! Error function - computed from the rational approximations of
!! W. J. Cody, Math. Comp. 22 (1969), pages 631-637.
!!
!! for abs(x) le 0.47 erf is calculated directly
!! for abs(x) gt 0.47 erf is calculated via erf(x)=1-erfc(x)
!
use w90_constants, only: dp
implicit none
real(kind=dp), intent(in) :: x
real(kind=dp) :: x2, p1(4), q1(4)
real(kind=dp) :: qe_erf
data p1/2.426679552305318E2_dp, 2.197926161829415E1_dp, &
6.996383488619136_dp, -3.560984370181538E-2_dp/
data q1/2.150588758698612E2_dp, 9.116490540451490E1_dp, &
1.508279763040779E1_dp, 1.000000000000000_dp/
!
if (abs(x) > 6.0_dp) then
!
! erf(6)=1-10^(-17) cannot be distinguished from 1
!
qe_erf = sign(1.0_dp, x)
else
if (abs(x) <= 0.47_dp) then
x2 = x**2
qe_erf = x*(p1(1) + x2*(p1(2) + x2*(p1(3) + x2*p1(4)))) &
/(q1(1) + x2*(q1(2) + x2*(q1(3) + x2*q1(4))))
else
qe_erf = 1.0_dp - qe_erfc(x)
endif
endif
!
return
end function qe_erf
```