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

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

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

using namespace std;

////////////////////////////////////////////////////////////////////////
// File manipulation
////////////////////////////////////////////////////////////////////////

FILE *Multipole_OpenOutputFile(CCTK_ARGUMENTS, const string &name)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  
  bool first_time = cctk_iteration == 0;
  const char *mode = first_time ? "w" : "a";
  CCTK_STRING *out_dir  = (CCTK_STRING *) CCTK_ParameterGet("out_dir", "IOUtil", NULL);

  if (*out_dir == 0)
  {
    CCTK_WARN(1,"Parameter IOUtil::out_dir not found");
    return 0;
  }

  string output_name(string(*out_dir) + "/" + name);

  FILE *fp = fopen(output_name.c_str(), mode);

  if (fp == 0)
  {
    CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, (string("Could not open output file ") + output_name).c_str());
  }

  return fp;
}

////////////////////////////////////////////////////////////////////////
// Unused
////////////////////////////////////////////////////////////////////////

void Multipole_OutputArray(CCTK_ARGUMENTS, FILE *f, int array_size, 
                           CCTK_REAL const th[], CCTK_REAL const ph[], 
                           CCTK_REAL const xs[], CCTK_REAL const ys[], CCTK_REAL const zs[],
                           CCTK_REAL const data[])
{
  DECLARE_CCTK_PARAMETERS;
  DECLARE_CCTK_ARGUMENTS;

  CCTK_REAL last_ph = ph[0];

  for (int i = 0; i < array_size; i++)
  {
    if (ph[i] != last_ph) // Separate blocks for gnuplot
      fprintf(f, "\n");
    fprintf(f, "%f %f %f %f %f %f %.19g\n", cctk_time, th[i], ph[i], xs[i], ys[i], zs[i], data[i]);
    last_ph = ph[i];
  }
}


void Multipole_OutputArrayToFile(CCTK_ARGUMENTS, const string &name, int array_size,
                                 CCTK_REAL const th[], CCTK_REAL const ph[],
                                 CCTK_REAL const xs[], CCTK_REAL const ys[], CCTK_REAL const zs[],
                                 CCTK_REAL const data[])
{
  DECLARE_CCTK_ARGUMENTS;

  if (FILE *fp = Multipole_OpenOutputFile(CCTK_PASS_CTOC, name))
  {
    Multipole_OutputArray(CCTK_PASS_CTOC, fp, array_size, th, ph, xs, ys, zs, data);
    fclose(fp);
  }
}

////////////////////////////////////////////////////////////////////////
// Misc
////////////////////////////////////////////////////////////////////////

void Multipole_Output1D(CCTK_ARGUMENTS, const string &name, int array_size,
                        CCTK_REAL const th[], CCTK_REAL const ph[], mp_coord coord,
                        CCTK_REAL const data[])
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;

  if (FILE *f = Multipole_OpenOutputFile(CCTK_PASS_CTOC, name))
  {
    fprintf(f, "\"Time = %.19g\n", cctk_time);

    if (coord == mp_theta)
    {
      for (int i = 0; i <= ntheta; i++)
      {
        int idx = Multipole_Index(i, 0, ntheta);
        fprintf(f, "%f %.19g\n", th[idx], data[idx]);
      }
    }
    else if (coord == mp_phi)
    {
      for (int i = 0; i <= nphi; i++)
      {
        int idx = Multipole_Index(ntheta / 4, i, ntheta);
        fprintf(f, "%f %.19g\n", ph[idx], data[idx]);
      }
    }
    fprintf(f, "\n\n");
    fclose(f);
  }
}

////////////////////////////////////////////////////////////////////////
// Complex IO
////////////////////////////////////////////////////////////////////////

void Multipole_OutputComplex(CCTK_ARGUMENTS, FILE *fp, CCTK_REAL redata, CCTK_REAL imdata)
{
  DECLARE_CCTK_PARAMETERS;
  DECLARE_CCTK_ARGUMENTS;
  fprintf(fp, "%f %.19g %.19g\n", cctk_time, redata, imdata);
}

