/*!\file Conductance.C
\brief  Conductance class definition.
*/
#include "Conductance.h"
#include <math.h>

//Range all;

/********************************************************************//**
 *  \author Slava Merkin (vgm at bu.edu)
 *  \since 08-23-2005 
 *
 *  \param[in] solver The pointer to the Solver class. This way we
 *  pass the necessary solver information, e.g. solver grid dimensions
 *  and coordinates
 *
 *  \param[in] params The pointer to the Params class. Mostly this is
 *  done to avoid passing a long list of conductance model parameters.
 *
 * The constructor gets dimensions of the Solver grid and redims the
 * necessary arrays accordingly. The conductance model parameters are
 * initialized from values contained in the Params class.
 ************************************************************************/
Conductance::Conductance(Solver *solver, const Params * const params)
  :euvModelType(params->euvModelType),   // INITIALIZE CONSTANTS
   alpha(params->alpha),               
   beta(params->beta),
   r(params->r),
   f107(params->f107),
   pedMin(params->pedmin),
   hallMin(params->hallmin),
   ped0(params->ped0),
   sigmaRatio(params->sigma_ratio),       
   constSigma(params->const_sigma)       // END INITIALIZE CONSTANTS
{
  nLat = solver->nLat;
  nLon = solver->nLon;

  pedersen.redim(nLat,nLon);
  hall.redim(nLat,nLon);
  avgEng.redim(nLat,nLon);
  numFlux.redim(nLat,nLon);
  engFlux.redim(nLat,nLon);
  deltaSigmaP.redim(nLat,nLon);
  deltaSigmaH.redim(nLat,nLon);
  euvSigmaP.redim(nLat,nLon);
  euvSigmaH.redim(nLat,nLon);
  sigmaP.redim(nLat,nLon);
  sigmaH.redim(nLat,nLon);

  jpara.redim(nLat,nLon);
  cs.redim(nLat,nLon);
  rho.redim(nLat,nLon);

  /** We initialize the necessary arrays to zero so that getTotal()
  * can be called without prior calling of precipitation(), and in
  * that case it's just supposed to give us the EUV background
  * contribution. So, we want #deltaSigmaP and
  * #deltaSigmaH to be initialized to 0.  Initialize
  * everything else just for consistency.
  */
  pedersen    = 0.;
  hall        = 0.;
  avgEng      = 0.;
  numFlux     = 0.;
  engFlux     = 0.;
  deltaSigmaP = 0.;
  deltaSigmaH = 0.;
  euvSigmaP   = 0.;
  euvSigmaH   = 0.;
  sigmaP      = 0.;
  sigmaH      = 0.;
  

  // Need to get coords from the solver grid
  x = solver->x;
  y = solver->y;

  cout << "Conductance constructor is here... " << endl;
}

/********************************************************************//**
 * \author Slava Merkin (vgm at bu.edu)
 * \since 08-23-2005 
 *
 * \param[in] tilt Magnetic dipole tilt
 *
 * \param[out] sigmaPFunction A grid function containing the Pedersen
 * conductance on the Solver grid
 *
 * \param[out] sigmaHFunction A grid function containing the Hall
 * conductance on the Solver grid
 *
 * This function is called from the main level program when it is
 * needed to calculate ionospheric condutances locally without
 * communication with an ITM model. The function precipitation is
 * assumed to have been called prior to calling getTotal(), or
 * otherwise the conductance will only contain the EUV contribution
 * (auroral contributions are initialized to zero in the contructor of
 * the class). The total conductances are calculated as follows. First
 * aurora() and euv() functions are called that fill in #euvSigmaP,
 * #euvSigmaH, #deltaSigmaP, and #deltaSigmaH arrays. Then,
 * \f{eqnarray*} \Sigma_P &= \sqrt{\Sigma_{P\_EUV}^2 + \delta
 * \Sigma_P^2}\\ \Sigma_H &= \sqrt{\Sigma_{H\_EUV}^2 + \delta
 * \Sigma_H^2} \f} where \f$\Sigma_{P\_EUV}\f$ and
 * \f$\Sigma_{H\_EUV}\f$ are EUV ionization contirbutions, while
 * \f$\delta \Sigma_P\f$ and \f$\delta \Sigma_H\f$ are contributions
 * due to auroral precipitation.
 ************************************************************************/
