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

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

#include "multipole.hh"
#include "interpolate.hh"
#include "utils.hh"
#include "sphericalharmonic.hh"

static const int var_name_length = 30;
static const int max_vars = 10;

static const int max_spin_weights = 3;
static const int max_l_modes = 10;
static const int max_m_modes = 2 * max_l_modes + 1;

typedef struct
{
  int index;
  int imag_index;
  int spin_weight;
  char name[var_name_length];
}
variable_desc;

typedef struct
{
  int n_vars;
  variable_desc *vars;
}
variables_desc;


static void fill_variable(int idx, const char *optstring, void *callback_arg)
{
  assert(idx >= 0);
  assert(callback_arg != 0);

  variables_desc *vs = (variables_desc * ) callback_arg;

  assert(vs->n_vars < max_vars); // Too many variables in the variables list
  variable_desc *v = &vs->vars[vs->n_vars];

  v->index = idx;

  // Default values if there is no option string or if the options are
  // not present
  v->imag_index = -1;
  v->spin_weight = 0;
  strcpy(v->name, CCTK_VarName(v->index));

  if (optstring != 0)
  {
    int table = Util_TableCreateFromString(optstring);

    if (table >= 0)
    {
      const int buffer_length = 256;
      char imag_name[buffer_length];
      
      Util_TableGetInt(table, &v->spin_weight , "sw");
/////////////////////////////////////////////////////////////
CCTK_VInfo(CCTK_THORNSTRING,"spinweight %d", v->spin_weight);
/////////////////////////////////////////////////////////////
      if (Util_TableGetString(table, buffer_length, imag_name , "cmplx") >= 0)
      {
        v->imag_index = CCTK_VarIndex(imag_name);
      }
      Util_TableGetString(table, var_name_length, v->name , "name");
    }
  }
  vs->n_vars++;
}

static void parse_variables_string(const char *var_string, variable_desc v[max_vars], int *n_variables)
{
  variables_desc  vars;

  vars.n_vars = 0;
  vars.vars = v;

  const char *actual_var_string = strlen(var_string) == 0 ? 
    "WeylScal4::Psi4r{sw=-2 cmplx=\'WeylScal4::Psi4i\' name=\'psi4\'}" : 
    var_string;

  if (strlen(var_string) == 0)
  {
    static bool have_warned = false;
    if (!have_warned)
    {
      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "WARNING: No Multipole::variables string has been specified.  For compatibility, the current default is to use \"%s\" but in future, this default will change to \"\".  Please update your parameter files to include this new parameter.", actual_var_string);
      have_warned = true;
    }
  }
    
  int ierr = CCTK_TraverseString(actual_var_string, fill_variable, &vars, CCTK_GROUP_OR_VAR);
  assert(ierr > 0);

  *n_variables = vars.n_vars;
}

static void output_mode(CCTK_ARGUMENTS, const variable_desc *v, CCTK_REAL rad, 
                     CCTK_REAL real_lm, CCTK_REAL imag_lm, int l, int m)
{
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS

  char name_tmp[1000];

  if (CCTK_MyProc(cctkGH) == 0)
  {
    sprintf(name_tmp, "mp_%s_l%d_m%d_r%1.2f.asc", v->name, l, m, rad);
    Multipole_OutputComplexToFile(CCTK_PASS_CTOC, name_tmp, real_lm, imag_lm);
  }
}

static void output_1D(CCTK_ARGUMENTS, const variable_desc *v, CCTK_REAL rad, 
                      CCTK_REAL *th, CCTK_REAL *ph,
                      CCTK_REAL *real, CCTK_REAL *imag, int array_size)
{
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS

  char name_tmp[1000];

  if (CCTK_MyProc(cctkGH) == 0)
  {
    if (out_1d_every != 0 && (cctk_iteration) % out_1d_every == 0)
    {
      const char *real_name = CCTK_VarName(v->index);
      const char *imag_name = CCTK_VarName(v->imag_index);

      sprintf(name_tmp, "mp_%s_r%1.2f.th.asc", real_name, rad);
      Multipole_Output1D(CCTK_PASS_CTOC, name_tmp, array_size, th, ph, mp_theta, real);
      sprintf(name_tmp, "mp_%s_r%1.2f.th.asc", imag_name, rad);
      Multipole_Output1D(CCTK_PASS_CTOC, name_tmp, array_size, th, ph, mp_theta, imag);

      sprintf(name_tmp, "mp_%s_r%1.2f.ph.asc", real_name, rad);
      Multipole_Output1D(CCTK_PASS_CTOC, name_tmp, array_size, th, ph, mp_phi, real);
      sprintf(name_tmp, "mp_%s_r%1.2f.ph.asc", imag_name, rad);
      Multipole_Output1D(CCTK_PASS_CTOC, name_tmp, array_size, th, ph, mp_phi, imag);
    }
  }
}

bool int_in_array(int a, const int array[], int len)
{
  for (int i = 0; i < len; i++)
  {
    if (array[i] == a)
      return true;
  }
  return false;
}

int find_int_in_array(int a, const int array[], int len)
{
  for (int i = 0; i < len; i++)
  {
    if (array[i] == a)
      return i;
  }
  return -1;
}


static
void get_spin_weights(variable_desc vars[], int n_vars, int spin_weights[max_spin_weights], int *n_weights)
{
  int n_spin_weights = 0;

  for (int i = 0; i < n_vars; i++)
  {
    if (!int_in_array(vars[i].spin_weight, spin_weights, n_spin_weights))
    {
      assert(n_spin_weights < max_spin_weights);
      spin_weights[n_spin_weights] = vars[i].spin_weight;
      n_spin_weights++;
    }
  }
  *n_weights = n_spin_weights;
}


