/** @file poisson.c
 *  @brief Poisson/Helmholtz solver, including the routines for
 *  inhomogeneous (point-plane) configurations.
 */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "cdr.h"
#include "cstream.h"
#include "fishpack.h"
#include "grid.h"
#include "interpol2.h"
#include "mapper.h"
#include "parameters.h"
#include "poisson.h"
#include "proto.h"
#include "rz_array.h"
#include "species.h"

static void get_bound_cond (pois_grid_t *grid, pois_problem_t *prob,
			    int *r_bound_cond, int *z_bound_cond, 
			    double lambda);
static int needs_refinement (pois_grid_t *grid, int ir, int iz, 
			     double threshold);
static int refine_in (pois_grid_t *grid, int cr0, int cz0, int cr1, int cz1);

decl_mapper_funcs(er);
decl_mapper_funcs(ez);
decl_mapper_funcs(etheta);
decl_mapper_funcs(charge);
/*!< Here we set three mapper_t objects that specify how we should do
 * the mapping between the potential (once it is calculated) and the three
 * components of the electric field.
 */

/* Global variables for the computations of inhomogeneous fields:
 *  initialized in pois_inhom_init. */

int pois_gridpoints_z, pois_gridpoints_r;
/*!< Number of gridpoints at level 0 of the poisson solver
 *
 * if !pois_inhom, this has to be equal to the parameter gridpoins_z.
 * We also allow a different extension in r, albeit this is not yet used.
 */

double pois_inhom_L;    /*!< Total length of the poisson domain. */
double pois_inhom_z;    /*!< Location of the floating charge. */
double pois_inhom_dphi;
/*!< Potential difference between the electrodes created by a unit charge 
 *  located at (0, pois_inhom_z)
 */

double pois_inhom_fixed_q_t;
/*!< Time-dependent fixed charge for the inhomogeneous poisson code
 *  It is updated in cstream_set_field_at_time acording to the parameter
 *  rise_time.
 */

extern double E_x, E_y, E_z;

static double inhom_phi_term (double r, double z, double b);
static double inhom_e_term (double r, double z, double b);
static double inhom_er_term (double r, double z, double b);
static double inhom_ez_term (double r, double z, double b);
/*!< To calculate the potential and electric fields created by a unit charge
 *  between two planar electrodes, we use the "mirror charges method"
 *  the contribution of each mirror charge is given by these functions.
 */

mapper_t er_mapper = mk_mapper_staggered (er, &interpol_bilin, 0, -1);
mapper_t ez_mapper = mk_mapper_staggered (ez, &interpol_bilin, -1, 0);
mapper_t etheta_mapper = mk_mapper (etheta, &interpol_quadratic);
mapper_t *e_mappers[] = {&er_mapper, &ez_mapper, &etheta_mapper, NULL};
/*!< These are the mapper structures to translate between cdr and poisson
 *  grids.  See mapper.c and mapper.h for more info.
 */

mapper_t charge_mapper = mk_mapper_down(charge, &interpol_wackers);
mapper_t *charge_mappers[] = {&charge_mapper, NULL};

pois_problem_t *pois_electrostatic;

/** @brief  All initialization of the Poisson solver has to be here. */
void
pois_init (void)
{
  if (pois_inhom) {
    pois_inhom_init ();
  }
  else {
    pois_gridpoints_z = gridpoints_z;
  }
  pois_gridpoints_r = gridpoints_r;

  pois_electrostatic = xmalloc (sizeof (pois_problem_t));

  pois_electrostatic->max_level = pois_max_level;
  pois_electrostatic->extra_levels = extra_pois_levels;
  pois_electrostatic->max_error = pois_max_error;

  pois_electrostatic->bnd_right = pois_bnd_right;
  pois_electrostatic->bnd_top = pois_bnd_top;
  pois_electrostatic->bnd_bottom = pois_bnd_bottom;
}

/** @brief Creates a 2D Poisson grid. */
pois_grid_t*
pois_new_a (int r0, int z0, int r1, int z1)
{
  pois_grid_t *grid;

  grid = pois_new_3d_a (r0, z0, r1, z1, 1);

  return grid;
}

/** @brief Creates a new 3D Poisson grid.
 *
 * @a r0, @a z0, @a r1, @a z1 are in GRID units
 */
pois_grid_t*
pois_new_3d_a (int r0, int z0, int r1, int z1, int ntheta)
{
  rz_array_t *phi, *charge, *error;
  pois_grid_t *grid;

  debug (2, "pois_new_a (r0 = %d, z0 = %d, r1 = %d, z1 = %d)\n", 
	 r0, z0, r1, z1);

  phi = rz_new_3d_a (r0, z0, r1, z1, ntheta);
  error = rz_new_3d_a (r0, z0, r1, z1, ntheta);

#ifdef F_PHI_SAME_AS_CHARGE
  charge = phi;
#else
  charge = rz_new_3d_a (r0, z0, r1, z1, ntheta);
#endif /* F_PHI_SAME_AS_CHARGE */

  grid = (pois_grid_t *) xmalloc (sizeof (pois_grid_t));

  grid->r0 = r0;
  grid->r1 = r1;
  grid->z0 = z0;
  grid->z1 = z1;
  grid->ntheta = ntheta;

  grid->phi = phi;
  grid->charge = charge;
  grid->error = error;

  /* By default, the array does not have any external boundary. */
  grid->ext_bound = BND_NONE;

  /* Initially, a grid does not have any relatives. */
  init_leaf (grid);

  return grid;
}

/** @brief  Frees the memory allocated by pois_new_a */
void
pois_free (pois_grid_t *grid)
{
  debug (3, "pois_free (...)\n");

  /* The decission of whether to use the same memory for the charge and
   * the potential can be made elsewhere.
   */
  if (grid->phi != grid->charge)
    rz_free (grid->charge);
  rz_free (grid->phi);
  rz_free (grid->error);

  free (grid);
}

/** @brief Recursively frees a tree of Poisson grids. */
void
pois_free_r (pois_grid_t *grid)
{
  pois_grid_t *leaf;

  debug (3, "pois_free_r (" grid_printf_str ")\n", grid_printf_args (grid));

  free_childs (grid, leaf, pois_free_r);

  pois_free (grid);
}

/** @brief Creates a new Poisson grid.
 *
 * @a r0, @a r1... are in global units (i.e.  independent of the grid level). 
 */
pois_grid_t*
pois_new_glob_a (int r0, int z0, int r1, int z1, int level)
{
  debug (3, "pois_new_glob_a (%d, %d, %d, %d)\n", r0, r1, z0, z1);
  pois_grid_t *grid;

  if (level >= 0) {
    grid = pois_new_a (r0 << level, z0 << level, 
		       r1 << level, z1 << level);
  } else {
    level = -level;
    grid = pois_new_a (r0 >> level, z0 >> level, 
		       r1 >> level, z1 >> level);
  }
  return grid;
}

/** @brief Starts the tree of Poisson grids with the two coarsest ones.
 *
 * Receives the grid dimensions at level 0.
 */
