35 double v2, v2max, v2min, acc = 1.e-13;
38 double lor, lor2, rho, sigma, th;
40 double fv2, dfv2, dv2, dv2old;
52 v2 = 0.5*(v2max+v2min);
56 lor2 = 1.0/(1.0 - v2);
67 th = sigma*x/sqrt(1.0 + 3.0*sigma*x);
68 h = 2.5*th + sqrt(2.25*th*th + 1.0);
71 scrh = -2.0*x/(3.0*lor)*sigma*(2.0*sigma*x - 3.0*th*th);
72 scrh /= 2.0*th*(1.0 + 3.0*sigma*x);
73 dW_dlor = W/lor + D*lor*(5.0*h - 8.0*th)/(2.0*h - 5.0*th)*
scrh;
80 dfv2 = 1.0 + m2/W3*dW_dlor*lor2*lor;
83 if ((((v2-v2max)*dfv2-fv2)*((v2-v2min)*dfv2-fv2) > 0.0)
84 || (fabs(2.*fv2) > fabs(dv2old*dfv2))) {
86 dv2 = 0.5*(v2max-v2min);
88 if (v2min == v2) done = 1;
94 if (temp == v2) done = 1;
97 if (fabs(dv2) < acc*v2 || fabs(fv2) < acc) done = 1;
99 lor2 = 1.0/(1.0 - v2);
110 x = pow(rho,2.0/3.0);
111 th = sigma*x/sqrt(1.0 + 3.0*sigma*x);
112 h = 2.5*th + sqrt(2.25*th*th + 1.0);
115 scrh = -2.0*x/(3.0*lor)*sigma*(2.0*sigma*x - 3.0*th*th);
116 scrh /= 2.0*th*(1.0 + 3.0*sigma*x);
117 dW_dlor = W/lor + D*lor*(5.0*h-8.0*th)/(2.0*h-5.0*th)*
scrh;
127 if (fv2 < 0.0) v2min = v2;
130 dfv2 = 1.0 + m2/W3*dW_dlor*lor2*lor;
134 par->
prs = sigma*x*rho;
141 print (
"! EntropySolve: too many iterations,%d, ",k);
144 par->
E = W - par->
prs;
148 if (par->
prs < 0.0 || sigma < 0.0 || W < 0.0) {
150 print (
"! EntropySolve: negative pressure, p = %12.6e\n",par->
prs);
153 par->
E = W - par->
prs;
157 par->
rho = par->
D/lor;
159 par->
E = W - par->
prs;
double m2
Square of total momentum (input).
double D
Lab density (input).
double g_smallPressure
Small value for pressure fix.
double rho
proper density (output)
int EntropySolve(Map_param *par)
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).