#include "udf.h"
/* const real R_Gas_C = 1.9872; Universal gas constant in cal/mole/K */ const real P_Atm = 101325.; /* unit atm in Pa */

DEFINE_NOX_RATE(user_nox, c, t, Pollut, Pollut_Par, NOx) {
  switch (Pollut_Par->pollut_io_pdf) {
  case IN_PDF:

POLLUT_FRATE(Pollut) = 0.0;
POLLUT_RRATE(Pollut) = 0.0;

switch (Pollut_Par->pollut_io_pdf) {

case IN_PDF:

 /* Included source terms other than those from char */

if (POLLUT_EQN(Pollut_Par) == EQ_NO) {

/* Thermal NOx */

if (NOx->thermal_nox && NOx->thermal_udf_replace) {

real kf1, kr1, kf2, kr2, kf3;

real o_eq=0.0, oh_eq=0.0;

real s1, s2, s3, rf, rr;

Rate_Const K_Z_F[3] = {{1.80e8, 0.0, 38370.0},

{1.80e4, 1.0, 4680.0},

{7.10e7, 0.0, 450.0}};

Rate_Const K_Z_R[3] = {{3.80e7, 0.0, 425.0},

{3.81e3, 1.0, 20820.0},

{1.70e8, 0.0, 24560.0}};

Rate_Const K_O[2] = {{36.64, 0.5, 27123.0},

{3.97e5, -0.5, 31090.0}};

Rate_Const K_OH = {2.129e2, -0.57, 4595.0};

kf1 = ARRH(Pollut, K_Z_F[0]);

kr1 = ARRH(Pollut, K_Z_R[0]);

kf2 = ARRH(Pollut, K_Z_F[1]);

kr2 = ARRH(Pollut, K_Z_R[1]);

kf3 = ARRH(Pollut, K_Z_F[2]);

/* determine O concentration */

if (NOx->o_mode == PARTIAL_EQUILIBRIUM) {

/* use partial equilibrium O model (Warnatz data) */

o_eq = ARRH(Pollut, K_O[0]);

}

else if (NOx->o_mode == FULL_EQUILIBRIUM) {

/* use full equilibrium O model (Westenburg) */

o_eq = ARRH(Pollut, K_O[1]);

}

else if (NOx->o_mode == LOCAL_MASS_FRACTION) {

/* use local [O] mass fraction */

o_eq = MOLECON(Pollut, O);

}

if ((NOx->o_mode == PARTIAL_EQUILIBRIUM || NOx->o_mode == FULL_EQUILIBRIUM))

o_eq *= sqrt(MOLECON(Pollut, O2));

/* determine OH concentration */

/* if local OH concentration is zero then calculate it from [O] and [H2O] */

if (NOx->oh_mode == PARTIAL_EQUILIBRIUM) {

/* (Westbrook&Baulch -- Andre Peters) */

oh_eq = ARRH(Pollut, K_OH);

oh_eq *= sqrt(o_eq*MOLECON(Pollut, H2O));

}

else if (NOx->oh_mode == LOCAL_MASS_FRACTION) {

/* use local OH mass fraction */

oh_eq = MOLECON(Pollut, OH);

}

else if (NOx->oh_mode == 0) {

/* do not include OH in NO calculation */

oh_eq = 0.0;

}

rf = rr = 0.0;

s2 = s3 = 0.0;

s1 = kf2*MOLECON(Pollut, O2);

if (s1 < 1e-6) {

rf = 2.*o_eq*kf1*MOLECON(Pollut, N2);

rr = 0.0;

POLLUT_FRATE(Pollut) = rf;
POLLUT_RRATE(Pollut) = rr;


}

else {

s2 = 1.+ kf3*oh_eq/s1;

s3 = s1 + kr1*MOLECON(Pollut, IDX(NO))/s2;

rf = 2.*o_eq*kf1*MOLECON(Pollut, N2)*s1/s3;

rr = 2.*o_eq*kr1*MOLECON(Pollut, IDX(NO))*kr2*MOLECON(Pollut, IDX(NO))/s3;


POLLUT_FRATE(Pollut) = rf;
POLLUT_RRATE(Pollut) = rr;


}
    }
    break;

    case OUT_PDF:
    /* Char Contributions, must be included here */
    break;

  default:
  /* Not used */
  break;
  }
 } 
break;

default:

;

}

}