erf_mod Module Reference

Functions/Subroutines
real(r8) function, public	erf (x)
real(r8) function, public	erfc (x)
real(r8) function, public	erfcc (x)
real(r8) function	gammln (xx)
real(r8) function, public	gammp (a, x)
real(r8) function, public	gammq (a, x)
subroutine	gser (gamser, a, x, gln)
subroutine	gcf (gammcf, a, x, gln)
real(r8) function	aerf (arg)

Function/Subroutine Documentation

aerf()

real(r8) function erf_mod::aerf ( real(r8), intent(in) arg )

private

Definition at line 409 of file erf.F.

!
!=======================================================================
!                                                                      !
!  The error function has been taken in approximate form from          !
!  Abramowitz and Stegun, page 299.                                    !
!                                                                      !
!  Reference:                                                          !
!                                                                      !
!  Abramowitz, M. and I.A. Stegun, 1972: Handbook of Mathematical      !
!       Functions with Formulas, Graphs, and Mathematical Tables,      !
!       National Bureau of Standards, Applied Mathematics Series       !
!       55, US Goverment Printing Office, Washington, D.C.             !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      real(r8), intent(in) :: arg
!
!  Local variable declarations.
!
      real(r8), parameter :: a1 = +0.254829592_r8
      real(r8), parameter :: a2 = -0.284496736_r8
      real(r8), parameter :: a3 = -1.421413741_r8
      real(r8), parameter :: a4 = -1.453152027_r8
      real(r8), parameter :: a5 = +1.061405429_r8
      real(r8), parameter :: p  = +0.3275911_r8
 
      real(r8) :: big, c, my_sign, t, y
      real(r8) :: value
!
!-----------------------------------------------------------------------
!  Compute approximated error function.
!-----------------------------------------------------------------------
!
      big=20.0_r8
      my_sign=1.0_r8
 
      IF (arg.lt.0.0_r8) THEN
        my_sign=-1.0
      ELSE
        my_sign=1.0_r8
      END IF
 
      IF (arg.eq.0.0_r8) THEN
        value=0.0_r8
        RETURN
      END IF
!
!  Polynomial approximation.
!
      y=abs(arg)
      t=1.0_r8/(1.0_r8+p*y)
      c=t*(a1+t*(a2+t*(a3+t*(a4+a5*t))))
      IF (y*y.lt.big) THEN
        value=my_sign*(1.0_r8-c*exp(-y*y))
      ELSE
        value=my_sign
      END IF
 
      RETURN

erf()

real(r8) function, public erf_mod::erf

(

real(r8), intent(in)

x

)

Definition at line 42 of file erf.F.

!
!=======================================================================
!                                                                      !
!  This routine computes the error function, ERF(x).                   !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      real(r8), intent(in) :: x
!
!  Local variable declarations.
!
      real(r8) :: value
!
!-----------------------------------------------------------------------
!  Compute error function in terms of the incomplete gamma function,
!  P(a,x).
!-----------------------------------------------------------------------
!
      IF (x.lt.0.0_r8) THEN
        value=-gammp(0.5_r8, x**2)
      ELSE
        value= gammp(0.5_r8, x**2)
      ENDIF
 
      RETURN

References gammp().

Referenced by analytical_mod::ana_nlminitial_tile().

Here is the call graph for this function:

Here is the caller graph for this function:

erfc()

real(r8) function, public erf_mod::erfc

(

real(r8), intent(in)

x

)

Definition at line 72 of file erf.F.

!
!=======================================================================
!                                                                      !
!  This routine computes the complementary error function, ERFC(x).    !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      real(r8), intent(in) :: x
!
!  Local variable declarations.
!
      real(r8) :: value
!
!-----------------------------------------------------------------------
!  Compute complementary error function in terms of the incomplete
!  gamma functions, P(a,x) and Q(a,x).
!-----------------------------------------------------------------------
!
      IF (x.lt.0.0_r8) THEN
        value=1.0_r8+gammp(0.5_r8, x**2)
      ELSE
        value=gammq(0.5_r8, x**2)
      endif
 
      RETURN

References gammp(), and gammq().

Here is the call graph for this function:

erfcc()

real(r8) function, public erf_mod::erfcc

(

real(r8), intent(in)

x

)

Definition at line 102 of file erf.F.

!
!=======================================================================
!                                                                      !
!  This routine computes the complementary error function, ERFCC(x),   !
!  with fractional error everywhere less than 1.2E-7.                  !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      real(r8), intent(in) :: x
!
!  Local variable declarations.
!
      real(r8) :: t, value, z