pois_grid_t*
pois_init_tree_a (int r0, int z0, int r1, int z1)
{
  pois_grid_t *coarsest, *sub_coarsest;
  int level;

  debug (3, "pois_init_tree_a (%d, %d, %d, %d)\n", r0, z0, r1, z1);

  level = -extra_pois_levels;

  coarsest = pois_new_glob_a (r0, z0, r1, z1, level);
  coarsest->ext_bound = BND_MASK_ALL;
  coarsest->level = level;

  sub_coarsest = pois_new_glob_a (r0, z0, r1, z1, level + 1);

  add_child (coarsest, sub_coarsest);
  grid_inherit_ext_bound ((grid_t *) sub_coarsest);

  return coarsest;
}

/** @brief Is it possible to use a single function to handle all four
 *  boundaries?
 *
 * Maybe this should be rewritten to use the interpol2 module? 
 */
REAL*
pois_boundary_a (pois_grid_t *grid, int boundary)
{
  pois_grid_t *parent;

  /* Forget the 0s, they are there to shut up the compiler. */
  int stride = 0, ortho_stride = 0, icells, invert, n;
  int rb, zb;
  double interp_array[3][3] = {{-4.0, 1.0, 18.0},
			       {-2.0, 78.0, 14.0},
  			       {-6.0, -7.0, 4.0}};

  const double prefactor = 1 / 96.0;
  REAL *bnd, *p, *lstart;

  debug (3, "pois_boundary_a (..., %d)\n", boundary);

  parent = grid->parent;

  /* First we check if the boundary is in r or in z. */
  if (boundary & BND_AT_Z) {
    icells = grid->r1 - grid->r0;
    zb = 1; rb = 0;
    if (NULL != parent) {
      stride = parent->phi->strides[R_INDX];
      ortho_stride = parent->phi->strides[Z_INDX];
    }
  } else {
    zb = 0; rb = 1;
    icells = grid->z1 - grid->z0;
    if (NULL != parent) {
      stride = parent->phi->strides[Z_INDX];
      ortho_stride = parent->phi->strides[R_INDX];
    }
  }

  debug (3, "   grid->ext_bound = %d\n", grid->ext_bound);

  if (grid->ext_bound & BND_MASK (boundary)) {
    /* If the boundary is external, we apply homogeneous Dirichlet/Neumann. */
    bnd = (REAL *) xcalloc (sizeof (REAL), icells);
    debug (3, "pois_boundary_a (...) -> [0.0] * %d\n", icells);

    assert (grid->r1 != 0);
    return bnd;
  } else {
    bnd = (REAL *) xmalloc (sizeof (REAL) * icells);
  }
  /* If the grid is root it should have returned in the previous if */
  assert (NULL != parent);

  /* Now we check if the boundary is at {r,z} minimum or {r,z} maximum */
  if (boundary & BND_MAX) {
    lstart = RZP (parent->phi, 
		(grid->r1 >> 1) - 1 + rb, 
		(grid->z1 >> 1) - 1 + zb);
    stride = -stride;
    ortho_stride = -ortho_stride;
    invert = TRUE;
  } else {
    lstart = RZP (parent->phi, (grid->r0 >> 1) - rb, (grid->z0 >> 1) - zb);
    invert = FALSE;
  }

  for (p = lstart, n = 0; n < icells; p += stride, n += 2) {
    /* p points to the cell in the parent grid and p + stride points to 
       an adjacent cell. */
    double s = 0.0, s2 = 0.0;
    int i, j;
    for (i = -1; i <= 1; i++) {
      for (j = -1; j <= 1; j++) {
	s += *(p + i * stride + j * ortho_stride) * interp_array[i + 1][j + 1];
	s2 += *(p - i * stride + j * ortho_stride) * interp_array[i + 1][j + 1];
      }
    }
    if (!invert) {
      bnd[n] = prefactor * s;
      if (n < icells - 1) bnd[n + 1] = prefactor * s2;
    } else {
      bnd[icells - n - 1] = prefactor * s;
      if (n < icells - 1) bnd[icells - n - 2] = prefactor * s2;
    }
  }
  return bnd;
}

/** @brief Solves the Poisson equation, calling the FISHPACK routine.
 *
 * Assumes that grid->charge is already set.
 */
void
pois_solve_grid (pois_grid_t *grid, pois_problem_t *prob, 
		 double lambda, double s) 
{
  REAL *boundaries[4];
  int i, z_bound_cond, r_bound_cond;
  double rmin, rmax, zmin, zmax;

  debug (3, "pois_solve_grid (..., lambda = %f, s = %f)\n", lambda, s);

  for (i = 0; i < 4; i++) {
    boundaries[i] = pois_boundary_a (grid, i);
    debug (3, "boundaries[%d] = %p\n", i, boundaries[i]);
  }

#ifndef F_PHI_SAME_AS_CHARGE
  rz_copy (grid->charge, grid->r0, grid->z0,
	   grid->phi, grid->r0, grid->z0,
	   grid->r1 - grid->r0, grid->z1 - grid->z0);
#endif

  rmin = grid->r0 * dr[grid->level];
  rmax = grid->r1 * dr[grid->level];
  zmin = grid->z0 * dz[grid->level];
  zmax = grid->z1 * dz[grid->level];

  debug (3, "HSTCYL/HSTCRT in a grid %d x %d\n", grid->r1 - grid->r0,
	 grid->z1 - grid->z0);

  get_bound_cond (grid, prob, &r_bound_cond, &z_bound_cond, lambda);

#ifndef TRUE2D
  fish_hstcyl (rmin, rmax, grid->r1 - grid->r0, 
	       r_bound_cond, boundaries[BND_LEFT], boundaries[BND_RIGHT],
	       zmin, zmax, grid->z1 - grid->z0,
	       z_bound_cond, boundaries[BND_BOTTOM], boundaries[BND_TOP],
	       lambda, s,
	       RZP (grid->phi, grid->r0, grid->z0), 
	       grid->phi->strides[Z_INDX]);
#else
  fish_hstcrt (rmin, rmax, grid->r1 - grid->r0, 
	       r_bound_cond, boundaries[BND_LEFT], boundaries[BND_RIGHT],
	       zmin, zmax, grid->z1 - grid->z0,
	       z_bound_cond, boundaries[BND_BOTTOM], boundaries[BND_TOP],
	       lambda,
	       RZP (grid->phi, grid->r0, grid->z0), 
	       grid->phi->strides[Z_INDX]);
#endif
  debug (3, "Finished HSTCYL/HSTCRT in a grid %d x %d\n", grid->r1 - grid->r0,
	 grid->z1 - grid->z0);

  pois_set_phi_boundaries (grid, boundaries, 
			   (grid->ext_bound & BND_MASK (BND_LEFT))
			    && 0. == lambda,
			   
			   prob->bnd_right == BND_CND_HNEUMANN
			   && (grid->ext_bound & BND_MASK (BND_RIGHT)),
			   
			   prob->bnd_bottom == BND_CND_HNEUMANN
			   && (grid->ext_bound & BND_MASK (BND_BOTTOM)),
			   
			   prob->bnd_top == BND_CND_HNEUMANN
			   && (grid->ext_bound & BND_MASK (BND_TOP)));

  for (i = 0; i < 4; i++) {
    debug (3, "boundaries[%d] = %p\n", i, boundaries[i]);
    free (boundaries[i]);
  }
  debug (3, "  <- pos_solve_grid (...)\n");
}

