/* See SOIntegral.h */
/* Last edited on 2003-05-07 16:33:01 by ra014852 */

#include <SOIntegral.h>
#include <SOBasic.h>
#include <js.h>
#include <math.h>
#include <stdlib.h>

int SOIntegral_GaussOrder = 10;
  /* Number of knots to be used by {SOIntegral_Gauss} below,
    along each coordinate axis. */

double SOIntegral_MaxSize = 0.01;
  /* Maximum size of an integral area size, considering that 
     in general integrals are used to calculate dot products
     over the root cell, that measures 1 along each side. */

/* Knots and weights for 2-point Gaussian quadrature: */
static double qpt_2[2] = 
  { -0.577350269189626, 
     0.577350269189626
  };
static double qwt_2[1] = 
  {  1.000000000000000 
  }; 

/* Knots and weights for 4-point Gaussian quadrature: */
static double qpt_4[4] = 
  { -0.861136311594053, 
     0.861136311594053, 
    -0.339981043584856, 
     0.339981043584856
  };
static double qwt_4[2] =
  { 0.347854845137454, 
    0.652145154862546
  };

/* Knots and weights for 6-point Gaussian quadrature: */
static double qpt_6[6] =  
  { -0.932469514203152, 
     0.932469514203152, 
    -0.661209386466265,
     0.661209386466265, 
    -0.238619186083197, 
     0.238619186083197
  };
static double qwt_6[3] = 
  { 0.171324492379170, 
    0.360761573048139, 
    0.467913934572691 
  };

/* Knots and weights for 8-point Gaussian quadrature: */
static double qpt_8[8] =
  { -0.960289856497536,
     0.960289856497536,
    -0.796666477413627,
     0.796666477413627,
    -0.525532409916329,
     0.525532409916329,
    -0.183434642495650,
     0.183434642495650
  };
static double qwt_8[4] = 
  { 0.101228536290376,
    0.222381034453374,
    0.313706645877887,
    0.362683783378362
  };

/* Knots and weights for 10-point Gaussian quadrature: */
static double qpt_10[10] =
  { -0.973906528517172,
     0.973906528517172,
    -0.865063366688985,    
     0.865063366688985, 
    -0.679409568299024,
     0.679409568299024,
    -0.433395394129247,
     0.433395394129247,
    -0.148874338981631,
     0.148874338981631
  };
static double qwt_10[5] = 
  { 0.066671344308688,
    0.149451349150581,
    0.219086362515982, 
    0.269266719309996,
    0.295524224714753
  };

/* Knots and weights for 12-point Gaussian quadrature: */
static double qpt_12[12] =
  { -0.981560634246719,
     0.981560634246719,
    -0.904117256370475,
     0.904117256370475,
    -0.769902674194305,
     0.769902674194305,
    -0.587317954286617,
     0.587317954286617,
    -0.367831498998180,
     0.367831498998180,
    -0.125233408511469,
     0.125233408511469
  };
static double qwt_12[6] = 
  { 0.047175336386512,
    0.106939325995318,
    0.160078328543346,
    0.203167426723066,
    0.233492536538355,
    0.249147045813403
  };


/* Pointers to knot and weight tables: */
static double *qpoints[6] = 
  { (double *)qpt_2,
    (double *)qpt_4,
    (double *)qpt_6,
    (double *)qpt_8,
    (double *)qpt_10,
    (double *)qpt_12
  };
static double *qweight[6] =
  { (double *)qwt_2,
    (double *)qwt_4,
    (double *)qwt_6,
    (double *)qwt_8,
    (double *)qwt_10,
    (double *)qwt_12
  };


