erfcc.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Function :  e r f c c
!
!> @file
!!
!! Computation of the complementary error function erfc(x) = 1-erf(x)
!! with a fractional error everywhere less than 1.2 x 10^(-7)
!! (formula by Press et al., 'Numerical Recipes in Fortran 77').
!!
!! @section Copyright
!!
!! Copyright 2009-2013 Reinhard Calov, Ralf Greve
!!
!! @section License
!!
!! This file is part of SICOPOLIS.
!!
!! SICOPOLIS is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! SICOPOLIS is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with SICOPOLIS.  If not, see <http://www.gnu.org/licenses/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Computation of the complementary error function erfc(x) = 1-erf(x)
!! with a fractional error everywhere less than 1.2 x 10^(-7)
!! (formula by Press et al., 'Numerical Recipes in Fortran 77').
!<------------------------------------------------------------------------------
function erfcc(x)

use sico_types

implicit none
real(dp)             :: erfcc
real(dp), intent(in) :: x

real(dp) :: t, z

z = abs(x)
t = 1.0_dp/(1.0_dp+0.5_dp*z)

erfcc = t * exp( -z*z     -1.26551223_dp &
                 + t  * (  1.00002368_dp &
                 + t  * (  0.37409196_dp &
                 + t  * (  0.09678418_dp &
                 + t  * ( -0.18628806_dp &
                 + t  * (  0.27886807_dp &
                 + t  * ( -1.13520398_dp &
                 + t  * (  1.48851587_dp &
                 + t  * ( -0.82215223_dp &
                 + t  *    0.17087277_dp ) ) ) ) ) ) ) ) )

if (x < 0.0_dp) erfcc = 2.0_dp-erfcc

end function erfcc
!