!**************************************************************
!* Program to demonstrate integration of a real function F(x) *
!*                    by Romberg's method                     *
!* ---------------------------------------------------------- *
!* REFERENCE: "Mathematiques en Turbo-Pascal (Part 1)  By     *
!*             Marc Ducamp and Alain Reverchon, Eyrolles,     *
!*             Paris, 1987" [BIBLI 03].                       *
!*                              F90 version by J-P Moreau     *
!*                                  (www.jpmoreau.fr)         *
!* ---------------------------------------------------------- *
!* SAMPLE RUN:                                                *
!* (Integrate sin(x) from x=0 to x=1)                         *
!*                                                            *
!* Integral of a function F(X) by Romberg's method            *
!*                                                            *
!* Input begin and end values for x variable:                 *
!*                                                            *
!*         X0 = 0                                             *
!*         X1 = 1                                             *
!*                                                            *
!* Input desired precision: 1e-10                             *
!*                                                            *
!*                                                            *
!* Value of integral   :  0.459697694131851                   *
!*                                                            *
!* Obtained precision  :  9.599082639866197E-011              *
!*                                                            *
!* Number of iterations:  4                                   *
!*                                                            *
!**************************************************************
PROGRAM TEST_ROMBERG

real*8  x0,        &   ! begin x value       
        x1,        &   ! end   x value       
        prec,      &   ! desired precision   
        integral,  &   ! result of integral  
        obtprec        ! obtained precision

integer niter,     &   ! number of actual iterations
        itermin,   &   ! minimal number of iterations
	itermax        ! maximal number of iterations


  print *,' '
  print *,'Integral of a function F(X) by Romberg''s method'
  print *,' '
  print *,'Input begin and end values for x variable:'
  print *,' '
  write(*,"('         X0 = ')",advance='no')
  read  *, x0
  write(*,"('         X1 = ')",advance='no')
  read  *, x1
  print *,' '
  write(*,"(' Input desired precision: ')",advance='no') 
  read  *, prec

  itermin=1; itermax=50

  ! call Romberg subroutine
  integral = RombergIntegral(x0,x1,prec,obtprec,niter,itermin,itermax) 

  print *,' '
  print *,' '
  print *,'Value of integral   : ', integral
  print *,' '
  print *,'Obtained precision  : ', obtprec

  write(*,10)  niter

  stop

10 format(/' Number of iterations: ',i3//)

END !of main program


  ! Given function to integrate
  real*8 Function FUNC(x)
    real*8 x
    FUNC = dsin(x)
  end


  !*******************************************************
  !* Integral of a function FUNC(X) by Romberg's method  *
  !* --------------------------------------------------- *
  !* INPUTS:                                             *
  !*          a      begin value of x variable           *
  !*          b      end value of x variable             *
  !*       prec      desired precision                   *
  !*    itermin      minimum number of ietrations        *
  !*    itermax      maximum number of iterations        *
  !*                                                     *
  !* OUTPUTS:                                            *
  !*    obtprec      obtained precision for integral     *
  !*      niter      number of iterations done           *
  !*   integral      the integral of FUNC(X) from a to b *
  !*                                                     *
  !*******************************************************
real*8 Function RombergIntegral(a,b,prec,obtprec,niter,itermin,itermax)
    parameter(MAXITER = 15)
    real*8 a,b,prec,obtprec,FUNC
    real*8 t(MAXITER,MAXITER)
    real*8 pas,r,s,ta
    if (itermax>MAXITER)  itermax=MAXITER
    r = FUNC(a)
    ta = (r + FUNC(b)) / 2.d0
    niter=0
    pas=b-a
    t(0,0)=ta*pas
100 niter=niter+1
    pas=pas/2.d0
    s=ta
    do i=1, 2**niter-1
      s = s + FUNC(a+pas*i)
    end do
    t(0,niter)=s*pas
    r=1.d0
    do i=1, niter
      r=r*4.d0
      j=niter-i
      t(i,j)=(r*t(i-1,j+1) - t(i-1,j))/(r-1.d0)
    end do
    obtprec = DABS(t(niter,0) - t(niter-1,0))
    if (niter > itermax)  goto 200
    if (niter<itermin)  goto 100
    if (obtprec>prec)   goto 100
200 RombergIntegral = t(niter,0)
    return
end

! End of file tromberg.f90