/* $Header$ */
 
#include <assert.h>
#include <stdio.h>

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

#define REFLEVEL round(log10((CCTK_REAL)(cctkGH->cctk_levfac[0]))/log10(2.0))

void PunctureTracker (CCTK_ARGUMENTS);
void PunctureTracker_init (CCTK_ARGUMENTS);


void PunctureTracker_init (CCTK_ARGUMENTS)
{
  DECLARE_CCTK_PARAMETERS;
  DECLARE_CCTK_ARGUMENTS;
  int i;

  if(pt_verbose>0)
    CCTK_INFO("Initializing PunctureTracker");

  if((pt_initial_x[0]!=0.0)||(pt_initial_y[0]!=0.0)||(pt_initial_z[0]!=0.0))
    for (i=0;i<pt_num_tracked;i++){
      pt_loc_x_p[i] = pt_initial_x[i];
      pt_loc_y_p[i] = pt_initial_y[i];
      pt_loc_z_p[i] = pt_initial_z[i];
    }
  else{
    pt_loc_x_p[0] = par_b;
    pt_loc_x_p[1] =-par_b;
    pt_loc_y_p[0] = 0.0;
    pt_loc_y_p[1] = 0.0;
    pt_loc_z_p[0] = 0.0;
    pt_loc_z_p[1] = 0.0;
  }
  if(pt_verbose>0)
    CCTK_INFO("done initializing PunctureTracker");
}

void PunctureTracker (CCTK_ARGUMENTS)
{
  DECLARE_CCTK_PARAMETERS;
  DECLARE_CCTK_ARGUMENTS;
  
  int n,ierror,i,numpoints;
  // portland doesnt like that
  //CCTK_REAL  pt_betax[pt_num_tracked], pt_betay[pt_num_tracked], pt_betaz[pt_num_tracked];
  //so I do that

  CCTK_REAL * restrict const pt_betax = malloc(pt_num_tracked*sizeof(CCTK_REAL));
  CCTK_REAL * restrict const pt_betay = malloc(pt_num_tracked*sizeof(CCTK_REAL));
  CCTK_REAL * restrict const pt_betaz = malloc(pt_num_tracked*sizeof(CCTK_REAL));

  int         operator_handle, param_table_handle, coordsys_handle;
  int         input_array_indices[3];
  const void *interp_coords[3];
  void       *output_arrays[3];

  const CCTK_INT input_array_type_codes[3] = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL};   
  const CCTK_INT output_array_type_codes[3]= { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL};
   
  if(pt_verbose >0)
    CCTK_INFO("tracking Punctures");
  
  interp_coords[0] = (const void *) pt_loc_x_p;
  interp_coords[1] = (const void *) pt_loc_y_p;
  interp_coords[2] = (const void *) pt_loc_z_p;
  
  operator_handle     = -1;
  coordsys_handle     = -1;
  param_table_handle  = -1;
  
  // set interpolation handles
  operator_handle = CCTK_InterpHandle("Lagrange polynomial interpolation");
  if (operator_handle < 0) CCTK_WARN(0, "can’t get interpolation handle!");

  param_table_handle = Util_TableCreateFromString("order=4");
  if (param_table_handle < 0) CCTK_WARN(0, "can't ge parameter table!");
  
  coordsys_handle = CCTK_CoordSystemHandle("cart3d");
  if (coordsys_handle < 0) CCTK_WARN(0,"can't get coordsys handle!");
  if(pt_verbose>0)
    for (i=0; i<pt_num_tracked; i++) {
      CCTK_VInfo(CCTK_THORNSTRING,
      "Puncture %d is at (x,y,z)=(%f,%f,%f)\n",i,pt_loc_x_p[i],pt_loc_y_p[i],pt_loc_z_p[i]);
  }
  
  output_arrays[0] = (void *) pt_betax;
  output_arrays[1] = (void *) pt_betay;
  output_arrays[2] = (void *) pt_betaz;

  input_array_indices[0] = CCTK_VarIndex ("ADMBase::betax");
  input_array_indices[1] = CCTK_VarIndex ("ADMBase::betay");
  input_array_indices[2] = CCTK_VarIndex ("ADMBase::betaz");

  //only proc 0 does interpolation
  numpoints = 0;
  if (CCTK_MyProc(cctkGH) == 0)
    numpoints = pt_num_tracked;

  ierror = CCTK_InterpGridArrays(cctkGH,3,
                                 operator_handle, param_table_handle,coordsys_handle,
                                 numpoints,
                                 CCTK_VARIABLE_REAL,
                                 interp_coords,
                                 3,input_array_indices,
                                 3,output_array_type_codes,
                                  output_arrays);
  if (ierror<0)
    CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
               "error return from interpolator ierror=%d",(int)ierror);
  

  for (n=0; n< pt_num_tracked; n++) {
    if(pt_verbose>1)
      printf("at the location of puncture %i:\n              is pt_betax=%f, pt_betay=%f,pt_betaz=%f \n",n,pt_betax[n],pt_betay[n],pt_betaz[n]);
        pt_shiftx[n] = pt_betax[n];
        pt_shifty[n] = pt_betay[n];
        pt_shiftz[n] = pt_betaz[n];
  }

    free(pt_betax,pt_betay,pt_betaz);
  if(pt_verbose>0)
    CCTK_INFO ("done");
}

void update_punc_loc(CCTK_ARGUMENTS){

  DECLARE_CCTK_PARAMETERS;
  DECLARE_CCTK_ARGUMENTS;

  CCTK_REAL pt_dt;
  int n;

  pt_dt = CCTK_DELTA_TIME;

  for (n=0; n< pt_num_tracked; n++) {
    //correct small errors
    //    if (pt_shiftx[n] < 10e-14) pt_shiftx[n] = 0.0;
    //    if (pt_shifty[n] < 10e-14) pt_shifty[n] = 0.0;
    //    if (pt_shiftz[n] < 10e-14) pt_shiftz[n] = 0.0;
    
    if(pt_verbose > 2){
      printf("updating puncture %i\n",n);
      printf("   pt_dt      is %f\n",pt_dt);
      printf("   pt_loc_x   is %f\n",pt_loc_x[n]);
      printf("   pt_loc_x_p is %f\n",pt_loc_x_p[n]);
      printf("   pt_shiftx  is %f\n",pt_shiftx[n]);
      printf("   pt_dt      is %f\n",pt_dt);
      fflush(stdout);
    }
     
     printf("WRITING TO SURFACE %d\n",pt_which_surface_to_take[n]);
    pt_loc_x[n] = pt_loc_x_p[n] - pt_shiftx[n] * pt_dt;
    pt_loc_y[n] = pt_loc_y_p[n] - pt_shifty[n] * pt_dt;
    pt_loc_z[n] = pt_loc_z_p[n] - pt_shiftz[n] * pt_dt;

    sf_centroid_x[pt_which_surface_to_take[n]]=pt_loc_x_p[n];    
    sf_centroid_y[pt_which_surface_to_take[n]]=pt_loc_y_p[n];
    sf_centroid_z[pt_which_surface_to_take[n]]=pt_loc_z_p[n];
     
    sf_valid[pt_which_surface_to_take[n]]=1.0;
    sf_valid[pt_which_surface_to_take[n]]=1.0;
    sf_valid[pt_which_surface_to_take[n]]=1.0;
  }
}