/** @file reaction.c

 *  @brief Functions to handle the "reaction" part of

 *  convection-diffusion-reaction equation.

 */

#include <stdlib.h>

#include <stdio.h>

#include <math.h>

#include "cdr.h"

#include "grid.h"

#include "parameters.h"

#include "photo.h"

#include "proto.h"

#include "react_table.h"

#include "species.h"

#include "reaction.h"

static void fill_react_gaps ();

reaction_t *reactions_list = NULL;

extern double z_cutoff;

double species_written[20];

/** @brief Returns the position in the species-array of a given species.

 *

 * Returns -1 if not found. */

int

find_species_by_name(const char *spec_name)

{

  int res = -1;

  int cnt;

  for (cnt = 0; cnt < no_species; cnt++)

  {

    if (strcmp(spec_index[cnt]->name, spec_name) == 0) 

      res = cnt;

  }

  if (res == -1) {

      printf("Species-lookup-failure: %s!\n",spec_name);

      assert(1 == 0);

  }

  return res;

}

/** @brief Adds a given reaction to the reaction list */

void

react_add (reaction_t *react)

{

  react_table *rt;

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

  printf("Adding reaction with table: %s\n",react->tablefile);

  rt = (react_table *) xmalloc (sizeof(react_table));

  if (react->tablefile != NULL) {

    react_table_read(react->tablefile, rt);

    react->rt = rt;

  }

  react->next = reactions_list;

  reactions_list = react;

}

/** @brief Applies a reaction to the given @a grid */

void 

react_apply (reaction_t *react, cdr_grid_t *grid, int overwrite)

{

  int i, ir, iz, itheta;

  rz_array_t *grid_in[REACTION_MAX_IN], *grid_out[REACTION_MAX_IN + REACTION_MAX_OUT],

    *grid_eabs;

  double *in = NULL, *out = NULL, eabs;

  double test;

  double rate;

  int pos;

  double e1, e2, val1, val2, log_e, res, r_mod;

  int cnt, curr_species;

  debug (3, "react_apply (..., " grid_printf_str ")\n",

         grid_printf_args (grid));

  grid_eabs = grid->dens[no_species];

  for (i = 0; i < react->nin; i++) {

    grid_in[i] = grid->dens[react->input[i]];

    grid_out[i] = grid->d_dens[react->input[i]];

  }

  for (i = react->nin; i < react->nin + react->nout; i++) {

    grid_out[i] = grid->d_dens[react->output[i - react->nin]];

  }

#pragma omp parallel private(ir, iz, i, in, out)

  {

    /* malloc(0) is legal, but I do not want to play with fire. */

    if (react->nin > 0) {

      in = (double *) xmalloc (sizeof(double) * react->nin);

    }

    /* Do not know what use nout == 0 may have, (perhaps some debugging?)

       but we leave it here as theoretically possible. */

    if (react->nout > 0) {

      out = (double *) xmalloc (sizeof(double) * (react->nout + react->nin));

    }

    #pragma omp for

    iter_grid_theta(grid, itheta) { //ITER3

      iter_grid_z(grid, iz) { //ITER2

        double back_dens;

        if (sprite_module) {

          back_dens = spr_density_at (z_at (iz, grid->level));

        } else {

          back_dens = 1.0;

        }

        iter_grid_r(grid, ir) { //ITER1

          if ( z_at (iz, grid->level) < z_cutoff ) { //IF1

            eabs = fabs (RZT(grid_eabs, ir, iz, itheta));

            for (i = 0; i < react->nin; i++) {

              in[i] = fabs (RZT(grid_in[i], ir, iz, itheta));

	    }

            log_e = log10(eabs);

            /* If the supplied fieldstrength falls outside the boundaries of the table,

              return predetermined under-/overflow values */

            if (log_e < react->rt->e_min) {

              rate = react->rt->underflow;

            } else if (log_e > react->rt->e_min + react->rt->e_step * react->rt->steps) {

              rate = react->rt->overflow;

            } else {

              pos = floor((log_e - react->rt->e_min) / react->rt->e_step);

              val1 = react->rt->values[pos];

              val2 = react->rt->values[pos+1];

              e1 = pow(10, react->rt->e_min + react->rt->e_step * (double) pos);

              e2 = e1 * pow(10, react->rt->e_step);

              rate = val1 + ((val2 - val1) / (e2 - e1)) * (eabs - e1);

            }

            for (i = 0; i < react->nin; i++){ rate *=  MYMAX(0, in[i]); }

            for (i = 0; i < react->nin + react->nout; i++) { //FOR1

              if (i < react->nin) {

                curr_species = react->input[i]; r_mod = -rate;

              } else { curr_species = react->output[i - react->nin]; r_mod = rate; }

              if (spec_index[curr_species]->charge != 0.0) {

                 RZT(grid_out[i], ir, iz, itheta) += r_mod;

              }

            } //FOR1

          } //IF1

        } //ITER1

      } //ITER2

    } //ITER3

    free (in);

    free (out);

  }

}

/** @brief Recursive version of @a react_apply  */

void

react_apply_r (reaction_t *react, cdr_grid_t *grid, int overwrite)

{

  cdr_grid_t *child;

  react_apply (react, grid, overwrite);

  iter_childs (grid, child) {

    react_apply_r (react, child, overwrite);

  }

}

/** @brief Sets the d_dens field of a grid to zero */

void 

zero_fill (cdr_grid_t* grid)

{

  int ir, iz, itheta, i;

  iter_grid_3d (grid, ir, iz, itheta)

  for (i = 0; i < no_species; i++) {

    RZT(grid->d_dens[i], ir, iz, itheta) = 0.0;

  }

}

/** @brief Applies all reactions to the given @a grid and his descendants */

void 

react_apply_all (cdr_grid_t *grid)

{

  reaction_t *react;

  int overwrite;

  int last = -1;

  int cnt;

  zero_fill(grid);

  overwrite = TRUE;

  for (react = reactions_list; react; react = react->next) {

    if (react->is_photo) {

      photo_calc (photo_terms, grid);

    } else {

      react_apply_r (react, grid, overwrite);

    }

    overwrite = FALSE;

  }

}

/** @brief Fill in the gaps that we left in the definitions */

static void

fill_react_gaps ()

{

}

/** @brief Initializes the list of reactions. */

void 

react_init ()

{

  /* Note that the reactions are applied in inverse order than listed here. */

  /* Rest of the kinetic model: */

  kinetic_init ();

}

/** Below this electric field, we do not waste time calculating anything.

   Besides, this avoid NaNs for eabs == 0. */

#define EPS_EABS 1e-6