void Conductance::getTotal(realCompositeGridFunction &solverFunctions,
			   const real &tilt, const bool &HEMISPHERE)
{
  int SIGMAP, SIGMAH, AVG_ENG, NUM_FLUX;
  if (HEMISPHERE==NORTH) {
    SIGMAP   = SOLVER_SIGMAP_NORTH;
    SIGMAH   = SOLVER_SIGMAH_NORTH;
    AVG_ENG  = SOLVER_AVG_ENG_NORTH;
    NUM_FLUX = SOLVER_NUM_FLUX_NORTH;
  } 
  else {
    SIGMAP = SOLVER_SIGMAP_SOUTH;
    SIGMAH = SOLVER_SIGMAH_SOUTH;
    AVG_ENG     = SOLVER_AVG_ENG_SOUTH;
    NUM_FLUX    = SOLVER_NUM_FLUX_SOUTH;
  }

  // October 8, 2008: Set the local members avgEng and numFlux from
  // solverFunctions to get these values from the correct pole.
  Index I1(0,nLat), I2(0,nLon);
  avgEng = solverFunctions[0](I1,I2,all,AVG_ENG);
  numFlux = solverFunctions[0](I1,I2,all,NUM_FLUX);

  if (!constSigma) 
    {
      // note, precipitation parameters should be defined by now 
      // (call the precipitation() function prior to calling calc
      aurora();
      euv(tilt);
      
      sigmaP = sqrt( pow(euvSigmaP,2.) + pow(deltaSigmaP,2.) );
      /**\note The sign of the Hall conductance is positive, unlike
       * the conventional LFM implementation. The positive Hall
       * conductance is implied in the Solver code */
      sigmaH = sqrt( pow(euvSigmaH,2.) + pow(deltaSigmaH,2.) ); 
    }
  else
    {
      sigmaP = ped0;
      sigmaH = 0.;
    }

  // Put the arrays calculated into the grid functions for output
  solverFunctions[0](I1,I2,all,SIGMAP) = sigmaP;
  solverFunctions[0](I1,I2,all,SIGMAH) = sigmaH;
}

/********************************************************************//**
 * Apply a cap on the Pedersen & Hall conductances.  This is done to
 * prevent the code from blowing up spectacularly due to high
 * conductances on the night side of the low Latitude boundary
 ************************************************************************/
void Conductance::applyCap(realCompositeGridFunction &solverFunctions)
{
  if (!constSigma) {
    Index I1(0,nLat), I2(0,nLon);
    const int PED_N = SOLVER_SIGMAP_NORTH;
    const int PED_S = SOLVER_SIGMAP_SOUTH;
    const int HALL_N = SOLVER_SIGMAH_NORTH;
    const int HALL_S = SOLVER_SIGMAH_SOUTH;
    // Pedersen cap
    solverFunctions[0](I1,I2,all,PED_N) = max(pedMin, solverFunctions[0](I1,I2,all,PED_N));
    solverFunctions[0](I1,I2,all,PED_S) = max(pedMin, solverFunctions[0](I1,I2,all,PED_S));
    // Hall cap
    solverFunctions[0](I1,I2,all,HALL_N) = min(max(hallMin, solverFunctions[0](I1,I2,all,HALL_N)),
					       solverFunctions[0](I1,I2,all,PED_N)*sigmaRatio);
    solverFunctions[0](I1,I2,all,HALL_S) = min(max(hallMin, solverFunctions[0](I1,I2,all,HALL_S)),
					       solverFunctions[0](I1,I2,all,PED_S)*sigmaRatio);
  }
}

/********************************************************************//**
 * \author Slava Merkin (vgm at bu.edu)
 * \since 06-2007
 *
 * \param[in] jparaFunction The field-aligned current grid function on the Solver grid
 *
 * \param[in] rhoFunction The density grid function on the Solver grid
 *
 * \param[in] csFunction The sound speed grid function on the Solver grid
 *
 * \param[in] HEMISPHERE A logical variable defining the hemisphere
 *
 * \param[out] avgEngFunction The average precipitation energy grid
 * function on the Solver grid
 *
 * \param[out] numFluxFunction The number flux grid function on the Solver grid
 *
 * This function converts the incoming grid functions into arrays
 * (0:#nLat,0:#nLon), where #nLat and #nLon are initialized from
 * Solver::nLat and Solver::nLon in the Conductance constructor. Once
 * fedder95() is called and arrays #avgEng and #numFlux are filled
 * in, the outcoming grid functions are filled using data in these
 * arrays. Essentially, this function is the upper level wrapper for
 * fedder95().
 ************************************************************************/