/** @brief Gets the boundary conditions that we pass to FISHPACK.
 *
 * Note that we use the same conditions for the Poisson equation
 * and for the Helmholtz equation used to find the photoionization
 * source.  From my point of view, it doesn't make sense to use different
 * conditions, so that will simply add an unneccesary complication.
 */
static void
get_bound_cond (pois_grid_t *grid, pois_problem_t *prob,
		int *r_bound_cond, int *z_bound_cond, double lambda)
{
  static int right_inside[] = {FISH_UNS_DIR, FISH_NEU_DIR, FISH_DIR_DIR};
  static int right_ext_neu[] = {FISH_UNS_NEU, FISH_NEU_NEU, FISH_DIR_NEU};
  static int *right_ext_dir = right_inside;
  static int fish_order[4] = { FISH_DIR_DIR, FISH_DIR_NEU, 
			       FISH_NEU_DIR, FISH_NEU_NEU};

  int bottom, top;
  int *right_sel;

  if ((grid->ext_bound & BND_MASK (BND_TOP))
      && prob->bnd_top == BND_CND_HNEUMANN)
    top = 1;
  else top = 0;

  if ((grid->ext_bound & BND_MASK (BND_BOTTOM))
      && prob->bnd_bottom == BND_CND_HNEUMANN)
    bottom = 1;
  else bottom = 0;

  *z_bound_cond = fish_order[(bottom << 1) + top];

  if ((grid->ext_bound & BND_MASK (BND_RIGHT))) {
    if (prob->bnd_right == BND_CND_HNEUMANN)
      right_sel = right_ext_neu;
    else 
      right_sel = right_ext_dir;
  } else {
    right_sel = right_inside;
  }

#ifndef TRUE2D
  *r_bound_cond = ((grid->ext_bound & BND_MASK (BND_LEFT)) && 0. == lambda)? 
		   right_sel[0]: right_sel[2];
#else
  *r_bound_cond = ((grid->ext_bound & BND_MASK (BND_LEFT))? 
		   right_sel[1]: right_sel[2]);
#endif

}

/** @brief Uses the @a boundaries[] vectors to set the boundaries by
 *  interpolation for the phi array of a grid.
 *
 *  (in the extra space we allocated when we created the array).
 *  xxx_neu tells us if the xxx boundary has Neumann b.c.
 *  (or unspecified, as FISHPACK calls them if they are applied in the axis)
 *
 *  TODO: Please rewrite this mess.
 */
void
pois_set_phi_boundaries (pois_grid_t *grid, REAL *boundaries[], 
			 int left_neu, int right_neu, int bottom_neu, 
			 int top_neu)
{
  int i;

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

  if ((grid->ext_bound & BND_MASK(BND_BOTTOM)) && bottom_neu) {
    for (i = grid->r0; i < grid->r1; i++) {
      RZ (grid->phi, i, grid->z0 - 1) = RZ (grid->phi, i, grid->z0);
    }
  } else {
    for (i = grid->r0; i < grid->r1; i++) {
      RZ (grid->phi, i, grid->z0 - 1) = 
	2 * boundaries[BND_BOTTOM][i - grid->r0] - RZ (grid->phi, i, grid->z0);
    }
  }

  if ((grid->ext_bound & BND_MASK(BND_TOP)) && top_neu) {
    for (i = grid->r0; i < grid->r1; i++) {
      RZ (grid->phi, i, grid->z1) = RZ (grid->phi, i, grid->z1 - 1);
    }
  } else {
    for (i = grid->r0; i < grid->r1; i++) {
      RZ (grid->phi, i, grid->z1) = 2 * boundaries[BND_TOP][i - grid->r0] 
	- RZ(grid->phi, i, grid->z1 - 1);
    }
  }

  if (left_neu) {
    for (i = grid->z0; i < grid->z1; i++) {
      RZ (grid->phi, grid->r0 - 1, i) = RZ(grid->phi, grid->r0, i);
    }
  } else {
    for (i = grid->z0; i < grid->z1; i++) {
      RZ(grid->phi, grid->r0 - 1, i) = 2 * boundaries[BND_LEFT][i - grid->z0] 
	- RZ(grid->phi, grid->r0, i);
    }
  }

  if (right_neu) {
    for (i = grid->z0; i < grid->z1; i++) {
      RZ (grid->phi, grid->r1, i) = RZ(grid->phi, grid->r1 - 1, i);
    }
  } else {
    for (i = grid->z0; i < grid->z1; i++) {      
      RZ(grid->phi, grid->r1, i) = 2 * boundaries[BND_RIGHT][i - grid->z0] 
	- RZ(grid->phi, grid->r1 - 1, i);
    }
  }

  /*  This is quite nightmarish.  Think about rewritting. */
  if ((grid->ext_bound & BND_MASK(BND_RIGHT)) && right_neu) {
    RZ(grid->phi, grid->r1, grid->z0 - 1) = RZ(grid->phi, grid->r1 - 1, 
						   grid->z0 - 1);
    RZ(grid->phi, grid->r1, grid->z1) = RZ(grid->phi, grid->r1 - 1, 
					       grid->z1);
  } else {
    RZ(grid->phi, grid->r1, grid->z0 - 1) =
      RZ(grid->phi, grid->r1 - 1, grid->z0 - 1) 
      + RZ(grid->phi, grid->r1, grid->z0)
      - 0.5 * (RZ(grid->phi, grid->r1 - 2, grid->z0 - 1) 
	       + RZ(grid->phi, grid->r1, grid->z0 + 1));

    RZ(grid->phi, grid->r1, grid->z1) =
      RZ(grid->phi, grid->r1 - 1, grid->z1) 
      + RZ(grid->phi, grid->r1, grid->z1 - 1)
      - 0.5 * (RZ(grid->phi, grid->r1 - 2, grid->z1) 
	       + RZ(grid->phi, grid->r1, grid->z1 - 2));
  }

  if ((grid->ext_bound & BND_MASK(BND_LEFT)) && left_neu) {
    RZ(grid->phi, grid->r0 - 1, grid->z0 - 1) = RZ(grid->phi, grid->r0, 
						   grid->z0 - 1);
    RZ(grid->phi, grid->r0 - 1, grid->z1) = RZ(grid->phi, grid->r0, 
					       grid->z1);
  } else {
    RZ(grid->phi, grid->r0 - 1, grid->z0 - 1) =  
      RZ(grid->phi, grid->r0, grid->z0 - 1)
      + RZ(grid->phi, grid->r0 - 1, grid->z0)
      - 0.5 * (RZ(grid->phi, grid->r0 + 1, grid->z0 - 1) 
	       + RZ(grid->phi, grid->r0 - 1, grid->z0 + 1));

    RZ(grid->phi, grid->r0 - 1, grid->z1) =  
      RZ(grid->phi, grid->r0, grid->z1)
      + RZ(grid->phi, grid->r0 - 1, grid->z1 - 1)
      - 0.5 * (RZ(grid->phi, grid->r0 + 1, grid->z1) 
	       + RZ(grid->phi, grid->r0 - 1, grid->z1 - 2));
  }
}

