/** @file grid.c

 *  @brief General grid functions common for cdr and poisson grids.

 */

#include <stdlib.h>

#include <stdio.h>

#include <math.h>

#include "grid.h"

#include "parameters.h"

#include "proto.h"

#include "species.h"

/** @brief Gets the maximum depth of a tree (the level of its deepest node) */

int

grid_max_depth_r (grid_t *grid)

{

  int l, lc;

  grid_t *child;

  l = grid->level;

  iter_childs (grid, child) {

    lc = grid_max_depth_r (child);

    l = (l > lc)? l: lc;

  }

  return l;

}

/** @brief Gets the minimum r of a grid and all its descendants */

double

grid_rmin_r (grid_t *grid)

{

  int lev;

  double minr, cminr;

  grid_t *child;

  lev = grid->level;

  minr = r_at (grid->r0, lev);

  iter_childs (grid, child) {

    cminr = grid_rmin_r (child);

    minr = (minr > cminr)? cminr: minr;

  }

  return minr;

}

/** @brief Recursively prints one grid and all its descendants. */

void

grid_print_r (grid_t *grid, int indent)

{

  int i;

  grid_t *child;

  for (i = 0; i < indent; i++) fputs ("  ", stdout);

  printf ("grid (" grid_printf_str") {\n", grid_printf_args(grid));

  iter_childs (grid, child) {

    grid_print_r (child, indent + 1);

  }

  for (i = 0; i < indent; i++) fputs ("  ", stdout);

  fputs ("}\n", stdout);

}

/** @brief Checks whether a point is in the interior or on a given boundary

 * of a grid.

 *

 * Does the grid contain the point (i, j) (in its own coordinates)?\n

 * @a check has the following meaning:\n

 *  GRID_INSIDE : Checks if the point is in the interior of the grid\n

 *  BND_MASK_*  : Checks if it is on the given boundary (BND_MASK_ALL checks for

 *               any boundary).\n

 */

int

grid_contains (grid_t *grid, int i, int j, int check)

{

  int r0, r1, z0, z1;

  int r;

  r0 = grid->r0;

  z0 = grid->z0;

  r1 = grid->r1 - 1;

  z1 = grid->z1 - 1;

  r = FALSE;

  if (check & BND_MASK_LEFT) {

    r0 -= 1;

  }

  if (check & BND_MASK_RIGHT) {

    r1 += 1;

  }

  if (check & BND_MASK_BOTTOM) {

    z0 -= 1;

  }

  if (check & BND_MASK_TOP) {

    z1 += 1;

  }

  if (r0 < i && i < r1 && z0 < j && j < z1) 

    r = TRUE;

  if (check & GRID_INSIDE) {

    return r;

  } else {

    return (r && !grid_contains (grid, i, j, GRID_INSIDE));

  }

}

/** @brief Checks whether two grids overlap and return the rectangle that they

 * in common.

 *

 * Checks whether two grids overlap and return the rectangle that they 

 * have in common in the coordinates of the finest grid in 

 * (*left, *bottom) - (*right, *top)\n

 *

 * Note that *left, *bottom... can be modified even if the two grids do

 * not overlap.

 * @a buff1 and @a buff2 are integers that specify the numbers of cells

 * that are taken as buffer in the boundaries for @a grid1 and @a grid2.

 */

int

grid_overlap (grid_t *grid1, grid_t *grid2, int buff1, int buff2, 

	      int *left, int *bottom, int *right, int *top, int *level_diff)

{

  return grid_overlap_with_shifts (grid1, grid2, buff1, buff2, 

				  left, bottom, right, top, level_diff,

				  0, 0);

}

/** @brief Checks the overlap between two grids.

 *

 * Checks the overlap between two grids before it shifts the

 * coarser grid by (shift_r << (level_diff - 1), shift_z << (level_diff - 1)

 * (in finer grid coordinates).

 *

 *  Q: WTF! WHY!!!? \n

 *  A: We use this to implement the mapping between electric fields.  Since

 *     The electric field is evaluated not in the cell centers but in the

 *     borders, when we interpolate we have to take a shift into account.

 *     See mapper.c and mapper.h for more info.

 */

int

grid_overlap_with_shifts (grid_t *grid1, grid_t *grid2, int buff1, int buff2, 

			 int *left, int *bottom, int *right, int *top, 

			 int *level_diff, int shift_r, int shift_z)

