diff --git a/src/photo.c b/src/photo.c
new file mode 100644
index 0000000000000000000000000000000000000000..ed3efc613f0fe041effffda9b22fb3f5d3c42c17
--- /dev/null
+++ b/src/photo.c
@@ -0,0 +1,400 @@
+/** @file photo.c
+ *  @brief Routines for the calculation of the photoionization source term.
+ *
+ * We use the Helmholtz approximation for the kernel in the photoionization
+ * integral and we use the routines of the Poisson solver for the Helmholtz
+ * equation (FISHPACK was modified to take a new inhomogeneous term into
+ * account).
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "cdr.h"
+#include "interpol2.h"
+#include "mapper.h"
+#include "parameters.h"
+#include "photo.h"
+#include "poisson.h"
+#include "proto.h"
+#include "rz_array.h"
+#include "species.h"
+
+void photo_copy (mapper_t *mapper, grid_t *source, grid_t *target, 
+		 int ir, int iz, int itheta);
+void photo_coarsen (mapper_t *mapper, grid_t *source, grid_t *target, 
+		    int ir, int iz, int itheta);
+int photo_interpol_set (mapper_t *mapper, grid_t *source, 
+			 interpol_t *interpol,
+			 int pr, int pz, int itheta);
+void photo_interpol (mapper_t *mapper, grid_t *source, grid_t *target, 
+		     interpol_t *interpol, 
+		     int ir, int iz, int itheta);
+
+photo_term_t *photo_terms = NULL;
+
+mapper_t photo_mapper = {&interpol_bilin, photo_coarsen, photo_copy,
+			 photo_interpol_set, photo_interpol};
+
+mapper_t *photo_mappers[] = {&photo_mapper, NULL};
+
+pois_problem_t *pois_photo_1;
+pois_problem_t *pois_photo_2;
+extern pois_problem_t *pois_electrostatic;
+
+/** @brief Initializes photoionization QQQQ? */
+void
+photo_init ()
+{
+  /* Since it makes no sense, we use photo_extra_levels < 0 to say that
+     the conditions for the photoionization are the same as for
+     the electrostatic problem.  This is needed for backwards compatibility. */
+  if (extra_photo_levels < 0) 
+  {
+    pois_photo_1 = pois_electrostatic;
+    pois_photo_2 = pois_electrostatic;
+    return;
+  }
+
+  /* Two different sets of refinement criteria for the two photoionization
+     terms. If extra_photo_levels_2 < 0, then pois_photo_2 = pois_photo_1 */
+  pois_photo_1 = xmalloc (sizeof (pois_problem_t));
+
+  pois_photo_1->max_level = photo_max_level;
+  pois_photo_1->extra_levels = extra_photo_levels;
+  pois_photo_1->max_error = photo_max_error;
+
+  pois_photo_1->bnd_right = photo_bnd_right;
+  pois_photo_1->bnd_top = photo_bnd_top;
+  pois_photo_1->bnd_bottom = photo_bnd_bottom;
+
+  if (extra_photo_levels_2 < 0)
+  {
+     pois_photo_2 = pois_photo_1;
+     return;
+  }
+
+  pois_photo_2 = xmalloc (sizeof (pois_problem_t));
+
+  pois_photo_2->max_level = photo_max_level_2;
+  pois_photo_2->extra_levels = extra_photo_levels_2;
+  pois_photo_2->max_error = photo_max_error_2;
+
+  pois_photo_2->bnd_right = photo_bnd_right_2;
+  pois_photo_2->bnd_top = photo_bnd_top_2;
+  pois_photo_2->bnd_bottom = photo_bnd_bottom_2;
+}
+
+
+/** @brief Copies the photoionization QQQQ? */
+void
+photo_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->photo, ir, iz, itheta) = RZ (pois->phi, ir, iz);
+}
+
+/** @brief Coarsens the photoionization QQQQ? */
+void 
+photo_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, z, r;
+
+  cdr = (cdr_grid_t*) target;
+  pois = (pois_grid_t*) source;
+
+  level_diff = pois->level - cdr->level;
+
+  z = (iz << level_diff) + (1 << (level_diff - 1));
+  r = (ir << level_diff) + (1 << (level_diff - 1));
+
+  if (grid_contains (source, r, z, GRID_INSIDE) && 
+      grid_contains (source, r - 1, z, GRID_INSIDE) &&
+      grid_contains (source, r, z - 1, GRID_INSIDE) &&
+      grid_contains (source, r - 1, z - 1, GRID_INSIDE)) {
+    
+    RZT (cdr->photo, ir, iz, itheta)  = 0.25 *
+      (RZ (pois->phi, r, z) 
+       + RZ (pois->phi, r - 1, z) 
+       + RZ (pois->phi, r, z - 1) 
+       + RZ (pois->phi, r - 1, z - 1));
+  }
+
+}
+
+/** @brief photo_interpol_set QQQQ? */
+int 
+photo_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;
+
+  interpol_set_stencil (interpol,
+			r_at (pr, pois->level),
+			z_at (pz, pois->level),
+			RZ (pois->phi, pr, pz), 
+			RZ (pois->phi, pr, pz + 1),
+			RZ (pois->phi, pr + 1, pz), 
+			RZ (pois->phi, pr + 1, pz + 1));
+
+  return TRUE;
+}
+
+/** @brief photo_interpol QQQQ? */
+void
+photo_interpol (mapper_t *mapper, grid_t *source, grid_t *target, 
+		interpol_t *interpol, int ir, int iz, int itheta)
+{
+  double r, z;
+  cdr_grid_t *cdr;
+
+  cdr = (cdr_grid_t *) target;
+
+  r = r_at (ir, cdr->level);
+  z = z_at (iz, cdr->level);
+
+  RZT (cdr->photo, ir, iz, itheta) = interpol_apply (interpol, r, z);
+}
+
+
+/** @brief Registers a photoionization term with given @a A and @a lambda. */
+void
+photo_register (double A, double lambda)
+{
+  photo_term_t *t;
+
+  t = (photo_term_t *) xmalloc (sizeof (photo_term_t));
+
+  t->A = A;
+  t->lambda = lambda;
+  t->next = photo_terms;
+
+  photo_terms = t;
+}
+
+/** @brief Copies a list of photo terms into @a *dest. */
+void
+photo_copy_list (photo_term_t *src, photo_term_t **dest)
+{
+  photo_term_t *ptr_term;
+  photo_term_t *p = NULL, *p1;
+
+  *dest = NULL;
+
+  while (src) {
+    p1 = (photo_term_t *) xmalloc (sizeof (photo_term_t));
+    if (p) {
+      p->next = p1;
+    } else {
+      *dest = p1;
+    }
+
+    p = p1;
+    p->A = src->A;
+    p->lambda = src->lambda;
+    src = src->next;
+  }
+
+  if (p) p->next = NULL;
+}
+
+
+/** @brief Unregisters all the photoionization terms and
+ *  frees the allocated space.
+ */
+void
+photo_unregister_all (void)
+{
+  photo_term_t *t, *next;
+
+  for (t = photo_terms; t; t = next) {
+    next = t->next;
+    free (t);
+  }
+
+  photo_terms = NULL;
+}
+
+/** @brief Transforms back the photoionization calculation into real space. */
+void
+photo_dft_r (cdr_grid_t *grid, int sign)
+{
+  cdr_grid_t *leaf;
+
+  debug (3, "photo_dft_r(" grid_printf_str ", %d)\n", 
+	 grid_printf_args(grid), sign);
+
+  iter_childs (grid, leaf) {
+    photo_dft_r (leaf, sign);
+  }
+
+  dft_transform (grid->photo, grid->photo, sign);
+}
+
+/** @brief Copies the derivative of the ion density into cdr->charge
+ *
+ * The ion density, which is supposed to be at this point, the impact
+ * ionization) into cdr->charge, which will be used as the source for
+ * the Poisson/Helmholtz solver. */
+void
+photo_copy_source (cdr_grid_t *grid)
+{
+  int ir, iz, itheta;
+
+   debug(3, "photo_copy_source (" grid_printf_str ")\n",
+	 grid_printf_args(grid));
+
+#pragma omp parallel
+  {
+#pragma omp for private (ir, iz)
+    iter_grid_3d_n (grid, ir, iz, itheta, 2) {
+      RZT (grid->charge, ir, iz, itheta) = 
+	RZT (grid->d_dens[ions], ir, iz, itheta) / grid->ntheta;
+    }
+  }
+}
+
+/** @brief Recursive version of @a photo_copy_source. */
+mk_recursive (photo_copy_source, cdr_grid_t)
+
+/** @brief Once a photoionization term is computed, we add it to d_dens.
+ *
+ * Note that we already copied the contents of d_dens[ions] into
+ * charge, so when we solve again the Helholtz equation, the source
+ * is still the same.
+ */
+void
+photo_add_term (photo_term_t *term, cdr_grid_t *cdr)
+{
+  int s, ir, iz, itheta;
+  int updated[2] = {electrons, photo_ions};
+  
+  debug (3, "photo_add_term (" photo_printf_str ", " grid_printf_str ")\n",
+	 photo_printf_args(term), grid_printf_args(cdr));
+    
+#pragma omp parallel
+    {
+#pragma omp for private (ir, iz, s)
+      iter_grid_3d_n (cdr, ir, iz, itheta, 2) {
+	for (s = 0; s < 2; s++) {
+	  RZT (cdr->d_dens[updated[s]], ir, iz, itheta) += 
+	    term->A * RZT (cdr->photo, ir, iz, itheta);
+	}
+      }
+    }  
+}
+
+/** @brief ...and the recursive version of @a photo_add_term */
+void
+photo_add_term_r (photo_term_t *term, cdr_grid_t *cdr)
+{
+  cdr_grid_t *child;
+
+  photo_add_term (term, cdr);
+
+  iter_childs (cdr, child) {
+    photo_add_term_r (term, child);
+  }
+}
+
+/** @brief photo_calc_term QQQQ */
+pois_grid_t **
+photo_calc_term (photo_term_t *term, cdr_grid_t *cdr, int i)
+{
+  /* Call to the generic Poisson/Helmholtz solver */
+  if (i == 1) {
+	return pois_gen_solve_a (cdr, pois_photo_2, photo_mappers, term->lambda);
+  } else {
+	return pois_gen_solve_a (cdr, pois_photo_1, photo_mappers, term->lambda);
+  }
+}
+
+/** @brief Calculates the photoionization and adds it to the derivatives
+ *  of the species densities. */
+void
+photo_calc (photo_term_t *terms, cdr_grid_t *cdr)
+{
+  pois_grid_t **pois_modes;
+  photo_term_t *term;
+  // photo_term_t *ptr_term;
+
+  photo_copy_source_r (cdr);
+  
+  if (cdr->ntheta != 1)
+    cdr_dft_charge_r (cdr, 1);
+  
+  // for (ptr_term = terms; ptr_term != NULL; ptr_term = ptr_term->next) {
+  //   printf("terms: A=%g lambda=%g ptr=%d\n", ptr_term->A, ptr_term->lambda, (int)ptr_term->next);
+  // }
+
+  int j = 0;
+  for (term = terms; term; term = term->next ) {
+    int i;
+    pois_modes = photo_calc_term (term, cdr, j);
+    j++;
+    if (cdr->ntheta != 1)
+      photo_dft_r (cdr, -1);
+
+    photo_add_term_r (term, cdr);
+    debug (3, "photo_calc (" photo_printf_str ", " grid_printf_str ")\n",
+        	 photo_printf_args(term), grid_printf_args(cdr));
+
+    /* Free the allocated memory. */
+    for (i = 0; i < max_ntheta; i++) {      
+      pois_free_r (pois_modes[i]);
+    }
+    free (pois_modes);
+  }
+}
+
+/** @brief Loads a photoionization file, consisting in a series of rows with
+   @a A and @a lambda. */
+void
+photo_load_file (char *fname)
+{
+  FILE *fp;
+  double A, lambda;
+  int c;
+  int i;
+
+  printf("\n");
+  printf ("Loading photoionization data from `%s'...\n", fname);
+  printf("\n");
+  fp = fopen (fname, "r");
+
+  if (NULL == fp) {
+    warning ("Unable to open photoionization file `%s'\n", fname);
+    exit (-1);
+    return;
+  }
+
+  
+  i=0;
+  do {
+    c = fscanf (fp, "%lf %lf", &A, &lambda);
+    if (c != 2) {
+      break;
+    }
+
+    printf ("Registring photoionization term A = %.5g, lambda = %.4g\n",
+	    A, lambda);
+
+    photo_register (A, lambda);
+    i++;
+  } while (TRUE);
+  printf("\n");
+
+  fclose (fp);
+}