/****************************************************************
* 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   5.417e-005           562         *
*   2    ( 0)   0.7853981634   -3.611e-006         2122         *
*   3    ( 0)   0.7853981634   -1.354e-006         4682         *
*   4    ( 0)   0.7853981634   5.732e-008          8242         * 
*   5    ( 0)   0.7853981634   3.090e-008         12802         *
*   6    ( 0)   0.7853981634   -7.586e-010        18362         *
*   7    ( 0)   0.7853981634   -4.549e-010        24922         *
* -----------------------------------------------------------   *
*                                                               *
* Reference:                                                    *
* "Numerical Algorithms with C, By Gisela Engeln-Muellges       *
*  and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].          *
*                                                               *
*                         C++ Release 1.0 By J-P Moreau, Paris. *
*                                  (www.jpmoreau.fr)            *
*****************************************************************
Note: to link with files: basis_r.cpp, kubnec.cpp, vmblock.cpp
      and f_beisp.cpp.
--------------------------------------------------------------- */
#include <basis.h>  //For REAL, etc.


int Kub4NeCn (      //See kubnec.cpp
              REAL, REAL, int,
              REAL, REAL, int,
              int,
              REAL  f (REAL,REAL),
              REAL*,
              int,
              REAL*
             );


void main()  {

  REAL   a, b, c, d;
  REAL   W, F;
  int    Nx = 10, Ny = 10;
  int    method;
  int    mSCH;                      // with error estimation, if <> 0
  REAL   exp2 (REAL,REAL);          // Function declared in f_beisp.cpp

  a = c = ZERO;                     //define limits of integration rectangle
  b = d = TEN;

  //print header
  printf("-----------------------------------------------------------\n");
  printf("  #    ERROR     RESULT        ERROR     NUMBER OF FUNCTION\n");
  printf("METHOD CODE                  ESTIMATE        EVALUATIONS\n");
  printf("-----------------------------------------------------------\n");
  
  //main loop
  for (mSCH = 0; mSCH <= 1; mSCH++) {
    if (mSCH) printf("With error estimation:\n");
	else printf("Without error estimation:\n");
    //integrate with 7 methods
	for (method = 1; method <= 7; method++)
    {
      int i, nF;

      i = Kub4NeCn (a, b, Nx, c, d, Ny, method,
                    exp2,
                    &W, mSCH, &F);

      printf ("%3d    (%2d)  %13.10"LZP"f",   method, i, W);

      if (mSCH) printf ("   %10.3"LZP"e", F);
      else      printf ("%13c",' ');

      nF = (Nx * method + 1) * (Ny * method + 1);

      if (mSCH) nF += (2 * Nx * method + 1)*(2 * Ny * method + 1);
      if (F>0.0) printf ("       %7d\n", nF);
	  else printf ("       %6d\n", nF);
    }
  }
  printf("-----------------------------------------------------------\n");

  return;
}

//end of file tkubnec.cpp