!*************************************************************************
!* 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].                  *
!*                                                                       *
!*                                 F90 Release 1.0 By J-P Moreau, Paris. *
!*                                          (www.jpmoreau.fr)            *
!*************************************************************************
!Note: to link with kubgauss.f90.
Program TK3GAN

  integer i, Verfahren, Kantenteilung
  real*8 WertDerKubatur
  real*8 Px, Py, Qx, Qy, Rx, Ry

  !define summits of triangle
  Px =  0.d0;  Py =  0.d0
  Qx = 10.d0;  Qy =  0.d0
  Rx =  0.d0;  Ry = 10.d0

  !print header
  print *,' # Method     Number       Value      Error Code  '
  print *,'           Subtriangles   Integral    (must be 0) '
  print *,' -------------------------------------------------'
  !main integration loop
  do Verfahren = 1, 7
    !only methods 1,2,3,7 are implemented
    if (Verfahren < 4.or.Verfahren > 6) then
      do Kantenteilung = 1, 5
        !call appropriate Gaussian cubature method
        i = Kub3GauN (Px, Py, Qx, Qy, Rx, Ry, Verfahren,  &
                      Kantenteilung,                      &
                      WertDerKubatur)
        !print one line of results
        write(*,10) Verfahren, Kantenteilung**2, WertDerKubatur, i  
      end do
    end if
  end do
  print *,' -------------------------------------------------'
  stop

10 format('      ',I1,'         ',I2,' ',F18.10,'      ',I1)

END

!user defined function(x,y)
real*8 Function func(x,y)
  real*8 x, y
  func = Exp(-(x*x+y*y))
  return
end

! end of file tk3gan.f90