38 double Y1, Y2,
scrh, th;
39 double W, W2, S2_W2, fW, dfW, dW;
40 double dv2_dW, chi, chi_p_rho;
41 double rho, p, lor, lor2;
42 double dp, dp_drho, dp_dchi;
43 double dchi_dW, drho_dW;
45 double one_m_v2,
vel2, Wp;
63 #if RMHD_REDUCED_ENERGY == YES
64 Y1 = -4.0*(E + D -
B2);
65 Y2 = m2 - 2.0*(E +
D)*B2 + B2*B2;
68 Y2 = m2 - 2.0*E*B2 + B2*
B2;
70 chi = Y1*Y1 - 12.0*Y2;
72 W = ( - Y1 + sqrt(chi))/6.0;
75 #if RMHD_REDUCED_ENERGY == YES
82 #if RMHD_REDUCED_ENERGY == YES
91 vel2 = S2_W2*Y1*(Y1*W + 1.0) + m2*Y2;
92 one_m_v2 = 1.0 -
vel2;
97 print (
"! EnergySolve: |v| = %f > 1 during iter # %d, ",
102 #if RMHD_REDUCED_ENERGY == YES
103 chi = Wp/lor2 - D*vel2/(lor + 1.0);
105 chi = (W - D*lor)*one_m_v2;
108 dv2_dW = -2.0*Y2*(3.0*S2_W2 + Y1*(S2_W2*B2*B2/W + m2));
116 dchi_dW = one_m_v2 - 0.5*lor*(D + 2.0*chi*lor)*dv2_dW;
117 drho_dW = -0.5*D*lor*dv2_dW;
127 chi_p_rho = chi + rho;
128 scrh = sqrt(9.0*chi*chi + 18.0*rho*chi + 25.0*rho*rho);
129 p = 2.0*chi*(chi_p_rho + rho)/(5.0*chi_p_rho + scrh);
131 scrh = 1.0/(5.0*chi_p_rho - 8.0*p);
132 dp_dchi = (2.0*chi_p_rho - 5.0*p)*scrh;
133 dp_drho = (2.0*chi - 5.0*p)*scrh;
139 dp = dp_dchi*dchi_dW + dp_drho*drho_dW;
140 #if RMHD_REDUCED_ENERGY == YES
141 fW = Wp + 0.5*(B2 + (B2*m2 - S2)*Y2) - (E + p);
142 dfW = 1.0 - dp - (B2*m2 - S2)*Y2*Y1;
145 if (fabs(dW) < acc*Wp || fabs(fW) < acc) done = 1;
147 fW = W + 0.5*(B2 + (B2*m2 - S2)*Y2) - (E + p);
148 dfW = 1.0 - dp - (B2*m2 - S2)*Y2*Y1;
151 if (fabs(dW) < acc*W || fabs(fW) < acc) done = 1;
158 print (
"! EnergySolve(): too many iterations, %12.6e %12.6e %12.6e, ",
165 print (
"! EnergySolve(): negative pressure, p = %12.6e, ",p);
172 #if RMHD_REDUCED_ENERGY == YES
189 par->
sigma_c = par->
prs*lor/pow(rho,2.0/3.0)*(1.5*th + sqrt(2.25*th*th + 1.0));
double m2
Square of total momentum (input).
double D
Lab density (input).
int EnergySolve(Map_param *par)
double S2
Square of S (input).
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 B2
Square of magnetic field (input).
double E
Total energy (input).