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 :: converged
   REAL(KIND=wp) :: error ! Error allowed in erf(x)

   ! This is a generic function interface
   ! One can call ErfSeries with either REAL or COMPLEX
   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/factor) ! 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/factor) ! 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=-CMPLX(2*k-1,0,wp)/CMPLX(2*k+1,0,wp)/CMPLX(k,0,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        

PROGRAM Erf_Module
   USE Precision
   USE Erf_Series, ONLY: error, converged, ErfSeries

   ! We are in position to choose whether we want REAL or COMPLEX x:
   ! Comment one of these out:
        REAL(KIND=wp) :: x, erf ! x and erf(x)  
        ! COMPLEX (KIND=wp) :: x, erf ! x and erf(x)
   INTEGER :: status ! I/O status
   
   MainLoop: DO
   
      WRITE(UNIT=*,FMT="(A)",ADVANCE="NO") "Enter the value of x and epsilon ('S' to stop): "
      READ(UNIT=*,FMT=*,IOSTAT=status) x,error      
      IF(status/=0) EXIT MainLoop

      erf=ErfSeries(x)
      IF(converged) THEN
        WRITE(UNIT=*,FMT=*) "Erf (",x,") = ",erf
      ELSE      
        WRITE(UNIT=*,FMT=*) "Convergence was not achieved! Last known answer:"
        WRITE(UNIT=*,FMT=*) "Erf (",x,") = ",erf  
      END IF     
      
   END DO MainLoop               

END PROGRAM Erf_Module