cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*RT_diffuse evaluate local diffuse emission for this zone, 00004 * fill in ConEmitLocal[depth][energy] with diffuse emission, 00005 * called by Cloudy, this routine adds energy to the outward beam 00006 * OTS rates for this zone were set in RT_OTS - not here */ 00007 #include "cddefines.h" 00008 #include "physconst.h" 00009 #include "taulines.h" 00010 #include "grains.h" 00011 #include "grainvar.h" 00012 #include "iso.h" 00013 #include "dense.h" 00014 #include "opacity.h" 00015 #include "trace.h" 00016 #include "coolheavy.h" 00017 #include "rfield.h" 00018 #include "phycon.h" 00019 #include "hmi.h" 00020 #include "radius.h" 00021 #include "atmdat.h" 00022 #include "heavy.h" 00023 #include "atomfeii.h" 00024 #include "lines_service.h" 00025 #include "h2.h" 00026 #include "ipoint.h" 00027 #include "rt.h" 00028 00029 #define TwoS (1+ipISO) 00030 00031 #if defined (__ICC) && defined(__ia64) && __INTEL_COMPILER < 910 00032 #pragma optimization_level 0 00033 #endif 00034 void RT_diffuse(void) 00035 { 00036 /* set this true to print a table of 2 photon emission coefficients 00037 * comparable to Brown & Mathews (1971) */ 00038 static long lgPrt2PhtEmisCoef = false; 00039 00040 /* arrays used in this routine 00041 * ARRAY rfield.ConEmitLocal[depth][energy] local emission per unit vol 00042 * ARRAY rfield.ConOTS_local_photons - local ots two-photon rates 00043 * ARRAY DiffuseEscape is the spectrum of diffuse emission that escapes this zone, 00044 * at end of this routine part is thrown into the outward beam 00045 * by adding to rfield.ConInterOut 00046 * units are photons s-1 cm-3 00047 * one-time init done on first call */ 00048 static realnum *DiffuseEscape; 00049 00050 /* DiffuseEscape and rfield.ConEmitLocal are same except that 00051 * rfield.ConEmitLocal is local emission, would be source function if div by opac 00052 * DiffuseEscape is part that escapes so has RT built into it 00053 * DiffuseEscape is used to define rfield.ConInterOut below as per this statement 00054 * rfield.ConInterOut[nu] += DiffuseEscape[nu]*(realnum)radius.dVolOutwrd; 00055 */ 00056 /* \todo 0 define only rfield.ConEmitLocal as it is now done, 00057 * do not define DiffuseEscape at all 00058 * at bottom of this routine use inward and outward optical depths to define 00059 * local and escaping parts 00060 * this routine only defines 00061 * rfield.ConInterOut - set to DiffuseEscape times vol element 00062 * rfield.ConOTS_local_photons (two photon only) 00063 * so this is only var that 00064 * needs to be set 00065 */ 00066 static bool lgMustInit=true; 00067 00068 long int i=-100000, 00069 ip=-100000, 00070 ipISO=-100000, 00071 ipHi=-100000, 00072 ipLo=-100000, 00073 ipla=-100000, 00074 nelem=-100000, 00075 ion=-100000, 00076 limit=-100000, 00077 n=-100000, 00078 nu=-10000; 00079 00080 double EnergyLimit; 00081 00082 double EdenAbund, 00083 difflya, 00084 xInWrd, 00085 arg, 00086 fac, 00087 factor, 00088 gamma, 00089 gion, 00090 gn, 00091 photon, 00092 sum, 00093 Sum1level, 00094 SumCaseB; 00095 00096 realnum Dilution , two_photon_emiss; 00097 00098 DEBUG_ENTRY( "RT_diffuse()" ); 00099 00100 /* one time creation of space for ThowOut array */ 00101 if( lgMustInit ) 00102 { 00103 lgMustInit = false; 00104 DiffuseEscape = (realnum*)MALLOC((unsigned)rfield.nupper*sizeof(realnum)); 00105 } 00106 00107 /* many arrays were malloced to nupper, and we will add unit flux to [nflux] - 00108 8 this must be true to work */ 00109 ASSERT(rfield.nflux < rfield.nupper ); 00110 00111 /* this routine evaluates the local diffuse fields 00112 * it fills in all of the following vectors */ 00113 memset(DiffuseEscape , 0 , (unsigned)rfield.nupper*sizeof(realnum) ); 00114 memset(rfield.ConEmitLocal[nzone] , 0 , (unsigned)rfield.nupper*sizeof(realnum) ); 00115 memset(rfield.ConOTS_local_photons , 0 , (unsigned)rfield.nupper*sizeof(realnum) ); 00116 memset(rfield.TotDiff2Pht , 0 , (unsigned)rfield.nupper*sizeof(realnum) ); 00117 00118 /* must abort after setting all of above to zero because some may be 00119 * used in various ways before abort is complete */ 00120 if( lgAbort ) 00121 { 00122 /* quit if we are aborting */ 00123 return; 00124 } 00125 00126 /* loop over iso-sequences of all elements 00127 * to add all recombination continua and lines*/ 00128 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00129 { 00130 /* >>chng 01 sep 23, rewrote for iso sequences */ 00131 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00132 { 00133 /* this will be the sum of recombinations to all excited levels */ 00134 SumCaseB = 0.; 00135 00136 /* the product of the densities of the parent ion and electrons */ 00137 EdenAbund = dense.eden*dense.xIonDense[nelem][nelem+1-ipISO]; 00138 00139 /* recombination continua for all iso seq - 00140 * if this stage of ionization exists */ 00141 if( dense.IonHigh[nelem] >= nelem+1-ipISO ) 00142 { 00143 /* loop over all levels to include recombination diffuse continua, 00144 * pick highest energy continuum point that opacities extend to 00145 * for ground continuum, go up to highest defined Boltzmann factor, 00146 * at bottom of loop will be reset to ground state photo edge */ 00147 ipHi = rfield.nflux; 00148 /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */ 00149 for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ ) 00150 { 00151 Sum1level = 0.; 00152 00153 /* >>chng 02 nov 20 - pull the plug on old he treatment */ 00154 /* the number is (2 pi me k/h^2) ^ -3/2 * 8 pi/c^2 / ge - it includes 00155 * the stat weight of the free electron in the demominator */ 00157 /*gamma = 2.0618684e11*StatesElem[ipISO][nelem][n].g/iso.stat_ion[ipISO]/phycon.te/phycon.sqrte;*/ 00158 gamma = 0.5*MILNE_CONST*StatesElem[ipISO][nelem][n].g/iso.stat_ion[ipISO]/phycon.te/phycon.sqrte; 00159 00160 /* loop over all recombination continua 00161 * escaping part of recombinations are added to rfield.ConEmitLocal 00162 * added to ConInterOut at end of routine */ 00163 for( nu=iso.ipIsoLevNIonCon[ipISO][nelem][n]-1; nu < ipHi; nu++ ) 00164 { 00165 /* dwid used to adjust where within WIDFLX exp is evaluated - 00166 * weighted to lower energy due to exp(-energy/T) */ 00167 double dwid = 0.2; 00168 00169 /* this is the term in the negative exponential Boltzmann factor 00170 * for continuum emission */ 00171 arg = (rfield.anu[nu]-iso.xIsoLevNIonRyd[ipISO][nelem][n]+ 00172 rfield.widflx[nu]*dwid)/phycon.te_ryd; 00173 arg = MAX2(0.,arg); 00174 /* don't bother evaluating for this or higher energies if 00175 * Boltzmann factor is tiny. */ 00176 if( arg > SEXP_LIMIT ) 00177 break; 00178 00179 /* photon is in photons cm^3 s^-1 per cell */ 00180 photon = gamma*exp(-arg)*rfield.widflx[nu]* 00181 opac.OpacStack[ nu-iso.ipIsoLevNIonCon[ipISO][nelem][n] + iso.ipOpac[ipISO][nelem][n] ] * 00182 rfield.anu2[nu]; 00183 00184 /*ASSERT( photon >= 0. );*/ 00185 Sum1level += photon; 00186 /* total local diffuse emission */ 00187 rfield.ConEmitLocal[nzone][nu] += (realnum)(photon*EdenAbund); 00188 /* iso.RadRecomb[ipISO][nelem][n][ipRecEsc] is escape probability 00189 * DiffuseEscape is local emission that escapes this zone */ 00190 DiffuseEscape[nu] += 00191 (realnum)(photon*EdenAbund*iso.RadRecomb[ipISO][nelem][n][ipRecEsc]); 00192 } 00193 /* this will be used below to confirm case B sum */ 00194 if( n > 0 ) 00195 { 00196 /* SumCaseB will be sum to all excited */ 00197 SumCaseB += Sum1level; 00198 } 00199 00200 /* on entry to this loop ipHi was set to the upper limit of the code, 00201 * and ground state recombination continua for all energies was added. For 00202 * excited states (which are the ones that will be done after 00203 * first pass through) only go to ground state threshold, 00204 * since that will be so much stronger than excited state recom*/ 00205 ipHi = iso.ipIsoLevNIonCon[ipISO][nelem][0]-1; 00206 } 00207 /* this is check on self-consistency */ 00208 iso.CaseBCheck[ipISO][nelem] = MAX2(iso.CaseBCheck[ipISO][nelem], 00209 (realnum)(SumCaseB/iso.RadRec_caseB[ipISO][nelem])); 00210 00211 /* this add line emission from the model atoms, 00212 * do not include very highest level since disturbed by topoff */ 00213 /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */ 00214 for( ipLo=0; ipLo < (iso.numLevels_local[ipISO][nelem] - 2); ipLo++ ) 00215 { 00216 /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */ 00217 for( ipHi=ipLo+1; ipHi < iso.numLevels_local[ipISO][nelem] - 1; ipHi++ ) 00218 { 00219 /* must not include fake transitions (2s-1s, two photon, or truly 00220 * forbidden transitions */ 00221 if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont < 1 ) 00222 continue; 00223 00224 /* pointer to line energy in continuum array */ 00225 ip = Transitions[ipISO][nelem][ipHi][ipLo].ipCont-1; 00226 00227 /* number of photons in the line has not been defined up until now, 00228 * do so now. this is redone in lines. */ 00229 Transitions[ipISO][nelem][ipHi][ipLo].Emis->phots = 00230 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul* 00231 StatesElem[ipISO][nelem][ipHi].Pop* 00232 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc* 00233 dense.xIonDense[nelem][nelem+1-ipISO]; 00234 00235 /* inward fraction of line */ 00236 xInWrd = Transitions[ipISO][nelem][ipHi][ipLo].Emis->phots* 00237 Transitions[ipISO][nelem][ipHi][ipLo].Emis->FracInwd; 00238 00239 /* reflected part of line */ 00240 rfield.reflin[ip] += (realnum)(xInWrd*radius.BeamInIn); 00241 00242 /* inward beam that goes out since sphere set */ 00243 /* in all this the ColOvTot term has been commented out, 00244 * since this is not meaningful for something like the hydrogen atom, 00245 * where most forms by recombination */ 00246 rfield.outlin[ip] += (realnum)(xInWrd*radius.BeamInOut*opac.tmn[ip]/* 00247 Transitions[ipISO][nelem][ipHi][ipLo].ColOvTot*/); 00248 00249 /* outward part of line */ 00250 rfield.outlin[ip] += 00251 (realnum)(Transitions[ipISO][nelem][ipHi][ipLo].Emis->phots* 00252 (1. - Transitions[ipISO][nelem][ipHi][ipLo].Emis->FracInwd)* 00253 radius.BeamOutOut*opac.tmn[ip]/* 00254 Transitions[ipISO][nelem][ipHi][ipLo].ColOvTot*/); 00255 } 00256 } 00257 00258 /*Iso treatment of two photon emission. */ 00259 /* NISO could in the future be increased, but we want this assert to blow 00260 * so that it is understood this may not be correct for other iso sequences, 00261 * probably should break since will not be present */ 00262 ASSERT( ipISO <= ipHE_LIKE ); 00263 00264 /* upper limit to 2-phot is energy of 2s to ground */ 00265 EnergyLimit = Transitions[ipISO][nelem][TwoS][0].EnergyWN/RYD_INF; 00266 00267 /* this could fail when pops very low and Pop2Ion is zero */ 00268 ASSERT( iso.ipTwoPhoE[ipISO][nelem]>0 && EnergyLimit>0. ); 00269 00270 sum = 0.; 00271 /* two photon emission, iso.ipTwoPhoE[ipISO][nelem] is 00272 * continuum array index for Lya energy */ 00273 for( nu=0; nu<iso.ipTwoPhoE[ipISO][nelem]; nu++ ) 00274 { 00275 /* We do not assert rfield.anu[nu]<=EnergyLimit because the maximum 00276 * index could be set to point to the next highest bin. */ 00277 ASSERT( rfield.anu[nu]/EnergyLimit < 1.01 || rfield.anu[nu-1]<EnergyLimit); 00278 00279 /* iso.As2nu[ipISO][nelem][nu] is transition probability A per bin 00280 * So sum is the total transition probability - this sum should 00281 * add up to the A two photon */ 00282 sum += iso.As2nu[ipISO][nelem][nu]; 00283 00284 /* flag saying whether to also do he triplets */ 00285 const bool lgDo2TripSToo = false; 00286 00287 if( ipISO == ipHE_LIKE && lgDo2TripSToo ) 00288 { 00289 two_photon_emiss = 2.f*dense.xIonDense[nelem][nelem+1-ipISO]*iso.As2nu[ipISO][nelem][nu]* 00290 ((realnum)StatesElem[ipISO][nelem][TwoS].Pop+0.0001f*(realnum)StatesElem[ipISO][nelem][ipHe2s3S].Pop); 00291 } 00292 else 00293 { 00294 /* StatesElem[ipISO][nelem][TwoS].Pop is dimensionless and represents the population 00295 * of the TwoS level relative to that of the ion. dense.xIonDense[nelem][nelem+1-ipISO] 00296 * is the density of the current ion (cm^-3). The factor of 2 is for two photons per 00297 * transition. Thus two_photon_emiss has dimension photons cm-3 s-1. */ 00298 two_photon_emiss = 2.f*(realnum)StatesElem[ipISO][nelem][TwoS].Pop*dense.xIonDense[nelem][nelem+1-ipISO] 00299 *iso.As2nu[ipISO][nelem][nu]; 00300 } 00301 00302 /* >>chng 03 mar 26, only do if induced processes turned on, 00303 * otherwise inconsistent with rate solver treatment. */ 00304 /* >>chng 02 aug 14, add induced two photon emission */ 00305 /* product of occupation numbers for induced emission */ 00306 if( iso.lgInd2nu_On) 00307 { 00308 /* this is the total rate */ 00309 two_photon_emiss *= (1.f + rfield.SummedOcc[nu]) * 00310 (1.f+rfield.SummedOcc[iso.ipSym2nu[ipISO][nelem][nu]-1]); 00311 } 00312 00313 /* information - only used in punch output */ 00314 rfield.TotDiff2Pht[nu] += two_photon_emiss; 00315 00316 /* total local diffuse emission */ 00317 rfield.ConEmitLocal[nzone][nu] += two_photon_emiss; 00318 00319 /* this is escaping part of two-photon emission, 00320 * as determined from optical depth to illuminated face */ 00321 DiffuseEscape[nu] += two_photon_emiss*opac.ExpmTau[nu]; 00322 00323 /* save locally destroyed OTS two-photon continuum - this is only 00324 * place this is set in this routine */ 00325 rfield.ConOTS_local_photons[nu] += two_photon_emiss*(1.f - opac.ExpmTau[nu]); 00326 } 00327 00328 /* a sanity check on the code, see Spitzer and Greenstein (1951), eqn 4. */ 00329 /* >>refer HI 2nu Spitzer, L., & Greenstein, J., 1951, ApJ, 114, 407 */ 00330 ASSERT( fabs( 1.f - (realnum)sum/Transitions[ipISO][nelem][TwoS][0].Emis->Aul ) < 0.01f ); 00331 } 00332 00333 /* option to print hydrogen and helium two-photon emission coefficients? */ 00334 if( lgPrt2PhtEmisCoef ) 00335 { 00336 long yTimes20; 00337 double y, E2nu; 00338 00339 lgPrt2PhtEmisCoef = false; 00340 00341 fprintf( ioQQQ, "\ny\tGammaNot(2q)\n"); 00342 00343 for( yTimes20=1; yTimes20<=10; yTimes20++ ) 00344 { 00345 y = yTimes20/20.; 00346 00347 fprintf( ioQQQ, "%.3e\t", (double)y ); 00348 00349 E2nu = Transitions[0][0][1][0].EnergyWN/RYD_INF; 00350 i = ipoint(y*E2nu); 00351 fprintf( ioQQQ, "%.3e\t", 00352 8./3.*HPLANCK*StatesElem[0][0][1].Pop/dense.eden 00353 *y*iso.As2nu[0][0][i]*E2nu/rfield.widflx[i] ); 00354 00355 E2nu = Transitions[1][1][2][0].EnergyWN/RYD_INF; 00356 i = ipoint(y*E2nu); 00357 fprintf( ioQQQ, "%.3e\n", 00358 8./3.*HPLANCK*StatesElem[1][1][2].Pop/dense.eden 00359 *y*iso.As2nu[1][1][i]*E2nu/rfield.widflx[i] ); 00360 00361 /* 00362 fprintf( ioQQQ, "%.3e\t%.3e\n", 00363 rfield.TotDiff2Pht[i]*rfield.anu[i]*EN1RYD 00364 /E2nu/rfield.widflx[i]/dense.eden/dense.xIonDense[nelem][nelem+1-ipISO]/FR1RYD, 00365 8./3.*HPLANCK*StatesElem[ipISO][nelem][TwoS].Pop/dense.eden 00366 *y*iso.As2nu[ipISO][nelem][i]*E2nu/rfield.widflx[i]); */ 00367 } 00368 fprintf( ioQQQ, "eden%.3e\n", dense.eden ); 00369 } 00370 } 00371 } 00372 00373 /* add recombination continua for elements heavier than those done with iso seq */ 00374 for( nelem=NISO; nelem < LIMELM; nelem++ ) 00375 { 00376 /* do not include species with iso-sequence in following */ 00377 /* >>chng 03 sep 09, upper bound was wrong, did not include NISO */ 00378 for( ion=dense.IonLow[nelem]; ion < nelem-NISO+1; ion++ ) 00379 { 00380 if( dense.xIonDense[nelem][ion+1] > 0. ) 00381 { 00382 long int ns, nshell,igRec , igIon, 00383 iplow , iphi , ipop; 00384 00385 ip = Heavy.ipHeavy[nelem][ion]-1; 00386 ASSERT( ip >= 0 ); 00387 00388 /* nflux was reset upward in ConvInitSolution to encompass all 00389 * possible line and continuum emission. this test should not 00390 * possibly fail. It could if the ionization were to increase with depth 00391 * although the continuum mesh is designed to deal with this. 00392 * This test is important because the nflux cell in ConInterOut 00393 * is used to carry out the unit integration, and if it gets 00394 * clobbered by diffuse emission the code will declare 00395 * insanity in PrtComment */ 00396 if( ip >= rfield.nflux ) 00397 continue; 00398 00399 /* get shell number, stat weights for this species */ 00400 atmdat_outer_shell( nelem+1 , nelem+1-ion , &nshell, &igRec , &igIon ); 00401 gn = (double)igRec; 00402 gion = (double)igIon; 00403 00404 /* shell number */ 00405 ns = Heavy.nsShells[nelem][ion]-1; 00406 ASSERT( ns == (nshell-1) ); 00407 00408 /* lower and upper energies, and offset for opacity stack */ 00409 iplow = opac.ipElement[nelem][ion][ns][0]-1; 00410 iphi = opac.ipElement[nelem][ion][ns][1]; 00411 iphi = MIN2( iphi , rfield.nflux ); 00412 ipop = opac.ipElement[nelem][ion][ns][2]; 00413 00414 /* now convert ipop to the offset in the opacity stack from threshold */ 00415 ipop = ipop - iplow; 00416 00417 gamma = 0.5*MILNE_CONST*gn/gion/phycon.te/phycon.sqrte*dense.eden*dense.xIonDense[nelem][ion+1]; 00418 00419 /* this is ground state continuum from stored opacities */ 00420 if( rfield.ContBoltz[iplow] > SMALLFLOAT ) 00421 { 00422 for( nu=iplow; nu < iphi; ++nu ) 00423 { 00424 photon = gamma*rfield.ContBoltz[nu]/rfield.ContBoltz[iplow]* 00425 rfield.widflx[nu]*opac.OpacStack[nu+ipop]*rfield.anu2[nu]; 00426 /* add heavy rec to ground in active beam,*/ 00431 rfield.ConEmitLocal[nzone][nu] += (realnum)photon; 00432 /*DiffuseEscape[i] += (realnum)photon;*/ 00433 /*>>chng 03 sep 08, index had been i, should have been nu */ 00434 DiffuseEscape[nu] += (realnum)photon*opac.ExpmTau[nu]; 00435 } 00436 } 00437 00438 /* now do the recombination Lya */ 00439 ipla = Heavy.ipLyHeavy[nelem][ion]-1; 00440 ASSERT( ipla >= 0 ); 00441 /* xLyaHeavy is set to a fraction of the total rad rec in ion_recomb, includes eden */ 00442 difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1]; 00443 /* >>chng 03 jul 10, here and below, use outlin_noplot */ 00444 rfield.outlin_noplot[ipla] += (realnum)(difflya*radius.dVolOutwrd*opac.tmn[ipla]*opac.ExpmTau[ipla]); 00445 00446 /* now do the recombination Balmer photons */ 00447 ipla = Heavy.ipBalHeavy[nelem][ion]-1; 00448 ASSERT( ipla >= 0 ); 00449 /* xLyaHeavy is set to fraction of total rad rec in ion_recomb, includes eden */ 00450 difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1]; 00451 rfield.outlin_noplot[ipla] += (realnum)(difflya*radius.dVolOutwrd*opac.tmn[ipla]*opac.ExpmTau[ipla]); 00452 } 00453 } 00454 } 00455 00456 /* free-free free free brems emission for all ions */ 00457 limit = MIN2( rfield.ipMaxBolt , rfield.nflux ); 00458 for( nu=0; nu < limit; nu++ ) 00459 { 00460 double TotBremsAllIons = 0., BremsThisIon; 00461 00462 /* do hydrogen first, before main loop since want to add on H- brems */ 00463 nelem = ipHYDROGEN; 00464 ion = 1; 00465 BremsThisIon = POW2( (realnum)ion )*dense.xIonDense[nelem][ion]*rfield.gff[ion][nu]; 00466 ASSERT( BremsThisIon >= 0. ); 00467 00468 /* for case of hydrogen, add H- brems - OpacStack contains the ratio 00469 * of the H- to H brems cross section - multiply by this and the fraction in ground */ 00470 BremsThisIon *= (1.+opac.OpacStack[nu-1+opac.iphmra]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop); 00471 TotBremsAllIons = BremsThisIon; 00472 00473 /* chng 02 may 16, by Ryan...do all brems for all ions in one fell swoop, 00474 * using gaunt factors from rfield.gff. */ 00475 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ ) 00476 { 00477 /* MAX2 occurs because we want to start at first ion (or above) 00478 * and do not want atom */ 00479 for( ion=MAX2(1,dense.IonLow[nelem]); ion<=dense.IonHigh[nelem]; ++ion ) 00480 { 00481 /* eff. charge is ion, so first rfield.gff argument must be "ion". */ 00482 BremsThisIon = POW2( (realnum)ion )*dense.xIonDense[nelem][ion]*rfield.gff[ion][nu]; 00483 /*ASSERT( BremsThisIon >= 0. );*/ 00484 00485 TotBremsAllIons += BremsThisIon; 00486 } 00487 } 00488 00490 /* >>chng 06 apr 05, no free free also turns off emission */ 00491 TotBremsAllIons *= dense.eden*1.032e-11*rfield.widflx[nu]*rfield.ContBoltz[nu]/rfield.anu[nu]/phycon.sqrte * 00492 CoolHeavy.lgFreeOn; 00493 ASSERT( TotBremsAllIons >= 0.); 00494 00495 /* >>chng 01 jul 01, move thick brems back to ConEmitLocal but do not add 00496 * to outward beam - ConLocNoInter array removed as result 00497 * if problems develop with very dense blr clouds, this may be reason */ 00498 /*rfield.ConLocNoInter[nu] += (realnum)fac;*/ 00499 /*rfield.ConEmitLocal[nzone][nu] += (realnum)TotBremsAllIons;*/ 00500 00501 if( nu >= rfield.ipEnergyBremsThin ) 00502 { 00503 /* >>chng 05 feb 20, move into this test on brems opacity - should not be 00504 * needed since would use expmtau to limit outward beam */ 00505 /* >>chng 01 jul 01, move thick brems back to ConEmitLocal but do not add 00506 * to outward beam - ConLocNoInter array removed as result 00507 * if problems develop with very dense BLR clouds, this may be reason */ 00508 /*rfield.ConLocNoInter[nu] += (realnum)fac;*/ 00509 rfield.ConEmitLocal[nzone][nu] += (realnum)TotBremsAllIons; 00510 00511 /* do not add optically thick part to outward beam since self absorbed 00512 * >>chng 96 feb 27, put back into outward beam since do not integrate 00513 * over it anyway. */ 00514 /* >>chng 99 may 28, take back out of beam since DO integrate over it 00515 * in very dense BLR clouds */ 00516 /* >>chng 01 jul 10, add here, in only one loop, where optically thin */ 00517 DiffuseEscape[nu] += (realnum)TotBremsAllIons; 00518 } 00519 } 00520 00521 /* grain dust emission */ 00522 /* >>chng 01 nov 22, moved calculation of grain flux to qheat.c, PvH */ 00523 if( gv.lgDustOn && gv.lgGrainPhysicsOn ) 00524 { 00525 /* this calculates diffuse emission from grains, 00526 * and stores the result in gv.GrainEmission */ 00527 GrainMakeDiffuse(); 00528 00529 for( nu=0; nu < rfield.nflux; nu++ ) 00530 { 00531 rfield.ConEmitLocal[nzone][nu] += gv.GrainEmission[nu]; 00532 DiffuseEscape[nu] += gv.GrainEmission[nu]; 00533 } 00534 } 00535 00536 /* hminus emission */ 00537 fac = dense.eden*(double)dense.xIonDense[ipHYDROGEN][0]; 00538 gn = 1.; 00539 gion = 2.; 00540 gamma = 0.5*MILNE_CONST*gn/gion/phycon.te/phycon.sqrte; 00541 /* >>chng 00 dec 15 change limit to -1 of H edge */ 00542 limit = MIN2(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1,rfield.nflux); 00543 00544 if( rfield.ContBoltz[hmi.iphmin-1] > 0. ) 00545 { 00546 for( nu=hmi.iphmin-1; nu < limit; nu++ ) 00547 { 00548 /* H- flux photons cm-3 s-1 00549 * ContBoltz is ratio of Boltzmann factor for each freq */ 00550 factor = gamma*rfield.ContBoltz[nu]/rfield.ContBoltz[hmi.iphmin-1]*rfield.widflx[nu]* 00551 opac.OpacStack[nu-hmi.iphmin+opac.iphmop]* 00552 rfield.anu2[nu]*fac; 00553 rfield.ConEmitLocal[nzone][nu] += (realnum)factor; 00554 DiffuseEscape[nu] += (realnum)factor; 00555 } 00556 } 00557 else 00558 { 00559 for( nu=hmi.iphmin-1; nu < limit; nu++ ) 00560 { 00561 arg = MAX2(0.,TE1RYD*(rfield.anu[nu]-0.05544)/phycon.te); 00562 /* this is the limit sexp normally uses */ 00563 if( arg > SEXP_LIMIT ) 00564 break; 00565 /* H- flux photons cm-3 s-1 00566 * flux is in photons per sec per ryd */ 00567 factor = gamma*exp(-arg)*rfield.widflx[nu]* 00568 opac.OpacStack[nu-hmi.iphmin+opac.iphmop]* 00569 rfield.anu2[nu]*fac; 00570 rfield.ConEmitLocal[nzone][nu] += (realnum)factor; 00571 DiffuseEscape[nu] += (realnum)factor; 00572 } 00573 } 00574 00575 /* this dilution is needed to conserve volume in spherical models. tests such 00576 * as parispn.in will fault if this is removed */ 00577 Dilution = (realnum)POW2( radius.rinner / (radius.Radius-radius.drad/2.) ); 00578 00579 /* this is a unit of energy that will be passed through the code as a test 00580 * that all integrations are carried out. A similar test is set in lineset1 00581 * and verified in PrtFinal. The opacity at this cell is zero so only 00582 * geometrical dilution will affect the integral 00583 * Radius is currently outer edge of zone, so radius-drad/2 is radius 00584 * of center of zone */ 00585 rfield.ConEmitLocal[nzone][rfield.nflux] = 1.e-10f * Dilution; 00586 DiffuseEscape[rfield.nflux] = 1.e-10f * Dilution; 00587 00588 /* opacity should be zero at this energy */ 00589 ASSERT( opac.opacity_abs[rfield.nflux] == 0. ); 00590 00591 /* many tmn added to conserve energy in very high Z models 00592 * rerun highZ qso model if any tmn ever deleted here 00593 * 00594 * tmn set in ZoneStart and includes both opacity and dilution 00595 * 00596 * use diffuse lines and continuum to add flux to outward beam 00597 * 00598 * NB!!! this routine adds flux to the outward beam 00599 * it can only be called once per zone 00600 */ 00601 00602 /* >>chng 96 nov 19, do not consider energies below plasma freq 00603 * ipPlasmaFreq points to lowest energy thin to brems and plasma freq */ 00604 00605 /* add DiffuseEscape continuum to ConInterOut, 00606 * lower limit of rfield.ipPlasma-1 since continuum is zero below rfield.ipPlasma-1 00607 * due to plasma frequency 00608 * note that upper range of sum is <= nflux, which contains the unit 00609 * verification token in the highest cell*/ 00610 for( nu=rfield.ipPlasma-1; nu <= rfield.nflux; nu++ ) 00611 { 00612 /* ConInterOut is the interactive continuum 00613 * DiffuseEscape contains all radiation thrown into outward beam */ 00614 /* radius.dVolOutwrd has factor or 1 or 0.5 for outward part */ 00615 /* >>chng 05 feb 23, activate this statement */ 00616 rfield.ConInterOut[nu] += DiffuseEscape[nu]*(realnum)radius.dVolOutwrd; 00617 /*rfield.ConInterOut[nu] += rfield.ConEmitLocal[nzone][nu]*opac.ExpmTau[nu]*(realnum)radius.dVolOutwrd;*/ 00618 ASSERT( rfield.ConInterOut[nu] >= 0.); 00619 } 00620 00621 { 00622 /*@-redef@*/ 00623 enum {DEBUG_LOC=false}; 00624 /*@+redef@*/ 00625 if( DEBUG_LOC ) 00626 { 00627 fprintf(ioQQQ,"rtdiffusedebugg %li\t%.2e\t%.2e\t%.2e\n", 00628 nzone, rfield.ConInterOut[1158] , DiffuseEscape[1158],radius.dVolOutwrd); 00629 } 00630 } 00631 00632 /* outward level 1 line photons, 0 is dummy line */ 00633 for( i=1; i <= nLevel1; i++ ) 00634 { 00635 outline( &TauLines[i] ); 00636 } 00637 00638 for( ipISO = ipHE_LIKE; ipISO < NISO; ipISO++ ) 00639 { 00640 for( nelem = ipISO; nelem < LIMELM; nelem++ ) 00641 { 00642 if( dense.lgElmtOn[nelem] ) 00643 { 00644 for( i=0; i<iso.numLevels_max[ipISO][nelem]; i++ ) 00645 { 00646 00647 SatelliteLines[ipISO][nelem][i].Emis->phots = 00648 SatelliteLines[ipISO][nelem][i].Emis->Aul* 00649 SatelliteLines[ipISO][nelem][i].Hi->Pop* 00650 SatelliteLines[ipISO][nelem][i].Emis->Pesc* 00651 dense.xIonDense[nelem][nelem+1-ipISO]; 00652 00653 outline( &SatelliteLines[ipISO][nelem][i] ); 00654 } 00655 } 00656 } 00657 } 00658 00659 /* outward level 2 line photons */ 00660 for( i=0; i < nWindLine; i++ ) 00661 { 00662 /* must not also do lines that were already done as part 00663 * of the isoelectronic sequences */ 00664 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 00665 { 00666 { 00667 /*@-redef@*/ 00668 enum {DEBUG_LOC=false}; 00669 /*@+redef@*/ 00670 if( DEBUG_LOC /*&& nzone > 10*/ && i==4821 ) 00671 { 00672 /* set up to dump the Fe 9 169A line */ 00673 fprintf(ioQQQ,"DEBUG dump lev2 line %li\n", i ); 00674 DumpLine( &TauLine2[i] ); 00675 fprintf(ioQQQ,"DEBUG dump %.3e %.3e %.3e\n", 00676 rfield.outlin[TauLine2[i].ipCont-1], 00677 TauLine2[i].Emis->phots*TauLine2[i].Emis->FracInwd*radius.BeamInOut*opac.tmn[i]*TauLine2[i].Emis->ColOvTot, 00678 TauLine2[i].Emis->phots*(1. - TauLine2[i].Emis->FracInwd)*radius.BeamOutOut* TauLine2[i].Emis->ColOvTot ); 00679 } 00680 } 00681 outline( &TauLine2[i] ); 00682 /*if( i==2576 ) fprintf(ioQQQ,"DEBUG dump %.3e %.3e \n", 00683 rfield.outlin[TauLine2[i].ipCont-1] , rfield.outlin_noplot[TauLine2[i].ipCont-1]);*/ 00684 } 00685 } 00686 00687 /* outward hyperfine structure line photons */ 00688 for( i=0; i < nHFLines; i++ ) 00689 { 00690 outline( &HFLines[i] ); 00691 } 00692 00693 /*Atoms & Molecules*/ 00694 for( i=0; i < linesAdded2; i++ ) 00695 { 00696 outline( atmolEmis[i].tran ); 00697 } 00698 00699 /* carbon monoxide */ 00700 for( i=0; i < nCORotate; i++ ) 00701 { 00702 outline( &C12O16Rotate[i] ); 00703 outline( &C13O16Rotate[i] ); 00704 } 00705 00706 /* H2 emission */ 00707 H2_RT_diffuse(); 00708 00709 /* do outward parts of FeII lines, if large atom is turned on */ 00710 FeII_RT_Out(); 00714 if( trace.lgTrace ) 00715 fprintf( ioQQQ, " RT_diffuse returns.\n" ); 00716 00717 /* >>chng 02 jul 25, zero out all light below plasma freq */ 00718 for( nu=0; nu<rfield.ipPlasma-1; nu++ ) 00719 { 00720 rfield.flux_beam_const[nu] = 0.; 00721 rfield.flux_beam_time[nu] = 0.; 00722 rfield.flux_isotropic[nu] = 0.; 00723 rfield.flux[nu] = 0.; 00724 rfield.ConEmitLocal[nzone][nu] = 0.; 00725 rfield.otscon[nu] =0.; 00726 rfield.otslin[nu] =0.; 00727 rfield.outlin[nu] =0.; 00728 rfield.outlin_noplot[nu] =0.; 00729 rfield.reflin[nu] = 0.; 00730 rfield.TotDiff2Pht[nu] = 0.; 00731 rfield.ConOTS_local_photons[nu] = 0.; 00732 rfield.ConInterOut[nu] = 0.; 00733 } 00734 00735 /* find occupation number, also assert that no continua are negative */ 00736 for( nu=0; nu < rfield.nflux; nu++ ) 00737 { 00738 /* >>chng 00 oct 03, add diffuse continua */ 00739 /* local diffuse continua */ 00740 rfield.OccNumbDiffCont[nu] = rfield.ConEmitLocal[nzone][nu]*rfield.convoc[nu]; 00741 00742 /* confirm that all are non-negative */ 00743 ASSERT( rfield.flux_beam_const[nu] >= 0.); 00744 ASSERT( rfield.flux_beam_time[nu] >= 0.); 00745 ASSERT( rfield.flux_isotropic[nu] >= 0.); 00746 ASSERT( rfield.flux[nu] >=0.); 00747 ASSERT( rfield.ConEmitLocal[nzone][nu] >= 0.); 00748 ASSERT( rfield.otscon[nu] >=0.); 00749 ASSERT( rfield.otslin[nu] >=0.); 00750 ASSERT( rfield.outlin[nu] >=0.); 00751 ASSERT( rfield.outlin_noplot[nu] >=0.); 00752 ASSERT( rfield.reflin[nu] >= 0.); 00753 ASSERT( rfield.TotDiff2Pht[nu] >= 0.); 00754 ASSERT( rfield.ConOTS_local_photons[nu] >= 0.); 00755 ASSERT( rfield.ConInterOut[nu] >=0.); 00756 } 00757 00758 /* option to kill outward lines with no outward lines command*/ 00759 if( rfield.lgKillOutLine ) 00760 { 00761 for( nu=0; nu < rfield.nflux; nu++ ) 00762 { 00763 rfield.outlin[nu] = 0.; 00764 rfield.outlin_noplot[nu] = 0.; 00765 } 00766 } 00767 00768 /* option to kill outward continua with no outward continua command*/ 00769 if( rfield.lgKillOutCont ) 00770 { 00771 for( nu=0; nu < rfield.nflux; nu++ ) 00772 { 00773 rfield.ConInterOut[nu] = 0.; 00774 } 00775 } 00776 return; 00777 }