erf.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE erf_mod
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  These routines compute the incomplete gamma function and error      !
!  function:                                                           !
!                                                                      !
!  ERF        Error function, ERF(x)                                   !
!  ERFC       Complementary error function, ERFC(x)                    !
!  ERFCC      Complementary error function, ERFCC(x): cheaper          !
!               Chebyshev fitting approximation.                       !
!  GAMMP      Incomplete gamma function, P(a,x)                        !
!  GAMMQ      Incomplete gamma function complement, Q(a,x)=1-P(a,x)    !
!                                                                      !
!  Adapted from Numerical Recipes:                                     !
!                                                                      !
!    Press, W.H, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery,   !
!     1986:  Numerical Recipes in Fortran 77,  The Art of Parallel     !
!     Scientific Computing, 1st Edition, Cambridge Univ. Press.        !
!                                                                      !
!=======================================================================
!
      USE mod_kinds
!
      implicit none
!
      PRIVATE
      PUBLIC :: erf
      PUBLIC :: erfc
      PUBLIC :: erfcc
      PUBLIC :: gammp
      PUBLIC :: gammq
!
      CONTAINS
!
      FUNCTION erf (x) RESULT (value)
!
!=======================================================================
!                                                                      !
!  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
      END FUNCTION erf
!
      FUNCTION erfc (x) RESULT (value)
!
!=======================================================================
!                                                                      !
!  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
      END FUNCTION erfc
!
      FUNCTION erfcc (x) RESULT (value)
!
!=======================================================================
!                                                                      !
!  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
      END FUNCTION erfcc
!
      FUNCTION gammln (xx) RESULT (value)
!
!=======================================================================
!                                                                      !
!  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
      END FUNCTION gammln
!
      FUNCTION gammp (a,x) RESULT (value)
!
!=======================================================================
!                                                                      !
!  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
      END FUNCTION gammp
!
      FUNCTION gammq (a,x) RESULT (value)
!
!=======================================================================
!                                                                      !
!  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
      END FUNCTION gammq
!
      SUBROUTINE gser (gamser, a, x, gln)
!
!=======================================================================
!                                                                      !
!  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
      END SUBROUTINE gser
!
      SUBROUTINE gcf (gammcf, a, x, gln)
!
!=======================================================================
!                                                                      !
!  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
      END SUBROUTINE gcf
!
      FUNCTION aerf (arg) RESULT (value)
!
!=======================================================================
!                                                                      !
!  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
      END FUNCTION aerf
 
      END MODULE erf_mod