/*RTDiffuse fill in DiffLocal and ThrowOut arrays, with diffuse emission for this zone */ #include #include "cddefines.h" #include "ip626.h" #include "showme.h" #include "taulines.h" #include "physconst.h" #include "nhe1lvl.h" #include "grainvar.h" #include "opacpoint.h" #include "hydrogenic.h" #include "opac.h" #include "tff.h" #include "ffsum.h" #include "maxcor.h" #include "trace.h" #include "gaunts.h" #include "iphmin.h" #include "em2nu.h" #include "he1blt.h" #include "he.h" #include "tepowers.h" #include "rfield.h" #include "hhe2phtots.h" #include "h.h" #include "phycon.h" #include "nsshells.h" #include "nh.h" #include "hcaseb.h" #include "nhe1.h" #include "hstat.h" #include "he1stat.h" #include "radius.h" #include "he1nionryd.h" #include "difher.h" #include "iphe1l.h" #include "outer.h" #include "phe1lv.h" #include "exptau.h" #include "grainemis.h" #include "diffheav.h" #include "ionfracs.h" #include "hydrophocs.h" #include "sexp.h" #include "ionrange.h" #include "heavy.h" #include "pbowen.h" #include "fekems.h" #include "dtrans.h" #include "heots.h" #include "sumoutinc.h" #include "he3lines.h" #include "plasnu.h" #include "dumpline.h" #include "pop371.h" #include "rt.h" void RTDiffuse(void) { long int i, ipPlasmaFreq, ip, ipHi, ipLo, ipZ, ipla, nelem, nion, limit, n, nd; double EdenAbund, difflya, esc, ots, xInWrd, arg, bfac, bolt, bolt2, cor, csThresh, dwid, f, fac, facff, fache2, factor, gamma, gion, gn, photon, sum1ex, sum1gr, sumh1[LIMELM], sumhex, sumtri; float Dilution , fach; /* *total local two photon emission */ float Diff2Pht[NCELL]; # ifdef DEBUG_FUN fputs( "<+>RTDiffuse()\n", debug_fp ); # endif /* this routine evaluates the local diffuse fields * it fills in the vector DiffLocal */ for( i=0; i < NCELL; i++ ) { rfield.ThrowOut[i] = 0.; rfield.ThrowOutNoInter[i] = 0.; rfield.DiffLocal[i] = 0.; Diff2Pht[i] = 0.; } /* hydrogen continuous recombination continuum (free-bound) emission * following used to adjust where within WIDFLX exp is evaluated */ dwid = 0.2; gion = 1.; sumhex = 0.; /* loop over all elements that have hydrogenic species*/ for( ipZ=0; ipZ < LIMELM; ipZ++ ) { sumh1[ipZ] = 0.; EdenAbund = phycon.eden*xIonFracs[ipZ+2][ipZ]; if( IonRange.IonHigh[ipZ] == ipZ + 2 ) { /* loop over all levels with explicit opacities */ limit = MIN2(NHPLPHOT-1,nhlevel); for( n=IP1S; n <= limit; n++ ) { /* reset this counter at 2s, lowest excited state, will use below * to confirm case B sum */ if( n == IP2S ) { /* sumh1 is sum rec to ground, sumhex will be sum to all excited */ sumh1[ipZ] = sumhex; sumhex = 0.; } gn = hstat.HStatWght[n]; gamma = 2.0618684e11*gn/gion/phycon.te/TePowers.sqrte; /* following is boltz fac for second cell, check not zero */ bolt2 = 1.57890e5*rfield.anu[nh.ipHn[ipZ][n]]/phycon.te; /* pick highest continuum to fill in */ if( n == IP1S ) { /* for Lyman continuum, go up to highest defined boltzmann factor */ ipHi = maxcor.ipMaxBolt; } else { /* for excited states go either up to Lyman continuum or highest * defined Boltzmann factor */ ipHi = MIN2(maxcor.ipMaxBolt,nh.ipHn[ipZ][IP1S]); } /* only escaping part of recombinations are thrown into ThrowOut here, * added to outcon. the lost photons are never put back in.*/ if( hydro.hcbolt[ipZ][n] > 0. ) { for( i=nh.ipHn[ipZ][n]-1; i < ipHi; i++ ) { /* flux is in photons per sec per ryd * ContBoltz is ratio of Boltz fac for each freq */ photon = gamma*rfield.ContBoltz[i]/hydro.hcbolt[ipZ][n]* rfield.widflx[i]*opac.OpacStack[i-nh.ipHn[ipZ][n]+ OpacPoint.ipHOpac[ipZ][n]]*rfield.anu2[i]; sumhex += photon; rfield.DiffLocal[i] += (float)(photon*EdenAbund); rfield.ThrowOut[i] += (float)(photon*EdenAbund*hydro.hrec[ipZ][n][ipRecEsc]/**/); } } else if( bolt2 < 3. ) { /* boltzmann factor is zero, but still emission within continuum */ for( i=nh.ipHn[ipZ][n]-1; i < ipHi; i++ ) { arg = MAX2(0.,1.57890e5* (rfield.anu[i]-hydro.HLevNIonRyd[ipZ][n]+ rfield.widflx[i]*dwid)/phycon.te); /* flux is in photons per sec per ryd */ photon = gamma*sexp(arg)*rfield.widflx[i]* opac.OpacStack[i-nh.ipHn[ipZ][n]+OpacPoint.ipHOpac[ipZ][n]]* rfield.anu2[i]; sumhex += photon; rfield.DiffLocal[i] += (float)(photon*EdenAbund); rfield.ThrowOut[i] += (float)(photon*EdenAbund*hydro.hrec[ipZ][n][ipRecEsc]); } } } for( n=NHPLPHOT; n <= nhlevel; n++ ) { /* reset this counter at 2s, lowest excited state, will use below * to confirm case B sum */ if( n == IP2S ) { /* sumh1 is sum rec to ground, sumhex will be sum to all excited */ sumh1[ipZ] = sumhex; sumhex = 0.; } gn = hstat.HStatWght[n]; csThresh = HydroPhoCS.STH[n]*rfield.anu3[nh.ipHn[ipZ][n]-1]/(POW2(ipZ+1.)); gamma = 2.0618684e11*gn/gion/phycon.te/TePowers.sqrte* csThresh; /* following is boltz fac for second cell, check not zero */ bolt2 = 1.57890e5*rfield.anu[nh.ipHn[ipZ][n]]/phycon.te; /* pick highest continuum to fill in */ if( n == IP1S ) { /* for Lyman continuum, go up to highest defined boltzmann factor */ ipHi = maxcor.ipMaxBolt; } else { /* for excited states go either up to Lyman continuum or highest * defined Boltzmann factor */ ipHi = MIN2(maxcor.ipMaxBolt,nh.ipHn[ipZ][IP1S]); } if( hydro.hcbolt[ipZ][n] > 0. ) { for( i=nh.ipHn[ipZ][n]-1; i < ipHi; i++ ) { /* flux is in photons per sec per ryd * ContBoltz is ratio of Boltz fac for each freq */ photon = gamma*rfield.ContBoltz[i]/hydro.hcbolt[ipZ][n]* rfield.widflx[i]/rfield.anu[i]; sumhex += photon; rfield.DiffLocal[i] += (float)(photon*EdenAbund); rfield.ThrowOut[i] += (float)(photon*EdenAbund*hydro.hrec[ipZ][n][ipRecEsc]); } } else if( bolt2 < 3. ) { /* boltzmann factor is zero, but still emission within continuum */ for( i=nh.ipHn[ipZ][n]-1; i < ipHi; i++ ) { arg = MAX2(0.,1.57890e5* (rfield.anu[i]-hydro.HLevNIonRyd[ipZ][n]+ rfield.widflx[i]*dwid)/phycon.te); /* flux is in photons per sec per ryd */ photon = gamma*sexp(arg)*rfield.widflx[i]/ rfield.anu[i]; sumhex += photon; rfield.DiffLocal[i] += (float)(photon*EdenAbund); rfield.ThrowOut[i] += (float)(photon*EdenAbund*hydro.hrec[ipZ][n][ipRecEsc]); } } } /* following is check on self-consistency */ difher.CaseBCheck[ipZ] = (float)(MAX2(difher.CaseBCheck[ipZ], sumhex/hcaseb.HRecRad[ipZ])); } } /* at this stage the array has only the hydrogenic recombination continua */ /* neutral helium triplets hei */ EdenAbund = phycon.eden*xIonFracs[2][ipHE]; sumtri = 0.; if( EdenAbund > 0. ) { gion = 2.; gn = 3.; gamma = 2.0618684e11*gn/gion/phycon.te/TePowers.sqrte; bolt2 = 1.57890e5*rfield.anu[he.nhei3-1]/phycon.te; if( rfield.ContBoltz[he.nhei3-1] > 0. ) { for( i=he.nhei3-1; i < nhe1Com.nhe1[0]; i++ ) { photon = gamma*rfield.ContBoltz[i]/rfield.ContBoltz[he.nhei3-1]* rfield.widflx[i]*opac.OpacStack[i-he.nhei3+OpacPoint.ioptri]* rfield.anu2[i]; sumtri += photon; rfield.DiffLocal[i] += (float)(photon*EdenAbund); rfield.ThrowOut[i] += (float)(photon*EdenAbund); } } else if( bolt2 < 4. ) { for( i=he.nhei3-1; i < nhe1Com.nhe1[0]; i++ ) { arg = MAX2(0.,1.57890e5*(rfield.anu[i]-rfield.anu[he.nhei3-1])/ phycon.te); photon = gamma*sexp(arg)*rfield.widflx[i]*opac.OpacStack[i-he.nhei3+OpacPoint.ioptri]* rfield.anu2[i]; sumtri += photon; rfield.DiffLocal[i] += (float)(photon*EdenAbund); rfield.ThrowOut[i] += (float)(photon*EdenAbund); } } } /* singlet neutral helium free-bound recombination continua * * HeI ground state from Milne relation */ EdenAbund = phycon.eden*xIonFracs[2][ipHE]; sum1gr = 0.; if( EdenAbund > 0. ) { gion = 2.; gn = 1.; gamma = 2.0618684e11*gn/gion/phycon.te/TePowers.sqrte; cor = he1blt.he1cbt[0]; if( cor > 0. ) { for( i=nhe1Com.nhe1[0]-1; i < maxcor.ipMaxBolt; i++ ) { /* flux is in photons per sec per ryd * ContBoltz is ratio of Boltz fac for each freq */ photon = gamma*rfield.ContBoltz[i]/cor*rfield.widflx[i]* opac.OpacStack[i-nhe1Com.nhe1[0]+OpacPoint.iophe1[0]]* rfield.anu2[i]; sum1gr += photon; rfield.DiffLocal[i] += (float)(photon*EdenAbund); rfield.ThrowOut[i] += (float)(photon*EdenAbund*phe1lv.he1rec[2][0]); } } else { for( i=nhe1Com.nhe1[0]-1; i < maxcor.ipMaxBolt; i++ ) { arg = MAX2(0.,1.57890e5*(rfield.anu[i]-He1NIonRyd.He1IonRyd[0])/ phycon.te); /* flux is in photons per sec per ryd */ photon = gamma*sexp(arg)*rfield.widflx[i]*opac.OpacStack[i-nhe1Com.nhe1[0]+ OpacPoint.iophe1[0]]*rfield.anu2[i]; sum1gr += photon; rfield.DiffLocal[i] += (float)(photon*EdenAbund); rfield.ThrowOut[i] += (float)(photon*EdenAbund*phe1lv.he1rec[2][0]); } } } /* do recombination from He+ */ sum1ex = 0.; EdenAbund = phycon.eden*xIonFracs[2][ipHE]; if( EdenAbund > 0. ) { gion = 2.; for( n=1; n < NHE1LVL; n++ ) { bolt2 = 1.57890e5*rfield.anu[nhe1Com.nhe1[n]-1]/phycon.te; /* this is trick to do 2p, 2s correctly */ gn = MAX2(4.,he1statCOM.he1stat[n]); gamma = 2.0618684e11*gn/gion/phycon.te/TePowers.sqrte; /* Milne relation for excited states */ limit = MIN2(nhe1Com.nhe1[0],rfield.nflux); if( rfield.ContBoltz[nhe1Com.nhe1[n]-1] > 0. ) { for( i=nhe1Com.nhe1[n]-1; i < limit; i++ ) { photon = gamma*rfield.ContBoltz[i]/rfield.ContBoltz[nhe1Com.nhe1[n]- 1]*rfield.widflx[i]*opac.OpacStack[i-nhe1Com.nhe1[n]+ OpacPoint.iophe1[n]]*rfield.anu2[i]; sum1ex += photon; rfield.DiffLocal[i] += (float)(photon*EdenAbund); rfield.ThrowOut[i] += (float)(photon*EdenAbund*phe1lv.he1rec[2][n]); } } else if( bolt2 < 3. ) { for( i=nhe1Com.nhe1[n]-1; i < limit; i++ ) { arg = MAX2(0.,1.57890e5*(rfield.anu[i]-He1NIonRyd.He1IonRyd[n])/ phycon.te); photon = gamma*sexp(arg)*rfield.widflx[i]* opac.OpacStack[i-nhe1Com.nhe1[n]+OpacPoint.iophe1[n]]* rfield.anu2[i]; sum1ex += photon; rfield.DiffLocal[i] += (float)(photon*EdenAbund); rfield.ThrowOut[i] += (float)(photon*EdenAbund*phe1lv.he1rec[2][n]); } } } } /* trace diffuse triggers this print */ if( trace.lgTrace && trace.lgTrDiff ) { fprintf( ioQQQ, " RTDiffuse: Hgrd;%10.2e Ex;%10.2e He3;%10.2e He1g;%10.2e Ex;%10.2e He2g;%10.2e Ex;\n", sumh1[0], sumhex, sumtri, sum1gr, sum1ex, sumh1[1] ); } /* hydrogen two photon emission * in all following the factor of two is because a single * decay produces two photons */ EdenAbund = 2.*hydro.hn[0][IP2S]*xIonFracs[2][ipH]; hhe2phtots.H12phtOTS = 0.; /* upper limit to H 2-phot is energy of Lya */ limit = HydroLines[0][IP2P][IP1S].ipCont; for( i=0; i < limit; i++ ) { fach = em2nu.Hy2nu[i]*(float)EdenAbund; rfield.DiffLocal[i] += fach; Diff2Pht[i] = fach; /* this is escaping part, * as determined from optical depth to illuminated face */ rfield.ThrowOut[i] += fach*exptau.ExpmTau[i]; /* not escaping part will be counted as ots at a single energy */ hhe2phtots.H12phtOTS += fach*(1.f - exptau.ExpmTau[i]); } { enum {DEBUG=FALSE}; if( DEBUG) { fprintf(ioQQQ," RTDiffuse expmtau %e %e\n" , exptau.ExpmTau[HydroLines[0][IP2P][IP1S].ipCont - 30 ], hhe2phtots.H12phtOTS); } } /* HeI two photon emission, first part that cannot ionize hydrogen */ EdenAbund = 2.*phe1lv.he12s*xIonFracs[2][ipHE]; hhe2phtots.He12phtOTS = 0.; /* following also protects from element being turned off */ if( EdenAbund > 0. ) { for( i=0; i < nh.ipHn[0][IP1S]; i++ ) { rfield.DiffLocal[i] += (float)(em2nu.he12nu[i]*EdenAbund); Diff2Pht[i] += (float)(em2nu.he12nu[i]*EdenAbund); rfield.ThrowOut[i] += (float)(em2nu.he12nu[i]*EdenAbund*exptau.ExpmTau[i]); } /* HeI two photon emission, energies that can photoionize hydrogen */ for( i=nh.ipHn[0][IP1S]; i < iphe1lCom.iphe1l[0][1]; i++ ) { rfield.DiffLocal[i] += (float)(em2nu.he12nu[i]*EdenAbund); Diff2Pht[i] += (float)(em2nu.he12nu[i]*EdenAbund); rfield.ThrowOut[i] += (float)(em2nu.he12nu[i]*EdenAbund*exptau.ExpmTau[i]); hhe2phtots.He12phtOTS += (float)(em2nu.he12nu[i]*EdenAbund*(1. - exptau.ExpmTau[i])); } } /* HeII two photon emission */ EdenAbund = 2.*hydro.hn[1][IP2S]*xIonFracs[3][ipHE]; hhe2phtots.He22phtOTS = 0.; if( EdenAbund > 0. ) { /* add HeII 2-nu for energies that cannot ionize H Lyman continuum */ for( i=0; i < nh.ipHn[0][IP1S]; i++ ) { rfield.DiffLocal[i] += (float)(em2nu.he22nu[i]*EdenAbund); Diff2Pht[i] += (float)(em2nu.he22nu[i]*EdenAbund); rfield.ThrowOut[i] += (float)(em2nu.he22nu[i]*EdenAbund*exptau.ExpmTau[i]); } /* add HeII 2-nu for energies that can ionize H Lyman continuum */ limit = HydroLines[1][IP2P][IP1S].ipCont; for( i=nh.ipHn[0][IP1S]; i < limit; i++ ) { rfield.DiffLocal[i] += (float)(em2nu.he22nu[i]*EdenAbund); Diff2Pht[i] += (float)(em2nu.he22nu[i]*EdenAbund); rfield.ThrowOut[i] += (float)(em2nu.he22nu[i]*EdenAbund*exptau.ExpmTau[i]); hhe2phtots.He22phtOTS += (float)(em2nu.he22nu[i]*EdenAbund*(1. - exptau.ExpmTau[i])); } } /* H, He, heavy element, brems free-free */ fach = (float)(phycon.eden*1.032e-11*(xIonFracs[2][ipH] + xIonFracs[2][ipHE])/ TePowers.sqrte); fache2 = phycon.eden*1.032e-11*4.*xIonFracs[3][ipHE]/TePowers.sqrte; facff = phycon.eden*1.032e-11*ffsum.EdenFFSum*5./TePowers.sqrte; /* ntff is pointer to energy where gas goes * optically thin to brems */ for( i=0; i < (tffCom.ntff - 1); i++ ) { fac = (fach*gaunts.gff[i]*(1. + opac.OpacStack[i-1+OpacPoint.iphmra]* hydro.hn[0][IP1S]) + fache2*gaunts.gffhe2[i] + facff*gaunts.gff[i])* rfield.widflx[i]*rfield.ContBoltz[i]/rfield.anu[i]; rfield.ThrowOutNoInter[i] += (float)fac; rfield.DiffLocal[i] += (float)fac; /* do not add to outward beam since self absorbed * >>chng 96 feb 27, put back into outward beam since do not integrate * over it anyway. */ /* >>chng 99 May 28, take back out of beam since DO integrate over it * in very dense BLR clouds */ /*rfield.ThrowOut[i] += (float)fac;*/ } for( i=tffCom.ntff-1; i < maxcor.ipMaxBolt; i++ ) { fac = (fach*gaunts.gff[i]*(1. + opac.OpacStack[i-1+OpacPoint.iphmra]* hydro.hn[0][IP1S]) + fache2*gaunts.gffhe2[i] + facff*gaunts.gff[i])* rfield.widflx[i]*rfield.ContBoltz[i]/rfield.anu[i]; rfield.DiffLocal[i] += (float)fac; rfield.ThrowOut[i] += (float)fac; } /* grain dust emission */ if( GrainVar.lgDustOn ) { /* add grain dust emission in IR */ for( i=0; i < rfield.nflux; i++ ) { GrainEmis.GrainEmission[i] = 0.; } for( nd=0; nd < NDUST; nd++ ) { if( GrainVar.lgDustOn1[nd] ) { factor = 2.*PI4*POW2(FR1RYD/SPEEDLIGHT)*FR1RYD; bolt = TE1RYD/GrainVar.tedust[nd]; f = 1.; i = 0; /* >>chng 96 apr 26 upper limt chng from 15 to 75 Peter van Hoof */ while( i < rfield.nflux && f < 75. ) { f = bolt*rfield.anu[i]; if( f > 75. ) { bfac = 1e33; } else if( f > .01 ) { bfac = exp(f) - 1.; } else { bfac = f; } /* >>chng 96 july 13, save grain emission array * >>chng 97 feb 28, GrainEmission was set to last evaluated grain, * zeroed out above, and now increment grain instead of setting it */ cor = factor*(GrainVar.dstab1[nd][i]*phycon.hden* GrainVar.dstAbund[nd])*rfield.anu2[i]*rfield.widflx[i]/ bfac; GrainEmis.GrainEmission[i] += (float)cor; rfield.DiffLocal[i] += (float)cor; rfield.ThrowOut[i] += (float)cor; ++i; } } } } /* hminus emission */ fac = phycon.eden*(double)h.hi; gn = 1.; gion = 2.; gamma = 2.0618684e11*gn/gion/phycon.te/TePowers.sqrte; cor = rfield.ContBoltz[iphminCom.iphmin-1]; limit = MIN2(nhe1Com.nhe1[0],rfield.nflux); if( cor > 0. ) { for( i=iphminCom.iphmin-1; i < limit; i++ ) { /* flux is in photons per sec per ryd * ContBoltz is ratio of Boltz fac for each freq */ factor = gamma*rfield.ContBoltz[i]/cor*rfield.widflx[i]* opac.OpacStack[i-iphminCom.iphmin+OpacPoint.iphmop]* rfield.anu2[i]*fac; rfield.DiffLocal[i] += (float)factor; rfield.ThrowOut[i] += (float)factor; } } else { for( i=iphminCom.iphmin-1; i < limit; i++ ) { arg = MAX2(0.,1.57890e5*(rfield.anu[i]-0.05544)/phycon.te); /* flux is in photons per sec per ryd */ factor = gamma*sexp(arg)*rfield.widflx[i]* opac.OpacStack[i-iphminCom.iphmin+OpacPoint.iphmop]* rfield.anu2[i]*fac; rfield.DiffLocal[i] += (float)factor; rfield.ThrowOut[i] += (float)factor; } } /* this is a unit of energy that will be passed through the code as a test * that all integrations are carried out. A similar test is set in lineset1 * and verified in PrtFinal. The opacity at this cell is zero so only * geometrical dilution will affect the integral * Radius is currently outer edge of zone, so radius-drad/2 is radius * of center of zone */ /* this dilution is needed to conserve volume in spherical models. tests such * as parispn.in will fault if this is removed */ Dilution = (float)POW2( radius.rinner / (radius.Radius-radius.drad/2.) );/**/ rfield.DiffLocal[rfield.nflux] = 1.e-10f * Dilution /* */ ; rfield.ThrowOut[rfield.nflux] = 1.e-10f * Dilution /* */; /* opacity should be zero at this energy */ assert( opac.opac[rfield.nflux] == 0. ); /* many tmn added to conserve energy in very high Z models * rerun highZ qso model if any tmn ever deleted here * * tmn set in StartZone and includes both opacity and dilution * * use duffuse lines and continuum to add flux to outward beam * * NB!!! this routine adds flux to the outward beam * it can only be called once per zone * * covrt is radiative transfer covering factor * covrt = 0 for open geometry * covrt = 1 for closed geometry * * outwrd is fraction that goes outward, of the escaping line radiation * default covrt is zero for open geometry, is 1 for closed * outwrd is 1/2 for opten, 1 for closed */ /*outwrd = (1. + cover.covrt)/2.;*/ /* radius is the outer edge of the current zone * dVolOutwrd = outwrd * dreff * >>chng 95 nov 16 * dVolOutwrd = outwrd * dreff * ( radius/(radius+drad) )**2 * dreff in upper since fill fac enters, but not in lower * >>chng 96 may 8, moved to zonsrt * dVolOutwrd = outwrd * radius * dreff / (radius + 2.*drad) */ /* save contents of outcon for checking how much was added */ for( i=0; i < rfield.nflux; i++ ) { SumOutInc.SavOutCon[i] = rfield.outcon[i] + rfield.outlin[i]; } /* >>chng 96 may 8, moved to zonsrt * dVolReflec = (1.-outwrd) * dReff * r1r0sq * dVolOutwrd = outwrd * radius * dReff / (radius + 2.*drad) * * continuum interaction arrays that occur here: * * ContNoInter - not counted in internal energy budget * of the cloud, used to get pretty plots afterwards * * outcon, outward only continuum, inc in heating rates * otscon, otslin, loc ots rates, inc in heating */ /* add HI hydrogenic lines to outward and reflected beams */ for( ipZ=0; ipZ < LIMELM; ipZ++ ) { if( (IonRange.IonHigh[ipZ] == ipZ + 2) ) { for( ipLo=IP1S; ipLo <= (nhlevel - 1); ipLo++ ) { /* must not include 2s-1s two photon in this since is actually a continuum */ for( ipHi=MAX2(IP2P,ipLo+1) ; ipHi <= nhlevel; ipHi++ ) { /* pointer to line energy in continuum array */ ip = HydroLines[ipZ][ipHi][ipLo].ipCont-1; /* number of photons in the line has not been defined up til now, * do so now. this is redone in lines. */ HydroLines[ipZ][ipHi][ipLo].phots = HydroLines[ipZ][ipHi][ipLo].Aul* HydroLines[ipZ][ipHi][ipLo].PopHi* HydroLines[ipZ][ipHi][ipLo].Pesc* xIonFracs[ipZ+2][ipZ]; /* inward fraction of line */ xInWrd = HydroLines[ipZ][ipHi][ipLo].phots* HydroLines[ipZ][ipHi][ipLo].FracInwd; /* reflected part of line */ rfield.reflin[ip] += (float)(xInWrd*radius.BeamInIn); /* inward beam that goes out since sphere set */ /* in all this the ColOvTot term has been commented out, * since this is not meaningful for something like the hydrogen atom, * where most forms by recombination */ rfield.outlin[ip] += (float)(xInWrd*radius.BeamInOut*opac.tmn[ip]/* HydroLines[ipZ][ipHi][ipLo].ColOvTot*/); /* outward part of line */ rfield.outlin[ip] += (float)(HydroLines[ipZ][ipHi][ipLo].phots* (1. - HydroLines[ipZ][ipHi][ipLo].FracInwd)* radius.BeamOutOut*opac.tmn[ip]/* HydroLines[ipZ][ipHi][ipLo].ColOvTot*/); } } } } /* HeI He I He I He I He I He I He I He I * outward part of escaping He I Ly alpha hei lya */ rfield.reflin[iphe1lCom.iphe1l[0][1]-1] += (float)radius.dVolReflec*heots.esc584; rfield.outlin[iphe1lCom.iphe1l[0][1]-1] += (float)radius.dVolOutwrd*heots.esc584; /* HeI higher lyman lines */ for( ipHi=2; ipHi < (NHE1LVL - 1); ipHi++ ) { rfield.outlin[iphe1lCom.iphe1l[0][ipHi]-1] += he.heii*phe1lv.he1n[ipHi]* phe1lv.he1mis[0][ipHi]*(float)radius.dVolOutwrd; rfield.reflin[iphe1lCom.iphe1l[0][ipHi]-1] += he.heii*phe1lv.he1n[ipHi]* phe1lv.he1mis[0][ipHi]*(float)radius.dVolReflec; } /* helium 1 triplet lines */ rfield.outlin[he3lines.ipHe3l[0]-1] += (float)(he3lines.p10830/1.84e-12* radius.dVolOutwrd); rfield.outlin[he3lines.ipHe3l[1]-1] += (float)(he3lines.p5876/3.39e-12* radius.dVolOutwrd); rfield.outlin[he3lines.ipHe3l[2]-1] += (float)(he3lines.p7065/2.82e-12* radius.dVolOutwrd); rfield.outlin[he3lines.ipHe3l[3]-1] += (float)(he3lines.p3889/5.12e-12* radius.dVolOutwrd); rfield.reflin[he3lines.ipHe3l[0]-1] += (float)(he3lines.p10830/1.84e-12* radius.dVolReflec); rfield.reflin[he3lines.ipHe3l[1]-1] += (float)(he3lines.p5876/3.39e-12* radius.dVolReflec); rfield.reflin[he3lines.ipHe3l[2]-1] += (float)(he3lines.p7065/2.82e-12* radius.dVolReflec); rfield.reflin[he3lines.ipHe3l[3]-1] += (float)(he3lines.p3889/5.12e-12* radius.dVolReflec); /* He I ionized helium singlet lines */ for( ipLo=1; ipLo < (NHE1LVL - 2); ipLo++ ) { for( ipHi=ipLo + 1; ipHi <= (NHE1LVL - 1); ipHi++ ) { rfield.outlin[iphe1lCom.iphe1l[ipLo][ipHi]-1] += he.heii* phe1lv.he1n[ipHi]*phe1lv.he1mis[ipLo][ipHi]*(float)radius.dVolOutwrd; rfield.reflin[iphe1lCom.iphe1l[ipLo][ipHi]-1] += he.heii* phe1lv.he1n[ipHi]*phe1lv.he1mis[ipLo][ipHi]*(float)radius.dVolReflec; } } rfield.outlin[ip626.ipHe23sGrnd-1] += (float)radius.dVolOutwrd*heots.out626; rfield.reflin[ip626.ipHe23sGrnd-1] += (float)radius.dVolReflec*heots.out626; rfield.outlin[pbowen.ip660-1] += (float)radius.dVolOutwrd*heots.esc660; rfield.reflin[pbowen.ip660-1] += (float)radius.dVolReflec*heots.esc660; /* >>chng 96 nov 19, do not consider energies below plasma freq * ipPlasmaFreq points to lowest energy thin to brems and plas freq */ ipPlasmaFreq = MAX2(plasnu.ipPlasma,1)-1; /* default dtrans.chDffTrns is "OU2" */ /* case where OTS, set with diffuse ots command*/ if( strcmp(dtrans.chDffTrns,"OTS") == 0 ) { /********************************************************************* * * OTS continuum OTS OTS OTS * *********************************************************************/ /* pick up these continua here if we are not doing outward only. * note that the upper bounds to the sum is <= nflux, the continuum * only goes up through nflux points. In RTDiffuse the [nflux] * cell was set to 1e-10f, which will be propagated through the * code and verified in PrtFinal. */ for( i=ipPlasmaFreq; i <= rfield.nflux; i++ ) { /* continuum not used for heating, just pretty pictures * add it here - outward approximation will do it in outcon */ /* rfield.ContNoInter is not in energy budget */ /* radius.dVolOutwrd has geo part that is 1 for sphere, 1/2 for open */ /* opac.tmn is needed for FIR brems emission in dense BLR clouds - * each zone is vastly optically thick to 1 cm radiation */ rfield.ConOutNoInter[i] += rfield.ThrowOut[i]*(float)radius.dVolOutwrd/* *opac.tmn[i]*/; /* reflected continuum, geo part of dVolReflec is 1/2 for open, 0 for sphere */ rfield.ConRefNoInter[i] += rfield.ThrowOut[i]*(float)radius.dVolReflec; } } else if( strncmp( dtrans.chDffTrns , "OU" , 2) == 0 ) { /********************************************************************* * * outward only continuum * *********************************************************************/ /* THIS IS THE DEFAULT option for outward only continuum - * only test first two char so this picks up all outward continua */ /* add rfield.ThrowOut continuum (set in RTDiffuse) to outcon, * lower limit of ipPlasmaFreq since continuum is zero below ipPlasmaFreq * due to plasma frequency * note that upper range of sum is <= nflux, which contains the unit * verification token */ for( i=ipPlasmaFreq; i <= rfield.nflux; i++ ) { /* outcon is the interactive continuum * tmn is attenuation across one zone * ThrowOut contains all radiation thrown into outward beam, * eval in RTDiffuse */ /* NB opac.tmn is needed for FIR brems emission in dense BLR clouds - * each zone is vastly optically thick to 1 cm radiation */ rfield.outcon[i] += rfield.ThrowOut[i]*(float)radius.dVolOutwrd/**opac.tmn[i]*/; } /* now save the outward continua to track changes */ for( i=0; i < rfield.nflux; i++ ) { /* >>chng 96 july 13, add extra grain emiss to save out, so it is * not counted in effective sum - benign source of radiation */ SumOutInc.SavOutCon[i] += GrainEmis.GrainEmission[i]* (float)radius.dVolOutwrd; } } else { /* dtrans was neither one nor other, this is impossible */ fprintf( ioQQQ, " chDffTrns=%3.3s in RTDiffuse. This is impossible.\n", dtrans.chDffTrns ); puts( "[Stop in RTDiffuse]" ); exit(1); } # if !defined(NDEBUG) /*begin sanity check for level 1 lines */ for( i=1; i <= nLevel1; i++ ) { if( TauLines[i].phots* TauLines[i].ColOvTot < 0. ) { fprintf( ioQQQ, " RTDiffuse finds insane NPohts*ColovTot for level 1 line, vals=%10.1e%10.1e\n", TauLines[i].phots, TauLines[i].ColOvTot ); DumpLine(&TauLines[i]); ShowMe(); puts( "[Stop in RTDiffuse]" ); exit(1); } } /*end sanity check */ # endif /* outward level 1 line photons, 0 is dummy line */ for( i=1; i <= nLevel1; i++ ) { /* pointer to line energy in continuum array */ ip = TauLines[i].ipCont-1; /* last factor does not accout for frac of lines pumped */ xInWrd = TauLines[i].phots*TauLines[i].FracInwd; rfield.reflin[ip] += (float)(xInWrd*radius.BeamInIn); /* inward beam that goes out since sphere set */ rfield.outlin[ip] += (float)(xInWrd*radius.BeamInOut*opac.tmn[ip]* TauLines[i].ColOvTot); rfield.outlin[ip] += (float)(TauLines[i].phots* (1. - TauLines[i].FracInwd)*radius.BeamOutOut* /*opac.tmn[ip]**/TauLines[i].ColOvTot); } # if !defined(NDEBUG) /*begin sanity check for level 2 lines */ for( i=0; i < nWindLine; i++ ) { if( TauLine2[i].phots*TauLine2[i].ColOvTot < 0. ) { fprintf( ioQQQ, " RTDiffuse finds insane NPohts*ColovTot for level 2 line, vals=%10.1e%10.1e\n", TauLine2[i].phots, TauLine2[i].ColOvTot ); DumpLine(&TauLine2[i]); ShowMe(); puts( "[Stop in RTDiffuse]" ); exit(1); } } /*end sanity check */ # endif /* outward level 2 line photons */ for( i=0; i < nWindLine; i++ ) { /* NB if this changes, then FeIIEmitOut in pop371 * probably needs to change too * pointer to line energy in continuum array */ ip = TauLine2[i].ipCont-1; /* last factor does not accout for frac of lines pumped */ xInWrd = TauLine2[i].phots*TauLine2[i].FracInwd; rfield.reflin[ip] += (float)(xInWrd*radius.BeamInIn); /* inward beam that goes out since sphere set */ rfield.outlin[ip] += (float)(xInWrd*radius.BeamInOut*opac.tmn[ip]* TauLine2[i].ColOvTot); rfield.outlin[ip] += (float)(TauLine2[i].phots* (1. - TauLine2[i].FracInwd)*radius.BeamOutOut* /*opac.tmn[ip]**/TauLine2[i].ColOvTot); /* >>chng 96 july 16, added ipLnColovTOT */ } /* do one of the FeII atoms */ if( FeII.lgFeIION ) { /* this will do the above, adding FeII to reflected, outward beam */ if( xIonFracs[2][25] > 0. && FeII.lgFeIION ) { FeIIEmitOut(radius.dVolOutwrd,radius.dVolReflec); } } /* add recombination continua for elements heavier than He */ /* nelem = 2 is Li */ for( nelem=2; nelem < LIMELM; nelem++ ) { /* do not include hydrogenic species in following */ for( nion=0; nion < nelem; nion++ ) { if( xIonFracs[nion+2][nelem] > 0. ) { long int ipNu, ns, nshell,igRec , igIon, iplow , iphi , ipop; ip = Heavy.ipHeavy[nion][nelem]-1; /* nflux was reset upward in ConvInitTemp to encompass all * possible line and continuum emission. this test should not * possibly fail. It could it the ionization were allowed to increase somehow. * This is important because the nflux cell in outcon is used to carry out the * unit integration, and if it gets clobbered by diffuse emission the code * will declare insanity in PrtComment */ assert( ip < rfield.nflux ); /* this is extinction to the illuminated face */ esc = exptau.ExpmTau[ip]; /* DiffusHeavy is rad rec to ground, with electron density, * it is evaluated in routine MakeRecomb */ EdenAbund = DiffHeav.DiffusHeavy[nion][nelem]*xIonFracs[nion+2][nelem]; rfield.outcon[ip] += (float)(EdenAbund*radius.dVolOutwrd*esc*opac.tmn[ip]); /* get shell number, stat weights for this species */ outer( nelem+1 , nelem+1-nion , &nshell, &igRec , &igIon ); gn = igRec; gion = igRec; ns = nsShellsCom.nsShells[nion][nelem]-1; assert( ns == (nshell-1) ); /* lower and upper energies, and offset for opacity stack */ iplow = OpacPoint.ipElement[0][ns][nion][nelem]-1; iphi = OpacPoint.ipElement[1][ns][nion][nelem]; iphi = MIN2( iphi , rfield.nflux ); ipop = OpacPoint.ipElement[2][ns][nion][nelem]; /* now convert ipop to the offset in the opacity stack from threshold */ ipop = ipop - iplow; bolt2 = 1.57890e5*rfield.anu[iplow]/phycon.te; gamma = 2.0618684e11*gn/gion/phycon.te/TePowers.sqrte*phycon.eden*xIonFracs[nion+2][nelem]; /* this is ground state continuum from stored opacities */ sum1ex = 0.; if( rfield.ContBoltz[iplow] > SMALLFLOAT ) { for( ipNu=iplow; ipNu < iphi; ++ipNu ) { photon = gamma*rfield.ContBoltz[ipNu]/rfield.ContBoltz[iplow]* rfield.widflx[ipNu]*opac.OpacStack[ipNu+ipop]*rfield.anu2[ipNu]; sum1ex += photon; rfield.ThrowOutNoInter[ipNu] += (float)photon; } } /* now do the recombination Lya */ ipla = Heavy.ipLyHeavy[nion][nelem]-1; esc = exptau.ExpmTau[ipla]; /* xLyaHeavy is set to a fraction of rad rec in MakeRecomb */ difflya = DiffHeav.xLyaHeavy[nion][nelem]*xIonFracs[nion+2][nelem]; assert( difflya >= 0. ); rfield.outlin[ipla] += (float)(difflya*radius.dVolOutwrd*opac.tmn[ipla]*esc); rfield.reflin[ipla] += (float)(difflya*radius.dVolReflec*opac.tmn[ipla]*esc); } } } /* Fe Auger K-alpha in outward only * first is rec plus hot auger * this is pointer to valence shell Lya * >>chng 97 apr 28, had been 26,20, so picked up vastly wrong energy for Ka * improve to real pointers for line * >>chng 96 nov 20, added tmn to two below, to prevent line from * overwhelming very thick models * >>chng 97 apr 28, put in reflected Ka, had been only outward * outlin(ip) = outlin(ip) + fekhot/1.11e-8 * dReff * */ rfield.outlin[fekems.ipkhot-1] += fekems.fekhot*(float)radius.dVolOutwrd* opac.tmn[fekems.ipkhot-1]; rfield.reflin[fekems.ipkhot-1] += fekems.fekhot*(float)radius.dVolReflec; /* >>chng 97 apr 28, put in reflected Ka, had been only outward * this is cold iron, one cell below hot * outlin(ip-1) = outlin(ip-1) +fekcld/1.03e-8*dReff* */ rfield.outlin[fekems.ipkcld-1] += fekems.fekcld*(float)radius.dVolOutwrd* opac.tmn[fekems.ipkcld-2]; rfield.reflin[fekems.ipkcld-1] += fekems.fekcld*(float)radius.dVolReflec; /********************************************************************** * * * now analyse what we have done to the outward beam * * * **********************************************************************/ /* now figure out how much we added to the beam * this is used for setting zone thickness - do not want to * add too much energy in outward beam */ SumOutInc.SumOutCon = 0.; SumOutInc.SumIncCon = 0.; ots = 0.; /* this will be photon energy where this peak occurs */ SumOutInc.SumOutMax = 0.; for( i=0; i < (nh.ipHn[0][IP1S] - 1); i++ ) { /* >>chng 96 jan 14, was photon number, now energy * >>chng 96 nov 19, save difference in SavOutCon array, for later punch */ SumOutInc.SavOutCon[i] = (float)(MAX2(0.,rfield.outlin[i]+rfield.outcon[i]- SumOutInc.SavOutCon[i])); esc = SumOutInc.SavOutCon[i]*opac.opac[i]*rfield.anu[i]; if( esc > ots ) { ots = esc; SumOutInc.SumOutMax = rfield.anu[i]; SumOutInc.ipSumOutMax = i+1; } /* the ratio of SumOutCon to SumIncCon is kept below 5% in nextdr */ SumOutInc.SumOutCon += (float)esc; /* >>chng 96 nov 18, add outcon, outlin to sum * SumIncCon = SumIncCon + (flux(i)+otscon(i)+otslin(i)) *opac(i)* */ SumOutInc.SumIncCon += (float)((rfield.flux[i] + rfield.otscon[i] + rfield.otslin[i] + rfield.outcon[i] + rfield.outlin[i])* opac.opac[i]*rfield.anu[i]); } SumOutInc.SumHOutCon = 0.; SumOutInc.SumHIncCon = 0.; /* this will be sum of H-ionizing continuum * >>chng 96 jan 14, was photon number, now energy */ for( i=nh.ipHn[ipH][IP1S]-1; i < rfield.nflux; i++ ) { /* esc = ( outlin(i)+outcon(i)-SavOutCon(i) )* opac(i)* anu(i) */ SumOutInc.SavOutCon[i] = (float)MAX2(0.,rfield.outlin[i]+rfield.outcon[i]- SumOutInc.SavOutCon[i]); esc = SumOutInc.SavOutCon[i]*opac.opac[i]*rfield.anu[i]; if( esc > ots ) { ots = esc; SumOutInc.SumOutMax = rfield.anu[i]; SumOutInc.ipSumOutMax = i+1; } SumOutInc.SumHOutCon += (float)esc; /* >>chng 96 nov 18, add outcon, outlin to sum */ SumOutInc.SumHIncCon += (float)((rfield.flux[i] + rfield.otscon[i] + rfield.otslin[i] + rfield.outcon[i] + rfield.outlin[i])* opac.opac[i]*rfield.anu[i]); } SumOutInc.SumOutCon += SumOutInc.SumHOutCon; SumOutInc.SumIncCon += SumOutInc.SumHIncCon; if( trace.lgTrace ) { fprintf( ioQQQ, " RTDiffuse returns.\n" ); } # ifdef DEBUG_FUN fputs( " <->RTDiffuse()\n", debug_fp ); # endif return; }