!************************************************************************
!*     Test program: Numerical differentation according to  Romberg     *
!* -------------------------------------------------------------------- *
!* SAMPLE RUN:                                                          *
!* (Find the first derivative of function f(x)=1/x, for x=0.12).        *
!*                                                                      *
!*  Numerical differentation according to Romberg                       *
!*                                                                      *
!*  Test function   f(x) = 1/x                                          *
!*                                                                      *
!*  Put in x-value at which you want to evaluate derivative: 0.12       *
!*  Put in desired accuracy: 1e-10                                      *
!*  Maximal number of columns in Romberg scheme: 4                      *
!*  Starting step size: 0.005                                           *
!*                                                                      *
!*  Results:                                                            *
!*                                                                      *
!*  er_app=   8.891731795301894E-011                                    *
!*  res   =   -69.4444444444433                                         *
!*  nend  =            3                                                *
!*  hend  =   6.250000000000000E-004                                    *
!*                                                                      *
!* -------------------------------------------------------------------- *
!* "Numerical Algorithms with C, By Gisela Engeln-Muellges              *
!*  and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].                 *
!*                                                                      *
!*                                F90 Release 1.0 By J-P Moreau, Paris. *
!*                                         (www.jpmoreau.fr)            *
!************************************************************************
Program Test_Difrom

integer error, n, nend
real*8  x,prec,h, res,schaetz, hend

   print *,' '
   print *,' Numerical differentation according to Romberg'
   print *,' '
   print *,' Test function   f(x) = 1/x'
   print *,' '

   write(*,10,advance='no'); read *, x
   write(*,20,advance='no'); read *, prec
   write(*,30,advance='no'); read *, n
   write(*,40,advance='no'); read *, h

   call difrom(x,prec,n,h,res,schaetz,nend,hend,error)

   call aus(error,schaetz,res,nend,hend)

   stop

10 format('  Put in x-value at which you want to evaluate derivative: ')
20 format('  Put in desired accuracy: ')
30 format('  Maximal number of columns in Romberg scheme: ')
40 format('  Starting step size: ')

END

!print results
Subroutine aus(error, schaetz, res, nend, hend)
integer error, nend
real*8 schaetz,res,hend
  if (error.ne.1) then
    print *,' '
    print *,' Results:'
    print *,' '
    print *,' er_app= ', schaetz
    print *,' res   = ', res
    print *,' nend  = ', nend
    print *,' hend  = ', hend
	print *,' '
  end if
return
End

!test function to derivate
real*8 Function func(x)
real*8 x
  func = 1.d0 / x
return
end

Subroutine difrom(x0, eps, n, h, res, er_app, nend, hend, error)
real*8 x0,eps,h,res,er_app,hend
integer error
!*************************************************************************
!*  Computes an approximation for the first derivative of func at x0     *
!*  using the ROMBERG method.                                            *
!* --------------------------------------------------------------------- *
!*  Input parameters :                                                   *
!*                                                                       *
!*    real*8 func(real*8)  (External) Name of function to be             *
!*                         differentiated.                               *
!*    real*8 x0            value of abscissa at which derivative is to   *
!*                         be found.                                     *
!*    real*8 eps           desired accuracy.                             *
!*    integer n            max. number of columns in the Romberg scheme  *
!*                         (n > 1).                                      *
!*    real*8 h             initial step size.                            *
!*                                                                       *
!*  Output parameters:                                                   *
!*                                                                       *
!*    real*8 res           approximate derivative                        *
!*    real*8 er_app        error estimate for  res                       *
!*    integer nend         number of columns actually used in scheme     *
!*    real*8 hend          final step size                               *
!*                                                                       *
!*  Return value (integer error):                                        *
!*                                                                       *
!*    0:   no error: er_app < eps                                        *
!*    1:   n < 1  or  eps <= 0  or  h < MACH_EPS                         *
!*    2:   desired accuracy not reached after n steps                    *
!*    3:   step size drooped below  MACH_EPS                             *
!*    4:   lack of sufficient memory (not used here)                     *
!*                                                                       *
!*************************************************************************
parameter(NMAX=10, MACH_EPS = 1.2d-16)
real*8 h2, d1, d2, func
real*8 d(0:NMAX)

  if (n <= 1.or.eps <= 0.0.or.h < MACH_EPS) then    !check input
    error=1
    return
  end if

  h2 = 2.d0 * h
  d(0) = (func(x0 + h) - func(x0 - h)) / h2

!************************************************************************
!* This loop runs until the maximum of Romberg rows is filled or until  *
!* the step size use drops below the machine constant or if the desired *
!* accuracy is reached.                                                 *
!************************************************************************
  error=2
  do j = 1, n - 1
    d(j) = 0.d0
    d1   = d(0)
    h2   = h
    h    = h * 0.5d0
    nend = j

    if (h < MACH_EPS) then
      error = 3                  !step size less than machine constant
      goto 10
    end if

    d(0) = (func(x0 + h) - func(x0 - h)) / h2

    m=4
    do i = 1, j
      d2 = d(i)
      d(i) = (m * d(i-1) - d1) / (m-1)
      d1 = d2
      m = m * 4
    end do

    er_app = DABS(d(j) - d(j-1))

    if (er_app < eps) then
      error = 0                          !desired accuracy reached
      goto 10
    end if
  end do

10 res = d(nend)                         !save final values
   hend = h

Return

END

!end of file tdifrom.f90