!*********************************************************
!* Integration by Gauss method of a real function F=F(X) *
!* or F=F(X,Y) or F=F(X,Y,Z). The integral is calculated *
!* by using from 2 to 10 Gauss points.                   *
!* ----------------------------------------------------- *
!* Ref.: "Mécanique des vibrations linéaires By          *
!*        M. Lalanne, P. Berthier and J. Der Hagopian,   *
!*        Masson, Paris, 1980" [BIBLI 16].               *
!* ----------------------------------------------------- *
!* SAMPLE RUN:                                           *
!*                                                       *
!* (Integrate F=sin(x) from x=0 to x=1).                 *
!*                                                       *
!* INTEGRATION OF A REAL FUNCTION BY GAUSS               *
!*        F(X), F(X,Y) or F(X,Y,Z)                       *
!*                                                       *
!* Number of variables (1 to 3): 1                       *
!*                                                       *
!* How many Gauss points (2 to 10): 4                    *
!*                                                       *
!* Minimum value of X: 0                                 *
!* Maximum value of X: 1                                 *
!*                                                       *
!* Value of integral =    0.4596977                      *
!*                                                       *
!*                    F90 version by J-P Moreau, Paris.  *
!*                           (www.jpmoreau.fr)           *
!*********************************************************
Program Test_Gauss

!Labels: 10,20,30,40

real    A(10), H(10)
real    x,x1,x2,y,y1,y2,z,z1,z2
integer i,j,k,n,n1,n2
real    a1,a2,a3,b1,b2,b3,xi9


  print *,' '
  print *,' INTEGRATION OF A REAL FUNCTION BY GAUSS'
  print *,'        F(X), F(X,Y) or F(X,Y,Z)'
  print *,' '
  write(*,50,advance='no'); read *, n2
  print *,' '
  write(*,60,advance='no'); read *, n
  print *,' '
  n1 = n - 1
  Select Case(n1)
    case(1)
        A(1) = -0.57735026919
        H(1) = 1.0
    case(2)
        A(1) = -0.774596669241
        A(2) = 0
        H(1) = 0.555555555556
        H(2) = 0.888888888889
    case(3)
        A(1) = -0.8611363115939999
        A(2) = -0.339981043585
        H(1) = 0.347854845137
        H(2) = 0.652145154863
    case(4)
        A(1) = -0.906179845939
        A(2) = -0.538469310106
        A(3) = 1e-12
        H(1) = 0.236926885056
        H(2) = 0.478628670499
        H(3) = 0.568888888889
    case(5)
        A(1) = -0.9324695142029999
        A(2) = -0.661209386466
        A(3) = -0.238619186083
        H(1) = 0.171324492379
        H(2) = 0.360761573048
        H(3) = 0.467913934573
    case(6)
        A(1) = -0.949107912343
        A(2) = -0.741531185599
        A(3) = -0.405845151377
        A(4) = 0
        H(1) = 0.129484966169
        H(2) = 0.279705391489
        H(3) = 0.381830050505
        H(4) = 0.417959183673
    case(7)
        A(1) = -0.949107912343
        A(2) = -0.741531185599
        A(3) = -0.405845151377
        A(4) = 0
        H(1) = 0.129484966169
        H(2) = 0.279705391489
        H(3) = 0.381830050505
        H(4) = 0.417959183673
    case(8)
        A(1) = -0.960289856497
        A(2) = -0.796666477414
        A(3) = -0.525532409916
        A(4) = -0.183434642496
        H(1) = 0.10122853629
        H(2) = 0.222381034453
        H(3) = 0.313706645878
        H(4) = 0.362683783378
    case(9)
        A(1) = -0.968160239508
        A(2) = -0.836031107327
        A(3) = -0.6133714327000001
        A(4) = -0.324253423404
        A(5) = 0
        H(1) = 0.0812743883616
        H(2) = 0.180648160695
        H(3) = 0.260610696403
        H(4) = 0.31234707704
        H(5) = 0.330239355001
    case(10)
        A(1) = -0.973906528517
        A(2) = -0.865063366689
        A(3) = -0.679409568299
        A(4) = -0.433395394129
        A(5) = -0.148874338982
        H(1) = 0.0666713443087
        H(2) = 0.149451349151
        H(3) = 0.219086362516
        H(4) = 0.26926671931
        H(5) = 0.295524224715
  End Select
  do i = 1, n/2
    j = n + 1 - i
    A(j) = -A(i)
    H(j) = H(i)
  end do
  print *,' '
  write(*,70,advance='no'); read *, x1
  write(*,71,advance='no'); read *, x2
  print *,' '
  if (n2 - 1 > 0) then
    write(*,80,advance='no'); read *, y1
    write(*,81,advance='no'); read *, y2
    print *,' '
  end if
  if (n2 - 2 > 0) then
    write(*,90,advance='no'); read *, z1
    write(*,91,advance='no'); read *, z2
  end if
  xi9 = 0
  do i = 1, n
    IF (n2 - 1 == 0) GOTO 10
    do j = 1, n
      IF (n2 - 2 == 0) GOTO 10
  	  do k = 1, n
10      a1 = -(x1 - x2) / 2
	    b1 = (x1 + x2) / 2
	    x = a1 * A(i) + b1
	    IF (n2 - 1 == 0) GOTO 20
	    a2 = -(y1 - y2) / 2
	    b2 = (y1 + y2) / 2
	    y = a2 * A(j) + b2
	    IF (n2 - 2 == 0) GOTO 20
	    a3 = -(z1 - z2) / 2
	    b3 = (z1 + z2) / 2
	    z = a3 * A(j) + b3
20      Select Case(n2)
          Case(1)
            xi9 = xi9 + F(x,y,z) * H(i) * a1
	        GOTO 40
          Case(2)
            xi9 = xi9 + F(x,y,z) * H(i) * H(j) * a1 * a2
	        GOTO 30
          Case(3)
		    xi9 = xi9 + F(x,y,z) * H(i) * H(j) * H(k) * a1 * a2 * a3
        End Select
	  end do
30   end do
40 end do 
  print *,' '
  print *,' Value of integral = ', xi9
  print *,' '
  print *,' '
  pause ' Any key to exit...'

50 format('  Number of variables (1 to 3): ')
60 format('  How many Gauss points (2 to 10): ')
70 format('  Minimum value of X: ')
71 format('  Maximum value of X: ')
80 format('  Minimum value of Y: ')
81 format('  Maximum value of Y: ')
90 format('  Minimum value of Z: ')
91 format('  Maximum value of Z: ')


END

!define here function F
real Function F(x,y,z)
real x,y,z
  F = SIN(x)
end

!end of file tgauss.f90