#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include "quadrature.h"

// http://en.wikipedia.org/wiki/Gaussian_quadrature

/**
@brief Take index and degree to integration point in Interval or [-1,1]
*/
float ixs(int ind, int N, Interval x){return x.lower+(x.upper-x.lower)*ind/(N+1);	}
float ix(int ind, int N){return ixs(ind,N,(Interval){-1.,1.});}

float gauss_to_interval(float a, Interval x)
{
	return .5*(x.lower+x.upper) + .5*(x.upper-x.lower)*a;
}

// x - integration points, w - integration weights
// for interval [-1,1]
// generated with python
static float x1[1] = {0.};
static float w1[1] = {2.};
static float x2[2] = {0.57735026919,-0.57735026919};
static float w2[2] = {1.,1.};
static float x3[3] = {0.,0.774596669241,-0.774596669241};
static float w3[3] = {0.888888888889,0.555555555556,0.555555555556};
static float x4[4] = {0.339981043585,-0.339981043585,0.861136311594,-0.861136311594};
static float w4[4] = {0.652145154863,0.652145154863,0.347854845137,0.347854845137};
static float x5[5] = {0.,0.538469310106,-0.538469310106,0.906179845939,-0.906179845939};
static float w5[5] = {0.501960784314,0.478628670499,0.478628670499,0.236926885056,0.236926885056};

/*
Return the gaussian quadrature absissa given by the ind point in the N point rule
*/
float gauss_pt(int ind, const int N)
{
	assert(ind<N);
	switch (N){
		case 1:
			return x1[ind];
		case 2:
			return x2[ind];
		case 3:
			return x3[ind];
		case 4:
			return x4[ind];
		case 5:
			return x5[ind];
	}
	return 0.;
}
float gauss_wt(int ind, const int N)
{
	assert(ind<N);
	switch (N){
		case 1:
			return w1[ind];
		case 2:
			return w2[ind];
		case 3:
			return w3[ind];
		case 4:
			return w4[ind];
		case 5:
			return w5[ind];
	}
	return 0.;
}