/***************************************************************************
 *   Copyright (C) 2009 by Douglas Castro,,,   *
 *   douglas@douglas-laptop   *
 *                                                                         *
 ***************************************************************************/

#include "interp_integ.h"


double Qx(double ul,double ur)
{
  return -(ur-ul)/8.0;
}

double Qy(double ul,double ur)
{
  return -(ur-ul)/8.0;
}
  /* ulr left pra x e right pra y, alusao a u_{i-1,j+1}*/
double Qxy(double ull,double ulr,double url, double urr)
{
  return (urr-url-ulr+ull)/64.0;
}

double integra(double xmin, double xmax, double ymin, double ymax, Funcao *f)
{
	double qsi[20], w[20];
	// raizes qsi[j] e pesos w[j] da integracao de gauss 
	qsi[0] = -0.906179846;
	qsi[1] = -0.53846931;
	qsi[2] = 0.0;
	qsi[3] = 0.53846931;
	qsi[4] = 0.906179846;
	
	w[0] = 0.236926885;
	w[1] = 0.47862867;
	w[2] = 0.568888889;
	w[3] = 0.47862867;
	w[4] = 0.236926885;
	// comeca a ingtegracao..

	// dominio de integracao [A,B]x[C,D]
	double A = xmin, B = xmax;

	// variaveis auxiliares
	double h1 = (B-A)/2.0, h2 = (B+A)/2.0;
	double J = 0.0;
	int i, j;
	for(i=0;i<5;i++)
	{
		double JX = 0.0;
		double x = h1*qsi[i] + h2;
		// caso C = C(x) e D = D(x) substitui ymin e ymax pelas funcoes 
		double C = ymin, D = ymax;
		double k1 = (D-C)/2.0, k2 = (D+C)/2.0;
		for(j=0;j<5;j++)
		{
			double y = k1*qsi[j] + k2;
			double Q = f(x,y);
			JX += w[j]*Q;
		}
		J += w[i]*k1*JX;
	}
	J = h1*J;
	return J;
	
}

double interpola_cubica(double x0, double fx0, double x1, double fx1, double x2, double fx2, double x3, double fx3, double x)
{
	double gamma0 = fx0;
	double gamma1 = (fx1 - gamma0)/(x1 - x0);
	double gamma2 = (fx2 - gamma0 - gamma1*(x2 - x0))/((x2-x0)*(x2-x1));
	double gamma3 = (fx3 - gamma0 - gamma1*(x3 - x0) - gamma2*(x3 - x0)*(x3 - x1))/((x3 - x0)*(x3 - x2)*(x3 - x1));
	double interp = (gamma0 + gamma1*(x - x0) + gamma2*(x - x0)*(x - x1) + gamma3*(x - x0)*(x - x1)*(x - x2));

	return(interp);
}

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 );
}

/*
double integra_20pontos(double xmin, double xmax, double ymin, double ymax, Funcao *f)
{
	double qsi[20], w[20];
	// raizes qsi[j] e pesos w[j] da integracao de gauss 
	qsi[0] = -0.993129;
	qsi[1] = -0.963972;
	qsi[2] = -0.912234;
	qsi[3] = -0.839117;
	qsi[4] = -0.746332;
	qsi[5] = -0.636054;
	qsi[6] = -0.510867;
	qsi[7] = -0.373706;
	qsi[8] = -0.227786;
	qsi[9] = -0.0765265;
	qsi[10] = 0.0765265;
	qsi[11] = 0.227786;
	qsi[12] = 0.373706;
	qsi[13] = 0.510867;
	qsi[14] = 0.636054;
	qsi[15] = 0.746332;
	qsi[16] = 0.839117;
	qsi[17] = 0.912234;
	qsi[18] = 0.963972;
	qsi[19] = 0.993129;
	w[0] = 0.017614;
	w[1] = 0.0406014;
	w[2] = 0.062672;
	w[3] = 0.0832767;
	w[4] = 0.10193;
	w[5] = 0.118195;
	w[6] = 0.131689;
	w[7] = 0.142096;
	w[8] = 0.149173;
	w[9] = 0.152753;
	w[10] = 0.152753;
	w[11] = 0.149173;
	w[12] = 0.142096;
	w[13] = 0.131689;
	w[14] = 0.118195;
	w[15] = 0.10193;
	w[16] = 0.0832767;
	w[17] = 0.062672;
	w[18] = 0.0406014;
	w[19] = 0.017614;
	// comeca a ingtegracao..

	// dominio de integracao [A,B]x[C,D]
	double A = xmin, B = xmax;

	// variaveis auxiliares
	double h1 = (B-A)/2.0, h2 = (B+A)/2.0;
	double J = 0.0;
	int i, j;
	for(i=0;i<20;i++)
{
		double JX = 0.0;
		double x = h1*qsi[i] + h2;
		// caso C = C(x) e D = D(x) substitui ymin e ymax pelas funcoes 
		double C = ymin, D = ymax;
		double k1 = (D-C)/2.0, k2 = (D+C)/2.0;
		for(j=0;j<20;j++)
{
			double y = k1*qsi[j] + k2;
			double Q = f(x,y);
			JX += w[j]*Q;
}
		J += w[i]*k1*JX;
}
	J = h1*J;
	return J;
	
}
*/