/** @brief  Estimates the error on this grid by comparing with the calculations
 *  for his parent and performing a Richardson extrapolation.
 */
void
pois_set_error (pois_grid_t *grid)
{
  pois_grid_t *parent;
  interpol_t *interpol;
  int pr, pz;

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

  parent = grid->parent;
  assert (NULL != parent);

  interpol = interpol_new_a (2.0, 2.0, &interpol_wackers);

  iter_grid_parent (grid, pr, pz) {
    int i, j;

    interpol_set_stencil_at ((grid_t *) grid, interpol, 
			     pr * 2.0 + 1, pz * 2.0 + 1,
			     parent->phi, pr, pz, 0);
    for (i = 0; i < 2; i++)
      for (j = 0; j < 2; j++) {
	int cr, cz;
	double v;

	cr = (pr << 1) + i;
	cz = (pz << 1) + j;
	
	v = interpol_apply (interpol, (double) cr + 0.5, (double) cz + 0.5);

	/* The error will only be used as its abs value.  Therefore
	   we calculate it here and forever. */
	RZ (grid->error, cr, cz) = fabs (RZ (grid->phi, cr, cz) - v);
	
      }
  }
  interpol_free (interpol);
}

/** @brief The refinement routine.
 *
 * Threshold is initially pois_max_error, but if we get too large refinement
 * it can be increased (after warning the user that we are not reaching
 * the accuracy that he asked for).
 */
int
pois_refine (pois_grid_t *grid, double threshold)
{
  int ir, iz;
  int zmin, zmax, rmin, rmax;
  /* Be careful: cz1 and cr1 are inclusive: when the child is actually
     created, one has to increse them by 1. */
  int cr0, cr1, cz0, cz1;
  int all_ok = TRUE;

  int tainted_line = FALSE;

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

  /* Initially we set values out of the grid, so the first comparisons
   * always change the values.
   */
  cr0 = grid->r1; cr1 = grid->r0 - 1;
  cz0 = grid->z1; cz1 = grid->z0 - 1;
 
  /* The boundary can only coincide with the parent's boundary if this
   * is an external boundary.
   */
  rmin = (grid->ext_bound & BND_MASK(BND_LEFT))? grid->r0: (grid->r0 + 1);
  rmax = (grid->ext_bound & BND_MASK(BND_RIGHT))? grid->r1: (grid->r1 - 1);
  zmin = (grid->ext_bound & BND_MASK(BND_BOTTOM))? grid->z0: (grid->z0 + 1);
  zmax = (grid->ext_bound & BND_MASK(BND_TOP))? grid->z1: (grid->z1 - 1);

  for (iz = zmin; iz < zmax && all_ok; iz++) {
    tainted_line = FALSE;
    for (ir = rmin; ir < rmax; ir++) {
      if (needs_refinement (grid, ir, iz, threshold)) {
	tainted_line = TRUE;
	if (ir < cr0) cr0 = ir;
	if (ir > cr1) cr1 = ir;
      }
    }
    if (tainted_line) {
      if (iz < cz0) cz0 = iz;
      if (iz > cz1) cz1 = iz;
    } else if (cr1 >= cr0 && cz1 >= cz0) {
      all_ok = all_ok && refine_in (grid, cr0, cz0, cr1, cz1);
      cr0 = grid->r1; cr1 = grid->r0 - 1;
      cz0 = grid->z1; cz1 = grid->z0 - 1;
    }
  }

  if (cr1 >= cr0 && cz1 >= cz0) {
    all_ok = all_ok && refine_in (grid, cr0, cz0, cr1, cz1);
  }

  /* If there was some error with some refinement, we delete the created
   * sub-grids and go back to the caller, which has to increase the refinement
   * threshold.
   */
  if (!all_ok) {
    pois_grid_t *leaf;
    free_childs (grid, leaf, pois_free_r);
    set_childless (grid);
  }
  return all_ok;
}

#define REFINE_IN_MAX_WARN_CNT 20;

/** @brief Refines and creates a new grid according to the coordinates given.
 *
 * (if the FISHPACK limit allows it). Returns TRUE if the refinement
 * succeeded, FALSE otherwise (if the FISHPACK limit was exceeded).
*/
static int
refine_in (pois_grid_t *grid, int cr0, int cz0, int cr1, int cz1)
{
  pois_grid_t *child;
  int nr0, nr1, nz0, nz1;
  static int warn_cnt = REFINE_IN_MAX_WARN_CNT;

  nr0 = cr0 << 1;
  nz0 =  cz0 << 1;
  nr1 = (++cr1) << 1;
  nz1 = (++cz1) << 1;

  if ((nr1 - nr0) > FISH_MAX_GRIDPOINTS ||  
      (nz1 - nz0) > FISH_MAX_GRIDPOINTS) {
    if (warn_cnt > 0) {
      warning ("FISHPACK limit exceeded.  Not refining grid "
	       grid_printf_str ".\n", nr0, nz0, nr1, nz1, grid->level + 1);
      warn_cnt --;
    }
    if (warn_cnt == 0) {
      warning ("No more warnings about FISHPACK limit.\n");
      warn_cnt --;
    }
    return FALSE;
  }

  child = pois_new_a (nr0, nz0, nr1, nz1);

  add_child (grid, child);
  grid_inherit_ext_bound ((grid_t *) child);
  debug (3, "new child created {r0 = %d, z0 = %d, r1 = %d, z1 = %d, "
	 "level = %d}\n", child->r0, child->z0, child->r1, child->z1, 
	 child->level);

  return TRUE;
}

/** @brief Calculates the error measure of this grid.
 *
 * For debug/testing purposes
 */
void
pois_error_measures (pois_grid_t *grid, double *L1, double *L2, double *Lmax)
{
  int ir, iz, i;
  double err;
  *L2 = 0.0;
  *L1 = 0.0;
  *Lmax = 0.0;
  i = 0;

  iter_grid (grid, ir, iz) {
    i++;
    err = RZ (grid->error, ir, iz);
    if (err > *Lmax) *Lmax = err;
    *L1 += err;
    *L2 += err * err;
  }

  *L1 /= i;
  *L2 = sqrt (*L2 / i);

}

/** @brief Writes some error measures of this grid and its descendants
 *  into @a fp
 */
void
pois_write_error_r (pois_grid_t *grid, FILE *fp)
{
  double L1, L2, Lmax;
  pois_grid_t *child;

  pois_error_measures (grid, &L1, &L2, &Lmax);

  fprintf (fp, "%d %g %g %g %g %g\n", grid->level, dr[grid->level],
	   dz[grid->level], L1, L2, Lmax);

  iter_childs (grid, child) {
    pois_write_error_r (child, fp);
  }
}