{

  int sign;

  debug (3, "grid_overlap(" grid_printf_str", " grid_printf_str 

	 ", buff1 = %d, buff2 = %d)\n",

	 grid_printf_args(grid1), grid_printf_args(grid2),

	 buff1, buff2);

  sign = 1;

  /* If level_diff < 0, grid2 is finer */

  *level_diff = grid1->level - grid2->level;

  if (*level_diff < 0) {

    sign = -1;

    *level_diff = - *level_diff;

    XCHG (grid1, grid2);

    XCHG (buff1, buff2);

  }

  /* At this point, level_diff >= 0 and grid1 is finer 

     or equivalent to grid2 */

  *left = MAX_AT_LEVEL_WITH_SHIFT (grid1->r0 - buff1, grid2->r0 - buff2, 

				   *level_diff, shift_r);

  *right = MIN_AT_LEVEL_WITH_SHIFT (grid1->r1 + buff1, grid2->r1 + buff2, 

				    *level_diff, shift_r);

  /* Note: *right is not contained in the grid.  Hence, if *left == *right

     the two grids do not overlap.  Same for *bottom and *top, below. */

  debug (3, "*left = %d, *right = %d\n", *left, *right);

  if (*left >= *right) {

    *level_diff *= sign;

    return FALSE;

  }

  *bottom = MAX_AT_LEVEL_WITH_SHIFT (grid1->z0 - buff1, grid2->z0 - buff2, 

				     *level_diff, shift_z);

  *top = MIN_AT_LEVEL_WITH_SHIFT (grid1->z1 + buff1, grid2->z1 + buff2, 

				  *level_diff, shift_z);

  debug (3, "*bottom = %d, *top = %d\n", *bottom, *top);

  if (*bottom >= *top) {

    *level_diff *= sign;

    return FALSE;

  }

  debug (3, "pois_overlap = TRUE\n");

  *level_diff *= sign;

  return TRUE;

}

/** @brief Makes a grid inherit the ext_bound field of its parent,

 *

 * with the appropiate modifications.

 */

void

grid_inherit_ext_bound (grid_t *grid)

{

  grid_t *parent = grid->parent;

  debug (3, "grid_inherit_ext_bound (...)\n");

  grid->ext_bound = BND_NONE;

  /* If the parent does not have any external boundary, the child will also

     not have */

  if (BND_NONE == parent->ext_bound) {

    return;

  }

  if (parent->ext_bound & BND_MASK_LEFT) {

    if ((parent->r0 << 1) == grid->r0) grid->ext_bound |= BND_MASK_LEFT;

  }

  if (parent->ext_bound & BND_MASK_RIGHT) {

    if ((parent->r1 << 1) == grid->r1) grid->ext_bound |= BND_MASK_RIGHT;

  }

  if (parent->ext_bound & BND_MASK_BOTTOM) {

    if ((parent->z0 << 1) == grid->z0) grid->ext_bound |= BND_MASK_BOTTOM;

  }

  if (parent->ext_bound & BND_MASK_TOP) {

    if ((parent->z1 << 1) == grid->z1) grid->ext_bound |= BND_MASK_TOP;

  }

}

/** @brief Finds the finest descendant of grid that contains the point given by

 * the (@a r, @a z) coordinates.

 *

 * Or NULL if the point is not contained in grid.

 * We assume that only ONE sub-grid can contain a given point.

 */

grid_t *

grid_finest_containing_r (grid_t *grid, double r, double z)

{

  grid_t *leaf, *ret;

  if (r < er_r_at(grid->r0 - 1, grid->level) || 

      r > er_r_at(grid->r1 - 1, grid->level))

    return NULL;

  if (z < ez_z_at(grid->z0 - 1, grid->level) || 

      z > ez_z_at(grid->z1 - 1, grid->level))

    return NULL;

  iter_childs (grid, leaf) {

    ret = grid_finest_containing_r (leaf, r, z);

    if (NULL != ret)

      return ret;

  }

  return grid;

}

/** @brief How many children does one grid have? */

int

grid_howmany_children (grid_t *grid)

{

  grid_t *leaf;

  int count;

  count = 0;

  iter_childs (grid, leaf) {

    count ++;

  }

  return count;

}

/** @brief Finds child number @a n of the given grid.

 *

 * If it doesn't have enough children just returns NULL

 */

grid_t*

grid_get_child (grid_t *grid, int n)

{

  grid_t *leaf;

  int count;

  count = 0;

  iter_childs (grid, leaf) {

    if (count == n) return leaf;

    count ++;

  }

  return NULL;

}