!
!_______________________________________________________________________
! Donev Aleksandar, PHY201, a Fortran 90 program
! that implements the trapezoid rule for numerical integration
! and integrates the function f(x)=exp(-x^2) on the interval [0,x].
! The analytical result is known to be sqrt(pi)/2*erf(x).
!_______________________________________________________________________
!
!_______________________________________________________________________
! This module implements the trapezoid rule for numerical integration
! of a smooth function f(x) over a finite interval [a,b] with
! a division with N intervals:
!
! Integrate[f(x),x=a..b] <> h*Sum[f(a+i*h),i=1..(N-1)]+(f(a)+f(b))/2
! where h=|b-a|/N
! Modeled after the solution of Khoa Lam for worksheet 8
!_______________________________________________________________________
MODULE Trapzd
      USE Precision, ONLY: wp
      PUBLIC
!
CONTAINS
   !___________________________________
   ! This function integrates the
   ! function f(x) over [a,b]
   ! using N intervals
   ! in the trapezoid rule:
   !___________________________________
      FUNCTION Trapezoid (f, a, b, N) RESULT (integral)
      !_______________
      ! Input arguments and result declarations:
      !_______________
         INTERFACE ! Interface for function to be integrated
            FUNCTION f (x) RESULT (f_x)
               USE Precision
               REAL (KIND=wp), INTENT (IN) :: x
               REAL (KIND=wp) :: f_x
            END FUNCTION f
         END INTERFACE
         REAL (KIND=wp), INTENT (IN) :: a, b ! The bounds
         INTEGER, INTENT (IN) :: N
         REAL (KIND=wp) :: integral ! The result
      !_______________
      ! Local variables declarations:
      !_______________
         REAL (KIND=wp) :: h, sum ! Step size and temporary for the sum
         INTEGER :: i ! Just a counter
      !_______________
      ! Calculation:
      !_______________
         h = (b-a) / N
         ! Calculate the sum:
         sum = 0.0_wp ! First initialize
         DO i = 1, N - 1
            sum = sum + f (a+i*h)! Or use x=x+h, starting at x=a
         END DO
         ! The result is:
         integral = h * sum + (h/2) * (f(a)+f(b))
      END FUNCTION Trapezoid
!
END MODULE Trapzd
!
!_______________________________________________________________________
! This module contains the function f(x)=exp(-x**2)
!_______________________________________________________________________
MODULE Function
      USE Precision, ONLY: wp
      IMPLICIT NONE
      PUBLIC
!
CONTAINS
  !___________________________________
  ! This is the function f(x)
  !___________________________________
      FUNCTION f (x)
         REAL (KIND=wp), INTENT (IN) :: x
         REAL (KIND=wp) :: f
         f = Exp (-x**2)
      END FUNCTION f
!
END MODULE Function
!_______________________________________________________________________
!
!
!_______________________________________________________________________
! The main program calls Trapzd passing f to it to obtain the integral:
! Integrate[exp(-t**2),t=0..x]=sqrt(pi)/2*erf(x)
! and writes the result to a file
!_______________________________________________________________________
PROGRAM Erf
      USE Precision, ONLY: wp
      USE Trapzd, ONLY: Trapezoid
      USE Function, ONLY: f
      USE Erf_Precise, ONLY: ErfPrecise
      IMPLICIT NONE
!
      INTEGER :: N ! The number of intervals
      REAL (KIND=wp) :: x ! We calculate erf(x)
!
      ! Open the file:
      OPEN (UNIT=10, FILE="Pi.dat", ACTION="write", STATUS="unknown")
!
      ! Do the calculation
      WRITE (UNIT=*, FMT=*) "The value of x?" ! Ask the user
      READ (UNIT=*, FMT=*) x ! The user enters N
      WRITE (UNIT=*, FMT=*) "How many intervals to use?" ! Ask the user
      READ (UNIT=*, FMT=*) N ! The user enters N
      ! Call Trapzd and write the result to the file:
      WRITE (UNIT=10, FMT=*) "The calculated value of erf(x) is: ", & 
             2.0_wp * Trapezoid (f, 0.0_wp, x, N) / Sqrt (4.0_wp*Atan(1.0_wp))
      WRITE (UNIT=10, FMT=*) "True value of erf(x) is: ", ErfPrecise (x)
!
      ! Close the file:
      CLOSE (UNIT=10)
!
END PROGRAM Erf
!