void Conductance::precipitation(realCompositeGridFunction &solverFunctions, 
				const real &tilt, const bool &HEMISPHERE)
{
  int JPARA, DENSITY, SOUND_SPEED, AVG_ENG, NUM_FLUX, jparaSign;

  /** 
   * Note that by convention, the precipitation calculation in
   * fedder95() assumes upward current to be positive. In the MIX
   * code, unlike the conventional LFM code, the field-aligned current
   * is what it is, \f$ \vec j \cdot \hat b\f$, and therefore the
   * downward current is positive in the northern hemisphere and
   * negative in the southern hemisphere. Thus we multiply the current
   * passed to this function by -1 in the northern hemisphere before
   * feeding it to fedder95().
   */
  if (HEMISPHERE==NORTH) {
    JPARA       = SOLVER_FAC_NORTH;
    DENSITY     = SOLVER_DENSITY_NORTH;
    SOUND_SPEED = SOLVER_SOUND_SPEED_NORTH;
    AVG_ENG     = SOLVER_AVG_ENG_NORTH;
    NUM_FLUX    = SOLVER_NUM_FLUX_NORTH;
    jparaSign   = -1;
  } 
  else {
    JPARA       = SOLVER_FAC_SOUTH;
    DENSITY     = SOLVER_DENSITY_SOUTH;
    SOUND_SPEED = SOLVER_SOUND_SPEED_SOUTH;
    AVG_ENG     = SOLVER_AVG_ENG_SOUTH;
    NUM_FLUX    = SOLVER_NUM_FLUX_SOUTH;
    jparaSign   = 1;
  }

  // We make some copies here
  // Could probably do it somewhat more efficient by using links
  Index I1(0,nLat), I2(0,nLon);
  jpara = jparaSign*solverFunctions[0](I1,I2,all,JPARA);
  rho = solverFunctions[0](I1,I2,all,DENSITY);
  cs = solverFunctions[0](I1,I2,all,SOUND_SPEED);

  fedder95(tilt, HEMISPHERE);

  solverFunctions[0](I1,I2,all,AVG_ENG) = avgEng;
  solverFunctions[0](I1,I2,all,NUM_FLUX) = numFlux;
}


