/*LineSet1 put energetics, H, and He lines into line intensity stack */ #include #include "cddefines.h" #include "taulines.h" /****************************************************************************** * note - order of routines in this file is important since hazy iv * gets a table from here. do not split this file into pieces. ****************************************************************************** */ #include "physconst.h" #include "nhe1lvl.h" #include "linesave.h" #include "coolheavy.h" #include "secondaries.h" #include "dclgas.h" #include "hydrogenic.h" #include "grainvar.h" #include "opac.h" #include "elmton.h" #include "totlum.h" #include "covib.h" #include "heots.h" #include "h.h" #include "heat.h" #include "heating.h" #include "heatlmax.h" #include "hsrate.h" #include "he1ex.h" #include "he.h" #include "smbeta.h" #include "phycon.h" #include "ionfracs.h" #include "logte.h" #include "dphoht.h" #include "hmihet.h" #include "compton.h" #include "phe1lv.h" #include "heinwd.h" #include "he3c.h" #include "he3tau.h" #include "he3lines.h" #include "dtrans.h" #include "sphere.h" #include "cionhe.h" #include "codish.h" #include "h2pcng.h" #include "h2dish.h" #include "hfbnet.h" #include "dhla.h" #include "rfield.h" #include "hahmin.h" #include "tepowers.h" #include "rin2ph.h" #include "hcaseb.h" #include "predcont.h" #include "trace.h" #include "photrate.h" #include "ip626.h" #include "iphe1l.h" #include "ionrange.h" #include "linadd.h" #include "lindst.h" #include "iwavlen.h" #include "sexp.h" #include "radius.h" #include "pntforline.h" #include "putline.h" #include "putextra.h" #include "getmaxhline.h" void LineSet1(void) { long int i, ipHi, ipLo, ipZ, ipnt, low, nd; double ConIn, ConOt, ConRef , cb206, cb5016, col, dhtot, e228, e504, e584, e626, e911, eff, em, EmisFac ,/* used for hydrogenic emissivity or intensity */ f10830, flow, hbetab, hbetac, HeatMetal , ee511, hlalph, pop, pump, rec, totl; # ifdef DEBUG_FUN fputs( "<+>LineSet1()\n", debug_fp ); # endif if( trace.lgTrace ) { fprintf( ioQQQ, " LineSet1 called\n" ); } /* initialize ExtraInten with a zero */ PutExtra(0.); /* put in something impossible in element 0 of line stack */ linadd(0.f,0,"zero",'i'); /* do we want intensity/luminosity or emissivity */ if( hydro.lgHydEmiss ) { /* option to print emissivity, not intensity */ EmisFac = 1./phycon.eden; } else { EmisFac = xIonFracs[2][0]; } /* total H-beta from multi-level atom */ ipZ = 0; ipHi = 4; ipLo = IP2P; hbetac = (HydroLines[ipZ][ipHi][IP2S].Aul * HydroLines[ipZ][ipHi][IP2S].Pesc + HydroLines[ipZ][ipHi][IP2P].Aul * HydroLines[ipZ][ipHi][IP2P].Pesc ) * HydroLines[ipZ][ipHi][ipLo].PopHi* HydroLines[ipZ][ipHi][ipLo].EnergyErg* EmisFac ; /* total H-beta from multi-level atom */ /*hbetac = HydroLines[0][4][IP2S].xIntensity + */ /*HydroLines[0][4][IP2P].xIntensity;*/ /* these lines added to outlin in metdif - following must be false */ lindst(hbetac,4861,"TOTL",HydroLines[0][4][2].ipCont, 'i',FALSE); /* total Ly-a from multi-level atom */ ipHi = IP2P; ipLo = IP1S; hlalph = HydroLines[ipZ][ipHi][ipLo].Aul* HydroLines[ipZ][ipHi][ipLo].PopHi* HydroLines[ipZ][ipHi][ipLo].Pesc* HydroLines[ipZ][ipHi][ipLo].EnergyErg* EmisFac ; /*hlalph = HydroLines[0][IP2P][IP1S].xIntensity;*/ linadd(hlalph,1216,"TOTL",'i'); totlum.totlsv = totlum.totlsv/(float)radius.dVeff*sphere.covgeo; /* someday, put in 21 cm, A from Gould, ApJ 423, 522, =2.88426e-15 s-1 */ /* total luminosity in incident continuum */ linadd(totlum.totlsv,0,"Inci",'i'); totlum.totlsv = totlum.totlsv*(float)radius.dVeff/sphere.covgeo; /* ipass is flag to indicate whether to only set up line array * (ipass=0) or actually evaluate lines intensities (ipass=1) */ if( LineSave.ipass > 0 ) { totlum.totlsv = 0.; } /* total heating, all forms, information since individuals added later */ linadd(heat.htot,0,"TotH",'i'); /* hydrogen photoionization heating, ground state only */ linadd(HeatingCom.heating[0][0],0,"BFH1",'h'); /* net hydrogen photoionization heating less rec cooling, all excited states * normally zero, positive if excited states are net heating */ linadd(HeatingCom.heating[1][0],0,"BFHx",'h'); if( heat.htot > 0. ) { if( HeatingCom.heating[22][0]/heat.htot > heatlmax.HeatLineMax ) { heatlmax.HeatLineMax = (float)(HeatingCom.heating[22][0]/heat.htot); GetMaxhLine(); } } /* heating due to induced lines absorption of continuum */ linadd(HeatingCom.heating[22][0],0,"Line",'h'); /* net cooling due to collisional ionization of Heo */ linadd(MAX2(0.,he1exCOM.he1ex*he.heii),0,"He1i",'c'); /* this is the heating due to 3-body recombination */ linadd(MAX2(-he1exCOM.he1ex*he.heii,0.),0,"3He1",'h'); /* total helium photoionization heating, all stages */ linadd(HeatingCom.heating[0][1]+HeatingCom.heating[1][1]+HeatingCom.heating[2][1], 0,"BFHe",'h'); HeatMetal = 0.; /* some sums that will be printed in the stack */ for( ipZ=2; ipZ 0 ) { smbeta.SimHBeta = 0.; } /* this is main loop creating hydrogenic intensity or emissivity */ for( ipZ=0; ipZ < LIMELM; ipZ++ ) { if( hydro.lgHydEmiss ) { /* option set with Hydrogen emissivity command * will produce emissivity not intensity - luminosity */ for( ipZ=0; ipZ < LIMELM; ipZ++ ) { if( elmton.lgElmtOn[ipZ] ) { for( ipHi=IP2S; ipHi <= nhlevel; ipHi++ ) { for( ipLo=IP1S; ipLo <= (ipHi - 1); ipLo++ ) { /* this is in real units not emissivit*/ HydroLines[ipZ][ipHi][ipLo].phots = (float)( HydroLines[ipZ][ipHi][ipLo].Aul* HydroLines[ipZ][ipHi][ipLo].PopHi* HydroLines[ipZ][ipHi][ipLo].Pesc* xIonFracs[ipZ+2][ipZ]); /* now find line intensity */ HydroLines[ipZ][ipHi][ipLo].xIntensity = (float)( HydroLines[ipZ][ipHi][ipLo].Aul* HydroLines[ipZ][ipHi][ipLo].PopHi* HydroLines[ipZ][ipHi][ipLo].Pesc* HydroLines[ipZ][ipHi][ipLo].EnergyErg/phycon.eden); } } } } } else { /* this is the default case - luminosity - intensity */ for( ipZ=0; ipZ < LIMELM; ipZ++ ) { if( elmton.lgElmtOn[ipZ] ) { for( ipHi=IP2S; ipHi <= nhlevel; ipHi++ ) { for( ipLo=IP1S; ipLo <= (ipHi - 1); ipLo++ ) { HydroLines[ipZ][ipHi][ipLo].phots = HydroLines[ipZ][ipHi][ipLo].Aul* HydroLines[ipZ][ipHi][ipLo].PopHi* HydroLines[ipZ][ipHi][ipLo].Pesc* xIonFracs[ipZ+2][ipZ]; /* now find line intensity and line ots rates */ /* set hydrogen line emission and ots fields */ HydroLines[ipZ][ipHi][ipLo].xIntensity = HydroLines[ipZ][ipHi][ipLo].Aul* HydroLines[ipZ][ipHi][ipLo].PopHi* HydroLines[ipZ][ipHi][ipLo].Pesc* xIonFracs[ipZ+2][ipZ]* HydroLines[ipZ][ipHi][ipLo].EnergyErg; } } } } } } /* create emissivity or intensity for hydrogenic species, * first combine/bring balmer series together */ for( ipZ=0; ipZ < LIMELM; ipZ++ ) { if( IonRange.IonHigh[ipZ] == ipZ + 2 ) { /* bring n=> 2s + 2p together into one line with 2s lower level * must be 2s since it comes first in line stack */ for( ipHi=3; ipHi <= nhlevel; ipHi++ ) { HydroLines[ipZ][ipHi][IP2S].xIntensity += HydroLines[ipZ][ipHi][IP2P].xIntensity; /* kill 2p */ HydroLines[ipZ][ipHi][IP2P].xIntensity = 0.; } } } /* H beta recombination, assuming old case B */ hbetab = (double)((pow(10.,-20.89 - 0.10612*POW2(logte.alogte - 4.4)))/ phycon.te); /* convert to emission per unit then dump into line stack */ if( !hydro.lgHydEmiss ) { hbetab *= xIonFracs[2][0]*phycon.eden; } /* this is old case b, had been in code for very long time */ linadd(hbetab,4861,"CaBo",'i'); /* 1640 1640 1640 */ em = 2.03e-20/(TePowers.te70*TePowers.te10*TePowers.te03); if( !hydro.lgHydEmiss ) { em *= xIonFracs[3][1]*phycon.eden; } /* old prediction of He II 1640, case B at low densities */ linadd(em,1640,"CaBo",'i'); /* hydrogenic helium */ /* old prediction of He II 4686, case B */ em = 2.52e-20/(pow(phycon.te,1.05881)); if( !hydro.lgHydEmiss ) { em *= xIonFracs[3][1]*phycon.eden; } linadd(em,4686,"CaBo",'i'); /* predict case b intensities of hydrogen lines */ ipZ = 0; for( ipZ=0; ipZ<2; ++ipZ ) { for( ipLo=IP2P; ipLo<6; ++ipLo ) { for( ipHi=ipLo+1; ipHi< ipLo+4; ++ipHi ) { /* Put case b predictions into line stack * NB NB NB each Hummer & Storey case b line must be * explicitly clobbered by hand in rutine final * if not ok flag is set since this indicates that we exceeded bounds of table, * DO NOT want to print lines in that case */ /* first do case b emssivity of balmer lines */ /* get HS predictions for Ha */ hbetab = (double)HSRate( ipHi,ipLo , ipZ+1, phycon.te , phycon.eden ); if( hbetab<=0. ) { /* routine returns -1 if outside bounds of table*/ if( ipZ==0 ) { /* hydrogen */ CaseBHS.lgHCaseBOK = FALSE; } else { /* helium */ CaseBHS.lgHeCaseBOK = FALSE; } hbetab = 0.f; } /* convert to emission per unit then dump into line stack */ if( !hydro.lgHydEmiss ) { hbetab *= xIonFracs[ipZ+2][ipZ]*phycon.eden; } /* this is interpolated case b from HS tables */ i = iWavLen(&HydroLines[ipZ][ipHi][ipLo] ); linadd(hbetab,i,"Ca B",'i'); } } } /* this is the main printout, where line intensities are entered into the stack */ for( ipZ=0; ipZ < LIMELM; ipZ++ ) { if( elmton.lgElmtOn[ipZ] ) { /* do not print more than first 50 levels, since my pc does not * have enough memory */ /* NB NB - low and high must be in this order so that all balmer, paschen, * etc series line up correctly in final printout */ for( ipLo=IP1S; ipLo < MIN2(50,nhlevel-1); ipLo++ ) { /* don't bother with decays to 2p ince we set them to zero above */ if( ipLo==IP2P ) continue; for( ipHi=MAX2(IP2P,ipLo+1); ipHi <= MIN2(50,nhlevel); ipHi++ ) { /* this will enter .xIntensity into the line stack */ PutLine(&HydroLines[ipZ][ipHi][ipLo]); } } } } /* total collisional cooling due to all hydrogen lines */ linadd(MAX2(0.,hydro.HLineTotCool[0]),912,"Clin",'c'); /* total collisional heating due to all hydrogen lines */ linadd(MAX2(0.,-hydro.HLineTotCool[0]),912,"Hlin",'h'); /* 2-photon two photon emission from multi-level atom */ linadd(h.hii*hydro.hn[0][IP2S]*8.23*1.634e-11,0,"2 NU",'r'); /* la contribution from suprathermal secondaries from ground */ linadd(Secondaries.x12tot*h.hii*hydro.hn[0][IP1S]*1.634e-11,1216,"LA X" ,'i'); /* H-alpha produce by H- mutual neutralization */ linadd(hahmin.HalphaHmin*3.032e-12,6563,"H-CT",'i'); /* "Ly alpha" produced by induced two photon */ linadd(hydro.hn[0][IP2S]*h.hii*rin2ph.ri2s1s*1.634e-11,1216, "Ind2" ,'i'); /* H-beta produced by continuum pumping in optically thin ld limit */ /* factor of 0.4836 is ratio of A(4-2)/(A(4-3)+A(4-2)) * the IPLNPUMP is the actual pumping rate per atom */ pump = (double)(HydroLines[0][4][IP1S].pump*hydro.hn[0][IP1S]* h.hii*4.09e-12*0.4836); linadd(pump,4861,"Pump",'r'); /* CHION is HEX*HII (in COOLR) and is net col ionz-3 body heat * collision ionization cooling of hydrogen */ linadd(MAX2(0.,hydro.chion[0]),0,"CION",'c'); /* this is the heating due to 3-body recombination */ linadd(MAX2(-hydro.chion[0],0.),0,"3bHt",'h'); /*fixit(); stark term not now added to total rad decay rate - * so energy not right in this printout - why is stark * not added */ /* Stark broadening contribution to line */ linadd(h.hii*hydro.hn[0][IP2P]*0.*hydro.pestrk[2][0]*1.634e-11, 1216,"Strk",'i'); /* Stark broadening contribution to line */ linadd(h.hii*hydro.hn[0][3]*hydro.pestrk[3][2]*3.025e-12, 6563,"Strk",'i'); /* Stark broadening contribution to line */ linadd(h.hii*hydro.hn[0][4]*hydro.pestrk[4][2]*4.084e-12, 4861,"Strk",'i'); /* Stark broadening contribution to line */ linadd(h.hii*hydro.hn[0][4]*hydro.pestrk[4][3]*1.059e-12, 18751,"Strk",'i'); /* pestrk[5,4] is A[4,5]*pest[4,5] * Stark broadening contribution to line */ linadd(h.hii*hydro.hn[0][5]*hydro.pestrk[5][4]*4.900e-13, 40512,"Strk",'i'); /* this can fail if RTMake never updates the ots rates, a logic error, * but only assert this during actual calculation (ipass>0), */ assert( LineSave.ipass <1 || HydroLines[0][IP2P][IP1S].ots>= 0.); /* portion of line lost due to absorp by background opacity */ linadd(HydroLines[0][IP2P][IP1S].ots*1.634e-11, 1216,"Dest",'i'); /* portion of line lost due to absorp by background opacity */ linadd(HydroLines[0][3][IP2S].ots*3.025e-12, 6563,"Dest",'i'); /* portion of line lost due to absorp by background opacity */ linadd(HydroLines[0][5][4].ots*4.915e-13,40516, "Dest",'i'); /* portion of line lost due to absorp by background opacity */ linadd(HydroLines[0][4][IP2S].ots*4.084e-12, 4861,"Dest",'i'); /* portion of line lost due to absorp by background opacity */ linadd(HydroLines[0][4][3].ots*1.059e-12,18751, "Dest",'i'); /* Ly-alpha destroyed by overlap with FeII */ linadd(hydro.hn[0][IP2P]*h.hii*HydroLines[0][IP2P][IP1S].Aul* hydro.dstfe2lya*1.637e-11,1216,"Fe 2",'i'); /* 511keV annihilation line */ ee511 = (phycon.hden + 4.*he.ahe)*PhotRate.PhotoRate[0][0][10][0]*2.*8.20e-7; PntForLine(2.427e-2,"e-e+",&ipnt); lindst(ee511,511,"e-e+",ipnt,'i',TRUE); /* predict emitted continuum at series of continuum points */ /* variables are located in predcont.h, * NPREDCONT number of continuum points, * next two arrays are defined in zerologic.c * EnrPredCont - energy (Ryd) where we want to predict the continuum, * these are set in zerologic.c, * ipPreCont - pointers to energy within continuum array, * these are set in CreatePoint * * the entry nFnu will only be printed if the command * print diffuse * is entered - */ for( i=0; i < NPREDCONT; i++ ) { /* reflected toal continuum (dif+incident emitted inward direction) */ ConRef = rfield.refcon[ipPredCont[i]-1]/rfield.widflx[ipPredCont[i]-1]* rfield.anu2[ipPredCont[i]-1]*EN1RYD/opac.tmn[ipPredCont[i]-1]/radius.GeoDil; /* reflected part of INCIDENT continuum (only incident) */ ConIn = rfield.reflux[ipPredCont[i]-1]/rfield.widflx[ipPredCont[i]-1]* rfield.anu2[ipPredCont[i]-1]*EN1RYD/opac.tmn[ipPredCont[i]-1]/radius.GeoDil; /* outward continuum */ ConOt = (rfield.outcon[ipPredCont[i]-1]+ rfield.ConOutNoInter[ipPredCont[i]-1])/ rfield.widflx[ipPredCont[i]-1]* rfield.anu2[ipPredCont[i]-1]*EN1RYD*radius.r1r0sq/ opac.tmn[ipPredCont[i]-1]/radius.GeoDil; /* this is wavelength in Angstroms, energies stored in EnrPredCont */ low = (long)(RYDLAM/EnrPredCont[i]); /* some are hundreds of microns, overflows in printout, cnvrt to microns */ if( low > 100000 ) { low /= 10000; } /* memory has not been allocated until ipass >= 0, so must not access this element, * this element will be used to save the following quantity */ if( LineSave.ipass > 0 ) { LineSv[LineSave.nsum].sumlin = 0.; } /* total continuum produced by cloud at selected energy points */ linadd((ConRef+ConOt)/radius.dVeff,low,"nFnu",'i'); /* here we must correct unit emissivity which is incorrect from above */ if( LineSave.ipass > 0 ) { LineSv[LineSave.nsum-1].emslin = (float)( rfield.DiffLocal[ipPredCont[i]-1]/rfield.widflx[ipPredCont[i]-1]* rfield.anu2[ipPredCont[i]-1]*EN1RYD); } /* memory has not been allocated until ipass >= 0 */ if( LineSave.ipass > 0 ) { LineSv[LineSave.nsum].sumlin = 0.; } /* total reflected continuum, total inward emission plus * reflected diffuse continuum */ linadd(ConRef/radius.dVeff,low,"InwT",'i'); if( LineSave.ipass > 0 ) { LineSv[LineSave.nsum-1].emslin = 0.; } /* memory has not been allocated until ipass >= 0 */ if( LineSave.ipass > 0 ) { LineSv[LineSave.nsum].sumlin = 0.; } /* reflected incident continuum (only incident) */ linadd(ConIn/radius.dVeff,low,"InwC",'i'); if( LineSave.ipass > 0 ) { LineSv[LineSave.nsum-1].emslin = 0.; } /*fprintf( ioQQQ,"%li\n",low );*/ } /* now find incident continuum */ smbeta.cn4861 /= (float)radius.dVeff; /* incident continuum nu*f_nu at H-beta, at illuminated face of cloud */ linadd(smbeta.cn4861,4860,"Inci",'i'); smbeta.cn1216 /= (float)radius.dVeff; /* incident continuum nu*f_nu near Ly-alpha, at illuminated face of cloud */ linadd(smbeta.cn1216,1215,"Inci",'i'); smbeta.cn1216 *= (float)radius.dVeff; smbeta.cn4861 *= (float)radius.dVeff; if( LineSave.ipass > 0 ) { smbeta.cn4861 = 0.; smbeta.cn1216 = 0.; } /* integrated Balmer continuum emission */ flow = (hydro.hrec[0][IP2P][ipRecRad] + hydro.hrec[0][IP2S][ipRecRad])* hydro.hrec[0][IP2P][ipRecEsc]*phycon.eden*xIonFracs[2][0]* 5.45e-12; linadd(flow,0,"Ba C",'i'); /* >>chng 96 june 13, had been hrec(3,1,ipRecEsc) * Paschen continuum emission */ flow = hydro.hrec[0][3][ipRecRad]*hydro.hrec[0][3][ipRecEsc]* phycon.eden*xIonFracs[2][0]*3.53e-12; linadd(flow,0,"PA C",'i'); /* total grain heating by all sources, lines, collisions, incident continuum */ linadd(dhla.dstotl,0,"GraT",'i'); /* grain heating by incident continuum */ linadd(dhla.dincon,0,"GraI",'i'); /* grain heating due to destruction of Ly alpha */ linadd(dhla.DustHeatLya,1216,"GraL",'i'); /* grain heating due to collisions with gas */ linadd(dhla.dhcol,0,"GraC",'i'); /* grain heating due to diffuse fields, may also have grain emission */ linadd(dhla.HDustDif,0,"GraD",'i'); /* part of H brems, in x-ray beyond 0.5KeV */ linadd((hydro.FreeFreeCool+CoolHeavy.bhe)*sexp(5.8e6/phycon.te),0,"FF X" ,'i'); /* total Compton cooling */ linadd(CoolHeavy.tccool,0,"ComC",'c'); /* expansion cooling, only non-zero for wind */ linadd(CoolHeavy.expans,0,"Expn",'c'); /* electron - electron brems */ linadd(CoolHeavy.eebrm,0,"eeff",'c'); /* H radiative recombination cooling */ linadd(hydro.HRadRecCool[0],0,"H FB",'i'); /* net free-bound cooling and heating */ linadd(MAX2(0.,HFBNet.HFreBndNet[0]),0,"HFBc",'c'); linadd(MAX2(0.,-HFBNet.HFreBndNet[0]),0,"HFBh",'h'); /* cooling due to induced rec of hydrogen */ linadd(HFBNet.crcind[0],0,"Hind",'c'); /* cooling due to induced rec of fully ionized helium */ linadd(CoolHeavy.c3ind,0,"3He2",'c'); /* cyclotron cooling */ linadd(CoolHeavy.cyntrn,0,"Cycn",'c'); /* changed from info to cooling june 25 95 to pick this up in primal.in * collisionally excited La cooling */ linadd(MAX2(0.,hydro.coola[0]),1216,"Cool",'i'); /* collisionally de-excited La heating */ linadd(MAX2(0.,-hydro.coola[0]),1216,"Heat",'i'); /* cooling due to n>2 Lyman lines */ linadd(MAX2(0.,hydro.restly[0]),960,"Crst",'i'); /* heating due to n>2 Lyman lines */ linadd(MAX2(0.,-hydro.restly[0]),960,"Hrst",'i'); /* cooling due to n>3 Balmer lines */ linadd(MAX2(0.,hydro.barest[0]),4861,"Crst",'i'); /* heating due to n>3 Balmer lines */ linadd(MAX2(0.,-hydro.barest[0]),4861,"Hrst",'i'); /* cooling due to higher Paschen lines */ linadd(MAX2(0.,hydro.hrest[0]),0,"Crst",'i'); /* heating due to higher Paschen lines */ linadd(MAX2(0.,-hydro.hrest[0]),0,"Hrst",'i'); /* H2 rotation lines from Lepp and Shull ApJ 270, 578. */ linadd(CoolHeavy.h2line,2,"H2 l",'c'); /* H2 dissociation by H atoms (not e) */ linadd(CoolHeavy.h2dis,0,"H2dC",'c'); /* molecular hydrogen heating */ h2dish.h2dtot += h2dish.HeatH2Dish*(float)radius.dVeff; h2dish.h2dfrc = (float)MAX2(h2dish.HeatH2Dish/heat.htot,h2dish.h2dfrc); /* heating by H2 dissociation by Lyman continuum */ linadd(h2dish.HeatH2Dish,0,"H2dH",'h'); /* neg H ion free-bound emission */ linadd(CoolHeavy.hmfb,0,"H-FB",'c'); /* H+ + H => H2+ + photon continuum cooling */ linadd(h2pcng.H2PlsCool,0,"H2+ ",'c'); /* HeH+ formation cooling */ linadd(MAX2(3.27e-12+phycon.te*BOLTZMANN,0.)*h.hii*he.hei*1e-20+ (1.76e-11+phycon.te*BOLTZMANN)*h.hi*he.heii*1e-16,0,"HEH+",'i'); /* carbon monoxide heating */ codish.codtot += codish.CODissHeat*(float)radius.dVeff; codish.codfrc = (float)MAX2(codish.CODissHeat/heat.htot,codish.codfrc); /* carbon monoxide co photodissociation */ linadd(codish.CODissHeat,0,"COdh",'h'); /* cooling due to coll of vib rot levels */ linadd(covib.COVibCool,0,"CO C",'c'); /* He brems emission */ linadd(CoolHeavy.bhe,0,"HeFF",'c'); /* He recombination cooling */ linadd(CoolHeavy.herec,0,"HeFB",'c'); /* heavy element recombination cooling */ linadd(CoolHeavy.heavfb,0,"MeFB",'c'); /* metal brems emission */ linadd(CoolHeavy.heavff,0,"MeFF",'c'); /* total brems emission */ linadd(hydro.FreeFreeCool+CoolHeavy.bhe+CoolHeavy.heavff,0,"ToFF",'i'); /* He I triplet net collisional ionization cooling */ linadd(MAX2(0.,cionhe.He3IonCool),0,"He3I",'c'); /* He I triplet net 3-body heating */ linadd(MAX2(0.,-cionhe.He3IonCool),0,"He3b",'h'); /* He I 584, collisional excitation cooling */ linadd(MAX2(0.,CoolHeavy.c584),584,"HeIC",'c'); /* He I 584, collisional de-excitation heating */ linadd(MAX2(0.,-CoolHeavy.c584),584,"HeIH",'h'); /* total 584 emergent from cloud in inner direction * now compute emergent Helium emergent flux */ if( sphere.lgSphere ) { e911 = 0.; e584 = 0.; e626 = 0.; e228 = 0.; e504 = 0.; } else { if( strncmp( dtrans.chDffTrns , "OU" , 2) == 0 ) { eff = 0.5; } else { eff = 1.; } e504 = heots.esc504/2.*3.95e-11*eff; /* this is inward looking tau but no fac of two in original defn */ e584 = heots.esc584/2.*3.41e-11/2.; /* this looked to illuminated face and does not need extra factor */ e626 = heots.esc626*3.18e-11; e911 = heots.esc911*2.18e-11*eff; e228 = heots.esc228*8.72e-11*eff; } /* He I rec to ground escaping cloud */ linadd(e504,504,"He I",'i'); /* He I Ly alpha escaping cloud * call linadd( e584 , 584, 'esc ' ,'i') * these lines added to outlin in metdif - following must be false */ lindst(e584,584,"esc ",iphe1lCom.iphe1l[0][1],'i',FALSE); /* He I 2^3S-ground emission escaping lluminated face of cloud * call linadd( e626 , 626 , 'esc ' ,'i') * these lines added to outlin in metdif - following must be false */ lindst(e626,626,"esc ",ip626.ipHe23sGrnd,'i',FALSE); /* He I 4471, */ if( phycon.te <= 1e4 ) { em = 2.42e-22/(phycon.te/TePowers.te10); } else { em = 8.79e-22/phycon.te/TePowers.te03/TePowers.te01; } /* >>chng 97 jan 4, label was HeI, changed to Ca B for symmetry with others * He I 4471 recombination only, fit to Brocklehurst '72 */ rec = phycon.eden*he.heii*em; PntForLine(4471.,"He 1",&ipnt); lindst(rec,4471,"Ca B",ipnt,'i',TRUE); /* He I total 5876, from n-level atom * these lines added to outlin in metdif - following must be false */ lindst(he3lines.p5876,5876,"TOTL",he3lines.ipHe3l[1],'i',FALSE); /* inward part of 5876 */ linadd(he3lines.p5876*he3tau[IPT5876-1].FracInwd, 5876,"Inwd",'i'); /* He I 5876 REC, simple fit to Brocklehurst */ linadd(phycon.eden*he.heii*4.22e-21/(phycon.te*TePowers.te10), 5876,"Ca B",'r'); /* He I 6678 REC, simple fit to Brocklehurst */ rec = phycon.eden*he.heii*1.32e-21/(phycon.te*TePowers.te10*TePowers.te01); PntForLine(6678.,"He 1",&ipnt); lindst(rec,6678,"Ca B",ipnt,'i',TRUE); /* 10830 from n-level atom * these lines added to outlin in metdif - following must be false */ lindst(he3lines.p10830,10830,"TOTL",he3lines.ipHe3l[0],'i',FALSE); f10830 = 1.02e7*he3tau[IPT10830-1].Pesc/(1.02e7* he3tau[IPT10830-1].Pesc + 2.51e-9*TePowers.sqrte* phycon.eden); /* He I 10830 produced by radiative recombination */ rec = he.heii*phycon.eden*5.13e-22/(TePowers.te70*TePowers.te10* TePowers.te01*TePowers.te01)*f10830; linadd(rec,10830,"reco",'i'); /* following expression from Clegg 1987 */ pop = he.heii*5.79e-6*(pow(phycon.te/1e4,-1.18))/(1. + 3110./phycon.eden* (pow(phycon.te/1e4,-0.51))); col = pop*he3c.c2s2p*1.84e-12*f10830; /* collisionally excited 10830 estimated from Clegg 1987 (not model atom) */ linadd(col,10830,"coll",'i'); /* inward escaping HeI 10830 */ linadd(heinwd.he10in,10830,"Inwd",'i'); /* He I 3889 from n-level atom * these lines added to outlin in metdif - following must be false */ lindst(he3lines.p3889,3889,"TOTL",he3lines.ipHe3l[3],'i',FALSE); /* He I 7065 from n-level atom * these lines added to outlin in metdif - following must be false */ lindst(he3lines.p7065,7065,"TOTL",he3lines.ipHe3l[2],'i',FALSE); /* total collisional He I triplet line cooling, from n-level atom */ linadd(MAX2(0.,he3lines.he3clx),0,"CcHE",'c'); /* total collisional de-exec He I heating, from n-level atom */ linadd(MAX2(0.,-he3lines.he3clx),0,"ChHE",'h'); em = phycon.eden*he.heii; cb206 = 1.064e-22/TePowers.te70/TePowers.te10/TePowers.te03/TePowers.te03* em; /* case B He I 2.06 micron */ linadd(cb206,2,"Ca B",'i'); /* He I 2.06 micron from model atom, all physics in */ totl = he.heii*phe1lv.he12p*t206.Aul*t206.Pesc* 9.66e-13; PntForLine(20600.,"He 1",&ipnt); lindst(totl,2,"TOTL",ipnt,'i',TRUE); /* Case B He I 5016 */ cb5016 = 5.43e-23/TePowers.te70/TePowers.te10*em; PntForLine(5016.,"He 1",&ipnt); lindst(cb5016,5016,"Ca B",ipnt,'i',TRUE); /* He II Lyman continuum */ linadd(e228,228,"HeII",'i'); /* He II Balmer continuum escaping from cloud */ linadd(e911,911,"He2C",'i'); /* ==================================================== * end helium * ====================================================*/ if( LineSave.ipass <= 0 ) { /* on first pass set flag saying that case b results are ok, * will be set false if we ever overstep bounds of table */ CaseBHS.lgHCaseBOK = TRUE; CaseBHS.lgHeCaseBOK = TRUE; } if( trace.lgTrace ) { fprintf( ioQQQ, " LineSet1 returns\n" ); } # ifdef DEBUG_FUN fputs( " <->LineSet1()\n", debug_fp ); # endif return; }