!___________________________________________________________________________________________________
! Aleksandar Donev
! PHY201 Solutions to Worksheet
! December, 2000, MSU
!___________________________________________________________________________________________________
MODULE Precision
PUBLIC
INTEGER, PARAMETER :: sp=KIND(0.0E0), dp=KIND(0.0D0)
INTEGER, PARAMETER :: wp=sp ! Use single precision in this run
END MODULE Precision
!
MODULE Erf_Series
USE Precision
!
IMPLICIT NONE
PUBLIC
INTEGER, PARAMETER, PRIVATE :: max_iterations = 100 ! Max # of iterations
LOGICAL, SAVE :: converged = .TRUE.
REAL (KIND=wp) :: error ! Error allowed in erf(x)
!
INTERFACE ErfSeries
MODULE PROCEDURE ErfReal
MODULE PROCEDURE ErfComplex
END INTERFACE
!
CONTAINS ! Do not forget this statement
!
! This is the error function for real arguments
FUNCTION ErfReal (x) RESULT (erf)
REAL (KIND=wp), INTENT (IN) :: x ! This is the function argument
REAL (KIND=wp) :: erf ! Function result erf(x)
!
INTEGER :: k ! Iteration counter
REAL (KIND=wp) :: sum, factor, expression ! A few auxiliary variables
REAL (KIND=wp) :: epsilon ! Error allowed in the sum
REAL (KIND=wp) :: pi ! The constant pi
!
pi = 4 * Atan (1.0_wp)! This is a good way of calculating pi to the needed precision
factor = 2 * x / Sqrt (pi)! The factor in front of the sum
epsilon = Abs (error) / (Abs(factor)+TINY(0.0_wp))! This is the error allowed *in* the sum
!
Magnitude: IF (Abs(x) < 5.0_wp) THEN ! We can use the series provided
sum = 1.0_wp ! This is E(0)
expression = 1.0_wp ! This is E(0) as well
k = 0 ! Initialize counter
Summation: DO
k = k + 1 ! Iteration counter
! Calculate the next term to be added
expression = - REAL (2*k-1, wp) / REAL (2*k+1, wp) / REAL (k, wp) * (x**2) * expression
! Now check for convergence:
IF (Abs(expression) < epsilon) THEN ! Converged
! converged=.TRUE.
erf = factor * sum
EXIT Summation
ELSE IF (k >= max_iterations) THEN ! Not converged
converged = .FALSE.
erf = factor * sum
EXIT Summation
ELSE ! Keep on adding terms
sum = sum + expression ! Add the term to the sum
CYCLE Summation ! This can be ommitted
END IF
END DO Summation
ELSE ! If x is large in magnitude, we can not use the above series easily
converged = .FALSE.
erf = 1.0_wp
END IF Magnitude
!
END FUNCTION ErfReal
!
! This is the error function for complex arguments
! It is an almost exact copy of the previous function
FUNCTION ErfComplex (x) RESULT (erf)
COMPLEX (KIND=wp), INTENT (IN) :: x ! This is the function argument
COMPLEX (KIND=wp) :: erf ! Function result erf(x)
!
INTEGER :: k ! Iteration counter
COMPLEX (KIND=wp) :: sum, factor, expression ! A few auxiliary variables
REAL (KIND=wp) :: epsilon ! Error allowed in the sum
REAL (KIND=wp) :: pi ! The constant pi
!
pi = 4 * Atan (1.0_wp)! This is a good way of calculating pi to the needed precision
factor = 2 * x / Sqrt (pi)! The factor in front of the sum
epsilon = Abs (error) / (Abs(factor)+TINY(0.0_wp))! This is the error allowed *in* the sum
!
Magnitude: IF (Abs(x) < 5.0_wp) THEN ! We can use the series provided
sum = (1.0_wp, 0.0_wp)! This is E(0)
expression = (1.0_wp, 0.0_wp)! This is E(0) as well
k = 0 ! Initialize counter
Summation: DO
k = k + 1 ! Iteration counter
! Calculate the next term to be added
expression = - REAL (2*k-1, wp) / REAL (2*k+1, wp) / REAL (k, wp) * (x**2) * expression
! Now check for convergence:
IF (Abs(expression) < epsilon) THEN ! Converged
! converged=.TRUE.
erf = factor * sum
EXIT Summation
ELSE IF (k >= max_iterations) THEN ! Not converged
converged = .FALSE.
erf = factor * sum
EXIT Summation
ELSE ! Keep on adding terms
sum = sum + expression ! Add the term to the sum
CYCLE Summation ! This can be ommitted
END IF
END DO Summation
ELSE ! If x is large in magnitude, we can not use the above series easily
converged = .FALSE.
erf = (0.0_wp, 0.0_wp)! There is no clear answer in this case
END IF Magnitude
!
END FUNCTION ErfComplex
!
END MODULE Erf_Series
!