/********************************************************************//**
 * \author Slava Merkin (vgm at bu.edu)
 * \since 08-23-2005 
 *
 * \param[in] tilt Magnetic dipole tilt
 * 
 * \retval #euvSigmaP The function fills in values in this array which
 * has the class scope
 *
 * \retval #euvSigmaH The function fills in values in this array which
 * has the class scope
 *
 * This function calculates the EUV contribution to the ionospheric
 * conductances either using the "AMIE" algorithm implemented in the
 * conventional LFM (H. Kroehl, private communication, 1990) or the
 * Moen-Brekke formulae (Moen and Brekke, 1993). This depends on the
 * value of #euvModelType initialized from Params::euvModelType
 * in the constructor. The Moen-Brekke formulae is not fully
 * functional yet (08-2007) because it is only valid on the dayside
 * and I need to think how to deal with it on the night side (smooth
 * gradent or something similar). Also the new (correct?) zenith
 * angle definition is used (see #coszen) which is different from
 * the LFM implementation of the AMIE model and should be checked.
 *
 * The AMIE calculation proceeds as follows. \f$\Sigma_{P,H}\f$
 * equals one of 
 *
 * \f{eqnarray*} 
 * \Sigma_{P,H}^{(1)}, & \theta_{zenith} \leq 65^o 
 * \\ \Sigma_{P,H}^{(2)}, & 65 < \theta_{zenith} \leq 100^o 
 * \\ \Sigma_{P,H}^{(3)}, & \theta_{zenith} > 100^o \f} where
 *
 * \f{eqnarray*} 
 * \Sigma_P^{(1)} &= & 0.5 F_{10.7}^{2/3} \cos^{2/3}(\theta_{zenith})
 * \\ \Sigma_H^{(1)} &=& 1.8 F_{10.7}^{1/2}\cos(\theta_{zenith})
 * \\ 
 * \\ \Sigma_P^{(2)} &=& \Sigma_P^{(1)}(65^o)- 0.24 
 * \cdot 250^{-2/3}F_{10.7}^{2/3} (\theta_{zenith}-65^o)
 * \\ \Sigma_H^{(2)} &=& \Sigma_H^{(1)}(65^o) - 0.27 
 * \cdot 250^{-1/2}F_{10.7}^{1/2} (\theta_{zenith}-65^o)
 * \\ 
 * \\ \Sigma_P^{(3)} &=& \Sigma_P^{(2)}(100^o)-0.26 
 * \cdot 250^{-2/3} (\theta_{zenith}-100^o)\\ \Sigma_H^{(3)} &=& \Sigma_H^{(2)} \f}
 * and \f$\theta_{zenith} = \frac{\pi}{2} - \left( asin(\frac{x}{R_{ion}})
 * + \theta_{tilt} \right)\f$ and \f$\theta_{tilt}\f$ is the dipole
 * tilt.
 *
 * Unlike the calculation above (implemented in the LFM code), for
 * the Moen-Brekke formulae (and possibly default future use) the
 * zenith angle is defined as \f$
 * \cos(\theta_{zenith})=x\cos(\theta_{tilt})+\sqrt{1-x^2-y^2}\sin(\theta_{tilt})\f$,
 * where coordinates are, of course, in SM coordinate system
 *
 * \note We use thresholds on both conductances once they are
 * calculated: \f$\Sigma_P \ge \f$ #pedMin, \f$\Sigma_H \ge \f$
 * #hallMin, and \f$\Sigma_H \ge \Sigma_P\f$*#sigmaRatio
 ************************************************************************/
void Conductance::euv(const real tilt) 
{
  static real PI2       = PI/2.;
  static real ang65     = PI/180.*65.;
  static real ang100    = PI*5./9.;
  static real pref      = 2.*pow(250.,-0.666666);
  static real href      = 1./(1.8*sqrt(250.));
  static real rad2deg   = 180./PI;
  static real shall     = 1.8*sqrt(f107);
  static real speder    = 0.5*pow((real)f107,0.666666);
  static real pedslope  = 0.24*pref*speder*rad2deg;
  static real pedslope2 = 0.13*pref*speder*rad2deg;
  static real hallslope = 0.27*href*shall*rad2deg;
  static real sigmap65  = speder*pow((real)cos(ang65),0.666666);
  static real sigmah65  = shall*cos(ang65);
  static real sigmap100 = sigmap65 - (ang100-ang65)*pedslope;

  realArray zenith = PI2 - ( asin(x) + tilt );
  // An alternative (correct) definition of the zenith angle (for Moen-Brekke)
  realArray coszen = x*cos(tilt)+sqrt(1.-x*x-y*y)*sin(tilt) ; //as it should be
  switch ( euvModelType ) {
  case AMIE:
    where (zenith <= ang65) {
      euvSigmaP = speder*pow(cos(zenith),0.666666);
      euvSigmaH = shall*cos(zenith);
    }
    elsewhere(zenith <= ang100) {
      euvSigmaP = sigmap65 - pedslope*(zenith - ang65);
      euvSigmaH = sigmah65 - hallslope*(zenith - ang65);
    }
    elsewhere(zenith > ang100) {
      euvSigmaP = sigmap100 - pedslope2*(zenith-ang100);
      euvSigmaH = sigmah65 - hallslope*(zenith-ang65);
    }
    break;
    
  case MOEN_BREKKE: 
    // This works only for the dayside (zenith <= pi/2)
    // Set it to pedMin (hallMin) on the nightside (may rethink it later)
    
    where (coszen >=0.) {
      //	euvSigmaP = pow((real)f107,0.49)*( 0.34*cos(zenith)+0.93*sqrt(cos(zenith)) );
      //	euvSigmaH = pow((real)f107,0.53)*( 0.81*cos(zenith)+0.54*sqrt(cos(zenith)) );
      euvSigmaP = pow((real)f107,0.49)*( 0.34*coszen+0.93*sqrt(coszen) );
      euvSigmaH = pow((real)f107,0.53)*( 0.81*coszen+0.54*sqrt(coszen) );
    }
    elsewhere() {
      euvSigmaP = pedMin;
      euvSigmaH = hallMin;
    }
    break;
    
  default:
    cout << "The EUV model type entered is not supported."<<endl;
    exit(1);
  }

   // MJW - Added a step to make conductance is never less than the floor

  // (VGM-070711 - moved it after the switch statement; otherwise this
  // was never executed). Also, using more compact form.
//    where(euvSigmaP < pedMin) {
//      euvSigmaP = pedMin;
//    }
//    where(euvSigmaH < hallMin) {
//      euvSigmaH = hallMin;
//    }
  euvSigmaP = max(euvSigmaP,pedMin);
  euvSigmaH = max(euvSigmaH,hallMin);


}


