# Paste: pi calc polished

Author: Rex fortran Sat, 2 May 2009 22:17:33
Plain Text |
```! Pi calculator, by Reginald Ford
! PLEDGE: The JMU Honor Code was obeyed during this work
!----------------------------------------------------------------------
! The following is a module for simpson's rule along with
! a function, f(x), which is the semicircle function.
! Simpson's rule is used to integrate over the domain of the function,
! and this provides pi*r^2/2 . After multiplying by (2/r^2), the approximation of
! pi is printed. Several approximations are done to see the convergence.
!----------------------------------------------------------------------
module simpsons_rule

! r2 should be a perfect square so that there is
! no error in computing r, which is used in the bounds of integration
! This value seems to have no effect on the accuracy of the calculations
real (kind=8)::r2=2**30

contains

! The semicircle function for radius sqrt(r2)
real (kind = 8) function f(x)
implicit none
real (kind=8),intent(in)::x

f=sqrt(r2-x**2.)
! sqrt(x) provides nan for x<0
! So this is just in case
if(isnan(f)) then
f=0.0
end if
return
end function f

! This is the simple simpson's rule, which is basically
! quadratic interpolation and integration from a to b
! using (a+b)/2 as the midpoint
real (kind=8) function simple(a,b)
implicit none
real (kind=8),intent(in)::a,b
simple=((b-a)/6.)*(f(a)+4*f((a+b)/2.)+f(b))
return
end function simple

! This will integrate f(x) from a to b using simpson's rule
! from a to b with 'accuracy' controlling the number of nodes used
real (kind=8) function composite(a,b,accuracy)
implicit none
real (kind=8) , intent(in)::a,b
integer, intent(in)::accuracy
integer::n,i
real (kind=8)::h,integral

integral=0.0
n=accuracy
! the number of nodes needs to be odd
if(mod(n,2)==0) then
n=n+1
end if
h=(b-a)/n
do i=0,n-1
integral=integral+simple(a+i*h,a+(i+1)*h)
end do
composite=integral
return
end function composite

end module simpsons_rule

! MAIN PROGRAM STARTS HERE !
program pi_calc
use simpsons_rule
implicit none

integer::max_t=18
integer::t,n
real (kind=8)::approx,r

r=sqrt(r2)

write(*,*)"Approximations of pi:"
do t=1,max_t,2
! This allows us to see how it converges, with a large range of n covered
n=3**t
approx=composite((-1.)*r,r,n)*(2./r2)
write(*,*)"n=",n,":",approx," error:",acos(-1.)-approx
end do
! The program takes too long to converge to epsilon(double) precision on my machine
write(*,*)"epsilon=",epsilon(r)
end program pi_calc```