/** @brief Takes a @a cdr tree and solves the Poisson equation for it.
 *
 * Returns a poisson tree.
 */
pois_grid_t **
pois_solve_a (cdr_grid_t *cdr, pois_problem_t *prob)
{
  /* Call to the generic Poisson/Helmholtz solver */
  return pois_gen_solve_a (cdr, prob, e_mappers, 0.0);
}

/** @brief Solves a Poisson/Helmholtz equation with FISHPACK and maps
 *  the result using a given set of mapper.
 *
 *  From this method downwards (in the sense of calling order, not file
 *  organization), all the code is shared by the Poisson solver and the
 *  Photoionization (Helmholtz) solver.
 */
pois_grid_t **
pois_gen_solve_a (cdr_grid_t *cdr, pois_problem_t *prob,
		  mapper_t **mappers, double es)
{
  cdr_grid_t *tree;
  pois_grid_t **pois_modes;
  int mode;

  pois_modes = (pois_grid_t**) xmalloc (sizeof(pois_grid_t*) * max_ntheta);
  assert (0 == cdr->level);
  tree = cdr_add_coarser_grids_a (cdr, prob->extra_levels);

  /* This is the high-level parallelizable loop! */
#pragma omp parallel
  {
#pragma omp for schedule(dynamic)
    for (mode = 0; mode < max_ntheta; mode++) {
      pois_modes[mode] = pois_solve_mode (tree, cdr, prob, mode, es);

      if (pois_inhom && es == 0.0 && mode == 0) {
	/* Add the inhomogeneous laplacian potential to phi.
	 * It has only to be added to the zero mode.
	 */
	double q_factor;
	q_factor = pois_inhom_q_factor (pois_modes[0]);
	
	/* Was:
	   pois_add_inhom_phi_r (pois_modes[0], q_factor);
	*/
	pois_add_inhom_phi_r (pois_modes[0], -q_factor);
      }

      /* For the interpolation, we require an extra buffer cell
       * in some boundaries and we check for overflows inside the
       * mapper methods.
       */
      map_trees_r (mappers, (grid_t*) pois_modes[mode], (grid_t*) cdr, 
		   mode, FALSE, TRUE, FALSE, 
		   /*s_buf = */ 2, 
		   /*t_buf = */ 4);
      map_trees_r (mappers, (grid_t*) pois_modes[mode], (grid_t*) cdr, 
		   mode, TRUE, FALSE, TRUE, 
		   /*s_buf = */ 1, 
		   /*t_buf = */ 2);
    }
  }
  cdr_free_coarser_grids (tree, extra_pois_levels);

  return pois_modes;
}

/** @brief Solves a single Fourier mode.
 *
 * @a tree is the root of the @a cdr tree with extra_pois_levels added.
 * @a cdr is the cdr grid at level 0, contained in tree:
 *
 *   cdr = tree->first_child->firts_child...->first_child
 *            \_________ extra_pois_levels __________/
 */
pois_grid_t *
pois_solve_mode (cdr_grid_t *tree, cdr_grid_t *cdr, pois_problem_t *prob,
		 int mode, double es)
{
  pois_grid_t *pois;

  pois = pois_init_tree_a (0, 0, pois_gridpoints_r, pois_gridpoints_z);

  /* Solves the coarsest Poisson grid. */
  /* Changed Wed Mar 19 13:50:23 2008: was s_buf = 1. */
  map_trees_r (charge_mappers, (grid_t *) tree, (grid_t *) pois, 
	       mode, TRUE, TRUE, FALSE, 0, 2);
  pois_solve_grid (pois, prob, -w2k[mode], es);

  /* Wed May 14 11:00:30 2008:
   * I divide pois_max_error by abs(wk[mode]) because in the calculation
   * of e_theta, the errors get multiplied by wk[mode].  Hence if we
   * want to get similar errors in all modes, we have to remove it from the
   * threshold.  The 1 is because wk[0] = 0, but we do not want to have
   * a \inf threshold.
   */
  pois_solve_r (pois->first_child, tree, prob, mode, es, 
		prob->max_error / (1 + abs(wk[mode])));

  return pois;
}

#define POIS_SOLVE_R_MAX_WARN 50

/** @brief  Recursively solves one Poisson grid. */
void
pois_solve_r (pois_grid_t *pois, cdr_grid_t *cdr, pois_problem_t *prob,
	      int mode, double es, double threshold)
{
  pois_grid_t *child;
  int succeeded = FALSE;
  static int warn_cnt = POIS_SOLVE_R_MAX_WARN;

  /* Changed Wed Mar 19 13:51:00 2008: was s_buf = 1 */
  map_grid_r (charge_mappers, (grid_t *) cdr, (grid_t *) pois, 
	      mode, TRUE, TRUE, FALSE, 0, 2);
  debug(3, "w2k[mode = %d] = %f\n", mode, w2k[mode]);

  pois_solve_grid (pois, prob, -w2k[mode], es);

  if (pois->level >= prob->max_level) return;

  pois_set_error (pois);

  while (!succeeded) {
    succeeded = pois_refine (pois, threshold);

    if (!succeeded) {
      threshold *= 2;
      if (warn_cnt > 0) {
	warning ("Poisson error threshold was too small: trying %g.\
 Consider increasing pois_max_error\n",
		 threshold);
	if (warn_cnt >= 0) warn_cnt --;
      }
    }
  }
  iter_childs (pois, child) {
    pois_solve_r (child, cdr, prob, mode, es, threshold);
  }
}

/** @brief Returns 1 if the given cell has to be refined.
 *
 * which means that it satisfies at least one of the following:
 * a. The error in the cell is larger than the threshold.
 * b. The cell has a neightbour with an error larger than the threshold
 */
static int
needs_refinement (pois_grid_t *grid, int ir, int iz, double threshold)
{
  int i, j;
  int r = FALSE;
  int irmin, irmax, izmin, izmax;

  irmin = ir > grid->r0 + 3? ir - 1: grid->r0;
  irmax = ir < grid->r1 - 4? ir + 1: grid->r1 - 1;
  izmin = iz > grid->z0 + 3? iz - 1: grid->z0;
  izmax = iz < grid->z1 - 4? iz + 1: grid->z1 - 1;

  for (i = irmin; i <= irmax && !r; i++)
    for (j = izmin; j <= izmax; j++)
      if (RZ(grid->error, i, j) > threshold) {
	r = TRUE;
	break;
      }

  return r;
}

/*****************************************************************************
 * Functions for the implemetation of a needle-plate discharge.              *
 * The method used here is that of the "floating charge":                    *
 * A charge is located at pois_inhom_z, which is used to keep the potential  *
 * at L_z constant.                                                          *
 ****************************************************************************/

