!___________________________________________________________________________________________________
! 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
!