/************************************************************************
* Test programm for cubature over triangles via summed Gaussian n-point *
* formula                                                               *
* --------------------------------------------------------------------- *
* SAMPLE RUN:                                                           *
* (Integrate function EXP(-(x*x + y*y)) over triangle defined by three  *
*  points (0,0), (10,0) and (0,10). )                                   *
*                                                                       *
* # Method     Number       Value      Error Code                       *
*           Subtriangles   Integral    (must be 0)                      *
* -------------------------------------------------                     *
*     1          1       0.0000000112      0                            *
*     1          4       0.0483240045      0                            *
*     1          9       0.4706075316      0                            *
*     1         16       0.7913533084      0                            *
*     1         25       0.8814598010      0                            *
*     2          1       0.0000000000      0                            *
*     2          4       0.0120496973      0                            *
*     2          9       0.2538493244      0                            *
*     2         16       0.5593099250      0                            *
*     2         25       0.7075499085      0                            * 
*     3          1       0.0644320023      0                            *
*     3          4       1.0390297430      0                            *
*     3          9       1.0188794444      0                            *
*     3         16       0.8590974900      0                            *
*     3         25       0.7974378946      0                            *
*     7          1       0.8091876107      0                            *
*     7          4       0.9654995491      0                            *
*     7          9       0.7922397459      0                            *
*     7         16       0.7849640612      0                            *
*     7         25       0.7862848602      0                            *
* -------------------------------------------------                     *
*                                                                       *
* --------------------------------------------------------------------- *
* 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 kubgauss.cpp.                                       */ 
#include <basis.h>          //for REAL, etc.

int Kub3GauN (              //see file kubgauss.cpp
              REAL, REAL,
              REAL, REAL,
              REAL, REAL,
              int,   int,
              REAL f(REAL,REAL),
              REAL *
             );


void main()  {
  int    i, Verfahren, Kantenteilung;
  REAL   WertDerKubatur;
  REAL   Px, Py, Qx, Qy, Rx, Ry;
  REAL   exp2 (REAL, REAL);                // Function from F_Beisp.cpp

  Px = ZERO;  Py = ZERO;
  Qx =  TEN;  Qy = ZERO;
  Rx = ZERO;  Ry =  TEN;

  //print header
  printf(" # Method     Number       Value      Error Code\n");
  printf("           Subtriangles   Integral    (must be 0)\n");
  printf(" -------------------------------------------------\n");
  for (Verfahren = 1; Verfahren <= 7; Verfahren++)
    if (Verfahren < 4 || Verfahren > 6)
    {
      for (Kantenteilung = 1; Kantenteilung <= 5; Kantenteilung++)
      {
        printf ("     %d", Verfahren);   
        printf ("         %2d", Kantenteilung * Kantenteilung);

        i = Kub3GauN (Px, Py, Qx, Qy, Rx, Ry, Verfahren,
                      Kantenteilung,
                      exp2,                // Function from F_Beisp.cpp
                      &WertDerKubatur);

        printf (" %18.10f      %d\n", WertDerKubatur, i);

      } /* for Kantenteilung */
    } /* if Verfahren */
    printf(" -------------------------------------------------\n");
  
}

//end of file tk3gan.cpp