#include <assert.h>
#include <stdlib.h>

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

#include "util_Table.h"



void
InterpToArray (CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  
  int interpolator;
  int options_table;
  int coord_handle;
  
  CCTK_REAL * restrict coordsx;
  CCTK_REAL * restrict coordsy;
  CCTK_REAL * restrict coordsz;
  CCTK_POINTER_TO_CONST coords[3];
  CCTK_INT * restrict inputs;
  CCTK_INT * restrict output_types;
  CCTK_POINTER * restrict outputs;
  CCTK_INT * restrict operation_codes;
  CCTK_INT * restrict time_deriv_order;
  
  int group;
  cGroupDynamicData dyndata;
  
  int nvars;
  int npoints;
  
  int n;
  int i, j, k;
  int ierr;
  
  
  
  interpolator = CCTK_InterpHandle (interpolator_name);
  assert (interpolator >= 0);
  
  options_table = Util_TableCreateFromString (interpolator_options);
  assert (options_table >= 0);
  
  coord_handle = CCTK_CoordSystemHandle (interpolator_coordinates);
  assert (coord_handle >= 0);
  
  
  
  /* Scalars */
  {
    nvars = nscalars;
    if (nvars > 0) {
      npoints = 1;
      
      coordsx = malloc (npoints * sizeof * coordsx);
      assert (npoints==0 || coordsx);
      coordsy = malloc (npoints * sizeof * coordsy);
      assert (npoints==0 || coordsy);
      coordsz = malloc (npoints * sizeof * coordsz);
      assert (npoints==0 || coordsz);
      coords[0] = coordsx;
      coords[1] = coordsy;
      coords[2] = coordsz;
      
      n = 0;
      assert (n <= npoints);
      coordsx[n] = scalar_x0;
      coordsy[n] = scalar_y0;
      coordsz[n] = scalar_z0;
      ++n;
      assert (n == npoints);
      
      inputs = malloc (nvars * sizeof * inputs);
      assert (inputs);
      
      for (n=0; n<nvars; ++n) {
        inputs[n] = CCTK_VarIndex (scalar_vars[n]);
        if (inputs[n] < 0) {
          inputs[n] = -1;
        }
      }
      
      output_types = malloc (nvars * sizeof * output_types);
      assert (output_types);
      
      for (n=0; n<nvars; ++n) {
        output_types[n] = CCTK_VARIABLE_REAL;
      }
      
      outputs = malloc (nvars * sizeof * outputs);
      assert (outputs);
      
      for (n=0; n<nvars; ++n) {
        outputs[n] = &scalars[npoints * n];
      }
      
      ierr = CCTK_InterpGridArrays
        (cctkGH, 3, interpolator, options_table, coord_handle,
         npoints, CCTK_VARIABLE_REAL, coords,
         nvars, inputs,
         nvars, output_types, outputs);
      assert (! ierr);
      
      free (coordsx);
      free (coordsy);
      free (coordsz);
      free (inputs);
      free (output_types);
      free (outputs);
    }
  }
  
  
  
  /* 1D Arrays */
  {
    nvars = narrays1d;
    if (nvars > 0) {
      npoints = array1d_npoints_i;
      
      coordsx = malloc (npoints * sizeof * coordsx);
      assert (npoints==0 || coordsx);
      coordsy = malloc (npoints * sizeof * coordsy);
      assert (npoints==0 || coordsy);
      coordsz = malloc (npoints * sizeof * coordsz);
      assert (npoints==0 || coordsz);
      coords[0] = coordsx;
      coords[1] = coordsy;
      coords[2] = coordsz;
      
      n = 0;
      for (i=0; i<array1d_npoints_i; ++i) {
        assert (n <= npoints);
        coordsx[n] = array1d_x0 + i * array1d_dx_i;
        coordsy[n] = array1d_y0 + i * array1d_dy_i;
        coordsz[n] = array1d_z0 + i * array1d_dz_i;
        ++n;
      }
      assert (n == npoints);
      
      inputs = malloc (nvars * sizeof * inputs);
      assert (inputs);
      
      for (n=0; n<nvars; ++n) {
        inputs[n] = CCTK_VarIndex (array1d_vars[n]);
        assert (inputs[n] >= 0);
        if (inputs[n] < 0) {
          inputs[n] = -1;
        }
      }
      
      output_types = malloc (nvars * sizeof * output_types);
      assert (output_types);
      
      for (n=0; n<nvars; ++n) {
        output_types[n] = CCTK_VARIABLE_REAL;
      }
      
      outputs = malloc (nvars * sizeof * outputs);
      assert (outputs);
      
      for (n=0; n<nvars; ++n) {
        outputs[n] = &arrays1d[npoints * n];
      }
      
      operation_codes = malloc (nvars * sizeof * operation_codes);
      assert (operation_codes);
      
      for (n=0; n<nvars; ++n) {
        operation_codes[n] = array1d_spacederivs[n];
      }
      
      ierr = Util_TableSetIntArray
        (options_table, nvars, operation_codes, "operation_codes");
      assert (! ierr);
      
      time_deriv_order = malloc (nvars * sizeof * time_deriv_order);
      assert (time_deriv_order);
      
      for (n=0; n<nvars; ++n) {
        time_deriv_order[n] = array1d_timederivs[n];
      }
      
      ierr = Util_TableSetIntArray
        (options_table, nvars, time_deriv_order, "time_deriv_order");
      assert (! ierr);
      
      ierr = CCTK_InterpGridArrays
        (cctkGH, 3, interpolator, options_table, coord_handle,
         npoints, CCTK_VARIABLE_REAL, coords,
         nvars, inputs,
         nvars, output_types, outputs);
      assert (! ierr);
      
      ierr = Util_TableDeleteKey (options_table, "operation_codes");
      assert (! ierr);
      
      ierr = Util_TableDeleteKey (options_table, "time_deriv_order");
      assert (! ierr);
      
      free (coordsx);
      free (coordsy);
      free (coordsz);
      free (inputs);
      free (output_types);
      free (outputs);
      free (operation_codes);
      free (time_deriv_order);
    }
  }
  
  
  
  /* 2D Arrays */
  {
    nvars = narrays2d;
    if (nvars > 0) {
      npoints = array2d_npoints_i * array2d_npoints_j;
      
      coordsx = malloc (npoints * sizeof * coordsx);
      assert (npoints==0 || coordsx);
      coordsy = malloc (npoints * sizeof * coordsy);
      assert (npoints==0 || coordsy);
      coordsz = malloc (npoints * sizeof * coordsz);
      assert (npoints==0 || coordsz);
      coords[0] = coordsx;
      coords[1] = coordsy;
      coords[2] = coordsz;
      
      n = 0;
      for (j=0; j<array2d_npoints_j; ++j) {
        for (i=0; i<array2d_npoints_i; ++i) {
          assert (n <= npoints);
          coordsx[n] = array2d_x0 + i * array2d_dx_i + j * array2d_dx_j;
          coordsy[n] = array2d_y0 + i * array2d_dy_i + j * array2d_dy_j;
          coordsz[n] = array2d_z0 + i * array2d_dz_i + j * array2d_dz_j;
          ++n;
        }
      }
      assert (n == npoints);
      
      inputs = malloc (nvars * sizeof * inputs);
      assert (inputs);
      
      for (n=0; n<nvars; ++n) {
        inputs[n] = CCTK_VarIndex (array2d_vars[n]);
        if (inputs[n] < 0) {
          inputs[n] = -1;
        }
      }
      
      output_types = malloc (nvars * sizeof * output_types);
      assert (output_types);
      
      for (n=0; n<nvars; ++n) {
        output_types[n] = CCTK_VARIABLE_REAL;
      }
      
      outputs = malloc (nvars * sizeof * outputs);
      assert (outputs);
      
      for (n=0; n<nvars; ++n) {
        outputs[n] = &arrays2d[npoints * n];
      }
      
      ierr = CCTK_InterpGridArrays
        (cctkGH, 3, interpolator, options_table, coord_handle,
         npoints, CCTK_VARIABLE_REAL, coords,
         nvars, inputs,
         nvars, output_types, outputs);
      assert (! ierr);
      
      free (coordsx);
      free (coordsy);
      free (coordsz);
      free (inputs);
      free (output_types);
      free (outputs);
    }
  }
  
  
  
  /* 3D Arrays */
  {
    nvars = narrays3d;
    if (nvars > 0) {
      npoints = array3d_npoints_i * array3d_npoints_j * array3d_npoints_k;
      
      coordsx = malloc (npoints * sizeof * coordsx);
      assert (npoints==0 || coordsx);
      coordsy = malloc (npoints * sizeof * coordsy);
      assert (npoints==0 || coordsy);
      coordsz = malloc (npoints * sizeof * coordsz);
      assert (npoints==0 || coordsz);
      coords[0] = coordsx;
      coords[1] = coordsy;
      coords[2] = coordsz;
      
      n = 0;
      for (k=0; k<array3d_npoints_k; ++k) {
        for (j=0; j<array3d_npoints_j; ++j) {
          for (i=0; i<array3d_npoints_i; ++i) {
            assert (n <= npoints);
            coordsx[n] = array3d_x0 + i * array3d_dx_i + j * array3d_dx_j + k * array3d_dx_k;
            coordsy[n] = array3d_y0 + i * array3d_dy_i + j * array3d_dy_j + k * array3d_dy_k;
            coordsz[n] = array3d_z0 + i * array3d_dz_i + j * array3d_dz_j + k * array3d_dz_k;
            ++n;
          }
        }
      }
      assert (n == npoints);
      
      inputs = malloc (nvars * sizeof * inputs);
      assert (inputs);
      
      for (n=0; n<nvars; ++n) {
        inputs[n] = CCTK_VarIndex (array3d_vars[n]);
        if (inputs[n] < 0) {
          inputs[n] = -1;
        }
      }
      
      output_types = malloc (nvars * sizeof * output_types);
      assert (output_types);
      
      for (n=0; n<nvars; ++n) {
        output_types[n] = CCTK_VARIABLE_REAL;
      }
      
      outputs = malloc (nvars * sizeof * outputs);
      assert (outputs);
      
      for (n=0; n<nvars; ++n) {
        outputs[n] = &arrays3d[npoints * n];
      }
      
      ierr = CCTK_InterpGridArrays
        (cctkGH, 3, interpolator, options_table, coord_handle,
         npoints, CCTK_VARIABLE_REAL, coords,
         nvars, inputs,
         nvars, output_types, outputs);
      assert (! ierr);
      
      free (coordsx);
      free (coordsy);
      free (coordsz);
      free (inputs);
      free (output_types);
      free (outputs);
    }
  }
  
  
  
  /* Parallel 1D Arrays */
  {
    nvars = nparrays1d;
    if (nvars > 0) {
      group = CCTK_GroupIndex ("InterpToArray::parrays1d");
      assert (group >= 0);
      ierr = CCTK_GroupDynamicData (cctkGH, group, &dyndata);
      assert (! ierr);
      assert (dyndata.dim == 1);
      npoints = dyndata.lsh[0];
      
      coordsx = malloc (npoints * sizeof * coordsx);
      assert (npoints==0 || coordsx);
      coordsy = malloc (npoints * sizeof * coordsy);
      assert (npoints==0 || coordsy);
      coordsz = malloc (npoints * sizeof * coordsz);
      assert (npoints==0 || coordsz);
      coords[0] = coordsx;
      coords[1] = coordsy;
      coords[2] = coordsz;
      
      n = 0;
      for (i=0; i<dyndata.lsh[0]; ++i) {
        assert (n <= npoints);
        coordsx[n] = (parray1d_x0 +
                      (dyndata.lbnd[0] + i) * parray1d_dx_i);
        coordsy[n] = (parray1d_y0 +
                      (dyndata.lbnd[0] + i) * parray1d_dy_i);
        coordsz[n] = (parray1d_z0 +
                      (dyndata.lbnd[0] + i) * parray1d_dz_i);
        ++n;
      }
      assert (n == npoints);
      
      inputs = malloc (nvars * sizeof * inputs);
      assert (inputs);
      
      for (n=0; n<nvars; ++n) {
        inputs[n] = CCTK_VarIndex (parray1d_vars[n]);
        assert (inputs[n] >= 0);
        if (inputs[n] < 0) {
          inputs[n] = -1;
        }
      }
      
      output_types = malloc (nvars * sizeof * output_types);
      assert (output_types);
      
      for (n=0; n<nvars; ++n) {
        output_types[n] = CCTK_VARIABLE_REAL;
      }
      
      outputs = malloc (nvars * sizeof * outputs);
      assert (outputs);
      
      for (n=0; n<nvars; ++n) {
        outputs[n] = &parrays1d[npoints * n];
      }
      
      operation_codes = malloc (nvars * sizeof * operation_codes);
      assert (operation_codes);
      
      for (n=0; n<nvars; ++n) {
        operation_codes[n] = parray1d_spacederivs[n];
      }
      
      ierr = Util_TableSetIntArray
        (options_table, nvars, operation_codes, "operation_codes");
      assert (! ierr);
      
      time_deriv_order = malloc (nvars * sizeof * time_deriv_order);
      assert (time_deriv_order);
      
      for (n=0; n<nvars; ++n) {
        time_deriv_order[n] = parray1d_timederivs[n];
      }
      
      ierr = Util_TableSetIntArray
        (options_table, nvars, time_deriv_order, "time_deriv_order");
      assert (! ierr);
      
      ierr = CCTK_InterpGridArrays
        (cctkGH, 3, interpolator, options_table, coord_handle,
         npoints, CCTK_VARIABLE_REAL, coords,
         nvars, inputs,
         nvars, output_types, outputs);
      assert (! ierr);
      
      ierr = Util_TableDeleteKey (options_table, "operation_codes");
      assert (! ierr);
      
      ierr = Util_TableDeleteKey (options_table, "time_deriv_order");
      assert (! ierr);
      
      free (coordsx);
      free (coordsy);
      free (coordsz);
      free (inputs);
      free (output_types);
      free (outputs);
      free (operation_codes);
      free (time_deriv_order);
    }
  }
  
  
  
  /* Parallel 2D Arrays */
  {
    nvars = nparrays2d;
    if (nvars > 0) {
      group = CCTK_GroupIndex ("InterpToArray::parrays2d");
      assert (group >= 0);
      ierr = CCTK_GroupDynamicData (cctkGH, group, &dyndata);
      assert (! ierr);
      assert (dyndata.dim == 2);
      npoints = dyndata.lsh[0] * dyndata.lsh[1];
      
      coordsx = malloc (npoints * sizeof * coordsx);
      assert (npoints==0 || coordsx);
      coordsy = malloc (npoints * sizeof * coordsy);
      assert (npoints==0 || coordsy);
      coordsz = malloc (npoints * sizeof * coordsz);
      assert (npoints==0 || coordsz);
      coords[0] = coordsx;
      coords[1] = coordsy;
      coords[2] = coordsz;
      
      n = 0;
      for (j=0; j<dyndata.lsh[1]; ++j) {
        for (i=0; i<dyndata.lsh[0]; ++i) {
          assert (n <= npoints);
          coordsx[n] = (parray2d_x0 +
                        (dyndata.lbnd[0] + i) * parray2d_dx_i +
                        (dyndata.lbnd[1] + j) * parray2d_dx_j);
          coordsy[n] = (parray2d_y0 +
                        (dyndata.lbnd[0] + i) * parray2d_dy_i +
                        (dyndata.lbnd[1] + j) * parray2d_dy_j);
          coordsz[n] = (parray2d_z0 +
                        (dyndata.lbnd[0] + i) * parray2d_dz_i +
                        (dyndata.lbnd[1] + j) * parray2d_dz_j);
          ++n;
        }
      }
      assert (n == npoints);
      
      inputs = malloc (nvars * sizeof * inputs);
      assert (inputs);
      
      for (n=0; n<nvars; ++n) {
        inputs[n] = CCTK_VarIndex (parray2d_vars[n]);
        if (inputs[n] < 0) {
          inputs[n] = -1;
        }
      }
      
      output_types = malloc (nvars * sizeof * output_types);
      assert (output_types);
      
      for (n=0; n<nvars; ++n) {
        output_types[n] = CCTK_VARIABLE_REAL;
      }
      
      outputs = malloc (nvars * sizeof * outputs);
      assert (outputs);
      
      for (n=0; n<nvars; ++n) {
        outputs[n] = &parrays2d[npoints * n];
      }
      
      ierr = CCTK_InterpGridArrays
        (cctkGH, 3, interpolator, options_table, coord_handle,
         npoints, CCTK_VARIABLE_REAL, coords,
         nvars, inputs,
         nvars, output_types, outputs);
      assert (! ierr);
      
      free (coordsx);
      free (coordsy);
      free (coordsz);
      free (inputs);
      free (output_types);
      free (outputs);
    }
  }
  
  
  
  /* Parallel 3D Arrays */
  {
    nvars = nparrays3d;
    if (nvars > 0) {
      group = CCTK_GroupIndex ("InterpToArray::parrays3d");
      assert (group >= 0);
      ierr = CCTK_GroupDynamicData (cctkGH, group, &dyndata);
      assert (! ierr);
      assert (dyndata.dim == 3);
      npoints = dyndata.lsh[0] * dyndata.lsh[1] * dyndata.lsh[2];
      
      coordsx = malloc (npoints * sizeof * coordsx);
      assert (npoints==0 || coordsx);
      coordsy = malloc (npoints * sizeof * coordsy);
      assert (npoints==0 || coordsy);
      coordsz = malloc (npoints * sizeof * coordsz);
      assert (npoints==0 || coordsz);
      coords[0] = coordsx;
      coords[1] = coordsy;
      coords[2] = coordsz;
      
      n = 0;
      for (k=0; k<dyndata.lsh[2]; ++k) {
        for (j=0; j<dyndata.lsh[1]; ++j) {
          for (i=0; i<dyndata.lsh[0]; ++i) {
            assert (n <= npoints);
            coordsx[n] = (parray3d_x0 +
                          (dyndata.lbnd[0] + i) * parray3d_dx_i +
                          (dyndata.lbnd[1] + j) * parray3d_dx_j +
                          (dyndata.lbnd[2] + k) * parray3d_dx_k);
            coordsy[n] = (parray3d_y0 +
                          (dyndata.lbnd[0] + i) * parray3d_dy_i +
                          (dyndata.lbnd[1] + j) * parray3d_dy_j +
                          (dyndata.lbnd[2] + k) * parray3d_dy_k);
            coordsz[n] = (parray3d_z0 +
                          (dyndata.lbnd[0] + i) * parray3d_dz_i +
                          (dyndata.lbnd[1] + j) * parray3d_dz_j +
                          (dyndata.lbnd[2] + k) * parray3d_dz_k);
            ++n;
          }
        }
      }
      assert (n == npoints);
      
      inputs = malloc (nvars * sizeof * inputs);
      assert (inputs);
      
      for (n=0; n<nvars; ++n) {
        inputs[n] = CCTK_VarIndex (parray3d_vars[n]);
        if (inputs[n] < 0) {
          inputs[n] = -1;
        }
      }
      
      output_types = malloc (nvars * sizeof * output_types);
      assert (output_types);
      
      for (n=0; n<nvars; ++n) {
        output_types[n] = CCTK_VARIABLE_REAL;
      }
      
      outputs = malloc (nvars * sizeof * outputs);
      assert (outputs);
      
      for (n=0; n<nvars; ++n) {
        outputs[n] = &parrays3d[npoints * n];
      }
      
      ierr = CCTK_InterpGridArrays
        (cctkGH, 3, interpolator, options_table, coord_handle,
         npoints, CCTK_VARIABLE_REAL, coords,
         nvars, inputs,
         nvars, output_types, outputs);
      assert (! ierr);
      
      free (coordsx);
      free (coordsy);
      free (coordsz);
      free (inputs);
      free (output_types);
      free (outputs);
    }
  }
  
  
  
  ierr = Util_TableDestroy (options_table);
  assert (! ierr);
}