/*

  Programa que integra uma funcao por varios metodos 
  numericos.  
  Adaptado de Numerical Methods Lectures Notes
              Natural Sciences Tripos 1B
              Stuart Dalziel

*/

#include <stdio.h>
#include <math.h>

/* funcao a ser integrada */
double f(double x) {
   return (sin(x));
}


double trapezio(double x0, double x1, int nx){ 
  double xa,xb,fa,fb;
  int i;
  double dx = (x1 - x0)/nx;
  double soma = 0.0;
  for (i=0; i<nx; i++){
     xa = x0 + i*dx;
     xb = x0 + (i+1)*dx;  
     fa = f(xa);
     fb = f(xb);
     soma = soma + fa + fb;
   }
   soma = soma*dx/2.0;
   return soma;
}

double pontomedio(double x0, double x1, int nx){
  double xa,fa;
  int i;
  double dx = (x1 - x0)/nx;
  double soma = 0.0;
  for (i=0; i<nx; i++){
     xa = x0 + ((double)i + 0.5)*dx;
     fa = f(xa);
     soma = soma + fa;
   }
   soma = soma*dx;
   return soma;
}

double simpson(double x0, double x1, int nx){ 
  double xa,xb,xc,fa,fb,fc;
  int i;
  double dx = (x1 - x0)/nx;
  double soma = 0.0;
  for (i=0; i<nx; i=i+2){
     xa = x0 + i*dx;
     xb = x0 + (i+1)*dx;  
     xc = x0 + (i+2)*dx;
     fa = f(xa);
     fb = f(xb);
     fc = f(xc);
     soma = soma + fa + 4*fb +fc;
    }
   soma = soma*dx/3.0;
   return soma;
}

double quad(double x0, double x1, int nx){
  double xa,xb,fa,fb;
  int i;
  double dx = (x1 - x0)/nx;
  double dxl = dx*(0.5 - sqrt(3.0)/6.0);
  double dxr = dx*(0.5 + sqrt(3.0)/6.0);
  double soma = 0.0;
  for (i=0; i<nx; i++){
     xa = x0 + i*dx + dxl;
     xb = x0 + i*dx + dxr;      
     fa = f(xa);
     fb = f(xb);
     soma = soma + fa + fb;
   }
   soma = soma*dx/2.0;
   return soma;
}


int main()
{  

  double PI = 4.0*atan2(1.0,1.0);

  double x0 = 0.0;
  double x1 = PI;
  double exato = 2.0;
  double integ;
  
  int nx;     /* numero de divisoes no dominio 
		 vamos de 1 ate 2^10 */
  int j;
  int jMax = 12;   

  char form[60] = "nx = %4d  integ = %18.15f exato-integ = %18.15f\n" ;

  printf("Regra do Trapezio \n");
  nx = 1;
  for (j=1; j<=jMax; j++){
   integ = trapezio(x0,x1,nx);
   printf(form, 
           nx, integ,  exato-integ);
   nx = nx*2;
  }

  printf("Regra do ponto medio \n");
  nx = 1;
  for (j=1; j<=jMax; j++){
   integ = pontomedio(x0,x1,nx);
   printf( form, 
           nx, integ,  exato-integ);
   nx=nx*2;
  }

  printf("Metodo de Simpson \n");
  nx=1;
  for (j=1; j<=jMax; j++){
   integ = simpson(x0,x1,nx);
   printf(form, 
           nx, integ,  exato-integ);
   nx = nx*2;
  }  

  nx = 1;
  printf("Quadratura de Gauss \n");
  for (j=1; j<=jMax; j++){
   integ = quad(x0,x1,nx);
   printf(form, 
           nx, integ,  exato-integ);
   nx = nx*2;
  }  
  return 0;
}