!
!-----------------------------------------------------------------------
!  Compute complementary error function based on Chebyshev fitting to
!  an inspired guess.
!-----------------------------------------------------------------------
!
      z=abs(x)
      t=1.0_r8/(1.0_r8+0.5_r8*z)
      value=t*exp(-z*z-1.26551223_r8+                                   &
     &            t*(1.00002368_r8+                                     &
     &               t*(0.37409196_r8+                                  &
     &                  t*(0.09678418_r8+                               &
     &                     t*(-0.18628806_r8+                           &
     &                        t*(0.27886807_r8+                         &
     &                           t*(-1.13520398_r8+                     &
     &                              t*(1.48851587_r8+                   &
     &                                 t*(-0.82215223_r8+               &
     &                                    t*.17087277)))))))))
      IF (x.lt.0.0_r8) value=2.0-value
 
      RETURN

gammln()

real(r8) function erf_mod::gammln ( real(r8), intent(in) xx )

private

Definition at line 141 of file erf.F.

!
!=======================================================================
!                                                                      !
!  This routine computes the natural logarithm of the gamma function,  !
!  LN(Gamma(xx)). The natural logarithm is computed to avoid floating- !
!  point overflow.                                                     !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      real(r8), intent(in) :: xx
!
!  Local variable declarations.
!
      integer :: j
 
      real(r8) :: ser, tmp, x, y
      real(r8) :: value
      real(r8), save, dimension(6) :: cof =                             &
     &          (/ 76.18009172947146_r8,                                &
     &            -86.50532032941677_r8,                                &
     &             24.01409824083091_r8,                                &
     &             -1.231739572450155_r8,                               &
     &              0.1208650973866179E-2_r8,                           &
     &             -0.5395239384953E-5_r8 /)
      real(r8), save :: stp = 2.5066282746310005_r8
!
!-----------------------------------------------------------------------
!  Compute natural logarithm of the gamma function.
!-----------------------------------------------------------------------
!
      x=xx
      y=x
      tmp=x+5.5_r8
      tmp=(x+0.5_r8)*log(tmp)-tmp
      ser=1.000000000190015_r8
      DO j=1,6
        y=y+1.d0
        ser=ser+cof(j)/y
      END DO
      value=tmp+log(stp*ser/x)
 
      RETURN

References mod_kinds::r8.

Referenced by gcf(), and gser().

Here is the caller graph for this function:

gammp()

real(r8) function, public erf_mod::gammp	(	real(r8), intent(in)	a,
		real(r8), intent(in)	x )

Definition at line 188 of file erf.F.

!
!=======================================================================
!                                                                      !
!  This routine computes the incomplete gamma function, P(a,x).        !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_scalars
      USE mod_iounits
!
!  Imported variable declarations.
!
      real(r8), intent(in) :: a, x
!
!  Local variable declarations.
!
      real(r8) :: gammcf, gamser, gln
      real(r8) :: value
!
!-----------------------------------------------------------------------
!  Compute incompleate gamma fucntion, P(a.x).
!-----------------------------------------------------------------------
!
      IF ((x.lt.0.0_r8).or.(a.le.0.0_r8)) THEN
        WRITE (stdout,10) a, x
  10    FORMAT (/,' GAMMAP - gamma function negative argument,',        &
     &          ' a = ',1pe13.6,'  x = ',1pe13.6)
        exit_flag=8
      END IF
 
      IF (x.lt.a+1.0_r8) THEN
        CALL gser (gamser, a, x, gln)
        value=gamser                 ! series representation
      ELSE
        CALL gcf (gammcf, a, x, gln)
        value=1.0_r8-gammcf          ! continued fraction representation
      ENDIF
 
      RETURN

References mod_scalars::exit_flag, gcf(), gser(), and mod_iounits::stdout.

Referenced by erf(), and erfc().

Here is the call graph for this function:

Here is the caller graph for this function:

gammq()

real(r8) function, public erf_mod::gammq	(	real(r8), intent(in)	a,
		real(r8), intent(in)	x )

Definition at line 231 of file erf.F.

!
!=======================================================================
!                                                                      !
!  This routine computes the incomplete gamma function complement,     !
!  Q(a,x)=1-P(a,x).                                                    !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_scalars
      USE mod_iounits
!
!  Imported variable declarations.
!
      real(r8), intent(in) :: a, x
!
!  Local variable declarations.
!
      real(r8) :: gammcf, gamser, gln
      real(r8) :: value
!
!-----------------------------------------------------------------------
!  Compute incompleate gamma fucntion, Q(a,x)=1-P(a.x).
!-----------------------------------------------------------------------
!
      IF ((x.lt.0.0_r8).or.(a.le.0.0_r8)) THEN
        WRITE (stdout,10) a, x
  10    FORMAT (/,' GAMMAQ - gamma function negative argument,',        &
     &          ' a = ',1pe13.6,'  x = ',1pe13.6)
        exit_flag=8
      END IF
 
      IF (x.lt.a+1.0_r8) THEN
        CALL gser (gamser, a, x, gln)
        value=1.0_r8-gamser          ! series representation
      ELSE
        CALL gcf (gammcf, a, x, gln)
        value=gammcf                 ! continued fraction representation
      ENDIF
 
      RETURN