void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, const string &name, CCTK_REAL redata, CCTK_REAL imdata)
{
  DECLARE_CCTK_ARGUMENTS;

  if (FILE *fp = Multipole_OpenOutputFile(CCTK_PASS_CTOC, name))
  {
    Multipole_OutputComplex(CCTK_PASS_CTOC, fp, redata, imdata);
    fclose(fp);
  }
}


////////////////////////////////////////////////////////////////////////
// Coordinates
////////////////////////////////////////////////////////////////////////

void Multipole_CoordSetup(int ntheta, int nphi, 
                          CCTK_REAL xhat[], CCTK_REAL yhat[], 
                          CCTK_REAL zhat[], CCTK_REAL th[], 
                          CCTK_REAL ph[])
{
  const CCTK_REAL PI = acos(-1.0);

  for (int it = 0; it <= ntheta; it++)
  {
    for (int ip = 0; ip <= nphi; ip++)
    {
      const int i = Multipole_Index(it, ip, ntheta);

      th[i] = it * PI / (ntheta);
      ph[i] = ip * 2 * PI / nphi;
      xhat[i] = cos(ph[i])*sin(th[i]);
      yhat[i] = sin(ph[i])*sin(th[i]);
      zhat[i] = cos(th[i]);
    }
  }
}

void Multipole_ScaleCartesian(int ntheta, int nphi, CCTK_REAL r,
                              CCTK_REAL const xhat[], CCTK_REAL const yhat[], CCTK_REAL const zhat[],
                              CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[])
{
  for (int it = 0; it <= ntheta; it++)
  {
    for (int ip = 0; ip <= nphi; ip++)
    {
      const int i = Multipole_Index(it, ip, ntheta);

      x[i] = r * xhat[i];
      y[i] = r * yhat[i];
      z[i] = r * zhat[i];
    }
  }
}

////////////////////////////////////////////////////////////////////////
// Integration
////////////////////////////////////////////////////////////////////////

void 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)
{
  DECLARE_CCTK_PARAMETERS

  int il = Multipole_Index(0,0,ntheta);
  int iu = Multipole_Index(1,0,ntheta);
  CCTK_REAL dth = th[iu] - th[il];
  iu = Multipole_Index(0,1,ntheta);
  CCTK_REAL dph = ph[iu] - ph[il];

  if (CCTK_Equals(integration_method, "midpoint"))
  {
    CCTK_REAL tempr = 0.0;
    CCTK_REAL tempi = 0.0;
      
    for (int i=0; i < array_size; i++) {
      // the below calculations take the integral of conj(array1)*array2
      tempr += ( array1r[i]*array2r[i] + array1i[i]*array2i[i] )
        *sin(th[i])*dth*dph;
      tempi += ( array1r[i]*array2i[i] - array1i[i]*array2r[i] )
        *sin(th[i])*dth*dph;
      
      *outre = tempr;
      *outim = tempi;
    } 
  }
  else if (CCTK_Equals(integration_method, "simpson"))
  {
    static CCTK_REAL *fr = 0;
    static CCTK_REAL *fi = 0;
    static bool allocated_memory = false;

    // Construct an array for the real integrand
    if (!allocated_memory)
    {
      fr = new CCTK_REAL[array_size];
      fi = new CCTK_REAL[array_size];
      allocated_memory = true;
    }
    
    for (int i = 0; i < array_size; i++)
    {
      fr[i] = (array1r[i] * array2r[i] + 
               array1i[i] * array2i[i] ) * sin(th[i]);
      fi[i] = (array1r[i] * array2i[i] - 
               array1i[i] * array2r[i] ) * sin(th[i]);
    }
    
    *outre = Simpson2DIntegral(fr, ntheta, nphi, dth, dph);
    *outim = Simpson2DIntegral(fi, ntheta, nphi, dth, dph);
  }
}