Athena: rsolvers/hlle_sr.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*===========================================================================*/
00003 /*! \file hlle_sr.c
00004  *  \brief Compute 1D fluxes using an HLLE-type relativistic Riemann solver.
00005  *
00006  * PURPOSE: Compute 1D fluxes using an HLLE-type relativistic Riemann solver.
00007  *   Works for both hydro and MHD, but is very diffusive.
00008  *
00009  * HISTORY:
00010  * - First version written by Kevin Tian, 2007
00011  * - Updated by Jonathan Fulton, February 2009
00012  *
00013  * CONTAINS PUBLIC FUNCTIONS: 
00014  * - fluxes() - all Riemann solvers in Athena must have this function name and
00015  *              use the same argument list as defined in rsolvers/prototypes.h*/
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 /* solves quartic equation defined by a and returns roots in root
00042  * returns the number of Real roots 
00043  * error specifies an accuracy
00044  * currently force four Real solutions b/c it's physical */
00045 int QUARTIC (Real b, Real c, Real d, Real e, Real z[]);
00046 
00047 /* solves cubic equation defined by a and stores roots in root
00048  * returns number of Real roots */
00049 int CUBIC(Real b, Real c, Real d, Real z[]);
00050 
00051 
00052 /*----------------------------------------------------------------------------*/
00053 /*! \fn void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00054  *            const Prim1DS Wl, const Prim1DS Wr, const Real Bx, Cons1DS *pFlux)
00055  *  \brief Computes 1D fluxes
00056  *
00057  *   Input Arguments:
00058  *   - Ul,Ur = L/R-states of CONSERVED variables at cell interface
00059  *   - Bx = B in direction of slice at cell interface
00060  *   Output Arguments:
00061  *   - pFlux = pointer to fluxes of CONSERVED variables at cell interface
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 /*--- Step 1. ------------------------------------------------------------------
00079  * Compute the max and min wave speeds used in Mignone 
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 /*--- Step 1a. -----------------------------------------------------------------
00111  * If PLUTO wavespeeds are bad, fall back to the estimate used in ECHO
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   /* compute L/R fluxes */
00145   flux_LR(Ul,Wl,&Fl,Bx,&Pl);
00146   flux_LR(Ur,Wr,&Fr,Bx,&Pr);
00147 
00148   if(Sl >= 0.0){
00149     /*printf("Flux_L\n");*/
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     /*printf("Flux_R\n");*/
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     /* Compute HLL average state */
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 /*! \fn void entropy_flux (const Cons1DS Ul, const Cons1DS Ur,
00216  *            const Prim1DS Wl, const Prim1DS Wr, const Real Bx, Real *pFlux)
00217  *  \brief Compute entropy flux. */
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 /*--- Step 1. ------------------------------------------------------------------
00234  * Compute the max and min wave speeds used in Mignone 
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 /*--- Step 1a. -----------------------------------------------------------------
00266  * If PLUTO wavespeeds are bad, fall back to the estimate used in ECHO
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   /* compute L/R fluxes */
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     /* Compute HLL average state */
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   /* calculate enthalpy */
00332 
00333   theta = W.P/W.d;
00334   gmmr = Gamma / Gamma_1;
00335 
00336   h = 1.0 + gmmr*theta;
00337 
00338   /* calculate gamma */
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;        /* smallest roots, Mignone Eq 55 */
00382   Real lpl,lpr;        /* largest roots, Mignone Eq 55 */
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;  /* -- this is bx/u0 -- */
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      Define coefficients for the quartic  
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     /*printPrim1D(W);*/
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     /*SHOW(uprim,i);*/
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;        /* smallest roots, Mignone Eq 55 */
00532   Real lpl,lpr;        /* largest roots, Mignone Eq 55 */
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 /*! \fn int QUARTIC (Real b, Real c, Real d, 
00596  *             Real e, Real z[]) 
00597  *  \brief Solve a quartic equation.
00598  *
00599  * PURPOSE:
00600  *
00601  *   Solve a quartic equation in the form 
00602  *
00603  *  -   z^4 + bz^3 + cz^2 + dz + e = 0
00604  *
00605  *   For its purpose, it is assumed that ALL 
00606  *   roots are real. This makes things faster.
00607  *
00608  *
00609  * ARGUMENTS
00610  *
00611  * - b, c,
00612  * - d, e  (IN)  = coefficient of the quartic
00613  *                 z^4 + bz^3 + cz^2 + dz + e = 0
00614  *
00615  * - z[]   (OUT) = a vector containing the 
00616  *                 (real) roots of the quartic
00617  *   
00618  *
00619  * REFERENCE:
00620  *
00621  * - http://www.1728.com/quartic2.htm 
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      verify that cmax and cmin satisfy original 
00674      equation
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       printf (" z: %f ; %f ; %f ; %f\n",z[0], z[1], z[2], z[3]);
00692   */
00693 }
00694 /* *************************************************** */
00695 /*! \fn int CUBIC(Real b, Real c, Real d, Real z[])
00696  *  \brief Solve a cubic equation.
00697  *
00698  * PURPOSE:
00699  *
00700  *   Solve a cubic equation in the form 
00701  *
00702  * -    z^3 + bz^2 + cz + d = 0
00703  *
00704  *   For its purpose, it is assumed that ALL 
00705  *   roots are real. This makes things faster.
00706  *
00707  *
00708  * ARGUMENTS
00709  *
00710  * - b, c, d (IN)  = coefficient of the cubic
00711  *                    z^3 + bz^2 + cz + d = 0
00712  *
00713  * - z[]   (OUT)   = a vector containing the 
00714  *                   (real) roots of the cubic.
00715  *                   Roots should be sorted
00716  *                   in increasing order.
00717  *   
00718  *
00719  * REFERENCE:
00720  *
00721  * - http://www.1728.com/cubic2.htm 
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       the expression for f should be 
00737       f = c - b*b/3.0; however, to avoid negative
00738       round-off making h > 0.0 or g^2/4 - h < 0.0
00739       we let c --> c(1- 1.1e-16)
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      Real roots are possible only when 
00750          
00751      h <= 0 
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       printf ("i2 < 0.0 %12.6e\n",i2);
00760       return(1);
00761     */
00762     i2 = 0.0;
00763   }
00764         
00765   /* --------------------------------------
00766      i^2 must be >= g2*0.25
00767      -------------------------------------- */
00768         
00769   i = sqrt(i2);       /*  > 0   */
00770   j = pow(i, one_3);  /*  > 0   */
00771   k = -0.5*g/i;
00772         
00773   /*  this is to prevent unpleseant situation 
00774       where both g and i are close to zero       */
00775         
00776   k = (k < -1.0 ? -1.0:k);
00777   k = (k >  1.0 ?  1.0:k);
00778         
00779   k = acos(k)*one_3;       /*  pi/3 < k < 0 */
00780         
00781   m = cos(k);              /*   > 0   */
00782   n = sqrt(3.0)*sin(k);    /*   > 0   */
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      Since j, m, n > 0, it should follow that from
00791     
00792      z0 = -jm - jn + p
00793      z1 = -jm + jn + p
00794      z2 = 2jm + p
00795          
00796      z2 is the greatest of the roots, while z0 is the 
00797      smallest one.
00798      ------------------------------------------------------ */
00799         
00800   return(0);
00801 }
00802 
00803 #endif /* SPECIAL_RELATIVITY */
00804 #endif /* HLLE_FLUX */