/* $Header$ */

#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 nvars;
  int npoints;
  
  int n;
  int i, j, k;
  int d;
  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 (coordsx);
      coordsy = malloc (npoints * sizeof * coordsy);
      assert (coordsy);
      coordsz = malloc (npoints * sizeof * coordsz);
      assert (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);
    
      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 (coordsx);
      coordsy = malloc (npoints * sizeof * coordsy);
      assert (coordsy);
      coordsz = malloc (npoints * sizeof * coordsz);
      assert (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 (coordsx);
      coordsy = malloc (npoints * sizeof * coordsy);
      assert (coordsy);
      coordsz = malloc (npoints * sizeof * coordsz);
      assert (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);
    
      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 (coordsx);
      coordsy = malloc (npoints * sizeof * coordsy);
      assert (coordsy);
      coordsz = malloc (npoints * sizeof * coordsz);
      assert (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);
    
      free (coordsx);
      free (coordsy);
      free (coordsz);
      free (inputs);
      free (output_types);
      free (outputs);
    }
  }
  
  
  
  ierr = Util_TableDestroy (options_table);
  assert (! ierr);
}