#include <stdio.h>
#include <string.h>
#include <math.h>
#include <assert.h>

#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"

#include "utils.hh"
#include "integrate.hh"

static CCTK_REAL test_simpson_integral(int n)
{
  const int nx = n;
  const int ny = n;
  const int array_size=(nx+1)*(ny+1);

  CCTK_REAL *f = new CCTK_REAL[array_size];

  const CCTK_REAL dx = 1./nx;
  const CCTK_REAL dy = 1./ny;

  for (int ix = 0; ix <= nx; ix++)
  {
    for (int iy = 0; iy <= ny; iy++)
    {
      const int i = Multipole_Index(ix, iy, nx);

      const CCTK_REAL x = ix*dx;
      const CCTK_REAL y = iy*dy;
      const CCTK_REAL PI = acos(-1.0);

      f[i] = x*pow(y,2)*pow(cos(2*PI*y),2)*pow(sin(2*PI*x),2);
    }
  }

  const CCTK_REAL result = Simpson2DIntegral(f, nx, ny, dx, dy);
  delete [] f;
  return result;
}


void Multipole_TestSimpson(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;

  const int n1 = 100;
  const int n2 = 200;
  const CCTK_REAL PI = acos(-1.0);
  const CCTK_REAL result1 = test_simpson_integral(100);
  const CCTK_REAL result2 = test_simpson_integral(200);
  const CCTK_REAL exact = 1./24 + 1./(64 * pow(PI,2));
  const CCTK_REAL error1 = fabs(result1 - exact);
  const CCTK_REAL error2 = fabs(result2 - exact);
  *test_simpson_convergence_order = log10(error1/error2) / log10((CCTK_REAL) n2/n1);
}

// void Multipole_TestIntegrate(CCTK_ARGUMENTS)
// {
//   const int n = 100;

//   const int nth = n;
//   const int nph = n;
//   const int array_size=(nth+1)*(nph+1);

//   CCTK_REAL *f = new CCTK_REAL[array_size];

//   const CCTK_REAL dth = 1/nx;
//   const CCTK_REAL dph = 1/ny;

//   for (int ix = 0; ix <= nx; ix++)
//   {
//     for (int iy = 0; iy <= ny; iy++)
//     {
//       const int i = Multipole_Index(ix, iy, nx);

//       const CCTK_REAL x = ix*dx;
//       const CCTK_REAL y = iy*dy;

//       f[i] = sin(2*PI*x)*cos(2*PI*y);
//     }
//   }




//   CCTK_REAL result = Multipole_Integrate(int array_size, int nthetap,
//                                          CCTK_REAL const array1r[], CCTK_REAL const array1i[],
//                                          CCTK_REAL const array2r[], CCTK_REAL const array2i[],
//                                          CCTK_REAL const th[], CCTK_REAL const ph[], 
//                                          CCTK_REAL *outre, CCTK_REAL *outim)



//   printf("Integration result: %.19g\n", result);
// }