/********************************************************
* 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.45969769                      *
*                                                       *
*                   C++ version by J-P Moreau, Paris.   *
*                          (www.jpmoreau.fr)            *
********************************************************/
#include <stdio.h>
#include <math.h>

// Labels: e10,e20,e30,e40

double A[11], H[11];
double x,x1,x2,y,yy1,y2,z,z1,z2;
int    i,j,k,n,n1,n2;
double a1,a2,a3,b1,b2,b3,xi9;

//define here function F
double F(double x,double y,double z) { 
  return sin(x);
}


void main() {
  printf("\n INTEGRATION OF A REAL FUNCTION BY GAUSS\n");
  printf("        F(X), F(X,Y) or F(X,Y,Z)\n");
  printf("\n Number of variables (1 to 3): "); scanf("%d",&n2);
  printf("\n How many Gauss points (2 to 10): "); scanf("%d",&n);
  printf("\n");
  n1 = n - 1;
  switch(n1) {
    case 1:A[1] = -0.57735026919;
           H[1] = 1.0;
		   break;
    case 2:A[1] = -0.774596669241;
           A[2] = 0.0;
           H[1] = 0.555555555556;
           H[2] = 0.888888888889;
		   break;
    case 3:A[1] = -0.8611363115939999;
           A[2] = -0.339981043585;
           H[1] = 0.347854845137;
           H[2] = 0.652145154863;
           break;
    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;
           break;
    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;
           break;   
    case 6:A[1] = -0.949107912343;
           A[2] = -0.741531185599;
           A[3] = -0.405845151377;
           A[4] = 0.0;
           H[1] = 0.129484966169;
           H[2] = 0.279705391489;
           H[3] = 0.381830050505;
           H[4] = 0.417959183673;
           break;
    case 7:A[1] = -0.949107912343;
           A[2] = -0.741531185599;
           A[3] = -0.405845151377;
           A[4] = 0.0;
           H[1] = 0.129484966169;
           H[2] = 0.279705391489;
           H[3] = 0.381830050505;
           H[4] = 0.417959183673;
           break;
    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;
           break;
    case 9:A[1] = -0.968160239508;
           A[2] = -0.836031107327;
           A[3] = -0.6133714327000001;
           A[4] = -0.324253423404;
           A[5] = 0.0;
           H[1] = 0.0812743883616;
           H[2] = 0.180648160695;
           H[3] = 0.260610696403;
           H[4] = 0.31234707704;
           H[5] = 0.330239355001;
           break;
    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;
  }
  for (i = 1; i<=n/2; i++) {
    j = n + 1 - i;
    A[j] = -A[i];
    H[j] = H[i];
  }
  
  printf("\n Minimum value of X: "); scanf("%lf", &x1);
  printf("\n Maximum value of X: "); scanf("%lf", &x2);
  printf("\n");
  if (n2 - 1 > 0) {
    printf("\n Minimum value of Y: "); scanf("%lf", &yy1);
    printf("\n Maximum value of Y: "); scanf("%lf", &y2);
    printf("\n");
  }
  if (n2 - 2 > 0) {
    printf("\n Minimum value of Z: "); scanf("%lf", &z1);
    printf("\n Maximum value of Z: "); scanf("%lf", &z2);
  }
  xi9 = 0;
  for (i = 1; i<=n; i++) {
    if (n2 - 1 == 0) goto e10;
    for (j = 1; j<=n; j++) {
      if (n2 - 2 == 0) goto e10;
	  for (k=1; k<=n; k++) {
e10:    a1 = -(x1 - x2) / 2;
	    b1 = (x1 + x2) / 2;
	    x = a1 * A[i] + b1;
	    if (n2 - 1 == 0) goto e20;
	    a2 = -(yy1 - y2) / 2;
	    b2 = (yy1 + y2) / 2;
	    y = a2 * A[j] + b2;
	    if (n2 - 2 == 0) goto e20;
	    a3 = -(z1 - z2) / 2;
	    b3 = (z1 + z2) / 2;
	    z = a3 * A[j] + b3;
e20:    switch(n2) {
          case 1:xi9 += F(x,y,z) * H[i] * a1;
	             goto e40;
                 break;
          case 2:xi9 += F(x,y,z) * H[i] * H[j] * a1 * a2;
	             goto e30;
                 break;
          case 3:xi9 += F(x,y,z) * H[i] * H[j] * H[k] * a1 * a2 * a3;
		}
	  }
e30: ;}   
e40:;} 
  printf("\n Value of integral = %12.8f\n\n", xi9);

}

// end of file tgauss.cpp