31 #define HELIUM_IONIZATION NO
43 static void GetHFracs (
double *,
double *);
62 double eH, eHe, eHpH, eHplus, eH2,
rhoe, eHeplus, eHeplusplus;
71 printf (
"! InternalEnergy: x or y are too large %f, %f\n",
92 #if HELIUM_IONIZATION == YES
96 eHeplus = eHeplusplus = 0.0;
100 eHpH = eHplus = eHeplus = eHeplusplus = 0.0;
108 rhoe = (eH2 + eH + eHe + eHpH + eHplus + eHeplus + eHeplusplus)*(
CONST_kB*T*v[RHO]/
CONST_mH);
115 void GetMu(
double T,
double rho,
double *mu)
132 #if HELIUM_IONIZATION == YES
157 double rhs1, rhs2, rhs3;
167 b = rhs1*rhs2*sqrt(rhs2)*rhs3;
169 fdeg[
DEG_y] = 2.0/(1.0 + sqrt(1.0 + 4.0/b));
174 b = rhs1*rhs2*sqrt(rhs2)*rhs3;
176 fdeg[
DEG_x] = 2.0/(1.0 + sqrt(1.0 + 4.0/b));
178 #if HELIUM_IONIZATION == YES
182 fdeg[
DEG_z1] = -2.0*c/(b + sqrt(b*b - 4.0*c));
187 fdeg[
DEG_z2] = -2.0*c/(b + sqrt(b*b - 4.0*c));
194 void GetHFracs(
double *
v,
double *f)
209 #elif COOLING == H2_COOL
224 print1 (
"! GetHFracs: not defined for this cooling\n");
267 double cv, mu, chirho = 1.0, chiT = 1.0, rho,
rhoe;
268 double ep, em, delta = 1.e-2;
269 double Tp, Tm, mup, mum, rhop, rhom;
270 double dmu_dT, dmu_drho, de_dT;
297 Tp = T*(1.0 + delta);
298 Tm = T*(1.0 - delta);
302 de_dT = (ep - em)/(2.0*delta*T);
308 GetMu(Tp, rho, &mup);
309 GetMu(Tm, rho, &mum);
310 dmu_dT = (mup - mum)/(2.0*delta*T);
311 chiT = 1.0 - T/mu*dmu_dT;
313 rhop = rho*(1.0 + delta);
314 rhom = rho*(1.0 - delta);
315 GetMu(T, rhop, &mup);
316 GetMu(T, rhom, &mum);
317 dmu_drho = (mup - mum)/(2.0*delta*rho);
318 chirho = 1.0 - rho/mu*dmu_drho;
325 gmm1 = v[PRS]/(cv*T*v[
RHO])*chiT*chiT + chirho;
332 void InternalEnergyBracket (
double rhoe,
double *
v,
double *Tlo,
double *Thi)
347 (*Thi) = q/(a*(1.0 + f[
DEG_x])*f[
DEG_y] + b + c*(1.0 - f[
DEG_y])*1.49);
#define CONST_h
Planck Constant.
void GetMu(double T, double rho, double *mu)
#define UNIT_DENSITY
Unit density in gr/cm^3.
double InternalEnergy(double *v, double T)
static void GetSahaHFracs(double, double, double *)
#define CONST_eV
Electron Volt in erg.
void print1(const char *fmt,...)
void GetFuncDum(double, double *)
#define CONST_me
Electron mass.
#define UNIT_VELOCITY
Unit velocity in cm/sec.
#define CONST_kB
Boltzmann constant.
#define CONST_mH
Hydrogen atom mass.
double MeanMolecularWeight(double *V)
int GetPV_Temperature(double *v, double *T)
#define QUIT_PLUTO(e_code)
double InternalEnergyFunc(double *v, double T)