program inegrate
implicit none
! simple integration program to calculate integrals.
! we will integrate using simpsons rule
real*8 tfun
external tfun
real*8 answer
real*8 pi

pi = 4.0d0*atan(1.0d0)

call qsimp(tfun,0.0d0,pi*0.5d0,answer)

write(6,*) answer

stop
end

real*8 function tfun(a)
implicit none
real*8 a
tfun = a*a
return
end



subroutine qsimp(func,a,b,s)
integer jmax
real*8 :: a,b,func,s,eps
external func
parameter( eps=1.0d-12, jmax = 20)
integer j
real*8 os,ost,st

ost = -1.0d30
os = 1.0d30
do j = 1, jmax
	call trapzd(func,a,b,st,j)
	s= (4.*st-ost)/3.
	if (abs(s-os).lt.eps*abs(os)) return
	write(6,*) s
	os = s
	ost = st
end do
write(6,*) 'too many steps in qsimp'
return
end


subroutine trapzd(func,a,b,s,n)
integer n
real*8 a,b,s,func
external func
integer it,j
real*8 del,sum,tnm,x
if (n.eq.1) then
	s=0.5*(b-a)*(func(a)+func(b))
else
	it=2**(n-2)
	tnm=it
	del=(b-a)/tnm
	x=a+0.5*del
	sum=0.0
	do j = 1 , it
		sum = sum+func(x)
		x=x+del
	end do
	s = 0.5*(s+(b-a)*sum/tnm)
end if

return
end