References mod_scalars::exit_flag, gcf(), gser(), and mod_iounits::stdout.

Referenced by erfc().

Here is the call graph for this function:

Here is the caller graph for this function:

gcf()

subroutine erf_mod::gcf ( real(r8), intent(out) gammcf, real(r8), intent(in) a, real(r8), intent(in) x, real(r8), intent(out) gln )

private

Definition at line 342 of file erf.F.

!
!=======================================================================
!                                                                      !
!  This routine computes the incomplete gamma function P(a,x)          !
!  evaliuated by its series representation as gamser.                  !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_scalars
      USE mod_iounits
!
!  Imported variable declarations.
!
      real(r8), intent(in) :: a, x
      real(r8), intent(out) :: gammcf, gln
!
!  Local variable declarations.
!
      logical :: Converged
 
      integer, parameter :: ITMAX = 100          ! number of iterations
      integer :: i
 
      real(r8), parameter :: EPS = 3.0e-7_r8     ! relative accuracy
      real(r8), parameter :: FPMIN = 1.0e-30_r8  ! smallest number
      real(r8) :: an, b, c, d, del, h
!
!-----------------------------------------------------------------------
!  Compue incomplete gamma function evaluated by its continued fraction
!  representation.
!-----------------------------------------------------------------------
!
      gln=gammln(a)
      b=x+1.0_r8-a
      c=1.0_r8/fpmin
      d=1.0_r8/b
      h=d
      converged=.false.
      DO i=1,itmax
        an=-real(i,r8)*(real(i,r8)-a)
        b=b+2.0_r8
        d=an*d+b
        IF (abs(d).lt.fpmin) d=fpmin
        c=b+an/c
        IF (abs(c).lt.fpmin) c=fpmin
        d=1.0_r8/d
        del=d*c
        h=h*del
        IF (abs(del-1.0_r8).lt.eps) THEN
          converged=.true.
          EXIT
        END IF
      END DO
      IF (converged) THEN
        gammcf=exp(-x+a*log(x)-gln)*h
      ELSE
        WRITE (stdout,10) itmax
  10    FORMAT (/,' GCF - Gamma function not converged, ITMAX = ',      &
     &          i4.4,/,8x,'a is too large')
        exit_flag=8
      END IF
 
      RETURN

References mod_scalars::exit_flag, gammln(), and mod_iounits::stdout.

Referenced by gammp(), and gammq().

Here is the call graph for this function:

Here is the caller graph for this function:

gser()

subroutine erf_mod::gser ( real(r8), intent(out) gamser, real(r8), intent(in) a, real(r8), intent(in) x, real(r8), intent(out) gln )

private

Definition at line 275 of file erf.F.

!
!=======================================================================
!                                                                      !
!  This routine computes the incomplete gamma function P(a,x)          !
!  evaliuated by its series representation as gamser.                  !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_scalars
      USE mod_iounits
!
!  Imported variable declarations.
!
      real(r8), intent(in) :: a, x
      real(r8), intent(out) :: gamser, gln
!
!  Local variable declarations.
!
      logical :: Converged
      integer, parameter :: ITMAX = 100
      integer :: i
 
      real(r8), parameter :: EPS =3.0e-7_r8
      real(r8) :: ap, del, my_sum
!
!-----------------------------------------------------------------------
!  Compue incomplete gamma function by its series representation.
!-----------------------------------------------------------------------
!
      gln=gammln(a)
      IF (x.le.0.0_r8) THEN
        IF (x.lt.0.0) THEN
          WRITE (stdout,10) x
  10      FORMAT (/,' GSER - gamma function negative argument, x = ',   &
     &            1pe13.6)
          exit_flag=8
        END IF
        gamser=0.0_r8
        RETURN
      END IF
      ap=a
      my_sum=1.0_r8/a
      del=my_sum
      converged=.false.
      DO i=1,itmax
        ap=ap+1.0_r8
        del=del*x/ap
        my_sum=my_sum+del
        IF (abs(del).lt.abs(my_sum)*eps) THEN
          converged=.true.
          EXIT
        END IF
      END DO
      IF (converged) THEN
        gamser=my_sum*exp(-x+a*log(x)-gln)
      ELSE
        WRITE (stdout,20) itmax
  20    FORMAT (/,' GSER - Gamma function not converged, ITMAX = ',     &
     &          i4.4,/,8x,'a is too large')
        exit_flag=8
      END IF
 
      RETURN

References mod_scalars::exit_flag, gammln(), and mod_iounits::stdout.

Referenced by gammp(), and gammq().

Here is the call graph for this function:

Here is the caller graph for this function: