Athena: convert_var.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file convert_var.c
00004  *  \brief Functions to convert conservative->primitive variables, and <->
00005  *
00006  * PURPOSE: Functions to convert conservative->primitive variables, and
00007  *   vice-versa. Routines for both Newtonian and special relativistic physics
00008  *   are included. Also contains function to compute fast magnetosonic speed for
00009  *   Newtonian flows only.
00010  *
00011  * CONTAINS PUBLIC FUNCTIONS: 
00012  * - Cons_to_Prim()     - converts Cons type to Prim type
00013  * - Cons1D_to_Prim1D() - converts 1D vector (Bx passed through arguments)
00014  * - Prim1D_to_Cons1D() - converts 1D vector (Bx passed through arguments)
00015  * - cfast()            - computes fast magnetosonic speed
00016  *
00017  * For special relativity, there are two versions of the Cons1D_to_Prim1D
00018  * functions, one for HYDRO and one for MHD. There is also a 'check_Prim'
00019  * and 'check_Prim1D' routine for SRMHD which returns a primitve state
00020  * without verifying physical status.
00021  *
00022  * REFERENCE:
00023  *   S. Noble et al., "Primitive Variable Solvers for Conservative General
00024  *   Relativistic Magnetohydrodynamics", ApJ. 641, 626 (2006)
00025  *
00026  *   A. Mignone & J. McKinney, "Equation of state in relativistic MHD: variable
00027  *   versus constant adiabatic index", MNRAS, 378, 1118 (2007)
00028  *
00029  * HISTORY: (functions for special relativity)
00030  * - First version written by Kevin Tian, summer 2007.
00031  * - Rewritten and revised by Jonathan Fulton, Senior Thesis 2009.
00032  * - Rewritten to use M&McK (2007) by Kris Beckwith, Fall 2009.               */
00033 /*============================================================================*/
00034 
00035 #include <math.h>
00036 #include "defs.h"
00037 #include "athena.h"
00038 #include "globals.h"
00039 #include "prototypes.h"
00040 
00041 #if defined(SPECIAL_RELATIVITY) && defined(MHD)
00042 /* prototypes for private functions needed with SR MHD */
00043 Real Gamma_1overGamma;
00044 
00045 /*! \fn Prim1DS vsq1D_fix (const Cons1DS *pU, const Real *pBx) 
00046  *  \brief Wrapper for the fix_vsq1D function, works 
00047  *          only for SPECIAL_RELATIVITY && MHD only */
00048 Prim1DS vsq1D_fix (const Cons1DS *pU, const Real *pBx);
00049 
00050 /*! \fn Prim1DS entropy_fix1D(const Cons1DS *U, const Real *Bx, 
00051  *                             const Real *ent) 
00052  *  \brief entropy_fix: SPECIAL RELATIVISTIC MHD VERSION */
00053 Prim1DS entropy_fix1D (const Cons1DS *pU, const Real *pBx, const Real *ent);
00054 
00055 /*! \fn static Real calc_func (Real Q, Real E, Real Bsq, Real Ssq, 
00056  *                             Real Vsq, Real pgas)
00057  *  \brief Evaluate equation A25 for E, rather than E' */
00058 static Real calc_func (Real Q, Real E, Real Bsq, Real Ssq, Real Vsq, Real pgas);
00059 
00060 /*! \fn static Real calc_dfunc(Real Q, Real Bsq, Real Msq, Real Ssq, Real d, 
00061  *                             Real Vsq, Real Gsq, Real Chi)
00062  *  \brief Evaluate A8 for E & Q, rather than E', W' */
00063 static Real calc_dfunc (Real Q, Real Bsq, Real Msq, Real Ssq, Real d, Real Vsq,
00064                         Real Gsq, Real Chi);
00065 
00066 /*! \fn static Real calc_ent_func(Real rho, Real pgas, Real d, Real ent)
00067  *  \brief Evaluate equation A25 for E, rather than E' */
00068 static Real calc_ent_func (Real rho, Real pgas, Real d, Real ent);
00069 
00070 /*! \fn static Real calc_ent_dfunc(Real Q, Real Bsq, Real Msq, Real Ssq,
00071  *                                 Real d, Real Vsq, Real Gsq, Real Chi,
00072  *                                 Real pgas, Real rho)
00073  *  \brief Evaluate A8 for E & Q, rather than E', W' */
00074 static Real calc_ent_dfunc(Real Q, Real Bsq, Real Msq, Real Ssq,
00075                            Real d, Real Vsq, Real Gsq, Real Chi,
00076                            Real pgas, Real rho);
00077 
00078 /*! \fn static Real calc_vsq(Real Bsq, Real Msq, Real Ssq, Real Q)
00079  *  \brief Obtain v^2 via eqn. A3 */
00080 static Real calc_vsq (Real Bsq, Real Msq, Real Ssq, Real Q);
00081 
00082 /*! \fn static Real calc_chi (Real d, Real Vsq, Real Gsq, Real Q)
00083  *  \brief Evaluate \Chi from eqn. A11 using Q, rathher than Q'*/
00084 static Real calc_chi (Real d, Real Vsq, Real Gsq, Real Q);
00085 
00086 void FUNV2(const Real d, const Real v2, 
00087            const Real p, const Real Bsq, const Real Msq, const Real Ssq,
00088            Real *W, Real *f);
00089 
00090 /* Parameter controlling accuracy of inversion scheme for SRMHD*/
00091 static Real tol=1.0e-10;
00092 #endif /* SPECIAL_RELATIVITY && MHD */
00093 
00094 
00095 /*----------------------------------------------------------------------------*/
00096 /*! \fn PrimS Cons_to_Prim(const ConsS *pCons)
00097  *  \brief Wrapper for the Cons1D_to_Prim1D function, works for both
00098  *   NEWTONIAN and SPECIAL_RELATIVITY
00099  *
00100  * - conserved variables = (d,M1,M2,M3,[E],[B1c,B2c,B3c],[s(n)])
00101  * - primitive variables = (d,V1,V2,V3,[P],[B1c,B2c,B3c],[r(n)])
00102  */
00103 
00104 PrimS Cons_to_Prim(const ConsS *pCons)
00105 {
00106   Cons1DS U;
00107   Prim1DS W;
00108   PrimS Prim;
00109 #if (NSCALARS > 0)
00110   int n;
00111 #endif
00112   Real Bx=0.0;
00113 
00114   U.d  = pCons->d;
00115   U.Mx = pCons->M1;
00116   U.My = pCons->M2;
00117   U.Mz = pCons->M3;
00118 #ifndef ISOTHERMAL
00119   U.E  = pCons->E;
00120 #endif /* ISOTHERMAL */
00121 #ifdef MHD
00122   Bx = pCons->B1c;
00123   U.By = pCons->B2c;
00124   U.Bz = pCons->B3c;
00125 #endif /* MHD */
00126 #if (NSCALARS > 0)
00127   for (n=0; n<NSCALARS; n++) U.s[n] = pCons->s[n];
00128 #endif
00129 
00130   W = Cons1D_to_Prim1D(&U, &Bx);
00131 
00132   Prim.d  = W.d;
00133   Prim.V1 = W.Vx;
00134   Prim.V2 = W.Vy;
00135   Prim.V3 = W.Vz;
00136 #ifndef ISOTHERMAL
00137   Prim.P = W.P;
00138 #endif /* ISOTHERMAL */
00139 #ifdef MHD
00140   Prim.B1c = Bx;
00141   Prim.B2c = W.By;
00142   Prim.B3c = W.Bz;
00143 #endif /* MHD */
00144 #if (NSCALARS > 0)
00145   for (n=0; n<NSCALARS; n++) Prim.r[n] = W.r[n];
00146 #endif
00147 
00148   return Prim;
00149 }
00150 
00151 /*----------------------------------------------------------------------------*/
00152 /*! \fn ConsS Prim_to_Cons(const PrimS *pW)
00153  *  \brief Wrapper for the Prim1D_to_Cons1D function, works for both
00154  *   NEWTONIAN and SPECIAL_RELATIVITY
00155  *
00156  * - conserved variables = (d,M1,M2,M3,[E],[B1c,B2c,B3c],[s(n)])
00157  * - primitive variables = (d,V1,V2,V3,[P],[B1c,B2c,B3c],[r(n)])
00158  */
00159 ConsS Prim_to_Cons (const PrimS *pW)
00160 {
00161   Cons1DS U;
00162   Prim1DS W;
00163   ConsS Cons;
00164 #if (NSCALARS > 0)
00165   int n;
00166 #endif
00167   Real Bx=0.0;
00168 
00169   W.d  = pW->d;
00170   W.Vx = pW->V1;
00171   W.Vy = pW->V2;
00172   W.Vz = pW->V3;
00173 #ifndef ISOTHERMAL
00174   W.P  = pW->P;
00175 #endif /* ISOTHERMAL */
00176 #ifdef MHD
00177   Bx = pW->B1c;
00178   W.By = pW->B2c;
00179   W.Bz = pW->B3c;
00180 #endif /* MHD */
00181 #if (NSCALARS > 0)
00182   for (n=0; n<NSCALARS; n++) W.r[n] = pW->r[n];
00183 #endif
00184  
00185   U = Prim1D_to_Cons1D(&W, &Bx);
00186 
00187   Cons.d  = U.d;
00188   Cons.M1 = U.Mx;
00189   Cons.M2 = U.My;
00190   Cons.M3 = U.Mz;
00191 #ifndef ISOTHERMAL
00192   Cons.E = U.E;
00193 #endif /* ISOTHERMAL */
00194 #ifdef MHD
00195   Cons.B1c = Bx;
00196   Cons.B2c = U.By;
00197   Cons.B3c = U.Bz;
00198 #endif /* MHD */
00199 #if (NSCALARS > 0)
00200   for (n=0; n<NSCALARS; n++) Cons.s[n] = U.s[n];
00201 #endif
00202 
00203   return Cons;
00204 }
00205 
00206 #ifdef SPECIAL_RELATIVITY /* special relativity only */
00207 #ifdef MHD /* MHD only */
00208 /*----------------------------------------------------------------------------*/
00209 /*! \fn PrimS fix_vsq(const ConsS *pCons)
00210  *  \brief Wrapper for the fix_vsq1D function, works 
00211  *          only for SPECIAL_RELATIVITY && MHD only
00212  *
00213  * - conserved variables = (d,M1,M2,M3,[E],[B1c,B2c,B3c],[s(n)])
00214  * - primitive variables = (d,V1,V2,V3,[P],[B1c,B2c,B3c],[r(n)])
00215  */
00216 PrimS fix_vsq (const ConsS *pCons)
00217 {
00218   Cons1DS U;
00219   Prim1DS W;
00220   PrimS Prim;
00221 #if (NSCALARS > 0)
00222   int n;
00223 #endif
00224   Real Bx=0.0;
00225 
00226   U.d  = pCons->d;
00227   U.Mx = pCons->M1;
00228   U.My = pCons->M2;
00229   U.Mz = pCons->M3;
00230 #ifndef ISOTHERMAL
00231   U.E  = pCons->E;
00232 #endif /* ISOTHERMAL */
00233 #ifdef MHD
00234   Bx = pCons->B1c;
00235   U.By = pCons->B2c;
00236   U.Bz = pCons->B3c;
00237 #endif /* MHD */
00238 #if (NSCALARS > 0)
00239   for (n=0; n<NSCALARS; n++) U.s[n] = pCons->s[n];
00240 #endif
00241 
00242   W = vsq1D_fix (&U, &Bx);
00243 
00244   Prim.d  = W.d;
00245   Prim.V1 = W.Vx;
00246   Prim.V2 = W.Vy;
00247   Prim.V3 = W.Vz;
00248 #ifndef ISOTHERMAL
00249   Prim.P = W.P;
00250 #endif /* ISOTHERMAL */
00251 #ifdef MHD
00252   Prim.B1c = Bx;
00253   Prim.B2c = W.By;
00254   Prim.B3c = W.Bz;
00255 #endif /* MHD */
00256 #if (NSCALARS > 0)
00257   for (n=0; n<NSCALARS; n++) Prim.r[n] = W.r[n];
00258 #endif
00259 
00260   return Prim;
00261 }
00262 
00263 /*----------------------------------------------------------------------------*/
00264 /*! \fn PrimS entropy_fix(const ConsS *pCons, const Real *ent)
00265  *  \brief Wrapper for the entropy_fix1D function, works 
00266  *          only for SPECIAL_RELATIVITY && MHD only
00267  *
00268  * - conserved variables = (d,M1,M2,M3,[E],[B1c,B2c,B3c],[s(n)])
00269  * - primitive variables = (d,V1,V2,V3,[P],[B1c,B2c,B3c],[r(n)])
00270  */
00271 PrimS entropy_fix (const ConsS *pCons, const Real *ent)
00272 {
00273   Cons1DS U;
00274   Prim1DS W;
00275   PrimS Prim;
00276 #if (NSCALARS > 0)
00277   int n;
00278 #endif
00279   Real Bx=0.0;
00280 
00281   U.d  = pCons->d;
00282   U.Mx = pCons->M1;
00283   U.My = pCons->M2;
00284   U.Mz = pCons->M3;
00285 #ifndef ISOTHERMAL
00286   U.E  = pCons->E;
00287 #endif /* ISOTHERMAL */
00288 #ifdef MHD
00289   Bx = pCons->B1c;
00290   U.By = pCons->B2c;
00291   U.Bz = pCons->B3c;
00292 #endif /* MHD */
00293 #if (NSCALARS > 0)
00294   for (n=0; n<NSCALARS; n++) U.s[n] = pCons->s[n];
00295 #endif
00296 
00297   W = entropy_fix1D (&U, &Bx, ent);
00298 
00299   Prim.d  = W.d;
00300   Prim.V1 = W.Vx;
00301   Prim.V2 = W.Vy;
00302   Prim.V3 = W.Vz;
00303 #ifndef ISOTHERMAL
00304   Prim.P = W.P;
00305 #endif /* ISOTHERMAL */
00306 #ifdef MHD
00307   Prim.B1c = Bx;
00308   Prim.B2c = W.By;
00309   Prim.B3c = W.Bz;
00310 #endif /* MHD */
00311 #if (NSCALARS > 0)
00312   for (n=0; n<NSCALARS; n++) Prim.r[n] = W.r[n];
00313 #endif
00314 
00315   return Prim;
00316 }
00317 #endif /* MHD */
00318 
00319 /*----------------------------------------------------------------------------*/
00320 /*! \fn PrimS check_Prim(const ConsS *pCons)
00321  *  \brief Wrapper for the check_Prim1D function, works 
00322  *          only for SPECIAL_RELATIVITY
00323  *
00324  * - conserved variables = (d,M1,M2,M3,[E],[B1c,B2c,B3c],[s(n)])
00325  * - primitive variables = (d,V1,V2,V3,[P],[B1c,B2c,B3c],[r(n)])
00326  */
00327 PrimS check_Prim(const ConsS *pCons)
00328 {
00329   Cons1DS U;
00330   Prim1DS W;
00331   PrimS Prim;
00332 #if (NSCALARS > 0)
00333   int n;
00334 #endif
00335   Real Bx=0.0;
00336 
00337   U.d  = pCons->d;
00338   U.Mx = pCons->M1;
00339   U.My = pCons->M2;
00340   U.Mz = pCons->M3;
00341 #ifndef ISOTHERMAL
00342   U.E  = pCons->E;
00343 #endif /* ISOTHERMAL */
00344 #ifdef MHD
00345   Bx = pCons->B1c;
00346   U.By = pCons->B2c;
00347   U.Bz = pCons->B3c;
00348 #endif /* MHD */
00349 #if (NSCALARS > 0)
00350   for (n=0; n<NSCALARS; n++) U.s[n] = pCons->s[n];
00351 #endif
00352 
00353 #ifdef MHD
00354   W = check_Prim1D(&U, &Bx);
00355 #else
00356   W = Cons1D_to_Prim1D(&U, &Bx);
00357 #endif
00358 
00359   Prim.d  = W.d;
00360   Prim.V1 = W.Vx;
00361   Prim.V2 = W.Vy;
00362   Prim.V3 = W.Vz;
00363 #ifndef ISOTHERMAL
00364   Prim.P = W.P;
00365 #endif /* ISOTHERMAL */
00366 #ifdef MHD
00367   Prim.B1c = Bx;
00368   Prim.B2c = W.By;
00369   Prim.B3c = W.Bz;
00370 #endif /* MHD */
00371 #if (NSCALARS > 0)
00372   for (n=0; n<NSCALARS; n++) Prim.r[n] = W.r[n];
00373 #endif
00374 
00375   return Prim;
00376 }
00377 #endif /* SPECIAL_RELATIVITY && MHD */
00378 
00379 #ifndef SPECIAL_RELATIVITY /* Following versions for Newtonian dynamics */
00380 /*----------------------------------------------------------------------------*/
00381 /*! \fn Prim1DS Cons1D_to_Prim1D(const Cons1DS *pU, const Real *pBx)
00382  *  \brief Cons1D_to_Prim1D: NEWTONIAN VERSION 
00383  *
00384  *   - conserved variables = (d,Mx,My,Mz,[E],[By,Bz],[s(n)])
00385  *   - primitive variables = (d,Vx,Vy,Vz,[P],[By,Bz],[r(n)])
00386  *   - Bx is passed in through the argument list.
00387  */
00388 
00389 Prim1DS Cons1D_to_Prim1D(const Cons1DS *pU, const Real *pBx)
00390 {
00391   Prim1DS Prim1D;
00392 #if (NSCALARS > 0)
00393   int n;
00394 #endif
00395   Real di = 1.0/pU->d;
00396 
00397   Prim1D.d  = pU->d;
00398   Prim1D.Vx = pU->Mx*di;
00399   Prim1D.Vy = pU->My*di;
00400   Prim1D.Vz = pU->Mz*di;
00401 
00402 #ifndef ISOTHERMAL
00403   Prim1D.P = pU->E - 0.5*(SQR(pU->Mx)+SQR(pU->My)+SQR(pU->Mz))*di;
00404 #ifdef MHD
00405   Prim1D.P -= 0.5*(SQR(*pBx) + SQR(pU->By) + SQR(pU->Bz));
00406 #endif /* MHD */
00407   Prim1D.P *= Gamma_1;
00408   Prim1D.P = MAX(Prim1D.P,TINY_NUMBER);
00409 #endif /* ISOTHERMAL */
00410 
00411 #ifdef MHD
00412   Prim1D.By = pU->By;
00413   Prim1D.Bz = pU->Bz;
00414 #endif /* MHD */
00415 
00416 #if (NSCALARS > 0)
00417   for (n=0; n<NSCALARS; n++) Prim1D.r[n] = pU->s[n]*di;
00418 #endif
00419 
00420   return Prim1D;
00421 }
00422 
00423 /*----------------------------------------------------------------------------*/
00424 /*! \fn Cons1DS Prim1D_to_Cons1D(const Prim1DS *pW, const Real *pBx)
00425  *  \brief Prim1D_to_Cons1D: NEWTONIAN VERSION
00426  *
00427  *   - primitive variables = (d,Vx,Vy,Vz,[P],[By,Bz],[r(n)])
00428  *   - conserved variables = (d,Mx,My,Mz,[E],[By,Bz],[s(n)])
00429  *   - Bx is passed in through the argument list.
00430  */
00431 
00432 Cons1DS Prim1D_to_Cons1D(const Prim1DS *pW, const Real *pBx)
00433 {
00434   Cons1DS Cons1D;
00435 #if (NSCALARS > 0)
00436   int n;
00437 #endif
00438 
00439   Cons1D.d  = pW->d;
00440   Cons1D.Mx = pW->d*pW->Vx;
00441   Cons1D.My = pW->d*pW->Vy;
00442   Cons1D.Mz = pW->d*pW->Vz;
00443 
00444 #ifndef ISOTHERMAL
00445   Cons1D.E = pW->P/Gamma_1 + 0.5*pW->d*(SQR(pW->Vx) +SQR(pW->Vy) +SQR(pW->Vz));
00446 #ifdef MHD
00447   Cons1D.E += 0.5*(SQR(*pBx) + SQR(pW->By) + SQR(pW->Bz));
00448 #endif /* MHD */
00449 #endif /* ISOTHERMAL */
00450 
00451 #ifdef MHD
00452   Cons1D.By = pW->By;
00453   Cons1D.Bz = pW->Bz;
00454 #endif /* MHD */
00455 
00456 #if (NSCALARS > 0)
00457   for (n=0; n<NSCALARS; n++) Cons1D.s[n] = pW->r[n]*pW->d;
00458 #endif
00459 
00460   return Cons1D;
00461 }
00462 
00463 /*----------------------------------------------------------------------------*/
00464 /*! \fn Real cfast(const Cons1DS *U, const Real *Bx)
00465  *
00466  *  \brief Returns fast magnetosonic speed given input 1D vector of conserved
00467  *   variables and Bx -- NEWTONIAN PHYSICS ONLY.
00468  */
00469 
00470 Real cfast(const Cons1DS *U, const Real *Bx)
00471 {
00472   Real asq;
00473 #ifndef ISOTHERMAL
00474   Real p,pb=0.0;
00475 #endif
00476 #ifdef MHD
00477   Real ctsq,casq,tmp,cfsq;
00478 #endif
00479 
00480 #ifdef ISOTHERMAL
00481   asq = Iso_csound2;
00482 #else
00483 #ifdef MHD
00484   pb = 0.5*(SQR(*Bx) + SQR(U->By) + SQR(U->Bz));
00485 #endif /* MHD */
00486   p = Gamma_1*(U->E - pb - 0.5*(SQR(U->Mx)+SQR(U->My)+SQR(U->Mz))/U->d);
00487   asq = Gamma*p/U->d;
00488 #endif /* ISOTHERMAL */
00489 
00490 #ifndef MHD
00491   return sqrt(asq);
00492 #else
00493   ctsq = (SQR(U->By) + SQR(U->Bz))/U->d;
00494   casq = SQR(*Bx)/U->d;
00495   tmp = casq + ctsq - asq;
00496   cfsq = 0.5*((asq+ctsq+casq) + sqrt(tmp*tmp + 4.0*asq*ctsq));
00497   return sqrt(cfsq);
00498 #endif
00499 }
00500 #endif /* not SPECIAL_RELATIVITY */
00501 
00502 #if defined(SPECIAL_RELATIVITY) && defined(HYDRO) /* special relativity only */
00503 /*----------------------------------------------------------------------------*/
00504 /*! \fn Prim1DS Cons1D_to_Prim1D(const Cons1DS *U, const Real *Bx)
00505  *  \brief Cons1D_to_Prim1D: SPECIAL RELATIVISTIC HYDRODYNAMICS VERSION
00506  *
00507  *   - conserved variables = (d,Mx,My,Mz,[E],[By,Bz])
00508  *   - primitive variables = (d,Vx,Vy,Vz,[P],[By,Bz])
00509  *   - Bx is passed in through the argument list.
00510  */
00511 
00512 Prim1DS Cons1D_to_Prim1D(const Cons1DS *U, const Real *Bx)
00513 {
00514   Prim1DS Prim1D;
00515   Real Msq, M, ME, Dsq, Gamma_1sq, denom;
00516   Real a3, a2, a1, a0;
00517   Real i1, i2, i3;
00518   Real iR, iS, iT;
00519   Real ix1=0.0, iB, iC, v, vOverM;
00520 
00521   /* Step One: Reduce equations to a quartic polynomial in v */
00522 
00523   Msq = SQR(U->Mx) + SQR(U->My) + SQR(U->Mz);
00524   M = sqrt(Msq);
00525 
00526   if (fabs(M) < TINY_NUMBER) {
00527 
00528     v = 0.0;   /* if M is zero, then v must be identically zero */
00529 
00530   } else {
00531 
00532     ME = M * U->E;
00533     Dsq = SQR(U->d);
00534 
00535     Gamma_1sq = SQR(Gamma_1);
00536 
00537     denom = 1.0 / (Gamma_1sq * (Msq + Dsq));
00538 
00539     a3 = (-2.0 * Gamma * Gamma_1 * ME) * denom;
00540     a2 = (SQR(Gamma) * SQR(U->E) + 2.0*Gamma_1*Msq - Gamma_1sq*Dsq)*denom;
00541     a1 = (-2.0 * Gamma * ME) * denom;
00542     a0 = Msq * denom;
00543 
00544     /* Step Two: Solve the polynomial analytically */
00545 
00546     i1 = -a2;
00547     i2 = a3 * a1 - 4.0 * a0;
00548     i3 = 4.0 * a2 * a0 - SQR(a1) - SQR(a3) * a0;
00549 
00550     iR = (9.0 * i1 * i2 - 27.0 * i3 - 2.0 * SQR(i1) * i1) / 54.0;
00551     iS = (3.0 * i2 - SQR(a2)) / 9.0;
00552     iT = SQR(iR) + SQR(iS) * iS;
00553 
00554     /* iT may be negative, but ix1 then adds a complex number and its conjugate,
00555      * giving us a real value, so ix1 is real no matter what */
00556 
00557     if (iT < 0) {
00558       ix1 = 2.0*pow(sqrt(iR*iR + iT),(ONE_3RD))*cos(atan2(sqrt(-iT),iR)/3.0)
00559         - i1/3.0;
00560     } else {
00561       ix1 = pow((iR + sqrt(iT)),(ONE_3RD)) + pow((iR - sqrt(iT)),(ONE_3RD)) 
00562         - i1/3.0;
00563     }
00564 
00565     iB = 0.5*(a3 + sqrt(SQR(a3) - 4.0*a2 + 4.0*ix1));
00566     iC = 0.5*(ix1 - sqrt(SQR(ix1) - 4.0*a0));
00567     v = (-iB + sqrt(SQR(iB) - 4.0*iC))/2.0;
00568 
00569   }
00570 
00571   /* Step Three: Resolve primitives with the value of v */
00572 
00573   /* ensure that v is physical */
00574   v = MAX(v, 0.0);
00575   v = MIN(v, 1.0 - 1.0e-15);
00576 
00577   vOverM = v / M;
00578 
00579   if (fabs(M) < TINY_NUMBER) {
00580     vOverM = 0.0;
00581   }
00582 
00583   Prim1D.d = sqrt(1.0 - SQR(v)) * U->d;
00584 
00585   Prim1D.Vx = U->Mx * vOverM;
00586   Prim1D.Vy = U->My * vOverM;
00587   Prim1D.Vz = U->Mz * vOverM;
00588 
00589   Prim1D.P = Gamma_1*
00590     ((U->E - U->Mx*Prim1D.Vx - U->My*Prim1D.Vy - U->Mz*Prim1D.Vz) - Prim1D.d);
00591 
00592   return Prim1D;
00593 }
00594 #endif /* SPECIAL_RELATIVITY && HYDRO */
00595 
00596 #if defined(SPECIAL_RELATIVITY) && defined(MHD) /* special relativity only */
00597 /*----------------------------------------------------------------------------*/
00598 /*! \fn Prim1DS Cons1D_to_Prim1D(const Cons1DS *U, const Real *Bx)
00599  *  \brief Cons1D_to_Prim1D: SPECIAL RELATIVISTIC MHD VERSION 
00600  *
00601  *   - conserved variables = (d,Mx,My,Mz,[E],[By,Bz])
00602  *   - primitive variables = (d,Vx,Vy,Vz,[P],[By,Bz])
00603  *   - Bx is passed in through the argument list.
00604  *
00605  * IMPORTANT: This algorithm uses an iterative (Newton-Raphson) root-finding
00606  * step, which requires an initial guess for W. This is provided by solving
00607  * a cubic equation for W based on passed values of U->E & U->d and assuming
00608  * that v^2 = 1, as in Appendix A3 of Mignone & McKinney. Note that the
00609  * conserved quantity is the total energy, 
00610  * - E = D*h*\gamma - p + 0.5*B^2 + 0.5*(v^2*B^2 - v \dot B)
00611  */
00612 
00613 Prim1DS Cons1D_to_Prim1D(const Cons1DS *U, const Real *Bx)
00614 {
00615   Prim1DS Prim1D;
00616   Real Bsq = 0.0, Msq = 0.0, S = 0.0, Ssq = 0.0;
00617   Real Qp = 0.0, Q = 0.0, Ep = 0.0, E = 0.0, d = 0.0;
00618   Real scrh1, scrh2, tmp1, tmp2;
00619   Real Usq, Vsq, Gsq, rho, Chi, pgas, ent;
00620   Real fQ, dfQ, dQstep;
00621 
00622   int nr_success;
00623   int q_incs;
00624 
00625   Gamma_1overGamma = Gamma_1/Gamma;
00626 
00627   /* Calculate Bsq = B^2, Msq = M^2, S = M \dot B, Ssq = S^2 */
00628   Bsq = SQR((*Bx)) + SQR(U->By) + SQR(U->Bz);
00629   Msq = SQR(U->Mx) + SQR(U->My) + SQR(U->Mz);
00630   S = U->Mx * (*Bx) + U->My * U->By + U->Mz * U->Bz;
00631   Ssq = SQR(S);
00632 
00633   /* Assign input energy & density to local variables */
00634   E = U->E;
00635   d = U->d;
00636 
00637   /* Starting guess for W, based on taking the +ve */
00638   /* root of Eqn. A27, guarantees that p is +ve    */
00639   scrh1 = -4.0*(E - Bsq);
00640   scrh2 = Msq - 2.0*E*Bsq + Bsq*Bsq;
00641   Q = ( - scrh1 + sqrt(fabs(scrh1*scrh1 - 12.0*scrh2)))/6.0;
00642 
00643   nr_success = 0;
00644   if (Q < 0.0) {
00645     Q = d;
00646   } else if (Q != Q) {
00647     nr_success = -4;
00648   }
00649 
00650   /* 1d NR algorithm to find W from A1, converges */
00651   /* when W'(k+1) / W(k) < tol, or max iterations reached */
00652   dQstep = 1.0;
00653   q_incs = 0;
00654   while (nr_success == 0 && q_incs < 100000) {
00655 
00656     if (fabs(dQstep) <= tol) nr_success = 1;
00657 
00658     /* Obtain scalar quantities based on current solution estimate */
00659     Vsq = calc_vsq (Bsq,Msq,Ssq,Q);
00660     if (Vsq != Vsq){
00661       nr_success = -3;
00662     }
00663     Gsq = 1.0/(1.0-Vsq);
00664     Chi = calc_chi (d,Vsq,Gsq,Q);
00665     rho = d / sqrt(fabs(Gsq));
00666     pgas = Gamma_1*Chi/Gamma;
00667 
00668     /* Evaluate Eqn. A24, A25 for the total energy density */ 
00669     fQ = calc_func (Q, E, Bsq, Ssq, Vsq, pgas);
00670     dfQ = calc_dfunc (Q, Bsq, Msq, Ssq, d, Vsq, Gsq, Chi);
00671 
00672     /* Check that we didn't get a NaN in the process */
00673     /*printf ("Function evaluations %12.6e %12.6e\n",fQ,dfQ);*/
00674     if (fQ != fQ) {
00675       nr_success = -1;
00676     }
00677     if (dfQ != dfQ) {
00678       nr_success = -2;
00679     }
00680 
00681     /* If we start out very close to the solution try and make sure we don't
00682      * overshoot. */
00683     if (fabs(fQ) < 0.1 && q_incs == 0) {
00684       Q *= 10;
00685       Vsq = calc_vsq (Bsq,Msq,Ssq,Q);
00686       Gsq = 1.0/(1.0-Vsq);
00687       Chi = calc_chi (d,Vsq,Gsq,Q);
00688       rho = d / sqrt(fabs(Gsq));
00689       pgas = Gamma_1*Chi/Gamma;
00690 
00691       fQ = calc_func (Q, E, Bsq, Ssq, Vsq, pgas);
00692       dfQ = calc_dfunc (Q, Bsq, Msq, Ssq, d, Vsq, Gsq, Chi);
00693     }
00694     
00695     /* Calculate step for NR update, check that it's not a Nan, update/check Q */
00696     dQstep = fQ / dfQ;
00697     if (dQstep != dQstep){
00698       nr_success = -2;
00699     }
00700     Q -= dQstep;
00701     if (Q != Q){
00702       nr_success = -1;
00703     }
00704     q_incs ++;
00705 
00706     /* Rinse and repeat */
00707   }
00708 
00709   /* If we convered (indicated by nr_success = 1) then check solution */
00710   if (nr_success == 1){
00711     /*Q = MAX(Q, d*(1.0e0+1.0e-6));*/
00712 
00713     Vsq = calc_vsq (Bsq,Msq,Ssq,Q);
00714     Gsq = 1.0/(1.0-Vsq);
00715     Chi = calc_chi (d,Vsq,Gsq,Q);
00716     rho = d / sqrt(fabs(Gsq));
00717     pgas = Gamma_1*Chi/Gamma;
00718 
00719     if (pgas < 0.0) nr_success = 2;
00720     if (Vsq > 1.0) nr_success =  3;
00721     if (Vsq < 0.0) nr_success =  4;
00722   }
00723 
00724   if (nr_success == 2) {
00725     /* Case of p < 0, in which case we return floor values */
00726     /* Should be amended to either solve the cold MHD eqn's, fall back on the
00727      * entropy or tell the integrator to switch to a more diffusive update */
00728 
00729     tmp1 = 1.0 / Q;
00730     tmp2 = 1.0 / (Q + Bsq);
00731     Prim1D.d = MAX(rho,1.0e-4);
00732     Prim1D.P = MAX(pgas,1.0e-5);
00733     Prim1D.Vx = (U->Mx + S*(*Bx)*tmp1)*tmp2;
00734     Prim1D.Vy = (U->My + S*U->By*tmp1)*tmp2;
00735     Prim1D.Vz = (U->Mz + S*U->Bz*tmp1)*tmp2;
00736 
00737     Prim1D.By = U->By;
00738     Prim1D.Bz = U->Bz;
00739   } else if (nr_success == 3) {
00740     /* Case of V^2 > 0, in which case we try again with some rescalings */
00741 
00742     tmp1 = 1.0 / Q;
00743     tmp2 = 1.0 / (Q + Bsq);
00744     Prim1D.Vx = (U->Mx + S*(*Bx)*tmp1)*tmp2;
00745     Prim1D.Vy = (U->My + S*U->By*tmp1)*tmp2;
00746     Prim1D.Vz = (U->Mz + S*U->Bz*tmp1)*tmp2;
00747 
00748     scrh1 = SQR(Prim1D.Vx) + SQR(Prim1D.Vy) + SQR(Prim1D.Vz);
00749 
00750     Prim1D.Vx *= 0.9999/scrh1;
00751     Prim1D.Vy *= 0.9999/scrh1;
00752     Prim1D.Vz *= 0.9999/scrh1;
00753     Vsq = SQR(Prim1D.Vx) + SQR(Prim1D.Vy) + SQR(Prim1D.Vz);
00754 
00755     Gsq = 1.0/(1.0-Vsq);
00756     Chi = calc_chi (d,Vsq,Gsq,Q);
00757     rho = d / sqrt(fabs(Gsq));
00758     pgas = Gamma_1*Chi/Gamma;
00759 
00760     Prim1D.d = MAX(rho,1.0e-4);
00761     Prim1D.P = MAX(pgas,1.0e-5);
00762     Prim1D.By = U->By;
00763     Prim1D.Bz = U->Bz;
00764 
00765   } else if (nr_success == 4) {
00766     /* Case of v^2 < 0, returns unphysical state */
00767     Prim1D.d = -1.0;
00768     Prim1D.P =  1.0;
00769     Prim1D.Vx = 1.0;
00770     Prim1D.Vy =1.0;
00771     Prim1D.Vz = 1.0;
00772 
00773     Prim1D.By = U->By;
00774     Prim1D.Bz = U->Bz;
00775   } else if (nr_success == 1) {
00776     /* It worked!!! Should have a valid solution, so now set up primitives */
00777     tmp1 = 1.0 / Q;
00778     tmp2 = 1.0 / (Q + Bsq);
00779     Prim1D.d = MAX(rho,1.0e-4);
00780     Prim1D.P = MAX(pgas,1.0e-5);
00781     Prim1D.Vx = (U->Mx + S*(*Bx)*tmp1)*tmp2;
00782     Prim1D.Vy = (U->My + S*U->By*tmp1)*tmp2;
00783     Prim1D.Vz = (U->Mz + S*U->Bz*tmp1)*tmp2;
00784 
00785     Prim1D.By = U->By;
00786     Prim1D.Bz = U->Bz;
00787   } else {
00788     /* NR iteration failed to converge, return unphysical values */
00789     Prim1D.d = -2.0;
00790     Prim1D.P =  2.0;
00791     Prim1D.Vx = 2.0;
00792     Prim1D.Vy = 2.0;
00793     Prim1D.Vz = 2.0;
00794 
00795     Prim1D.By = U->By;
00796     Prim1D.Bz = U->Bz;
00797   }
00798 
00799   return Prim1D;
00800 }
00801 
00802 /*----------------------------------------------------------------------------*/
00803 /*! \fn Prim1DS check_Prim1D(const Cons1DS *U, const Real *Bx) 
00804  *  \brief check_Prim1D: SPECIAL RELATIVISTIC MHD VERSION 
00805  *
00806  *   - conserved variables = (d,Mx,My,Mz,[E],[By,Bz])
00807  *   - primitive variables = (d,Vx,Vy,Vz,[P],[By,Bz])
00808  *   - Bx is passed in through the argument list.
00809  *
00810  * IMPORTANT: This algorithm uses an iterative (Newton-Raphson) root-finding
00811  * step, which requires an initial guess for W. This is provided by solving
00812  * a cubic equation for W based on passed values of U->E & U->d and assuming
00813  * that v^2 = 1, as in Appendix A3 of Mignone & McKinney. Note that the
00814  * conserved quantity is the total energy, 
00815  * - E = D*h*\gamma - p + 0.5*B^2 + 0.5*(v^2*B^2 - v \dot B)
00816  */
00817 
00818 Prim1DS check_Prim1D (const Cons1DS *U, const Real *Bx)
00819 {
00820   Prim1DS Prim1D;
00821   Real Bsq = 0.0, Msq = 0.0, S = 0.0, Ssq = 0.0;
00822   Real Qp = 0.0, Q = 0.0, Ep = 0.0, E = 0.0, d = 0.0;
00823   Real scrh1, scrh2, tmp1, tmp2;
00824   Real Usq, Vsq, Gsq, rho, Chi, pgas, ent;
00825   Real fQ, dfQ, dQstep;
00826 
00827   int nr_success;
00828   int q_incs;
00829 
00830   Gamma_1overGamma = Gamma_1/Gamma;
00831 
00832   /* Calculate Bsq = B^2, Msq = M^2, S = M \dot B, Ssq = S^2 */
00833   Bsq = SQR((*Bx)) + SQR(U->By) + SQR(U->Bz);
00834   Msq = SQR(U->Mx) + SQR(U->My) + SQR(U->Mz);
00835   S = U->Mx * (*Bx) + U->My * U->By + U->Mz * U->Bz;
00836   Ssq = SQR(S);
00837 
00838   /* Assign input energy & density to local variables */
00839   E = U->E;
00840   d = U->d;
00841 
00842   /* Starting guess for W, based on taking the +ve */
00843   /* root of Eqn. A27, guarantees that p is +ve    */
00844   scrh1 = -4.0*(E - Bsq);
00845   scrh2 = Msq - 2.0*E*Bsq + Bsq*Bsq;
00846   Q = ( - scrh1 + sqrt(fabs(scrh1*scrh1 - 12.0*scrh2)))/6.0;
00847 
00848   nr_success = 0;       
00849   if (Q < 0.0) {
00850     Q = d;
00851   } else if (Q != Q) {
00852     nr_success = 9;
00853   }
00854 
00855   /* 1d NR algorithm to find W from A1, converges */
00856   /* when W'(k+1) / W(k) < tol, or max iterations reached */
00857   dQstep = 1.0;
00858   q_incs = 0;
00859   while (nr_success == 0 && q_incs < 1000) {
00860 
00861     if (fabs(dQstep) <= tol) nr_success = 1;
00862 
00863     /* Obtain scalar quantities based on current solution estimate */
00864     Vsq = calc_vsq (Bsq,Msq,Ssq,Q);
00865     Gsq = 1.0/(1.0-Vsq);
00866     Chi = calc_chi (d,Vsq,Gsq,Q);
00867     rho = d / sqrt(fabs(Gsq));
00868     pgas = Gamma_1*Chi/Gamma;
00869 
00870     /* Evaluate Eqn. A24, A25 for the total energy density */ 
00871     fQ = calc_func (Q, E, Bsq, Ssq, Vsq, pgas);
00872     dfQ = calc_dfunc (Q, Bsq, Msq, Ssq, d, Vsq, Gsq, Chi);
00873 
00874     /* Check that we didn't get a NaN in the process */
00875     if (fQ != fQ) {
00876       nr_success = 8;
00877     }
00878     if (dfQ != dfQ) {
00879       nr_success = 7;
00880     }
00881 
00882     /* If we start out very close to the solution try and make sure we don't
00883      * overshoot.--- Not actually clear that this does anything, so it can probably
00884      * be removed */
00885     if (fabs(fQ) < 0.1 && q_incs == 0) {
00886       Q *= 10;
00887       Vsq = calc_vsq (Bsq,Msq,Ssq,Q);
00888       Gsq = 1.0/(1.0-Vsq);
00889       Chi = calc_chi (d,Vsq,Gsq,Q);
00890       rho = d / sqrt(fabs(Gsq));
00891       pgas = Gamma_1*Chi/Gamma;
00892 
00893       fQ = calc_func (Q, E, Bsq, Ssq, Vsq, pgas);
00894       dfQ = calc_dfunc (Q, Bsq, Msq, Ssq, d, Vsq, Gsq, Chi);
00895     }
00896     
00897     /* Calculate step for NR update, check that it's not a Nan, update/check Q */
00898     dQstep = fQ / dfQ;
00899     if (dQstep != dQstep){
00900       nr_success = 6;   
00901     }
00902     Q -= dQstep;
00903     if (Q != Q){
00904       nr_success = 5;
00905     }
00906     q_incs ++;
00907 
00908     /* Rinse and repeat */
00909 
00910   }
00911 
00912   /* If we convered (indicated by nr_success = 1) then check solution */
00913   if (nr_success == 1){
00914 
00915     Vsq = calc_vsq (Bsq,Msq,Ssq,Q);
00916     Gsq = 1.0/(1.0-Vsq);
00917     Chi = calc_chi (d,Vsq,Gsq,Q);
00918     rho = d / sqrt(fabs(Gsq));
00919     pgas = Gamma_1*Chi/Gamma;
00920  
00921     tmp1 = 1.0 / Q;
00922     tmp2 = 1.0 / (Q + Bsq);
00923     Prim1D.d = rho;
00924     Prim1D.P = pgas;
00925     Prim1D.Vx = (U->Mx + S*(*Bx)*tmp1)*tmp2;
00926     Prim1D.Vy = (U->My + S*U->By*tmp1)*tmp2;
00927     Prim1D.Vz = (U->Mz + S*U->Bz*tmp1)*tmp2;
00928         
00929     Prim1D.By = U->By;
00930     Prim1D.Bz = U->Bz;
00931           
00932   } else {
00933     Prim1D.d = -1.0;
00934     Prim1D.P = -1.0;
00935     Prim1D.Vx = 1.0;
00936     Prim1D.Vy = 1.0;
00937     Prim1D.Vz = 1.0;
00938           
00939     Prim1D.By = U->By;
00940     Prim1D.Bz = U->Bz;
00941   }
00942           
00943   return Prim1D;
00944 }
00945 #endif /* SPECIAL_RELATIVITY && MHD */
00946 
00947 #ifdef SPECIAL_RELATIVITY /* special relativity only */
00948 /*----------------------------------------------------------------------------*/
00949 /*! \fn Cons1DS Prim1D_to_Cons1D(const Prim1DS *W, const Real *Bx)
00950  *  \brief Prim1D_to_Cons1D: SPECIAL RELATIVITY VERSION
00951  *
00952  *   - primitive variables = (d,Vx,Vy,Vz,[P],[By,Bz])
00953  *   - conserved variables = (d,Mx,My,Mz,[E],[By,Bz])
00954  *   - Bx is passed in through the argument list.
00955  */
00956 
00957 Cons1DS Prim1D_to_Cons1D(const Prim1DS *W, const Real *Bx)
00958 {
00959   Cons1DS Cons1D;
00960   Real vsq, U0, Bsq = 0.0, vDotB = 0.0, wU0sq;
00961 
00962   vsq = SQR(W->Vx) + SQR(W->Vy) + SQR(W->Vz);
00963   U0 = 1.0 / (1.0 - vsq);
00964 
00965 #ifdef MHD
00966   Bsq = (*Bx)*(*Bx) + SQR(W->By) + SQR(W->Bz);
00967   vDotB = (*Bx)*W->Vx + W->By*W->Vy + W->Bz*W->Vz;
00968 #endif
00969 
00970   wU0sq = (W->d + Gamma/Gamma_1 * W->P)*U0;
00971 
00972   Cons1D.d  = sqrt(U0) * W->d;
00973   Cons1D.Mx = wU0sq * W->Vx;
00974   Cons1D.My = wU0sq * W->Vy;
00975   Cons1D.Mz = wU0sq * W->Vz;
00976 
00977   Cons1D.E  = wU0sq - W->P;
00978 
00979 #ifdef MHD
00980   Cons1D.Mx += Bsq*W->Vx - vDotB*(*Bx);
00981   Cons1D.My += Bsq*W->Vy - vDotB*W->By;
00982   Cons1D.Mz += Bsq*W->Vz - vDotB*W->Bz;
00983 
00984   Cons1D.E  += (1.0 + vsq) * Bsq/2.0 - SQR(vDotB) / 2.0;
00985 
00986   Cons1D.By = W->By;
00987   Cons1D.Bz = W->Bz;
00988 #endif /* MHD */
00989 
00990   return Cons1D;
00991 }
00992 #endif /* SPECIAL_RELATIVITY */
00993 
00994 #if defined(SPECIAL_RELATIVITY) && defined(MHD)
00995 /*=========================== PRIVATE FUNCTIONS ==============================*/
00996 
00997 /*----------------------------------------------------------------------------*/
00998 /*! \fn Prim1DS entropy_fix1D(const Cons1DS *U, const Real *Bx, 
00999  *                             const Real *ent) 
01000  *  \brief entropy_fix: SPECIAL RELATIVISTIC MHD VERSION 
01001  *
01002  *  - conserved variables = (d,Mx,My,Mz,[E],[By,Bz])
01003  *  - primitive variables = (d,Vx,Vy,Vz,[P],[By,Bz])
01004  *  - Bx is passed in through the argument list.
01005  */
01006 
01007 Prim1DS entropy_fix1D (const Cons1DS *U, const Real *Bx, const Real *ent)
01008 {
01009   Prim1DS Prim1D;
01010   Real Bsq = 0.0, Msq = 0.0, S = 0.0, Ssq = 0.0, Ent = 0.0;
01011   Real Qp = 0.0, Q = 0.0, Ep = 0.0, E = 0.0, d = 0.0;
01012   Real scrh1, scrh2, tmp1, tmp2;
01013   Real Usq, Vsq, Gsq, rho, Chi, pgas;
01014   Real fQ, dfQ, dQstep;
01015 
01016   int nr_success;
01017   int q_incs;
01018 
01019   Gamma_1overGamma = Gamma_1/Gamma;
01020 
01021   /* Calculate Bsq = B^2, Msq = M^2, S = M \dot B, Ssq = S^2 */
01022   Bsq = SQR((*Bx)) + SQR(U->By) + SQR(U->Bz);
01023   Msq = SQR(U->Mx) + SQR(U->My) + SQR(U->Mz);
01024   S = U->Mx * (*Bx) + U->My * U->By + U->Mz * U->Bz;
01025   Ssq = SQR(S);
01026 
01027   /* Assign input energy & density to local variables */
01028   E = U->E;
01029   d = U->d;
01030   Ent = *ent;
01031 
01032   /* Starting guess for W, based on taking the +ve */
01033   /* root of Eqn. A27, guarantees that p is +ve    */
01034   scrh1 = -4.0*(E - Bsq);
01035   scrh2 = Msq - 2.0*E*Bsq + Bsq*Bsq;
01036   Q = ( - scrh1 + sqrt(fabs(scrh1*scrh1 - 12.0*scrh2)))/6.0;
01037 
01038   nr_success = 0;       
01039   if (Q < 0.0) {
01040     Q = d;
01041   } else if (Q != Q) {
01042     nr_success = 9;
01043   }
01044 
01045   /* 1d NR algorithm to find W from A1, converges */
01046   /* when W'(k+1) / W(k) < tol, or max iterations reached */
01047   dQstep = 1.0;
01048   q_incs = 0;
01049   while (nr_success == 0 && q_incs < 1000) {
01050 
01051     if (fabs(dQstep) <= tol) nr_success = 1;
01052 
01053     /* Obtain scalar quantities based on current solution estimate */
01054     Vsq = calc_vsq (Bsq,Msq,Ssq,Q);
01055     Gsq = 1.0/(1.0-Vsq);
01056     Chi = calc_chi (d,Vsq,Gsq,Q);
01057     rho = d / sqrt(fabs(Gsq));
01058     pgas = Gamma_1*Chi/Gamma;
01059 
01060     /* Evaluate Eqn. A24, A25 for the total energy density */ 
01061     fQ = calc_ent_func (rho, pgas, d, Ent);
01062     dfQ = calc_ent_dfunc (Q, Bsq, Msq, Ssq, d, Vsq, Gsq, Chi, pgas, rho);
01063 
01064     /* Check that we didn't get a NaN in the process */
01065     /*printf ("Function evaluations %12.6e %12.6e\n",fQ,dfQ);*/
01066     if (fQ != fQ) {
01067       nr_success = 8;
01068     }
01069     if (dfQ != dfQ) {
01070       nr_success = 7;
01071     }
01072     
01073     /* Calculate step for NR update, check that it's not a Nan, update/check Q */
01074     dQstep = fQ / dfQ;
01075     if (dQstep != dQstep){
01076       nr_success = 6;       }
01077     Q -= dQstep;
01078     if (Q != Q){
01079       nr_success = 5;
01080     }
01081     q_incs ++;
01082 
01083     /* Rinse and repeat */
01084 
01085   }
01086 
01087   /* If we convered (indicated by nr_success = 1) then check solution */
01088   if (nr_success == 1){
01089 
01090     Vsq = calc_vsq (Bsq,Msq,Ssq,Q);
01091     Gsq = 1.0/(1.0-Vsq);
01092     Chi = calc_chi (d,Vsq,Gsq,Q);
01093     rho = d / sqrt(fabs(Gsq));
01094     pgas = Gamma_1*Chi/Gamma;
01095  
01096     tmp1 = 1.0 / Q;
01097     tmp2 = 1.0 / (Q + Bsq);
01098     Prim1D.d = rho;
01099     Prim1D.P = pgas;
01100     Prim1D.Vx = (U->Mx + S*(*Bx)*tmp1)*tmp2;
01101     Prim1D.Vy = (U->My + S*U->By*tmp1)*tmp2;
01102     Prim1D.Vz = (U->Mz + S*U->Bz*tmp1)*tmp2;
01103         
01104     Prim1D.By = U->By;
01105     Prim1D.Bz = U->Bz;
01106           
01107   } else {
01108     Prim1D.d = -1.0;
01109     Prim1D.P = -1.0;
01110     Prim1D.Vx = 1.0;
01111     Prim1D.Vy = 1.0;
01112     Prim1D.Vz = 1.0;
01113           
01114     Prim1D.By = U->By;
01115     Prim1D.Bz = U->Bz;
01116   }
01117           
01118   return Prim1D;
01119 }
01120 
01121 /*----------------------------------------------------------------------------*/
01122 /*! \fn Prim1DS vsq1D_fix(const Cons1DS *U, const Real *Bx) 
01123  *  \brief vsq1D_fix: SPECIAL RELATIVISTIC MHD VERSION  
01124  *
01125  *  - conserved variables = (d,Mx,My,Mz,[E],[By,Bz])
01126  *  - primitive variables = (d,Vx,Vy,Vz,[P],[By,Bz])
01127  *  - Bx is passed in through the argument list.
01128  */
01129 
01130 Prim1DS vsq1D_fix (const Cons1DS *U, const Real *Bx)
01131 {
01132   Prim1DS Prim1D;
01133   Cons1DS Con1D;
01134   int    k, done=0;
01135   Real Bsq = 0.0, Msq = 0.0, S = 0.0, Ssq = 0.0;
01136   Real v2, v2c, fc, f, W, dW, lor2, lor;
01137   Real fmin, fmax, v2min, v2max, p, d;
01138         
01139   d = 1.0e0;/*MAX(U->d,1.0e-2);*/
01140   p = 1.0e-1;/*0.000255457;*/
01141         
01142   /* Calculate Bsq = B^2, Msq = M^2, S = M \dot B, Ssq = S^2 */
01143   Bsq = SQR((*Bx)) + SQR(U->By) + SQR(U->Bz);
01144   Msq = SQR(U->Mx) + SQR(U->My) + SQR(U->Mz);
01145   S = U->Mx * (*Bx) + U->My * U->By + U->Mz * U->Bz;
01146   Ssq = SQR(S);
01147         
01148   v2max = 1.0-1.e-8;
01149   v2c = 0.95;
01150   FUNV2(d, v2c, p, Bsq, Msq, Ssq, &W, &fc);
01151   v2  = 0.96;
01152   for (k = 1; k < 100; k++){
01153     FUNV2(d, v2, p, Bsq, Msq, Ssq, &W, &f);
01154     if (done == 1) break;
01155     dW  = (v2 - v2c)/(f - fc)*f;
01156     v2c = v2; fc = f;
01157     v2 -= dW;
01158     v2 = MIN(v2max,v2);
01159     v2 = MAX(v2, 0.0);
01160     if (fabs(v2) < 1.e-9) done = 1;
01161     if (fabs(f) < 1.e-9) done = 1;
01162   }
01163   if (v2 > 1.0 || k >= 100) {
01164     ath_error("[fix_vsq1D]: too many iter while fixing p , v^2 = %f\n", v2);
01165   }
01166         
01167   lor2  = 1.0/(1.0 - v2);
01168   lor   = sqrt(lor2);
01169         
01170   Con1D = *U;
01171   Con1D.d = d;
01172   Con1D.E = W - p + 0.5*(1.0 + v2)*Bsq - 0.5*Ssq/(W*W);
01173   Prim1D = Cons1D_to_Prim1D(&Con1D, Bx);
01174   v2 = SQR(Prim1D.Vx) + SQR(Prim1D.Vy) + SQR(Prim1D.Vz);
01175         
01176   return Prim1D;
01177 }
01178 
01179 /*! \fn static Real calc_func (Real Q, Real E, Real Bsq, Real Ssq, 
01180  *                             Real Vsq, Real pgas)
01181  *  \brief Evaluate equation A25 for E, rather than E'
01182  */
01183 static Real calc_func (Real Q, Real E, Real Bsq, Real Ssq, Real Vsq, Real pgas)
01184 {
01185   Real tmp1;
01186   /* Evaluate equation A25 for E, rather than E' */
01187   return Q - pgas + 0.5*(1.0+Vsq)*Bsq - (0.5*Ssq/Q/Q) - E;
01188 }
01189 
01190 /*! \fn static Real calc_ent_func(Real rho, Real pgas, Real d, Real ent)
01191  *  \brief Evaluate equation A25 for E, rather than E' */
01192 static Real calc_ent_func (Real rho, Real pgas, Real d, Real ent)
01193 {
01194   Real tmp1;
01195   /* Evaluate equation A25 for E, rather than E' */
01196   tmp1 = pgas * pow(rho,-Gamma);
01197   return d*tmp1 - ent;
01198 }
01199 
01200 /*! \fn static Real calc_dfunc(Real Q, Real Bsq, Real Msq, Real Ssq, Real d, 
01201  *                             Real Vsq, Real Gsq, Real Chi)
01202  *  \brief Evaluate A8 for E & Q, rather than E', W' */
01203 static Real calc_dfunc(Real Q, Real Bsq, Real Msq, Real Ssq, Real d, Real Vsq,
01204                        Real Gsq, Real Chi)
01205 {
01206   /* Evaluate A8 for E & Q, rather than E', W' */
01207   Real dp_dQ, dp_dchi, dchi_dQ, dp_drho, drho_dQ, dVsq_dQ;
01208   Real G, Qsq, Qth, scrh1;
01209 
01210   Qsq = Q*Q;
01211   Qth = Qsq*Q;
01212   G = sqrt(fabs(Gsq));
01213 
01214   scrh1 = Q + Bsq;
01215   dVsq_dQ  = Ssq*(3.0*Q*scrh1 + Bsq*Bsq) + Msq*Qth;
01216   dVsq_dQ *= -2.0/Qth/(scrh1*scrh1*scrh1);
01217 
01218   /* -- kinematical terms -- */
01219   /* see eqn. A14 - A16 */
01220   dchi_dQ =  1.0 - Vsq - 0.5*G*(d + 2.0*Chi*G)*dVsq_dQ;
01221   drho_dQ = -0.5*d*G*dVsq_dQ;
01222 
01223   /* -- thermo terms, change for different EOS --*/
01224   /* see Section A2 of M&M */
01225   dp_dchi = (Gamma - 1.0)/Gamma;
01226   dp_drho = 0.0;
01227 
01228   dp_dQ = dp_dchi*dchi_dQ + dp_drho*drho_dQ;
01229 
01230   return 1.0 - dp_dQ +0.5*Bsq*dVsq_dQ + Ssq / Qth;
01231 }
01232 
01233 /*! \fn static Real calc_ent_dfunc(Real Q, Real Bsq, Real Msq, Real Ssq,
01234  *                                 Real d, Real Vsq, Real Gsq, Real Chi,
01235  *                                 Real pgas, Real rho)
01236  *  \brief Evaluate A8 for E & Q, rather than E', W' */
01237 static Real calc_ent_dfunc(Real Q, Real Bsq, Real Msq, Real Ssq,
01238                            Real d, Real Vsq, Real Gsq, Real Chi,
01239                            Real pgas, Real rho)
01240 {
01241   /* Evaluate A8 for E & Q, rather than E', W' */
01242   Real dp_dQ, dp_dchi, dchi_dQ, dp_drho, drho_dQ, dVsq_dQ;
01243   Real G, Qsq, Qth, scrh1;
01244 
01245   Qsq = Q*Q;
01246   Qth = Qsq*Q;
01247   G = sqrt(fabs(Gsq));
01248 
01249   scrh1 = Q + Bsq;
01250   dVsq_dQ  = Ssq*(3.0*Q*scrh1 + Bsq*Bsq) + Msq*Qth;
01251   dVsq_dQ *= -2.0/Qth/(scrh1*scrh1*scrh1);
01252 
01253   /* -- kinematical terms -- */
01254   /* see eqn. A14 - A16 */
01255   dchi_dQ =  1.0 - Vsq - 0.5*G*(d + 2.0*Chi*G)*dVsq_dQ;
01256   drho_dQ = -0.5*d*G*dVsq_dQ;
01257 
01258   /* -- thermo terms, change for different EOS --*/
01259   /* see Section A2 of M&M */
01260   dp_dchi = (Gamma - 1.0)/Gamma;
01261   dp_drho = 0.0;
01262 
01263   dp_dQ = dp_dchi*dchi_dQ + dp_drho*drho_dQ;
01264 
01265   return d*pow(rho,-Gamma)*dp_dQ - Gamma*pgas*pow(rho,Gamma+1.0)*drho_dQ;
01266 }
01267 
01268 /*! \fn static Real calc_vsq(Real Bsq, Real Msq, Real Ssq, Real Q)
01269  *  \brief Obtain v^2 via eqn. A3 */
01270 static Real calc_vsq (Real Bsq, Real Msq, Real Ssq, Real Q)
01271 {
01272   /* Obtain v^2 via eqn. A3 */
01273   Real Qsq, Ssq_Qsq, scrh1;
01274   Qsq = Q*Q;
01275   Ssq_Qsq = Ssq / Qsq;
01276   scrh1 = Q + Bsq;
01277   return (Msq + Ssq_Qsq*(scrh1 + Q))/(scrh1*scrh1);
01278 }
01279 
01280 /*! \fn static Real calc_chi (Real d, Real Vsq, Real Gsq, Real Q)
01281  *  \brief Evaluate \Chi from eqn. A11 using Q, rathher than Q' */
01282 static Real calc_chi (Real d, Real Vsq, Real Gsq, Real Q)
01283 {
01284   /* Evaluate \Chi from eqn. A11 using Q, rathher than Q' */
01285   Real G, tmp1, tmp2;
01286   G = sqrt(fabs(Gsq));
01287   tmp1 = 1.0 - Vsq;
01288   tmp2 = Q - d*G;
01289   return tmp2*tmp1;
01290 }
01291 
01292 void FUNV2(const Real d, const Real v2, 
01293            const Real p, const Real Bsq, const Real Msq, const Real Ssq,
01294            Real *W, Real *f)
01295 {
01296   Real lor2, lor, W2, pg;
01297         
01298   lor2 = 1.0/(1.0 - v2);
01299   lor  = sqrt(lor2);
01300   pg   = p*lor;
01301   *W  = (d + pg*Gamma/(Gamma - 1.0))*lor;
01302   W2    = SQR(*W);
01303         
01304   *f  =  Ssq*(2.0*(*W) + Bsq) + Msq*W2;
01305   *f /= ((*W) + Bsq)*((*W) + Bsq)*W2;
01306   *f -= v2;
01307 }
01308 #endif /* SPECIAL_RELATIVITY && MHD */