void SOIntegral_Gauss
  ( SOIntegral_Func f, 
    SOGrid_Dim d, 
    SOGrid_Dim fd, 
    double *sum, 
    double *corr
  )
  { int ktable = (SOIntegral_GaussOrder/2)-1;   /* Which tables to use */
    double *pt = qpoints[ktable]; /* Knot positions in [-1_+1]. */
    double *wt = qweight[ktable];
    double volfactor = (double)1.0/((double)(1<<d)); /* {2^{-d}}: */
    static int ix[MAX_PDIM];    /* Counters for the grid. */
    static double x[MAX_PDIM];  /* Coordinates of sample points, in [0_1]. */
    static double fx[MAX_FDIM]; /* Function value at {x}. */
    int i, j;  
    
    /* Initialize counter variables: */
    for (i = 0; i < d; i++) { ix[i] = 0;  x[i] = (double)(1 + pt[0])/2; }
    
    /* Simulate {d} nested loops, with {ix[i]} varying from 0 to {SOIntegral_GaussOrder-1}: */
    do 
      { 
        /* Compute weight {w} of sample point, including the box measure: */
        double w = volfactor;
	for (i = 0; i < d; i++) { w = w * wt[ix[i]/2]; }
       
        /* Evaluate the function at the sample point: */
        f(&(x[0]), &(fx[0]));
        
        /* Accumulate {w*f(x)} on the integral: */
         for (j = 0; j < fd; j++)
          { double term = w * fx[j];
            /* Kahan's summation formula: */
            double sumj = sum[j];
            double tcorr = term - corr[j];
            double newSum = sumj + tcorr;
            corr[j] = (newSum - sumj) - tcorr;
            sum[j] = newSum;
          }
        
        /* Increment the loop counters and update sample point: */
        i = d;
        do
          { i--;
            ix[i] = (ix[i] + 1) % SOIntegral_GaussOrder;
            x[i] = (double)(1 + pt[ix[i]])/2; 
          }
        while ((i > 0) && (ix[i] == 0));
        /* After the last point, {i = ix[0] = 0}. */
      }
    while ((i > 0) || (ix[i] != 0));  

  }


void SOIntegral_Gauss_Limits
  ( SOIntegral_Func f, 
    SOGrid_Dim d, 
    SOGrid_Dim fd, 
    double *sum, 
    double *corr,
    double *min,    /* Lower coordinates of integration limits. */
    double *max,    /* Higher coordinates of integration limits. */
    bool verbose
  )
  { int ktable = (SOIntegral_GaussOrder/2)-1;   /* Which tables to use */
    double *pt = qpoints[ktable]; /* Knot positions in [-1_+1]. */
    double *wt = qweight[ktable];
    double volfactor = 1;
    static int ix[MAX_PDIM];    /* Counters for the grid. */
    static double x[MAX_PDIM];  /* Coordinates of sample points, in [0_1]. */
    static double fx[MAX_FDIM]; /* Function value at {x}. */
    int i, j;  
    
    /* Initialize counter variables: */
    for (i = 0; i < d; i++) 
      { 
        volfactor *= (double)(max[i]-min[i]) / 2;

	//        volfactor = volfactor/pow(2.0, (double)i/d);   // Canonical geometry.

        ix[i] = 0;  
        x[i] = (double)((max[i]+min[i])+(double)(max[i]-min[i])*pt[0])/2; 
      }
    
    /* Simulate {d} nested loops, with {ix[i]} varying from 0 to {SOIntegral_GaussOrder-1}: */
    do 
      { 
        /* Compute weight {w} of sample point, including the box measure: */
        double w = volfactor;
	for (i = 0; i < d; i++) { w = (double)w * wt[ix[i]/2]; }
       
        /* Evaluate the function at the sample point: */
        f(&(x[0]), &(fx[0]));
    
        /* Accumulate {w*f(x)} on the integral: */
         for (j = 0; j < fd; j++)
          { double term = (double)w * fx[j];
            /* Kahan's summation formula: */
            double sumj = sum[j];
            double tcorr = term - corr[j];
            double newSum = sumj + tcorr;
            corr[j] = (newSum - sumj) - tcorr;
            sum[j] = newSum;
          }
        
        /* Increment the loop counters and update sample point: */
        i = d;
        do
          { i--;
            ix[i] = (ix[i] + 1) % SOIntegral_GaussOrder;
            x[i] = (double)((max[i]+min[i])+(double)(max[i]-min[i])*pt[ix[i]])/2;
          }
        while ((i > 0) && (ix[i] == 0));
        /* After the last point, {i = ix[0] = 0}. */
      }
    while ((i > 0) || (ix[i] != 0));  

  }


void SOIntegral_OnRootCell
  ( SOIntegral_Func f, 
    SOGrid_Dim d,       /* Dimension {m} of domain. */
    SOGrid_Dim fd,      /* Dimension {n} of the function's range. */
    SOGrid_Tree *tree,  /* Finite cell grid to use. */
    double *sum, 
    double *corr
  )
  {
    /* If single cell, use SOIntegral_Gauss, else recurse on
       each sub-cell and add results. */
    //    if (tree == NULL)
    //      { SOIntegral_Gauss(f, d, fd, sum, corr); }
    //    else
    { assert(FALSE, "not implemented yet"); }
  }