/*coolr main routine to call others, to evaluate total cooling */ #include "cddefines.h" #include "showme.h" #include "physconst.h" #include "hydrogenic.h" #include "nhe1lvl.h" #include "taulines.h" #include "tff.h" #include "wind.h" #include "coolheavy.h" #include "tseton.h" #include "radius.h" #include "dclgas.h" #include "ffsum.h" #include "freeon.h" #include "tcool.h" #include "trace.h" #include "rh2dis.h" #include "cnegtk.h" #include "colneg.h" #include "called.h" #include "numderiv.h" #include "nomole.h" #include "doppvel.h" #include "cionhe.h" #include "he3lines.h" #include "pressure.h" #include "densty.h" #include "teinv.h" #include "compton.h" #include "logte.h" #include "he3tau.h" #include "hmi.h" #include "tepowers.h" #include "he.h" #include "cyclot.h" #include "h.h" #include "phycon.h" #include "he1ex.h" #include "zonecnt.h" #include "hmihet.h" #include "dndt.h" #include "covib.h" #include "cextra.h" #include "heating.h" #include "h2pcng.h" #include "he1deltat.h" #include "phe1lv.h" #include "he1stat.h" #include "he1blt.h" #include "ionfracs.h" #include "bit32.h" #include "he1crt.h" #include "sexp.h" #include "tfidle.h" #include "coladd.h" #include "coolmetals.h" #include "coolsum.h" #include "colzro.h" #include "widlin.h" #include "fndstr.h" #include "fndneg.h" #include "coolr.h" void coolr(double *tot) { long int ion, ipZ, iu, l, nelem; static long int nhit = 0, nzSave=0; double brems, BremsFactor , deriv, dhe1l, factor, he1ion, he1line, he1lmax, he1lone, pgrd, qn, recem, rothi, rotlow, tripcl, x; static double oltcool=0., oldtemp=0.; # ifdef DEBUG_FUN fputs( "<+>coolr()\n", debug_fp ); # endif /* returns tot, the total cooling, * and dc, the derivative of the cooling */ /* routine level2( t10 ) * routine level3( abund , t10,t21,t20) * tsq1 = 1. / (te**2) * POPEXC( O12,g1,g2,A21,excit,abund) ; result already*a21 * POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) * beseq(cs23,cs24,cs34,tarray,a41) * FIVEL( G(1-5) , ex(wn,1-5), cs12,cs13,14,15,23,24,25,34,35,45, * A21,31,41,51,32,42,52,43,53,54, pop(1-5), abund) */ if( trace.lgTrace ) { fprintf( ioQQQ, " COOLR called, TE=%12.4e\n", phycon.te ); } /* possibly update temperature variables * edensqte = eden/sqrte * cdsqte = 8.629e-6*eden/sqrte */ tfidle(); /* now zero out the cooling stack */ colzro(); /* grain recombination cooling, evaluated elsewhere * can either heat or cool the gas, do cooling here */ coladd("dust",0,MAX2(0.,dclgas.DustCoolGas)); /* grain cooling proportional to temperature ^3/2 */ tcool.dCooldT += MAX2(0.,dclgas.DustCoolGas)*3./(2.*phycon.te); /* molecular cooling */ if( nomole.lgNoH2Mole ) { CoolHeavy.hmfb = 0.; CoolHeavy.h2line = 0.; CoolHeavy.h2dis = 0.; /* H + H+ => H2+ cooling */ h2pcng.H2PlsCool = 0.; } else { /* H + H+ => H2+ cooling */ h2pcng.H2PlsCool = (float)(MAX2((2.325*phycon.te-1875.)*1e-20,0.)* h.hi*h.hii*1.66e-11); /* H- FB */ CoolHeavy.hmfb = hmihetCom.hmicol; /* H- FF is in with H ff * rate for rotation lines from * >>refer Lepp, S., & Shull, J.M. 1983, ApJ, 270, 578 */ x = logte.alogte - 4.; if( phycon.te > 1087. ) { rothi = 3.90e-19*sexp(6118./phycon.te); } else { rothi = pow(10.,-19.24 + 0.474*x - 1.247*x*x); } /* low density rotation cooling */ qn = pow(MAX2(hmi.htwo,1e-37),0.77) + 1.2*pow(MAX2(h.hi,1e-37),0.77); if( phycon.te > 4031. ) { rotlow = 1.38e-22*sexp(9243./phycon.te)*qn; } else { rotlow = pow(10.,-22.90 - 0.553*x - 1.148*x*x)*qn; } if( rotlow > 0. ) { CoolHeavy.h2line = hmi.htwo*rothi/(1. + rothi/rotlow); } else { CoolHeavy.h2line = 0.; } CoolHeavy.h2dis = (float)(rh2disCom.rh2dis*h.hi*hmi.htwo*7.19e-12* 0.); } coladd("H-fb",0,CoolHeavy.hmfb); /* >>chng 96 nov 15, fac of 2 in deriv to help convergence in very dense * models where H- is important, this takes change in eden into * partial account */ tcool.dCooldT += 2.*CoolHeavy.hmfb*teinvCom.teinv; coladd("H2ln",0,CoolHeavy.h2line); tcool.dCooldT += CoolHeavy.h2line*teinvCom.teinv; coladd("H2ds",0,CoolHeavy.h2dis); tcool.dCooldT += CoolHeavy.h2dis*teinvCom.teinv; coladd("H2+ ",0,h2pcng.H2PlsCool); tcool.dCooldT += h2pcng.H2PlsCool*teinvCom.teinv; /* carbon monoxide (CO) vib-rot cooling */ coladd("CO C",0,covib.COVibCool); tcool.dCooldT += covib.dCOVdT; /* heating due to three-body, will be incremented in HydroCool*/ HeatingCom.heating[3][0] = 0.; /* heating due to hydrogen lines */ HeatingCom.heating[23][0] = 0.; /* heating due to photoionization of all excited states of hydrogen species */ HeatingCom.heating[1][0] = 0.; /* hydrogenic species cooling, mainly lines, and level ionization */ for( ipZ=0; ipZ < LIMELM; ipZ++ ) { /* must always call HydroCool since we must zero those variables * that would have been set had the species been present */ HydroCool(ipZ); } /* brems cooling */ /* Hydrogen BremsFactor free free cooling, free-free heating */ /* NB major change in code with extension of continuum down to 1e-8 ryd * firstdr does not check opacity below 1e-5 - if free--free heating runs * away could be strong fir continuum hitting free free opacity */ factor = 0.9 + 1.2/(2.2 + fabs(logte.alogte-5.398)); /* >>chng 96 jan 8, for stability always use corrected free free */ /* tffCom.tff is freq (Ryd) where gas becomes optically thin to absorption */ BremsFactor = 1.42e-27*phycon.eden*TePowers.sqrte*factor*sexp(tffCom.tff* TE1RYD/phycon.te); /* charge squared is hydrogenic brems factor */ hydro.FreeFreeCool = 0.; /* following set FALSE with no free-free command */ if( freeon.lgFreeOn ) { for( ipZ=0; ipZ>chng 96 oct 30, from HFFNet to just FreeFreeCool, * since SumHeat picks up FreeFreeHeat */ tcool.dCooldT += hydro.FreeFreeCool*tcool.halfte; /* net cooling for helium, HeI * HE1X evaluated in HE1LEV */ he1ion = he1exCOM.he1ex*he.heii; coladd("He1i",0,he1ion); tcool.dCooldT += he1ion*3.9e4*tcool.tsq1; /* all cooling due to hei lines */ if( he.heii > 0. ) { he1line = 0.; he1lmax = 0.; dhe1l = 0.; /* first do 1s-2s and 1s-2p collisions */ he1line = phycon.eden*he.heiii*BOLTZMANN*He1deltaT.He1deltaTe[0][1]* (phe1lv.he1n[0]*he1blt.he1bol[1][0]*(he1crt.he1c2s1s + 3.* phe1lv.he1cll[0][1]) - phe1lv.he12s*he1crt.he1c2s1s - phe1lv.he12p* phe1lv.he1cll[0][1]); he1lmax = he1line; he1lone = he1line; if( he1lone > 0. ) { /* line is coolant, note derivative counted from ground */ dhe1l = he1lone*(He1deltaT.He1deltaTe[0][1]*tcool.tsq1 - tcool.halfte); } else { /* line is a heating agent */ dhe1l = -he1lone*tcool.halfte; } /* upper limits changed to 5 and 6 due to strong cooling from * highest psuedo state */ for( l=0; l < 5; l++ ) { for( iu=MAX2(2,l+1); iu < 6; iu++ ) { he1lone = phe1lv.he1cll[l][iu]*(phe1lv.he1n[l]* he1blt.he1bol[iu][l]*he1statCOM.he1stat[iu]/he1statCOM.he1stat[l] - phe1lv.he1n[iu])*phycon.eden*he.heii*(BOLTZMANN* He1deltaT.He1deltaTe[l][iu]); /* 3 eden*heii* (1.384e-16 * HdeltaT(iu,l)*4. ) * >>chng 97 mar 4, above heii was heiii */ if( fabs(he1lone) > fabs(he1lmax) ) { he1lmax = he1lone; } he1line += he1lone; /* >>chng 97 july 1, add in cooling derivative */ if( he1lone > 0. ) { /* line is coolant, note derivative counted from ground */ dhe1l += he1lone*(He1deltaT.He1deltaTe[0][iu]* tcool.tsq1 - tcool.halfte); } else { /* line is a heating agent */ dhe1l += -he1lone*tcool.halfte; } } } } else { dhe1l = 0.; he1line = 0.; } /* all HeI line collisional excitation */ coladd("He1l",0,he1line); tcool.dCooldT += dhe1l; /* total Compton cooling */ CoolHeavy.tccool = Compton.cmcool*phycon.te; coladd("Comp",0,CoolHeavy.tccool); tcool.dCooldT += Compton.cmcool; /* option for "extra" cooling, expressed as power-law */ if( cextra.lgCExtraOn ) { CoolHeavy.cextxx = (float)(cextra.CoolExtra*pow(phycon.te/1e4,cextra.cextpw)); } else { CoolHeavy.cextxx = 0.; } coladd("Extr",0,CoolHeavy.cextxx); /* cooling due to wind expansion, only for winds expansion cooling */ if( wind.windv > 0. ) { dndt.dDensityDT = (float)(wind.AccelTot/wind.windv + 2.*wind.windv/ radius.Radius); CoolHeavy.expans = densty.pden*phycon.te*BOLTZMANN*dndt.dDensityDT; } else { dndt.dDensityDT = 0.; CoolHeavy.expans = 0.; } coladd("Expn",0,CoolHeavy.expans); tcool.dCooldT += CoolHeavy.expans/phycon.te; /* cyclotron cooling */ CoolHeavy.cyntrn = cyclot.CycloCoolCoef*phycon.eden*phycon.te; coladd("Cycl",0,CoolHeavy.cyntrn); tcool.dCooldT += cyclot.CycloCoolCoef*phycon.eden; /* electron-electron brems, approx form from * >>refer Stepney and Guilbert, MNRAS 204, 1269 (1983) * ok for T<10**9 */ CoolHeavy.eebrm = POW2(phycon.eden*phycon.te*1.84e-21); /* >>chng 97 mar 12, added deriv */ tcool.dCooldT += CoolHeavy.eebrm*tcool.halfte; coladd("eeff",0,CoolHeavy.eebrm); /* collisional ionization of triplet helium */ coladd("HeI3",0,MAX2(0.,cionhe.He3IonCool)); tcool.dCooldT += MAX2(0.,cionhe.He3IonCool)*5.53e4*tcool.tsq1; /* net three body heating * >>chng 96 july 15, had always been cooling */ HeatingCom.heating[4][1] = MAX2(0.,-cionhe.He3IonCool); /* collisional excitation of all triplet lines, done in he3lev * can heat or cool */ if( he.hei3 > 0. ) { tripcl = he3lines.he3clx; } else { tripcl = 0.; } coladd("Hetr",0,MAX2(0.,tripcl)); HeatingCom.heating[3][1] = MAX2(0.,-tripcl); /* >>chng 96 nov 19, to following expression for HeI Lya cooling *666 error! break out 2s and 2p, really should write entire he cooling rout */ CoolHeavy.c584 = phycon.eden*(double)he.heii*3.41e-11*phe1lv.he1cll[0][1]* (he1blt.he1bol[1][0]*phe1lv.he1n[0]*4. - phe1lv.he1n[1]); tcool.dCooldT += CoolHeavy.c584*2.46e5*tcool.tsq1; coladd("He I",584,CoolHeavy.c584); /* radiation pressure due to HE I 10830 */ pgrd = he3lines.p2s - he3lines.p2p/3.; if( pgrd != 0. ) { pressure.plines[3] = (float)( 2.69e-10*he3lines.p2p/(he3lines.p2s - he3lines.p2p/3.)* widlin(he3tau[IPT10830-1].TauIn, he3tau[IPT10830-1].TauTot, 1.37e-4,DoppVel.doppler[1]) ); } else { pressure.plines[3] = 0.; } /* He I brems (free-free) cooling */ factor = 0.9 + 1.2/(2.2 + fabs(logte.alogte-5.398)); /* >>chng 95 dec 23 always do ots for brems for stability * ots continuum transport - do not count optically thick part */ brems = 1.42e-27*phycon.eden*TePowers.sqrte*factor*sexp(tffCom.tff* TE1RYD/phycon.te); /* do not include ion since was done in hydrogenic above */ /* following set FALSE with no free-free command */ if( freeon.lgFreeOn ) { CoolHeavy.bhe = brems*he.heii; } else { CoolHeavy.bhe = 0.; } coladd("Heff",0,CoolHeavy.bhe); /* recombination cooling for He+ to Heo */ recem = 9.52e-26*TePowers.te10*phycon.eden; CoolHeavy.herec = recem*he.heii; tcool.dCooldT += CoolHeavy.bhe*tcool.halfte + CoolHeavy.herec* 0.27*teinvCom.teinv; coladd("Hefb",0,CoolHeavy.herec); /* heavy element recombination cooling, do not count hydrogenic since * already done above, also helium singlets have been done * EdenFFSum set in ESUM, given by Harrington et al. 82 * 666 error! put in charge dependent rec term */ CoolHeavy.heavfb = 0.; for( nelem=2; nelem < LIMELM; nelem++ ) { for( ion=2; ion <= nelem ; ion++ ) { CoolHeavy.heavfb += xIonFracs[ion][nelem]; } } CoolHeavy.heavfb *= recem; coladd("hvFB",0,CoolHeavy.heavfb); /* heavy element brems (free-free) cooling */ if( freeon.lgFreeOn ) { CoolHeavy.heavff = ffsum.EdenFFSum*1.44e-27*TePowers.sqrte*phycon.eden* factor; } else { CoolHeavy.heavff = 0.; } tcool.dCooldT += CoolHeavy.heavfb*.113*teinvCom.teinv; coladd("Hvff",0,CoolHeavy.heavff); /* heavy element collisional ionization * derivative should be zero since increased coll ion rate * decreases neutral fraction by proportional amount */ coladd("Hvin",0,CoolHeavy.colmet); /* Carbon cooling */ CoolCarb(); /* Nitrogen cooling */ CoolNitr(); /* Oxygen cooling */ CoolOxyg(); /* Fluorine cooling */ CoolFluo(); /* Neon cooling */ CoolNeon(); /* Magnesium cooling */ CoolMagn(); /* Sodium cooling */ CoolSodi(); /* Aluminum cooling */ CoolAlum(); /* Silicon cooling */ CoolSili(); /* Phosphorus */ CoolPhos(); /* Sulphur cooling */ CoolSulf(); /* Chlorine cooling */ CoolChlo(); /* Argon cooling */ CoolArgo(); /* Potasium cooling */ CoolPota(); /* Calcium cooling */ CoolCalc(); /* Scandium cooling */ CoolScan(); /* Titanium cooling */ CoolTita(); /* Vanadium cooling */ CoolVana(); /* Chromium cooling */ CoolChro(); /* Iron cooling */ CoolIron(); /* Manganese cooling */ CoolMang(); /* Cobalt cooling */ CoolCoba(); /* Nickel cooling */ CoolNick(); /* do all the thousands of lines Dima added with g-bar approximation */ CoolDima(); /* now add up all the coolants */ CoolSum(tot); /* negative cooling */ if( *tot <= 0. ) { fprintf( ioQQQ, " COOLR; cooling is <=0, this is impossible.\n" ); ShowMe(); puts( "[Stop in coolr]" ); exit(1); } /* bad derivative */ if( tcool.dCooldT == 0. ) { fprintf( ioQQQ, " COOLR; cooling slope <=0, this is impossible.\n" ); if( (*tot > 0. && bit32.lgBit32) && phycon.hden < 1e-4 ) { fprintf( ioQQQ, " Probably due to very low density.\n" ); } ShowMe(); puts( "[Stop in coolr]" ); exit(1); } if( trace.lgTrace ) { fndstr(*tot,tcool.dCooldT); } /* lgTSetOn true for constant temperature model */ if( (((((!tseton.lgTSetOn) && *tot < 0.) && called.lgTalk) && colneg.lgColNeg) && cnegtk.lgCNegChk) && ZoneCnt.nzone > 0 ) { fprintf( ioQQQ, " Negative cooling, zone%4ld, =%10.2e coola=%10.2e CHION=%10.2e Te=%10.2e\n", ZoneCnt.nzone, *tot, hydro.coola[0], hydro.chion[0], phycon.te ); fndneg(); } /* possibility of getting emperical cooling derivative * normally false, set true with 'set numerical derivatives' command */ if( NumDeriv.lgNumDeriv ) { if( ((ZoneCnt.nzone > 2 && ZoneCnt.nzone == nzSave) && oldtemp != phycon.te) && nhit > 4 ) { /* hnit is number of trys on this zone - use to stop numerical problems * do not evaluate numerical deriv until well into solution */ deriv = (oltcool - *tot)/(oldtemp - phycon.te); tcool.dCooldT = deriv; } else { deriv = tcool.dCooldT; } if( ZoneCnt.nzone != nzSave ) nhit = 0; nzSave = ZoneCnt.nzone; nhit += 1; oltcool = *tot; oldtemp = phycon.te; } # ifdef DEBUG_FUN fputs( " <->coolr()\n", debug_fp ); # endif return; }