#ifndef QUADRATURE_H
#define QUADRATURE_H

typedef struct 
{
      int size;
      double* points;
      double* weights;
} Quadrature;

Quadrature SIMPSONS_RULE;
Quadrature FIVE_POINT_GAUSS;

void Create_Quadratures()
{
   SIMPSONS_RULE.size = 3;
   SIMPSONS_RULE.points = new double[3];
   SIMPSONS_RULE.points[0]=0.0; SIMPSONS_RULE.points[1]=0.5; SIMPSONS_RULE.points[2]=1.0;
   SIMPSONS_RULE.weights = new double[3];
   SIMPSONS_RULE.weights[0]=1.0/6; SIMPSONS_RULE.weights[1]=4.0/6; SIMPSONS_RULE.weights[2]=1.0/6;

   //NOTE: Unlike the Simpson's rule above (which is given exact in machine precision),
   //  the points for the 5 point Gaussian rule are given manually and are exact up to the
   //  8-9th decimal point. Thus it'll evaluate the integrals with an error ~10^(-9) even
   //  for a low degree polynomials (for which the Simpson's rule will give an error ~10^(-16)).
   FIVE_POINT_GAUSS.size = 5;
   FIVE_POINT_GAUSS.points = new double[5];
   FIVE_POINT_GAUSS.points[0] = 0.046910075;  FIVE_POINT_GAUSS.points[1] = 0.230765345;
   FIVE_POINT_GAUSS.points[2] = 0.5;
   FIVE_POINT_GAUSS.points[3] = 0.769234655;  FIVE_POINT_GAUSS.points[4] = 0.953089925; 
   FIVE_POINT_GAUSS.weights = new double[5];
   FIVE_POINT_GAUSS.weights[0] = 0.118463445;  FIVE_POINT_GAUSS.weights[1] = 0.239314335;
   FIVE_POINT_GAUSS.weights[2] = 0.284444445;
   FIVE_POINT_GAUSS.weights[3] = 0.239314335;  FIVE_POINT_GAUSS.weights[4] = 0.118463445; 
}

void Destroy_Quadratures()
{
   delete SIMPSONS_RULE.points;
   delete SIMPSONS_RULE.weights;

   delete FIVE_POINT_GAUSS.points;
   delete FIVE_POINT_GAUSS.weights;
}

#endif