/** @brief pois_inhom_init QQQQ */
void
pois_inhom_init (void)
{
  /* The domain of solution of the Poisson equation will be different
   * from that of the CDR equation.  The difference between them is
   * given by the parameter needle_length, but has to be rounded
   * to a multiple of the grid size at level 0.
   */
  pois_gridpoints_z = (gridpoints_z
		       + (int) nearbyint (needle_length / dz[0]));
  pois_inhom_L = pois_gridpoints_z * dz[0];
		       
  pois_inhom_z = L_z + needle_radius;
  pois_inhom_dphi = pois_inhom_phi (0, pois_inhom_L) - pois_inhom_phi (0, L_z);

}

/* We use the following functions to calculate the fields and potential
 * created by a unit charge located at (r = 0, z = b).
 * These are created by the method of mirror charges, with a total of
 * pois_inhom_reflections charges.
 */

/* Single term of phi. */
static double
inhom_phi_term (double r, double z, double b)
{
  double insqrt;

  insqrt = (SQ(z - b) + SQ(r));

  return invfourpi / sqrt (insqrt);
}

/** @brief Common decay factor for the electric fields. */
static double
inhom_e_term (double r, double z, double b)
{
  double insqrt;

  insqrt = (SQ(z - b) + SQ(r));
  insqrt = insqrt * insqrt * insqrt;

  return invfourpi / sqrt (insqrt);
}

/** @brief Single term of E_r. */
static double
inhom_er_term (double r, double z, double b)
{
  return r * inhom_e_term (r, z, b);
}

/** @brief Single term of E_z */
static double
inhom_ez_term (double r, double z, double b)
{
  return (z - b) * inhom_e_term (r, z, b);
}

/* Since we have to sum over the same locations of the mirror charges
 * but with different "field functions", we use the following macro.
 *
 * Note that using a single function for inhom_phi, inhom_er, and inhom_ez
 * would not be a good idea since they are evaluated at different points.
 * Passing that function a pointer to another function would have been somewhat
 * slower and, I think, not so readable.
 */
#define sum_reflections(total_, field_, r_, z_)					\
       do {									\
	 int i__;								\
	 total_ = (field_ (r_, z_, pois_inhom_z) -				\
		   field_ (r_, z_, -pois_inhom_z));				\
										\
	 for (i__ = 2; i__ <= pois_inhom_reflections; i__ += 2) {		\
	   total_ += (field_ (r_, z_, pois_inhom_z + i__ * pois_inhom_L) +	\
		      field_ (r_, z_, pois_inhom_z - i__ * pois_inhom_L));	\
	 									\
	   total_ -= (field_ (r_, z_, -pois_inhom_z + i__ * pois_inhom_L) +	\
		      field_ (r_, z_, -pois_inhom_z - i__ * pois_inhom_L));	\
	 }								\
       } while (0)

/** @brief Returns the potential created by a unit charge located at pois_inhom_z
 *  between an electrode at 0 and a second one at pois_inhom_L.
 */
double
pois_inhom_phi (double r, double z)
{
  double res;

  sum_reflections (res, inhom_phi_term, r, z);

  return res;
}

/** @brief Radial component of the field. */
double
pois_inhom_er (double r, double z)
{
  double res;

  sum_reflections (res, inhom_er_term, r, z);

  return res;
}

/** @brief Axial component of the field. */
double
pois_inhom_ez (double r, double z)
{
  double res;

  sum_reflections (res, inhom_ez_term, r, z);

  return res;
}

/** @brief Once we have the electrostatic potential created by the space charges,
 * we use this function to compute the multiplying factor of the floating
 * charge.
 *
 * One can also forget about keeping the potential fixed somewhere and impose
 * a constant (non-floating, but may depend on time) charge pois_inhom_fixed_q,
 * which is used if nonzero. The use of that is to simulate a charged cloud
 * close to the earth, which acts as an electrode: usually you would also combine
 * that with the sprite module, that allows varing neutral densities.
 */
double
pois_inhom_q_factor (pois_grid_t *pois)
{
  double u;
  /* Note that the signs preceding the numerical potentials are inverted since
   * in our notation, we use the opposite of the physical one, in order
   * to use the charge as the source of the Poisson equation.
   */

  if (pois_inhom_fixed_q != 0.0) {
    return pois_inhom_fixed_q_t;
  }

  /* Was: 
  u = - E_z * (pois_inhom_L - L_z) + pois_phi_at (pois, 0.0, L_z, 0.0);
  */

  u = E_z * (L_z - pois_inhom_L) + pois_phi_at (pois, 0.0, L_z, 0.0);

  return  -u / pois_inhom_dphi;
}

/** @brief Adds the potential created by the needle to a Poisson grid and their
 *  descendants.
 *
 * Note that \f$\phi\f$ here is in Fourier space and the Laplacian
 * potential is always axi-symmetric, so you only have to call this function
 * with the zero-mode of the Poisson grids.
 *
 * The reason that it is better to add the inhomogeneous field to the potential
 * and not later to the electric fields is that if we do that, we will
 * have to substract two large and very curved functions that give something
 * close to zero: the error we make then will be comparable to the result.
 */
void
pois_add_inhom_phi_r (pois_grid_t *grid, double q)
{
  int ir, iz;
  pois_grid_t *child;
  double res;

  debug (2, "pois_add_inhom_phi_r (" grid_printf_str ", q = %f)\n",
	 grid_printf_args(grid), q);

  /* Note that there is no factor grid->ntheta coming from thr FFT 
     un-normalization.  The reason is that actually q is sqrt(N) smaller
     than it has to be and later, in the inverse FFT we will multiply by
     sqrt(N).  So everything is correct, I believe. */
  iter_grid_n (grid, ir, iz, 2) {

    sum_reflections (res, inhom_phi_term, r_at(ir, grid->level), z_at(iz, grid->level));


    RZ(grid->phi, ir, iz) += q * res;
  }

  iter_childs (grid, child) {
    pois_add_inhom_phi_r (child, q);
  }
}

/***************************************************************
 *  Mapper functions for the components of the electric field. *
 *  See mapper.h for an explanation of each of these methods.  *
 ***************************************************************/

/** @brief er_copy QQQQ */
void
er_copy (mapper_t *mapper, grid_t *source, grid_t *target, 
	 int ir, int iz, int itheta)
{
  cdr_grid_t *cdr;
  pois_grid_t *pois;

  cdr = (cdr_grid_t*) target;
  pois = (pois_grid_t*) source;

  if (ir < pois->r1)
    RZT (cdr->er, ir, iz, itheta) = ER_RZ (pois, ir, iz);
}

/** @brief ez_copy QQQQ */
void
ez_copy (mapper_t *mapper, grid_t *source, grid_t *target, 
	 int ir, int iz, int itheta)
{
  cdr_grid_t *cdr;
  pois_grid_t *pois;

  cdr = (cdr_grid_t*) target;
  pois = (pois_grid_t*) source;

  if (iz < pois->z1)
    RZT (cdr->ez, ir, iz, itheta) = EZ_RZ (pois, ir, iz);
}

