!*****************************************************************
!* Test program for cubature over rectangles using Newton-Cotes  *
!* ------------------------------------------------------------- *
!* SAMPLE RUN:                                                   *
!* (Integrate function EXP(-(x*x + y*y)) in rectangle defined by *
!*  x1=0, x2=10, y1=0, y2=10).                                   *
!*                                                               *
!* -----------------------------------------------------------   *
!*   #    ERROR     RESULT        ERROR     NUMBER OF FUNCTION   *
!* METHOD CODE                  ESTIMATE        EVALUATIONS      *
!* -----------------------------------------------------------   *
!* Without error estimation:                                     *
!*   1    ( 0)   0.7855606650                        121         *
!*   2    ( 0)   0.7853439999                        441         *
!*   3    ( 0)   0.7853778519                        961         *
!*   4    ( 0)   0.7854017744                       1681         *
!*   5    ( 0)   0.7854001104                       2601         *
!*   6    ( 0)   0.7853979700                       3721         *
!*   7    ( 0)   0.7853980474                       5041         *
!* With error estimation:                                        *
!*   1    ( 0)   0.7853981634   0.541672E-04         562         *
!*   2    ( 0)   0.7853981634  -0.361090E-05        2122         *
!*   3    ( 0)   0.7853981634  -0.135410E-05        4682         *
!*   4    ( 0)   0.7853981634   0.573169E-07        8242         * 
!*   5    ( 0)   0.7853981634   0.309050E-07       12802         *
!*   6    ( 0)   0.7853981634  -0.758605E-09       18362         *
!*   7    ( 0)   0.7853981634  -0.454919E-09       24922         *
!* -----------------------------------------------------------   *
!*                                                               *
!* Reference:                                                    *
!* "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)             *
!*****************************************************************
! To link with Kubnec.f90
Program TKUBNEC

IMPLICIT NONE

  real*8  a, b, c, d, W, F
  integer i, Nx, Ny, method, nF 
  integer mSCH                          !with error estimation, if <>0

  integer Kub4NeCn

       !Function to integrate func(x,y) is declared in file kubnec.f90

  a = 0.d0; c = 0.d0           !define limits of integration rectangle
  b = 10.d0; d = 10.d0

  Nx=10; Ny=10 

  !print header
  print *, '-----------------------------------------------------------'
  print *, '  #    ERROR     RESULT        ERROR     NUMBER OF FUNCTION'
  print *, 'METHOD CODE                  ESTIMATE        EVALUATIONS   '
  print *, '-----------------------------------------------------------'
  
  !main loop
  do mSCH = 0, 1
    if (mSCH.ne.0) then
      print *, 'With error estimation:'
    else
      print *, 'Without error estimation:'
    end if

    !integrate with 7 methods
    do method = 1, 7

      ! return value i of Kub4NeCn must be zero 
      i = Kub4NeCn (a, b, Nx, c, d, Ny, method, W, mSCH, F)  !see kubnec.f90

      nF = (Nx * method + 1) * (Ny * method + 1)
      if (mSCH.ne.0)  nF = nF + (2 * Nx * method + 1)*(2 * Ny * method + 1)
      
	  if (mSCH.ne.0) then   !with error estimation 
	    write(*,10) method, i, W, F, nF
      else                  !without error estimation
	    write(*,20) method, i, W, nF
      end if

    end do
  end do
  print *, '-----------------------------------------------------------'

  stop

10 format(I4,'    (',I2,') ',F13.10,'   ',E13.6,'      ',I6)
20 format(I4,'    (',I2,') ',F13.10,'                      ',I6)

END

!end of file tkubnec.f90