/** @file cstream.c
 *  @brief General initialization/termination functions.
 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#define ALLOC_PARAMS

#include "cstream.h"
#include "parameters.h"
#include "proto.h"
#include "species.h"

double E_x, E_y, E_z;
/**< For time-dependent electric fields E_x, E_y, E_z there is a difference between
   E0 (the highest field) and E (the field at a given time). */

extern double pois_inhom_fixed_q_t;
/**< The same is true for time-dependent fixed charges (useful to model
   the charging of a cloud in modeling sprites. */

double *dr, *dz, dtheta;
/**< dr[i] and dz[i] represent the grid size at level @a i.  Note that
   @a i can go from -extra_pois_levels...max_levels. */
double *dr_start, *dz_start;

double *w2k, *wk;
/**< w2k[k] = 2 / dtheta^2 (1 - cos(k dtheta)). This numbers appear due to the
   finite-differences derivatives of a Fourier transform.  If we would make
   a continuous fourier transform, they will approach k^2.
*/


char *invok_name = "cstream";
/**< Name of this program. */

const double twopi = 6.283185307179586477L;
/**< \f$ 2 \times \pi \f$. */

const double invfourpi = 0.079577471545947667884441882L;
/**< \f$ 1 / (4 \pi )\f$ */

const double invpi32 = 0.179587122125166561689081L;
/**< \f$ 1 / (\pi ^ {3/2} )\f$ */

/** @brief Initializes the parameters.
 *
 *  The global parameters will get their default values.
 *  These values correspond with the values in input/defaults.cfg
 */