/** @brief etheta_copy QQQQ */
void
etheta_copy (mapper_t *mapper, grid_t *source, grid_t *target, 
	     int ir, int iz, int itheta)
{
  cdr_grid_t *cdr;
  pois_grid_t *pois;

  cdr = (cdr_grid_t*) target;
  pois = (pois_grid_t*) source;

  RZT (cdr->etheta, ir, iz, itheta) = RZ (pois->phi, ir, iz);
}

/** @brief er_coarsen QQQQ */
void 
er_coarsen (mapper_t *mapper, grid_t *source, grid_t *target, 
	    int ir, int iz, int itheta)
{
  cdr_grid_t *cdr;
  pois_grid_t *pois;
  int level_diff, er_ir, er_iz;

  cdr = (cdr_grid_t*) target;
  pois = (pois_grid_t*) source;

  level_diff = pois->level - cdr->level;

  er_iz = (iz << level_diff) + (1 << (level_diff - 1));
  er_ir = ((ir + 1) << level_diff) - 1;

  if (grid_contains (source, er_ir, er_iz - 1, GRID_INSIDE) && 
      grid_contains (source, er_ir, er_iz, GRID_INSIDE)) {

    RZT (cdr->er, ir, iz, itheta) = 0.5 * (ER_RZ(pois, er_ir, er_iz) +
					   ER_RZ(pois, er_ir, er_iz - 1));
  }
}

/** @brief ez_coarsen QQQQ */
void 
ez_coarsen (mapper_t *mapper, grid_t *source, grid_t *target, 
	    int ir, int iz, int itheta)
{
  cdr_grid_t *cdr;
  pois_grid_t *pois;
  int level_diff, ez_ir, ez_iz;

  cdr = (cdr_grid_t*) target;
  pois = (pois_grid_t*) source;

  level_diff = pois->level - cdr->level;

  ez_iz = ((iz + 1) << level_diff) - 1;
  ez_ir = (ir << level_diff) + (1 << (level_diff - 1));

  if (grid_contains (source, ez_ir - 1, ez_iz, GRID_INSIDE) && 
      grid_contains (source, ez_ir, ez_iz, GRID_INSIDE)) {

    RZT(cdr->ez, ir, iz, itheta) = 0.5 * (EZ_RZ(pois, ez_ir, ez_iz) +
					  EZ_RZ(pois, ez_ir - 1, ez_iz));
  }
}

/** @brief etheta_coarsen QQQQ */
void 
etheta_coarsen (mapper_t *mapper, grid_t *source, grid_t *target, 
		int ir, int iz, int itheta)
{
  cdr_grid_t *cdr;
  pois_grid_t *pois;
  int level_diff, ez_ir, er_iz;

  cdr = (cdr_grid_t*) target;
  pois = (pois_grid_t*) source;

  level_diff = pois->level - cdr->level;

  er_iz = (iz << level_diff) + (1 << (level_diff - 1));
  ez_ir = (ir << level_diff) + (1 << (level_diff - 1));

  if (grid_contains (source, ez_ir, er_iz, GRID_INSIDE) && 
      grid_contains (source, ez_ir - 1, er_iz, GRID_INSIDE) &&
      grid_contains (source, ez_ir, er_iz - 1, GRID_INSIDE) &&
      grid_contains (source, ez_ir - 1, er_iz - 1, GRID_INSIDE)) {

    RZT (cdr->etheta, ir, iz, itheta)  = 0.25 *
      (RZ (pois->phi, ez_ir, er_iz) 
       + RZ (pois->phi, ez_ir - 1, er_iz) 
       + RZ (pois->phi, ez_ir, er_iz - 1) 
       + RZ (pois->phi, ez_ir - 1, er_iz - 1));
  }
}

/** @brief er_interpol QQQQ */
int 
er_interpol_set (mapper_t *mapper, grid_t *source, interpol_t *interpol,
		      int pr, int pz, int itheta)
{
  pois_grid_t *pois;

  pois = (pois_grid_t*) source;

  /* When we set Neumann b.c. in one boundary and Dirichlet in the opposite
   * a larger-than-usual error appears in the boundaries because we are
   * matching potetntials in one side and fields in the other.  The only way
   * I see of avoiding this is to forget about the extrapolated values of
   * the grids.  But we still need those values on the external boundaries.
   */
  if (pr < pois->r0 || pz < pois->z0
      || pr > pois->r1 - 1 || pz > pois->z1)
    return FALSE;

  interpol_set_stencil (interpol, 
			er_r_at (pr - 1, pois->level), 
			er_z_at (pz - 1, pois->level),
			ER_RZ (pois, pr - 1, pz - 1), 
			ER_RZ (pois, pr - 1, pz),
			ER_RZ (pois, pr, pz - 1), 
			ER_RZ (pois, pr, pz));

  return TRUE;
}

/** @brief ez_interpol_set QQQQ */
int
ez_interpol_set (mapper_t *mapper, grid_t *source, interpol_t *interpol,
		      int pr, int pz, int itheta)
{
  pois_grid_t *pois;

  pois = (pois_grid_t*) source;

  if (pr < pois->r0 || pz < pois->z0
      || pr > pois->r1 || pz > pois->z1 - 1)
    return FALSE;

  interpol_set_stencil (interpol, 
			ez_r_at (pr - 1, pois->level), 
			ez_z_at (pz - 1, pois->level),
			EZ_RZ (pois, pr - 1, pz - 1), 
			EZ_RZ (pois, pr - 1, pz),
			EZ_RZ (pois, pr, pz - 1), 
			EZ_RZ (pois, pr, pz));
  return TRUE;
}

/** @brief etheta_interpol_set QQQQ */
int 
etheta_interpol_set (mapper_t *mapper, grid_t *source, interpol_t *interpol,
		      int pr, int pz, int itheta)
{
  pois_grid_t *pois;

  pois = (pois_grid_t*) source;

  if (pr < pois->r0 || pz < pois->z0
      || pr > pois->r1 - 1|| pz > pois->z1 - 1)
    return FALSE;

  interpol_set_stencil_at (source, interpol, 
			   r_at (pr, pois->level),
			   z_at (pz, pois->level),
			   pois->phi, pr, pz, 0);
  return TRUE;
}

/** @brief er_interpol QQQQ */
void
er_interpol (mapper_t *mapper, grid_t *source, grid_t *target, 
	     interpol_t *interpol, int ir, int iz, int itheta)
{
  double er_r, er_z;
  cdr_grid_t *cdr;
  pois_grid_t *pois;
  int level_diff;

  cdr = (cdr_grid_t *) target;
  pois = (pois_grid_t *) source;

  level_diff = cdr->level - pois->level;

  /* Note that when pr = pois->r1, the electric field is calculated
   * from the extrapolation of phi and it is hence less accurate.
   * The same applies below to pz and pois->z1.
   */
  if ((ir >> level_diff) >= pois->r1) return;

  er_r = er_r_at (ir, cdr->level);
  er_z = er_z_at (iz, cdr->level);

  RZT (cdr->er, ir, iz, itheta) = interpol_apply (interpol, er_r, er_z);
}

