'*********************************************************
'* 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.                   *
'* ----------------------------------------------------- *
'* 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.4596976938639075               *
'*                                                       *
'* ----------------------------------------------------- *
'* Ref.: "Mécanique des vibrations linéaires By          *
'*        M. Lalanne, P. Berthier and J. Der Hagopian,   *
'*        Masson, Paris, 1980" [BIBLI 16].               *
'*********************************************************
DEFDBL A-H, O-Z
DEFINT I-N

DIM H(10), A(10)

CLS
PRINT " INTEGRATION OF A REAL FUNCTION BY GAUSS"
PRINT "        F(X), F(X,Y) or F(X,Y,Z)"
PRINT
A$ = "deb"
INPUT " Number of variables (1 to 3): ", n2
PRINT
INPUT " How many Gauss points (2 to 10): ", n
PRINT
n1 = n - 1
ON n1 GOTO 380, 410, 460, 510, 580, 650, 740, 830, 940
380 A(1) = -.57735026919#
H(1) = 1#
GOTO 1040
410 A(1) = -.774596669241#
A(2) = 0#
H(1) = .555555555556#
H(2) = .888888888889#
GOTO 1040
460 A(1) = -.8611363115939999#
A(2) = -.339981043585#
H(1) = .347854845137#
H(2) = .652145154863#
GOTO 1040
510 A(1) = -.906179845939#
A(2) = -.538469310106#
A(3) = 1E-12
H(1) = .236926885056#
H(2) = .478628670499#
H(3) = .568888888889#
GOTO 1040
580 A(1) = -.9324695142029999#
A(2) = -.661209386466#
A(3) = -.238619186083#
H(1) = .171324492379#
H(2) = .360761573048#
H(3) = .467913934573#
GOTO 1040
650 A(1) = -.949107912343#
A(2) = -.741531185599#
A(3) = -.405845151377#
A(4) = 0#
H(1) = .129484966169#
H(2) = .279705391489#
H(3) = .381830050505#
H(4) = .417959183673#
GOTO 1040
740 A(1) = -.960289856497#
A(2) = -.796666477414#
A(3) = -.525532409916#
A(4) = -.183434642496#
H(1) = .10122853629#
H(2) = .222381034453#
H(3) = .313706645878#
H(4) = .362683783378#
GOTO 1040
830 A(1) = -.968160239508#
A(2) = -.836031107327#
A(3) = -.6133714327000001#
A(4) = -.324253423404#
A(5) = 0#
H(1) = .0812743883616#
H(2) = .180648160695#
H(3) = .260610696403#
H(4) = .31234707704#
H(5) = .330239355001#
GOTO 1040
940 A(1) = -.973906528517#
A(2) = -.865063366689#
A(3) = -.679409568299#
A(4) = -.433395394129#
A(5) = -.148874338982#
H(1) = .0666713443087#
H(2) = .149451349151#
H(3) = .219086362516#
H(4) = .26926671931#
H(5) = .295524224715#
1040 FOR i = 1 TO INT(n / 2)
  j = n + 1 - i
  A(j) = -A(i)
  H(j) = H(i)
NEXT i
PRINT
INPUT " Minimum value of X: ", x1
INPUT " Maximum value of X: ", x2
PRINT
IF n2 - 1 = 0 THEN GOTO 1280
INPUT " Minimum value of Y: ", y1
INPUT " Maximum value of Y: ", y2
PRINT
IF n2 - 2 = 0 THEN GOTO 1280
INPUT " Minimum value of Z: ", z1
INPUT " Maximum value of Z: ", z2
1280 xi9 = 0
FOR i = 1 TO n
  IF n2 - 1 = 0 THEN GOTO 1340
  FOR j = 1 TO n
    IF n2 - 2 = 0 THEN GOTO 1340
	FOR k = 1 TO n
1340      a1 = -(x1 - x2) / 2#
	  b1 = (x1 + x2) / 2#
	  x = a1 * A(i) + b1
	  IF n2 - 1 = 0 THEN GOTO 1450
	  a2 = -(y1 - y2) / 2#
	  b2 = (y1 + y2) / 2#
	  y = a2 * A(j) + b2
	  IF n2 - 2 = 0 THEN GOTO 1450
	  a3 = -(z1 - z2) / 2#
	  b3 = (z1 + z2) / 2#
	  z = a3 * A(j) + b3
1450      GOSUB 2000 'call function F
	  ON n2 GOTO 1470, 1490, 1510
1470      xi9 = xi9 + F * H(i) * a1
	  GOTO 1540
1490      xi9 = xi9 + F * H(i) * H(j) * a1 * a2
	  GOTO 1530
1510      xi9 = xi9 + F * H(i) * H(j) * H(k) * a1 * a2 * a3
	NEXT k
1530  NEXT j
1540 NEXT i
PRINT
PRINT " Value of integral = "; xi9
PRINT
END

2000 'define here function F
  F = SIN(x)
RETURN

'end of file tgauss.bas