Athena: reconstruction/esystem

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file esystem_prim.c
00004  *  \brief Functions to evaluate the eigenvalues, and left- and
00005  * right-eigenvectors for the linearized system in the 
00006  * PRIMITIVE variables.
00007  *
00008  * PURPOSE: Functions to evaluate the eigenvalues, and left- and
00009  * right-eigenvectors for the linearized system in the 
00010  * PRIMITIVE variables, i.e. W,t = AW,x, where W=(d,vx,vy,vz,[P],[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 (R_i[*] = right_eigenmatrix[*][i]),
00014  * and left-eigenvectors stored as ROWS (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  * REFERENCES:
00022  * - J. Stone, T. Gardiner, P. Teuben, J. Hawley, & J. Simon "Athena: A new
00023  *   code for astrophysical MHD", ApJS, (2008), Appendix A
00024  *   Equation numbers refer to this paper.
00025  *
00026  * CONTAINS PUBLIC FUNCTIONS:
00027  * - esys_prim_iso_hyd() - isothermal hydrodynamics
00028  * - esys_prim_adb_hyd() - adiabatic hydrodynamics
00029  * - esys_prim_iso_mhd() - isothermal MHD
00030  * - esys_prim_adb_mhd() - adiabatic MHD
00031  *============================================================================*/
00032 
00033 #include <math.h>
00034 #include <float.h>
00035 #include "../defs.h"
00036 #include "../athena.h"
00037 #include "../globals.h"
00038 #include "prototypes.h"
00039 #include "../prototypes.h"
00040 
00041 /*----------------------------------------------------------------------------*/
00042 /*! \fn void esys_prim_iso_hyd(const Real d, const Real v1, 
00043  *  Real eigenvalues[],
00044  *  Real right_eigenmatrix[][4], Real left_eigenmatrix[][4])
00045  *  \brief ISOTHERMAL HYDRO
00046  *   Input: d, v1 = primitive variables
00047  *   Output: eigenvalues[4], right_eigenmatrix[4,4], left_eigenmatrix[4,4];
00048  */
00049 
00050 #if defined(BAROTROPIC) && defined(HYDRO)
00051 void esys_prim_iso_hyd(const Real d, const Real v1, 
00052   Real eigenvalues[],
00053   Real right_eigenmatrix[][4], Real left_eigenmatrix[][4])
00054 {
00055 
00056 /* Compute eigenvalues (eq. A6) */
00057 
00058   eigenvalues[0] = v1 - Iso_csound;
00059   eigenvalues[1] = v1;
00060   eigenvalues[2] = v1;
00061   eigenvalues[3] = v1 + Iso_csound;
00062 
00063 /* Right-eigenvectors, stored as COLUMNS (eq. A3) */
00064 
00065   right_eigenmatrix[0][0] = 1.0;
00066   right_eigenmatrix[1][0] = -Iso_csound/d;
00067 /*right_eigenmatrix[2][0] = 0.0; */
00068 /*right_eigenmatrix[3][0] = 0.0; */
00069 
00070 /*right_eigenmatrix[0][1] = 0.0; */
00071 /*right_eigenmatrix[1][1] = 0.0; */
00072   right_eigenmatrix[2][1] = 1.0;
00073 /*right_eigenmatrix[3][1] = 0.0; */
00074 
00075 /*right_eigenmatrix[0][2] = 0.0; */
00076 /*right_eigenmatrix[1][2] = 0.0; */
00077 /*right_eigenmatrix[2][2] = 0.0; */
00078   right_eigenmatrix[3][2] = 1.0;
00079 
00080   right_eigenmatrix[0][3] = 1.0;
00081   right_eigenmatrix[1][3] = Iso_csound/d;
00082 /*right_eigenmatrix[2][3] = 0.0; */
00083 /*right_eigenmatrix[3][3] = 0.0; */
00084 
00085 /* Left-eigenvectors, stored as ROWS (eq. A7) */
00086 
00087   left_eigenmatrix[0][0] = 0.5;
00088   left_eigenmatrix[0][1] = -0.5*d/Iso_csound;
00089 /*left_eigenmatrix[0][2] = 0.0; */
00090 /*left_eigenmatrix[0][3] = 0.0; */
00091 
00092 /*left_eigenmatrix[1][0] = 0.0; */
00093 /*left_eigenmatrix[1][1] = 0.0; */
00094   left_eigenmatrix[1][2] = 1.0;
00095 /*left_eigenmatrix[1][3] = 0.0; */
00096 
00097 /*left_eigenmatrix[2][0] = 0.0; */
00098 /*left_eigenmatrix[2][1] = 0.0; */
00099 /*left_eigenmatrix[2][2] = 0.0; */
00100   left_eigenmatrix[2][3] = 1.0;
00101 
00102   left_eigenmatrix[3][0] = 0.5;
00103   left_eigenmatrix[3][1] = 0.5*d/Iso_csound;
00104 /*left_eigenmatrix[3][2] = 0.0; */
00105 /*left_eigenmatrix[3][3] = 0.0; */
00106 
00107 }
00108 #endif
00109 
00110 /*----------------------------------------------------------------------------*/
00111 /*! \fn void esys_prim_adb_hyd(const Real d, const Real v1, const Real rho_a2,
00112  *  Real eigenvalues[],
00113  *  Real right_eigenmatrix[][5], Real left_eigenmatrix[][5])
00114  *  \brief ADIABATIC HYDRO
00115  *   Input: d, v1, p = primitive variables
00116  *   Output: eigenvalues[5], right_eigenmatrix[5,5], left_eigenmatrix[5,5];
00117  */
00118 
00119 #if !defined(BAROTROPIC) && defined(HYDRO)
00120 void esys_prim_adb_hyd(const Real d, const Real v1, const Real rho_a2,
00121   Real eigenvalues[],
00122   Real right_eigenmatrix[][5], Real left_eigenmatrix[][5])
00123 {
00124   Real asq,a;
00125   asq = rho_a2/d;
00126   a = sqrt(asq);
00127 
00128 /* Compute eigenvalues (eq. A2) */
00129 
00130   eigenvalues[0] = v1 - a;
00131   eigenvalues[1] = v1;
00132   eigenvalues[2] = v1;
00133   eigenvalues[3] = v1;
00134   eigenvalues[4] = v1 + a;
00135 
00136 /* Right-eigenvectors, stored as COLUMNS (eq. A3) */
00137 
00138   right_eigenmatrix[0][0] = 1.0;
00139   right_eigenmatrix[1][0] = -a/d;
00140 /*right_eigenmatrix[2][0] = 0.0; */
00141 /*right_eigenmatrix[3][0] = 0.0; */
00142   right_eigenmatrix[4][0] = asq;
00143 
00144   right_eigenmatrix[0][1] = 1.0;
00145 /*right_eigenmatrix[1][1] = 0.0; */
00146 /*right_eigenmatrix[2][1] = 0.0; */
00147 /*right_eigenmatrix[3][1] = 0.0; */
00148 /*right_eigenmatrix[4][1] = 0.0; */
00149 
00150 /*right_eigenmatrix[0][2] = 0.0; */
00151 /*right_eigenmatrix[1][2] = 0.0; */
00152   right_eigenmatrix[2][2] = 1.0;
00153 /*right_eigenmatrix[3][2] = 0.0; */
00154 /*right_eigenmatrix[4][2] = 0.0; */
00155 
00156 /*right_eigenmatrix[0][3] = 0.0; */
00157 /*right_eigenmatrix[1][3] = 0.0; */
00158 /*right_eigenmatrix[2][3] = 0.0; */
00159   right_eigenmatrix[3][3] = 1.0;
00160 /*right_eigenmatrix[4][3] = 0.0; */
00161 
00162   right_eigenmatrix[0][4] = 1.0;
00163   right_eigenmatrix[1][4] = -right_eigenmatrix[1][0];
00164 /*right_eigenmatrix[2][4] = 0.0; */
00165 /*right_eigenmatrix[3][4] = 0.0; */
00166   right_eigenmatrix[4][4] = asq;
00167 
00168 /* Left-eigenvectors, stored as ROWS (eq. A4) */
00169 
00170 /*left_eigenmatrix[0][0] = 0.0; */
00171   left_eigenmatrix[0][1] = -0.5*d/a;
00172 /*left_eigenmatrix[0][2] = 0.0; */
00173 /*left_eigenmatrix[0][3] = 0.0; */
00174   left_eigenmatrix[0][4] = 0.5/asq;
00175 
00176   left_eigenmatrix[1][0] = 1.0;
00177 /*left_eigenmatrix[1][1] = 0.0; */
00178 /*left_eigenmatrix[1][2] = 0.0; */
00179 /*left_eigenmatrix[1][3] = 0.0; */
00180   left_eigenmatrix[1][4] = -1.0/asq;
00181 
00182 /*left_eigenmatrix[2][0] = 0.0; */
00183 /*left_eigenmatrix[2][1] = 0.0; */
00184   left_eigenmatrix[2][2] = 1.0;
00185 /*left_eigenmatrix[2][3] = 0.0; */
00186 /*left_eigenmatrix[2][4] = 0.0; */
00187 
00188 /*left_eigenmatrix[3][0] = 0.0; */
00189 /*left_eigenmatrix[3][1] = 0.0; */
00190 /*left_eigenmatrix[3][2] = 0.0; */
00191   left_eigenmatrix[3][3] = 1.0;
00192 /*left_eigenmatrix[3][4] = 0.0; */
00193 
00194 /*left_eigenmatrix[4][0] = 0.0; */
00195   left_eigenmatrix[4][1] = -left_eigenmatrix[0][1];
00196 /*left_eigenmatrix[4][2] = 0.0; */
00197 /*left_eigenmatrix[4][3] = 0.0; */
00198   left_eigenmatrix[4][4] = left_eigenmatrix[0][4];
00199 }
00200 #endif
00201 
00202 /*----------------------------------------------------------------------------*/
00203 /*! \fn void esys_prim_iso_mhd(const Real d, const Real v1, const Real b1, 
00204  *  const Real b2, const Real b3, Real eigenvalues[],
00205  *  Real right_eigenmatrix[][6], Real left_eigenmatrix[][6])
00206  *  \brief ISOTHERMAL MHD
00207  *   Input: d, v1, b1, b2, b3 = density, velocities, and B field
00208  *   Output: eigenvalues[6], right_eigenmatrix[6,6], left_eigenmatrix[6,6];
00209  */
00210 
00211 #if defined(BAROTROPIC) && defined(MHD)
00212 void esys_prim_iso_mhd(const Real d, const Real v1, const Real b1,
00213   const Real b2, const Real b3, Real eigenvalues[],
00214   Real right_eigenmatrix[][6], Real left_eigenmatrix[][6])
00215 {
00216   Real btsq,vaxsq,cfsq,cf,cssq,cs,bt,bet2,bet3,alpha_f,alpha_s;
00217   Real sqrtd,s,qf,qs,af,as,vax,norm,af_prime,as_prime;
00218   Real ct2,tsum,tdif,cf2_cs2;
00219   Real di = 1.0/d;
00220   btsq  = b2*b2 + b3*b3;
00221   vaxsq = b1*b1*di;
00222 
00223 /* Compute fast- and slow-magnetosonic speeds (eq. A10) */
00224 
00225   ct2 = btsq*di;
00226   tsum = vaxsq + ct2 + (Iso_csound2);
00227   tdif = vaxsq + ct2 - (Iso_csound2);
00228   cf2_cs2 = sqrt((double)(tdif*tdif + 4.0*(Iso_csound2)*ct2));
00229 
00230   cfsq = 0.5*(tsum + cf2_cs2);
00231   cf = sqrt((double)cfsq);
00232 
00233   cssq = (Iso_csound2)*vaxsq/cfsq;
00234   cs = sqrt((double)cssq);
00235 
00236 /* Compute beta(s) (eq A17) */
00237 
00238   bt  = sqrt(btsq);
00239   if (bt == 0.0) {
00240     bet2 = 1.0;
00241     bet3 = 0.0;
00242   } else {
00243     bet2 = b2/bt;
00244     bet3 = b3/bt;
00245   }
00246 
00247 /* Compute alpha(s) (eq A16) */
00248 
00249   if ((cfsq-cssq) == 0.0) {
00250     alpha_f = 1.0;
00251     alpha_s = 0.0;
00252   } else if ( (Iso_csound2 - cssq) <= 0.0) {
00253     alpha_f = 0.0;
00254     alpha_s = 1.0;
00255   } else if ( (cfsq - Iso_csound2) <= 0.0) {
00256     alpha_f = 1.0;
00257     alpha_s = 0.0;
00258   } else {
00259     alpha_f = sqrt((Iso_csound2 - cssq)/(cfsq - cssq));
00260     alpha_s = sqrt((cfsq - Iso_csound2)/(cfsq - cssq));
00261   }
00262 
00263 /* Compute Q(s) and A(s) (eq. A14-15), etc. */
00264 
00265   sqrtd = sqrt(d);
00266   s = SIGN(b1);
00267   qf = cf*alpha_f*s;
00268   qs = cs*alpha_s*s;
00269   af = Iso_csound*alpha_f*sqrtd;
00270   as = Iso_csound*alpha_s*sqrtd;
00271 
00272 /* Compute eigenvalues (eq. A21) */
00273 
00274   vax = sqrt(vaxsq);
00275   eigenvalues[0] = v1 - cf;
00276   eigenvalues[1] = v1 - vax;
00277   eigenvalues[2] = v1 - cs;
00278   eigenvalues[3] = v1 + cs;
00279   eigenvalues[4] = v1 + vax;
00280   eigenvalues[5] = v1 + cf;
00281 
00282 /* Right-eigenvectors, stored as COLUMNS (eq A12) */
00283 
00284   right_eigenmatrix[0][0] = d*alpha_f;
00285   right_eigenmatrix[1][0] = -cf*alpha_f;
00286   right_eigenmatrix[2][0] = qs*bet2;
00287   right_eigenmatrix[3][0] = qs*bet3;
00288   right_eigenmatrix[4][0] = as*bet2;
00289   right_eigenmatrix[5][0] = as*bet3;
00290 
00291 /*right_eigenmatrix[0][1] = 0.0; */
00292 /*right_eigenmatrix[1][1] = 0.0; */
00293   right_eigenmatrix[2][1] = -bet3;
00294   right_eigenmatrix[3][1] = bet2;
00295   right_eigenmatrix[4][1] = -bet3*s*sqrtd;
00296   right_eigenmatrix[5][1] = bet2*s*sqrtd;
00297 
00298   right_eigenmatrix[0][2] = d*alpha_s;
00299   right_eigenmatrix[1][2] = -cs*alpha_s;
00300   right_eigenmatrix[2][2] = -qf*bet2;
00301   right_eigenmatrix[3][2] = -qf*bet3;
00302   right_eigenmatrix[4][2] = -af*bet2;
00303   right_eigenmatrix[5][2] = -af*bet3;
00304 
00305   right_eigenmatrix[0][3] = right_eigenmatrix[0][2];
00306   right_eigenmatrix[1][3] = -right_eigenmatrix[1][2];
00307   right_eigenmatrix[2][3] = -right_eigenmatrix[2][2];
00308   right_eigenmatrix[3][3] = -right_eigenmatrix[3][2];
00309   right_eigenmatrix[4][3] = right_eigenmatrix[4][2];
00310   right_eigenmatrix[5][3] = right_eigenmatrix[5][2];
00311 
00312 /*right_eigenmatrix[0][4] = 0.0; */
00313 /*right_eigenmatrix[1][4] = 0.0; */
00314   right_eigenmatrix[2][4] = bet3;
00315   right_eigenmatrix[3][4] = -bet2;
00316   right_eigenmatrix[4][4] = right_eigenmatrix[4][1];
00317   right_eigenmatrix[5][4] = right_eigenmatrix[5][1];
00318 
00319   right_eigenmatrix[0][5] = right_eigenmatrix[0][0];
00320   right_eigenmatrix[1][5] = -right_eigenmatrix[1][0];
00321   right_eigenmatrix[2][5] = -right_eigenmatrix[2][0];
00322   right_eigenmatrix[3][5] = -right_eigenmatrix[3][0];
00323   right_eigenmatrix[4][5] = right_eigenmatrix[4][0];
00324   right_eigenmatrix[5][5] = right_eigenmatrix[5][0];
00325 
00326 /* Left-eigenvectors, stored as ROWS (eq A22) */
00327 
00328   norm = 0.5/Iso_csound2;
00329   qf = norm*qf;
00330   qs = norm*qs;
00331   af_prime = norm*af*di;
00332   as_prime = norm*as*di;
00333 
00334   left_eigenmatrix[0][0] = norm*alpha_f*Iso_csound2*di;
00335   left_eigenmatrix[0][1] = -norm*cf*alpha_f;
00336   left_eigenmatrix[0][2] = qs*bet2;
00337   left_eigenmatrix[0][3] = qs*bet3;
00338   left_eigenmatrix[0][4] = as_prime*bet2;
00339   left_eigenmatrix[0][5] = as_prime*bet3;
00340 
00341 /*left_eigenmatrix[1][0] = 0.0; */
00342 /*left_eigenmatrix[1][1] = 0.0; */
00343   left_eigenmatrix[1][2] = -0.5*bet3;
00344   left_eigenmatrix[1][3] = 0.5*bet2;
00345   left_eigenmatrix[1][4] = -0.5*bet3*s/sqrtd;
00346   left_eigenmatrix[1][5] = 0.5*bet2*s/sqrtd;
00347 
00348   left_eigenmatrix[2][0] = norm*alpha_s*Iso_csound2*di;
00349   left_eigenmatrix[2][1] = -norm*cs*alpha_s;
00350   left_eigenmatrix[2][2] = -qf*bet2;
00351   left_eigenmatrix[2][3] = -qf*bet3;
00352   left_eigenmatrix[2][4] = -af_prime*bet2;
00353   left_eigenmatrix[2][5] = -af_prime*bet3;
00354 
00355   left_eigenmatrix[3][0] = left_eigenmatrix[2][0];
00356   left_eigenmatrix[3][1] = -left_eigenmatrix[2][1];
00357   left_eigenmatrix[3][2] = -left_eigenmatrix[2][2];
00358   left_eigenmatrix[3][3] = -left_eigenmatrix[2][3];
00359   left_eigenmatrix[3][4] = left_eigenmatrix[2][4];
00360   left_eigenmatrix[3][5] = left_eigenmatrix[2][5];
00361 
00362 /*left_eigenmatrix[4][0] = 0.0; */
00363 /*left_eigenmatrix[4][1] = 0.0; */
00364   left_eigenmatrix[4][2] = -left_eigenmatrix[1][2];
00365   left_eigenmatrix[4][3] = -left_eigenmatrix[1][3];
00366   left_eigenmatrix[4][4] = left_eigenmatrix[1][4];
00367   left_eigenmatrix[4][5] = left_eigenmatrix[1][5];
00368 
00369   left_eigenmatrix[5][0] = left_eigenmatrix[0][0];
00370   left_eigenmatrix[5][1] = -left_eigenmatrix[0][1];
00371   left_eigenmatrix[5][2] = -left_eigenmatrix[0][2];
00372   left_eigenmatrix[5][3] = -left_eigenmatrix[0][3];
00373   left_eigenmatrix[5][4] = left_eigenmatrix[0][4];
00374   left_eigenmatrix[5][5] = left_eigenmatrix[0][5];
00375 }
00376 #endif
00377 
00378 /*----------------------------------------------------------------------------*/
00379 /*! \fn void esys_prim_adb_mhd(const Real d, const Real v1, const Real rho_a2,
00380  *  const Real b1, const Real b2, const Real b3, Real eigenvalues[],
00381  *  Real right_eigenmatrix[][7], Real left_eigenmatrix[][7])
00382  *  \brief ADIABATIC MHD
00383  *   Input: d, v1, p, b1, b2, b3 = density, velocities, pressure, and B field
00384  *   Output: eigenvalues[7], right_eigenmatrix[7,7], left_eigenmatrix[7,7];
00385  */
00386 
00387 #if !defined(BAROTROPIC) && defined(MHD)
00388 void esys_prim_adb_mhd(const Real d, const Real v1, const Real rho_a2,
00389   const Real b1, const Real b2, const Real b3, Real eigenvalues[],
00390   Real right_eigenmatrix[][7], Real left_eigenmatrix[][7])
00391 {
00392   Real di,cfsq,cf,cssq,cs,bt,bet2,bet3,alpha_f,alpha_s;
00393   Real sqrtd,s,a,qf,qs,af,as,vax,na,af_prime,as_prime;
00394   Real tsum,tdif,cf2_cs2,ct2;
00395   Real btsq,vaxsq,asq;
00396   di = 1.0/d;
00397   btsq  = b2*b2 + b3*b3;
00398   vaxsq = b1*b1*di;
00399   asq   = rho_a2*di;
00400 
00401 /* Compute fast- and slow-magnetosonic speeds (eq. A10) */
00402 
00403   ct2 = btsq*di;
00404   tsum = vaxsq + ct2 + asq;
00405   tdif = vaxsq + ct2 - asq;
00406   cf2_cs2 = sqrt((double)(tdif*tdif + 4.0*asq*ct2));
00407 
00408   cfsq = 0.5*(tsum + cf2_cs2);
00409   cf = sqrt((double)cfsq);
00410 
00411   cssq = asq*vaxsq/cfsq;
00412   cs = sqrt((double)cssq);
00413 
00414 /* Compute beta(s) (eq A17) */
00415 
00416   bt  = sqrt(btsq);
00417   if (bt == 0.0) {
00418     bet2 = 1.0;
00419     bet3 = 0.0;
00420   } else {
00421     bet2 = b2/bt;
00422     bet3 = b3/bt;
00423   }
00424 
00425 /* Compute alpha(s) (eq A16) */
00426 
00427   if (cf2_cs2 == 0.0) {
00428     alpha_f = 1.0;
00429     alpha_s = 0.0;
00430   } else if ( (asq - cssq) <= 0.0) {
00431     alpha_f = 0.0;
00432     alpha_s = 1.0;
00433   } else if ( (cfsq - asq) <= 0.0) {
00434     alpha_f = 1.0;
00435     alpha_s = 0.0;
00436   } else {
00437     alpha_f = sqrt((asq - cssq)/cf2_cs2);
00438     alpha_s = sqrt((cfsq - asq)/cf2_cs2);
00439   }
00440 
00441 /* Compute Q(s) and A(s) (eq. A14-15), etc. */
00442 
00443   sqrtd = sqrt(d);
00444   s = SIGN(b1);
00445   a = sqrt(asq);
00446   qf = cf*alpha_f*s;
00447   qs = cs*alpha_s*s;
00448   af = a*alpha_f*sqrtd;
00449   as = a*alpha_s*sqrtd;
00450 
00451 /* Compute eigenvalues (eq. A9) */
00452 
00453   vax = sqrt(vaxsq);
00454   eigenvalues[0] = v1 - cf;
00455   eigenvalues[1] = v1 - vax;
00456   eigenvalues[2] = v1 - cs;
00457   eigenvalues[3] = v1;
00458   eigenvalues[4] = v1 + cs;
00459   eigenvalues[5] = v1 + vax;
00460   eigenvalues[6] = v1 + cf;
00461 
00462 /* Right-eigenvectors, stored as COLUMNS (eq. A12) */
00463 /* Note statements are grouped in ROWS for optimization, even though rem[*][n]
00464  * is the nth right eigenvector */
00465 
00466 
00467   right_eigenmatrix[0][0] = d*alpha_f;
00468 /*right_eigenmatrix[0][1] = 0.0; */
00469   right_eigenmatrix[0][2] = d*alpha_s;
00470   right_eigenmatrix[0][3] = 1.0;
00471   right_eigenmatrix[0][4] = right_eigenmatrix[0][2];
00472 /*right_eigenmatrix[0][5] = 0.0; */
00473   right_eigenmatrix[0][6] = right_eigenmatrix[0][0];
00474 
00475   right_eigenmatrix[1][0] = -cf*alpha_f;
00476 /*right_eigenmatrix[1][1] = 0.0; */
00477   right_eigenmatrix[1][2] = -cs*alpha_s;
00478 /*right_eigenmatrix[1][3] = 0.0; */
00479   right_eigenmatrix[1][4] = -right_eigenmatrix[1][2];
00480 /*right_eigenmatrix[1][5] = 0.0; */
00481   right_eigenmatrix[1][6] = -right_eigenmatrix[1][0];
00482 
00483   right_eigenmatrix[2][0] = qs*bet2;
00484   right_eigenmatrix[2][1] = -bet3;
00485   right_eigenmatrix[2][2] = -qf*bet2;
00486 /*right_eigenmatrix[2][3] = 0.0; */
00487   right_eigenmatrix[2][4] = -right_eigenmatrix[2][2];
00488   right_eigenmatrix[2][5] = bet3;
00489   right_eigenmatrix[2][6] = -right_eigenmatrix[2][0];
00490 
00491   right_eigenmatrix[3][0] = qs*bet3;
00492   right_eigenmatrix[3][1] = bet2;
00493   right_eigenmatrix[3][2] = -qf*bet3;
00494 /*right_eigenmatrix[3][3] = 0.0; */
00495   right_eigenmatrix[3][4] = -right_eigenmatrix[3][2];
00496   right_eigenmatrix[3][5] = -bet2;
00497   right_eigenmatrix[3][6] = -right_eigenmatrix[3][0];
00498 
00499   right_eigenmatrix[4][0] = d*asq*alpha_f;
00500 /*right_eigenmatrix[4][1] = 0.0; */
00501   right_eigenmatrix[4][2] = d*asq*alpha_s;
00502 /*right_eigenmatrix[4][3] = 0.0; */
00503   right_eigenmatrix[4][4] = right_eigenmatrix[4][2];
00504 /*right_eigenmatrix[4][5] = 0.0; */
00505   right_eigenmatrix[4][6] = right_eigenmatrix[4][0];
00506 
00507   right_eigenmatrix[5][0] = as*bet2;
00508   right_eigenmatrix[5][1] = -bet3*s*sqrtd;
00509   right_eigenmatrix[5][2] = -af*bet2;
00510 /*right_eigenmatrix[5][3] = 0.0; */
00511   right_eigenmatrix[5][4] = right_eigenmatrix[5][2];
00512   right_eigenmatrix[5][5] = right_eigenmatrix[5][1];
00513   right_eigenmatrix[5][6] = right_eigenmatrix[5][0];
00514 
00515   right_eigenmatrix[6][0] = as*bet3;
00516   right_eigenmatrix[6][1] = bet2*s*sqrtd;
00517   right_eigenmatrix[6][2] = -af*bet3;
00518 /*right_eigenmatrix[6][3] = 0.0; */
00519   right_eigenmatrix[6][4] = right_eigenmatrix[6][2];
00520   right_eigenmatrix[6][5] = right_eigenmatrix[6][1];
00521   right_eigenmatrix[6][6] = right_eigenmatrix[6][0];
00522 
00523 /* Left-eigenvectors, stored as ROWS (eq. A18) */
00524 
00525   na = 0.5/asq;
00526   qf = na*qf;
00527   qs = na*qs;
00528   af_prime = na*af*di;
00529   as_prime = na*as*di;
00530 
00531 /*left_eigenmatrix[0][0] = 0.0; */
00532   left_eigenmatrix[0][1] = -na*cf*alpha_f;
00533   left_eigenmatrix[0][2] = qs*bet2;
00534   left_eigenmatrix[0][3] = qs*bet3;
00535   left_eigenmatrix[0][4] = na*alpha_f*di;
00536   left_eigenmatrix[0][5] = as_prime*bet2;
00537   left_eigenmatrix[0][6] = as_prime*bet3;
00538 
00539 /*left_eigenmatrix[1][0] = 0.0; */
00540 /*left_eigenmatrix[1][1] = 0.0; */
00541   left_eigenmatrix[1][2] = -0.5*bet3;
00542   left_eigenmatrix[1][3] = 0.5*bet2;
00543 /*left_eigenmatrix[1][4] = 0.0; */
00544   left_eigenmatrix[1][5] = -0.5*bet3*s/sqrtd;
00545   left_eigenmatrix[1][6] = 0.5*bet2*s/sqrtd;
00546 
00547 /*left_eigenmatrix[2][0] = 0.0; */
00548   left_eigenmatrix[2][1] = -na*cs*alpha_s;
00549   left_eigenmatrix[2][2] = -qf*bet2;
00550   left_eigenmatrix[2][3] = -qf*bet3;
00551   left_eigenmatrix[2][4] = na*alpha_s*di;
00552   left_eigenmatrix[2][5] = -af_prime*bet2;
00553   left_eigenmatrix[2][6] = -af_prime*bet3;
00554 
00555   left_eigenmatrix[3][0] = 1.0;
00556 /*left_eigenmatrix[3][1] = 0.0; */
00557 /*left_eigenmatrix[3][2] = 0.0; */
00558 /*left_eigenmatrix[3][3] = 0.0; */
00559   left_eigenmatrix[3][4] = -1.0/asq;
00560 /*left_eigenmatrix[3][5] = 0.0; */
00561 /*left_eigenmatrix[3][6] = 0.0; */
00562 
00563 /*left_eigenmatrix[4][0] = 0.0; */
00564   left_eigenmatrix[4][1] = -left_eigenmatrix[2][1];
00565   left_eigenmatrix[4][2] = -left_eigenmatrix[2][2];
00566   left_eigenmatrix[4][3] = -left_eigenmatrix[2][3];
00567   left_eigenmatrix[4][4] = left_eigenmatrix[2][4];
00568   left_eigenmatrix[4][5] = left_eigenmatrix[2][5];
00569   left_eigenmatrix[4][6] = left_eigenmatrix[2][6];
00570 
00571 /*left_eigenmatrix[5][0] = 0.0; */
00572 /*left_eigenmatrix[5][1] = 0.0; */
00573   left_eigenmatrix[5][2] = -left_eigenmatrix[1][2];
00574   left_eigenmatrix[5][3] = -left_eigenmatrix[1][3];
00575 /*left_eigenmatrix[5][4] = 0.0; */
00576   left_eigenmatrix[5][5] = left_eigenmatrix[1][5];
00577   left_eigenmatrix[5][6] = left_eigenmatrix[1][6];
00578 
00579 /*left_eigenmatrix[6][0] = 0.0; */
00580   left_eigenmatrix[6][1] = -left_eigenmatrix[0][1];
00581   left_eigenmatrix[6][2] = -left_eigenmatrix[0][2];
00582   left_eigenmatrix[6][3] = -left_eigenmatrix[0][3];
00583   left_eigenmatrix[6][4] = left_eigenmatrix[0][4];
00584   left_eigenmatrix[6][5] = left_eigenmatrix[0][5];
00585   left_eigenmatrix[6][6] = left_eigenmatrix[0][6];
00586 }
00587 #endif
reconstruction/esystem_prim.c
