/***************************************************************************
 *   Copyright (C) 2008 by Douglas Castro,,,   *
 *   douglas@douglas-laptop   *
 *                                                                         *
 ***************************************************************************/
/* Last edited on 2008-08-21 13:15:51 by stolfi */

#include "interp_integ.h"

/**
 * Dados limites do dominio Q=[xmin,xmax]x[ymin,ymax] calcula-se a integral de f em Q
 * @param xmin [in] abscissa inferior do dominio de integracao
 * @param xmax [in] abscissa superior do dominio de integracao
 * @param ymin [in] ordenda inferior do dominio de integracao
 * @param ymax [in] ordenada superior do dominio de integracao
 * @param (f) funcao a ser integrada
 */
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;
	
}


/**
 * Dados limites do dominio Q=[xmin,xmax]x[ymin,ymax] calcula-se a integral de f em Q
 * @param xmin [in] abscissa inferior do dominio de integracao
 * @param xmax [in] abscissa superior do dominio de integracao
 * @param ymin [in] ordenda inferior do dominio de integracao
 * @param ymax [in] ordenada superior do dominio de integracao
 * @param (f) funcao a ser integrada
 */
/*
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.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;
	
}
*/


/**
 * Funcao de interpolacao cubica. Recebe quatro valores {fx0,fx1,fx2,fx3}
 * de uma funcao {f} arbitrária em quatro pontos igualmente espacados {x0,x1,x2,x3},
 * e uma posicao {x}, e interpola {f(x)}.
 * @param x0 [in] abscissa 
 * @param fx0 [in] ordenada
 * @param x1 [in] abscissa 
 * @param fx1 [in] ordenada
 * @param x2 [in] abscissa 
 * @param fx2 [in] ordenada
 * @param x3 [in] abscissa 
 * @param fx3 [in] ordenada
 * @param x [in] abscissa entre x1 e x2, onde sera feita interpolacao
 */
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 )
{
  return ( x1 + (x0 - x2)/8.0 );
}