/*he1lev evaluate ionization balance level populations for helium singlets */ #include #include #include "cddefines.h" #include "physconst.h" #include "taulines.h" #include "nhe3lvl.h" #include "nhe1lvl.h" #include "showme.h" #include "rfield.h" #include "herec.h" #include "he1blt.h" #include "he1crt.h" #include "secondaries.h" #include "he1ind.h" #include "trace.h" #include "tepowers.h" #include "h.h" #include "he1tau.h" #include "phycon.h" #include "phe1lv.h" #include "he1ex.h" #include "he1cbr.h" #include "he1bn.h" #include "he1rsp.h" #include "he3n.h" #include "che2sn.h" #include "he1rb.h" #include "the1lo.h" #include "zonecnt.h" #include "nhe1.h" #include "hevmolec.h" #include "he1net.h" #include "destcrt.h" #include "he1stat.h" #include "he1rate.h" #include "he3as.h" #include "he3tau.h" #include "he3pht.h" #include "he1nionryd.h" #include "typmatrx.h" #include "linpack.h" #include "matin1.h" #include "veclib.h" #include "negcon.h" #include "he1.h" void he1lev(double *he2ov1, float *sm2ov1) { long int i, ipiv[NHE1LVL + 1], j, job, nerror; double chk2ov1, constant, edp, gdist, he1exx[NHE1LVL + 1], hex12, rdest, rfac, totcap[10]; static double pop3n1 = 0., rtrip2p, rtrip2s, rtripg; double amat[NHE1LVL + 1][NHE1LVL + 1], bvec[NHE1LVL + 1], c, r, rcond, work[NHE1LVL + 1], z[NHE1LVL + 2][NHE1LVL + 2]; # ifdef DEBUG_FUN fputs( "<+>he1lev()\n", debug_fp ); # endif /* nhe1lvl is 9 */ /* following used for linpac matrix inversion */ /* real a2sg * this is in he3as.com * parameter (a2sg=1.13e-4) */ /* real fudge */ /* 2S LEVEL IS LEVEL TEN; Z(2,I) IS 2P, Z(10,I) = 2S, HN(2) IS SUM * level 7 is actually 7 - 10; 8 is 11 - 20; 9 is 21 - 100 * free statement labels >=12 * * last term is frac of trips that get photoionized * He1RecB is case b singlet rec coef * he1rtt is effect total rec coef */ herecCom.reci = he1rb.he1rtt + herecCom.rectri*he3pht.he3relax; /* do simple estimate of ionization balance * following assumes all trip recomb go to singlets */ *sm2ov1 = (float)((phe1lv.he1gam[0] + Secondaries.csupra*0.5 + phe1lv.he1clc[0]* phycon.eden)/(phycon.eden*herecCom.reci + hevmolec.hevmol[IPCO-1]* 1.10e-9)); /* >>chng 96 dec 17, for HeI oscillations when pumped Lyman lines *resulted in photoionization, simple estimate no good, so use *previous solution */ if( ZoneCnt.nzone > 1 ) { /* use previous matrix solution if it exists */ chk2ov1 = he1rate.clhe2ov1; } else { /* use simple estimate for beginning of calculation */ chk2ov1 = *sm2ov1; } he1rate.smhe2ov1 = *sm2ov1; if( trace.lgHe1Bug && trace.lgTrace ) { fprintf( ioQQQ, " HE1LEV HE1GAM%10.2e csupra%10.2e HE1CLC(1)%10.2e rec*re%10.2e smn2ov1%10.2e\n", phe1lv.he1gam[0], Secondaries.csupra, phe1lv.he1clc[0], phycon.eden* herecCom.reci, *sm2ov1 ); } if( trace.lgHeBug && trace.lgTrace ) { fprintf( ioQQQ, " HE1LEV HE1GAM tot ioniz%10.2e total reco%10.2e SM2OV1=%10.2e\n", phe1lv.he1gam[0] + Secondaries.csupra*0.5 + phe1lv.he1clc[0]* phycon.eden, phycon.eden*herecCom.reci + hevmolec.hevmol[IPCO-1]* 1.10e-9, *sm2ov1 ); } DestCrt.create[0][1] = phycon.eden*herecCom.reci + hevmolec.hevmol[IPCO-1]* 1.10e-9; DestCrt.destroy[0][1] = phe1lv.he1gam[0] + Secondaries.csupra*0.5 + phe1lv.he1clc[0]*phycon.eden; /* check how low temperature is; departure coef tend to infinity * as TE goes to zero; matrix formalism won't work on a 32-bit machine * if EXP(H NU/KT) > 1E38 for 13.6EV. * also use this if very low ionization * >>chng 96 dec 17, from sm2ov1 to chk2ov1 */ if( ((phycon.te <= the1loCom.the1lo || he1crt.he1lte[0] == 0.) || he1crt.he1lte[1] == 0.) || chk2ov1 < 1e-8 ) { /* remember that things are too cool, for later analysis */ the1loCom.nhe1lo += 1; /* last part above added to prevent code from using following * when helium is highly COLLISIONALLY ionized */ he1exCOM.he1ex = 0.; he1bnCOM.he1bn[0] = 1.; for( i=1; i < (NHE1LVL + 1); i++ ) { he1bnCOM.he1bn[i] = 1.; phe1lv.he1n[i] = 0.; } /* He1RecB is evaluated with rec coef, is case b sum * recexc = He1RecB + * 1 (he3n(1)*(c2ss2p+c2ss2s) + he3n(2)*(c2ps2s+c2ps2p) ) / eden * reci = he1rec(1,1)*he1rec(1,2) + recexc + * 1 (he3n(1)*(c2sg+a2sg) + he3n(2)*c2pg ) / eden * * do simple two level for helium, then fill in HN(i) for rest */ rdest = phe1lv.he1gam[0] + Secondaries.csupra*0.5 + phe1lv.he1clc[0]* phycon.eden; if( rdest > 0. ) { /* charge transfer with CO in here */ phe1lv.he1n[0] = (float)((phycon.eden*herecCom.reci + hevmolec.hevmol[IPCO-1]* 1.10e-9)/rdest); } else { phe1lv.he1n[0] = FLT_MAX; } /* effective case B rec coef for low den limit Martin 1988 */ phe1lv.he1n[3] = (float)(1.44e-11/(TePowers.te70/TePowers.te03/TePowers.te03* TePowers.te01)*phycon.eden/he1netCOM.he1net[1][3]); constant = phe1lv.he1n[3]*he1netCOM.he1net[1][3]; phe1lv.he1n[2] = (float)(constant*14.5/TePowers.te10/TePowers.te03/TePowers.te03* TePowers.te01/he1netCOM.he1net[1][2]); phe1lv.he1n[4] = (float)(constant*0.291*TePowers.te03*TePowers.te01/ he1netCOM.he1net[1][4]); phe1lv.he1n[5] = (float)(constant*0.138*TePowers.te03*TePowers.te01* TePowers.te01/he1netCOM.he1net[1][5]); /* He1RecB is evaluated with rec coef, is case b sum * this needs to have a term for exchange collisions from trips, * as in recexc which is commented out above. It is commented * out because of oscillations it caused. */ phe1lv.he12s = (float)(he1rb.He1RecB*phycon.eden/(8.23 + che2sn.cs2s2p + he1crt.he1c2s1s*phycon.eden)*0.25); phe1lv.he12p = (float)(he1rb.He1RecB*phycon.eden/(phe1lv.he1cll[0][1]* phycon.eden + phe1lv.he1dst[0][1] + he1netCOM.he1net[0][1])*0.75); phe1lv.he1n[1] = phe1lv.he12s + phe1lv.he12p; *he2ov1 = 1./phe1lv.he1n[0]; he1rate.clhe2ov1 = (float)*he2ov1; if( trace.lgHeBug && trace.lgTrace ) { fprintf( ioQQQ, " LOW TE, HE1N(1)=%12.3e rec=%12.3eHE1gam(1)=%12.3e\n", phe1lv.he1n[0], herecCom.reci, phe1lv.he1gam[0] ); } # ifdef DEBUG_FUN fputs( " <->he1lev()\n", debug_fp ); # endif return; } /*------------------------------------------------------------------- * * radiative and three body rec, fix 2s and 2p first * remember that HE1LTE(i) is only LTE pop when mult by Ne Nh*/ /* get rates from triplets */ if( he3nCom.he3n[1] > 0. ) { /* 666 error! when going over to one eden change per iteration, get rid * of only updates when flag set * use this if we have a valid triplet solution, hn3n(2)>0 * he3n is pop relative to heii */ rtrip2s = (he3nCom.he3n[0]*che2sn.c2ss2s + he3nCom.he3n[1]* che2sn.c2ps2s)/phycon.eden; } else { /* no triplet solution, so assume no recs come here * >>chng 96 feb 14, have been half of triplets here */ rtrip2s = 0.; } z[10][9] = (he1rsp.he1r2s*phe1lv.he1rec[1][1] + rtrip2s + he1ind.he1rin[1]* 0.25)/he1crt.he1lte[9] + phe1lv.he1clc[1]*phycon.eden; totcap[9] = z[10][9]*he1crt.he1lte[9]; herecCom.reci = (float)totcap[9]; if( he3nCom.he3n[1] > 0. ) { /* use this if we have a valid triplet solution */ rtrip2p = (he3nCom.he3n[0]*che2sn.c2ss2p + he3nCom.he3n[1]* che2sn.c2ps2p)/phycon.eden; } else { /* no triplet solution, so assume all recs come here * >>chng 96 feb 14, have been half of triplets here */ rtrip2p = 0.; } z[10][1] = (he1rsp.he1r2p*phe1lv.he1rec[1][1] + rtrip2p + he1ind.he1rin[1]* 0.75)/he1crt.he1lte[1] + phe1lv.he1clc[1]*phycon.eden; totcap[1] = z[10][1]*he1crt.he1lte[1]; herecCom.reci += (float)totcap[1]; if( ZoneCnt.nzone > 1 ) { /* this atom couples with the triplets - * damp out oscillations here */ pop3n1 = (pop3n1 + he3nCom.he3n[0])/2.; } else { pop3n1 = he3nCom.he3n[0]; } pop3n1 = he3nCom.he3n[0]; if( he3nCom.he3n[1] > 0. ) { rtripg = (pop3n1*(che2sn.c2sg + he3as.a2sg) + he3nCom.he3n[1]* che2sn.c2pg)/phycon.eden; } else { /* no triplet solution, so assume all recs come here * >>chng 96 feb 14, had been zero */ rtripg = herecCom.rectri*he3pht.he3relax; } z[10][0] = (phe1lv.he1rec[0][0]*phe1lv.he1rec[1][0] + rtripg + he1ind.he1rin[0])/he1crt.he1lte[0] + phe1lv.he1clc[0]*phycon.eden; totcap[0] = z[10][0]*he1crt.he1lte[0]; herecCom.reci += (float)totcap[0]; for( i=2; i < NHE1LVL; i++ ) { z[10][i] = (phe1lv.he1rec[0][i]*phe1lv.he1rec[1][i] + he1ind.he1rin[i])/ he1crt.he1lte[i] + phe1lv.he1clc[i]*phycon.eden; totcap[i] = z[10][i]*he1crt.he1lte[i]; herecCom.reci += (float)totcap[i]; } /* Janev et al show HeI 100 eV cross section is about 6 times * smaller than that for hydrogen */ hex12 = Secondaries.x12tot/6.; /*--------------------------------------------------------------------- * now set up matrix for dreaded b(i)'s*/ /* level 1 balance equation * HE1CLL(u,l), etc, are rates for downward collisions (cm^3 s^-1) * corr of hi*1.7e-4 accounts for col ion by HI; Drawin Zs Phys 225, 483. */ z[0][0] = Secondaries.csupra*0.5 + phe1lv.he1gam[0] + hex12 + he1cbr.he1coldn[0]* phycon.eden + phe1lv.he1clc[0]*(phycon.eden + h.hi*1.7e-4) + he1cbr.he1bul[0] + (che2sn.csg2s + che2sn.csg2p); z[1][0] = -(he1netCOM.he1net[0][1] + he1tauCOM.he1con[0][1] + phe1lv.he1cll[0][1]*phycon.eden + hex12/3.)*he1crt.he1lte[1]/ he1crt.he1lte[0]; z[2][0] = -(he1netCOM.he1net[0][2] + he1tauCOM.he1con[0][2] + phe1lv.he1cll[0][2]*phycon.eden)*he1crt.he1lte[2]/he1crt.he1lte[0]; z[3][0] = -(he1netCOM.he1net[0][3] + he1tauCOM.he1con[0][3] + phe1lv.he1cll[0][3]*phycon.eden)*he1crt.he1lte[3]/he1crt.he1lte[0]; z[4][0] = -(he1netCOM.he1net[0][4] + he1tauCOM.he1con[0][4] + phe1lv.he1cll[0][4]*phycon.eden)*he1crt.he1lte[4]/he1crt.he1lte[0]; z[5][0] = -(he1netCOM.he1net[0][5] + he1tauCOM.he1con[0][5] + phe1lv.he1cll[0][5]*phycon.eden)*he1crt.he1lte[5]/he1crt.he1lte[0]; z[6][0] = -(he1netCOM.he1net[0][6] + he1tauCOM.he1con[0][6] + phe1lv.he1cll[0][6]*phycon.eden)*he1crt.he1lte[6]/he1crt.he1lte[0]; z[7][0] = -(he1netCOM.he1net[0][7] + he1tauCOM.he1con[0][7] + phe1lv.he1cll[0][7]*phycon.eden)*he1crt.he1lte[7]/he1crt.he1lte[0]; z[8][0] = -(he1netCOM.he1net[0][8] + he1tauCOM.he1con[0][8] + phe1lv.he1cll[0][8]*phycon.eden)*he1crt.he1lte[8]/he1crt.he1lte[0]; z[9][0] = -(8.23 + hex12 + he1crt.he1c2s1s*phycon.eden)*he1crt.he1lte[9]/ he1crt.he1lte[0]; gdist = z[0][0]; /* LEVEL 2P BALANCE EQUATION - 2p is higher energy than 2s */ edp = phycon.eden*0.1 + h.hii*0.9; c = 0.75/he1crt.he1lte[1]; /* R is low density limit branching ratio for 2s-p (NOT l-mixed) * this is to simulate mixing of the n=3 level */ rfac = 0.32 - 0.07*phycon.eden/(phycon.eden + 1e7); r = (1. - rfac)/he1crt.he1lte[1]; z[0][1] = -(hex12*0.33 + phe1lv.he1cll[0][1]*phycon.eden + he1tauCOM.he1con[0][1]/ he1blt.he1bol[1][0]); z[1][1] = phe1lv.he1gam[1] + he1cbr.he1radn[1] + hex12/3. + he1cbr.he1coldn[1]* phycon.eden + phe1lv.he1clc[1]*phycon.eden + he1crt.he1c2p2s* edp + (he1cbr.he1bul[1] + he1tauCOM.he1con[0][1]) + (che2sn.cs2p2s + che2sn.cs2p2p); z[2][1] = -(he1netCOM.he1net[1][2]*r + (he1tauCOM.he1con[1][2] + phe1lv.he1cll[1][2]*phycon.eden)*c)*he1crt.he1lte[2]; z[3][1] = -(he1netCOM.he1net[1][3]*r + (he1tauCOM.he1con[1][3] + phe1lv.he1cll[1][3]*phycon.eden)*c)*he1crt.he1lte[3]; z[4][1] = -(he1netCOM.he1net[1][4]*r + (he1tauCOM.he1con[1][4] + phe1lv.he1cll[1][4]*phycon.eden)*c)*he1crt.he1lte[4]; z[5][1] = -(he1netCOM.he1net[1][5]*r + (he1tauCOM.he1con[1][5] + phe1lv.he1cll[1][5]*phycon.eden)*c)*he1crt.he1lte[5]; z[6][1] = -(he1netCOM.he1net[1][6]*r + (he1tauCOM.he1con[1][6] + phe1lv.he1cll[1][6]*phycon.eden)*c)*he1crt.he1lte[6]; z[7][1] = -(he1netCOM.he1net[1][7]*r + (he1tauCOM.he1con[1][7] + phe1lv.he1cll[1][7]*phycon.eden)*c)*he1crt.he1lte[7]; z[8][1] = -(he1netCOM.he1net[1][8]*r + (he1tauCOM.he1con[1][8] + phe1lv.he1cll[1][8]*phycon.eden)*c)*he1crt.he1lte[8]; z[9][1] = -he1crt.he1c2p2s*edp; /* LEVEL 2S BALANCE EQUATION - define 2s as lower energy than 2p */ c = 0.25/he1crt.he1lte[9]; r = rfac/he1crt.he1lte[9]; z[0][9] = -(hex12*0.33 + he1crt.he1c2s1s*phycon.eden); z[1][9] = -(he1crt.he1c2p2s*edp + t206.Pesc*t206.Aul)* he1crt.he1lte[1]/he1crt.he1lte[9]; z[2][9] = -(he1netCOM.he1net[1][2]*r + (he1tauCOM.he1con[1][2] + phe1lv.he1cll[1][2]*phycon.eden)*c)*he1crt.he1lte[2]; z[3][9] = -(he1netCOM.he1net[1][3]*r + (he1tauCOM.he1con[1][3] + phe1lv.he1cll[1][3]*phycon.eden)*c)*he1crt.he1lte[3]; z[4][9] = -(he1netCOM.he1net[1][4]*r + (he1tauCOM.he1con[1][4] + phe1lv.he1cll[1][4]*phycon.eden)*c)*he1crt.he1lte[4]; z[5][9] = -(he1netCOM.he1net[1][5]*r + (he1tauCOM.he1con[1][5] + phe1lv.he1cll[1][5]*phycon.eden)*c)*he1crt.he1lte[5]; z[6][9] = -(he1netCOM.he1net[1][6]*r + (he1tauCOM.he1con[1][6] + phe1lv.he1cll[1][6]*phycon.eden)*c)*he1crt.he1lte[6]; z[7][9] = -(he1netCOM.he1net[1][7]*r + (he1tauCOM.he1con[1][7] + phe1lv.he1cll[1][7]*phycon.eden)*c)*he1crt.he1lte[7]; z[8][9] = -(he1netCOM.he1net[1][8]*r + (he1tauCOM.he1con[1][8] + phe1lv.he1cll[1][8]*phycon.eden)*c)*he1crt.he1lte[8]; z[9][9] = he1cbr.he1radn[9] + phe1lv.he1gam[1] + (he1cbr.he1coldn[9] + phe1lv.he1clc[1])*phycon.eden + hex12 + he1crt.he1c2p2s*3.*edp + he1cbr.he1bul[1] + (che2sn.cs2s2s + che2sn.cs2s2p); /* LEVEL 3 */ z[0][2] = -(phe1lv.he1cll[0][2]*phycon.eden + (he1tauCOM.he1con[0][2]* he1statCOM.he1stat[2]/he1statCOM.he1stat[0])*he1crt.he1lte[0]/ he1crt.he1lte[2]); z[1][2] = -(phe1lv.he1cll[1][2]*phycon.eden + he1tauCOM.he1con[1][2]/ he1blt.he1bol[2][1])*0.75; z[2][2] = phe1lv.he1gam[2] + he1cbr.he1radn[2] + (phe1lv.he1clc[2] + he1cbr.he1coldn[2])*phycon.eden + he1cbr.he1bul[2]; z[3][2] = -(he1netCOM.he1net[2][3] + he1tauCOM.he1con[2][3] + phe1lv.he1cll[2][3]*phycon.eden)*he1crt.he1lte[3]/he1crt.he1lte[2]; z[4][2] = -(he1netCOM.he1net[2][4] + he1tauCOM.he1con[2][4] + phe1lv.he1cll[2][4]*phycon.eden)*he1crt.he1lte[4]/he1crt.he1lte[2]; z[5][2] = -(he1netCOM.he1net[2][5] + he1tauCOM.he1con[2][5] + phe1lv.he1cll[2][5]*phycon.eden)*he1crt.he1lte[5]/he1crt.he1lte[2]; z[6][2] = -(he1netCOM.he1net[2][6] + he1tauCOM.he1con[2][6] + phe1lv.he1cll[2][6]*phycon.eden)*he1crt.he1lte[6]/he1crt.he1lte[2]; z[7][2] = -(he1netCOM.he1net[2][7] + he1tauCOM.he1con[2][7] + phe1lv.he1cll[2][7]*phycon.eden)*he1crt.he1lte[7]/he1crt.he1lte[2]; z[8][2] = -(he1netCOM.he1net[2][8] + he1tauCOM.he1con[2][8] + phe1lv.he1cll[2][8]*phycon.eden)*he1crt.he1lte[8]/he1crt.he1lte[2]; z[9][2] = -(phe1lv.he1cll[1][2]*phycon.eden + he1tauCOM.he1con[1][2]/ he1blt.he1bol[2][1])*0.25; /* level 4 */ z[0][3] = -phe1lv.he1cll[0][3]*phycon.eden - he1tauCOM.he1con[0][3]/ he1blt.he1bol[3][0]; z[1][3] = -(phe1lv.he1cll[1][3]*phycon.eden + he1tauCOM.he1con[1][3]/ he1blt.he1bol[3][1])*0.75; z[2][3] = -phe1lv.he1cll[2][3]*phycon.eden - he1tauCOM.he1con[2][3]/ he1blt.he1bol[3][2]; z[3][3] = (he1cbr.he1coldn[3] + phe1lv.he1clc[3])*phycon.eden + he1cbr.he1radn[3] + phe1lv.he1gam[3] + he1cbr.he1bul[3]; z[4][3] = -(he1netCOM.he1net[3][4] + he1tauCOM.he1con[3][4] + phe1lv.he1cll[3][4]*phycon.eden)*he1crt.he1lte[4]/he1crt.he1lte[3]; z[5][3] = -(he1netCOM.he1net[3][5] + he1tauCOM.he1con[3][5] + phe1lv.he1cll[3][5]*phycon.eden)*he1crt.he1lte[5]/he1crt.he1lte[3]; z[6][3] = -(he1netCOM.he1net[3][6] + he1tauCOM.he1con[3][6] + phe1lv.he1cll[3][6]*phycon.eden)*he1crt.he1lte[6]/he1crt.he1lte[3]; z[7][3] = -(he1netCOM.he1net[3][7] + he1tauCOM.he1con[3][7] + phe1lv.he1cll[3][7]*phycon.eden)*he1crt.he1lte[7]/he1crt.he1lte[3]; z[8][3] = -(he1netCOM.he1net[3][8] + he1tauCOM.he1con[3][8] + phe1lv.he1cll[3][8]*phycon.eden)*he1crt.he1lte[8]/he1crt.he1lte[3]; z[9][3] = -(phe1lv.he1cll[1][3]*phycon.eden + he1tauCOM.he1con[1][3]/ he1blt.he1bol[3][1])*0.25; /* level 5 */ z[0][4] = -phe1lv.he1cll[0][4]*phycon.eden - he1tauCOM.he1con[0][4]/ he1blt.he1bol[4][0]; z[1][4] = -(phe1lv.he1cll[1][4]*phycon.eden + he1tauCOM.he1con[1][4]/ he1blt.he1bol[4][1])*0.75; z[2][4] = -phe1lv.he1cll[2][4]*phycon.eden - he1tauCOM.he1con[2][4]/ he1blt.he1bol[4][2]; z[3][4] = -phe1lv.he1cll[3][4]*phycon.eden - he1tauCOM.he1con[3][4]/ he1blt.he1bol[4][3]; z[4][4] = phe1lv.he1gam[4] + he1cbr.he1radn[4] + (he1cbr.he1coldn[4] + phe1lv.he1clc[4])*phycon.eden + he1cbr.he1bul[4]; z[5][4] = -(he1netCOM.he1net[4][5] + he1tauCOM.he1con[4][5] + phe1lv.he1cll[4][5]*phycon.eden)*he1crt.he1lte[5]/he1crt.he1lte[4]; z[6][4] = -(he1netCOM.he1net[4][6] + he1tauCOM.he1con[4][6] + phe1lv.he1cll[4][6]*phycon.eden)*he1crt.he1lte[6]/he1crt.he1lte[4]; z[7][4] = -(he1netCOM.he1net[4][7] + he1tauCOM.he1con[4][7] + phe1lv.he1cll[4][7]*phycon.eden)*he1crt.he1lte[7]/he1crt.he1lte[4]; z[8][4] = -(he1netCOM.he1net[4][8] + he1tauCOM.he1con[4][8] + phe1lv.he1cll[4][8]*phycon.eden)*he1crt.he1lte[8]/he1crt.he1lte[4]; z[9][4] = -(phe1lv.he1cll[1][4]*phycon.eden + he1tauCOM.he1con[1][4]/ he1blt.he1bol[4][1])*0.25; /* level 6 */ z[0][5] = -phe1lv.he1cll[0][5]*phycon.eden - he1tauCOM.he1con[0][5]/ he1blt.he1bol[5][0]; z[1][5] = -(phe1lv.he1cll[1][5]*phycon.eden + he1tauCOM.he1con[1][5]/ he1blt.he1bol[5][1])*0.75; z[2][5] = -phe1lv.he1cll[2][5]*phycon.eden - he1tauCOM.he1con[2][5]/ he1blt.he1bol[5][2]; z[3][5] = -phe1lv.he1cll[3][5]*phycon.eden - he1tauCOM.he1con[3][5]/ he1blt.he1bol[5][3]; z[4][5] = -phe1lv.he1cll[4][5]*phycon.eden - he1tauCOM.he1con[4][5]/ he1blt.he1bol[5][4]; z[5][5] = (he1cbr.he1coldn[5] + phe1lv.he1clc[5])*phycon.eden + phe1lv.he1gam[5] + he1cbr.he1radn[5] + he1cbr.he1bul[5]; z[6][5] = -(he1netCOM.he1net[5][6] + he1tauCOM.he1con[5][6] + phe1lv.he1cll[5][6]*phycon.eden)*he1crt.he1lte[6]/he1crt.he1lte[5]; z[7][5] = -(he1netCOM.he1net[5][7] + he1tauCOM.he1con[5][7] + phe1lv.he1cll[5][7]*phycon.eden)*he1crt.he1lte[7]/he1crt.he1lte[5]; z[8][5] = -(he1netCOM.he1net[5][8] + he1tauCOM.he1con[5][8] + phe1lv.he1cll[5][8]*phycon.eden)*he1crt.he1lte[8]/he1crt.he1lte[5]; z[9][5] = -(phe1lv.he1cll[1][5]*phycon.eden + he1tauCOM.he1con[1][5]/ he1blt.he1bol[5][1])*0.25; /* level 7 */ z[0][6] = -phe1lv.he1cll[0][6]*phycon.eden - he1tauCOM.he1con[0][6]/ he1blt.he1bol[6][0]; z[1][6] = -(phe1lv.he1cll[1][6]*phycon.eden + he1tauCOM.he1con[1][6]/ he1blt.he1bol[6][1])*0.75; z[2][6] = -phe1lv.he1cll[2][6]*phycon.eden - he1tauCOM.he1con[2][6]/ he1blt.he1bol[6][2]; z[3][6] = -phe1lv.he1cll[3][6]*phycon.eden - he1tauCOM.he1con[3][6]/ he1blt.he1bol[6][3]; z[4][6] = -phe1lv.he1cll[4][6]*phycon.eden - he1tauCOM.he1con[4][6]/ he1blt.he1bol[6][4]; z[5][6] = -phe1lv.he1cll[5][6]*phycon.eden - he1tauCOM.he1con[5][6]/ he1blt.he1bol[6][5]; z[6][6] = (he1cbr.he1coldn[6] + phe1lv.he1clc[6])*phycon.eden + phe1lv.he1gam[6] + he1cbr.he1radn[6] + he1cbr.he1bul[6]; z[7][6] = -(he1netCOM.he1net[6][7] + he1tauCOM.he1con[6][7] + phe1lv.he1cll[6][7]*phycon.eden)*he1crt.he1lte[7]/he1crt.he1lte[6]; z[8][6] = -(he1netCOM.he1net[6][8] + he1tauCOM.he1con[6][8] + phe1lv.he1cll[6][8]*phycon.eden)*he1crt.he1lte[8]/he1crt.he1lte[6]; z[9][6] = -(phe1lv.he1cll[1][6]*phycon.eden + he1tauCOM.he1con[1][6]/ he1blt.he1bol[6][1])*0.25; /* level 8 */ z[0][7] = -phe1lv.he1cll[0][7]*phycon.eden - he1tauCOM.he1con[0][7]/ he1blt.he1bol[7][0]; z[1][7] = -(phe1lv.he1cll[1][7]*phycon.eden + he1tauCOM.he1con[1][7]/ he1blt.he1bol[7][1])*0.75; z[2][7] = -phe1lv.he1cll[2][7]*phycon.eden - he1tauCOM.he1con[2][7]/ he1blt.he1bol[7][2]; z[3][7] = -phe1lv.he1cll[3][7]*phycon.eden - he1tauCOM.he1con[3][7]/ he1blt.he1bol[7][3]; z[4][7] = -phe1lv.he1cll[4][7]*phycon.eden - he1tauCOM.he1con[4][7]/ he1blt.he1bol[7][4]; z[5][7] = -phe1lv.he1cll[5][7]*phycon.eden - he1tauCOM.he1con[5][7]/ he1blt.he1bol[7][5]; z[6][7] = -phe1lv.he1cll[6][7]*phycon.eden - he1tauCOM.he1con[6][7]/ he1blt.he1bol[7][6]; z[7][7] = (he1cbr.he1coldn[7] + phe1lv.he1clc[7])*phycon.eden + phe1lv.he1gam[7] + he1cbr.he1radn[7] + he1cbr.he1bul[7]; z[8][7] = -(he1netCOM.he1net[7][8] + he1tauCOM.he1con[7][8] + phe1lv.he1cll[7][8]*phycon.eden)*he1crt.he1lte[8]/he1crt.he1lte[7]; z[9][7] = -(phe1lv.he1cll[1][7]*phycon.eden + he1tauCOM.he1con[1][7]/ he1blt.he1bol[7][1])*0.25; /* level 9 */ z[0][8] = -phe1lv.he1cll[0][8]*phycon.eden - he1tauCOM.he1con[0][8]/ he1blt.he1bol[8][0]; z[1][8] = -(phe1lv.he1cll[1][8]*phycon.eden + he1tauCOM.he1con[1][8]/ he1blt.he1bol[8][1])*0.75; z[2][8] = -phe1lv.he1cll[2][8]*phycon.eden - he1tauCOM.he1con[2][8]/ he1blt.he1bol[8][2]; z[3][8] = -phe1lv.he1cll[3][8]*phycon.eden - he1tauCOM.he1con[3][8]/ he1blt.he1bol[8][3]; z[4][8] = -phe1lv.he1cll[4][8]*phycon.eden - he1tauCOM.he1con[4][8]/ he1blt.he1bol[8][4]; z[5][8] = -phe1lv.he1cll[5][8]*phycon.eden - he1tauCOM.he1con[5][8]/ he1blt.he1bol[8][5]; z[6][8] = -phe1lv.he1cll[6][8]*phycon.eden - he1tauCOM.he1con[6][8]/ he1blt.he1bol[8][6]; z[7][8] = -phe1lv.he1cll[7][8]*phycon.eden - he1tauCOM.he1con[7][8]/ he1blt.he1bol[8][7]; z[8][8] = (he1cbr.he1coldn[8] + phe1lv.he1clc[8])*phycon.eden + phe1lv.he1gam[8] + he1cbr.he1radn[8] + he1cbr.he1bul[8]; z[9][8] = -(phe1lv.he1cll[1][8]*phycon.eden + he1tauCOM.he1con[1][8]/ he1blt.he1bol[8][1])*0.25; if( trace.lgTrace && trace.lgHe1Bug ) { for( i=0; i < (NHE1LVL + 1); i++ ) { fprintf( ioQQQ, " HeIs%2ld", i ); for( j=0; j < (NHE1LVL + 2); j++ ) { fprintf( ioQQQ, "%10.2e", z[j][i] ); } fprintf( ioQQQ, "\n" ); } } /* scale to largest number of order unity */ for( j=0; j < (NHE1LVL + 1); j++ ) { for( i=0; i < (NHE1LVL + 2); i++ ) { z[i][j] *= he1crt.he1lte[j]; } } /* which matrix solver? */ if( strcmp(TypMatrx.chMatrix,"matin1 ") == 0 ) { /*matin1();*/ nerror = 1; if( nerror != 0 ) { fprintf( ioQQQ, " HE1LEV MATRIX ERROR, STOP.\n" ); puts( "[Stop in he1lev]" ); exit(1); } } else if( strcmp(TypMatrx.chMatrix,"linpack") == 0 ) { /* this is the default matrix inverter */ for( j=0; j < (NHE1LVL + 1); j++ ) { for( i=0; i < (NHE1LVL + 1); i++ ) { amat[i][j] = z[i][j]; } bvec[j] = z[NHE1LVL+1][j]; } DGETRF(NHE1LVL+1,NHE1LVL+1,(double*)amat,NHE1LVL+1,ipiv,&nerror); if( nerror != 0 ) { fprintf( ioQQQ, " he1lvl dgetrf error\n" ); puts( "[Stop in he1lev]" ); exit(1); } DGETRS('N',NHE1LVL+1,1,(double*)amat,NHE1LVL+1,ipiv,bvec,NHE1LVL+1,&nerror); if( nerror != 0 ) { fprintf( ioQQQ, " he1lvl dgetrs error\n" ); puts( "[Stop in he1lev]" ); exit(1); } /* now put results back into z so rest of code treates only * one case - as if matin1 had been used */ for( i=0; i < (NHE1LVL + 1); i++ ) { z[NHE1LVL+1][i] = bvec[i]; } } else if( strcmp(TypMatrx.chMatrix,"veclib ") == 0 ) { /* Jason found this one on the Exemplar, distributed source just stops */ fprintf( ioQQQ, " this has not been checked since H atom conv\n" ); for( j=0; j < (NHE1LVL + 1); j++ ) { for( i=0; i < (NHE1LVL + 1); i++ ) { amat[i][j] = z[i][j]; } bvec[j] = z[NHE1LVL][j]; } job = 0; rcond = 0.; dgeco((double*)amat,NHE1LVL+1,NHE1LVL+1,ipiv,rcond,work); dgesl((double*)amat,NHE1LVL+1,NHE1LVL+1,ipiv,bvec,job); /* now put results back into z so rest of code treates only * one case - as if matin1 had been used */ for( i=0; i < (NHE1LVL + 1); i++ ) { z[NHE1LVL+1][i] = bvec[i]; } puts( "[Stop in he1lev]" ); exit(1); } else { fprintf( ioQQQ, " chMatrix type insane in hydrogen, was%7.7s\n", TypMatrx.chMatrix ); puts( "[Stop in he1lev]" ); exit(1); } /* convert from departure coef in Z(I) to level pop rel to HeII * HE1EX is net cooling rate due to col ionization-three body rec (s-1 here) */ he1bnCOM.he1bn[0] = z[10][0]; phe1lv.he1n[0] = (float)(he1bnCOM.he1bn[0]*he1crt.he1lte[0]*phycon.eden); he1bnCOM.he1bn[1] = z[10][1]; phe1lv.he1n[1] = (float)(he1bnCOM.he1bn[1]*he1crt.he1lte[1]*phycon.eden); for( i=2; i < NHE1LVL; i++ ) { he1bnCOM.he1bn[i] = z[10][i]; /* HE1N(I) are pop per level rel to He+; i.e., N(n)/N(He+) */ phe1lv.he1n[i] = (float)(he1bnCOM.he1bn[i]*he1crt.he1lte[i]*phycon.eden); } he1exCOM.he1ex = 0.; for( i=0; i < (NHE1LVL - 2); i++ ) { /* units of HE1EX are erg/cm^3/s */ he1exx[i] = He1NIonRyd.He1IonRyd[i]*phe1lv.he1clc[i]*phycon.eden* (phe1lv.he1n[i] - he1crt.he1lte[i]*phycon.eden); he1exCOM.he1ex += (float)he1exx[i]; } he1exCOM.he1ex *= (float)EN1RYD; phe1lv.he12p = phe1lv.he1n[1]; he1bnCOM.he1bn[NHE1LVL] = z[10][NHE1LVL]; phe1lv.he12s = (float)(he1bnCOM.he1bn[NHE1LVL]*he1crt.he1lte[NHE1LVL]* phycon.eden); phe1lv.he1n[NHE1LVL] = phe1lv.he12s; phe1lv.he1n[1] += phe1lv.he12s; *he2ov1 = phe1lv.he1n[0] + phe1lv.he1n[1] + phe1lv.he1n[2] + phe1lv.he1n[3] + phe1lv.he1n[4] + phe1lv.he1n[5]; /* HE2OV1 is He+ / He sin */ *he2ov1 = 1./ *he2ov1; he1rate.clhe2ov1 = (float)*he2ov1; he1exCOM.he1ex += (float)(EN1RYD*rfield.anu[nhe1Com.nhe1[1]-1]*phe1lv.he1clc[1]* phycon.eden*(phe1lv.he12s - he1crt.he1lte[NHE1LVL]*phycon.eden)); if( trace.lgHe1Bug && trace.lgTrace ) { fprintf( ioQQQ, " HE1LEV HE1GAM" ); for(i=0L; i < NHE1LVL + 1; i++) fprintf( ioQQQ, "%10.2e", phe1lv.he1gam[i] ); fprintf( ioQQQ, "\n" ); fprintf( ioQQQ, " HE1LEV TOTCAP" ); for(i=0L; i < 10; i++) fprintf( ioQQQ, "%10.2e", totcap[i] ); fprintf( ioQQQ, "%10.2e\n", herecCom.reci ); fprintf( ioQQQ, " HE1LEV IND Rc" ); for(i=0L; i < NHE1LVL + 1; i++) fprintf( ioQQQ, "%10.2e", he1ind.he1rin[i] ); fprintf( ioQQQ, "\n" ); fprintf( ioQQQ, " HE1LEV IND As" ); for(i=0L; i < NHE1LVL + 1; i++) fprintf( ioQQQ, "%10.2e", he1cbr.he1bul[i] ); fprintf( ioQQQ, "\n" ); fprintf( ioQQQ, " HE1LEV HE1LTE" ); for(i=0L; i < NHE1LVL + 1; i++) fprintf( ioQQQ, "%10.2e", he1crt.he1lte[i] ); fprintf( ioQQQ, "\n" ); fprintf( ioQQQ, " HE1LEV HE1N" ); for(i=0L; i < NHE1LVL + 1; i++) fprintf( ioQQQ, "%10.2e", phe1lv.he1n[i] ); fprintf( ioQQQ, "%10.2e%10.2e\n", phe1lv.he12s, phe1lv.he12p ); fprintf( ioQQQ, " HE1LEV b(n)" ); for( i=0; i < NHE1LVL; i++ ) { fprintf( ioQQQ, "%10.2e", he1bnCOM.he1bn[i] ); } fprintf( ioQQQ, "%10.2e%10.2e\n", he1bnCOM.he1bn[9], he1bnCOM.he1bn[1] ); fprintf( ioQQQ, " HE1LEV X12tot%10.2e Heat Exc%10.2e Grn dest:%10.2e\n", Secondaries.x12tot, he1exCOM.he1ex, gdist ); } for( i=0; i < (NHE1LVL + 1); i++ ) { if( phe1lv.he1n[i] <= 0. ) { fprintf( ioQQQ, " HE1LEV finds negative population,%3ld%10.2e simple=%" "10.2e Te=%10.2e NZONE=%4ld\n", i, phe1lv.he1n[i], *sm2ov1, phycon.te, ZoneCnt.nzone ); fprintf( ioQQQ, " HE1LEV pops," ); for( j=0; j < (NHE1LVL + 1); j++ ) { fprintf( ioQQQ, "%9.1e", phe1lv.he1n[j] ); } fprintf( ioQQQ, "\n" ); negcon(); ShowMe(); puts( "[Stop in he1lev]" ); exit(1); } } if( trace.lgTrace ) { fprintf( ioQQQ, " HE1LEV returns, mtrx finds HE2OV1 =%10.2e simple=%10.2e Gamma1=%10.2e rec1%10.2e rec3%10.2e\n", *he2ov1, *sm2ov1, phe1lv.he1gam[0], herecCom.reci, herecCom.rectri ); } # ifdef DEBUG_FUN fputs( " <->he1lev()\n", debug_fp ); # endif return; }