00001 #include "../copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include <math.h>
00019 #include <stdio.h>
00020 #include <stdlib.h>
00021 #include "../defs.h"
00022 #include "../athena.h"
00023 #include "../globals.h"
00024 #include "prototypes.h"
00025 #include "../prototypes.h"
00026
00027 #ifdef HLLE_FLUX
00028 #ifdef SPECIAL_RELATIVITY
00029
00030 #if (NSCALARS > 0)
00031 #error : The SR HLLE flux does not work with passive scalars.
00032 #endif
00033
00034 void flux_LR(Cons1DS U, Prim1DS W, Cons1DS *flux, Real Bx, Real* p);
00035 void getMaxSignalSpeeds_pluto(const Prim1DS Wl, const Prim1DS Wr,
00036 const Real Bx, Real* low, Real* high);
00037 void getMaxSignalSpeeds_echo(const Prim1DS Wl, const Prim1DS Wr,
00038 const Real Bx, Real* low, Real* high);
00039 void getVChar_echo(const Prim1DS W, const Real Bx, Real* lml, Real* lmr);
00040 void getVChar_pluto(const Prim1DS W, const Real Bx, Real* lml, Real* lmr);
00041
00042
00043
00044
00045 int QUARTIC (Real b, Real c, Real d, Real e, Real z[]);
00046
00047
00048
00049 int CUBIC(Real b, Real c, Real d, Real z[]);
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00065 const Prim1DS Wl, const Prim1DS Wr, const Real Bx, Cons1DS *pFlux)
00066 {
00067 Cons1DS Fl, Fr;
00068 Cons1DS Uhll, Fhll;
00069 Prim1DS Whll;
00070 Real Pl, Pr;
00071 Real Sl, Sr;
00072 Real Sla, Sra;
00073 Real dS_1, V2l;
00074 int wave_speed_fail;
00075
00076 wave_speed_fail = 0;
00077
00078
00079
00080
00081 getMaxSignalSpeeds_pluto(Wl,Wr,Bx,&Sl,&Sr);
00082
00083 if (Sl != Sl) {
00084 wave_speed_fail = 1;
00085 printf("[hlle_sr_mhd]: NaN in Sl %10.4e %10.4e\n",Sl,Sr);
00086 Sl = -1.0;
00087 Sr = 1.0;
00088 }
00089
00090 if (Sr != Sr) {
00091 wave_speed_fail = 1;
00092 printf("[hlle_sr_mhd]: NaN in Sr %10.4e %10.4e\n",Sl,Sr);
00093 Sl = -1.0;
00094 Sr = 1.0;
00095 }
00096
00097 if (Sl < -1.0) {
00098 wave_speed_fail = 1;
00099 printf("[hlle_sr_mhd]: Superluminal Sl %10.4e %10.4e\n",Sl,Sr);
00100 Sl = -1.0;
00101 Sr = 1.0;
00102 }
00103 if (Sr > 1.0) {
00104 wave_speed_fail = 1;
00105 printf("[hlle_sr_mhd]: Superluminal Sr %10.4e %10.4e\n",Sl,Sr);
00106 Sl = -1.0;
00107 Sr = 1.0;
00108 }
00109
00110
00111
00112
00113 if (wave_speed_fail){
00114 getMaxSignalSpeeds_echo (Wl,Wr,Bx,&Sla,&Sra);
00115
00116 if (Sla != Sla) {
00117 printf("[hlle_sr_mhd]: NaN in Sl %10.4e %10.4e\n",Sl,Sr);
00118 Sla = -1.0;
00119 Sra = 1.0;
00120 }
00121
00122 if (Sra != Sra) {
00123 printf("[hlle_sr_mhd]: NaN in Sr %10.4e %10.4e\n",Sl,Sr);
00124 Sla = -1.0;
00125 Sra = 1.0;
00126 }
00127
00128 if (Sla < -1.0) {
00129 printf("[hlle_sr_mhd]: Superluminal Sl %10.4e %10.4e\n",Sl,Sr);
00130 Sla = -1.0;
00131 Sra = 1.0;
00132 }
00133 if (Sra > 1.0) {
00134 printf("[hlle_sr_mhd]: Superluminal Sr %10.4e %10.4e\n",Sl,Sr);
00135 Sla = -1.0;
00136 Sra = 1.0;
00137 }
00138
00139 Sl = Sla;
00140 Sr = Sra;
00141
00142 }
00143
00144
00145 flux_LR(Ul,Wl,&Fl,Bx,&Pl);
00146 flux_LR(Ur,Wr,&Fr,Bx,&Pr);
00147
00148 if(Sl >= 0.0){
00149
00150 pFlux->d = Fl.d;
00151 pFlux->Mx = Fl.Mx;
00152 pFlux->My = Fl.My;
00153 pFlux->Mz = Fl.Mz;
00154 pFlux->E = Fl.E;
00155 #ifdef MHD
00156 pFlux->By = Fl.By;
00157 pFlux->Bz = Fl.Bz;
00158 #endif
00159
00160 return;
00161 }
00162 else if(Sr <= 0.0){
00163
00164 pFlux->d = Fr.d;
00165 pFlux->Mx = Fr.Mx;
00166 pFlux->My = Fr.My;
00167 pFlux->Mz = Fr.Mz;
00168 pFlux->E = Fr.E;
00169 #ifdef MHD
00170 pFlux->By = Fr.By;
00171 pFlux->Bz = Fr.Bz;
00172 #endif
00173
00174 return;
00175 }
00176 else{
00177
00178
00179 dS_1 = 1.0/(Sr - Sl);
00180
00181 Uhll.d = (Sr*Ur.d - Sl*Ul.d + Fl.d - Fr.d ) * dS_1;
00182 Uhll.Mx = (Sr*Ur.Mx - Sl*Ul.Mx + Fl.Mx - Fr.Mx) * dS_1;
00183 Uhll.My = (Sr*Ur.My - Sl*Ul.My + Fl.My - Fr.My) * dS_1;
00184 Uhll.Mz = (Sr*Ur.Mz - Sl*Ul.Mz + Fl.Mz - Fr.Mz) * dS_1;
00185 Uhll.E = (Sr*Ur.E - Sl*Ul.E + Fl.E - Fr.E ) * dS_1;
00186 #ifdef MHD
00187 Uhll.By = (Sr*Ur.By - Sl*Ul.By + Fl.By - Fr.By) * dS_1;
00188 Uhll.Bz = (Sr*Ur.Bz - Sl*Ul.Bz + Fl.Bz - Fr.Bz) * dS_1;
00189 #endif
00190
00191 Fhll.d = (Sr*Fl.d - Sl*Fr.d + Sl*Sr*(Ur.d - Ul.d )) * dS_1;
00192 Fhll.Mx = (Sr*Fl.Mx - Sl*Fr.Mx + Sl*Sr*(Ur.Mx - Ul.Mx)) * dS_1;
00193 Fhll.My = (Sr*Fl.My - Sl*Fr.My + Sl*Sr*(Ur.My - Ul.My)) * dS_1;
00194 Fhll.Mz = (Sr*Fl.Mz - Sl*Fr.Mz + Sl*Sr*(Ur.Mz - Ul.Mz)) * dS_1;
00195 Fhll.E = (Sr*Fl.E - Sl*Fr.E + Sl*Sr*(Ur.E - Ul.E )) * dS_1;
00196 #ifdef MHD
00197 Fhll.By = (Sr*Fl.By - Sl*Fr.By + Sl*Sr*(Ur.By - Ul.By)) * dS_1;
00198 Fhll.Bz = (Sr*Fl.Bz - Sl*Fr.Bz + Sl*Sr*(Ur.Bz - Ul.Bz)) * dS_1;
00199 #endif
00200
00201 pFlux->d = Fhll.d;
00202 pFlux->Mx = Fhll.Mx;
00203 pFlux->My = Fhll.My;
00204 pFlux->Mz = Fhll.Mz;
00205 pFlux->E = Fhll.E;
00206 #ifdef MHD
00207 pFlux->By = Fhll.By;
00208 pFlux->Bz = Fhll.Bz;
00209 #endif
00210
00211 return;
00212 }
00213 }
00214
00215
00216
00217
00218 void entropy_flux (const Cons1DS Ul, const Cons1DS Ur,
00219 const Prim1DS Wl, const Prim1DS Wr, const Real Bx, Real *pFlux)
00220 {
00221 Real Fl, Fr;
00222 Real USl, USr;
00223 Real WSl, WSr;
00224 Real Uhll, Fhll;
00225 Real Pl, Pr;
00226 Real Sl, Sr;
00227 Real Sla, Sra;
00228 Real dS_1;
00229 int wave_speed_fail;
00230
00231 wave_speed_fail = 0;
00232
00233
00234
00235
00236 getMaxSignalSpeeds_pluto(Wl,Wr,Bx,&Sl,&Sr);
00237
00238 if (Sl != Sl) {
00239 wave_speed_fail = 1;
00240 printf("[hlle_sr_mhd]: NaN in Sl %10.4e %10.4e\n",Sl,Sr);
00241 Sl = -1.0;
00242 Sr = 1.0;
00243 }
00244
00245 if (Sr != Sr) {
00246 wave_speed_fail = 1;
00247 printf("[hlle_sr_mhd]: NaN in Sr %10.4e %10.4e\n",Sl,Sr);
00248 Sl = -1.0;
00249 Sr = 1.0;
00250 }
00251
00252 if (Sl < -1.0) {
00253 wave_speed_fail = 1;
00254 printf("[hlle_sr_mhd]: Superluminal Sl %10.4e %10.4e\n",Sl,Sr);
00255 Sl = -1.0;
00256 Sr = 1.0;
00257 }
00258 if (Sr > 1.0) {
00259 wave_speed_fail = 1;
00260 printf("[hlle_sr_mhd]: Superluminal Sr %10.4e %10.4e\n",Sl,Sr);
00261 Sl = -1.0;
00262 Sr = 1.0;
00263 }
00264
00265
00266
00267
00268 if (wave_speed_fail){
00269 getMaxSignalSpeeds_echo (Wl,Wr,Bx,&Sla,&Sra);
00270
00271 if (Sla != Sla) {
00272 printf("[hlle_sr_mhd]: NaN in Sl %10.4e %10.4e\n",Sl,Sr);
00273 Sla = -1.0;
00274 Sra = 1.0;
00275 }
00276
00277 if (Sra != Sra) {
00278 printf("[hlle_sr_mhd]: NaN in Sr %10.4e %10.4e\n",Sl,Sr);
00279 Sla = -1.0;
00280 Sra = 1.0;
00281 }
00282
00283 if (Sla < -1.0) {
00284 printf("[hlle_sr_mhd]: Superluminal Sl %10.4e %10.4e\n",Sl,Sr);
00285 Sla = -1.0;
00286 Sra = 1.0;
00287 }
00288 if (Sra > 1.0) {
00289 printf("[hlle_sr_mhd]: Superluminal Sr %10.4e %10.4e\n",Sl,Sr);
00290 Sla = -1.0;
00291 Sra = 1.0;
00292 }
00293
00294 Sl = Sla;
00295 Sr = Sra;
00296
00297 }
00298
00299
00300 WSl = Wl.P*pow(Wl.d,1.0-Gamma);
00301 WSr = Wr.P*pow(Wr.d,1.0-Gamma);
00302 USl = WSl * Ul.d/Wl.d;
00303 USr = WSr * Ur.d/Wr.d;
00304 Fl = USl * Wl.Vx;
00305 Fr = USr * Wr.Vx;
00306
00307 if(Sl >= 0.0){
00308 *pFlux = Fl;
00309 return;
00310 }
00311 else if(Sr <= 0.0){
00312 *pFlux = Fr;
00313 return;
00314 }
00315 else{
00316
00317 dS_1 = 1.0/(Sr - Sl);
00318 *pFlux = (Sr*Fl - Sl*Fr + Sl*Sr*(USr - USl)) * dS_1;
00319 return;
00320 }
00321 }
00322
00323
00324 void flux_LR(Cons1DS U, Prim1DS W, Cons1DS *flux, Real Bx, Real* p)
00325 {
00326 Real wtg2, pt, g, g2, g_2, h, gmmr, theta;
00327 #ifdef MHD
00328 Real bx, by, bz, vB, b2, Bmag2;
00329 #endif
00330
00331
00332
00333 theta = W.P/W.d;
00334 gmmr = Gamma / Gamma_1;
00335
00336 h = 1.0 + gmmr*theta;
00337
00338
00339
00340 g = U.d/W.d;
00341 g2 = SQR(g);
00342 g_2 = 1.0/g2;
00343
00344 pt = W.P;
00345 wtg2 = W.d*h*g2;
00346
00347 #ifdef MHD
00348 vB = W.Vx*Bx + W.Vy*W.By + W.Vz*W.Bz;
00349 Bmag2 = SQR(Bx) + SQR(W.By) + SQR(W.Bz);
00350
00351 bx = g*( Bx*g_2 + vB*W.Vx);
00352 by = g*(W.By*g_2 + vB*W.Vy);
00353 bz = g*(W.Bz*g_2 + vB*W.Vz);
00354
00355 b2 = Bmag2*g_2 + vB*vB;
00356
00357 pt += 0.5*b2;
00358 wtg2 += b2*g2;
00359 #endif
00360
00361 flux->d = U.d*W.Vx;
00362 flux->Mx = wtg2*W.Vx*W.Vx + pt;
00363 flux->My = wtg2*W.Vy*W.Vx;
00364 flux->Mz = wtg2*W.Vz*W.Vx;
00365 flux->E = U.Mx;
00366 #ifdef MHD
00367 flux->Mx -= bx*bx;
00368 flux->My -= by*bx;
00369 flux->Mz -= bz*bx;
00370 flux->By = W.Vx*W.By - Bx*W.Vy;
00371 flux->Bz = W.Vx*W.Bz - Bx*W.Vz;
00372 #endif
00373
00374 *p = pt;
00375 }
00376
00377 void getMaxSignalSpeeds_pluto(const Prim1DS Wl, const Prim1DS Wr,
00378 const Real Bx, Real* low, Real* high)
00379 {
00380
00381 Real lml,lmr;
00382 Real lpl,lpr;
00383 Real al,ar;
00384
00385 getVChar_pluto(Wl,Bx,&lml,&lpl);
00386 getVChar_pluto(Wr,Bx,&lmr,&lpr);
00387
00388 *low = MIN(lml, lmr);
00389 *high = MAX(lpl, lpr);
00390 }
00391
00392 void getVChar_pluto(const Prim1DS W, const Real Bx, Real* lm, Real* lp)
00393 {
00394 Real rhoh,vsq,bsq;
00395 Real cssq,vasq,asq;
00396 Real Vx2,gamma,gamma2;
00397 Real Bx2,Bsq,vDotB,vDotBsq,b0,bx;
00398 Real w_1,a0,a1,a2,a3,a4,Q;
00399 Real scrh,scrh1,scrh2,eps2;
00400 Real a2_w,one_m_eps2,lambda[5];
00401 int iflag;
00402
00403 rhoh = W.d + (Gamma/Gamma_1) * (W.P);
00404
00405 Vx2 = SQR(W.Vx);
00406 vsq = Vx2 + SQR(W.Vy) + SQR(W.Vz);
00407 if (vsq > 1.0){
00408 printf("[getVChar]: |v|= %f > 1\n",vsq);
00409 *lm = -1.0;
00410 *lp = 1.0;
00411 return;
00412 }
00413 gamma = 1.0 / sqrt(1 - vsq);
00414 gamma2 = SQR(gamma);
00415
00416 Bx2 = SQR(Bx);
00417 Bsq = Bx2 + SQR(W.By) + SQR(W.Bz);
00418 vDotB = W.Vx*Bx + W.Vy*W.By + W.Vz*W.Bz;
00419 vDotBsq = SQR(vDotB);
00420 b0 = gamma * vDotB;
00421 bx = Bx/gamma2 + W.Vx*vDotB;
00422 bsq = Bsq / gamma2 + SQR(vDotB);
00423
00424 cssq = (Gamma * W.P) / (rhoh);
00425 vasq = bsq / (rhoh + bsq);
00426
00427 if (bsq < 0.0) bsq = 0.0;
00428 if (cssq < 0.0) cssq = 0.0;
00429 if (cssq > 1.0) cssq = 1.0;
00430 if (vasq > 1.0) bsq = rhoh + bsq;
00431
00432 if (vsq < 1.0e-12) {
00433 w_1 = 1.0/(rhoh + bsq);
00434 eps2 = cssq + bsq*w_1*(1.0 - cssq);
00435 a0 = cssq*Bx*Bx*w_1;
00436 a1 = - a0 - eps2;
00437 scrh = a1*a1 - 4.0*a0;
00438 if (scrh < 0.0) scrh = 0.0;
00439
00440 scrh = sqrt(0.5*(-a1 + sqrt(scrh)));
00441 *lp = scrh;
00442 *lm = -scrh;
00443 return;
00444 }
00445
00446 w_1 = 1.0/(rhoh + bsq);
00447
00448 if (Bx < 1.0e-14) {
00449
00450 eps2 = cssq + bsq*w_1*(1.0 - cssq);
00451
00452 scrh1 = (1.0 - eps2)*gamma2;
00453 scrh2 = cssq*vDotBsq*w_1 - eps2;
00454
00455 a2 = scrh1 - scrh2;
00456 a1 = -2.0*W.Vx*scrh1;
00457 a0 = Vx2*scrh1 + scrh2;
00458
00459 *lp = 0.5*(-a1 + sqrt(a1*a1 - 4.0*a2*a0))/a2;
00460 *lm = 0.5*(-a1 - sqrt(a1*a1 - 4.0*a2*a0))/a2;
00461
00462 return;
00463 }
00464
00465 scrh1 = bx;
00466 scrh2 = scrh1*scrh1;
00467
00468 a2_w = cssq*w_1;
00469 eps2 = (cssq*rhoh + bsq)*w_1;
00470 one_m_eps2 = gamma2*rhoh*(1.0 - cssq)*w_1;
00471
00472
00473
00474
00475
00476 scrh = 2.0*(a2_w*vDotB*scrh1 - eps2*W.Vx);
00477 a4 = one_m_eps2 - a2_w*vDotBsq + eps2;
00478 a3 = - 4.0*W.Vx*one_m_eps2 + scrh;
00479 a2 = 6.0*Vx2*one_m_eps2 + a2_w*(vDotBsq - scrh2) + eps2*(Vx2 - 1.0);
00480 a1 = - 4.0*W.Vx*Vx2*one_m_eps2 - scrh;
00481 a0 = Vx2*Vx2*one_m_eps2 + a2_w*scrh2 - eps2*Vx2;
00482
00483 if (a4 < 1.e-12){
00484
00485 printf("[MAX_CH_SPEED]: Can not divide by a4 in MAX_CH_SPEED\n");
00486
00487 *lm = -1.0;
00488 *lp = 1.0;
00489
00490 return;
00491 }
00492
00493 scrh = 1.0/a4;
00494
00495 a3 *= scrh;
00496 a2 *= scrh;
00497 a1 *= scrh;
00498 a0 *= scrh;
00499 iflag = QUARTIC(a3, a2, a1, a0, lambda);
00500
00501 if (iflag){
00502 printf ("Can not find max speed:\n");
00503
00504 printf("QUARTIC: f(x) = %12.6e + x*(%12.6e + x*(%12.6e ",
00505 a0*a4, a1*a4, a2*a4);
00506 printf("+ x*(%12.6e + x*%12.6e)))\n", a3*a4, a4);
00507 printf("[MAX_CH_SPEED]: Failed to find wave speeds");
00508
00509 *lm = -1.0;
00510 *lp = 1.0;
00511
00512 return;
00513 }
00514
00515 *lp = MIN(1.0,MAX(lambda[3], lambda[2]));
00516 *lp = MIN(1.0,MAX(*lp, lambda[1]));
00517 *lp = MIN(1.0,MAX(*lp, lambda[0]));
00518
00519 *lm = MAX(-1.0,MIN(lambda[3], lambda[2]));
00520 *lm = MAX(-1.0,MIN(*lm, lambda[1]));
00521 *lm = MAX(-1.0,MIN(*lm, lambda[0]));
00522
00523 return;
00524
00525 }
00526
00527 void getMaxSignalSpeeds_echo (const Prim1DS Wl, const Prim1DS Wr,
00528 const Real Bx, Real* low, Real* high)
00529 {
00530
00531 Real lml,lmr;
00532 Real lpl,lpr;
00533 Real al,ar;
00534
00535 getVChar_echo(Wl,Bx,&lml,&lpl);
00536 getVChar_echo(Wr,Bx,&lmr,&lpr);
00537
00538 *low = MIN(lml, lmr);
00539 *high = MAX(lpl, lpr);
00540 }
00541
00542 void getVChar_echo(const Prim1DS W, const Real Bx, Real* lm, Real* lp)
00543 {
00544 Real rhoh,vsq,bsq;
00545 Real cssq,vasq,asq;
00546 Real gamma,gamma2;
00547 Real Bsq,vDotB,b0,bx;
00548 Real tmp1,tmp2,tmp3,tmp4,tmp5;
00549 Real vm,vp;
00550
00551 rhoh = W.d + (Gamma/Gamma_1) * (W.P);
00552
00553 vsq = SQR(W.Vx) + SQR(W.Vy) + SQR(W.Vz);
00554 gamma = 1.0 / sqrt(1 - vsq);
00555 gamma2 = SQR(gamma);
00556
00557 Bsq = SQR(Bx) + SQR(W.By) + SQR(W.Bz);
00558 vDotB = W.Vx*Bx + W.Vy*W.By + W.Vz*W.Bz;
00559 b0 = gamma * vDotB;
00560 bx = Bx/gamma2 + W.Vx*vDotB;
00561 bsq = Bsq / gamma2 + SQR(vDotB);
00562
00563 cssq = (Gamma * W.P) / (rhoh);
00564 vasq = bsq / (rhoh + bsq);
00565 asq = cssq + vasq - (cssq*vasq);
00566
00567 if (cssq < 0.0) cssq = 0.0;
00568 if (vasq > 0.0) vasq = 0.0;
00569 if (asq < 0.0) asq = 0.0;
00570 if (cssq > 1.0) cssq = 1.0;
00571 if (vasq > 1.0) vasq = 1.0;
00572 if (asq > 1.0) asq = 1.0;
00573
00574 tmp1 = (1.0 - asq);
00575 tmp2 = (1.0 - vsq);
00576 tmp3 = (1.0 - vsq*asq);
00577 tmp4 = SQR(W.Vx);
00578 tmp5 = 1.0 / tmp3;
00579
00580 vm = tmp1*W.Vx - sqrt(asq*tmp2*(tmp3 - tmp1*tmp4));
00581 vp = tmp1*W.Vx + sqrt(asq*tmp2*(tmp3 - tmp1*tmp4));
00582 vm *=tmp5;
00583 vp *=tmp5;
00584
00585 if (vp > vm) {
00586 *lm = vm;
00587 *lp = vp;
00588 } else {
00589 *lm = vp;
00590 *lp = vm;
00591 }
00592 }
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626 int QUARTIC (Real b, Real c, Real d,
00627 Real e, Real z[])
00628 {
00629 int n, ifail;
00630 Real b2, f, g, h;
00631 Real a2, a1, a0, u[4];
00632 Real p, q, r, s;
00633 static Real three_256 = 3.0/256.0;
00634 static Real one_64 = 1.0/64.0;
00635
00636 b2 = b*b;
00637
00638 f = c - b2*0.375;
00639 g = d + b2*b*0.125 - b*c*0.5;
00640 h = e - b2*b2*three_256 + 0.0625*b2*c - 0.25*b*d;
00641
00642 a2 = 0.5*f;
00643 a1 = (f*f - 4.0*h)*0.0625;
00644 a0 = -g*g*one_64;
00645
00646 ifail = CUBIC(a2, a1, a0, u);
00647
00648 if (ifail)return(1);
00649
00650 if (u[1] < 1.e-14){
00651
00652 p = sqrt(u[2]);
00653 s = 0.25*b;
00654 z[0] = z[2] = - p - s;
00655 z[1] = z[3] = + p - s;
00656
00657 }else{
00658
00659 p = sqrt(u[1]);
00660 q = sqrt(u[2]);
00661
00662 r = -0.125*g/(p*q);
00663 s = 0.25*b;
00664
00665 z[0] = - p - q + r - s;
00666 z[1] = p - q - r - s;
00667 z[2] = - p + q - r - s;
00668 z[3] = p + q + r - s;
00669
00670 }
00671
00672
00673
00674
00675
00676
00677 for (n = 0; n < 4; n++){
00678 s = e + z[n]*(d + z[n]*(c + z[n]*(b + z[n])));
00679 if (s != s) {
00680 printf ("Nan found in QUARTIC \n");
00681 return(1);
00682 }
00683 if (fabs(s) > 1.e-6) {
00684 printf ("Solution does not satisfy f(z) = 0; f(z) = %12.6e\n",s);
00685 return(1);
00686 }
00687 }
00688
00689 return(0);
00690
00691
00692
00693 }
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726 int CUBIC(Real b, Real c, Real d, Real z[])
00727 {
00728 Real b2, g2;
00729 Real f, g, h;
00730 Real i, i2, j, k, m, n, p;
00731 static Real one_3 = 1.0/3.0, one_27=1.0/27.0;
00732
00733 b2 = b*b;
00734
00735
00736
00737
00738
00739
00740
00741
00742 f = c*(1.0 - 1.e-16) - b2*one_3;
00743 g = b*(2.0*b2 - 9.0*c)*one_27 + d;
00744 g2 = g*g;
00745 i2 = -f*f*f*one_27;
00746 h = g2*0.25 - i2;
00747
00748
00749
00750
00751
00752
00753
00754 if (h > 1.e-12){
00755 printf ("Only one real root (%12.6e)!\n", h);
00756 }
00757 if (i2 < 0.0){
00758
00759
00760
00761
00762 i2 = 0.0;
00763 }
00764
00765
00766
00767
00768
00769 i = sqrt(i2);
00770 j = pow(i, one_3);
00771 k = -0.5*g/i;
00772
00773
00774
00775
00776 k = (k < -1.0 ? -1.0:k);
00777 k = (k > 1.0 ? 1.0:k);
00778
00779 k = acos(k)*one_3;
00780
00781 m = cos(k);
00782 n = sqrt(3.0)*sin(k);
00783 p = -b*one_3;
00784
00785 z[0] = -j*(m + n) + p;
00786 z[1] = -j*(m - n) + p;
00787 z[2] = 2.0*j*m + p;
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800 return(0);
00801 }
00802
00803 #endif
00804 #endif