Athena: rsolvers/esystem_roe.c Source File

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file esystem_roe.c
00004  *  \brief Functions to evaluate the eigenvalues, and left- and
00005  * right-eigenvectors of "Roe's matrix A" for the linearized system in the 
00006  * CONSERVED variables.
00007  *
00008  * PURPOSE: Functions to evaluate the eigenvalues, and left- and
00009  * right-eigenvectors of "Roe's matrix A" for the linearized system in the 
00010  * CONSERVED variables, i.e. U,t = AU,x, where U=(d,d*vx,d*vy,d*vz,[E],[By,Bz]).
00011  * The eigenvalues are returned through the argument list as a vector of length
00012  * NWAVE.  The eigenvectors are returned as matrices of size (NWAVE)x(NWAVE),
00013  * with right-eigenvectors stored as COLUMNS (so R_i = right_eigenmatrix[*][i]),
00014  * and left-eigenvectors stored as ROWS (so L_i = left_eigenmatrix[i][*]).
00015  *
00016  * To improve performance components of the eigenvectors which are zero
00017  * are not set here (eigenmatrices must be initialized to zero in calling
00018  * routine).  However, for completeness, statements which set these values
00019  * are included but are commented out.
00020  *
00021  * If the input argument for the left- or right-eigenmatrix is set to the NULL
00022  * pointer, only the eigenvalues are returned
00023  *
00024  * The "Roe-averaging" of the L/R states must be performed in the calling funct
00025  *
00026  * REFERENCES:
00027  * - P. Cargo & G. Gallice, "Roe matrices for ideal MHD and systematic
00028  *   construction of Roe matrices for systems of conservation laws",
00029  *   JCP, 136, 446 (1997)
00030  *
00031  * - J. Stone, T. Gardiner, P. Teuben, J. Hawley, & J. Simon "Athena: A new
00032  *   code for astrophysical MHD", ApJS, (2008), Appendix B
00033  *   Equation numbers refer to this paper.
00034  *
00035  * CONTAINS PUBLIC FUNCTIONS:
00036  * - esys_roe_iso_hyd()
00037  * - esys_roe_adb_hyd()
00038  * - esys_roe_iso_mhd()
00039  * - esys_roe_adb_mhd()                                                       */
00040 /*============================================================================*/
00041 
00042 #include <stdlib.h>
00043 #include <math.h>
00044 #include <float.h>
00045 #include "../defs.h"
00046 #include "../athena.h"
00047 #include "../globals.h"
00048 #include "prototypes.h"
00049 #include "../prototypes.h"
00050 
00051 /*----------------------------------------------------------------------------*/
00052 /*! \fn  void esys_roe_iso_hyd(const Real v1, const Real v2, const Real v3,
00053  *  Real eigenvalues[],
00054  *  Real right_eigenmatrix[][4], Real left_eigenmatrix[][4])
00055  *  \brief ISOTHERMAL HYDRO
00056  *
00057  *  - Input: v1,v2,v3 = Roe averaged components of velocity
00058  *  - Output: eigenvalues[4],right_eigenmatrix[4,4],left_eigenmatrix[4,4]
00059  */
00060 
00061 #if defined(ISOTHERMAL) && defined(HYDRO)
00062 void esys_roe_iso_hyd(const Real v1, const Real v2, const Real v3,
00063   Real eigenvalues[],
00064   Real right_eigenmatrix[][4], Real left_eigenmatrix[][4])
00065 {
00066 
00067 /* Compute eigenvalues (eq. B6) */
00068 
00069   eigenvalues[0] = v1 - Iso_csound;
00070   eigenvalues[1] = v1;
00071   eigenvalues[2] = v1;
00072   eigenvalues[3] = v1 + Iso_csound;
00073   if (right_eigenmatrix == NULL || left_eigenmatrix == NULL) return;
00074 
00075 /* Right-eigenvectors, stored as COLUMNS (eq. B3) */
00076 
00077   right_eigenmatrix[0][0] = 1.0;
00078   right_eigenmatrix[1][0] = v1 - Iso_csound;
00079   right_eigenmatrix[2][0] = v2;
00080   right_eigenmatrix[3][0] = v3;
00081 
00082 /*right_eigenmatrix[0][1] = 0.0; */
00083 /*right_eigenmatrix[1][1] = 0.0; */
00084   right_eigenmatrix[2][1] = 1.0;
00085 /*right_eigenmatrix[3][1] = 0.0; */
00086 
00087 /*right_eigenmatrix[0][2] = 0.0; */
00088 /*right_eigenmatrix[1][2] = 0.0; */
00089 /*right_eigenmatrix[2][2] = 0.0; */
00090   right_eigenmatrix[3][2] = 1.0;
00091 
00092   right_eigenmatrix[0][3] = 1.0;
00093   right_eigenmatrix[1][3] = v1 + Iso_csound;
00094   right_eigenmatrix[2][3] = v2;
00095   right_eigenmatrix[3][3] = v3;
00096 
00097 /* Left-eigenvectors, stored as ROWS (eq. B7) */
00098 
00099   left_eigenmatrix[0][0] = 0.5*(1.0 + v1/Iso_csound);
00100   left_eigenmatrix[0][1] = -0.5/Iso_csound;
00101 /*left_eigenmatrix[0][2] = 0.0; */
00102 /*left_eigenmatrix[0][3] = 0.0; */
00103 
00104   left_eigenmatrix[1][0] = -v2;
00105 /*left_eigenmatrix[1][1] = 0.0; */
00106   left_eigenmatrix[1][2] = 1.0;
00107 /*left_eigenmatrix[1][3] = 0.0; */
00108 
00109   left_eigenmatrix[2][0] = -v3;
00110 /*left_eigenmatrix[2][1] = 0.0; */
00111 /*left_eigenmatrix[2][2] = 0.0; */
00112   left_eigenmatrix[2][3] = 1.0;
00113 
00114   left_eigenmatrix[3][0] = 0.5*(1.0 - v1/Iso_csound);
00115   left_eigenmatrix[3][1] = 0.5/Iso_csound;
00116 /*left_eigenmatrix[3][2] = 0.0; */
00117 /*left_eigenmatrix[3][3] = 0.0; */
00118 }
00119 #endif
00120 
00121 /*----------------------------------------------------------------------------*/
00122 /*! \fn void esys_roe_adb_hyd(const Real v1, const Real v2, const Real v3, 
00123  *                            const Real h, Real right_eigenmatrix[][5], 
00124  *                            Real left_eigenmatrix[][5])
00125  *  \brief ADIABATIC HYDRO
00126  *
00127  * - Input: v1,v2,v3,h = Roe averaged velocities and enthalpy
00128  * - Output: eigenvalues[5], right_eigenmatrix[5,5], left_eigenmatrix[5,5];
00129  */
00130 
00131 #if defined(ADIABATIC) && defined(HYDRO)
00132 void esys_roe_adb_hyd(const Real v1, const Real v2, const Real v3, const Real h,
00133   Real eigenvalues[],
00134   Real right_eigenmatrix[][5], Real left_eigenmatrix[][5])
00135 {
00136   Real vsq,asq,a,na,qa;
00137   vsq = v1*v1 + v2*v2 + v3*v3;
00138   asq = Gamma_1*MAX((h-0.5*vsq), TINY_NUMBER);
00139   a = sqrt(asq);
00140 
00141 /* Compute eigenvalues (eq. B2) */
00142 
00143   eigenvalues[0] = v1 - a;
00144   eigenvalues[1] = v1;
00145   eigenvalues[2] = v1;
00146   eigenvalues[3] = v1;
00147   eigenvalues[4] = v1 + a;
00148   if (right_eigenmatrix == NULL || left_eigenmatrix == NULL) return;
00149 
00150 /* Right-eigenvectors, stored as COLUMNS (eq. B3) */
00151 
00152   right_eigenmatrix[0][0] = 1.0;
00153   right_eigenmatrix[1][0] = v1 - a;
00154   right_eigenmatrix[2][0] = v2;
00155   right_eigenmatrix[3][0] = v3;
00156   right_eigenmatrix[4][0] = h - v1*a;
00157 
00158 /*right_eigenmatrix[0][1] = 0.0; */
00159 /*right_eigenmatrix[1][1] = 0.0; */
00160   right_eigenmatrix[2][1] = 1.0;
00161 /*right_eigenmatrix[3][1] = 0.0; */
00162   right_eigenmatrix[4][1] = v2;
00163 
00164 /*right_eigenmatrix[0][2] = 0.0; */
00165 /*right_eigenmatrix[1][2] = 0.0; */
00166 /*right_eigenmatrix[2][2] = 0.0; */
00167   right_eigenmatrix[3][2] = 1.0;
00168   right_eigenmatrix[4][2] = v3;
00169 
00170   right_eigenmatrix[0][3] = 1.0;
00171   right_eigenmatrix[1][3] = v1;
00172   right_eigenmatrix[2][3] = v2;
00173   right_eigenmatrix[3][3] = v3;
00174   right_eigenmatrix[4][3] = 0.5*vsq;
00175 
00176   right_eigenmatrix[0][4] = 1.0;
00177   right_eigenmatrix[1][4] = v1 + a;
00178   right_eigenmatrix[2][4] = v2;
00179   right_eigenmatrix[3][4] = v3;
00180   right_eigenmatrix[4][4] = h + v1*a;
00181 
00182 /* Left-eigenvectors, stored as ROWS (eq. B4) */
00183 
00184   na = 0.5/asq;
00185   left_eigenmatrix[0][0] = na*(0.5*Gamma_1*vsq + v1*a);
00186   left_eigenmatrix[0][1] = -na*(Gamma_1*v1 + a);
00187   left_eigenmatrix[0][2] = -na*Gamma_1*v2;
00188   left_eigenmatrix[0][3] = -na*Gamma_1*v3;
00189   left_eigenmatrix[0][4] = na*Gamma_1;
00190 
00191   left_eigenmatrix[1][0] = -v2;
00192 /*left_eigenmatrix[1][1] = 0.0; */
00193   left_eigenmatrix[1][2] = 1.0;
00194 /*left_eigenmatrix[1][3] = 0.0; */
00195 /*left_eigenmatrix[1][4] = 0.0; */
00196 
00197   left_eigenmatrix[2][0] = -v3;
00198 /*left_eigenmatrix[2][1] = 0.0; */
00199 /*left_eigenmatrix[2][2] = 0.0; */
00200   left_eigenmatrix[2][3] = 1.0;
00201 /*left_eigenmatrix[2][4] = 0.0; */
00202 
00203   qa = Gamma_1/asq;
00204   left_eigenmatrix[3][0] = 1.0 - na*Gamma_1*vsq;
00205   left_eigenmatrix[3][1] = qa*v1;
00206   left_eigenmatrix[3][2] = qa*v2;
00207   left_eigenmatrix[3][3] = qa*v3;
00208   left_eigenmatrix[3][4] = -qa;
00209 
00210   left_eigenmatrix[4][0] = na*(0.5*Gamma_1*vsq - v1*a);
00211   left_eigenmatrix[4][1] = -na*(Gamma_1*v1 - a);
00212   left_eigenmatrix[4][2] = left_eigenmatrix[0][2];
00213   left_eigenmatrix[4][3] = left_eigenmatrix[0][3];
00214   left_eigenmatrix[4][4] = left_eigenmatrix[0][4];
00215 }
00216 #endif
00217 
00218 /*----------------------------------------------------------------------------*/
00219 /*! \fn void esys_roe_iso_mhd(const Real d, const Real v1, const Real v2, 
00220  *                      const Real v3, const Real b1, const Real b2, 
00221  *                      const Real b3, const Real x, const Real y, 
00222  *                      Real eigenvalues[],
00223  *                      Real right_eigenmatrix[][6], Real left_eigenmatrix[][6])
00224  *  \brief ISOTHERMAL MHD
00225  *
00226  * - Input: d,v1,v2,v3,b1,b2,b3 = Roe averaged density, velocities, and B field
00227  *          x,y = numerical factors (eqs. )
00228  * - Output: eigenvalues[6], right_eigenmatrix[6,6], left_eigenmatrix[6,6];
00229  */
00230 
00231 #if defined(ISOTHERMAL) && defined(MHD)
00232 void esys_roe_iso_mhd(const Real d, const Real v1, const Real v2, const Real v3,
00233   const Real b1, const Real b2, const Real b3, const Real x, const Real y, 
00234   Real eigenvalues[],
00235   Real right_eigenmatrix[][6], Real left_eigenmatrix[][6])
00236 {
00237   Real btsq,bt_starsq,vaxsq,twid_csq,cfsq,cf,cssq,cs;
00238   Real bt,bt_star,bet2,bet3,bet2_star,bet3_star,bet_starsq,alpha_f,alpha_s;
00239   Real sqrtd,s,twid_c,qf,qs,af_prime,as_prime,vax;
00240   Real norm,cff,css,af,as,afpb,aspb,q2_star,q3_star,vqstr;
00241   Real ct2,tsum,tdif,cf2_cs2;
00242   Real di = 1.0/d;
00243   btsq = b2*b2 + b3*b3;
00244   bt_starsq = btsq*y;
00245   vaxsq = b1*b1*di;
00246   twid_csq = Iso_csound2 + x;
00247 
00248 /* Compute fast- and slow-magnetosonic speeds (eq. B39) */
00249 
00250   ct2 = bt_starsq*di;
00251   tsum = vaxsq + ct2 + twid_csq;
00252   tdif = vaxsq + ct2 - twid_csq;
00253   cf2_cs2 = sqrt((double)(tdif*tdif + 4.0*twid_csq*ct2));
00254 
00255   cfsq = 0.5*(tsum + cf2_cs2);
00256   cf = sqrt((double)cfsq);
00257 
00258   cssq = twid_csq*vaxsq/cfsq;
00259   cs = sqrt((double)cssq);
00260 
00261 /* Compute beta's (eqs. A17, B28, B40) */
00262 
00263   bt = sqrt(btsq);
00264   bt_star = sqrt(bt_starsq);
00265   if (bt == 0.0) {
00266     bet2 = 1.0;
00267     bet3 = 0.0;
00268   } 
00269   else {
00270     bet2 = b2/bt;
00271     bet3 = b3/bt;
00272   }
00273   bet2_star = bet2/sqrt(y);
00274   bet3_star = bet3/sqrt(y);
00275   bet_starsq = bet2_star*bet2_star + bet3_star*bet3_star;
00276 
00277 /* Compute alpha's (eq. A16) */
00278 
00279   if ((cfsq-cssq) == 0.0) {
00280     alpha_f = 1.0;
00281     alpha_s = 0.0;
00282   } else if ((twid_csq - cssq) <= 0.0) {
00283     alpha_f = 0.0;
00284     alpha_s = 1.0;
00285   } else if ((cfsq - twid_csq) <= 0.0) {
00286     alpha_f = 1.0;
00287     alpha_s = 0.0;
00288   } else {
00289     alpha_f = sqrt((twid_csq - cssq)/(cfsq - cssq));
00290     alpha_s = sqrt((cfsq - twid_csq)/(cfsq - cssq));
00291   }
00292 
00293 /* Compute Q's (eq. A14-15), etc. */
00294 
00295   sqrtd = sqrt(d);
00296   s = SIGN(b1);
00297   twid_c = sqrt(twid_csq);
00298   qf = cf*alpha_f*s;
00299   qs = cs*alpha_s*s;
00300   af_prime = twid_c*alpha_f/sqrtd;
00301   as_prime = twid_c*alpha_s/sqrtd;
00302 
00303 /* Compute eigenvalues (eq. B38) */
00304 
00305   vax  = sqrt(vaxsq);
00306   eigenvalues[0] = v1 - cf;
00307   eigenvalues[1] = v1 - vax;
00308   eigenvalues[2] = v1 - cs;
00309   eigenvalues[3] = v1 + cs;
00310   eigenvalues[4] = v1 + vax;
00311   eigenvalues[5] = v1 + cf;
00312   if (right_eigenmatrix == NULL || left_eigenmatrix == NULL) return;
00313 
00314 /* Right-eigenvectors, stored as COLUMNS (eq. B21) */
00315 
00316   right_eigenmatrix[0][0] = alpha_f;
00317   right_eigenmatrix[1][0] = alpha_f*(v1 - cf);
00318   right_eigenmatrix[2][0] = alpha_f*v2 + qs*bet2_star;
00319   right_eigenmatrix[3][0] = alpha_f*v3 + qs*bet3_star;
00320   right_eigenmatrix[4][0] = as_prime*bet2_star;
00321   right_eigenmatrix[5][0] = as_prime*bet3_star;
00322 
00323 /*right_eigenmatrix[0][1] = 0.0; */
00324 /*right_eigenmatrix[1][1] = 0.0; */
00325   right_eigenmatrix[2][1] = -bet3;
00326   right_eigenmatrix[3][1] = bet2;
00327   right_eigenmatrix[4][1] = -bet3*s/sqrtd;
00328   right_eigenmatrix[5][1] = bet2*s/sqrtd;
00329 
00330   right_eigenmatrix[0][2] = alpha_s;
00331   right_eigenmatrix[1][2] = alpha_s*(v1 - cs);
00332   right_eigenmatrix[2][2] = alpha_s*v2 - qf*bet2_star;
00333   right_eigenmatrix[3][2] = alpha_s*v3 - qf*bet3_star;
00334   right_eigenmatrix[4][2] = -af_prime*bet2_star;
00335   right_eigenmatrix[5][2] = -af_prime*bet3_star;
00336 
00337   right_eigenmatrix[0][3] = alpha_s;
00338   right_eigenmatrix[1][3] = alpha_s*(v1 + cs);
00339   right_eigenmatrix[2][3] = alpha_s*v2 + qf*bet2_star;
00340   right_eigenmatrix[3][3] = alpha_s*v3 + qf*bet3_star;
00341   right_eigenmatrix[4][3] = right_eigenmatrix[4][2];
00342   right_eigenmatrix[5][3] = right_eigenmatrix[5][2];
00343 
00344 /*right_eigenmatrix[0][4] = 0.0; */
00345 /*right_eigenmatrix[1][4] = 0.0; */
00346   right_eigenmatrix[2][4] = bet3;
00347   right_eigenmatrix[3][4] = -bet2;
00348   right_eigenmatrix[4][4] = right_eigenmatrix[4][1];
00349   right_eigenmatrix[5][4] = right_eigenmatrix[5][1];
00350 
00351   right_eigenmatrix[0][5] = alpha_f;
00352   right_eigenmatrix[1][5] = alpha_f*(v1 + cf);
00353   right_eigenmatrix[2][5] = alpha_f*v2 - qs*bet2_star;
00354   right_eigenmatrix[3][5] = alpha_f*v3 - qs*bet3_star;
00355   right_eigenmatrix[4][5] = right_eigenmatrix[4][0];
00356   right_eigenmatrix[5][5] = right_eigenmatrix[5][0];
00357 
00358 /* Left-eigenvectors, stored as ROWS (eq. B41) */
00359 /* Normalize by 1/2a^{2}: quantities denoted by \hat{f} */
00360 
00361   norm = 0.5/twid_csq;
00362   cff = norm*alpha_f*cf;
00363   css = norm*alpha_s*cs;
00364   qf *= norm;
00365   qs *= norm;
00366   af = norm*af_prime*d;
00367   as = norm*as_prime*d;
00368   afpb = norm*af_prime*bt_star;
00369   aspb = norm*as_prime*bt_star;
00370 
00371   q2_star = bet2_star/bet_starsq;
00372   q3_star = bet3_star/bet_starsq;
00373   vqstr = (v2*q2_star + v3*q3_star);
00374 
00375   left_eigenmatrix[0][0] = cff*(cf+v1) - qs*vqstr - aspb;
00376   left_eigenmatrix[0][1] = -cff;
00377   left_eigenmatrix[0][2] = qs*q2_star;
00378   left_eigenmatrix[0][3] = qs*q3_star;
00379   left_eigenmatrix[0][4] = as*q2_star;
00380   left_eigenmatrix[0][5] = as*q3_star;
00381 
00382   left_eigenmatrix[1][0] = 0.5*(v2*bet3 - v3*bet2);
00383 /*left_eigenmatrix[1][1] = 0.0; */
00384   left_eigenmatrix[1][2] = -0.5*bet3;
00385   left_eigenmatrix[1][3] = 0.5*bet2;
00386   left_eigenmatrix[1][4] = -0.5*sqrtd*bet3*s;
00387   left_eigenmatrix[1][5] = 0.5*sqrtd*bet2*s;
00388 
00389   left_eigenmatrix[2][0] = css*(cs+v1) + qf*vqstr + afpb;
00390   left_eigenmatrix[2][1] = -css;
00391   left_eigenmatrix[2][2] = -qf*q2_star;
00392   left_eigenmatrix[2][3] = -qf*q3_star;
00393   left_eigenmatrix[2][4] = -af*q2_star;
00394   left_eigenmatrix[2][5] = -af*q3_star;
00395 
00396   left_eigenmatrix[3][0] = css*(cs-v1) - qf*vqstr + afpb;
00397   left_eigenmatrix[3][1] = css;
00398   left_eigenmatrix[3][2] = -left_eigenmatrix[2][2];
00399   left_eigenmatrix[3][3] = -left_eigenmatrix[2][3];
00400   left_eigenmatrix[3][4] = left_eigenmatrix[2][4];
00401   left_eigenmatrix[3][5] = left_eigenmatrix[2][5];
00402 
00403   left_eigenmatrix[4][0] = -left_eigenmatrix[1][0];
00404 /*left_eigenmatrix[4][1] = 0.0; */
00405   left_eigenmatrix[4][2] = -left_eigenmatrix[1][2];
00406   left_eigenmatrix[4][3] = -left_eigenmatrix[1][3];
00407   left_eigenmatrix[4][4] = left_eigenmatrix[1][4];
00408   left_eigenmatrix[4][5] = left_eigenmatrix[1][5];
00409 
00410   left_eigenmatrix[5][0] = cff*(cf-v1) + qs*vqstr - aspb;
00411   left_eigenmatrix[5][1] = cff;
00412   left_eigenmatrix[5][2] = -left_eigenmatrix[0][2];
00413   left_eigenmatrix[5][3] = -left_eigenmatrix[0][3];
00414   left_eigenmatrix[5][4] = left_eigenmatrix[0][4];
00415   left_eigenmatrix[5][5] = left_eigenmatrix[0][5];
00416 }
00417 #endif
00418 
00419 /*----------------------------------------------------------------------------*/
00420 /*! \fn void esys_roe_adb_mhd(const Real d, const Real v1, const Real v2, 
00421  *                      const Real v3, const Real h, const Real b1, 
00422  *                      const Real b2, const Real b3, Real eigenvalues[], 
00423  *                      Real right_eigenmatrix[][7], Real left_eigenmatrix[][7])
00424  *  \brief ADIABATIC MHD
00425  *
00426  * - Input: d,v1,v2,v3,h,b1,b2,b3=Roe averaged density, velocities, enthalpy, B
00427  *          x,y = numerical factors (see eqn XX)
00428  * - Output: eigenvalues[7], right_eigenmatrix[7,7], left_eigenmatrix[7,7];
00429  */
00430 
00431 #if defined(ADIABATIC) && defined(MHD)
00432 void esys_roe_adb_mhd(const Real d, const Real v1, const Real v2, const Real v3,
00433   const Real h, const Real b1, const Real b2, const Real b3, 
00434   const Real x, const Real y,
00435   Real eigenvalues[],
00436   Real right_eigenmatrix[][7], Real left_eigenmatrix[][7])
00437 {
00438   Real di,vsq,btsq,bt_starsq,vaxsq,hp,twid_asq,cfsq,cf,cssq,cs;
00439   Real bt,bt_star,bet2,bet3,bet2_star,bet3_star,bet_starsq,vbet,alpha_f,alpha_s;
00440   Real isqrtd,sqrtd,s,twid_a,qf,qs,af_prime,as_prime,afpbb,aspbb,vax;
00441   Real norm,cff,css,af,as,afpb,aspb,q2_star,q3_star,vqstr;
00442   Real ct2,tsum,tdif,cf2_cs2;
00443   Real qa,qb,qc,qd;
00444   di = 1.0/d;
00445   vsq = v1*v1 + v2*v2 + v3*v3;
00446   btsq = b2*b2 + b3*b3;
00447   bt_starsq = (Gamma_1 - Gamma_2*y)*btsq;
00448   vaxsq = b1*b1*di;
00449   hp = h - (vaxsq + btsq*di);
00450   twid_asq = MAX((Gamma_1*(hp-0.5*vsq)-Gamma_2*x), TINY_NUMBER);
00451 
00452 /* Compute fast- and slow-magnetosonic speeds (eq. B18) */
00453 
00454   ct2 = bt_starsq*di;
00455   tsum = vaxsq + ct2 + twid_asq;
00456   tdif = vaxsq + ct2 - twid_asq;
00457   cf2_cs2 = sqrt((double)(tdif*tdif + 4.0*twid_asq*ct2));
00458 
00459   cfsq = 0.5*(tsum + cf2_cs2);
00460   cf = sqrt((double)cfsq);
00461 
00462   cssq = twid_asq*vaxsq/cfsq;
00463   cs = sqrt((double)cssq);
00464 
00465 /* Compute beta(s) (eqs. A17, B20, B28) */
00466 
00467   bt = sqrt(btsq);
00468   bt_star = sqrt(bt_starsq);
00469   if (bt == 0.0) {
00470     bet2 = 1.0;
00471     bet3 = 0.0;
00472   } else {
00473     bet2 = b2/bt;
00474     bet3 = b3/bt;
00475   }
00476   bet2_star = bet2/sqrt(Gamma_1 - Gamma_2*y);
00477   bet3_star = bet3/sqrt(Gamma_1 - Gamma_2*y);
00478   bet_starsq = bet2_star*bet2_star + bet3_star*bet3_star;
00479   vbet = v2*bet2_star + v3*bet3_star;
00480 
00481 /* Compute alpha(s) (eq. A16) */
00482 
00483   if ((cfsq-cssq) == 0.0) {
00484     alpha_f = 1.0;
00485     alpha_s = 0.0;
00486   } else if ( (twid_asq - cssq) <= 0.0) {
00487     alpha_f = 0.0;
00488     alpha_s = 1.0;
00489   } else if ( (cfsq - twid_asq) <= 0.0) {
00490     alpha_f = 1.0;
00491     alpha_s = 0.0;
00492   } else {
00493     alpha_f = sqrt((twid_asq - cssq)/(cfsq - cssq));
00494     alpha_s = sqrt((cfsq - twid_asq)/(cfsq - cssq));
00495   }
00496 
00497 /* Compute Q(s) and A(s) (eq. A14-15), etc. */
00498 
00499   sqrtd = sqrt(d);
00500   isqrtd = 1.0/sqrtd;
00501   s = SIGN(b1);
00502   twid_a = sqrt(twid_asq);
00503   qf = cf*alpha_f*s;
00504   qs = cs*alpha_s*s;
00505   af_prime = twid_a*alpha_f*isqrtd;
00506   as_prime = twid_a*alpha_s*isqrtd;
00507   afpbb = af_prime*bt_star*bet_starsq;
00508   aspbb = as_prime*bt_star*bet_starsq;
00509 
00510 /* Compute eigenvalues (eq. B17) */
00511 
00512   vax = sqrt(vaxsq);
00513   eigenvalues[0] = v1 - cf;
00514   eigenvalues[1] = v1 - vax;
00515   eigenvalues[2] = v1 - cs;
00516   eigenvalues[3] = v1;
00517   eigenvalues[4] = v1 + cs;
00518   eigenvalues[5] = v1 + vax;
00519   eigenvalues[6] = v1 + cf;
00520   if (right_eigenmatrix == NULL || left_eigenmatrix == NULL) return;
00521 
00522 /* Right-eigenvectors, stored as COLUMNS (eq. B21) */
00523 /* Note statements are grouped in ROWS for optimization, even though rem[*][n]
00524  * is the nth right eigenvector */
00525 
00526   right_eigenmatrix[0][0] = alpha_f;
00527 /*right_eigenmatrix[0][1] = 0.0; */
00528   right_eigenmatrix[0][2] = alpha_s;
00529   right_eigenmatrix[0][3] = 1.0;
00530   right_eigenmatrix[0][4] = alpha_s;
00531 /*right_eigenmatrix[0][5] = 0.0; */
00532   right_eigenmatrix[0][6] = alpha_f;
00533 
00534   right_eigenmatrix[1][0] = alpha_f*eigenvalues[0];
00535 /*right_eigenmatrix[1][1] = 0.0; */
00536   right_eigenmatrix[1][2] = alpha_s*eigenvalues[2];
00537   right_eigenmatrix[1][3] = v1;
00538   right_eigenmatrix[1][4] = alpha_s*eigenvalues[4];
00539 /*right_eigenmatrix[1][5] = 0.0; */
00540   right_eigenmatrix[1][6] = alpha_f*eigenvalues[6];
00541 
00542   qa = alpha_f*v2;
00543   qb = alpha_s*v2;
00544   qc = qs*bet2_star;
00545   qd = qf*bet2_star;
00546   right_eigenmatrix[2][0] = qa + qc;
00547   right_eigenmatrix[2][1] = -bet3;
00548   right_eigenmatrix[2][2] = qb - qd;
00549   right_eigenmatrix[2][3] = v2;
00550   right_eigenmatrix[2][4] = qb + qd;
00551   right_eigenmatrix[2][5] = bet3;
00552   right_eigenmatrix[2][6] = qa - qc;
00553 
00554   qa = alpha_f*v3;
00555   qb = alpha_s*v3;
00556   qc = qs*bet3_star;
00557   qd = qf*bet3_star;
00558   right_eigenmatrix[3][0] = qa + qc;
00559   right_eigenmatrix[3][1] = bet2;
00560   right_eigenmatrix[3][2] = qb - qd;
00561   right_eigenmatrix[3][3] = v3;
00562   right_eigenmatrix[3][4] = qb + qd;
00563   right_eigenmatrix[3][5] = -bet2;
00564   right_eigenmatrix[3][6] = qa - qc;
00565 
00566   right_eigenmatrix[4][0] = alpha_f*(hp - v1*cf) + qs*vbet + aspbb;
00567   right_eigenmatrix[4][1] = -(v2*bet3 - v3*bet2);
00568   right_eigenmatrix[4][2] = alpha_s*(hp - v1*cs) - qf*vbet - afpbb;
00569   right_eigenmatrix[4][3] = 0.5*vsq + Gamma_2*x/Gamma_1;
00570   right_eigenmatrix[4][4] = alpha_s*(hp + v1*cs) + qf*vbet - afpbb;
00571   right_eigenmatrix[4][5] = -right_eigenmatrix[4][1];
00572   right_eigenmatrix[4][6] = alpha_f*(hp + v1*cf) - qs*vbet + aspbb;
00573 
00574   right_eigenmatrix[5][0] = as_prime*bet2_star;
00575   right_eigenmatrix[5][1] = -bet3*s*isqrtd;
00576   right_eigenmatrix[5][2] = -af_prime*bet2_star;
00577 /*right_eigenmatrix[5][3] = 0.0; */
00578   right_eigenmatrix[5][4] = right_eigenmatrix[5][2];
00579   right_eigenmatrix[5][5] = right_eigenmatrix[5][1];
00580   right_eigenmatrix[5][6] = right_eigenmatrix[5][0];
00581 
00582   right_eigenmatrix[6][0] = as_prime*bet3_star;
00583   right_eigenmatrix[6][1] = bet2*s*isqrtd;
00584   right_eigenmatrix[6][2] = -af_prime*bet3_star;
00585 /*right_eigenmatrix[6][3] = 0.0; */
00586   right_eigenmatrix[6][4] = right_eigenmatrix[6][2];
00587   right_eigenmatrix[6][5] = right_eigenmatrix[6][1];
00588   right_eigenmatrix[6][6] = right_eigenmatrix[6][0];
00589 
00590 /* Left-eigenvectors, stored as ROWS (eq. B29) */
00591 
00592 /* Normalize by 1/2a^{2}: quantities denoted by \hat{f} */
00593   norm = 0.5/twid_asq;
00594   cff = norm*alpha_f*cf;
00595   css = norm*alpha_s*cs;
00596   qf *= norm;
00597   qs *= norm;
00598   af = norm*af_prime*d;
00599   as = norm*as_prime*d;
00600   afpb = norm*af_prime*bt_star;
00601   aspb = norm*as_prime*bt_star;
00602 
00603 /* Normalize by (gamma-1)/2a^{2}: quantities denoted by \bar{f} */
00604   norm *= Gamma_1;
00605   alpha_f *= norm;
00606   alpha_s *= norm;
00607   q2_star = bet2_star/bet_starsq;
00608   q3_star = bet3_star/bet_starsq;
00609   vqstr = (v2*q2_star + v3*q3_star);
00610   norm *= 2.0;
00611 
00612   left_eigenmatrix[0][0] = alpha_f*(vsq-hp) + cff*(cf+v1) - qs*vqstr - aspb;
00613   left_eigenmatrix[0][1] = -alpha_f*v1 - cff;
00614   left_eigenmatrix[0][2] = -alpha_f*v2 + qs*q2_star;
00615   left_eigenmatrix[0][3] = -alpha_f*v3 + qs*q3_star;
00616   left_eigenmatrix[0][4] = alpha_f;
00617   left_eigenmatrix[0][5] = as*q2_star - alpha_f*b2;
00618   left_eigenmatrix[0][6] = as*q3_star - alpha_f*b3;
00619 
00620   left_eigenmatrix[1][0] = 0.5*(v2*bet3 - v3*bet2);
00621 /*left_eigenmatrix[1][1] = 0.0; */
00622   left_eigenmatrix[1][2] = -0.5*bet3;
00623   left_eigenmatrix[1][3] = 0.5*bet2;
00624 /*left_eigenmatrix[1][4] = 0.0; */
00625   left_eigenmatrix[1][5] = -0.5*sqrtd*bet3*s;
00626   left_eigenmatrix[1][6] = 0.5*sqrtd*bet2*s;
00627 
00628   left_eigenmatrix[2][0] = alpha_s*(vsq-hp) + css*(cs+v1) + qf*vqstr + afpb;
00629   left_eigenmatrix[2][1] = -alpha_s*v1 - css;
00630   left_eigenmatrix[2][2] = -alpha_s*v2 - qf*q2_star;
00631   left_eigenmatrix[2][3] = -alpha_s*v3 - qf*q3_star;
00632   left_eigenmatrix[2][4] = alpha_s;
00633   left_eigenmatrix[2][5] = -af*q2_star - alpha_s*b2;
00634   left_eigenmatrix[2][6] = -af*q3_star - alpha_s*b3;
00635 
00636   left_eigenmatrix[3][0] = 1.0 - norm*(0.5*vsq - Gamma_2*x/Gamma_1); 
00637   left_eigenmatrix[3][1] = norm*v1;
00638   left_eigenmatrix[3][2] = norm*v2;
00639   left_eigenmatrix[3][3] = norm*v3;
00640   left_eigenmatrix[3][4] = -norm;
00641   left_eigenmatrix[3][5] = norm*b2;
00642   left_eigenmatrix[3][6] = norm*b3;
00643 
00644   left_eigenmatrix[4][0] = alpha_s*(vsq-hp) + css*(cs-v1) - qf*vqstr + afpb;
00645   left_eigenmatrix[4][1] = -alpha_s*v1 + css;
00646   left_eigenmatrix[4][2] = -alpha_s*v2 + qf*q2_star;
00647   left_eigenmatrix[4][3] = -alpha_s*v3 + qf*q3_star;
00648   left_eigenmatrix[4][4] = alpha_s;
00649   left_eigenmatrix[4][5] = left_eigenmatrix[2][5];
00650   left_eigenmatrix[4][6] = left_eigenmatrix[2][6];
00651 
00652   left_eigenmatrix[5][0] = -left_eigenmatrix[1][0];
00653 /*left_eigenmatrix[5][1] = 0.0; */
00654   left_eigenmatrix[5][2] = -left_eigenmatrix[1][2];
00655   left_eigenmatrix[5][3] = -left_eigenmatrix[1][3];
00656 /*left_eigenmatrix[5][4] = 0.0; */
00657   left_eigenmatrix[5][5] = left_eigenmatrix[1][5];
00658   left_eigenmatrix[5][6] = left_eigenmatrix[1][6];
00659 
00660   left_eigenmatrix[6][0] = alpha_f*(vsq-hp) + cff*(cf-v1) + qs*vqstr - aspb;
00661   left_eigenmatrix[6][1] = -alpha_f*v1 + cff;
00662   left_eigenmatrix[6][2] = -alpha_f*v2 - qs*q2_star;
00663   left_eigenmatrix[6][3] = -alpha_f*v3 - qs*q3_star;
00664   left_eigenmatrix[6][4] = alpha_f;
00665   left_eigenmatrix[6][5] = left_eigenmatrix[0][5];
00666   left_eigenmatrix[6][6] = left_eigenmatrix[0][6];
00667 }
00668 #endif