sprites.c

Go to the documentation of this file.

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>

#include "cdr.h"
#include "parameters.h"
#include "photo.h"
#include "proto.h"
#include "rz_array.h"
#include "species.h"

static int var_max_axis (cdr_grid_t *grid, rz_array_t *var, int sign);

double spr_nu_a, spr_nu_d, spr_nu_rt;
photo_term_t *spr_photo_terms;
extern photo_term_t *photo_terms;

double
spr_density_at (double altitude)
{
  return exp (-altitude / dens_decay_len);
}


void
spr_init ()
{
  /* We make the Courant numbers depend on the altitude: this is because
     when we apply the Courant-Levy criterium we assume v = -E, which is not
     true if the mobility now is != 1.0 */
  spr_nu_a = nu_a;
  spr_nu_d = nu_d;
  spr_nu_rt = nu_rt;

  photo_copy_list (photo_terms, &spr_photo_terms);
}

void
spr_hook (cdr_grid_t *grid)
{
  double altitude;
  
  altitude = spr_head_altitude (grid, sprite_sign);

  spr_update (altitude);
}

void
spr_update (double altitude)
{
  double back_dens = spr_density_at (altitude);
  photo_term_t *p0, *p;

  nu_a = spr_nu_a;
  nu_d = spr_nu_d * back_dens;
  nu_rt = spr_nu_rt;

  for (p0 = spr_photo_terms, p = photo_terms; p0; p0 = p0->next, p = p->next) {
    p->lambda = p0->lambda * back_dens * back_dens;
    p->A = p0->A * (sprite_dens_0 + sprite_dens_q) 
      / (back_dens + sprite_dens_q);
  }
  
}

double
spr_head_altitude (cdr_grid_t *grid, int sign)
{
  int ih;
  ih = var_max_axis (grid, grid->charge, sign);
  if (0 == ih) ih = var_max_axis (grid, grid->dens[electrons], 1);
  return z_at (ih, grid->level);
}

static int
var_max_axis (cdr_grid_t *grid, rz_array_t *var, int sign)
{
  int iz, iz_max = 0;
  double var_max = 0.0;

  iter_grid_z (grid, iz) {
    if (sign * RZT (var, 0, iz, 0) > var_max) {
      iz_max = iz;
      var_max = RZT (var, 0, iz, 0);
    }
  }
  
  return iz_max;
}