Paste: pi calc polished

! Pi calculator, by Reginald Ford
! James Madison University
! 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