/********************************************************************//**
 * \author Slava Merkin (vgm at bu.edu)
 * \since 08-23-2005 
 *
 * \retval #engFlux The energy flux in ergs/cm^2 is calculated by this function
 * but it is never used outside so far.
 *
 * \retval #deltaSigmaP array is filled in
 *
 * \retval #deltaSigmaH array is filled in
 *
 * Based on the average energy (#avgEng) and number flux (#numFlux)
 * calculated by #fedder95(), which must be called in advance to fill
 * in these arrays, this function calculates the auroral contribution
 * to the conductances according to Robinson et al. [1987, eqs. 3 and
 * 4]:
 *
 * \f{eqnarray*} \delta \Sigma_P &=&
 * \frac{40.\bar E}{16.+\bar E^2}\Phi_E^{1/2}
 * \\ \delta \Sigma_H &=&0.45 \delta \Sigma_P \bar E^{0.85} \f} 
 *
 * where \f$\bar E\f$ is the average energy in keV and \f$ \Phi_E \f$
 * is the energy flux in ergs/cm^2/s.
 *
 * \note In the conventional LFM code a slight modification to
 * Robinson is used \f$ \delta \Sigma_H = 0.45 \delta \Sigma_P
 * \frac{\bar E^{0.85}}{1.+(0.05\bar E)^2} \f$.
 ************************************************************************/
void Conductance::aurora() {
  engFlux = kev2erg*avgEng*numFlux;  //Energy flux in ergs/cm^2/s
  realArray factor1(nLat,nLon);

 /* --------- Lower influence of flux near lower boundary. Contact S. Slinker ---------- */
  realArray R = sqrt( x*x + y*y );
  real Redge = max(R);
  real Rmin = 0.85*Redge;
  real Rmin2= 0.45*Redge;
  real Rmax = 1.0*Redge;
  real rfac = 0.1;

  where     ( R <= Rmin ) { 
    factor1 = 1.; 
  }
  elsewhere ( R <= Rmax ) { 
    factor1 = 1. + (rfac - 1.)*(R - Rmin)/(Rmax - Rmin);
  }

  deltaSigmaP = 40.*factor1*avgEng*sqrt(engFlux)/(16.+avgEng*avgEng);
  deltaSigmaH = 0.45*deltaSigmaP*pow(avgEng,0.85)/(1.+0.0025*avgEng*avgEng);
}


/********************************************************************//**
 * \author Slava Merkin (vgm at bu.edu)
 * \since 08-23-2005 
 *
 *\retval #avgEng array is filled in
 *
 *\retval #numFlux array is filled in
 *
 * Based on the sound speed (#cs), density (#rho) and parallel current
 * (#jpara) received from the MHD code, calculates the average energy
 * (#avgEng) and number flux (#numFlux) of precipitating electrons
 * according to [Fedder el al., 1995].
 ************************************************************************/
