/*hmole determine populations of hydrogen molecules */ #include #include #include "cddefines.h" #define CONST_ 9.940622e26 #ifdef NMOLE # undef NMOLE #endif #define NMOLE 5 #include "powi.h" #include "showme.h" #include "physconst.h" #include "nhe1lvl.h" #include "grainvar.h" #include "opacpoint.h" #include "heat.h" #include "secondaries.h" #include "colden.h" #include "opac.h" #include "logte.h" #include "heavy.h" #include "h.h" #include "rfield.h" #include "tepowers.h" #include "hmi.h" #include "heating.h" #include "iph2.h" #include "timesc.h" #include "hydrogenic.h" #include "trace.h" #include "nomole.h" #include "hehp.h" #include "he.h" #include "phycon.h" #include "rh2dis.h" #include "hmihet.h" #include "iphmin.h" #include "hmidep.h" #include "nh.h" #include "h2max.h" #include "thlo.h" #include "ghabng.h" #include "gionrc.h" #include "h2dish.h" #include "zonecnt.h" #include "doppvel.h" #include "hahmin.h" #include "hmollte.h" #include "h2opac.h" #include "nhe1.h" #include "hphoto.h" #include "typmatrx.h" #include "linpack.h" #include "sexp.h" #include "rtescprob.h" #include "matin1.h" #include "gammas.h" #include "veclib.h" #include "hmrate.h" #include "hmirat.h" #include "hmole.h" void hmole( /* hmovhn1 is ratio of molecular to neutral H, * is product of this routine*/ float *hmovh1, /* this is the ratio of ion to atom for hydrogen, as produced * in the ionization balance */ double h2ovh1) { char chLab[NMOLE][5]; /* will punch debug output to this file */ FILE*ioFile; long int i, ih2, ih2p, ih3p, ihmi, iho, ipConServ, ipOffSet, ipiv[NMOLE], j, job, limit, merror, nd; double batach, bh2dis, bh2h2p, bhneut, c3bod, cionhm, corr, damper, desh2p, drate, eh2hh, ex3hp, exph2, exphmi, exphp, f, fa, fac, faneut, fhneut, fmol, g, gamhd, gamheh, h1fnd, h1good, h2pcin, h2ph3p, h2phhp, h2pion, h32h2, h3petc, h3ph2p, hnnew, ph3lte, radasc, radath, ratach, rate, renorm, rh2h2p, tdlow; static double colind, eh2hhm, eh2old, gamhmi, gamtwo, gh2dis, h2esc, h2phet, hmiht, hneut, ratind, rgrain, rspon, th2; double amat[NMOLE][NMOLE], b2pcin, bvec[NMOLE], c[NMOLE + 1][NMOLE + 1], plte, rcond, work[NMOLE]; # ifdef DEBUG_FUN fputs( "<+>hmole()\n", debug_fp ); # endif /* following used for linpac matrix inversion */ /* CONST is 8\pi * nu(912)**3 / c**2 */ /* there are two "no molecules" options, the no co, which turns off everything, * and the no n2, which only turns off the h2. in order to not kill the co * part we still need to compute the hydrogen network here, and then set h2 to * small values */ if( nomole.lgNoCOMole || phycon.te > 1e5 ) { /* thtmol = 0. */ *hmovh1 = 0.; hmi.hminus = 0.; hmihetCom.hmihet = 0.; h2opac.H2Opacity = 0.; hmiht = 0.; hmihetCom.hmicol = 0.; HeatingCom.heating[15][0] = 0.; HeatingCom.heating[16][0] = 0.; HeatingCom.heating[17][0] = 0.; h2dish.HeatH2Dish = 0.; hmidepCom.hmidep = 1.; gamhmi = 0.; hmi.htwo = 0.; hmi.h2plus = 0.; hehpCom.hehp = 0.; rh2disCom.rh2dis = 0.; gionrc.GrnIonRec = 0.; hahmin.HalphaHmin = 0.; # ifdef DEBUG_FUN fputs( " <->hmole()\n", debug_fp ); # endif return; } /* make sure ion ratio is positive */ assert( h2ovh1 >= 0. ); /* population of H- in LTE * IP is 0.754 eV * LTE population of H minus - cm^3 */ exphmi = sexp(8.745e3/phycon.te); if( exphmi > 0. ) { /* these are ratio n*(H-)/[ n*(ne) n*(Ho) ] */ HMolLTE.phmlte = SAHA/(TePowers.te32*exphmi)*(1./(2.*2.)); } else { HMolLTE.phmlte = 0.; } /* population of H2 in LTE * dissociation energy is 4.477eV */ exph2 = sexp(5.195e4/phycon.te); if( exph2 > 0. ) { /* extra factor accounts for mass of H instead of e- in SAHA * last factor was put in ver 85.23, missing before */ if( ZoneCnt.nzone <= 0 ) { damper = 0.; } else { damper = 0.5; } HMolLTE.ph2lte = (1. - damper)*SAHA/(TePowers.te32*exph2)* (1./(2.*2.))*3.634e-5 + damper*HMolLTE.ph2lte; } else { HMolLTE.ph2lte = 0.; } /* population of H2+ in LTE, phplte is H_2+/H / H+ * dissociation energy is 2.647 */ exphp = sexp(3.072e4/phycon.te); if( exphp > 0. ) { /* stat weight of H2+ is 4 * last factor was put in ver 85.23, missing before */ HMolLTE.phplte = SAHA/(TePowers.te32*exphp)*(4./(2.*1.))*3.634e-5; } else { HMolLTE.phplte = 0.; } /* population of H3+ in LTE, ph3lte is H_3+/( H2+ H+ ) * dissociation energy is 2.647 */ ex3hp = sexp(1.882e4/phycon.te); if( ex3hp > 0. ) { /* stat weight of H2+ is 4 * last factor was put in ver 85.23, missing before */ ph3lte = SAHA/(TePowers.te32*ex3hp)*(4./(2.*1.))*3.634e-5; } else { ph3lte = 0.; } /* first do gamma function for H minus, with spon+induc rec */ rspon = 0.; ratind = 0.; gamhmi = 0.; hmiht = 0.; colind = 0.; hmihetCom.hmicol = 0.; ipOffSet = iphminCom.iphmin + OpacPoint.iphmop; limit = MIN2(nhe1Com.nhe1[0],rfield.nflux); /* convert this into GammaBn someday, instead of done in a loop */ for( i=iphminCom.iphmin-1; i < limit; i++ ) { f = rfield.SummedCon[i]*opac.OpacStack[i+ipOffSet]; gamhmi += f; hmiht += f*rfield.anu[i]; ratind += f*rfield.ContBoltz[i]; /* induced cooling */ colind += f*rfield.ContBoltz[i]*rfield.anu[i]; g = opac.OpacStack[i-iphminCom.iphmin+OpacPoint.iphmop]* rfield.ContBoltz[i]*rfield.widflx[i]*rfield.anu2[i]; rspon += g; hmihetCom.hmicol += g*rfield.anu[i]; } /* correct for work function, convert from Ryd to erg */ hmiht = (hmiht - gamhmi*0.055502)*EN1RYD; colind = (colind - ratind*0.055502)*EN1RYD; hmihetCom.hmicol = (hmihetCom.hmicol - rspon*0.055502)*EN1RYD; /* add on high energy ionization, assume hydrogen cross section * n.b.; HGAMNC with secondaries */ gamhmi += HPhoto.hgamnc[0][IP1S]; hmihetCom.hmicol = (hmihetCom.hmicol*CONST_)*(HMolLTE.phmlte* phycon.eden)*h.hi; ratind *= HMolLTE.phmlte*phycon.eden; colind *= HMolLTE.phmlte*phycon.eden*h.hi; /* radiative association rate, spon and induced */ if( phycon.te >= 1e4 ) { rspon = (rspon*CONST_)*(HMolLTE.phmlte*phycon.eden); } else { rspon = hmirat(phycon.te)*phycon.eden; } /* photodissociation by Lyman band absorption: esc prob treatment, * treatment based on Tielens and Hollenbach 1985 ApJ 291, 722. */ ghabng.GammaHabing = 0.; /* find total intensity over carbon-ionizing continuum */ for( i=Heavy.ipHeavy[0][5]-1; i < nh.ipHn[0][IP1S]; i++ ) { ghabng.GammaHabing += (rfield.flux[i] + rfield.outcon[i] + rfield.outlin[i])*rfield.anu[i]; } /* now convert to Habing ISM units * GammaHabing is FUV continuum relative to Habing value */ ghabng.GammaHabing = (float)(ghabng.GammaHabing*2.18e-11/1.6e-3); /* cross section is eqn A8 of Teilens and Hollenbach 85a * branching ratio of 10% in, so like their table 5 */ gh2dis = ghabng.GammaHabing*3.4e-11; /* escape prob is line shielding, eqn A10 */ h2opac.H2Opacity = (float)(1.2e-14*(1e5/DoppVel.doppler[0])); th2 = coldenCom.colden[IPCOLH2-1]*1.2e-14*(1e5/DoppVel.doppler[0]); h2esc = escinc(th2,1e-4); /* H_2 + h nu => 2H */ gh2dis *= h2esc; /* collisional ionization of H-, rate from Janev, Langer et al. */ if( phycon.te < 3074. ) { cionhm = 1.46e-32*(powi(phycon.te,6))*TePowers.sqrte*exphmi; } else if( phycon.te >= 3074. && phycon.te < 30000. ) { cionhm = 5.9e-19*phycon.te*phycon.te*TePowers.sqrte*TePowers.te03* TePowers.te01*TePowers.te01; } else { cionhm = 3e-7; } /* H2 formation on grains; * rate from Hollenback and McKee Ap.J.Sup 41, 555 eq 3.4 3.8 */ if( GrainVar.lgDustOn ) { drate = 0.; tdlow = FLT_MAX; for( nd=0; nd < NDUST; nd++ ) { if( GrainVar.lgDustOn1[nd] && GrainVar.tedust[nd] > 0. ) { tdlow = MIN2(tdlow,GrainVar.tedust[nd]); drate = GrainVar.dstfactor[nd]; } } fa = 1./(1. + 1e4*sexp(600./tdlow)); rgrain = 3e-18*TePowers.sqrte*fa/(1. + 0.04*sqrt(tdlow+ phycon.te) + 0.002*phycon.te + 8e-6*phycon.te*phycon.te)* drate; /* rate ions recombine on grain surface - used elsewhere * based on Drain and Sutin 1987 ApJ 320, 803 eqn 5.15 * GRECON is usually 1, option to turn this process off * >>chng 97 feb 24, following had SQRT( tdlow ) not sqrte, * caught by Jon Slavin */ gionrc.GrnIonRec = (float)(5.8e-13/TePowers.sqrte*drate*phycon.hden* gionrc.grecon); } else { rgrain = 0.; gionrc.GrnIonRec = 0.; } /* collisional dissociation, rate from Dove and Mandy, Ap.J.(Let) 311, L93. * corr is correction for approach to high density limit * H2 + H => 3H - rate very uncertain */ corr = MIN2(6.,14.44-logte.alogte*3.08); corr = pow(10.,MAX2(0.,corr)/(1. + 1.6e4/h.hi)); rh2disCom.rh2dis = (float)(1.55e-8/TePowers.sqrte*sexp(65107./phycon.te)* corr); /* H2 + H+ => H2+ + H *>>chng 98 jan 02, from 2.12e4 to 2.123e4 */ rh2h2p = 1.8e-12*TePowers.sqrte*TePowers.te10/TePowers.te01*sexp(2.123e4/ phycon.te); /*>>chng 98 jan 02, following had used ratio of lte pops to cancel above sexp * inverse rate */ bh2h2p = 1.8e-12*TePowers.sqrte*TePowers.te10/TePowers.te01*2./16.; /* H2+ + HNU => H+ + H */ gamtwo = GammaK(OpacPoint.ih2pnt[0],OpacPoint.ih2pnt[1],OpacPoint.ih2pof,1.); h2phet = heat.HeatNet; /* molecule formation */ for( i=0; i < NMOLE; i++ ) { c[NMOLE][i] = 0.; for( j=0; j < NMOLE; j++ ) { c[j][i] = 0.; } } hneut = h.hi + hmi.hminus + 2.*(hmi.htwo + hmi.h2plus) + 3.*hmi.h3plus; /* these should be totally equivalent?? * >>chng 96 june 3, hneut floated above hden, ended negative */ hneut = phycon.hden - h.hii; /* protect against roundoff in ionized gas */ if( hneut/phycon.hden < 1e-3 ) { hneut = h.hi + hmi.hminus + 2.*(hmi.htwo + hmi.h2plus) + 3.* hmi.h3plus; } iho = 1; strcpy( chLab[iho-1], "HO " ); ihmi = 2; strcpy( chLab[ihmi-1], "H- " ); ih2 = 3; strcpy( chLab[ih2-1], "H2 " ); ih2p = 4; strcpy( chLab[ih2p-1], "H2+ " ); ih3p = 5; strcpy( chLab[ih3p-1], "H3+ " ); /*-------------------------------------------------------------------- */ /* H- H minus hminus balance equations * (IHMI,IHO) == processes making H- from Ho =+sign * radiative attachment: HI + NE => H-; */ c[iho-1][ihmi-1] += rspon + ratind; c[iho-1][iho-1] += -rspon - ratind; /* (IHMI,IHMI) = processes destroying H- =-sign * photodissociation, H- + H NU => H + NE */ c[ihmi-1][ihmi-1] -= gamhmi; c[ihmi-1][iho-1] += gamhmi; /* mutual neutralization with heavies, rate from Dalgarno and Mcray * all charged ions contribute equally */ faneut = 4e-6/TePowers.sqrte*MAX2(0.,phycon.eden-h.hii-he.heii- 2.*he.heiii); c[ihmi-1][ihmi-1] -= faneut; /* electron collisional ionization of H- */ c[ihmi-1][ihmi-1] += -cionhm*phycon.eden; c[ihmi-1][iho-1] += cionhm*phycon.eden; /* inverse process; three body rec */ c3bod = cionhm*(HMolLTE.phmlte*phycon.eden)*phycon.eden; c[iho-1][ihmi-1] += c3bod; c[iho-1][iho-1] -= c3bod; /* associative detachment: H- + H => H2 + E */ ratach = h.hi*1.35e-9; c[ihmi-1][ihmi-1] -= ratach; c[iho-1][iho-1] += -hmi.hminus*1.35e-9; /* (1,IH2) convert H2 into H- = +sign * the back reaction, H2 + e => H- + Ho */ if( HMolLTE.ph2lte > 0. ) { batach = (ratach/h.hi)*HMolLTE.phmlte/HMolLTE.ph2lte*phycon.eden; } else { batach = 0.; } c[ih2-1][ihmi-1] += batach; c[ih2-1][iho-1] += batach; /* mutual neut, mostly into n=3; rates from Janev et al * H- + H+ => H + H(n=3) */ fhneut = h.hii*7e-8; c[ihmi-1][ihmi-1] -= fhneut; c[ihmi-1][iho-1] += fhneut; /* back reaction from excited state H */ if( phycon.te > thlo.HydTempLimit ) { /* HBN(3,1) is defined; when =4) + e- => H + H-, * from Janev et al, quoted in lenzuni etal * density dep is for non-lte, guess from dalgarno and roberge apl 233, 25 * extra expo factor added for low temps */ if( ZoneCnt.nzone <= 1 ) { /* this is initial setup of code, so set rate coef to actual val */ eh2hhm = 2.7e-8*pow((double)10,-0.7*POW2(logte.alogte - 3.615))* phycon.eden*(phycon.hden/(1e7 + phycon.hden))*sexp(52000./ phycon.te); eh2old = eh2hhm; } else { /* this is deeper into the cloud, and there is danger of oscillation */ eh2old = eh2hhm; eh2hhm = 2.7e-8*pow((double)10,-0.7*POW2(logte.alogte - 3.615))* phycon.eden*(phycon.hden/(1e7 + phycon.hden))*sexp(52000./ phycon.te); fac = 0.5; eh2hhm = eh2hhm*fac + eh2old*(1. - fac); } c[ih2-1][ihmi-1] += eh2hhm; c[ih2-1][iho-1] += eh2hhm; /*-------------------------------------------------------------------- * * molecular hydrogen H2 htwo balance equation * (IH2,IHO)==create H2 from Ho =+ * * H2 formation on grains; * rate coefficient from Hollenback and McKee Ap.J.Sup 41, 555 eq 3.4 3.8 * really is rate coefficient, needs two HIs */ c[iho-1][ih2-1] += rgrain*h.hi; c[iho-1][iho-1] += -rgrain*h.hi; /* excited atom radiative association, * H(n=2) + H(n=1) => H2 + hnu * from Latter, W.B., and black, J.H., 1991, Ap.J. 372, 161 */ radasc = ((hydro.hn[0][IP2P] + hydro.hn[0][IP2S])*h.hii)*3e-14; c[iho-1][ih2-1] += radasc; c[iho-1][iho-1] -= radasc; /* make H2 from H- =+ sign * associative detachment; H- + H => H2: Browne and Dalgarno J PHys B 2, 885 */ c[ihmi-1][ih2-1] += ratach; /* photo-destroy H2 */ c[ih2-1][ih2-1] -= gh2dis; c[ih2-1][iho-1] += gh2dis; /* the process H2 + e- => H + H + e * from Lenzuni et al. apj sup 76, 759, quoted from Janev et al. */ eh2hh = 1.3e-18*phycon.te*phycon.te*sexp(52000./phycon.te)*phycon.eden; c[ih2-1][ih2-1] -= eh2hh; c[ih2-1][iho-1] += eh2hh; /* the processes H2(v>=4) + e- => H + H- * from Lenzuni et al. apj sup 76, 759, quoted from Janev et al. * density dep is for non-lte, guess from dalgarno and roberge apl 233, 25 */ c[ih2-1][ih2-1] -= eh2hhm; /* add on inverse of radiative attachment */ c[ih2-1][ih2-1] -= batach; /* collisional dissociation, rate from Dove and Mandy, Ap.J.(Let) 311, L93. * H_2 + S => 2H + S */ c[ih2-1][ih2-1] += -rh2disCom.rh2dis*h.hi; c[ih2-1][iho-1] += rh2disCom.rh2dis*h.hi; /* back rate, three body recombination, 2H + S => H_2 + S */ bh2dis = rh2disCom.rh2dis*HMolLTE.ph2lte*h.hi*h.hi; c[iho-1][ih2-1] += bh2dis; c[iho-1][iho-1] -= bh2dis; /* H2 + HNU=> H2+ + E * photoionization by hard photons, crossection=3*HI */ gamhd = HPhoto.hgamnc[0][IP1S]; c[ih2-1][ih2-1] -= gamhd; /* assume cosmic rays do the same thing */ c[ih2-1][ih2-1] -= Secondaries.csupra; /* H2 + H+ => H2+ + H */ c[ih2-1][ih2-1] += -rh2h2p*h.hii; c[ih2-1][iho-1] += rh2h2p*h.hii; /* H2+ + H => H+ + H2 */ c[ih2p-1][ih2-1] += bh2h2p*h.hi; c[iho-1][iho-1] += -bh2h2p*hmi.h2plus; /* H + H3+ => H2 + H2+ */ h3ph2p = hmrate(2.08e-9,0.,1.88e4)*h.hi; c[ih3p-1][ih2-1] += h3ph2p; c[iho-1][iho-1] += -hmrate(2.08e-9,0.,1.88e4)*hmi.h3plus; /* H2 + H3+ => H2 + H2+ + H */ h3petc = hmrate(3.41e-11,0.5,7.16e4)*hmi.htwo; c[ih3p-1][ih2-1] += h3petc; c[ih2-1][iho-1] += h3petc; /* H2 + H3+ => H2 + H+ + H2 */ h32h2 = hmrate(3.41e-11,0.5,5.04e4)*hmi.htwo; c[ih3p-1][ih2-1] += 2.*h32h2; /* e + H3+ => H2 + H */ c[ih3p-1][ih2-1] += hmrate(5.00e-9,-0.5,0.)*phycon.eden; c[ih3p-1][iho-1] += hmrate(5.00e-9,-0.5,0.)*phycon.eden; if( (trace.lgTrace && trace.lgTrMole) || iph2.lgPunH2 ) { if( iph2.lgPunH2 ) { ioFile = iph2.ipPunH2; } else { ioFile = ioQQQ; } if( c[ih2-1][ih2-1] != 0. ) { rate = -c[ih2-1][ih2-1]; fprintf( ioFile, " Destroy H2: rate=%10.2eDIS;%5.3f bat;%5.3f h2dis;%5.3f gamhd;%5.3f h2h2p;%5.3f E-h;%5.3f E-h-;%5.2f\n", rate, gh2dis/rate, batach/rate, rh2disCom.rh2dis*h.hi/ rate, gamhd/rate, rh2h2p*h.hii/rate, eh2hh/rate, eh2hhm/ rate ); } else { fprintf( ioFile, " Destroy H2: rate=0\n" ); } } /*------------------------------------------------------------------- */ /* h2plus H2+ balance equations */ /* make H2+ from Ho * H+ + H => H2+ + HNU * approximation was from Kurucz thesis, not meant for hot gas */ radath = MAX2(0.,2.325*MIN2(5000.,phycon.te)-1375.)*1e-20; c[iho-1][ih2p-1] += radath*h.hii; c[iho-1][iho-1] += -radath*h.hii; /* make H2+ from H2 =+sign * H2 + HNU => H2+ + E * also cosmic rays */ c[ih2-1][ih2p-1] += gamhd + Secondaries.csupra; /* H2 + H+ => H2+ + H */ c[ih2-1][ih2p-1] += rh2h2p*h.hii; /* (3,IH2P) == destroy H2+ = -sign * H + H2+ => H+ + H2 */ c[ih2p-1][ih2p-1] += -bh2h2p*h.hi; c[iho-1][iho-1] += -bh2h2p*hmi.h2plus; /* H2+ + p => H + H+ + H+; Janev et al. 3.2.6 */ h2pion = 2.4e-27*POW3(phycon.te)*h.hii; c[ih2p-1][ih2p-1] -= h2pion; c[ih2p-1][iho-1] += h2pion; /* H2+ + E => H + H+ + e-; Janev et al. */ h2pcin = 2e-7*sexp(30720./phycon.te); c[ih2p-1][ih2p-1] += -h2pcin*phycon.eden; c[ih2p-1][iho-1] += h2pcin*phycon.eden; /* back reaction, H + H+ + e => h2+ + e */ b2pcin = h2pcin*HMolLTE.phplte; /* this is the hot reaction at high densities */ c[iho-1][ih2p-1] += b2pcin*phycon.eden*h.hii; c[iho-1][iho-1] += -b2pcin*phycon.eden*h.hii; /* H2+ + HNU => H+ + H */ c[ih2p-1][ih2p-1] -= gamtwo; /* photoionization by hard photons, crossection =2*HI (wild guess) */ c[ih2p-1][ih2p-1] += -2.*HPhoto.hgamnc[0][IP1S]; /* H2 + H2+ => H + H3+ */ h2ph3p = 1.40e-9*hmi.htwo*(1. - sexp(9940./phycon.te)); c[ih2p-1][ih2p-1] -= h2ph3p; c[ih2p-1][iho-1] += h2ph3p; /* h + h3+ => h2 + h2+ */ c[ih3p-1][ih2p-1] += h3ph2p; /* H2 + H3+ => H2 + H2+ + H */ c[ih3p-1][ih2p-1] += h3petc; /* destroy H2+ via H2+ + H2 => H + H+ + H2 */ h2phhp = 2.41e-12*TePowers.sqrte*sexp(30720./phycon.te)*hmi.htwo; c[ih2p-1][ih2p-1] -= h2phhp; c[ih2p-1][iho-1] += h2phhp; /* save total H2P destruction rate for possible later printout: * NB this must come last */ desh2p = -c[ih2p-1][ih2p-1]; /*------------------------------------------------------------------ */ /* H3+ balance equations*/ /* H2 + H2+ => H + H3+ */ c[ih2p-1][ih3p-1] += h2ph3p; /* H + H3+ => H2 + H2+ */ c[ih3p-1][ih3p-1] -= h3ph2p; /* H2 + H3+ => H2 + H2+ + H */ c[ih3p-1][ih3p-1] -= h3petc; /* H2 + H3+ => H2 + H+ + H2 */ c[ih3p-1][ih3p-1] -= h32h2; /* e + H3+ => 3H * e + H3+ => H2 + H */ c[ih3p-1][ih3p-1] += -2.*hmrate(5.00e-9,-0.5,0.)*phycon.eden; /* photoionization by hard photons, crossection =2*HI (wild guess) */ c[ih3p-1][ih3p-1] += -2.*HPhoto.hgamnc[0][IP1S]; /*------------------------------------------------------------------ */ /* first equation is hydrogen conservation * remember that htwo/hden = 0.5 when totally molecular */ if( hmi.htwo/phycon.hden > 0.1 ) { ipConServ = ih2; } else { /* little molecular hydrogen, leave out the atomic */ ipConServ = iho; } /* this was not the cause of the problem, go back to old way * until I have time to debug this logic */ ipConServ = 1; c[iho-1][ipConServ-1] = 1.; c[ihmi-1][ipConServ-1] = 1.; c[ih2-1][ipConServ-1] = 2.; c[ih2p-1][ipConServ-1] = 2.; c[ih3p-1][ipConServ-1] = 3.; /* the above will add up to the following */ c[NMOLE][ipConServ-1] = hneut; /*------------------------------------------------------------------ */ if( trace.lgTrace && trace.lgTrMole ) { for( i=0; i < NMOLE; i++ ) { fprintf( ioQQQ, " MOLE%2ld", i ); for( j=0; j < (NMOLE + 1); j++ ) { fprintf( ioQQQ, "%10.2e", c[j][i] ); } fprintf( ioQQQ, "\n" ); } } /* establish age timescale for molecule formation */ if( -c[ih2-1][ih2-1] > SMALLFLOAT ) { timesc.AgeH2MoleDest = -1./c[ih2-1][ih2-1]; } else { timesc.AgeH2MoleDest = 0.; } /* invert matrix * which matrix solver? */ if( strcmp(TypMatrx.chMatrix,"matin1 ") == 0 ) { /*matin1();*/ merror = 0; if( merror != 0 ) { fprintf( ioQQQ, " HMOLE matrix error, stop.\n" ); ShowMe(); puts( "[Stop in hmole]" ); exit(1); } } else if( strcmp(TypMatrx.chMatrix,"linpack") == 0 ) { /* this one may be more robust */ for( j=0; j < NMOLE; j++ ) { for( i=0; i < NMOLE; i++ ) { amat[i][j] = c[i][j]; } bvec[j] = c[NMOLE][j]; } DGETRF(NMOLE,NMOLE,(double*)amat,NMOLE,ipiv, &merror); DGETRS('N',NMOLE,1,(double*)amat,NMOLE,ipiv,bvec,NMOLE,&merror); if( merror != 0 ) { fprintf( ioQQQ, " hmole dgetrs error\n" ); puts( "[Stop in hmole]" ); 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 < NMOLE; i++ ) { c[NMOLE][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 < NMOLE; j++ ) { for( i=0; i < NMOLE; i++ ) { amat[i][j] = c[i][j]; } bvec[j] = c[NMOLE][j]; } job = 0; rcond = 0.; dgeco((double*)amat,NMOLE,NMOLE,ipiv,rcond,work); dgesl((double*)amat,NMOLE,NMOLE,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 < NMOLE; i++ ) { c[NMOLE][i] = bvec[i]; } puts( "[Stop in hmole]" ); exit(1); } else { fprintf( ioQQQ, " chMatrix type insane in hydrogen, was%7.7s\n", TypMatrx.chMatrix ); puts( "[Stop in hmole]" ); exit(1); } h1fnd = c[NMOLE][iho-1]; hmi.hminus = (float)c[NMOLE][ihmi-1]; hmi.h2plus = (float)c[NMOLE][ih2p-1]; hmi.h3plus = (float)c[NMOLE][ih3p-1]; hmi.htwo = (float)c[NMOLE][ih2-1]; /* rate H-alpha is created by H- ct */ hahmin.HalphaHmin = (float)(fhneut*hmi.hminus); /* identify creation and destruction processes for H2+ */ if( trace.lgTrace && trace.lgTrMole ) { rate = desh2p; if( rate != 0. ) { fprintf( ioQQQ, " Destroy H2+: rate=%10.2e e-;%5.3f phot;%5.3f hard gam;%5.3f H2col;%5.3f h2phhp;%5.3f pion;%5.3f bh2h2p:%5.3f\n", rate, h2pcin*phycon.eden/rate, gamtwo/rate, 2.*HPhoto.hgamnc[0][IP1S]/ rate, h2ph3p/rate, h2phhp/rate, h2pion/rate, bh2h2p* h.hi/rate ); rate *= hmi.h2plus; if( rate > 0. ) { fprintf( ioQQQ, " Create H2+: rate=%10.2e HII HI;%5.3f Col H2;%5.3f HII H2;%5.3f HI HI;%5.3f\n", rate, radath*h.hii*h.hi/rate, (gamhd + Secondaries.csupra)* hmi.htwo/rate, rh2h2p*h.hii*hmi.htwo/rate, b2pcin* h.hi*h.hii*phycon.eden/rate ); } else { fprintf( ioQQQ, " Create H2+: rate= is zero\n" ); } } } /* heating due to H2 dissociation */ if( nomole.lgNoH2Mole ) { h2dish.HeatH2Dish = 0.; HeatingCom.heating[17][0] = 0.; } else { h2dish.HeatH2Dish = (float)(1.36e-23*hmi.htwo*h2esc*ghabng.GammaHabing); HeatingCom.heating[17][0] = h2dish.HeatH2Dish; } *hmovh1 = (float)(((hmi.htwo + hmi.h2plus)*2. + 3.*hmi.h3plus + hmi.hminus) / h.hi); if( *hmovh1 < 0. ) { fprintf( ioQQQ, " HMOLE: I found a negative molecular fraction=%10.2e\n", *hmovh1 ); fprintf( ioQQQ, " HMOLE: htwo,h2plus,h3plus,hminus,hi,eden=%9.1e%9.1e%9.1e%9.1e%9.1e%9.1e\n", hmi.htwo, hmi.h2plus, hmi.h3plus, hmi.hminus, h.hi, phycon.eden ); fprintf( ioQQQ, " HMOLE: h1fnd=%10.2e\n", h1fnd ); ShowMe(); puts( "[Stop in hmole]" ); exit(1); } h1good = phycon.hden/(1.e0 + (double)(*hmovh1) + h2ovh1); hnnew = h1good*(1.e0 + (double)(*hmovh1)); renorm = hnnew/hneut; if( trace.lgTrace && trace.lgTrMole ) { fprintf( ioQQQ, " raw; " ); for( i=1; i <= NMOLE; i++ ) { fprintf( ioQQQ, "%4.4s:%10.2e", chLab[i-1], c[NMOLE][i-1] ); } fprintf( ioQQQ, " \n" ); fprintf( ioQQQ, " ren; " ); for( i=1; i <= NMOLE; i++ ) { fprintf( ioQQQ, "%4.4s:%10.2e", chLab[i-1], c[NMOLE][i-1]* renorm ); } fprintf( ioQQQ, " \n" ); fprintf( ioQQQ, " RENORM factor was%10.2e\n", renorm ); } h1fnd *= renorm; hmi.hminus *= (float)renorm; hmi.htwo *= (float)renorm; hmi.h2plus *= (float)renorm; hmi.h3plus *= (float)renorm; /* following not used any more */ /*hatmic = (h.hi + h.hii)/(h.hi + h.hii + hmi.hminus + 2.*(hmi.htwo + hmi.h2plus) + 3*hmi.h3plus);*/ if( (trace.lgTrace && trace.lgTrMole) || iph2.lgPunH2 ) { if( iph2.lgPunH2 ) { ioFile = iph2.ipPunH2; } else { ioFile = ioQQQ; } rate = rgrain*h.hi + ratach*hmi.hminus + bh2dis*h.hi + bh2h2p* h.hi*hmi.h2plus + radasc*h.hi + h3ph2p*hmi.h3plus + h3petc* hmi.h3plus; fprintf( ioFile, " Create H2, rate=%10.2e grain;%5.3f hmin;%5.3f bhedis;%5.3f h2+;%5.3f radasc:%5.3f h3ph2p:%5.3f h3petc:%5.3f\n", rate, rgrain*h.hi/rate, ratach*hmi.hminus/rate, bh2dis/rate, bh2h2p*h.hi*hmi.h2plus/rate, radasc*h.hi/rate, h3ph2p*hmi.h3plus/ rate, h3petc*hmi.h3plus/rate ); } if( trace.lgTrace && trace.lgTrMole ) { rate = rh2h2p*hmi.htwo*h.hii + b2pcin*h.hii*phycon.eden*h.hi + h3ph2p*hmi.h3plus + h3petc*hmi.h3plus; if( rate > 1e-25 ) { fprintf( ioQQQ, " Create H2+, rate=%10.2e rh2h2p;%5.3f b2pcin;%5.3f h3ph2p;%5.3f h3petc+;%5.3f\n", rate, rh2h2p*h.hii*hmi.htwo/rate, b2pcin*h.hi*h.hii* phycon.eden/rate, h3ph2p*hmi.h3plus/rate, h3petc*hmi.h3plus/ rate ); } else { fprintf( ioQQQ, " Create H2+, rate=0\n" ); } } if( hmi.hminus > 0. && HMolLTE.phmlte > 0. ) { hmidepCom.hmidep = hmi.hminus/h.hi/phycon.eden/HMolLTE.phmlte; } else { hmidepCom.hmidep = 1.; } /* this will be net volume heating rate, photo heat - induc cool */ hmihetCom.hmihet = hmiht*hmi.hminus - colind; /* save it */ HeatingCom.heating[15][0] = hmihetCom.hmihet; HeatingCom.heating[16][0] = h2phet*hmi.h2plus; /* THTMOL is total heating due to absorption of Balmer continuum * thtmol = hmihet + h2phet * h2plus * departure coef for H2 defined rel to N(1s)**2 * (see Littes and Mihalas Solar Phys 93, 23) */ plte = (double)(h.hi) * (double)(h.hi) * HMolLTE.ph2lte; if( plte > 0. ) { hmidepCom.h2dep = hmi.htwo/plte; } else { hmidepCom.h2dep = 1.; } /* departure coef of H2+ defined rel to N(1s) N(p) * sec den was HI before 85.34 */ plte = h.hi*h.hii*HMolLTE.phplte; if( plte > 0. ) { hmidepCom.h2pdep = hmi.h2plus/plte; } else { hmidepCom.h2pdep = 1.; } /* departure coef of H3+ defined rel to N(H2+) N(p) */ if( ph3lte > 0. ) { hmidepCom.h3pdep = hmi.h3plus/ph3lte; } else { hmidepCom.h3pdep = 1.; } if( trace.lgTrace && trace.lgTrMole ) { fprintf( ioQQQ, " HMOLE, Dep Coef, H-:%10.2e H2:%10.2e H2+:%10.2e\n", hmidepCom.hmidep, hmidepCom.h2dep, hmidepCom.h2pdep ); fprintf( ioQQQ, " H- creat: Rad atch%10.3e Induc%10.3e bHneut%10.2e 3bod%10.2e b=>H2%10.2e N(H-);%10.2e b(H-);%10.2e\n", rspon, ratind, bhneut, c3bod, batach, hmi.hminus, hmidepCom.hmidep ); fprintf( ioQQQ, " H- destr: Photo;%10.3e mut neut%10.2e e- coll ion%10.2e =>H2%10.2e x-ray%10.2e p+H-%10.2e\n", gamhmi, faneut, cionhm, ratach, HPhoto.hgamnc[0][IP1S], fhneut ); fprintf( ioQQQ, " H- heating:%10.3e Ind cooling%10.2e Spon cooling%10.2e\n", hmihetCom.hmihet, colind, hmihetCom.hmicol ); } /* fmol is fraction of hydrgen in form of any molecule or ion */ fmol = (hmi.hminus + 2.*(hmi.htwo + hmi.h2plus))/phycon.hden; /* remember the largest fraction that occurs in the model */ h2max.BiggestH2 = (float)(MAX2(h2max.BiggestH2,fmol)); /*---------------------------------------------------------------- */ /* He H+ formation * rates taken from Flower+Roueff, Black * photodissociation through 1.6->2.3 continuum */ /* why is this in a look instead of GammaK? ] * to fix must set opacities into stack */ gamheh = 0.; for( i=hehpCom.iheh1-1; i < hehpCom.iheh2; i++ ) { gamheh += rfield.flux[i] + rfield.outcon[i] + rfield.outlin[i]; } gamheh *= 4e-18; /* both He+ + H and He + H+ * first rate is radiative associatino from * >>refer Zygelman, B., and Dalgarno, A. 1990, ApJ 365, 239 */ hehpCom.hehp = (float)(1e-15*h.hi*he.heii + 1e-20*h.hii*he.hei); /* H2+ + HE => HEH+ + H0 */ hehpCom.hehp += (float)(3e-10*exp(-6717./phycon.te)*he.hei*hmi.h2plus); /* hard radiation */ gamheh += 3.*HPhoto.hgamnc[0][IP1S]; /* HeH+ + H => H2+ + He */ gamheh += 1e-10*h.hi; /* recombination, HeH+ + e => He + H */ if( gamheh > hehpCom.hehp*1e-15 ) { gamheh += phycon.eden*1e-9; } else { gamheh = -1; } hehpCom.hehp /= (float)gamheh; # ifdef DEBUG_FUN fputs( " <->hmole()\n", debug_fp ); # endif return; }