00001 #include "../copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include <math.h>
00018 #include <stdio.h>
00019 #include <stdlib.h>
00020 #include "../defs.h"
00021 #include "../athena.h"
00022 #include "../globals.h"
00023 #include "prototypes.h"
00024 #include "../prototypes.h"
00025
00026 #ifdef HLLD_FLUX
00027 #ifdef SPECIAL_RELATIVITY
00028
00029 #ifndef MHD
00030 #error : The HLLD flux only works for mhd.
00031 #endif
00032
00033 #if (NSCALARS > 0)
00034 #error : The SR HLLD flux does not work with passive scalars.
00035 #endif
00036
00037 #if HYDRO
00038 #error : The SR HLLD flux does not work for gas=hydro
00039 #endif
00040
00041 #define MAX_ITER 100
00042
00043 typedef struct CONS_STATE{
00044 Real DN,EN;
00045 Real M1,M2,M3;
00046 Real B1,B2,B3;
00047 } CONS_STATE;
00048
00049 typedef struct RIEMANN_STATE{
00050 int fail;
00051 Real vx, vy, vz;
00052 Real Bx, By, Bz;
00053 Real Kx, Ky, Kz, K2;
00054 Real w, eta, p, rho;
00055 CONS_STATE u, R;
00056 Real S, Sa, sw;
00057 } Riemann_State;
00058
00059 void flux_LR(Cons1DS U, Prim1DS W, Cons1DS *flux, Real Bx, Real* p);
00060
00061 Real Fstar (Riemann_State *PaL, Riemann_State *PaR,
00062 Real *Sc, Real p, const Real Bx);
00063 int GET_RIEMANN_STATE (Riemann_State *Pv, Real p, int side, const Real Bx);
00064 void GET_ASTATE (Riemann_State *Pa, Real p, const Real Bx);
00065 void GET_CSTATE (Riemann_State *PaL, Riemann_State *PaR, Real p,
00066 CONS_STATE *Uc, const Real Bx);
00067 void getPtot (const Real Bx, const Prim1DS W, Real *pt);
00068 void getMaxSignalSpeeds_pluto(const Prim1DS Wl, const Prim1DS Wr,
00069 const Real Bx, Real* low, Real* high);
00070 void getMaxSignalSpeeds_echo(const Prim1DS Wl, const Prim1DS Wr,
00071 const Real Bx, Real* low, Real* high);
00072 void getVChar_echo(const Prim1DS W, const Real Bx, Real* lml, Real* lmr);
00073 void getVChar_pluto(const Prim1DS W, const Real Bx, Real* lml, Real* lmr);
00074
00075
00076
00077
00078
00079 int QUARTIC (Real b, Real c, Real d, Real e, Real z[]);
00080
00081
00082
00083 int CUBIC(Real b, Real c, Real d, Real z[]);
00084
00085 #ifdef MHD
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00099 const Prim1DS Wl, const Prim1DS Wr, const Real Bxi, Cons1DS *pFlux)
00100 {
00101 CONS_STATE Uhll,Fhll,Uc;
00102 Cons1DS Fl,Fr,Usl,Usr,Utmp,Ftmp;
00103 Prim1DS Wsl,Wsr,Whll,Wavg;
00104 Real Sl, Sr, Pl, Pr;
00105 Real Sla, Sra;
00106 Real dS_1,scrh,Sc;
00107 Real Bx,p0,pguess,f0,p,f,dp;
00108 Riemann_State PaL,PaR;
00109 int switch_to_hll,wave_speed_fail,k,ISTEP;
00110 int current_sheet_fail;
00111
00112 wave_speed_fail = 0;
00113 switch_to_hll = 0;
00114
00115
00116
00117
00118 getMaxSignalSpeeds_pluto(Wl,Wr,Bxi,&Sl,&Sr);
00119
00120 if (Sl != Sl) {
00121 wave_speed_fail = 1;
00122 printf("[hllc_sr_mhd]: NaN in Sl %10.4e %10.4e\n",Sl,Sr);
00123 Sl = -1.0;
00124 Sr = 1.0;
00125 }
00126
00127 if (Sr != Sr) {
00128 wave_speed_fail = 1;
00129 printf("[hllc_sr_mhd]: NaN in Sr %10.4e %10.4e\n",Sl,Sr);
00130 Sl = -1.0;
00131 Sr = 1.0;
00132 }
00133
00134 if (Sl < -1.0) {
00135 wave_speed_fail = 1;
00136 printf("[hllc_sr_mhd]: Superluminal Sl %10.4e %10.4e\n",Sl,Sr);
00137 Sl = -1.0;
00138 Sr = 1.0;
00139 }
00140 if (Sr > 1.0) {
00141 wave_speed_fail = 1;
00142 printf("[hllc_sr_mhd]: Superluminal Sr %10.4e %10.4e\n",Sl,Sr);
00143 Sl = -1.0;
00144 Sr = 1.0;
00145 }
00146
00147
00148
00149
00150 if (wave_speed_fail){
00151 getMaxSignalSpeeds_echo (Wl,Wr,Bxi,&Sla,&Sra);
00152
00153 if (Sla != Sla) {
00154 switch_to_hll = 1;
00155 printf("[hllc_sr_mhd]: NaN in Sl %10.4e %10.4e\n",Sl,Sr);
00156 Sla = -1.0;
00157 Sra = 1.0;
00158 }
00159
00160 if (Sra != Sra) {
00161 switch_to_hll = 1;
00162 printf("[hllc_sr_mhd]: NaN in Sr %10.4e %10.4e\n",Sl,Sr);
00163 Sla = -1.0;
00164 Sra = 1.0;
00165 }
00166
00167 if (Sla < -1.0) {
00168 switch_to_hll = 1;
00169 printf("[hllc_sr_mhd]: Superluminal Sl %10.4e %10.4e\n",Sl,Sr);
00170 Sla = -1.0;
00171 Sra = 1.0;
00172 }
00173 if (Sra > 1.0) {
00174 switch_to_hll = 1;
00175 printf("[hllc_sr_mhd]: Superluminal Sr %10.4e %10.4e\n",Sl,Sr);
00176 Sla = -1.0;
00177 Sra = 1.0;
00178 }
00179
00180 Sl = Sla;
00181 Sr = Sra;
00182
00183 }
00184
00185
00186 flux_LR(Ul,Wl,&Fl,Bxi,&Pl);
00187 flux_LR(Ur,Wr,&Fr,Bxi,&Pr);
00188
00189
00190
00191
00192 dS_1 = 1.0/(Sr - Sl);
00193
00194 Uhll.DN = (Sr*Ur.d - Sl*Ul.d + Fl.d - Fr.d ) * dS_1;
00195 Uhll.M1 = (Sr*Ur.Mx - Sl*Ul.Mx + Fl.Mx - Fr.Mx) * dS_1;
00196 Uhll.M2 = (Sr*Ur.My - Sl*Ul.My + Fl.My - Fr.My) * dS_1;
00197 Uhll.M3 = (Sr*Ur.Mz - Sl*Ul.Mz + Fl.Mz - Fr.Mz) * dS_1;
00198 Uhll.EN = (Sr*Ur.E - Sl*Ul.E + Fl.E - Fr.E ) * dS_1;
00199 Uhll.B1 = (Sr* Bxi- Sl* Bxi ) * dS_1;
00200 Uhll.B2 = (Sr*Ur.By - Sl*Ul.By + Fl.By - Fr.By) * dS_1;
00201 Uhll.B3 = (Sr*Ur.Bz - Sl*Ul.Bz + Fl.Bz - Fr.Bz) * dS_1;
00202
00203 Fhll.DN = (Sr*Fl.d - Sl*Fr.d + Sl*Sr*(Ur.d - Ul.d )) * dS_1;
00204 Fhll.M1 = (Sr*Fl.Mx - Sl*Fr.Mx + Sl*Sr*(Ur.Mx - Ul.Mx)) * dS_1;
00205 Fhll.M2 = (Sr*Fl.My - Sl*Fr.My + Sl*Sr*(Ur.My - Ul.My)) * dS_1;
00206 Fhll.M3 = (Sr*Fl.Mz - Sl*Fr.Mz + Sl*Sr*(Ur.Mz - Ul.Mz)) * dS_1;
00207 Fhll.EN = (Sr*Fl.E - Sl*Fr.E + Sl*Sr*(Ur.E - Ul.E )) * dS_1;
00208 Fhll.B1 = ( Sl*Sr*( Bxi- Bxi)) * dS_1;
00209 Fhll.B2 = (Sr*Fl.By - Sl*Fr.By + Sl*Sr*(Ur.By - Ul.By)) * dS_1;
00210 Fhll.B3 = (Sr*Fl.Bz - Sl*Fr.Bz + Sl*Sr*(Ur.Bz - Ul.Bz)) * dS_1;
00211
00212
00213 if (switch_to_hll) {
00214
00215 pFlux->d = Fhll.DN;
00216 pFlux->Mx = Fhll.M1;
00217 pFlux->My = Fhll.M2;
00218 pFlux->Mz = Fhll.M3;
00219 pFlux->E = Fhll.EN;
00220 pFlux->By = Fhll.B2;
00221 pFlux->Bz = Fhll.B3;
00222
00223 return;
00224 }
00225
00226
00227
00228
00229 if(Sl >= 0.0){
00230
00231 pFlux->d = Fl.d;
00232 pFlux->Mx = Fl.Mx;
00233 pFlux->My = Fl.My;
00234 pFlux->Mz = Fl.Mz;
00235 pFlux->E = Fl.E;
00236 pFlux->By = Fl.By;
00237 pFlux->Bz = Fl.Bz;
00238
00239 return;
00240
00241 }
00242 else if(Sr <= 0.0){
00243
00244 pFlux->d = Fr.d;
00245 pFlux->Mx = Fr.Mx;
00246 pFlux->My = Fr.My;
00247 pFlux->Mz = Fr.Mz;
00248 pFlux->E = Fr.E;
00249 pFlux->By = Fr.By;
00250 pFlux->Bz = Fr.Bz;
00251
00252 return;
00253 }
00254 else {
00255
00256 PaL.S = Sl;
00257 PaR.S = Sr;
00258 Bx = Uhll.B1;
00259
00260 PaL.R.DN = Sl*Ul.d - Fl.d;
00261 PaL.R.EN = Sl*Ul.E - Fl.E;
00262 PaL.R.M1 = Sl*Ul.Mx - Fl.Mx;
00263 PaL.R.M2 = Sl*Ul.My - Fl.My;
00264 PaL.R.M3 = Sl*Ul.Mz - Fl.Mz;
00265 PaL.R.B1 = Sl* Bxi ;
00266 PaL.R.B2 = Sl*Ul.By - Fl.By;
00267 PaL.R.B3 = Sl*Ul.Bz - Fl.Bz;
00268
00269 PaR.R.DN = Sr*Ur.d - Fr.d;
00270 PaR.R.EN = Sr*Ur.E - Fr.E;
00271 PaR.R.M1 = Sr*Ur.Mx - Fr.Mx;
00272 PaR.R.M2 = Sr*Ur.My - Fr.My;
00273 PaR.R.M3 = Sr*Ur.Mz - Fr.Mz;
00274 PaR.R.B1 = Sr* Bxi ;
00275 PaR.R.B2 = Sr*Ur.By - Fr.By;
00276 PaR.R.B3 = Sr*Ur.Bz - Fr.Bz;
00277
00278
00279
00280 Utmp.d = Uhll.DN;
00281 Utmp.Mx = Uhll.M1;
00282 Utmp.My = Uhll.M2;
00283 Utmp.Mz = Uhll.M3;
00284 Utmp.E = Uhll.EN;
00285 Utmp.By = Uhll.B2;
00286 Utmp.Bz = Uhll.B3;
00287
00288 Ftmp.d = Fhll.DN;
00289 Ftmp.Mx = Fhll.M1;
00290 Ftmp.My = Fhll.M2;
00291 Ftmp.Mz = Fhll.M3;
00292 Ftmp.E = Fhll.EN;
00293 Ftmp.By = Fhll.B2;
00294 Ftmp.Bz = Fhll.B3;
00295
00296
00297
00298
00299 scrh = MAX(Wl.P, Wr.P);
00300 if (Bx*Bx/scrh < 0.01) {
00301
00302 Real a,b,c;
00303 a = Sr - Sl;
00304 b = PaR.R.EN - PaL.R.EN + Sr*PaL.R.M1 - Sl*PaR.R.M1;
00305 c = PaL.R.M1*PaR.R.EN - PaR.R.M1*PaL.R.EN;
00306 scrh = b*b - 4.0*a*c;
00307 scrh = MAX(scrh,0.0);
00308 p0 = 0.5*(- b + sqrt(scrh))*dS_1;
00309
00310 } else {
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320 Whll = check_Prim1D(&Utmp, &Bxi);
00321 getPtot(Bxi,Whll,&p0);
00322 }
00323
00324 pguess = p0;
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340 pguess = p0;
00341 switch_to_hll = 0;
00342 f0 = Fstar(&PaL, &PaR, &Sc, p0, Bx);
00343 if (f0 != f0 || PaL.fail) switch_to_hll = 1;
00344
00345
00346
00347
00348 k = 0;
00349 if (fabs(f0) > 1.e-12 && !switch_to_hll){
00350 p = 1.025*p0; f = f0;
00351 for (k = 1; k < MAX_ITER; k++){
00352
00353 f = Fstar(&PaL, &PaR, &Sc, p, Bx);
00354 if ( f != f || PaL.fail || (k > 7) ||
00355 (fabs(f) > fabs(f0) && k > 4)) {
00356 switch_to_hll = 1;
00357 break;
00358 }
00359
00360 dp = (p - p0)/(f - f0)*f;
00361
00362 p0 = p; f0 = f;
00363 p -= dp;
00364 if (p < 0.0) p = 1.e-6;
00365 if (fabs(dp) < 1.e-5*p || fabs(f) < 1.e-6) break;
00366 }
00367 }else p = p0;
00368
00369 if (PaL.fail) switch_to_hll = 1;
00370
00371 if (switch_to_hll) {
00372
00373 *pFlux = Ftmp;
00374
00375 return;
00376 }
00377
00378 if (PaL.Sa >= -1.e-6){
00379
00380 GET_ASTATE (&PaL, p, Bx);
00381
00382 pFlux->d = Fl.d + Sl*(PaL.u.DN - Ul.d );
00383 pFlux->Mx = Fl.Mx + Sl*(PaL.u.M1 - Ul.Mx);
00384 pFlux->My = Fl.My + Sl*(PaL.u.M2 - Ul.My);
00385 pFlux->Mz = Fl.Mz + Sl*(PaL.u.M3 - Ul.Mz);
00386 pFlux->E = Fl.E + Sl*(PaL.u.EN - Ul.E );
00387 pFlux->By = Fl.By + Sl*(PaL.u.B2 - Ul.By);
00388 pFlux->Bz = Fl.Bz + Sl*(PaL.u.B3 - Ul.Bz);
00389
00390 if (pFlux->d != pFlux->d || pFlux->E != pFlux->E ||
00391 pFlux->Mx != pFlux->Mx || pFlux->My != pFlux->My ||
00392 pFlux->Mz != pFlux->Mz || pFlux->By != pFlux->By ||
00393 pFlux->Bz != pFlux->Bz) {
00394
00395 pFlux->d = Ftmp.d;
00396 pFlux->Mx = Ftmp.Mx;
00397 pFlux->My = Ftmp.My;
00398 pFlux->Mz = Ftmp.Mz;
00399 pFlux->E = Ftmp.E;
00400 pFlux->By = Ftmp.By;
00401 pFlux->Bz = Ftmp.Bz;
00402
00403 }
00404
00405 return;
00406
00407 }else if (PaR.Sa <= 1.e-6){
00408
00409 GET_ASTATE (&PaR, p, Bx);
00410
00411 pFlux->d = Fr.d + Sr*(PaR.u.DN - Ur.d );
00412 pFlux->Mx = Fr.Mx + Sr*(PaR.u.M1 - Ur.Mx);
00413 pFlux->My = Fr.My + Sr*(PaR.u.M2 - Ur.My);
00414 pFlux->Mz = Fr.Mz + Sr*(PaR.u.M3 - Ur.Mz);
00415 pFlux->E = Fr.E + Sr*(PaR.u.EN - Ur.E );
00416 pFlux->By = Fr.By + Sr*(PaR.u.B2 - Ur.By);
00417 pFlux->Bz = Fr.Bz + Sr*(PaR.u.B3 - Ur.Bz);
00418
00419 if (pFlux->d != pFlux->d || pFlux->E != pFlux->E ||
00420 pFlux->Mx != pFlux->Mx || pFlux->My != pFlux->My ||
00421 pFlux->Mz != pFlux->Mz || pFlux->By != pFlux->By ||
00422 pFlux->Bz != pFlux->Bz) {
00423
00424 pFlux->d = Ftmp.d;
00425 pFlux->Mx = Ftmp.Mx;
00426 pFlux->My = Ftmp.My;
00427 pFlux->Mz = Ftmp.Mz;
00428 pFlux->E = Ftmp.E;
00429 pFlux->By = Ftmp.By;
00430 pFlux->Bz = Ftmp.Bz;
00431
00432 }
00433
00434 return;
00435
00436 }else{
00437
00438 GET_CSTATE (&PaL, &PaR, p, &Uc, Bx);
00439
00440 if (Sc > 0.0){
00441
00442 pFlux->d = Fl.d + Sl*(PaL.u.DN - Ul.d ) + PaL.Sa*(Uc.DN - PaL.u.DN);
00443 pFlux->Mx = Fl.Mx + Sl*(PaL.u.M1 - Ul.Mx) + PaL.Sa*(Uc.M1 - PaL.u.M1);
00444 pFlux->My = Fl.My + Sl*(PaL.u.M2 - Ul.My) + PaL.Sa*(Uc.M2 - PaL.u.M2);
00445 pFlux->Mz = Fl.Mz + Sl*(PaL.u.M3 - Ul.Mz) + PaL.Sa*(Uc.M3 - PaL.u.M3);
00446 pFlux->E = Fl.E + Sl*(PaL.u.EN - Ul.E ) + PaL.Sa*(Uc.EN - PaL.u.EN);
00447 pFlux->By = Fl.By + Sl*(PaL.u.B2 - Ul.By) + PaL.Sa*(Uc.B2 - PaL.u.B2);
00448 pFlux->Bz = Fl.Bz + Sl*(PaL.u.B3 - Ul.Bz) + PaL.Sa*(Uc.B3 - PaL.u.B3);
00449
00450 if (pFlux->d != pFlux->d || pFlux->E != pFlux->E ||
00451 pFlux->Mx != pFlux->Mx || pFlux->My != pFlux->My ||
00452 pFlux->Mz != pFlux->Mz || pFlux->By != pFlux->By ||
00453 pFlux->Bz != pFlux->Bz) {
00454
00455 pFlux->d = Ftmp.d;
00456 pFlux->Mx = Ftmp.Mx;
00457 pFlux->My = Ftmp.My;
00458 pFlux->Mz = Ftmp.Mz;
00459 pFlux->E = Ftmp.E;
00460 pFlux->By = Ftmp.By;
00461 pFlux->Bz = Ftmp.Bz;
00462
00463 }
00464
00465 return;
00466
00467
00468 }else{
00469
00470 pFlux->d = Fr.d + Sr*(PaR.u.DN - Ur.d ) + PaR.Sa*(Uc.DN - PaR.u.DN);
00471 pFlux->Mx = Fr.Mx + Sr*(PaR.u.M1 - Ur.Mx) + PaR.Sa*(Uc.M1 - PaR.u.M1);
00472 pFlux->My = Fr.My + Sr*(PaR.u.M2 - Ur.My) + PaR.Sa*(Uc.M2 - PaR.u.M2);
00473 pFlux->Mz = Fr.Mz + Sr*(PaR.u.M3 - Ur.Mz) + PaR.Sa*(Uc.M3 - PaR.u.M3);
00474 pFlux->E = Fr.E + Sr*(PaR.u.EN - Ur.E ) + PaR.Sa*(Uc.EN - PaR.u.EN);
00475 pFlux->By = Fr.By + Sr*(PaR.u.B2 - Ur.By) + PaR.Sa*(Uc.B2 - PaR.u.B2);
00476 pFlux->Bz = Fr.Bz + Sr*(PaR.u.B3 - Ur.Bz) + PaR.Sa*(Uc.B3 - PaR.u.B3);
00477
00478 if (pFlux->d != pFlux->d || pFlux->E != pFlux->E ||
00479 pFlux->Mx != pFlux->Mx || pFlux->My != pFlux->My ||
00480 pFlux->Mz != pFlux->Mz || pFlux->By != pFlux->By ||
00481 pFlux->Bz != pFlux->Bz) {
00482
00483 pFlux->d = Ftmp.d;
00484 pFlux->Mx = Ftmp.Mx;
00485 pFlux->My = Ftmp.My;
00486 pFlux->Mz = Ftmp.Mz;
00487 pFlux->E = Ftmp.E;
00488 pFlux->By = Ftmp.By;
00489 pFlux->Bz = Ftmp.Bz;
00490
00491 }
00492
00493 return;
00494
00495 }
00496 }
00497 }
00498 }
00499
00500
00501
00502
00503
00504
00505 Real Fstar (Riemann_State *PaL, Riemann_State *PaR,
00506 Real *Sc, Real p, const Real Bx)
00507 {
00508 int success = 1;
00509 Real dK, Bxc, Byc, Bzc;
00510 Real sBx, fun;
00511 Real vxcL, KLBc;
00512 Real vxcR, KRBc;
00513
00514 sBx = Bx > 0.0 ? 1.0:-1.0;
00515
00516 success *= GET_RIEMANN_STATE (PaL, p, -1, Bx);
00517 success *= GET_RIEMANN_STATE (PaR, p, 1, Bx);
00518
00519
00520
00521 dK = PaR->Kx - PaL->Kx + 1.e-12;
00522
00523 Bxc = Bx*dK;
00524 Byc = PaR->By*(PaR->Kx - PaR->vx)
00525 - PaL->By*(PaL->Kx - PaL->vx)
00526 + Bx*(PaR->vy - PaL->vy);
00527 Bzc = PaR->Bz*(PaR->Kx - PaR->vx)
00528 - PaL->Bz*(PaL->Kx - PaL->vx)
00529 + Bx*(PaR->vz - PaL->vz);
00530
00531 KLBc = PaL->Kx*Bxc + PaL->Ky*Byc + PaL->Kz*Bzc;
00532 KRBc = PaR->Kx*Bxc + PaR->Ky*Byc + PaR->Kz*Bzc;
00533
00534 vxcL = PaL->Kx - dK*Bx*(1.0 - PaL->K2)/(PaL->sw*dK - KLBc);
00535 vxcR = PaR->Kx - dK*Bx*(1.0 - PaR->K2)/(PaR->sw*dK - KRBc);
00536
00537 PaL->Sa = PaL->Kx;
00538 PaR->Sa = PaR->Kx;
00539 *Sc = 0.5*(vxcL + vxcR);
00540 fun = vxcL - vxcR;
00541
00542
00543
00544
00545
00546
00547
00548
00549 success = (vxcL - PaL->Kx) > -1.e-6;
00550 success *= (PaR->Kx - vxcR) > -1.e-6;
00551
00552 success *= (PaL->S - PaL->vx) < 0.0;
00553 success *= (PaR->S - PaR->vx) > 0.0;
00554
00555 success *= (PaR->w - p) > 0.0;
00556 success *= (PaL->w - p) > 0.0;
00557 success *= (PaL->Sa - PaL->S) > -1.e-6;
00558 success *= (PaR->S - PaR->Sa) > -1.e-6;
00559
00560 PaL->fail = !success;
00561
00562 return (fun);
00563 }
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574 int GET_RIEMANN_STATE (Riemann_State *Pv, Real p, int side, const Real Bx)
00575 {
00576 double A, C, G, X, s;
00577 double vx, vy, vz, scrh;
00578
00579 A = Pv->R.M1 + p*(1.0 - Pv->S*Pv->S) - Pv->S*Pv->R.EN;
00580 C = Pv->R.B2*Pv->R.M2 + Pv->R.B3*Pv->R.M3;
00581 G = Pv->R.B2*Pv->R.B2 + Pv->R.B3*Pv->R.B3;
00582 X = Bx*(A*Pv->S*Bx + C) - (A + G)*(Pv->S*p + Pv->R.EN);
00583
00584 vx = ( Bx*(A*Bx + C*Pv->S) - (Pv->R.M1 + p)*(G + A) );
00585 vy = ( - (A + G - Bx*Bx*(1.0 - Pv->S*Pv->S))*Pv->R.M2
00586 + Pv->R.B2*(C + Bx*(Pv->S*Pv->R.M1 - Pv->R.EN)) );
00587 vz = ( - (A + G - Bx*Bx*(1.0 - Pv->S*Pv->S))*Pv->R.M3
00588 + Pv->R.B3*(C + Bx*(Pv->S*Pv->R.M1 - Pv->R.EN)) );
00589
00590 scrh = vx*Pv->R.M1 + vy*Pv->R.M2 + vz*Pv->R.M3;
00591 scrh = X*Pv->R.EN - scrh;
00592 Pv->w = p + scrh/(X*Pv->S - vx);
00593
00594 if (Pv->w < 0.0) return(0);
00595
00596 Pv->vx = vx/X;
00597 Pv->vy = vy/X;
00598 Pv->vz = vz/X;
00599
00600 Pv->Bx = Bx;
00601 Pv->By = -(Pv->R.B2*(Pv->S*p + Pv->R.EN) - Bx*Pv->R.M2)/A;
00602 Pv->Bz = -(Pv->R.B3*(Pv->S*p + Pv->R.EN) - Bx*Pv->R.M3)/A;
00603
00604 s = Bx > 0.0 ? 1.0:-1.0;
00605 if (side < 0) s *= -1.0;
00606
00607 Pv->sw = s*sqrt(Pv->w);
00608
00609 scrh = 1.0/(Pv->S*p + Pv->R.EN + Bx*Pv->sw);
00610 Pv->Kx = scrh*(Pv->R.M1 + p + Pv->R.B1*Pv->sw);
00611 Pv->Ky = scrh*(Pv->R.M2 + Pv->R.B2*Pv->sw);
00612 Pv->Kz = scrh*(Pv->R.M3 + Pv->R.B3*Pv->sw);
00613
00614 Pv->K2 = Pv->Kx*Pv->Kx + Pv->Ky*Pv->Ky + Pv->Kz*Pv->Kz;
00615 return(1);
00616 }
00617
00618
00619
00620
00621 void GET_ASTATE (Riemann_State *Pa, Real p, const Real Bx)
00622 {
00623 Real vB;
00624 Real scrh;
00625
00626 scrh = 1.0/(Pa->S - Pa->vx);
00627
00628 Pa->u.DN = Pa->R.DN*scrh;
00629 Pa->u.B1 = Bx;
00630 Pa->u.B2 = (Pa->R.B2 - Bx*Pa->vy)*scrh;
00631 Pa->u.B3 = (Pa->R.B3 - Bx*Pa->vz)*scrh;
00632
00633 vB = Pa->vx*Pa->u.B1 + Pa->vy*Pa->u.B2 + Pa->vz*Pa->u.B3;
00634 Pa->u.EN = (Pa->R.EN + p*Pa->vx - vB*Bx)*scrh;
00635
00636 Pa->u.M1 = (Pa->u.EN + p)*Pa->vx - vB*Pa->u.B1;
00637 Pa->u.M2 = (Pa->u.EN + p)*Pa->vy - vB*Pa->u.B2;
00638 Pa->u.M3 = (Pa->u.EN + p)*Pa->vz - vB*Pa->u.B3;
00639 }
00640
00641
00642
00643
00644
00645
00646 void GET_CSTATE (Riemann_State *PaL, Riemann_State *PaR, Real p,
00647 CONS_STATE *Uc, const Real Bx)
00648 {
00649 CONS_STATE ua;
00650 double dK;
00651 double vxcL, vycL, vzcL, KLBc;
00652 double vxcR, vycR, vzcR, KRBc;
00653 double vxc, vyc, vzc, vBc;
00654 double Bxc, Byc, Bzc, Sa, vxa;
00655 double scrhL, scrhR;
00656
00657 GET_ASTATE (PaL, p, Bx);
00658 GET_ASTATE (PaR, p, Bx);
00659 dK = (PaR->Kx - PaL->Kx) + 1.e-12;
00660
00661 Bxc = Bx*dK;
00662 Byc = PaR->By*(PaR->Kx - PaR->vx)
00663 - PaL->By*(PaL->Kx - PaL->vx)
00664 + Bx*(PaR->vy - PaL->vy);
00665 Bzc = PaR->Bz*(PaR->Kx - PaR->vx)
00666 - PaL->Bz*(PaL->Kx - PaL->vx)
00667 + Bx*(PaR->vz - PaL->vz);
00668
00669 Bxc = Bx;
00670 Byc /= dK;
00671 Bzc /= dK;
00672
00673 Uc->B1 = Bxc;
00674 Uc->B2 = Byc;
00675 Uc->B3 = Bzc;
00676
00677 KLBc = PaL->Kx*Bxc + PaL->Ky*Byc + PaL->Kz*Bzc;
00678 KRBc = PaR->Kx*Bxc + PaR->Ky*Byc + PaR->Kz*Bzc;
00679
00680 scrhL = (1.0 - PaL->K2)/(PaL->sw - KLBc);
00681 scrhR = (1.0 - PaR->K2)/(PaR->sw - KRBc);
00682
00683 vxcL = PaL->Kx - Uc->B1*scrhL;
00684 vxcR = PaR->Kx - Uc->B1*scrhR;
00685
00686 vycL = PaL->Ky - Uc->B2*scrhL;
00687 vycR = PaR->Ky - Uc->B2*scrhR;
00688
00689 vzcL = PaL->Kz - Uc->B3*scrhL;
00690 vzcR = PaR->Kz - Uc->B3*scrhR;
00691
00692 vxc = 0.5*(vxcL + vxcR);
00693 vyc = 0.5*(vycL + vycR);
00694 vzc = 0.5*(vzcL + vzcR);
00695
00696 if (vxc > 0.0) {
00697 GET_ASTATE (PaL, p, Bx);
00698 ua = PaL->u;
00699 Sa = PaL->Sa;
00700 vxa = PaL->vx;
00701 } else {
00702 GET_ASTATE (PaR, p, Bx);
00703 ua = PaR->u;
00704 Sa = PaR->Sa;
00705 vxa = PaR->vx;
00706 }
00707
00708 vBc = vxc*Uc->B1 + vyc*Uc->B2 + vzc*Uc->B3;
00709
00710 Uc->DN = ua.DN*(Sa - vxa)/(Sa - vxc);
00711 Uc->EN = (Sa*ua.EN - ua.M1 + p*vxc - vBc*Bx)/(Sa - vxc);
00712
00713 Uc->M1 = (Uc->EN + p)*vxc - vBc*Bx;
00714 Uc->M2 = (Uc->EN + p)*vyc - vBc*Uc->B2;
00715 Uc->M3 = (Uc->EN + p)*vzc - vBc*Uc->B3;
00716 }
00717
00718
00719
00720
00721
00722
00723 void getPtot (const Real Bx, const Prim1DS W, Real *pt)
00724 {
00725 Real vel2, Bmag2, vB;
00726
00727 vel2 = SQR(W.Vx) + SQR(W.Vy) + SQR(W.Vz);
00728 Bmag2 = SQR(Bx) + SQR(W.By) + SQR(W.Bz);
00729 vB = W.Vx*Bx + W.Vy*W.By + W.Vz*W.Bz;
00730
00731 *pt = W.P + 0.5*(Bmag2*(1.0 - vel2) + vB*vB);
00732 }
00733
00734
00735
00736
00737 void entropy_flux (const Cons1DS Ul, const Cons1DS Ur,
00738 const Prim1DS Wl, const Prim1DS Wr, const Real Bx, Real *pFlux)
00739 {
00740 Real Fl, Fr;
00741 Real USl, USr;
00742 Real WSl, WSr;
00743 Real Uhll, Fhll;
00744 Real Pl, Pr;
00745 Real Sl, Sr;
00746 Real Sla, Sra;
00747 Real dS_1;
00748 int wave_speed_fail;
00749
00750 wave_speed_fail = 0;
00751
00752
00753
00754
00755 getMaxSignalSpeeds_pluto(Wl,Wr,Bx,&Sl,&Sr);
00756
00757 if (Sl != Sl) {
00758 wave_speed_fail = 1;
00759 printf("[hlle_sr_mhd]: NaN in Sl %10.4e %10.4e\n",Sl,Sr);
00760 Sl = -1.0;
00761 Sr = 1.0;
00762 }
00763
00764 if (Sr != Sr) {
00765 wave_speed_fail = 1;
00766 printf("[hlle_sr_mhd]: NaN in Sr %10.4e %10.4e\n",Sl,Sr);
00767 Sl = -1.0;
00768 Sr = 1.0;
00769 }
00770
00771 if (Sl < -1.0) {
00772 wave_speed_fail = 1;
00773 printf("[hlle_sr_mhd]: Superluminal Sl %10.4e %10.4e\n",Sl,Sr);
00774 Sl = -1.0;
00775 Sr = 1.0;
00776 }
00777 if (Sr > 1.0) {
00778 wave_speed_fail = 1;
00779 printf("[hlle_sr_mhd]: Superluminal Sr %10.4e %10.4e\n",Sl,Sr);
00780 Sl = -1.0;
00781 Sr = 1.0;
00782 }
00783
00784
00785
00786
00787 if (wave_speed_fail){
00788 getMaxSignalSpeeds_echo (Wl,Wr,Bx,&Sla,&Sra);
00789
00790 if (Sla != Sla) {
00791 printf("[hlle_sr_mhd]: NaN in Sl %10.4e %10.4e\n",Sl,Sr);
00792 Sla = -1.0;
00793 Sra = 1.0;
00794 }
00795
00796 if (Sra != Sra) {
00797 printf("[hlle_sr_mhd]: NaN in Sr %10.4e %10.4e\n",Sl,Sr);
00798 Sla = -1.0;
00799 Sra = 1.0;
00800 }
00801
00802 if (Sla < -1.0) {
00803 printf("[hlle_sr_mhd]: Superluminal Sl %10.4e %10.4e\n",Sl,Sr);
00804 Sla = -1.0;
00805 Sra = 1.0;
00806 }
00807 if (Sra > 1.0) {
00808 printf("[hlle_sr_mhd]: Superluminal Sr %10.4e %10.4e\n",Sl,Sr);
00809 Sla = -1.0;
00810 Sra = 1.0;
00811 }
00812
00813 Sl = Sla;
00814 Sr = Sra;
00815
00816 }
00817
00818
00819 WSl = Wl.P*pow(Wl.d,1.0-Gamma);
00820 WSr = Wr.P*pow(Wr.d,1.0-Gamma);
00821 USl = WSl * Ul.d/Wl.d;
00822 USr = WSr * Ur.d/Wr.d;
00823 Fl = USl * Wl.Vx;
00824 Fr = USr * Wr.Vx;
00825
00826 if(Sl >= 0.0){
00827 *pFlux = Fl;
00828 return;
00829 }
00830 else if(Sr <= 0.0){
00831 *pFlux = Fr;
00832 return;
00833 }
00834 else{
00835
00836 dS_1 = 1.0/(Sr - Sl);
00837 *pFlux = (Sr*Fl - Sl*Fr + Sl*Sr*(USr - USl)) * dS_1;
00838 return;
00839 }
00840 }
00841
00842 void flux_LR(Cons1DS U, Prim1DS W, Cons1DS *flux, Real Bx, Real* p)
00843 {
00844 Real wtg2, pt, g, g2, g_1,g_2, h, gmmr, theta;
00845 Real bx, by, bz, vB, b2, Bmag2;
00846
00847
00848
00849 theta = W.P/W.d;
00850 gmmr = Gamma / Gamma_1;
00851
00852 h = 1.0 + gmmr*theta;
00853
00854
00855 g = U.d/W.d;
00856 g2 = SQR(g);
00857 g_1 = 1.0/g;
00858 g_2 = 1.0/g2;
00859
00860 pt = W.P;
00861 wtg2 = W.d*h*g2;
00862
00863 vB = W.Vx*Bx + W.Vy*W.By + W.Vz*W.Bz;
00864 Bmag2 = SQR(Bx) + SQR(W.By) + SQR(W.Bz);
00865
00866 bx = g*( Bx*g_2 + vB*W.Vx);
00867 by = g*(W.By*g_2 + vB*W.Vy);
00868 bz = g*(W.Bz*g_2 + vB*W.Vz);
00869
00870 b2 = Bmag2*g_2 + vB*vB;
00871
00872 pt += 0.5*b2;
00873 wtg2 += b2*g2;
00874
00875 flux->d = U.d*W.Vx;
00876 flux->Mx = wtg2*W.Vx*W.Vx + pt;
00877 flux->My = wtg2*W.Vy*W.Vx;
00878 flux->Mz = wtg2*W.Vz*W.Vx;
00879 flux->E = U.Mx;
00880
00881 flux->Mx -= bx*bx;
00882 flux->My -= by*bx;
00883 flux->Mz -= bz*bx;
00884 flux->By = W.Vx*W.By - Bx*W.Vy;
00885 flux->Bz = W.Vx*W.Bz - Bx*W.Vz;
00886
00887 *p = pt;
00888 }
00889
00890 void getMaxSignalSpeeds_pluto(const Prim1DS Wl, const Prim1DS Wr,
00891 const Real Bx, Real* low, Real* high)
00892 {
00893
00894 Real lml,lmr;
00895 Real lpl,lpr;
00896 Real al,ar;
00897
00898 getVChar_pluto(Wl,Bx,&lml,&lpl);
00899 getVChar_pluto(Wr,Bx,&lmr,&lpr);
00900
00901 *low = MIN(lml, lmr);
00902 *high = MAX(lpl, lpr);
00903 }
00904
00905 void getVChar_pluto(const Prim1DS W, const Real Bx, Real* lm, Real* lp)
00906 {
00907 Real rhoh,vsq,bsq;
00908 Real cssq,vasq,asq;
00909 Real Vx2,gamma,gamma2;
00910 Real Bx2,Bsq,vDotB,vDotBsq,b0,bx;
00911 Real w_1,a0,a1,a2,a3,a4,Q;
00912 Real scrh,scrh1,scrh2,eps2;
00913 Real a2_w,one_m_eps2,lambda[5];
00914 int iflag;
00915
00916 rhoh = W.d + (Gamma/Gamma_1) * (W.P);
00917
00918 Vx2 = SQR(W.Vx);
00919 vsq = Vx2 + SQR(W.Vy) + SQR(W.Vz);
00920 if (vsq > 1.0){
00921
00922 *lm = -1.0;
00923 *lp = 1.0;
00924 return;
00925 }
00926 gamma = 1.0 / sqrt(1 - vsq);
00927 gamma2 = SQR(gamma);
00928
00929 Bx2 = SQR(Bx);
00930 Bsq = Bx2 + SQR(W.By) + SQR(W.Bz);
00931 vDotB = W.Vx*Bx + W.Vy*W.By + W.Vz*W.Bz;
00932 vDotBsq = SQR(vDotB);
00933 b0 = gamma * vDotB;
00934 bx = Bx/gamma2 + W.Vx*vDotB;
00935 bsq = Bsq / gamma2 + SQR(vDotB);
00936
00937 cssq = (Gamma * W.P) / (rhoh);
00938 vasq = bsq / (rhoh + bsq);
00939
00940 if (bsq < 0.0) bsq = 0.0;
00941 if (cssq < 0.0) cssq = 0.0;
00942 if (cssq > 1.0) cssq = 1.0;
00943 if (vasq > 1.0) bsq = rhoh + bsq;
00944
00945 if (vsq < 1.0e-12) {
00946 w_1 = 1.0/(rhoh + bsq);
00947 eps2 = cssq + bsq*w_1*(1.0 - cssq);
00948 a0 = cssq*Bx*Bx*w_1;
00949 a1 = - a0 - eps2;
00950 scrh = a1*a1 - 4.0*a0;
00951 if (scrh < 0.0) scrh = 0.0;
00952
00953 scrh = sqrt(0.5*(-a1 + sqrt(scrh)));
00954 *lp = scrh;
00955 *lm = -scrh;
00956 return;
00957 }
00958
00959 w_1 = 1.0/(rhoh + bsq);
00960
00961 if (Bx < 1.0e-14) {
00962
00963 eps2 = cssq + bsq*w_1*(1.0 - cssq);
00964
00965 scrh1 = (1.0 - eps2)*gamma2;
00966 scrh2 = cssq*vDotBsq*w_1 - eps2;
00967
00968 a2 = scrh1 - scrh2;
00969 a1 = -2.0*W.Vx*scrh1;
00970 a0 = Vx2*scrh1 + scrh2;
00971
00972 *lp = 0.5*(-a1 + sqrt(a1*a1 - 4.0*a2*a0))/a2;
00973 *lm = 0.5*(-a1 - sqrt(a1*a1 - 4.0*a2*a0))/a2;
00974
00975 return;
00976 }
00977
00978 scrh1 = bx;
00979 scrh2 = scrh1*scrh1;
00980
00981 a2_w = cssq*w_1;
00982 eps2 = (cssq*rhoh + bsq)*w_1;
00983 one_m_eps2 = gamma2*rhoh*(1.0 - cssq)*w_1;
00984
00985
00986
00987
00988
00989 scrh = 2.0*(a2_w*vDotB*scrh1 - eps2*W.Vx);
00990 a4 = one_m_eps2 - a2_w*vDotBsq + eps2;
00991 a3 = - 4.0*W.Vx*one_m_eps2 + scrh;
00992 a2 = 6.0*Vx2*one_m_eps2 + a2_w*(vDotBsq - scrh2) + eps2*(Vx2 - 1.0);
00993 a1 = - 4.0*W.Vx*Vx2*one_m_eps2 - scrh;
00994 a0 = Vx2*Vx2*one_m_eps2 + a2_w*scrh2 - eps2*Vx2;
00995
00996 if (a4 < 1.e-12){
00997
00998 printf("[MAX_CH_SPEED]: Can not divide by a4 in MAX_CH_SPEED\n");
00999
01000 *lm = -1.0;
01001 *lp = 1.0;
01002
01003 return;
01004 }
01005
01006 scrh = 1.0/a4;
01007
01008 a3 *= scrh;
01009 a2 *= scrh;
01010 a1 *= scrh;
01011 a0 *= scrh;
01012 iflag = QUARTIC(a3, a2, a1, a0, lambda);
01013
01014 if (iflag){
01015 printf ("Can not find max speed:\n");
01016
01017 printf("QUARTIC: f(x) = %12.6e + x*(%12.6e + x*(%12.6e ",
01018 a0*a4, a1*a4, a2*a4);
01019 printf("+ x*(%12.6e + x*%12.6e)))\n", a3*a4, a4);
01020 printf("[MAX_CH_SPEED]: Failed to find wave speeds");
01021
01022 *lm = -1.0;
01023 *lp = 1.0;
01024
01025 return;
01026 }
01027
01028 *lp = MIN(1.0,MAX(lambda[3], lambda[2]));
01029 *lp = MIN(1.0,MAX(*lp, lambda[1]));
01030 *lp = MIN(1.0,MAX(*lp, lambda[0]));
01031
01032 *lm = MAX(-1.0,MIN(lambda[3], lambda[2]));
01033 *lm = MAX(-1.0,MIN(*lm, lambda[1]));
01034 *lm = MAX(-1.0,MIN(*lm, lambda[0]));
01035
01036 return;
01037
01038 }
01039
01040 void getMaxSignalSpeeds_echo (const Prim1DS Wl, const Prim1DS Wr,
01041 const Real Bx, Real* low, Real* high)
01042 {
01043
01044 Real lml,lmr;
01045 Real lpl,lpr;
01046 Real al,ar;
01047
01048 getVChar_echo(Wl,Bx,&lml,&lpl);
01049 getVChar_echo(Wr,Bx,&lmr,&lpr);
01050
01051 *low = MIN(lml, lmr);
01052 *high = MAX(lpl, lpr);
01053 }
01054
01055 void getVChar_echo(const Prim1DS W, const Real Bx, Real* lm, Real* lp)
01056 {
01057 Real rhoh,vsq,bsq;
01058 Real cssq,vasq,asq;
01059 Real gamma,gamma2;
01060 Real Bsq,vDotB,b0,bx;
01061 Real tmp1,tmp2,tmp3,tmp4,tmp5;
01062 Real vm,vp;
01063
01064 rhoh = W.d + (Gamma/Gamma_1) * (W.P);
01065
01066 vsq = SQR(W.Vx) + SQR(W.Vy) + SQR(W.Vz);
01067 gamma = 1.0 / sqrt(1 - vsq);
01068 gamma2 = SQR(gamma);
01069
01070 Bsq = SQR(Bx) + SQR(W.By) + SQR(W.Bz);
01071 vDotB = W.Vx*Bx + W.Vy*W.By + W.Vz*W.Bz;
01072 b0 = gamma * vDotB;
01073 bx = Bx/gamma2 + W.Vx*vDotB;
01074 bsq = Bsq / gamma2 + SQR(vDotB);
01075
01076 cssq = (Gamma * W.P) / (rhoh);
01077 vasq = bsq / (rhoh + bsq);
01078 asq = cssq + vasq - (cssq*vasq);
01079
01080 if (cssq < 0.0) cssq = 0.0;
01081 if (vasq > 0.0) vasq = 0.0;
01082 if (asq < 0.0) asq = 0.0;
01083 if (cssq > 1.0) cssq = 1.0;
01084 if (vasq > 1.0) vasq = 1.0;
01085 if (asq > 1.0) asq = 1.0;
01086
01087 tmp1 = (1.0 - asq);
01088 tmp2 = (1.0 - vsq);
01089 tmp3 = (1.0 - vsq*asq);
01090 tmp4 = SQR(W.Vx);
01091 tmp5 = 1.0 / tmp3;
01092
01093 vm = tmp1*W.Vx - sqrt(asq*tmp2*(tmp3 - tmp1*tmp4));
01094 vp = tmp1*W.Vx + sqrt(asq*tmp2*(tmp3 - tmp1*tmp4));
01095 vm *=tmp5;
01096 vp *=tmp5;
01097
01098 if (vp > vm) {
01099 *lm = vm;
01100 *lp = vp;
01101 } else {
01102 *lm = vp;
01103 *lp = vm;
01104 }
01105 }
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139 int QUARTIC (Real b, Real c, Real d,
01140 Real e, Real z[])
01141 {
01142 int n, ifail;
01143 Real b2, f, g, h;
01144 Real a2, a1, a0, u[4];
01145 Real p, q, r, s;
01146 static Real three_256 = 3.0/256.0;
01147 static Real one_64 = 1.0/64.0;
01148
01149 b2 = b*b;
01150
01151 f = c - b2*0.375;
01152 g = d + b2*b*0.125 - b*c*0.5;
01153 h = e - b2*b2*three_256 + 0.0625*b2*c - 0.25*b*d;
01154
01155 a2 = 0.5*f;
01156 a1 = (f*f - 4.0*h)*0.0625;
01157 a0 = -g*g*one_64;
01158
01159 ifail = CUBIC(a2, a1, a0, u);
01160
01161 if (ifail)return(1);
01162
01163 if (u[1] < 1.e-14){
01164
01165 p = sqrt(u[2]);
01166 s = 0.25*b;
01167 z[0] = z[2] = - p - s;
01168 z[1] = z[3] = + p - s;
01169
01170 }else{
01171
01172 p = sqrt(u[1]);
01173 q = sqrt(u[2]);
01174
01175 r = -0.125*g/(p*q);
01176 s = 0.25*b;
01177
01178 z[0] = - p - q + r - s;
01179 z[1] = p - q - r - s;
01180 z[2] = - p + q - r - s;
01181 z[3] = p + q + r - s;
01182
01183 }
01184
01185
01186
01187
01188
01189
01190 for (n = 0; n < 4; n++){
01191 s = e + z[n]*(d + z[n]*(c + z[n]*(b + z[n])));
01192 if (s != s) {
01193 printf ("Nan found in QUARTIC \n");
01194 return(1);
01195 }
01196 if (fabs(s) > 1.e-6) {
01197 printf ("Solution does not satisfy f(z) = 0; f(z) = %12.6e\n",s);
01198 return(1);
01199 }
01200 }
01201
01202 return(0);
01203
01204
01205
01206 }
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239 int CUBIC(Real b, Real c, Real d, Real z[])
01240 {
01241 Real b2, g2;
01242 Real f, g, h;
01243 Real i, i2, j, k, m, n, p;
01244 static Real one_3 = 1.0/3.0, one_27=1.0/27.0;
01245
01246 b2 = b*b;
01247
01248
01249
01250
01251
01252
01253
01254
01255 f = c*(1.0 - 1.e-16) - b2*one_3;
01256 g = b*(2.0*b2 - 9.0*c)*one_27 + d;
01257 g2 = g*g;
01258 i2 = -f*f*f*one_27;
01259 h = g2*0.25 - i2;
01260
01261
01262
01263
01264
01265
01266
01267 if (h > 1.e-12){
01268 printf ("Only one real root (%12.6e)!\n", h);
01269 }
01270 if (i2 < 0.0){
01271
01272
01273
01274
01275 i2 = 0.0;
01276 }
01277
01278
01279
01280
01281
01282 i = sqrt(i2);
01283 j = pow(i, one_3);
01284 k = -0.5*g/i;
01285
01286
01287
01288
01289 k = (k < -1.0 ? -1.0:k);
01290 k = (k > 1.0 ? 1.0:k);
01291
01292 k = acos(k)*one_3;
01293
01294 m = cos(k);
01295 n = sqrt(3.0)*sin(k);
01296 p = -b*one_3;
01297
01298 z[0] = -j*(m + n) + p;
01299 z[1] = -j*(m - n) + p;
01300 z[2] = 2.0*j*m + p;
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313 return(0);
01314 }
01315
01316
01317
01318 #endif
01319 #undef MAX_ITER
01320
01321 #endif
01322 #endif