static
void setup_harmonics(const int spin_weights[max_spin_weights], int n_spin_weights, int lmin, int lmax, 
                     int m_mode, bool allmodes, CCTK_REAL th[], CCTK_REAL ph[], int array_size, 
                     CCTK_REAL *reY[max_spin_weights][max_l_modes][max_m_modes],
                     CCTK_REAL *imY[max_spin_weights][max_l_modes][max_m_modes])
{
  for (int si = 0; si < max_spin_weights; si++)
  {
    for (int l = 0; l < max_l_modes; l++)
    {
      for (int mi = 0; mi < max_m_modes; mi++)
      {
        reY[si][l][mi] = 0;
        imY[si][l][mi] = 0;
      }
    }
  }
  for (int si = 0; si < n_spin_weights; si++)
  {
    int sw = spin_weights[si];

    for (int l=lmin; l <= lmax; l++)
    {
      int mmin = (allmodes ? (m_mode > l ? -l : -m_mode) : m_mode);
      int mmax = (allmodes ? (m_mode > l ? l : m_mode) : m_mode);

      for (int m=mmin; m <= mmax; m++)
      {
        reY[si][l-lmin][m-mmin] = new CCTK_REAL[array_size];
        imY[si][l-lmin][m-mmin] = new CCTK_REAL[array_size];

        Multipole_HarmonicSetup(sw, l, m, array_size, th, ph, reY[si][l-lmin][m-mmin], imY[si][l-lmin][m-mmin]);
      }
    }
  } 
}

extern "C" void Multipole_Calc(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS

  static CCTK_REAL *xs, *ys, *zs;
  static CCTK_REAL *xhat, *yhat, *zhat;
  static CCTK_REAL *th, *ph;
  static CCTK_REAL *real = 0, *imag = 0;
  static CCTK_REAL *reY[max_spin_weights][max_l_modes][max_m_modes];
  static CCTK_REAL *imY[max_spin_weights][max_l_modes][max_m_modes];
  static variable_desc   vars[max_vars];
  static int n_variables = 0;
  static int spin_weights[max_spin_weights];
  static int n_spin_weights = 0;

  static bool initialized = false;

  int i, array_size=(ntheta+1)*(nphi+1);
  CCTK_REAL real_lm = 0.0, imag_lm = 0.0;
  int lmin, lmax, mmin, mmax, allmodes=0;

  if (cctk_iteration % out_every != 0)
    return;

  // Determine modes to compute
  if (CCTK_EQUALS(mode_type, "all modes"))
  {
    allmodes=1;
  }
  //if mode_type is set to "all modes", all modes with l <= l_mode, 
  //abs(m) <= min{l_mode, m_mode} will be output.  Otherwise, only
  //the l_mode, m_mode mode will be output.
  lmin = (allmodes ? l_min : l_mode);
  lmax = l_mode;

  assert(lmax + 1 <= max_l_modes);

  if (!initialized)
  {
    printf("Multipole: Allocating memory in Multipole_Calc\n");
    real = new CCTK_REAL[array_size];
    imag = new CCTK_REAL[array_size];
    th = new CCTK_REAL[array_size];  
    ph = new CCTK_REAL[array_size];
    xs = new CCTK_REAL[array_size];
    ys = new CCTK_REAL[array_size];
    zs = new CCTK_REAL[array_size];
    xhat = new CCTK_REAL[array_size];
    yhat = new CCTK_REAL[array_size];
    zhat = new CCTK_REAL[array_size];
  
    printf("Multipole: parsing variables string.\n");
    parse_variables_string(variables, vars, &n_variables);
    get_spin_weights(vars, n_variables, spin_weights, &n_spin_weights);
    Multipole_CoordSetup(ntheta, nphi, xhat, yhat, zhat, th, ph);
    setup_harmonics(spin_weights, n_spin_weights, lmin, lmax, m_mode, allmodes, th, ph, array_size, reY, imY);
    initialized = true;
  }

  for (int v = 0; v < n_variables; v++)
  {
    //assert(vars[v].spin_weight == -2);

    int si = find_int_in_array(vars[v].spin_weight, spin_weights, n_spin_weights);
    assert(si != -1);

    for(i=0; i<nradii; i++)
    {
      // Compute x^i = r * \hat x^i
      Multipole_ScaleCartesian(ntheta, nphi, radius[i], xhat, yhat, zhat, xs, ys, zs);
      
      // Interpolate Psi4r and Psi4i
      Multipole_Interp(CCTK_PASS_CTOC, xs, ys, zs, vars[v].index, vars[v].imag_index, 
                       real, imag);
      for (int l=lmin; l <= lmax; l++)
      {
        mmin = (allmodes ? (m_mode > l ? -l : -m_mode) : m_mode) ;
        mmax = (allmodes ? (m_mode > l ? l : m_mode) : m_mode);
        
        for (int m=mmin; m <= mmax; m++)
        {
          // Integrate sYlm (real + i imag) over the sphere at radius r
          Multipole_Integrate(array_size, ntheta, reY[si][l-lmin][m-mmin], imY[si][l-lmin][m-mmin], real, imag, th, ph,
                              &real_lm, &imag_lm);
    
          output_mode(CCTK_PASS_CTOC, &vars[v], radius[i], real_lm, imag_lm, l, m);
        }//loop over m
      }//loop over l
      output_1D(CCTK_PASS_CTOC, &vars[v], radius[i], th, ph, real, imag, array_size);
    }//loop over radii
  }//loop over variables
}