/*hmole determine populations of hydrogen molecules */ /*hmirat compute radiative association rate for H- */ #include "cddefines.h" #ifdef NMOLE # undef NMOLE #endif #define NMOLE 6 #include "physconst.h" #include "iso.h" #include "nhe1lvl.h" #include "grainvar.h" #include "abundances.h" #include "secondaries.h" #include "colden.h" #include "opacity.h" #include "heavy.h" #include "ionfracs.h" #include "rfield.h" #include "converge.h" #include "hmi.h" #include "recom.h" #include "heat.h" #include "punch.h" #include "timesc.h" #include "trace.h" #include "nomole.h" #include "phycon.h" #include "iphmin.h" #include "ghabng.h" #include "doppvel.h" #include "nhe1.h" #include "typmatrx.h" #include "linpack.h" #include "rtescprob.h" #include "gammas.h" #include "veclib.h" #include "hmrate.h" #include "hmole.h" /*hmirat compute radiative association rate for H- */ static double hmirat(double te); 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, ih2_star, ih2p, ih3p, ihmi, iho, ipConServ, ipiv[NMOLE], ipLo, j, job, limit , merror, nd; int lgNegPop; double batach, bh2dis, bh2h2p, bhneut, c3bod, cionhm, corr, damper, deexc, desh2p, eh2hh, ex3hp, excit, exph2, exphmi, exphp, fa, fac, faneut, fhneut, fmol, gamhd, gamheh, h1fnd, h1good, h2pcin, h2ph3p, h2phhp, h2pion, h32h2, h3petc, h3ph2p, hnnew, ph3lte, radasc, radath, ratach, rate, renorm, rh2h2p, saha; /* lte populations of H-, H2, and H2+ */ static double phmlte, phplte; static double HMinus_induc_rec_cooling, eh2hhm, eh2old, gamtwo, gh2dis, gh2exc_dissoc, h2esc, h2phet, HMinus_photo_heat, hneut, HMinus_induc_rec_rate, rgrain, rspon, th2; static long int nzoneEval=9999; double amat[NMOLE][NMOLE], b2pcin, bvec[NMOLE], c[NMOLE+1][NMOLE+1], plte, rcond, eh3_h2h , work[NMOLE]; # ifdef DEBUG_FUN fputs( "<+>hmole()\n", debug_fp ); # endif /* CONST is 8\pi * nu(912)**3 / c**2 */ # define CONST_ 9.940622e26 /* 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.; hmi.htwo = 0.; /* this is where the Emline struc expects to find the abundance */ xIonFracs[LIMELM+2][0] = 0.; hmi.h2plus = 0.; hmi.hmihet = 0.; hmi.H2Opacity = 0.; hmi.htwo_star = 0.; hmi.hmicol = 0.; hmi.HeatH2Dish = 0.; hmi.HeatH2Dexc = 0.; hmi.hmidep = 1.; hmi.hehp = 0.; hmi.rh2dis = 0.; hmi.HalphaHmin = 0.; heat.heating[0][15] = 0.; heat.heating[0][16] = 0.; heat.heating[0][17] = 0.; recom.GrnIonRec = 0.; # ifdef DEBUG_FUN fputs( " <->hmole()\n", debug_fp ); # endif return; } /* >>chng 02 may 27, add option to continue in full H2 limit */ if( hmi.lgAll_H2 ) { hmi.hminus = 0.; hmi.h2plus = 0.; hmi.h3plus = 0.; hmi.htwo = (float)(phycon.hden/2.); /* NB this must be kept parallel with nelem and ionstag in H2Lines EmLine struc, * since that struc expects to find the abundances here */ xIonFracs[LIMELM+2][0] = hmi.htwo; hmi.htwo_star = 0.; hmi.hmihet = 0.; hmi.H2Opacity = 0.; hmi.hmicol = 0.; hmi.HeatH2Dish = 0.; hmi.HeatH2Dexc = 0.; hmi.hmidep = 1.; hmi.hehp = 0.; hmi.HalphaHmin = 0.; heat.heating[0][15] = 0.; heat.heating[0][16] = 0.; heat.heating[0][17] = 0.; recom.GrnIonRec = 0.; # ifdef DEBUG_FUN fputs( " <->hmole()\n", debug_fp ); # endif return; } /* this checks whether very first call in this model */ if( !conv.nTotalIoniz ) { phmlte = 0.; hmi.ph2lte = 0.; phplte = 0.; } /* 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 */ saha = sqrt(SAHA2); exphmi = sexp(8.745e3/phycon.te); if( exphmi > 0. ) { /* these are ratio n*(H-)/[ n*(ne) n*(Ho) ] */ phmlte = saha/(phycon.te32*exphmi)*(1./(2.*2.)); } else { phmlte = 0.; } /* population of H2 in LTE * dissociation energy is 4.477eV */ exph2 = sexp(5.195e4/phycon.te); if( exph2 > 0. ) { if( nzone <= 0 ) { damper = 0.; } else { damper = 0.5; } /* extra factor accounts for mass of H instead of e- in SAHA * last factor was put in ver 85.23, missing before */ hmi.ph2lte = (1. - damper)*saha/(phycon.te32*exph2)* (1./(2.*2.))*3.634e-5 + damper*hmi.ph2lte; } else { hmi.ph2lte = 0.; } { /*@-redef@*/ /* often the H- route is the most efficient formation mechanism for H2, * will be through rate called ratach * this debug print statement is to trace h2 oscillations */ enum {DEBUG=FALSE}; /*@+redef@*/ if( DEBUG && nzone>187&& iteration > 1/**/) { /* rapid increase in H2 density caused by rapid increase in hmi.ph2lte */ fprintf(ioQQQ,"ph2lteee\t%.2e\t%.1e\t%.1e\n", hmi.ph2lte, exph2, phycon.te); } } /* 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 */ phplte = saha/(phycon.te32*exphp)*(4./(2.*1.))*3.634e-5; } else { 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/(phycon.te32*ex3hp)*(4./(2.*1.))*3.634e-5; } else { ph3lte = 0.; } rspon = hmirat(phycon.te); hmi.hmicol = rspon*EN1RYD*phycon.te*1.15e-5; /*fprintf(ioQQQ,"%.2e %.2e %.2e %.2e\n", phycon.te, rspon , hmi.hmicol, hmi.hmicol/(rspon*EN1RYD*phycon.te*1.15e-5) );*/ /* get per unit vol */ rspon *= phycon.eden; hmi.hmicol *= phycon.eden*xIonFracs[ipHYDROGEN][0]; /* ================================================================= */ /* evaluate H- photodissociation rate, induc rec and rec cooling rates */ /* >>chng 00 dec 24, add test so that photo rate only reevaluated two times per zone. * in grain-free models this was sometimes dominated by Lya and so oscillated. * especially bad in primal.in - change 2 to 4 and primal.in will stop due to Lya oscil */ if( /**/(conv.nPres2Ioniz < 2) || (nzone==0) || (nzoneEval!=nzone) ) { /* >>chng 02 feb 16, add damper on H- photo rate, wild oscillations in Lya photo rate in * grain free models */ static double hm_damper=0.25; static double phot_new=0. , phot_old; nzoneEval = nzone; /*hmi.HMinus_photo_rate = GammaBn( iphminCom.iphmin-1 , nhe1Com.nhe1[0] , opac.iphmop , 0.055502 , &HMinus_induc_rec_rate , &HMinus_induc_rec_cooling );*/ phot_old = phot_new; phot_new = GammaBn( iphminCom.iphmin-1 , nhe1Com.nhe1[0] , opac.iphmop , 0.055502 , &HMinus_induc_rec_rate , &HMinus_induc_rec_cooling ); if( nzone==0 ) { hmi.HMinus_photo_rate = phot_new; } else { hmi.HMinus_photo_rate = hm_damper*phot_new+ (1.-hm_damper)*phot_old; } { /* following should be set true to print populations */ /*@-redef@*/ enum {DEBUG=FALSE}; /*@+redef@*/ if( DEBUG) { fprintf(ioQQQ,"hminphoto\t%li\t%li\t%.2e\n", nzone, conv.nPres2Ioniz , hmi.HMinus_photo_rate ); } } /* save H- photodissociation heating */ HMinus_photo_heat = heat.HeatNet; /* induced recombination */ HMinus_induc_rec_rate *= phmlte*phycon.eden; /* induced recombination cooling per unit volume */ HMinus_induc_rec_cooling *= phmlte*phycon.eden*xIonFracs[ipHYDROGEN][0]; { /* following should be set true to debug H- photoionization rates */ /*@-redef@*/ enum {DEBUG=FALSE}; /*@+redef@*/ if( DEBUG && nzone>187&& iteration > 1) { fprintf(ioQQQ,"hmoledebugg %li ",nzone); GammaPrt( iphminCom.iphmin-1 , nhe1Com.nhe1[0] , opac.iphmop , /* io unit we will write to */ ioQQQ, /* total photo rate from previous call to GammaK */ hmi.HMinus_photo_rate, /* we will print contributors that are more than this rate */ hmi.HMinus_photo_rate*0.05); } } } /* add on high energy ionization, assume hydrogen cross section * n.b.; HGAMNC with secondaries */ /* >>chng 00 dec 24, above goes to hei edge, no need for this, and was not important*/ /*hmi.HMinus_photo_rate += iso.gamnc[ipH_LIKE][0][ipH1s];*/ /* ================================================================= */ /* photodissociation by Lyman band absorption: esc prob treatment, * treatment based on * >>refer HI abs Tielens and Hollenbach 1985 ApJ 291, 722. */ /* do up to carbon photo edge if carbon is turned on */ /* >>>chng 00 apr 07, add test for whether element is turned on */ if( abundances.lgElmtOn[5] ) { /* carbon is turned on, use carbon 1 edge */ ipLo = Heavy.ipHeavy[ipCARBON][0] - 1; } else { /* carbon truned off, use hydrogen balmer continuum */ ipLo = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2s]-1; } ghabng.GammaHabing = 0.; /* >>chng 00 apr 07 from explicit ipHeavy to ipLo */ /* find total intensity over carbon-ionizing continuum */ for( i=ipLo; i < iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]; i++ ) { ghabng.GammaHabing += (rfield.flux[i] + rfield.ConInterOut[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); /* escape prob is line shielding, eqn A10 of TH85*/ hmi.H2Opacity = (float)(1.2e-14*(1e5/DoppVel.doppler[0])); th2 = coldenCom.colden[ipCOLH2]*1.2e-14*(1e5/DoppVel.doppler[0]); h2esc = esc_PRD_1side(th2,1e-4); /* cross section is eqn A8 of Teilens and Hollenbach 85a * branching ratio of 10% in, so like their table 5 */ /* the total H2 pump rate, eqn A8 of TH85, * 10% lead to dissociation through H_2 + h nu => 2H */ gh2dis = ghabng.GammaHabing*3.4e-11 * h2esc; /* collisional ionization of H-, rate from Janev, Langer et al. */ if( phycon.te < 3074. ) { cionhm = 1.46e-32*(powi(phycon.te,6))*phycon.sqrte*exphmi; } else if( phycon.te >= 3074. && phycon.te < 30000. ) { cionhm = 5.9e-19*phycon.te*phycon.te*phycon.sqrte*phycon.te03* phycon.te01*phycon.te01; } else { cionhm = 3e-7; } /* ion recomb on grain surfaces */ recom.GrnIonRec = 0.; /* H2 formation on grains; * rate from * >>refer H2 grain formation Hollenback, D., & McKee, C.F., ApJS, 41, 555 eq 3.4 3.8 */ if( gv.lgDustOn ) { double RateCoef; RateCoef = 0.; recom.GrnIonRec = 0.; for( nd=0; nd < gv.nBin; nd++ ) { /* >>chng 02 feb 15, removed check tedust > 1.01, change in InitGrains * guarantees that all relevant parameters are initialized, PvH */ /* sticking probability, equation 3.7 of * >>refer Hollenbach, D.J., & McKee, C.F., 1979, ApJS, 41, 555, * fraction of H impacts on grain surface that stick */ double s = 1./(1. + 0.04*sqrt(gv.bin[nd]->tedust+phycon.te) + 0.002*phycon.te + 8e-6*phycon.te*phycon.te); /* fraction of impacts that produce H2 before evaporation from grain surface. * this is equation 3.4 of * >>refer Hollenbach, D.J., & McKee, C.F., 1979, ApJS, 41, 555 * 1e4 is ratio of total absorption sites to apropriate sites * 920 is D_H and chosen to get f_a = 0.5 at 100 K. * factor of 0.6252 needed to obtain std ism rate to be 3e-17 at 100 K, * the value deduced by * >>refer H2 grain form Jura, M., ApJ, 197, 581 */ fa = 0.6252/(1. + 1e4*sexp(920./gv.bin[nd]->tedust)); /* NB IntArea is total, not projected, area, must div by 4 */ RateCoef += 0.5 * DoppVel.AveVel[ipHYDROGEN]* gv.bin[nd]->IntArea/4. * gv.bin[nd]->cnv_H_pCM3 * s * fa ; ASSERT( RateCoef > 0. ); /* rate ions recombine on grain surface - used elsewhere * based on * >>refer grain rec Draine 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 */ recom.GrnIonRec +=(float)(5.8e-13/phycon.sqrte*gv.bin[nd]->cnv_H_pCM3*recom.grecon); ASSERT( recom.GrnIonRec > 0. ); } /* >>chng 02 jan 05, use resolved grains for H_2 formation */ rgrain = RateCoef; /*fprintf(ioQQQ," H2 grain form rate %.2e\n", rgrain );*/ } else { rgrain = 0.; recom.GrnIonRec = 0.; } /* collisional dissociation, rate from * >>refer H2 coll disc 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-phycon.alogte*3.08); corr = pow(10.,MAX2(0.,corr)/(1. + 1.6e4/xIonFracs[ipHYDROGEN][0])); hmi.rh2dis = (float)(1.55e-8/phycon.sqrte*sexp(65107./phycon.te)* corr); /* H2 + H+ => H2+ + H *>>chng 98 jan 02, from 2.12e4 to 2.123e4 */ rh2h2p = 1.8e-12*phycon.sqrte*phycon.te10/phycon.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*phycon.sqrte*phycon.te10/phycon.te01*2./16.; /* H2+ + HNU => H+ + H */ gamtwo = GammaK(opac.ih2pnt[0],opac.ih2pnt[1],opac.ih2pof,1.); /*GammaPrt(opac.ih2pnt[0],opac.ih2pnt[1],opac.ih2pof,ioQQQ,gamtwo,0.01);*/ h2phet = heat.HeatNet; /* set up matrix that will be used for the molecule formation rates */ for( i=0; i < NMOLE; i++ ) { c[NMOLE][i] = 0.; for( j=0; j < NMOLE; j++ ) { c[j][i] = 0.; } } hneut = xIonFracs[ipHYDROGEN][0] + hmi.hminus + 2.*(hmi.htwo + hmi.h2plus) + 3.*hmi.h3plus; /* these should be totally equivalent?? * >>chng 96 jun 3, hneut floated above hden, ended negative */ hneut = phycon.hden - xIonFracs[ipHYDROGEN][1]; /* protect against roundoff in ionized gas */ if( hneut/phycon.hden < 1e-3 ) { hneut = xIonFracs[ipHYDROGEN][0] + hmi.hminus + 2.*(hmi.htwo + hmi.htwo_star + hmi.h2plus) + 3.* hmi.h3plus; } iho = 0; strcpy( chLab[iho], "HO " ); ihmi = 1; strcpy( chLab[ihmi], "H- " ); ih2 = 2; strcpy( chLab[ih2], "H2 " ); ih2p = 3; strcpy( chLab[ih2p], "H2+ " ); ih3p = 4; strcpy( chLab[ih3p], "H3+ " ); ih2_star = 5; strcpy( chLab[ih2_star], "H2* " ); /*-------------------------------------------------------------------- */ /* H- H minus hminus balance equations * (IHMI,IHO) == processes making H- from Ho =+sign * radiative attachment: HI + NE => H-; */ c[iho][ihmi] += rspon + HMinus_induc_rec_rate; c[iho][iho] += -rspon - HMinus_induc_rec_rate; /* (IHMI,IHMI) = processes destroying H- =-sign * photodissociation, H- + H NU => H + NE */ c[ihmi][ihmi] -= hmi.HMinus_photo_rate; c[ihmi][iho] += hmi.HMinus_photo_rate; /* mutual neutralization with heavies, rate from Dalgarno and Mcray * all charged ions contribute equally */ faneut = 4e-6/phycon.sqrte*MAX2(0.,phycon.eden-xIonFracs[ipHYDROGEN][1]-xIonFracs[ipHELIUM][1]- 2.*xIonFracs[ipHELIUM][2]); c[ihmi][ihmi] -= faneut; /* electron collisional ionization of H- */ cionhm *= phycon.eden; c[ihmi][ihmi] += -cionhm; c[ihmi][iho] += cionhm; /* inverse process; three body rec */ c3bod = cionhm*(phmlte*phycon.eden); c[iho][ihmi] += c3bod; c[iho][iho] -= c3bod; /* associative detachment: H- + H => H2 + E */ ratach = xIonFracs[ipHYDROGEN][0]*1.35e-9; c[ihmi][ihmi] -= ratach; c[iho][iho] -= hmi.hminus*1.35e-9; /* (1,IH2) convert H2 into H- = +sign * the back reaction, H2 + e => H- + Ho */ if( hmi.ph2lte > 0. ) { batach = (ratach/xIonFracs[ipHYDROGEN][0])*phmlte/hmi.ph2lte*phycon.eden; } else { batach = 0.; } { /*@-redef@*/ /* often the H- route is the most efficient formation mechanism for H2, * will be through rate called ratach * this debug print statement is to trace h2 oscillations */ enum {DEBUG=FALSE}; /*@+redef@*/ if( DEBUG && nzone>187&& iteration > 1/**/) { /* rapid increase in H2 density caused by rapid increase in hmi.ph2lte */ fprintf(ioQQQ,"batach\t%.2e\t%.1e\t%.1e\t%.1e\t%.1e\t%.1e\t%.1e\n", batach, ratach, xIonFracs[ipHYDROGEN][0], phmlte, hmi.ph2lte, phycon.eden, hmi.hminus ); } } c[ih2][ihmi] += batach; c[ih2][iho] += batach; /* mutual neut, mostly into n=3; rates from Janev et al * H- + H+ => H + H(n=3) */ fhneut = xIonFracs[ipHYDROGEN][1]*7e-8; c[ihmi][ihmi] -= fhneut; c[ihmi][iho] += fhneut; /* back reaction from excited state H */ if( phycon.te > 1000. ) { /* 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( nzone <= 1 ) { /* this is initial setup of code, so set rate coef to actual val */ eh2hhm = 2.7e-8*pow(10.,-0.7*POW2(phycon.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(phycon.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][ihmi] += eh2hhm; c[ih2][iho] += eh2hhm; /*-------------------------------------------------------------------- * * molecular hydrogen H2 htwo balance equation * (IH2,IHO)==create H2 from Ho =+ * * H2 formation on grains; * rate coefficient from * >>refer H2 rates Hollenback, D.J., and McKee, C.F., 1979, ApJS, 41, 555 eq 3.4 3.8 * really is rate coefficient, needs two HIs */ /* >>chng 01 jan 05, remove from matrix part and add hden to rgrain, */ c[iho][ih2] += rgrain/* *xIonFracs[ipHYDROGEN][0]*/; c[iho][iho] += -rgrain/* *xIonFracs[ipHYDROGEN][0]*/; /* excited atom radiative association, * H(n=2) + H(n=1) => H2 + hnu * >>refer H2 rates Latter, W.B., & Black, J.H., 1991, Ap.J. 372, 161 */ radasc = ((iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2p] + iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2s])*xIonFracs[ipHYDROGEN][1])*3e-14; c[iho][ih2] += radasc; c[iho][iho] -= radasc; /* make H2 from H- =+ sign * associative detachment; H- + H => H2: * >>refer H2 rates Browne & Dalgarno J PHys B 2, 885 */ c[ihmi][ih2] += ratach; /* photo-destroy H2 */ /* >>chng 00 nov 25 factor of 0.1, assume pump is total, and 10% distroy H2 */ c[ih2][ih2] -= gh2dis*0.1; c[ih2][iho] += gh2dis*0.1; /* the process H2 + e- => H + H + e * >>refer H2 rates 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][ih2] -= eh2hh; c[ih2][iho] += eh2hh; /* the processes H2(v>=4) + e- => H + H- * >>refer H2 rates Lenzuni et al. apj sup 76, 759 * >>refer H2 rates Janev et al. * density dep is for non-lte, guess from dalgarno and roberge apl 233, 25 */ c[ih2][ih2] -= eh2hhm; /* add on inverse of radiative attachment */ c[ih2][ih2] -= batach; /* collisional dissociation, rate from * >>refer H2 rates Dove and Mandy, Ap.J.(Let) 311, L93. * H_2 + S => 2H + S */ c[ih2][ih2] += -hmi.rh2dis*xIonFracs[ipHYDROGEN][0]; c[ih2][iho] += hmi.rh2dis*xIonFracs[ipHYDROGEN][0]; /* back rate, three body recombination, 2H + S => H_2 + S */ bh2dis = hmi.rh2dis*hmi.ph2lte*xIonFracs[ipHYDROGEN][0]*xIonFracs[ipHYDROGEN][0]; c[iho][ih2] += bh2dis; c[iho][iho] -= bh2dis; /* H2 + HNU=> H2+ + E * photoionization by hard photons, crossection=3*HI */ gamhd = iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]; c[ih2][ih2] -= gamhd; /* assume cosmic rays do the same thing */ /* >>chng 00 nov 28, factor of 0.93 from >>refer cosmic ray ionization rate Tielens, A.G.G.M., & Hollenbach, D., 1998, ApJ, 291, 722 */ c[ih2][ih2] -= Secondaries.csupra*0.93; /* H2 + H+ => H2+ + H */ c[ih2][ih2] += -rh2h2p*xIonFracs[ipHYDROGEN][1]; c[ih2][iho] += rh2h2p*xIonFracs[ipHYDROGEN][1]; /* H2+ + H => H+ + H2 */ c[ih2p][ih2] += bh2h2p*xIonFracs[ipHYDROGEN][0]; c[iho][iho] += -bh2h2p*hmi.h2plus; /*TODO this rate drives numerical instability in such models as secondary1 and 2.in */ /*h3ph2p = hmrate(2.08e-9,0.,1.88e4)*xIonFracs[ipHYDROGEN][0];*/ /* this rate couples H2+ and H3+, and tends to destabalize the matrix in both highly * ionized and fully molecular conditions. Setting this to zero had no effect - the th85 * predictions were identical. * The long range solution is to break the H2+ and H3+ parts of the network out of the * main H2 solver. For now it is on the medium range to do list since it had no effects * whatsoever on pdr calculations. */ h3ph2p = 0.; c[ih3p][ih2] += h3ph2p; c[iho][iho] += -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][ih2] += h3petc; c[ih2][iho] += h3petc; /* H2 + H3+ => H2 + H+ + H2 */ h32h2 = hmrate(3.41e-11,0.5,5.04e4)*hmi.htwo; c[ih3p][ih2] += 2.*h32h2; /* e + H3+ => H2 + H */ eh3_h2h = hmrate(5.00e-9,-0.5,0.)*phycon.eden; c[ih3p][ih2] += hmrate(5.00e-9,-0.5,0.)*phycon.eden; c[ih3p][iho] += hmrate(5.00e-9,-0.5,0.)*phycon.eden; if( (trace.lgTrace && trace.lgTrMole) || punch.lgPunH2 ) { if( punch.lgPunH2 ) { ioFile = punch.ipPunH2; } else { ioFile = ioQQQ; } if( c[ih2][ih2] != 0. ) { rate = -c[ih2][ih2]; fprintf( ioFile, " Destroy H2: rate=%10.2e DIS;%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*0.1/rate, batach/rate, hmi.rh2dis*xIonFracs[ipHYDROGEN][0]/ rate, gamhd/rate, rh2h2p*xIonFracs[ipHYDROGEN][1]/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][ih2p] += radath*xIonFracs[ipHYDROGEN][1]; c[iho][iho] += -radath*xIonFracs[ipHYDROGEN][1]; /* make H2+ from H2 =+sign * H2 + HNU => H2+ + E * also cosmic rays */ /* >>chng 00 nov 28, factor of 0.93 from >>refer cosmic ray ionization rate Tielens, A.G.G.M., & Hollenbach, D., 1998, ApJ, 291, 722 */ c[ih2][ih2p] += gamhd + Secondaries.csupra*0.93; /* H2 + H+ => H2+ + H */ c[ih2][ih2p] += rh2h2p*xIonFracs[ipHYDROGEN][1]; /* (3,IH2P) == destroy H2+ = -sign * H + H2+ => H+ + H2 */ c[ih2p][ih2p] += -bh2h2p*xIonFracs[ipHYDROGEN][0]; c[iho][iho] += -bh2h2p*hmi.h2plus; /* H2+ + p => H + H+ + H+; Janev et al. 3.2.6 */ h2pion = 2.4e-27*POW3(phycon.te)*xIonFracs[ipHYDROGEN][1]; c[ih2p][ih2p] -= h2pion; c[ih2p][iho] += h2pion; /* H2+ + E => H + H+ + e-; Janev et al. */ h2pcin = 2e-7*sexp(30720./phycon.te); c[ih2p][ih2p] += -h2pcin*phycon.eden; c[ih2p][iho] += h2pcin*phycon.eden; /* back reaction, H + H+ + e => h2+ + e */ b2pcin = h2pcin*phplte; /* this is the hot reaction at high densities */ c[iho][ih2p] += b2pcin*phycon.eden*xIonFracs[ipHYDROGEN][1]; c[iho][iho] += -b2pcin*phycon.eden*xIonFracs[ipHYDROGEN][1]; /* H2+ + HNU => H+ + H */ c[ih2p][ih2p] -= gamtwo; /* photoionization by hard photons, crossection =2*HI (wild guess) */ c[ih2p][ih2p] += -2.*iso.gamnc[ipH_LIKE][0][ipH1s]; /* H2 + H2+ => H + H3+ */ h2ph3p = 1.40e-9*hmi.htwo*(1. - sexp(9940./phycon.te)); c[ih2p][ih2p] -= h2ph3p; c[ih2p][iho] += h2ph3p; /* h + h3+ => h2 + h2+ */ c[ih3p][ih2p] += h3ph2p; /* H2 + H3+ => H2 + H2+ + H */ c[ih3p][ih2p] += h3petc; /* destroy H2+ via H2+ + H2 => H + H+ + H2 */ h2phhp = 2.41e-12*phycon.sqrte*sexp(30720./phycon.te)*hmi.htwo; c[ih2p][ih2p] -= h2phhp; c[ih2p][iho] += h2phhp; /* save total H2P destruction rate for possible later printout: * NB this must come last */ desh2p = -c[ih2p][ih2p]; /*------------------------------------------------------------------ */ /* H3+ balance equations*/ /* H2 + H2+ => H + H3+ */ c[ih2p][ih3p] += h2ph3p; /* H + H3+ => H2 + H2+ */ c[ih3p][ih3p] -= h3ph2p; /* H2 + H3+ => H2 + H2+ + H */ c[ih3p][ih3p] -= h3petc; /* H2 + H3+ => H2 + H+ + H2 */ c[ih3p][ih3p] -= h32h2; /* e + H3+ => 3H * e + H3+ => H2 + H */ c[ih3p][ih3p] -= 2.*hmrate(5.00e-9,-0.5,0.)*phycon.eden; /* photoionization by hard photons, crossection =2*HI (wild guess) */ c[ih3p][ih3p] -= 2.*iso.gamnc[ipH_LIKE][0][ipH1s]; /*------------------------------------------------------------------ */ /* vib excited H2, called H2* balance equations, these closely follow * >>refer mh2 fits Tielens, A.G.G.M., & Hollenbach, D., 1985a, ApJ 291, 722 */ /* population of vib-excited H2, from discussion on pp 736-737 of TH85 */ /* deexcitation rate from upper level, H2* => H2 */ deexc = hmi.htwo*1.4e-12*phycon.sqrte * sexp( 18100./(phycon.te + 1200.) ) + xIonFracs[ipHYDROGEN][0] * 1e-12*phycon.sqrte * sexp(1000./phycon.te ); /* depopulate H2_star, 2e-7 is spontaneous deexcitation rate, * which also appears in lines where intensity of vib lines is entered into line stack */ c[ih2_star][ih2_star] -= (2e-7+ deexc); /* H2_star making producing H2 */ c[ih2_star][ih2] += (2e-7+ deexc); /* collisional excitation of vib from ground, * stat weight of ground 1, excit 6, as per TH discussion * this must normally be zero */ excit = 6. * deexc * sexp( 30172./phycon.te); c[ih2][ih2] -= excit; /* H2 producing H2_star */ c[ih2][ih2_star] += excit; /* assume that 0.9 of H2 dissociations lead to H2_star, * H2 + 0.9*gh2dis => h2_star */ c[ih2][ih2] -= 0.9*gh2dis ; c[ih2][ih2_star] += 0.9*gh2dis ; /* rate of photodissoc of vib-excit H2, A12 of TH85 */ gh2exc_dissoc = ghabng.GammaHabing*1e-11; c[ih2_star][ih2_star] -= gh2exc_dissoc ; /* it goes to neutral H */ c[ih2_star][iho] += gh2exc_dissoc; { /* following should be set true to print populations */ /*@-redef@*/ enum {DEBUG=FALSE}; /*@+redef@*/ if( DEBUG ) { /* these are the raw results */ fprintf( ioQQQ, " HMOLE h2 %.2e h2* %.2e\n" , hmi.htwo ,hmi.htwo_star ); } } /*------------------------------------------------------------------ */ /* 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 = 0; c[iho][ipConServ] = 1.; c[ihmi][ipConServ] = 1.; c[ih2][ipConServ] = 2.; c[ih2p][ipConServ] = 2.; c[ih3p][ipConServ] = 3.; c[ih2_star][ipConServ] = 2.; /* the above will add up to the following */ c[NMOLE][ipConServ] = hneut; /*------------------------------------------------------------------ */ if( trace.lgTrace && trace.lgTrMole ) { /* print the full matrix */ fprintf( ioQQQ, " "); for( i=0; i < NMOLE; i++ ) { fprintf( ioQQQ, " %s", chLab[i] ); } fprintf( ioQQQ, " \n" ); /* [0][0] [0][1] [0][2] [0][3] [0][4] [0][5] [1][0] [1][1] [1][2] [1][3] [1][4] [1][5] [2][0] [2][1] [2][2] [2][3] [2][4] [2][5] [3][0] [3][1] [3][2] [3][3] [3][4] [3][5] [4][0] [4][1] [4][2] [4][3] [4][4] [4][5] [5][0] [5][1] [5][2] [5][3] [5][4] [5][5] [iho][iho] [iho][ihmi] [iho][ih2] [iho][ih2p] [iho][ih3p] [iho][ih2_star] [ihmi][iho] [ihmi][ihmi] [ihmi][ih2] [ihmi][ih2p] [ihmi][ih3p] [ihmi][ih2_star] [ih2][iho] [ih2][ihmi] [ih2][ih2] [ih2][ih2p] [ih2][ih3p] [ih2][ih2_star] [ih2p][iho] [ih2p][ihmi] [ih2p][ih2] [ih2p][ih2p] [ih2p][ih3p] [ih2p][ih2_star] [ih3p][iho] [ih3p][ihmi] [ih3p][ih2] [ih3p][ih2p] [ih3p][ih3p] [ih3p][ih2_star] [ih2_star][iho] [ih2_star][ihmi] [ih2_star][ih2] [ih2_star][ih2p] [ih2_star][ih3p] [ih2_star][ih2_star] */ for( i=0; i < NMOLE; i++ ) { fprintf( ioQQQ, " MOLE%2ld %s", i ,chLab[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][ih2] > SMALLFLOAT ) { timesc.AgeH2MoleDest = -1./c[ih2][ih2]; } 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]" ); cdEXIT(EXIT_FAILURE); } } 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]" ); cdEXIT(EXIT_FAILURE); } /* 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]" ); cdEXIT(EXIT_FAILURE); } else { fprintf( ioQQQ, " chMatrix type insane in hydrogen, was%7.7s\n", TypMatrx.chMatrix ); puts( "[Stop in hmole]" ); cdEXIT(EXIT_FAILURE); } /* put derived abundances back into appropriate molecular species */ h1fnd = c[NMOLE][iho]; hmi.hminus = (float)c[NMOLE][ihmi]; hmi.h2plus = (float)c[NMOLE][ih2p]; hmi.h3plus = (float)c[NMOLE][ih3p]; hmi.htwo = (float)c[NMOLE][ih2]; /* NB this must be kept parallel with nelem and ionstag in H2Lines EmLine struc, * since that struc expects to find the abundances here */ xIonFracs[LIMELM+2][0] = hmi.htwo; hmi.htwo_star = (float)c[NMOLE][ih2_star]; lgNegPop = FALSE; for( i=0; i < NMOLE; i++ ) { if( c[NMOLE][i] < 0. ) lgNegPop = TRUE; } /* >>>chng 00 apr 04 add test for neg pops - had not been included before!! */ if( lgNegPop ) { enum {FIXUP=TRUE}; static long int nFixup=0; if( hmi.htwo/phycon.hden > 0.49 ) { /* we hit fully molecular limit - assume all H is in H2 */ fprintf( ioQQQ, " HMOLE finds fully H2 molecular limt in zone %li - H molecules network turned off.\n",nzone ); /* this is flag saying we hit this limit, in future will not attempt solution */ hmi.lgAll_H2 = TRUE; hmi.hminus = 0.; hmi.h2plus = 0.; hmi.h3plus = 0.; hmi.htwo = (float)(phycon.hden/2.); /* NB this must be kept parallel with nelem and ionstag in H2Lines EmLine struc, * since that struc expects to find the abundances here */ xIonFracs[LIMELM+2][0] = hmi.htwo; hmi.htwo_star = 0.; } else { fprintf( ioQQQ, " HMOLE finds negative H molecule, in zone %li.\n",nzone ); fprintf( ioQQQ, " The populations follow:"); for( i=0; i < NMOLE; i++ ) { fprintf(ioQQQ," %s %.2e", chLab[i], c[NMOLE][i]); } fprintf( ioQQQ, "\n" ); if( nzone == 0 && FIXUP && nFixup < 10 ) { for( i=0; i < NMOLE; i++ ) { c[NMOLE][i] = MAX2(0., c[NMOLE][i]); } ++nFixup; fprintf(ioQQQ," FIXUP taken %li time.\n\n", nFixup); } else { ShowMe(); puts( "[Stop in hmole]" ); cdEXIT(EXIT_FAILURE); } } } { /* following should be set true to print populations */ /*@-redef@*/ enum {DEBUG=FALSE}; /*@+redef@*/ if( DEBUG && (nzone>187) && (iteration > 1) ) { double create_from_Ho,create_3body_Ho,create_batach,create_eh2hhm,destroy_photo, destroy_coll_heavies,destroy_coll_electrons,destroy_Hattach,destroy_fhneut, destsum , createsum; create_from_Ho = (rspon + HMinus_induc_rec_rate); create_3body_Ho = c3bod; create_batach = batach; create_eh2hhm = eh2hhm; destroy_photo = hmi.HMinus_photo_rate; destroy_coll_heavies = faneut; destroy_coll_electrons = cionhm; destroy_Hattach = ratach; destroy_fhneut = fhneut; destsum = destroy_photo + destroy_coll_heavies + destroy_coll_electrons + destroy_Hattach + destroy_fhneut; fprintf(ioQQQ,"hminusdebugg\t%.2e\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f", destsum, destroy_photo/destsum , destroy_coll_heavies/destsum, destroy_coll_electrons/destsum, destroy_Hattach/destsum, destroy_fhneut/destsum ); fprintf(ioQQQ,"\n"); createsum = create_from_Ho+create_3body_Ho+create_batach+create_eh2hhm; fprintf(ioQQQ,"createsum\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", createsum, create_from_Ho/createsum, create_3body_Ho/createsum, create_batach/createsum, create_eh2hhm/createsum); } } /* rate H-alpha is created by H- ct */ hmi.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.*iso.gamnc[ipH_LIKE][0][ipH1s]/ rate, h2ph3p/rate, h2phhp/rate, h2pion/rate, bh2h2p* xIonFracs[ipHYDROGEN][0]/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*xIonFracs[ipHYDROGEN][1]*xIonFracs[ipHYDROGEN][0]/rate, (gamhd + Secondaries.csupra*0.93)*hmi.htwo/rate, rh2h2p*xIonFracs[ipHYDROGEN][1]*hmi.htwo/rate, b2pcin*xIonFracs[ipHYDROGEN][0]*xIonFracs[ipHYDROGEN][1]*phycon.eden/rate ); } else { fprintf( ioQQQ, " Create H2+: rate= is zero\n" ); } } } /* heating due to H2 dissociation */ if( nomole.lgNoH2Mole ) { hmi.HeatH2Dish = 0.; heat.heating[0][17] = 0.; hmi.HeatH2Dexc = 0.; heat.heating[0][8] = 0.; } else { /* H2 photodissociation heating, eqn A9 of Tielens & Hollenbach 1985a */ hmi.HeatH2Dish = (float)(1.36e-23*hmi.htwo*h2esc*ghabng.GammaHabing); heat.heating[0][17] = hmi.HeatH2Dish; /* >>chng 00 nov 25, explictly break out this heat agent */ /* 2.6 eV of heat per deexcitation, consider difference * between deexcitation (heat) and excitation (cooling) */ hmi.HeatH2Dexc = (float)( (hmi.htwo_star * deexc -hmi.htwo*excit) * 4.17e-12); heat.heating[0][8] = hmi.HeatH2Dexc; } *hmovh1 = (float)(((hmi.htwo + hmi.h2plus)*2. + 3.*hmi.h3plus + hmi.hminus) / xIonFracs[ipHYDROGEN][0]); 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, xIonFracs[ipHYDROGEN][0], phycon.eden ); fprintf( ioQQQ, " HMOLE: h1fnd=%10.2e\n", h1fnd ); ShowMe(); puts( "[Stop in hmole]" ); cdEXIT(EXIT_FAILURE); } h1good = phycon.hden/(1.e0 + (double)(*hmovh1) + h2ovh1); hnnew = h1good*(1.e0 + (double)(*hmovh1)); renorm = hnnew/hneut; { /* following should be set true to print populations */ /*@-redef@*/ enum {DEBUG=FALSE}; /*@+redef@*/ if( DEBUG ) { /* these are the raw results */ fprintf( ioQQQ, " HMOLE raw; hi\t%.2e" , xIonFracs[ipHYDROGEN][0]); for( i=0; i < NMOLE; i++ ) { fprintf( ioQQQ, "\t%s\t%.2e", chLab[i], c[NMOLE][i] ); } fprintf( ioQQQ, " \n" ); } } if( trace.lgTrace && trace.lgTrMole ) { /* these are the raw results */ fprintf( ioQQQ, " raw; " ); for( i=0; i < NMOLE; i++ ) { fprintf( ioQQQ, " %s:%.2e", chLab[i], c[NMOLE][i] ); } fprintf( ioQQQ, " \n" ); /* these are the renormalixed */ fprintf( ioQQQ, " ren; " ); for( i=0; i < NMOLE; i++ ) { fprintf( ioQQQ, " %s:%.2e", chLab[i], c[NMOLE][i]* 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 = (xIonFracs[ipHYDROGEN][1] + xIonFracs[ipHYDROGEN][1])/(xIonFracs[ipHYDROGEN][1] + xIonFracs[ipHYDROGEN][1] + hmi.hminus + 2.*(hmi.htwo + hmi.h2plus) + 3*hmi.h3plus);*/ /* trace.lgTrMole is trace molecules option, * punch htwo */ if( (trace.lgTrace && trace.lgTrMole) || punch.lgPunH2 ) { if( punch.lgPunH2 ) { ioFile = punch.ipPunH2; } else { ioFile = ioQQQ; } rate = rgrain/* *xIonFracs[ipHYDROGEN][0]*/ + ratach*hmi.hminus + bh2dis*xIonFracs[ipHYDROGEN][0] + bh2h2p* xIonFracs[ipHYDROGEN][0]*hmi.h2plus + radasc*xIonFracs[ipHYDROGEN][0] + h3ph2p*hmi.h3plus + h3petc* hmi.h3plus; if( rate > 0. ) { 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/rate, ratach*hmi.hminus/rate, bh2dis/rate, bh2h2p*xIonFracs[ipHYDROGEN][0]*hmi.h2plus/rate, radasc*xIonFracs[ipHYDROGEN][0]/rate, h3ph2p*hmi.h3plus/rate, h3petc*hmi.h3plus/rate ); } else { fprintf( ioFile, " Create H2, rate=0\n" ); } } /* this is H2+ */ if( trace.lgTrace && trace.lgTrMole ) { rate = rh2h2p*hmi.htwo*xIonFracs[ipHYDROGEN][1] + b2pcin*xIonFracs[ipHYDROGEN][1]*phycon.eden*xIonFracs[ipHYDROGEN][0] + 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*xIonFracs[ipHYDROGEN][1]*hmi.htwo/rate, b2pcin*xIonFracs[ipHYDROGEN][1]*xIonFracs[ipHYDROGEN][1]* phycon.eden/rate, h3ph2p*hmi.h3plus/rate, h3petc*hmi.h3plus/ rate ); } else { fprintf( ioQQQ, " Create H2+, rate=0\n" ); } } if( hmi.hminus > 0. && phmlte > 0. ) { hmi.hmidep = (double)hmi.hminus/((double)xIonFracs[ipHYDROGEN][0]*phycon.eden*phmlte); } else { hmi.hmidep = 1.; } /* this will be net volume heating rate, photo heat - induc cool */ hmi.hmihet = HMinus_photo_heat*hmi.hminus - HMinus_induc_rec_cooling; /* save it */ heat.heating[0][15] = hmi.hmihet; heat.heating[0][16] = 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)(xIonFracs[ipHYDROGEN][0]) * (double)(xIonFracs[ipHYDROGEN][0]) * hmi.ph2lte; if( plte > 0. ) { hmi.h2dep = hmi.htwo/plte; } else { hmi.h2dep = 1.; } /* departure coef of H2+ defined rel to N(1s) N(p) * sec den was HI before 85.34 */ plte = xIonFracs[ipHYDROGEN][0]*xIonFracs[ipHYDROGEN][1]*phplte; if( plte > 0. ) { hmi.h2pdep = hmi.h2plus/plte; } else { hmi.h2pdep = 1.; } /* departure coef of H3+ defined rel to N(H2+) N(p) */ if( ph3lte > 0. ) { hmi.h3pdep = hmi.h3plus/ph3lte; } else { hmi.h3pdep = 1.; } if( trace.lgTrace && trace.lgTrMole ) { fprintf( ioQQQ, " HMOLE, Dep Coef, H-:%10.2e H2:%10.2e H2+:%10.2e\n", hmi.hmidep, hmi.h2dep, hmi.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, HMinus_induc_rec_rate, bhneut, c3bod, batach, hmi.hminus, hmi.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", hmi.HMinus_photo_rate, faneut, cionhm, ratach, iso.gamnc[ipH_LIKE][0][ipH1s], fhneut ); fprintf( ioQQQ, " H- heating:%10.3e Ind cooling%10.2e Spon cooling%10.2e\n", hmi.hmihet, HMinus_induc_rec_cooling, hmi.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 */ hmi.BiggestH2 = MAX2(hmi.BiggestH2,(float)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.; limit = MIN2(hmi.iheh2-1 , rfield.nflux ); for( i=hmi.iheh1-1; i < limit; i++ ) { gamheh += rfield.flux[i] + rfield.ConInterOut[i]+ rfield.outlin[i]; } gamheh *= 4e-18; /* both He+ + H and He + H+ * first rate is radiative associatino from * >>refer mheh+ rate Zygelman, B., and Dalgarno, A. 1990, ApJ 365, 239 */ hmi.hehp = (float)(1e-15*xIonFracs[ipHYDROGEN][0]*xIonFracs[ipHELIUM][1] + 1e-20*xIonFracs[ipHYDROGEN][1]*xIonFracs[ipHELIUM][0]); /* H2+ + HE => HEH+ + H0 */ hmi.hehp += (float)(3e-10*exp(-6717./phycon.te)*xIonFracs[ipHELIUM][0]*hmi.h2plus); /* hard radiation */ gamheh += 3.*iso.gamnc[ipH_LIKE][0][ipH1s]; /* HeH+ + H => H2+ + He */ gamheh += 1e-10*xIonFracs[ipHYDROGEN][0]; /* recombination, HeH+ + e => He + H */ if( gamheh > hmi.hehp*1e-15 ) { gamheh += phycon.eden*1e-9; } else { gamheh = -1; } hmi.hehp /= (float)gamheh; # ifdef DEBUG_FUN fputs( " <->hmole()\n", debug_fp ); # endif return; } /*hmirat compute radiative association rate for H- */ static double hmirat(double te) { double hmirat_v; # ifdef DEBUG_FUN fputs( "<+>hmirat()\n", debug_fp ); # endif /* fits to radiative association rate coefficients */ if( te < 31.62 ) { hmirat_v = 8.934e-18*phycon.sqrte*phycon.te003*phycon.te001* phycon.te001; } else if( te < 90. ) { hmirat_v = 5.159e-18*phycon.sqrte*phycon.te10*phycon.te03* phycon.te03*phycon.te003*phycon.te001; } else if( te < 1200. ) { hmirat_v = 2.042e-18*te/phycon.te10/phycon.te03; } else if( te < 3800. ) { hmirat_v = 8.861e-18*phycon.te70/phycon.te03/phycon.te01* phycon.te003; } else if( te <= 1.4e4 ) { /* following really only optimized up to 1e4 */ hmirat_v = 8.204e-17*phycon.sqrte/phycon.te10/phycon.te01* phycon.te003; } else { /* >>chng 00 sep 28, add this branch */ hmirat_v = 5.44e-16*phycon.te20/phycon.te01; } # ifdef DEBUG_FUN fputs( " <->hmirat()\n", debug_fp ); # endif return( hmirat_v ); }