107 int ii,
k, nv, status;
108 double T, mu, t3, rho, prs, fn, gn, hn, x;
109 double N_H, n_el, cr, ci, rlosst, src_pr, crold,crnew;
110 double kr1, kr2, kr3, kr4, em[20], krvalues[9];
112 static int first_call = 1;
122 N_tot = N_H_rho*(1.0 + frac_He +
frac_Z);
132 v[nv] =
MAX(v[nv], 0.0);
133 v[nv] =
MIN(v[nv], 1.0);
144 print1 (
"! Radiat: mu = %f < 0 \n",mu);
170 print1 (
"! Radiat: negative electron density\n");
204 rhs[
X_HI] = Unit_Time*N_H*( gn*(kr2*fn + kr3*gn + kr4*x) + cr*hn*x
205 -fn*(ci*x + kr1*fn));
206 rhs[
X_H2] = 0.5*Unit_Time*N_H*(kr1*fn*fn - gn*(kr2*fn + kr3*gn + kr4*x));
207 rhs[
X_HII] = Unit_Time*N_H*(ci*fn - cr*hn)*x;
218 double logt3 = log(t3);
229 em[3] = em[1] + em[2];
232 em[4] = -16.818342 + logt3*(37.383713 + logt3*(58.145166 + logt3*(48.656103 + logt3*(20.159831 + logt3*(3.8479610)))));
234 em[4] = -24.311209 + logt3*(3.5692468 + logt3*(-11.332860 + logt3*(-27.850082 + logt3*(-21.328264 + logt3*(-4.2519023)))));
237 em[4] = -24.311209 + logt3*(4.6450521 + logt3*(-3.7209846 + logt3*(5.9369081 + logt3*(-5.5108047 + logt3*(1.5538288)))));
240 em[4] =
POW10(em[4]);
242 em[5] = -23.962112 + logt3*(2.09433740 + logt3*(-0.77151436 + logt3*(0.43693353 + logt3*(-0.14913216 + logt3*(-0.033638326)))));
244 em[5] =
POW10(em[5]);
245 em[6] = (fn*em[4] + gn*em[5])*N_H;
249 em[7] = gn*N_H*em[3]/(1.0 + (em[3]/em[6]));
253 em[8] = ci*13.6*1.6e-12*fn;
257 em[9] = cr*0.67*1.6e-12*hn*
T/11590.0;
261 em[10] = 7.18e-12*(kr2*fn*gn + kr3*gn*gn + kr4*gn*(n_el/N_H))*N_H*N_H;
265 em[11] = krvalues[8]*N_H*N_H*fn*fn;
267 em[13] = em[11] + em[10] + em[7] + (em[8] + em[9])*n_el*N_H;
278 rhs[RHOE] = -E_cost*rlosst;
#define UNIT_DENSITY
Unit density in gr/cm^3.
#define T_CUT_RHOE
Sets the lowest cut-off temperature used in the PVTE_LAW equation of state.
double InternalEnergy(double *v, double T)
void print1(const char *fmt,...)
#define CONST_amu
Atomic mass unit.
double g_smallPressure
Small value for pressure fix.
#define UNIT_VELOCITY
Unit velocity in cm/sec.
int GetEV_Temperature(double rhoe, double *v, double *T)
#define UNIT_LENGTH
Unit Length in cm.
#define CONST_AZ
Mean atomic weight of heavy elements.
double g_minCoolingTemp
The minimum temperature (in K) below which cooling is suppressed.
double MeanMolecularWeight(double *V)
#define CONST_AHe
Atomic weight of Helium.
void H2RateTables(double T, double *krvals)
#define QUIT_PLUTO(e_code)