31 double p, D_1,
alpha, alpha2, lor2, lor, m, tol=1.e-14;
32 double tau, theta, h, dh_dp, dh_dtau,
gmmr;
33 double yp, dyp, dp,
scrh;
47 for (iter = 0; iter <
MAX_ITER; iter++) {
51 lor2 = 1.0 - m2/alpha2;
53 lor2 =
MAX(lor2,1.e-9);
65 h = 2.5*theta + sqrt(2.25*theta*theta + 1.0);
66 scrh = (5.0*h - 8.0*theta)/(2.0*h - 5.0*theta);
72 dyp = D*lor*dh_dp - m2*lor2*lor/(alpha2*
alpha)*(lor*dh_dtau + D*h) - 1.0;
75 if (fabs (dp) < tol*p)
break;
84 print(
"! EnergySolve: NaN found while recovering pressure, ");
91 print(
"! EnergySolve: negative pressure (%8.2e), ", p);
96 par->
rho = par->
D/lor;
108 par->
sigma_c = p*lor/pow(par->
rho,2.0/3.0)*(1.5*theta + sqrt(2.25*theta*theta + 1.0);
double m2
Square of total momentum (input).
double D
Lab density (input).
int EnergySolve(Map_param *par)
double rho
proper density (output)
double W
D*h*lor (output).
void print(const char *fmt,...)
double lor
Lorentz factor (output).
double prs
Thermal pressure (output).
double sigma_c
Conserved entropy (input).
double E
Total energy (input).