/** @brief ez_interpol QQQQ */
void
ez_interpol (mapper_t *mapper, grid_t *source, grid_t *target, 
	     interpol_t *interpol, int ir, int iz, int itheta)
{
  double ez_r, ez_z;
  cdr_grid_t *cdr;
  pois_grid_t *pois;
  int level_diff;

  cdr = (cdr_grid_t *) target;
  pois = (pois_grid_t *) source;

  level_diff = cdr->level - pois->level;

  if ((iz >> level_diff) >= pois->z1) return;

  ez_r = ez_r_at (ir, cdr->level);
  ez_z = ez_z_at (iz, cdr->level);

  RZT (cdr->ez, ir, iz, itheta) = interpol_apply (interpol, ez_r, ez_z);
}

/** @brief etheta_interpol QQQQ */
void
etheta_interpol (mapper_t *mapper, grid_t *source, grid_t *target, 
		 interpol_t *interpol, int ir, int iz, int itheta)
{
  double eth_r, eth_z;
  cdr_grid_t *cdr;

  cdr = (cdr_grid_t *) target;

  eth_r = r_at (ir, cdr->level);
  eth_z = z_at (iz, cdr->level);

  RZT (cdr->etheta, ir, iz, itheta) = 
    interpol_apply (interpol, eth_r, eth_z);
}

/** @brief Mapping of the charge. */
void
charge_copy (mapper_t *mapper, grid_t *source, grid_t *target, 
	     int ir, int iz, int itheta)
{
  cdr_grid_t *cdr;
  pois_grid_t *pois;

  pois = (pois_grid_t*) target;
  cdr = (cdr_grid_t*) source;

  RZ (pois->charge, ir, iz) = RZT (cdr->charge, ir, iz, itheta);
}

/** @brief charge_interpol_set QQQQ. */
int
charge_interpol_set (mapper_t *mapper, grid_t *source, interpol_t *interpol,
		     int pr, int pz, int itheta)
{
  cdr_grid_t *cdr;

  cdr = (cdr_grid_t*) source;

  if (pr < cdr->r0 || pz < cdr->z0
      || pr > cdr->r1 - 1|| pz > cdr->z1 - 1)
    return FALSE;

  interpol_set_stencil_at (source, interpol, 
			   r_at (pr, cdr->level),
			   z_at (pz, cdr->level),
			   cdr->charge, pr, pz, itheta);
  return TRUE;
}

/** @brief charge_interpol QQQQ. */
void
charge_interpol (mapper_t *mapper, grid_t *source, grid_t *target, 
		 interpol_t *interpol, 
		 int ir, int iz, int itheta)
{
  double r, z;
  pois_grid_t *pois;

  pois = (pois_grid_t *) target;

  r = r_at (ir, pois->level);
  z = z_at (iz, pois->level);

  RZT (pois->charge, ir, iz, 0) = interpol_apply (interpol, r, z);
}

/** @brief Finds the best-possible approximation for the potential at a given
 * point (@a r, @a z, @a theta) by interpolating from the finest grid that covers
 * that point.
 */
double
pois_phi_at (pois_grid_t *grid, double r, double z, double theta)
{
  int ir, iz, itheta, it, itn;
  REAL phi_it[2];
  interpol_t *phi_at_interpol = NULL;

  debug (3, "pois_phi_at (" grid_printf_str ", r = %g, z = %g, theta = %g)\n",
	 grid_printf_args (grid), r, z, theta);

  grid = (pois_grid_t *) grid_finest_containing_r ((grid_t*) grid, r, z);

  if (NULL == grid) {
    fatal ("In pois_phi_at: potential outside the grid");
  }

  phi_at_interpol = interpol_new_a (dr[grid->level], dz[grid->level],
				    &interpol_quadratic);

  ir = (int) (r / dr[grid->level]);
  iz = (int) (z / dz[grid->level]);

  /* If r == 0, it doesn't matter which theta we use, so we use 0.
   * also if we are in a 2D simulation. 
   */
  if (grid->ntheta > 1 && r > 0.0) {
    itheta = (int) (theta / dtheta);
    itn = 2;
  } else {
    itheta = 0;
    itn = 1;
  }

  for (it = 0; it < itn; it++) {
    interpol_set_stencil_at ((grid_t *) grid, phi_at_interpol,
			     r_at (ir, grid->level),
			     z_at (iz, grid->level), 
			     grid->phi, ir, iz, itheta + it);

    phi_it[it] = interpol_apply (phi_at_interpol, r, z);
  }

  interpol_free (phi_at_interpol);

  if (grid->ntheta > 1 && r > 0.0) {
    return (phi_it[0] * (theta - theta_at (itheta)) +
	    phi_it[1] * (theta_at (itheta + 1) - theta)) / dtheta;
  } else {
    return phi_it[0];
  }
}

/** @brief Dumps all the contents of the given grid into filenames given by
 * prefix and name
 */
void
pois_dump (pois_grid_t *grid, const char *prefix, const char *name)
{
  char *fname;
  int m = pois_output_margin;

  asprintf (&fname, "%s/r.%s.tsv", prefix, name);
  rz_axis_dump (fname, grid->r0 - m, grid->r1 + m, dr[grid->level]);
  free (fname);

  asprintf (&fname, "%s/z.%s.tsv", prefix, name);
  rz_axis_dump (fname, grid->z0 - m, grid->z1 + m, dz[grid->level]);
  free (fname);

  asprintf (&fname, "%s/phi.%s.tsv", prefix, name);
  rz_dump (grid->phi, fname, "w", grid->r0 - m, grid->z0 - m, 
	   grid->r1 + m, grid->z1 + m);
  free (fname);

  asprintf (&fname, "%s/charge.%s.tsv", prefix, name);
  rz_dump (grid->charge, fname, "w", grid->r0 - m, grid->z0 - m, 
	   grid->r1 + m, grid->z1 + m);
  free (fname);

  asprintf (&fname, "%s/error.%s.tsv", prefix, name);
  rz_dump (grid->error, fname, "w", grid->r0 - m, grid->z0 - m, 
	   grid->r1 + m, grid->z1 + m);
  free (fname);
}

/** @brief Dumps all the contents of the given grid into filenames given by
 * prefix and name
 */
void
pois_dump_r (pois_grid_t *grid, const char *prefix, const char *name)
{
  pois_grid_t *child;
  char *cname;
  int i = 0, nchilds;
  char codes[] = "abcdefghijklmnopqrstuvwxyz";

  pois_dump (grid, prefix, name);

  nchilds = grid_howmany_children ((grid_t *) grid);

  for (i = 0; i < nchilds; i++) {
    child = (pois_grid_t *) grid_get_child ((grid_t*) grid, nchilds - i - 1);

    assert (NULL != child);

    if (i > sizeof(codes) - 2)
      asprintf (&cname, "%s{%d}", name, i);
    else
      asprintf (&cname, "%s%c", name, codes[i]);

    pois_dump_r (child, prefix, cname);
    free (cname);
  }
}