/***************************************************************************
 *   Copyright (C) 2009 by Douglas Castro   *
 *   douglas@ime.unicamp.br   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include "interp_integ.h"


double integra(int d, double xmin[], double xmax[], Funcao *f)
{
  int m = 5; //ordem de integracao.
	// raizes qsi[j] e pesos w[j] da integracao de gauss 
  double qsi[ ] = {-0.906179846, -0.53846931, 0.0, 0.53846931, 0.906179846};
	double w[ ] = { 0.236926885, 0.47862867, 0.568888889, 0.47862867, 0.236926885};
	// comeca a ingtegracao..
  double J = 0.0;
  double A, B;
	double h1, h2;
  int i, j, k;
  
  switch(d)
  {
    case 1:
      // dominio de integracao [A,B]
      A = xmin[0];
      B = xmax[0];

		  // variaveis auxiliares
      h1 = (B-A)/2.0;
      h2 = (B+A)/2.0;
      for(i=0;i<m;i++)
      {
        double x = h1*qsi[i] + h2;
        double y,z;
        double Q = f(d,x,y,z); // aqui y,z serao descartados...isso f que ira fazer
        J += w[i]*Q;
      }
      J = h1*J;
      break;
    case 2:
      // dominio de integracao [A,B]x[C,D]
      A = xmin[0];
      B = xmax[0];

		  // variaveis auxiliares
      double h1 = (B-A)/2.0;
      double h2 = (B+A)/2.0;
      for(i=0;i<m;i++)
      {
        double JX = 0.0;
        double x = h1*qsi[i] + h2;
			  // caso C = C(x) e D = D(x) substitui xmin e xmax pelas funcoes 
        double C = xmin[1];
        double D = xmax[1];
        double k1 = (D-C)/2.0;
        double k2 = (D+C)/2.0;
        for(j=0;j<m;j++)
        {
          double y = k1*qsi[j] + k2;
          double z;
          double Q = f(d,x,y,z); // aqui z sera descartado...isso f que ira fazer
          JX += w[j]*Q;
        }
        J += w[i]*k1*JX;
      }
      J = h1*J;
      break;
    case 3:
      // dominio de integracao [A,B]x[C,D]x[E,F]
      A = xmin[0];
      B = xmax[0];
      h1 = (B - A)/2.0;
      h2 = (B + A)/2.0;
      for(i=0;i<m;i++)
      {
        double JX =0.0;
        double x = h1*qsi[i] + h2;
        double D = xmax[1]; // caso D= D(x)  troca xmax por D(x)
        double C = xmin[1]; // caso C= C(x)  troca xmin por C(x)
        double k1 = (D - C)/2.0;
        double k2 = (D + C)/2.0;
        for(j=0;j<m;j++)
        {
          double JY =0.0;
          double y = k1*qsi[j] + k2;
          double F = xmax[2]; // caso F= F(x)  troca xmax por F(x)
          double E = xmin[2]; // caso E= E(x)  troca xmin por E(x)
          double l1 = (D - C)/2.0;
          double l2 = (D + C)/2.0;
          for(k=0;k<m;k++)
          {
            double z = l1*qsi[k] + l2;
            double Q = f(d,x,y,z); //supondo que aqui sei que a f esta definida em R^3
            JY += w[k]*Q;
          }
          JX += w[j]*l1*JY;
        }
        J += w[i]*k1*JX;
      }
      J = h1*J;
      break;
    default:
      fprintf(stderr," Nao ha regra de integacao para essa dimensao =/");
  }
  return J;
}

double interpola_quad( double x0, double x1, double x2 )
{
	// preve media do filho esquerdo 
  return ( x1 - (x2 - x0)/8.0 );
	// preve media do filho direito 
// 	return ( x1 + (x2 - x0)/8.0 );
}

// regras de integracao...
/*  if(d==1)
{
		// dominio de integracao [A,B]
    double A = xmin[0];
    double B = xmax[0];

		// variaveis auxiliares
    double h1 = (B-A)/2.0;
    double h2 = (B+A)/2.0;
    int i, j;
    for(i=0;i<m;i++)
{
      double JX = 0.0;
      double x = h1*qsi[i] + h2;
      double y,z;
      double Q = f(d,x,y,z); // aqui y,z serao descartados...isso f que ira fazer
      J += w[i]*Q;
}
    J = h1*J;
}
  else
    if(d==2)
{
		// dominio de integracao [A,B]x[C,D]
    double A = xmin[0];
    double B = xmax[0];

		// variaveis auxiliares
    double h1 = (B-A)/2.0;
    double h2 = (B+A)/2.0;
    int i, j;
    for(i=0;i<m;i++)
{
      double JX = 0.0;
      double x = h1*qsi[i] + h2;
			// caso C = C(x) e D = D(x) substitui xmin e xmax pelas funcoes 
      double C = xmin[1];
      double D = xmax[1];
      double k1 = (D-C)/2.0;
      double k2 = (D+C)/2.0;
      for(j=0;j<m;j++)
{
        double y = k1*qsi[j] + k2;
        double z;
        double Q = f(d,x,y,z); // aqui z sera descartado...isso f que ira fazer
        JX += w[j]*Q;
}
      J += w[i]*k1*JX;
}
    J = h1*J;
}
  else if(d==3)
{
		// dominio de integracao [A,B]x[C,D]x[E,F]
    int i,j,k;
    double A = xmin[0];
    double B = xmax[0];
    double h1 = (B - A)/2.0;
    double h2 = (B + A)/2.0;
    for(i=0;i<m;i++)
{
      double JX =0.0;
      double x = h1*qsi[i] + h2;
      double D = xmax[1]; // caso D= D(x)  troca xmax por D(x)
      double C = xmin[1]; // caso C= C(x)  troca xmin por C(x)
      double k1 = (D - C)/2.0;
      double k2 = (D + C)/2.0;
      for(j=0;j<m;j++)
{
        double JY =0.0;
        double y = k1*qsi[j] + k2;
        double F = xmax[2]; // caso F= F(x)  troca xmax por F(x)
        double E = xmin[2]; // caso E= E(x)  troca xmin por E(x)
        double l1 = (D - C)/2.0;
        double l2 = (D + C)/2.0;
        for(k=0;k<m;k++)
{
          double z = l1*qsi[k] + l2;
          double Q = f(d,x,y,z); //supondo que aqui sei que a f esta definida em R^3
          JY += w[k]*Q;
}
        JX += w[j]*l1*JY;
}
      J += w[i]*k1*JX;
}
    J = h1*J;
}
  else
{
    fprintf(stderr," Nao ha regra de integacao para essa dimensao =/");
}
  return J;*/