function adapt12(fun,a,b,acc,ncalls,err)
!
! Adaptive integrator which uses a combination of 1- and 2- point
! methods. Parameters:
! 	fun    -- name of the function to be integrated
!     a      -- lower limit
!     b      -- upper limit
!     acc    -- required relative accuracy
!     ncalls -- the number of times the integrand has been called
!
	parameter(max=900)
	dimension x(max),f(max)
	external fun

x(1)=a ; x(2)=(a+b)/2 ; x(3)=b
f(1)=fun(x(1)) ; f(2)=fun(x(2)) ; f(3)=fun(x(3))

i=1 ; ncalls=3 ; sum=0 ; err=0

10 	continue
	step=x(i+2)-x(i)
	a1=f(i+1)*step                               ! 1-point rule.
	a2=(f(i)+f(i+2))/2.*step                     ! 2-point rule.
	if(abs(a1-a2).gt.2*acc*(abs(a1)+abs(a2)))then  ! accurate enough?
		do j=ncalls,i+2,-1
			x(j+2)=x(j)                      ! If not, give space for
			f(j+2)=f(j)                      ! two more points.
		end do
		x(i+2)=x(i+1)
		f(i+2)=f(i+1)
		x(i+3)=x(i+2)+(x(i+4)-x(i+2))/2.       ! Add a point.
		f(i+3)=fun(x(i+3))
		x(i+1)=x(i)+(x(i+2)-x(i))/2.           ! Add a point.
		f(i+1)=fun(x(i+1))
		ncalls=ncalls+2
		if( ncalls.ge.max ) stop 'max steps reached'
		goto 10
	else                                         ! Step accepted.
		sum=sum+(2*a1+a2)/3                    ! Simpson's rule.
		err=err+abs(a1-a2)/2
		i=i+2                                  ! Integrate further.
		if( i.eq.ncalls ) goto 20              ! All points counted?
		goto 10
	end if	

20	continue
	
	adapt12=sum
	return
	end