void
init_parameters(void)
{
    debug (1, "entry init_parameters\n");
    /* An identification name for this run */
    prog_id = "example";

    /* Output directory.  */
    output_dir = "output/";

    /* Kinetics input file */
    kin_input = "input/kinetic_example.cfg";

    /* If restart is TRUE, the simulation will continue with data from a previous run */
    restart = 0;

    /* If restart is TRUE, the name of the file with data from previous run, otherwise empty */
    load_file = "";

    /* Time interval for output to be written to disk */
    output_dt = 0.100;

    /* Output of the Poisson grids, including the potential? */
    pois_output = 0;

    /* Margin outside the grids in the output of the cdr equation */
    cdr_output_margin =  0;

    /* Margin outside the grids in the output of the poisson equation */
    pois_output_margin = 1;

    /* If the time steps are smaller than this number, the program issues a warning */
    warn_min_timestep = 1e-06;

    /* Maximum disk space, in Mb, to use */
    max_disk_space_mb = 1048576;

    /* Number of R and Z gridpoints at level 0 */
    gridpoints_r = 600;
    gridpoints_z = 600;

    /* Number of azimuthal gridcells and modes */
    max_ntheta = 1;

    /* Initial and end time */
    start_t = 0.0;
    end_t = 0.12;

    /* Attempted timestep.  The actual timestep may be larger */
    attempt_dt = 50.0;

    /* Extra levels for the Poisson solver */
    extra_pois_levels = 2;

    /* Maximum level of refinement. Use a big number here */
    max_levels = 64;

    /* Error threshold that leads to refinement in the Poisson code. */
    pois_max_error = 0.001;

    /* Maximum level of refinement in the Poisson equation. */
    pois_max_level = 3;

    /* Extra levels for the photo-ionization solver */
    extra_photo_levels = -1;

    /* Maximum level of refinement in the photo-ionization solver. */
    photo_max_level = 4;

    /* Error threshold that leads to refinement in the photo-ionization code. */
    photo_max_error = 0.01;

    /* Photo-ionization boundary condition at r = L_r, z = 0, z = L_z.  1 for Hom. Neumann, -1 for Hom. Dirichlet */
    photo_bnd_right  = -1;
    photo_bnd_bottom = -1;
    photo_bnd_top    = -1;

    /* Extra levels for the photo-ionization solver */
    extra_photo_levels_2 = -1;

    /* Maximum level of refinement in the photo-ionization solver. */
    photo_max_level_2 = 4;

    /* Error threshold that leads to refinement in the photo-ionization code. */
    photo_max_error_2 = 0.01;

    /* Photo-ionization boundary condition at r = L_r, z = 0, z = L_z.  1 for Hom. Neumann, -1 for Hom. Dirichlet */
    photo_bnd_right_2  = -1;
    photo_bnd_bottom_2 = -1;
    photo_bnd_top_2    = -1;

    /* Particles boundary condition at z = 0, z = L_z, r = L_r.  1 for Hom. Neumann, -1 for Hom. Dirichlet */
    cdr_bnd_bottom = 1;
    cdr_bnd_top    = 1;
    cdr_bnd_right  = 1;

    /* Potential boundary condition at r = L_r, z = 0,  z = L_z.  1 for Hom. Neumann, -1 for Hom. Dirichlet */
    pois_bnd_right  = -1;
    pois_bnd_bottom = -1;
    pois_bnd_top    = -1;

    /* Maximum advection and diffusion Courant number */
    nu_a = 0.2;
    nu_d = 0.2;

    /* Maximum ratio of dt/relaxation time */
    nu_rt = 0.2;

    /* Maximum ratio of change of the densities (set to a very large number to ignore) */
    nu_f = 1e+20;

    /* Refinement threshold for the electric field */
    ref_threshold_eabs = 0.2;

    /* Maximum refinement level reached through ref_threshold_eabs */
    ref_level_eabs = 4;

    /* Refinement threshold for the curvature of the charge, densities */
    ref_threshold_charge = 0.004;
    ref_threshold_dens = 0.004;

    /* Refinement threshold for the densities in the leading edge */
    ref_threshold_edge = 10000.0;

    /* r-length and  z-length  of the minimal refinement area in the cdr equation */
    cdr_brick_dr = 8;
    cdr_brick_dz = 8;

    /* Maximum level of refinement in the Fluid equation. */
    cdr_max_level = 3;

    /* Interpolation method for the grid interior, and grid boundaries (0=zero_masses, 1=quadratic_masses [default], 2=wackers_masses, 3=quadlog */
    cdr_interp_in = 1;
    cdr_interp_bnd = 1;

    /* Length in r and z of the complete domain */
    L_r = 13044.0;
    L_z = 13044.0;

    /* Isotropic difussion coefficient */
    diffusion_coeff = 0.1;

    /* Whether the code includes photoionization or not */
    has_photoionization = 1;

    /* The name of a file from which we can read the photoionization parameters */
    photoionization_file = "input/air760torr.photo";

    /* Rate of dissociative attachment */
    attachment_rate = 0.0;

    /* E0 in the exp(-E0/E) factor in the attachment expression. */
    attachment_E0 = 0.0;

    /* x-, y- and z-component of the external electric field */
    E0_x = 0.0;
    E0_y = 0.0;
    E0_z = -0.06;

    /* Rise time of the electric field (0 for instantaneous rise) */
    rise_time = 0.0;

    /* Time to switch off the electric field (0.0 means never) */
    off_time = 0.0;

    /* x-, y- and z-width of the initial seed */
    seed_sigma_x = 0.0;
    seed_sigma_y = 0.0;
    seed_sigma_z = 0.0;

    /* Number of electrons in the initial seed */
    seed_N = 0.0;

    /* Initial at z=0 densities of electrons and ions */
    background_ionization = 0.0;

    /* Length of exponential increase of the pre-ionization (for atmospherical models) */
    background_increase_length = 0.0;

    /* Use the point-plane geometry? */
    pois_inhom = 1;

    /* Number of mirror charges to use */
    pois_inhom_reflections = 4;

    /* Length and radius of the needle */
    needle_length = 2500.0;
    needle_radius = 400.0;

    /* If nonzero, the charge is fixed, not floating (simulation of charged clouds close to the earth surface) */
    pois_inhom_fixed_q = 0.0;

    /* Constant ionization rate */
    constant_source = 0.0;

    /* Initial perturbation to the axisymmetric configuration */
    perturb_epsilon = 0.0;

    /* Perturb only modes up to perturb_max_k (large number to perturb all) */
    perturb_max_k = 1024;

    /* 1 if the sprite module is activated, 0 otherwise */
    sprite_module = 0;

    /* Lenght of exponential decay of the density w/r to altitude */
    dens_decay_len = 0.0;

    /* Density at z = 0 */
    sprite_dens_0 = 0.0;

    /* Quenching density */
    sprite_dens_q = 0.0;

    /* Sign of the sprite head that we are following (the other will not be reliable */
    sprite_sign = -1;
    debug (1, "exit init_parameters\n");
}

  int n, i;
/** @brief Initializes the grid sizes.  The parameters dr_root and dz_root
 *  specify the grid size at level 0 (dz_root = dz[0], dr_root[0] = dr[0])
 */
