00001 #include "../copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
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
00053
00054
00055
00056
00057
00058
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
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
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
00083
00084 right_eigenmatrix[2][1] = 1.0;
00085
00086
00087
00088
00089
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
00098
00099 left_eigenmatrix[0][0] = 0.5*(1.0 + v1/Iso_csound);
00100 left_eigenmatrix[0][1] = -0.5/Iso_csound;
00101
00102
00103
00104 left_eigenmatrix[1][0] = -v2;
00105
00106 left_eigenmatrix[1][2] = 1.0;
00107
00108
00109 left_eigenmatrix[2][0] = -v3;
00110
00111
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
00117
00118 }
00119 #endif
00120
00121
00122
00123
00124
00125
00126
00127
00128
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
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
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
00159
00160 right_eigenmatrix[2][1] = 1.0;
00161
00162 right_eigenmatrix[4][1] = v2;
00163
00164
00165
00166
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
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
00193 left_eigenmatrix[1][2] = 1.0;
00194
00195
00196
00197 left_eigenmatrix[2][0] = -v3;
00198
00199
00200 left_eigenmatrix[2][3] = 1.0;
00201
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
00220
00221
00222
00223
00224
00225
00226
00227
00228
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
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
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
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
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
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
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
00324
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
00345
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
00359
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
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
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
00421
00422
00423
00424
00425
00426
00427
00428
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
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
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
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
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
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
00523
00524
00525
00526 right_eigenmatrix[0][0] = alpha_f;
00527
00528 right_eigenmatrix[0][2] = alpha_s;
00529 right_eigenmatrix[0][3] = 1.0;
00530 right_eigenmatrix[0][4] = alpha_s;
00531
00532 right_eigenmatrix[0][6] = alpha_f;
00533
00534 right_eigenmatrix[1][0] = alpha_f*eigenvalues[0];
00535
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
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
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
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
00591
00592
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
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
00622 left_eigenmatrix[1][2] = -0.5*bet3;
00623 left_eigenmatrix[1][3] = 0.5*bet2;
00624
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
00654 left_eigenmatrix[5][2] = -left_eigenmatrix[1][2];
00655 left_eigenmatrix[5][3] = -left_eigenmatrix[1][3];
00656
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