00001 #include "../copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include <math.h>
00022 #include <stdio.h>
00023 #include <stdlib.h>
00024 #include <float.h>
00025 #include "../defs.h"
00026 #include "../athena.h"
00027 #include "../globals.h"
00028 #include "prototypes.h"
00029 #include "../prototypes.h"
00030
00031 #ifdef EXACT_FLUX
00032 #ifdef SPECIAL_RELATIVITY
00033
00034 #ifdef MHD
00035 #error : The exact flux for MHD has not been implemented.
00036 #endif
00037
00038 #if (NSCALARS > 0)
00039 #error : Passive scalars have not been implemented in the exact flux.
00040 #endif
00041
00042 #define EPS 3.0e-11
00043 #define JMAX 40
00044 enum WaveType {Two_S, RS, SR, Two_R, None};
00045
00046 static enum WaveType wave;
00047
00048 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00049 const Prim1DS Wl, const Prim1DS Wr, const Real Bx, Cons1DS *pF);
00050
00051 void gauleg(double x1, double x2, Real x[], Real w[], int n);
00052
00053 Real rtbis_vel(Real (*func)(const Prim1DS Wl, const Prim1DS Wr, Real p), const Prim1DS Wl,
00054 const Prim1DS Wr, Real x1, Real x2, Real xacc);
00055
00056 Real getDelVRel(const Prim1DS Wl, const Prim1DS Wr, Real p);
00057
00058 Real getVRel_2R(const Prim1DS Wl, const Prim1DS Wr, Real p);
00059
00060 Real getVRel_RS(const Prim1DS Wl, const Prim1DS Wr, Real p);
00061
00062 Real getVRel_2S(const Prim1DS Wl, const Prim1DS Wr, Real p);
00063
00064 Real getVlim_2S(const Prim1DS Wl, const Prim1DS Wr);
00065
00066 Real getVlim_RS(const Prim1DS Wl, const Prim1DS Wr);
00067
00068 Real getVlim_2R(const Prim1DS Wl, const Prim1DS Wr);
00069
00070 Real integrateRaref(const Prim1DS W, Real a, Real b);
00071
00072 Real getP(const Prim1DS Wl, const Prim1DS Wr);
00073
00074 void getShockVars(const Prim1DS Wa, Real Pb, char dir, Real *pJ, Real *pV_shock, Real *pd);
00075
00076 Real getVb_Shock(const Prim1DS Wa, Real Pb, char dir);
00077
00078 Real getVb_Raref(const Prim1DS Wa, Real Pb, char dir);
00079
00080 Real getXi(const Prim1DS Wa, Real P, Real vxc, char dir);
00081
00082 Real rtbis_xi(Real (*func)(const Prim1DS Wa, Real P, Real vx, char), const Prim1DS Wa,
00083 char dir, Real f_x1, Real f_x2, Real x1, Real x2, Real xacc);
00084
00085 void getVelT_Raref(const Prim1DS Wa, Real P, Real vxb, Real *pVy, Real *pVz);
00086
00087 void getVelT_Shock(const Prim1DS Wa, Real P, Real vxb, Real *pVy, Real *pVz);
00088
00089 void setFluxes(Real vx, Real vy, Real vz, Real P, Real d,Cons1DS *pF );
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00104 const Prim1DS Wl, const Prim1DS Wr, const Real Bx, Cons1DS *pF)
00105 {
00106 Real pc, vxc, dcl, dcr;
00107 Real vl_shock, vr_shock;
00108 Real vr_hd, vr_tl, vl_hd, vl_tl;
00109 Real eps, xacc;
00110 Real vt, p, vx, vy, vz, C, d;
00111 Real signY = 1.0;
00112 Real signZ = 1.0;
00113 Real TOL = 1e-5;
00114
00115
00116 if ((fabs(Wl.P-Wr.P) <= TOL) && (fabs(Wl.Vx-Wr.Vx) <= TOL)) {
00117 if (vxc >= 0.0) {
00118 setFluxes(Wl.Vx, Wl.Vy, Wl.Vz, Wl.P, Wl.d, pF);
00119 return;
00120 }
00121 else {
00122 setFluxes(Wl.Vx, Wr.Vy, Wr.Vz, Wl.P, Wr.d, pF);
00123 return;
00124 }
00125 }
00126
00127
00128 pc = getP(Wl, Wr);
00129
00130
00131 if (wave == RS) {
00132 vxc = getVb_Shock(Wr, pc, 'R');
00133 getShockVars(Wr, pc, 'R', NULL, &vr_shock, &dcr);
00134 dcl = Wl.d*pow((double)(pc/Wl.P), (double) (1.0/Gamma));
00135 }
00136 else if (wave == SR) {
00137 vxc = getVb_Shock(Wl, pc, 'L');
00138 getShockVars(Wl, pc, 'L', NULL, &vl_shock, &dcl);
00139 dcr = Wr.d*pow((double)(pc/Wr.P), (double) (1.0/Gamma));
00140 }
00141 else if (wave == Two_R) {
00142 vxc = getVb_Raref(Wl, pc, 'L');
00143 dcl = Wl.d*pow((double)(pc/Wl.P), (double) (1.0/Gamma));
00144 dcr = Wr.d*pow((double)(pc/Wr.P), (double) (1.0/Gamma));
00145 }
00146 else if (wave == Two_S) {
00147 vxc = getVb_Shock(Wl, pc, 'L');
00148 getShockVars(Wl, pc, 'L', NULL, &vl_shock, &dcl);
00149 getShockVars(Wr, pc, 'R', NULL, &vr_shock, &dcr);
00150 }
00151
00152
00153
00154
00155
00156 if (pc > Wl.P) {
00157
00158
00159 if (vl_shock >= 0.0) {
00160
00161 setFluxes(Wl.Vx, Wl.Vy, Wl.Vz, Wl.P, Wl.d, pF);
00162 return;
00163 }
00164 }
00165 else {
00166
00167
00168
00169 vl_hd = getXi(Wl, Wl.P, Wl.Vx, 'L');
00170 vl_tl = getXi(Wl, pc, vxc, 'L');
00171
00172 if (vl_hd >= 0.0) {
00173
00174 setFluxes(Wl.Vx, Wl.Vy, Wl.Vz, Wl.P, Wl.d, pF);
00175 return;
00176 }
00177 else if (vl_tl >= 0.0) {
00178
00179
00180
00181 eps = 1.0f;
00182 do eps /= 2.0f;
00183 while ((float)(1.0 + (eps/2.0)) != 1.0);
00184 xacc = eps*0.5*(Wl.P + pc);
00185
00186
00187 p = rtbis_xi(getXi, Wl, 'L', vl_hd, vl_tl, Wl.P, pc, xacc);
00188 vx = getVb_Raref(Wl, p, 'L');
00189 d = Wl.d*pow((double)(p/Wl.P), (double) (1.0/Gamma));
00190
00191 getVelT_Raref(Wl, p, vx, &vy, &vz);
00192 setFluxes(vx, vy, vz, p, d, pF);
00193 return;
00194 }
00195 }
00196
00197 if (pc > Wr.P) {
00198
00199
00200
00201
00202 if (vr_shock <= 0.0) {
00203
00204
00205 setFluxes(Wr.Vx, Wr.Vy, Wr.Vz, Wr.P, Wr.d, pF);
00206 return;
00207 }
00208 }
00209 else {
00210
00211
00212
00213 vr_hd = getXi(Wr, Wr.P, Wr.Vx, 'R');
00214 vr_tl = getXi(Wr, pc, vxc, 'R');
00215
00216 if (vr_hd <= 0.0) {
00217
00218
00219 setFluxes(Wr.Vx, Wr.Vy, Wr.Vz, Wr.P, Wr.d, pF);
00220 return;
00221
00222 }
00223 else if (vr_tl <= 0.0) {
00224
00225
00226
00227 eps = 1.0f;
00228 do eps /= 2.0f;
00229 while ((float)(1.0 + (eps/2.0)) != 1.0);
00230 xacc = eps*0.5*(Wr.P + pc);
00231
00232
00233 p = rtbis_xi(getXi, Wr, 'R', vr_hd, vr_tl, Wr.P, pc, xacc);
00234 vx = getVb_Raref(Wr, p, 'R');
00235 d = Wr.d*pow((double)(p/Wr.P), (double) (1.0/Gamma));
00236
00237 getVelT_Raref(Wr, p, vx, &vy, &vz);
00238 setFluxes(vx, vy, vz, p, d, pF);
00239 return;
00240 }
00241 }
00242
00243
00244
00245
00246
00247 if (vxc >= 0.0) {
00248
00249 if (Wl.P > pc) {
00250
00251 getVelT_Raref(Wl, pc, vxc, &vy, &vz);
00252 setFluxes(vxc, vy, vz, pc, dcl, pF);
00253 return;
00254 }
00255 else {
00256
00257 getVelT_Shock(Wl, pc, vxc, &vy, &vz);
00258 setFluxes(vxc, vy, vz, pc, dcl, pF);
00259 return;
00260 }
00261
00262 }
00263 else {
00264
00265 if (pc > Wr.P) {
00266
00267 getVelT_Shock(Wr, pc, vxc, &vy, &vz);
00268 setFluxes(vxc, vy, vz, pc, dcr, pF);
00269 return;
00270 }
00271 else {
00272
00273
00274 getVelT_Raref(Wr, pc, vxc, &vy, &vz);
00275 setFluxes(vxc, vy, vz, pc, dcr, pF);
00276 return;
00277 }
00278 }
00279 }
00280
00281
00282
00283
00284 Real integrateRaref(const Prim1DS Wa, Real a, Real b)
00285 {
00286
00287 int i;
00288 Real integral, sum, diff, f, xx;
00289 Real * x;
00290 Real * w;
00291 int n = 10;
00292 Real k, h, vt, G, A;
00293 Real dd, ccs2, hh;
00294
00295
00296
00297 x = (Real *)malloc((size_t) ((n+1)*sizeof(Real)));
00298 w = (Real *)malloc((size_t) ((n+1)*sizeof(Real)));
00299
00300
00301 gauleg(-1.0, 1.0, x, w, n);
00302
00303 sum = 0.5 * (b + a);
00304 diff = 0.5 * (b - a);
00305 integral = 0.0;
00306
00307 for (i = 1; i <= n; i++) {
00308 xx = diff*x[i] + sum;
00309
00310 h = 1.0 + Gamma*Wa.P/((Gamma-1.0)*Wa.d);
00311 vt = sqrt((double)(Wa.Vy*Wa.Vy + Wa.Vz*Wa.Vz));
00312
00313 G = 1.0/sqrt((double)(1.0 - Wa.Vx*Wa.Vx - Wa.Vy*Wa.Vy - Wa.Vz*Wa.Vz));
00314 A = h*G*vt;
00315
00316 dd = Wa.d*pow((double)(xx/Wa.P), (double) (1.0/Gamma));
00317 ccs2 = Gamma*(Gamma - 1.)*xx / (Gamma*xx + (Gamma - 1.)*dd);
00318 hh = 1.0 + xx*Gamma/(dd*(Gamma - 1.0));
00319
00320 f = sqrt((double) (hh*hh + A*A*(1.0-ccs2)))/(dd*sqrt((double)ccs2)*(hh*hh + A*A));
00321
00322 integral = integral + diff*w[i]*f;
00323 }
00324
00325 free(x);
00326 free(w);
00327
00328 return integral;
00329 }
00330
00331
00332
00333
00334
00335
00336 void getShockVars(const Prim1DS Wa, Real Pb, char dir, Real *pJ, Real *pV_s,
00337 Real *pD)
00338 {
00339 Real sign, ha, Ga;
00340 Real db, hb;
00341 Real A, B, C,D;
00342 Real d, J, v_s;
00343 Real TOL = 1e-5;
00344
00345 if (dir == 'L') sign = -1.0;
00346 if (dir == 'R') sign = 1.0;
00347
00348
00349 if (fabs(Wa.P - Pb) <= TOL) {
00350 if (pJ != NULL)
00351 *pJ = 0.0;
00352 if (pD != NULL)
00353 *pD = Wa.d;
00354 if (pV_s != NULL)
00355 *pV_s = Wa.Vx;
00356
00357 return;
00358 }
00359
00360 ha = 1.0 + Gamma*Wa.P/((Gamma-1.0)*Wa.d);
00361 Ga = 1.0/sqrt((double)(1.0 - Wa.Vx*Wa.Vx - Wa.Vy*Wa.Vy - Wa.Vz*Wa.Vz));
00362
00363
00364 A = 1.0 + (Gamma - 1.0)*(Wa.P-Pb)/(Gamma*Pb);
00365 B = 1.0 - A;
00366 C = ha*(Wa.P - Pb)/Wa.d - ha*ha;
00367
00368
00369 if (C > (B*B/(4.0*A)))
00370 ath_error("[exact flux]: Unphysical specific enthalpy in intermediate state");
00371
00372 hb = (-B + sqrt((double)(B*B - 4.0*A*C)))/(2.0*A);
00373
00374 d = Gamma*Pb/((Gamma - 1.0)*(hb - 1.0));
00375 if (pD != NULL)
00376 *pD = d;
00377
00378 J = sign*sqrt((double)((Pb - Wa.P)/(ha/Wa.d - hb/d)));
00379 if (pJ != NULL)
00380 *pJ = J;
00381
00382 A = Wa.d*Wa.d*Ga*Ga;
00383 if (pV_s != NULL)
00384 *pV_s = (A*Wa.Vx +sign*fabs((float) J)*
00385 sqrt((double)(J*J + A*(1.0 - Wa.Vx*Wa.Vx))))/(A + J*J);
00386
00387 return;
00388 }
00389
00390
00391
00392
00393 Real getVlim_2R(const Prim1DS Wl, const Prim1DS Wr)
00394 {
00395
00396 Real vxlc, vxrc;
00397 Real integral;
00398 Real vlb, vrb;
00399
00400
00401 vlb = getVb_Raref(Wl, 0.0, 'L');
00402 vxlc = (Wl.Vx - vlb)/(1.0 - Wl.Vx*vlb);
00403
00404
00405 vrb = getVb_Raref(Wr, 0.0, 'R');
00406 vxrc = (Wr.Vx - vrb)/(1.0 - Wr.Vx*vrb);
00407
00408 return (vxlc - vxrc) / (1.0 - vxlc*vxrc);
00409
00410 }
00411
00412
00413
00414
00415 Real getVlim_RS(const Prim1DS Wl, const Prim1DS Wr)
00416 {
00417
00418 Real a, b, integral;
00419 Real vlb, vrb, vlc, vrc;
00420
00421 if (Wl.P < Wr.P) {
00422
00423
00424 vlb = Wl.Vx;
00425 vrb = getVb_Raref(Wr, Wl.P, 'R');
00426 }
00427 else {
00428
00429
00430 vlb = getVb_Raref(Wl, Wr.P, 'L');
00431 vrb = Wr.Vx;
00432 }
00433
00434 vlc = (Wl.Vx - vlb)/(1.0 - Wl.Vx*vlb);
00435 vrc = (Wr.Vx - vrb)/(1.0 - Wr.Vx*vrb);
00436
00437 return (vlc - vrc)/(1.0 - vlc*vrc);
00438 }
00439
00440
00441
00442
00443 Real getVlim_2S(const Prim1DS Wl, const Prim1DS Wr)
00444 {
00445 if (Wl.P > Wr.P) {
00446 Real vlc, vrc, vrb;
00447
00448 vrb = getVb_Shock(Wr, Wl.P, 'R');
00449 vrc = (Wr.Vx - vrb)/(1.0 - vrb*Wr.Vx);
00450 vlc = 0.0;
00451 return (vlc - vrc)/(1.0 - vlc*vrc);
00452 }
00453 else {
00454 Real vlc, vrc, vlb;
00455
00456 vlb = getVb_Shock(Wl, Wr.P, 'L');
00457 vlc = (Wl.Vx - vlb)/(1.0 - vlb*Wl.Vx);
00458 vrc = 0.0;
00459 return (vlc - vrc)/(1.0 - vlc*vrc);
00460
00461 }
00462 }
00463
00464
00465
00466
00467 Real getDelVRel(const Prim1DS Wl, const Prim1DS Wr, Real p)
00468 {
00469 Real vRel;
00470 Real vRel_0 = (Wl.Vx - Wr.Vx)/(1.0 - Wl.Vx*Wr.Vx);
00471
00472 if (wave == Two_R) {
00473 vRel = getVRel_2R(Wl, Wr, p);
00474 return vRel - vRel_0;
00475 }
00476
00477 else if (wave == RS || wave == SR) {
00478 vRel = getVRel_RS(Wl, Wr, p);
00479 return vRel - vRel_0;
00480 }
00481 else if (wave == Two_S) {
00482 vRel = getVRel_2S(Wl, Wr, p);
00483 return vRel - vRel_0;
00484 }
00485
00486 return -1.0;
00487 }
00488
00489
00490
00491
00492 Real getP(const Prim1DS Wl, const Prim1DS Wr)
00493 {
00494 Real vRel, vRS, vSS, vRR;
00495 Real pmin, pmax;
00496 Real xacc;
00497 Real pc;
00498 float tol, eps;
00499
00500
00501 vRel = (Wl.Vx - Wr.Vx) / (1.0 - Wl.Vx*Wr.Vx);
00502
00503
00504
00505 vSS = getVlim_2S(Wl, Wr);
00506 vRS = getVlim_RS(Wl, Wr);
00507
00508 if (vRel <= vRS) {
00509 wave = Two_R;
00510 pmin = 0.0;
00511 if (Wl.P < Wr.P) pmax = Wl.P;
00512 else pmax = Wr.P;
00513 }
00514 else if (vRel > vRS && vRel <= vSS) {
00515 if (Wl.P > Wr.P) {
00516 wave = RS;
00517 pmin = Wr.P;
00518 pmax = Wl.P;
00519 }
00520 else {
00521 wave = SR;
00522 pmin = Wl.P;
00523 pmax = Wr.P;
00524 }
00525 }
00526 else {
00527 wave = Two_S;
00528 if (Wl.P > Wr.P) pmin = Wl.P;
00529 else pmin = Wr.P;
00530 pmax = 1000 * 0.5*(Wl.P + Wr.P);
00531 }
00532
00533
00534 eps = 1.0f;
00535 do eps /= 2.0f;
00536 while ((float)(1.0 + (eps/2.0)) != 1.0);
00537 xacc = eps*0.5*(pmin + pmax);
00538
00539
00540 pc = rtbis_vel(getDelVRel, Wl, Wr, pmin, pmax, xacc);
00541 return pc;
00542
00543 }
00544
00545
00546
00547
00548 Real getVRel_2R(const Prim1DS Wl, const Prim1DS Wr, Real p)
00549 {
00550
00551 Real vlb, vrb;
00552 Real integral;
00553 Real vlc, vrc;
00554
00555
00556 vlb = getVb_Raref(Wl, p, 'L');
00557
00558 vrb = getVb_Raref(Wr, p, 'R');
00559
00560
00561 vlc = (Wl.Vx - vlb) / (1.0 - Wl.Vx*vlb);
00562 vrc = (Wr.Vx - vrb) / (1.0 - Wr.Vx*vrb);
00563
00564 return (vlc - vrc)/(1.0 - vrc*vlc);
00565 }
00566
00567
00568
00569
00570 Real getVRel_RS(const Prim1DS Wl, const Prim1DS Wr, Real p)
00571 {
00572 Real vlb, vrb;
00573 Real vlc, vrc;
00574 Real v_shock, G_shock, J, ha, Ga;
00575 Real num, denom, integral;
00576
00577 if (wave == SR) {
00578
00579
00580
00581 vlb = getVb_Shock(Wl, p, 'L');
00582
00583
00584 vrb = getVb_Raref(Wr, p, 'R');
00585
00586 vlc = (Wl.Vx - vlb) / (1.0 - Wl.Vx*vlb);
00587 vrc = (Wr.Vx - vrb) / (1.0 - Wr.Vx*vrb);
00588
00589 return (vlc - vrc)/(1.0 - vrc*vlc);
00590 }
00591 else {
00592
00593
00594
00595 vrb = getVb_Shock(Wr, p, 'R');
00596
00597
00598 vlb = getVb_Raref(Wl, p, 'L');
00599
00600 vlc = (Wl.Vx - vlb) / (1.0 - Wl.Vx*vlb);
00601 vrc = (Wr.Vx - vrb) / (1.0 - Wr.Vx*vrb);
00602
00603 return (vlc - vrc)/(1.0 - vrc*vlc);
00604 }
00605 }
00606
00607
00608
00609 Real getVRel_2S(const Prim1DS Wl, const Prim1DS Wr, Real p)
00610 {
00611 Real vlb, vrb;
00612 Real vlc, vrc;
00613 char dir;
00614 Real v_shock, G_shock, J, ha, Ga;
00615 Real num, denom;
00616
00617
00618 vlb = getVb_Shock(Wl, p, 'L');
00619
00620
00621 vrb = getVb_Shock(Wr, p, 'R');
00622
00623
00624 vlc = (Wl.Vx - vlb) / (1.0 - Wl.Vx*vlb);
00625 vrc = (Wr.Vx - vrb) / (1.0 - Wr.Vx*vrb);
00626
00627 return (vlc - vrc)/(1.0 - vlc*vrc);
00628 }
00629
00630
00631
00632 Real getVb_Shock(const Prim1DS Wa, Real Pb, char dir)
00633 {
00634 Real num, denom;
00635 Real G_shock, ha, Ga;
00636 Real J, v_shock;
00637
00638 getShockVars(Wa, Pb, dir, &J, &v_shock, NULL);
00639 if (fabs(J) == 0.0) return Wa.Vx;
00640
00641 G_shock = 1.0 / sqrt((double) (1.0 - v_shock*v_shock));
00642
00643
00644 ha = 1.0 + Gamma*Wa.P/((Gamma-1.0)*Wa.d);
00645 Ga = 1.0/sqrt((double)(1.0 - Wa.Vx*Wa.Vx - Wa.Vy*Wa.Vy - Wa.Vz*Wa.Vz));
00646
00647 num = ha*Ga*Wa.Vx + G_shock*(Pb - Wa.P)/J;
00648 denom = ha*Ga + (Pb - Wa.P)*(G_shock*Wa.Vx/J + 1.0/(Wa.d*Ga));
00649
00650 return num/denom;
00651 }
00652
00653
00654
00655 Real getVb_Raref(const Prim1DS Wa, Real Pb, char dir)
00656 {
00657 Real integral, sign, B;
00658
00659 if (dir == 'L') sign = -1.0;
00660 if (dir == 'R') sign = 1.0;
00661
00662 integral = integrateRaref(Wa, Wa.P, Pb);
00663 B = 0.5*log((double)((1. + Wa.Vx)/(1. - Wa.Vx))) + sign*integral;
00664
00665 return tanh(B);
00666 }
00667
00668
00669
00670
00671 Real getXi(const Prim1DS Wa, Real P, Real vxc, char dir)
00672 {
00673 Real vt, h, G, A;
00674
00675 Real dc, hc, vtc, cs, v2, num, denom;
00676 Real sign;
00677
00678 if (dir == 'L') sign = -1.0;
00679 if (dir == 'R') sign = 1.0;
00680
00681 vt = sqrt((double)(Wa.Vy*Wa.Vy + Wa.Vz*Wa.Vz));
00682 h = 1.0 + Gamma*Wa.P/((Gamma-1.0)*Wa.d);
00683 G = 1.0/sqrt((double)(1.0 - Wa.Vx*Wa.Vx - Wa.Vy*Wa.Vy - Wa.Vz*Wa.Vz));
00684 A = h*G*vt;
00685
00686 dc = Wa.d*pow((double)(P/Wa.P), (double) (1.0/Gamma));
00687 hc = 1.0 + Gamma*P/((Gamma-1.0)*dc);
00688 vtc = A*sqrt((double)((1.0 - vxc*vxc)/(hc*hc + A*A)));
00689 cs = sqrt((double)(Gamma*(Gamma-1.0)*P/((Gamma -1.0)*dc + Gamma*P)));
00690
00691 v2 = vxc*vxc + vtc*vtc;
00692 num = vxc*(1-cs*cs) + sign*cs*sqrt((double)((1-v2)*(1.-v2*cs*cs-vxc*vxc*(1.-cs*cs))));
00693 denom = 1.0 - v2*cs*cs;
00694
00695 return num/denom;
00696 }
00697
00698
00699
00700
00701
00702 void getVelT_Raref(const Prim1DS Wa, Real P, Real vxb, Real *pVy, Real *pVz)
00703 {
00704 Real va_t, ha, Ga, A;
00705 Real db, hb, vb_t, Gb, C;
00706 Real signY = 1.0;
00707 Real signZ = 1.0;
00708
00709 if (Wa.Vy == 0.0 && Wa.Vz == 0.0) {
00710 *pVy = 0.0;
00711 *pVz = 0.0;
00712 return;
00713 }
00714
00715 va_t = sqrt((double)(Wa.Vy*Wa.Vy + Wa.Vz*Wa.Vz));
00716 ha = 1.0 + Gamma*Wa.P/((Gamma-1.0)*Wa.d);
00717 Ga = 1.0/sqrt((double)(1.0 - Wa.Vx*Wa.Vx - Wa.Vy*Wa.Vy - Wa.Vz*Wa.Vz));
00718 A = ha*Ga*va_t;
00719
00720 db = Wa.d*pow((double)(P/Wa.P), (double) (1.0/Gamma));
00721 hb = 1.0 + Gamma*P/((Gamma-1.0)*db);
00722
00723 vb_t = A*sqrt((double)((1.0 - vxb*vxb)/(hb*hb + A*A)));
00724 Gb = 1.0/sqrt((double)(1.0 -vxb*vxb - vb_t*vb_t));
00725
00726
00727 if (Wa.Vy == 0.0) {
00728 if (Wa.Vz < 0.0) signZ = -1.0;
00729 *pVy = 0.0;
00730 *pVz = signZ*vb_t;
00731 }
00732 else if (Wa.Vz == 0.0) {
00733 if (Wa.Vy < 0.0) signY = -1.0;
00734 *pVz = 0.0;
00735 *pVy = signY*vb_t;
00736 }
00737 else {
00738 C = (Wa.Vy*Wa.Vy) / (Wa.Vz*Wa.Vz);
00739 if (Wa.Vy < 0.0) signY = -1.0;
00740 if (Wa.Vz < 0.0) signZ = -1.0;
00741
00742 *pVz = signZ*sqrt((double)(vb_t*vb_t / (1.0 + C)));
00743 *pVy = signY*sqrt((double)(vb_t*vb_t / (1.0 + 1.0/C)));
00744 }
00745 }
00746
00747
00748
00749
00750
00751
00752
00753 void getVelT_Shock(const Prim1DS Wa, Real P, Real vxb, Real *pVy, Real *pVz)
00754 {
00755 Real va_t, ha, Ga, Ay, Az;
00756 Real hb, db, vy, vz, Gb;
00757 Real Cy, Cz;
00758 Real D;
00759
00760 if (Wa.Vy == 0.0 && Wa.Vz == 0.0) {
00761 *pVy = 0.0;
00762 *pVz = 0.0;
00763 return;
00764 }
00765
00766 ha = 1.0 + Gamma*Wa.P/((Gamma-1.0)*Wa.d);
00767 Ga = 1.0/sqrt((double)(1.0 - Wa.Vx*Wa.Vx - Wa.Vy*Wa.Vy - Wa.Vz*Wa.Vz));
00768
00769 Ay = ha*Ga*Wa.Vy;
00770 Az = ha*Ga*Wa.Vz;
00771 db = Wa.d*pow((double)(P/Wa.P), (double) (1.0/Gamma));
00772 hb = 1.0 + Gamma*P/((Gamma-1.0)*db);
00773
00774 Cz = Az*Az/(hb*hb+Az*Az);
00775 Cy = Ay*Ay/(hb*hb+Ay*Ay);
00776 D = (1.0 - Cz*Cy);
00777
00778 vz = sqrt((double)(Cz*(1-vxb*vxb)*(1-Cy)/D));
00779 vy = sqrt((double)(Cy*(1-vxb*vxb)*(1-Cz)/D));
00780
00781 Gb = 1.0/sqrt((double)(1.0 - vxb*vxb - vy*vy - vz*vz));
00782
00783 if (Wa.Vy >= 0.0)
00784 *pVy = vy;
00785 else
00786 *pVy = -vy;
00787
00788 if (Wa.Vz >= 0.0)
00789 *pVz = vz;
00790 else
00791 *pVz = -vz;
00792 }
00793
00794
00795
00796 void setFluxes(Real vx, Real vy, Real vz, Real P, Real d, Cons1DS *pF)
00797 {
00798 Real G, h, D, Sx, Sy, Sz;
00799
00800 if ((vx*vx + vy*vy + vz*vz) >= 1.0)
00801 ath_error("[exact_flux]: Superluminal velocities vx = %f, vy = %f, vz = %f\n",
00802 vx, vy, vz);
00803
00804 G = 1.0/sqrt((double)(1.0 - vx*vx - vy*vy - vz*vz));
00805 h = 1.0 + Gamma*P/((Gamma-1.0)*d);
00806 D = d*G;
00807 Sx = d*h*G*G*vx;
00808 Sy = d*h*G*G*vy;
00809 Sz = d*h*G*G*vz;
00810
00811 pF->d = D*vx;
00812 pF->Mx = Sx*vx + P;
00813 pF->My = Sy*vx;
00814 pF->Mz = Sz*vx;
00815 pF->E = Sx;
00816 }
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826 void gauleg(double x1, double x2, Real x[], Real w[], int n)
00827 {
00828 int m,j,i;
00829 Real z1,z,xm,xl,pp,p3,p2,p1;
00830 m=(n+1)/2;
00831 xm=0.5*(x2+x1);
00832 xl=0.5*(x2-x1);
00833
00834 for (i=1; i<=m; i++) {
00835 z=cos(3.141592654*(i-0.25)/(n+0.5));
00836
00837
00838 do {
00839 p1=1.0;
00840 p2=0.0;
00841 for (j=1;j<=n;j++) {
00842 p3=p2;
00843 p2=p1;
00844 p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
00845 }
00846
00847
00848
00849 pp=n*(z*p1-p2)/(z*z-1.0);
00850 z1=z;
00851 z=z1-p1/pp;
00852 } while (fabs(z-z1) > EPS);
00853 x[i]=xm-xl*z;
00854 x[n+1-i]=xm+xl*z;
00855 w[i]=2.0*xl/((1.0-z*z)*pp*pp);
00856 w[n+1-i]=w[i];
00857 }
00858 }
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868 Real rtbis_vel(Real (*func)(const Prim1DS Wl, const Prim1DS Wr, Real p),
00869 const Prim1DS Wl, const Prim1DS Wr, Real x1, Real x2, Real xacc)
00870 {
00871 Real j;
00872 Real dx, f, fmid, xmid, rtb;
00873 Real vRel_0 = (Wl.Vx-Wr.Vx)/(1.0-Wl.Vx*Wr.Vx);
00874
00875 if (wave == RS || wave == SR) {
00876 f = getVlim_RS(Wl, Wr) - vRel_0;
00877 fmid = getVlim_2S(Wl, Wr) - vRel_0;
00878 }
00879
00880 if (wave == Two_R) {
00881 f = getVlim_2R(Wl, Wr) - vRel_0;
00882 fmid = getVlim_RS(Wl, Wr) - vRel_0;
00883 }
00884
00885 if (wave == Two_S) {
00886 f = getVlim_2S(Wl, Wr) - vRel_0;
00887 fmid = (*func)(Wl, Wr, x2);
00888 }
00889 if (f < 0.0) {
00890 dx = x2-x1;
00891 rtb =x1;
00892 }
00893 else {
00894 dx = x1-x2;
00895 rtb = x2;
00896 }
00897
00898 for (j = 1; j <= JMAX; j++) {
00899 dx = dx / 2.0;
00900 xmid = rtb + dx;
00901 fmid = (*func)(Wl, Wr, xmid);
00902 if (fmid <= 0.0) rtb = xmid;
00903 if (fabs(dx) < xacc || fmid == 0.0) {
00904 return rtb;
00905 }
00906 }
00907 }
00908
00909
00910
00911
00912
00913
00914
00915
00916 Real rtbis_xi(Real (*func)(const Prim1DS Wa, Real P, Real vx, char),
00917 const Prim1DS Wa, char dir, Real f_x1, Real f_x2, Real x1,
00918 Real x2, Real xacc)
00919 {
00920 Real j, vx;
00921 Real dx, f, fmid, xmid, rtb;
00922
00923 f = f_x1;
00924 fmid = f_x2;
00925
00926 if (f < 0.0) {
00927 dx = x2-x1;
00928 rtb =x1;
00929 }
00930 else {
00931 dx = x1-x2;
00932 rtb = x2;
00933 }
00934
00935 for (j = 1; j <= JMAX; j++) {
00936 dx = dx / 2.0;
00937 xmid = rtb + dx;
00938 vx = getVb_Raref(Wa, xmid, dir);
00939 fmid = (*func)(Wa, xmid, vx, dir);
00940
00941
00942 if (fmid <= 0.0) rtb = xmid;
00943 if (fabs(dx) < xacc || fmid == 0.0) {
00944 return rtb;
00945 }
00946 }
00947
00948
00949 }
00950 #endif
00951 #endif