void Conductance::fedder95(const real &tilt, const bool &HEMISPHERE) {
  realArray E0(nLat,nLon);             // initial guess for energy
  realArray phi0(nLat,nLon);           // initial guess for flux
  realArray deltaE(nLat,nLon);         // potential drop along FAC
  realArray aRes(nLat,nLon);           // anomalous resistivity
  realArray factor1(nLat,nLon), factor2(nLat,nLon);
  
 /* --------- Lower influence of flux near lower boundary. Contact S. Slinker ---------- */
  realArray R = sqrt( x*x + y*y );
  real Redge = max(R);
  real Rmin = 0.85*Redge;
  real Rmin2= 0.45*Redge;
  real Rmax = 1.0*Redge;
  real rfac = 0.1;

  where     ( R <= Rmin ) { factor1 = 1.; }
  elsewhere ( R <= Rmax ) { factor1 = 1. + (rfac - 1.)*(R - Rmin)/(Rmax - Rmin); }

  // The y-axis points from dawn to dusk in the northern hemisphere,
  // but from dusk to dawn in the southern one. So the sign of "y"
  // there should be switched to "+" for the southern hemisphere
  // within the "elsewhere" below.
  real signOfY = 1.0;
  if (HEMISPHERE == NORTH){
    signOfY = -1.0;
  }

  where     ( R <= Rmin2 ) { factor2 = 1.; }
  // 11-05-05 **VGM** Commented out the mask below in the argument of
  // elsewhere due to a bug(?). As long as Rmax = Redge, it's
  // fine. With the mask, at least for a double resolution situation,
  // factor2 is set to 1.e+07 somewhere at the boundary for an uknown
  // reason
  elsewhere ( /* R <= Rmax */  ) { factor2 = 1. + 0.33*signOfY*y/R*(R - Rmin2)/(Rmax - Rmin2); }

  /* This version of the Ramp function overrides the implentation used
 *   in the OMP version of the LFM with a new version.  This version
 *   is desgined to have the function which has the loss cone full
 *   inside the polar cap and nearly empty below the nominal
 *   equatorward edge of the auoral oval.  See [P Janhunen, 1996] for
 *   a good justification of this factor
 */
  factor1 = 1.0;
  factor2 = 1.0;
  real rDeg2Rad = PI/180.0;
  real rPolarBound = 20.0*rDeg2Rad;
  real rEquatBound = 30.0*rDeg2Rad;
  real rLowLimit = 0.02;
  realArray RampFactor(nLat,nLon);
  where (R < rPolarBound) {RampFactor = 1.0;}
  where (R > rPolarBound) {RampFactor = 1.0+(rLowLimit - 1.0)*
                         (R - rPolarBound)/(rEquatBound-rPolarBound);}
  where (R > rEquatBound) {RampFactor = rLowLimit;}

  /* ---------------------- End corrections near lower boundary --------------------- */

  /**
   * 1. Calculate \f$ E_0 \f$ according to [Fedder et al., 1995, Eq. 3]:
   * \f[
   * E_0 = \alpha \cdot C_{sound}^2
   * \f]
   * which becomes \f$ \alpha \cdot M \cdot erg2kev \cdot C_{sound}^2 \f$ in 
   * order to have \f$ E_0 \f$ have units of keV. Here, \f$ M = m_p \cdot heFrac \f$ 
   * includes the helium correction (in g) and \f$ C_{sound} \f$ is in cm/s.
   * Evaluated this is 
   * \f[
   * E_0 = 12.11 \cdot 10^{-16} \cdot \alpha \cdot C_{sound}^2 \mbox{ [keV]}
   * \f]
   */

  /*
    Take into account the factor at the lower boundary
     Eo = alpha*mp*hefrac*erg2kev*cs**2*factor2
  */
  E0 = alpha*mp*heFrac*erg2kev*cs*cs*factor2;
  E0 = E0*RampFactor;


  /**
   * 2. Calculate \f$ \Phi_0 \f$ according to [Fedder et al., 1995, Eq. 4]:
   * \f[
   * \Phi_0 = \beta \cdot \rho E^{1/2}_0
   * \f]
   * which becomes \f$ \beta \cdot M^{-1} \cdot \rho \cdot kev2erg^{1/2} 
   * \cdot M^{-1/2} \cdot E_0^{1/2} \f$
   * in order to have \f$ \Phi_0 \f$ in units of \f$ [cm^{-2} s^{-1} ] \f$.
   * Evalulated this is 
   * \f[
   * \Phi_0 = 1.48 \cdot 10^{31} \cdot \beta \cdot \rho \cdot E_0^{1/2}\: [cm^{-2} s^{-1}]
   * \f]
   */

  phi0 = sqrt(kev2erg)/pow(heFrac*mp,1.5)*beta*rho*sqrt(E0);
  phi0 = phi0*RampFactor;

  /**
   * 3. According to [Chiu et al., 1981, Eq. 2] the parallel
   * electric potential energy is given by
   * \f[
   * E_\parallel = e \Delta \phi \sim \frac{e j_\parallel M \bar v}{N e^2}
   * \f]
   * Considering, that \f$ \rho = MN \f$ and \f$ E_0 \sim M \bar v^2 \f$ we obtain
   * \f[
   * E_\parallel \sim \frac{j_\parallel M^{3/2} E_0^{1/2}}{\rho e}
   * \f]
   * i.e. Eq. 5 from [Fedder et al., 1995]:
   * \f[
   * E_\parallel = R \frac{j_\parallel E_0^{1/2}}{\rho} 
   * \mbox{, where } R = \frac{cM^{3/2}}{e} 
   * \f]
   * and \e c is a dimensionless constant.
   *
   *  The coefficient R is taken bo be r*aRes, where #r is read from
   *  the parameter file and aRes is taken to be 5 times larger for
   *  current out of the ionosphere than for current into the
   *  ionosphere.
   */

  /*  
      resistence out of the ionosphere is 2*rout resistence into the
      ionosphere is 2*rin outward current is positive
  */
  real Rout = 6.;
  real Rin  = 1.2;
  where     ( jpara >=0. ) { aRes = 2.*Rout*factor1; }
  elsewhere( jpara < 0.)        { aRes = 2.*Rin*factor1; }

  /**
   * The coefficient in Eq. 5 from [Fedder et al., 1995] becomes \f$ R
   * \cdot M^{3/2} \cdot e^{-1} \cdot 10^{-4} \cdot ergs2kev^{1/2} \f$
   * in order for R to be dimensionless.

   * \note The factor of 1e-4 comes from the m to cm conversion
   * required to move the current density from \f$ A/m^2 \f$ to \f$
   * A/cm^2 \f$.
   */
  /* Impose floor on density to limit characteristic energy.  See Wiltberger et
   * al. 2009 for details.
   */
  real rhoFactor = 3.3e-24*0.5;
  euv(tilt);
  where ( rho < rhoFactor*euvSigmaP ) { rho = rhoFactor*euvSigmaP; }
  deltaE = (pow(heFrac*mp,1.5)/eCharge)*1.e-4*sqrt(erg2kev)*r*aRes*jpara*sqrt(E0)/rho;

  /*
	  Remove potential drops for cells with  nj > nj - nkill (Fortran indexing)
	  deltaE(J,all) = 0.; Note, the range J depends on the array base
  */

  /*
    07-06-06-VGM Removed the nkill code below because I've made this code LFM grid independent
    int nkill = 4;
    int last_j =  deltaE.getBound(0);
    Range J(last_j-nkill-1,last_j);
    deltaE(J,all) = 0.;           
   */
  
  /*      limit the max potential energy drop to 20 [keV]  */
  where ( deltaE >= 20. ) { deltaE = 20.; }

  /**
   * 4. Use [Fedder et al., 1995, Eq. 8] to calculate the
   * total energy of precipitating electrons:
   * \f[
   * E = E_0 + E_\parallel
   * \f]
   */
  avgEng = E0 + deltaE;

  /*       do not allow the avgEng to be less then 1.e-8 */
  where ( avgEng <= 1.e-8 ) { avgEng = 1.e-8; }

  /**
   * Finally apply [Fedder et al., 1995, Eq. 6, 7] to adjust the number flux
   * \f{eqnarray*}
   * \Phi &=& \Phi_0 \left( 8-7 \exp{\frac{-E_\parallel}{7E_0}} \right) 
   * \mbox{ for } E_\parallel > 0 \\
   * \Phi &=& \Phi_0 \exp{\frac{E_\parallel}{E_0}} \mbox{ for } E_\parallel < 0
   * \f}
   *
   * The resulting energy units are keV and number flux units are  \f$ cm^{-2} s^{-1}\f$. 
   */
  where     ( deltaE > 0 ) { numFlux = phi0*(8.-7.*exp(-deltaE/7./E0)); }
  elsewhere ( )            { numFlux = phi0*exp(deltaE/E0); }
}