static void
init_gridsizes_a (void)
{
  int n, i;
  double root_dr, root_dz;

  root_dr = L_r / gridpoints_r;
  root_dz = L_z / gridpoints_z;

  n = max_levels + extra_pois_levels + 1;

  dr_start = (double *) xmalloc(sizeof(double) * n);
  dz_start = (double *) xmalloc(sizeof(double) * n);

  dr = dr_start + extra_pois_levels;
  dz = dz_start + extra_pois_levels;

  dr_start[0] = root_dr * (1 << extra_pois_levels);
  dz_start[0] = root_dz * (1 << extra_pois_levels);

  for (i = 1; i < n; i++) {
    dr_start[i] = dr_start[i - 1] / 2.0;
    dz_start[i] = dz_start[i - 1] / 2.0;

    debug (3, "dr[%d] = %e\n", i - extra_pois_levels,
	   dr[i - extra_pois_levels]);
    debug (3, "dz[%d] = %e\n", i - extra_pois_levels,
	   dz[i - extra_pois_levels]);
  }

  dtheta = twopi / max_ntheta;
}

/** @brief Initializes the vector of wk's for the Helmholtz equation. */
static void
init_wk_a (void)
{
  int k;
  double twobydtheta2;

  debug (2, "init_w2k_a()\n");

  w2k = (double*) xmalloc (sizeof(double) * max_ntheta);
  wk  = (double*) xmalloc (sizeof(double) * max_ntheta);

  twobydtheta2 = 2. / (dtheta * dtheta);

  for (k = 0; k < max_ntheta / 2 + (max_ntheta % 2); k++) {
    double w2k_ = twobydtheta2 * (1 - cos (k * dtheta));
    double re_wk, im_wk;

    re_wk = (cos (k * dtheta) - 1) / dtheta;
    im_wk = sin (k * dtheta) / dtheta;

    wk[k] = re_wk;
    w2k[k]= w2k_;

    if (k != 0) {
      wk[max_ntheta - k] = im_wk;
      w2k[max_ntheta - k] = w2k_;
    }
  }

  if ((max_ntheta % 2) == 0) {
    w2k[max_ntheta / 2] = twobydtheta2 * (1 - cos (k * dtheta));
    wk[max_ntheta / 2] = (cos (k * dtheta) - 1) / dtheta;
  }

  assert (wk[0] == 0.);
}

/** @brief Frees the space allocated for dr and dz. */
static void
free_gridsizes (void)
{
  free (dr_start);
  free (dz_start);
}

/** @brief Frees the allocated space for \f$ \lambda \f$. */
static void
free_wk (void)
{
  free (wk);
  free (w2k);
}

/** @brief  Here all the initialization calls. */
void
cstream_init (void)
{
  /* For the perturbations, it is better to start with different seeds at
     each run. */
  srand (time (0));

  init_gridsizes_a ();
  init_wk_a ();
  react_init ();
  cdr_init ();
  pois_init ();
  photo_init ();
  if (has_photoionization) {
    printf("Photoionization\n");
    photo_load_file (photoionization_file);
  } else {
    printf("No photoionization\n");
  }

  if (sprite_module) spr_init ();
}

/** @brief Here all the cleaning and finishing calls. */
void
cstream_end (void)
{
  photo_unregister_all ();
  free_gridsizes ();
  free_wk ();
  cdr_end ();
}


/** @brief When we use a time-varying external field, we update the field
 *  components at each timesteps. */
void
cstream_set_field_at_time (double t)
{
  static int rise_reached = FALSE, off_reached = FALSE;
  double factor;

  if (rise_reached && off_reached)
    return;

  if (t >= rise_time) {
    rise_reached = TRUE;
  }

  if (t >= off_time || off_time == 0.0) {
    off_reached = TRUE;
  }

  /* We make sure that we never use a field _larger_ than E0 and that
     if rise_time is 0 we do not calculate t / 0.0.  */
  factor = t >= rise_time? 1.0: (t / rise_time);

  factor = (off_time > 0.0 && t >= off_time)? 0.0: factor;

  E_x = factor * E0_x;
  E_y = factor * E0_y;
  E_z = factor * E0_z;

  if (pois_inhom_fixed_q != 0.0) {
    pois_inhom_fixed_q_t = factor * pois_inhom_fixed_q;
  }

  return;
}

/** @brief Functions for a constant external electric field given in (x, y, z) components. */
double
e0_r (double r, double z, double theta)
{
  return E_x * cos (theta) + E_y * sin (theta);
}

/** @brief Functions for a constant external electric field given in (x, y, z) components. */
double
e0_z (double r, double z, double theta)
{
  return E_z;
}

/** @brief Functions for a constant external electric field given in (x, y, z) components. */
double
e0_theta (double r, double z, double theta)
{
  return -E_x * cos (theta) + E_y * sin (theta);
}

/** Initializes the component of the external electric field in \f$r\f$-direction. */
decl_field_comp(r) = e0_r;
/** Initializes the component of the external electric field in \f$z\f$-direction. */
decl_field_comp(z) = e0_z;
/** Initializes the component of the external electric field in \f$\theta\f$-direction. */
decl_field_comp(theta) = e0_theta;