Solution to the Test Exam

!_______________________________________________________________________
! Donev Aleksandar, PHY201, a Fortran 90 program
! that implements the trapezoid rule for numerical integration
! and integrates the function f(x)=4/(1+x^2) on the interval [0,1].
! The analytical result is known to be Pi.
!_______________________________________________________________________

!_______________________________________________________________________
! This module contains single and double precision kind constants for reals
!_______________________________________________________________________
module Precision
  implicit none
  public
  integer, parameter :: sp=kind(0.0e0), dp=kind(0.0d0) ! Single and double
  integer, parameter :: wp=sp ! Or dp
end module Precision  

!_______________________________________________________________________
! 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
  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)=4/(1+x^2)
!_______________________________________________________________________
module Function
  use Precision
  implicit none
  public

  contains
  !___________________________________
  ! This is the function f(x)
  !___________________________________
  function f(x)
    real(kind=wp), intent(in) :: x
    real(kind=wp) :: f
    f=4.0_wp/(1.0_wp+x**2)    
  end function f

end module Function
!_______________________________________________________________________


!_______________________________________________________________________
! The main program calls Trapzd passing f to it to obtain the integral:
! Integrate[1/(1+x^2),x=0..1]=Pi and writes the result to the file Pi.dat
!_______________________________________________________________________
program Pi
  use Precision
  use Trapzd
  use Function 
  implicit none
  
  integer :: N ! The number of intervals

  ! Open the file:
  open(unit=10,file="Pi.dat",action="write",status="unknown")
  
  ! Do the calculation
  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 integral value is: ",Trapezoid(f,0.0_wp,1.0_wp,N)  
  write(unit=10,fmt=*) "True value of Pi is: ",4.0_wp*atan(1.0_wp)
  
  ! Close the file:
  close(unit=10)

end program Pi