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_tau_init set initial outward optical depths at start of first iteration, 00004 * it is only called by cloudy one time per complete calculation, just after 00005 * continuum set up and before start of convergence attempts. 00006 * 00007 * RT_tau_reset after first iteration, updates the optical depths, mirroring this 00008 * routine but with the previous iteration's variables */ 00009 #include "cddefines.h" 00010 #define TAULIM 1e8 00011 #include "taulines.h" 00012 #include "doppvel.h" 00013 #include "iso.h" 00014 #include "h2.h" 00015 #include "lines_service.h" 00016 #include "rfield.h" 00017 #include "dense.h" 00018 #include "opacity.h" 00019 #include "thermal.h" 00020 #include "geometry.h" 00021 #include "stopcalc.h" 00022 #include "ipoint.h" 00023 #include "abund.h" 00024 #include "conv.h" 00025 #include "atomfeii.h" 00026 #include "rt.h" 00027 #include "trace.h" 00028 00029 void RT_tau_init(void) 00030 { 00031 long int i, 00032 nelem, 00033 ipISO, 00034 ipHi, 00035 ipLo, 00036 nHi; 00037 long lgHit; /* used for logic check */ 00038 00039 double AbunRatio, 00040 balc, 00041 coleff; 00042 00043 realnum f, BalmerTauOn; 00044 00045 DEBUG_ENTRY( "RT_tau_init()" ); 00046 00047 ASSERT( dense.eden > 0. ); 00048 00049 /* Zero lines first */ 00050 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00051 { 00052 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00053 { 00054 if( dense.lgElmtOn[nelem] ) 00055 { 00056 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00057 { 00058 if( iso.lgDielRecom[ipISO] ) 00059 TransitionZero( &SatelliteLines[ipISO][nelem][ipHi] ); 00060 00061 for( ipLo=0; ipLo < ipHi; ipLo++ ) 00062 { 00063 TransitionZero( &Transitions[ipISO][nelem][ipHi][ipLo] ); 00064 } 00065 } 00066 for( ipHi=2; ipHi <iso.nLyman[ipISO]; ipHi++ ) 00067 { 00068 TransitionZero( &ExtraLymanLines[ipISO][nelem][ipHi] ); 00069 } 00070 } 00071 } 00072 } 00073 00074 /* initialize heavy element line array */ 00075 for( i=1; i <= nLevel1; i++ ) 00076 { 00077 TransitionZero( &TauLines[i] ); 00078 } 00079 /* this is a dummy optical depth array for non-existant lines 00080 * when this goes over to struc, make sure all are set to zero here since 00081 * init in grid may depend on it */ 00082 TransitionZero( &TauDummy ); 00083 00084 /* lines in cooling function with Mewe approximate collision strengths */ 00085 for( i=0; i < nWindLine; i++ ) 00086 { 00087 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 00088 { 00089 /* inward optical depth */ 00090 TransitionZero( &TauLine2[i] ); 00091 } 00092 } 00093 00094 /* inner shell lines */ 00095 for( i=0; i < nUTA; i++ ) 00096 { 00097 /* these are line optical depth arrays 00098 * inward optical depth */ 00099 /* heat is special for this array - it is heat per pump */ 00100 double hsave = UTALines[i].Coll.heat; 00101 TransitionZero( &UTALines[i] ); 00102 UTALines[i].Coll.heat = hsave; 00103 } 00104 00105 /* hyperfine structure lines */ 00106 for( i=0; i < nHFLines; i++ ) 00107 { 00108 TransitionZero( &HFLines[i] ); 00109 } 00110 00111 /*Atoms & Molecules*/ 00112 for( i=0; i < linesAdded2; i++) 00113 { 00114 TransitionZero( atmolEmis[i].tran ); 00115 } 00116 00117 /* CO carbon monoxide lines */ 00118 for( i=0; i < nCORotate; i++ ) 00119 { 00120 /* >>chng 03 feb 14, use main zero function */ 00121 TransitionZero( &C12O16Rotate[i] ); 00122 TransitionZero( &C13O16Rotate[i] ); 00123 } 00124 00125 /* initialize optical depths in H2 */ 00126 H2_LineZero(); 00127 00128 /* initialize large atom FeII arrays */ 00129 FeII_LineZero(); 00130 00131 /* ==================================================================*/ 00132 /* end setting lines to zero */ 00133 00134 /* >>chng 02 feb 08, moved to here from opacitycreateall, which was called later. 00135 * bug inhibited convergence in some models. Caught by PvH */ 00136 /* this is option to stop at certain optical depth */ 00137 if( StopCalc.taunu > 0. ) 00138 { 00139 StopCalc.iptnu = ipoint(StopCalc.taunu); 00140 StopCalc.iptnu = MIN2(StopCalc.iptnu,rfield.nupper-1); 00141 } 00142 else 00143 { 00144 StopCalc.iptnu = rfield.nupper; 00145 /* >>chng 04 dec 21, remove from here and init to 1e30 in zero */ 00146 /*StopCalc.tauend = 1e30f;*/ 00147 } 00148 00149 /* if taunu was set with a stop optical depth command, then iptnu must 00150 * have also been set to a continuum index before this code is reaced */ 00151 ASSERT( StopCalc.taunu == 0. || StopCalc.iptnu >= 0 ); 00152 00153 /* set initial and total optical depths in arrays 00154 * TAUNU is set when lines read in, sets stopping radius */ 00155 if( StopCalc.taunu > 0. ) 00156 { 00157 /* an optical depth has been specified */ 00158 if( StopCalc.iptnu >= iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] ) 00159 { 00160 /* at ionizing energies */ 00161 for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ ) 00162 { 00163 /* taumin set to 1e-20, can be reset with taumin command */ 00164 opac.TauAbsGeo[1][i] = opac.taumin; 00165 opac.TauScatGeo[1][i] = opac.taumin; 00166 opac.TauTotalGeo[1][i] = opac.taumin; 00167 } 00168 00169 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ ) 00170 { 00171 /* TauAbsGeo(i,2) = tauend * (anu(i)/anu(iptnu))**(-2.43) */ 00172 opac.TauAbsGeo[1][i] = StopCalc.tauend; 00173 opac.TauScatGeo[1][i] = opac.taumin; 00174 opac.TauTotalGeo[1][i] = opac.TauAbsGeo[1][i] + opac.TauScatGeo[1][i]; 00175 } 00176 } 00177 00178 else 00179 { 00180 /* not specified at ionizing energies */ 00181 for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ ) 00182 { 00183 opac.TauAbsGeo[1][i] = StopCalc.tauend; 00184 opac.TauScatGeo[1][i] = StopCalc.tauend; 00185 opac.TauTotalGeo[1][i] = StopCalc.tauend; 00186 } 00187 00188 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ ) 00189 { 00190 opac.TauAbsGeo[1][i] = (realnum)(TAULIM*pow(rfield.anu[i],(realnum)-2.43f)); 00191 opac.TauScatGeo[1][i] = opac.taumin; 00192 opac.TauTotalGeo[1][i] = opac.TauAbsGeo[1][i] + opac.TauScatGeo[1][i]; 00193 } 00194 00195 } 00196 } 00197 00198 else 00199 { 00200 /* ending optical depth not specified, assume 1E8 at 1 Ryd */ 00201 for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ ) 00202 { 00203 opac.TauAbsGeo[1][i] = opac.taumin; 00204 opac.TauScatGeo[1][i] = opac.taumin; 00205 opac.TauTotalGeo[1][i] = opac.taumin; 00206 } 00207 00208 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ ) 00209 { 00210 opac.TauAbsGeo[1][i] = (realnum)(TAULIM*pow(rfield.anu[i],(realnum)-2.43f)); 00211 opac.TauScatGeo[1][i] = opac.taumin; 00212 opac.TauTotalGeo[1][i] = opac.TauAbsGeo[1][i] + opac.TauScatGeo[1][i]; 00213 } 00214 } 00215 00216 /* if lgSphere then double outer, set inner to half 00217 * assume will be optically thin at sub-ionizing energies */ 00218 if( geometry.lgSphere || opac.lgCaseB ) 00219 { 00220 for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ ) 00221 { 00222 opac.TauAbsGeo[0][i] = opac.taumin; 00223 opac.TauAbsGeo[1][i] = opac.taumin*2.f; 00224 opac.TauScatGeo[0][i] = opac.taumin; 00225 opac.TauScatGeo[1][i] = opac.taumin*2.f; 00226 opac.TauTotalGeo[0][i] = 2.f*opac.taumin; 00227 opac.TauTotalGeo[1][i] = 4.f*opac.taumin; 00228 } 00229 00230 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ ) 00231 { 00232 opac.TauAbsGeo[0][i] = opac.TauAbsGeo[1][i]; 00233 opac.TauAbsGeo[1][i] *= 2.; 00234 opac.TauScatGeo[0][i] = opac.TauScatGeo[1][i]; 00235 opac.TauScatGeo[1][i] *= 2.; 00236 opac.TauTotalGeo[0][i] = opac.TauTotalGeo[1][i]; 00237 opac.TauTotalGeo[1][i] *= 2.; 00238 } 00239 00240 if( StopCalc.taunu > 0. ) 00241 { 00242 /* ending optical depth specified, and lgSphere also */ 00243 StopCalc.tauend *= 2.; 00244 } 00245 } 00246 00247 /* fix escape prob for He, metals, first set log of Te, needed by RECEFF 00248 * do not do this if temperature has been set by command */ 00249 if( !thermal.lgTemperatureConstant ) 00250 { 00251 double TeNew; 00252 /* this is a typical temperature for the H+ zone, and will use it is 00253 * the line widths & opt depth in the following. this is not meant to be the first 00254 * temperature during the search phase. */ 00255 TeNew = 2e4; 00256 /* >>chng 05 jul 19, in PDR models the guess of the optical depth in Lya could 00257 * be very good since often limiting column density is set, but must use 00258 * the temperature of H0 gas. in a dense BLR this is roughly 7000K and 00259 * closer to 100K for a low-density PDR - use these guesses */ 00260 if( dense.gas_phase[ipHYDROGEN] >= 1e9 ) 00261 { 00262 /* this is a typical BLR H0 temp */ 00263 TeNew = 7000.; 00264 } 00265 else if( dense.gas_phase[ipHYDROGEN] <= 1e5 ) 00266 { 00267 /* this is a typical PDR H0 temp */ 00268 TeNew = 100.; 00269 } 00270 else 00271 { 00272 /* power law interpolation between them */ 00273 TeNew = 0.5012 * pow( (double)dense.gas_phase[ipHYDROGEN], 0.46 ); 00274 } 00275 00276 /* propagate this temperature through the code */ 00277 /* must not call tfidle at this stage since not all vectors have been allocated */ 00278 TempChange(TeNew ); 00279 } 00280 00281 /* set inward optical depths for hydrogenic ions to small number proportional to abundance */ 00282 for( nelem=0; nelem < LIMELM; nelem++ ) 00283 { 00284 if( dense.lgElmtOn[nelem] ) 00285 { 00286 /* now get actual optical depths */ 00287 AbunRatio = abund.solar[nelem]/abund.solar[0]; 00288 for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ ) 00289 { 00290 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ ) 00291 { 00292 /* set all inward optical depths to taumin, regardless of abundance 00293 * this is a very small number, 1e-20 */ 00294 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn = opac.taumin; 00295 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot = 2.f*opac.taumin; 00296 } 00297 } 00298 00299 /* La may be case B, tlamin set to 1e9 by default with case b command */ 00300 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn = opac.tlamin; 00301 00302 /* scale factor so that all other lyman lines are appropriate for this Lya optical depth*/ 00303 f = opac.tlamin/Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity; 00304 fixit(); /* this appears to be redundant to code below. */ 00305 00306 for( nHi=3; nHi<=iso.n_HighestResolved_max[ipH_LIKE][nelem]; nHi++ ) 00307 { 00308 ipHi = iso.QuantumNumbers2Index[ipH_LIKE][nelem][nHi][1][2]; 00309 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn = MAX2( opac.taumin, 00310 f*Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity ); 00311 } 00312 for( ipHi=iso.numLevels_max[ipH_LIKE][nelem] - iso.nCollapsed_max[ipH_LIKE][nelem]; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ ) 00313 { 00314 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn = MAX2( opac.taumin, 00315 f*Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity ); 00316 } 00317 00318 /* after this set of if's the total lya optical depth will be known, 00319 * common code later on will set rest of lyman lines 00320 * if case b then set total lyman to twice inner */ 00321 if( opac.lgCaseB ) 00322 { 00323 /* force outer optical depth to twice inner if case B */ 00324 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot = 00325 (realnum)(2.*Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn); 00326 /* force off Balmer et al optical depths */ 00327 BalmerTauOn = 0.f; 00328 } 00329 00330 else 00331 { 00332 /* usual case for H LyA; try to guess outer optical depth */ 00333 if( StopCalc.colnut < 6e29 ) 00334 { 00335 /* neutral column is defined */ 00336 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot = (realnum)(StopCalc.colnut* 00337 rt.DoubleTau*Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity/ 00338 DoppVel.doppler[nelem]*AbunRatio); 00339 ASSERT( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot > 0. ); 00340 } 00341 /* has optical depth at some energy where we want to stop been specified? 00342 * taunu is energy where 00343 * this has been specified - this checks on Lyman continuum, taunu = 1 */ 00344 else if( StopCalc.taunu < 3. && StopCalc.taunu >= 0.99 ) 00345 { 00346 /* Lyman continuum optical depth defined */ 00347 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot = (realnum)(StopCalc.tauend* 00348 1.2e4*1.28e6/DoppVel.doppler[nelem]*rt.DoubleTau*AbunRatio); 00349 ASSERT( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot > 0. ); 00350 } 00351 else if( StopCalc.HColStop < 6e29 ) 00352 { 00353 /* neutral column not defined, guess from total col and U */ 00354 coleff = StopCalc.HColStop - 00355 MIN2(MIN2(rfield.qhtot/dense.eden,1e24)/2.6e-13,0.6*StopCalc.HColStop); 00356 00357 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot = (realnum)(coleff* 00358 7.6e-14*AbunRatio); 00359 ASSERT( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot > 0. ); 00360 } 00361 else 00362 { 00363 /* no way to estimate 912 optical depth, set to large number */ 00364 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot = (realnum)(1e20* 00365 AbunRatio); 00366 ASSERT( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot > 0. ); 00367 } 00368 /* allow Balmer et al. optical depths */ 00369 BalmerTauOn = 1.f; 00370 } 00371 00372 /* Lya total optical depth now known, is it optically thick?*/ 00373 if( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot > 1. ) 00374 { 00375 rt.TAddHLya = (realnum)MIN2(1.,Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot/ 00376 1e4); 00377 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn += rt.TAddHLya; 00378 } 00379 00380 else 00381 { 00382 rt.TAddHLya = opac.tlamin; 00383 } 00384 00385 /* this scale factor is to set other lyman lines, given the Lya optical depth */ 00386 f = Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot/ 00387 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity; 00388 00389 ipISO = ipH_LIKE; 00390 ASSERT( ipISO<NISO && nelem < LIMELM ); 00391 for( nHi=3; nHi<=iso.n_HighestResolved_max[ipISO][nelem]; nHi++ ) 00392 { 00393 ipHi = iso.QuantumNumbers2Index[ipISO][nelem][nHi][1][2]; 00394 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauTot = MAX2( opac.taumin, 00395 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity * f ); 00396 00397 /* increment inward optical depths by rt all lyman lines, inward 00398 * optical depth was set above, this adds to it. this was originally 00399 * set some place above */ 00400 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn += rt.TAddHLya* 00401 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity/ 00402 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity; 00403 00404 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn = MAX2( 00405 opac.taumin, Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn ); 00406 } 00407 for( ipHi=iso.numLevels_max[ipH_LIKE][nelem] - iso.nCollapsed_max[ipH_LIKE][nelem]; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ ) 00408 { 00409 /* set total optical depth for higher lyman lines */ 00410 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauTot = MAX2( opac.taumin, 00411 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity * f ); 00412 00413 /* increment inward optical depths by rt all lyman lines, inward 00414 * optical depth was set above, this adds to it. this was originally 00415 * set some place above */ 00416 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn += rt.TAddHLya* 00417 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity/ 00418 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity; 00419 00420 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn = MAX2( 00421 opac.taumin, Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn ); 00422 } 00423 00424 /* try to guess what Balmer cont optical guess, 00425 * first branch is case where we will stop at a balmer continuum optical 00426 * depth - taunu is energy where tauend was set */ 00427 if( StopCalc.taunu > 0.24 && StopCalc.taunu < 0.7 ) 00428 { 00429 Transitions[ipH_LIKE][nelem][ipH3p][ipH2s].Emis->TauTot = (realnum)(StopCalc.tauend* 00430 3.7e4*BalmerTauOn*AbunRatio + 1e-20); 00431 00432 Transitions[ipH_LIKE][nelem][ipH3s][ipH2p].Emis->TauTot = (realnum)(StopCalc.tauend* 00433 3.7e4*BalmerTauOn*AbunRatio + 1e-20); 00434 00435 Transitions[ipH_LIKE][nelem][ipH3d][ipH2p].Emis->TauTot = (realnum)(StopCalc.tauend* 00436 3.7e4*BalmerTauOn*AbunRatio + 1e-20); 00437 } 00438 00439 else 00440 { 00441 /* this is a guess based on Ferland&Netzer 1979, but it gets very large */ 00442 balc = rfield.qhtot*2.1e-19*BalmerTauOn*AbunRatio + 1e-20; 00443 00444 Transitions[ipH_LIKE][nelem][ipH3p][ipH2s].Emis->TauTot = 00445 (realnum)(MIN2(2e5, balc*3.7e4*BalmerTauOn+1e-20)); 00446 00447 Transitions[ipH_LIKE][nelem][ipH3s][ipH2p].Emis->TauTot = 00448 (realnum)(MIN2(2e5, balc*3.7e4*BalmerTauOn+1e-20)); 00449 00450 Transitions[ipH_LIKE][nelem][ipH3d][ipH2p].Emis->TauTot = 00451 (realnum)(MIN2(2e5, balc*3.7e4*BalmerTauOn+1e-20)); 00452 00453 ASSERT( Transitions[ipH_LIKE][nelem][ipH3p][ipH2s].Emis->TauTot > 0.); 00454 ASSERT( Transitions[ipH_LIKE][nelem][ipH3s][ipH2p].Emis->TauTot > 0.); 00455 ASSERT( Transitions[ipH_LIKE][nelem][ipH3d][ipH2p].Emis->TauTot > 0.); 00456 } 00457 00458 /* 2s - 1s strongly forbidden */ 00459 Transitions[ipH_LIKE][nelem][ipH2s][ipH1s].Emis->TauTot = 2.f*opac.taumin; 00460 Transitions[ipH_LIKE][nelem][ipH2s][ipH1s].Emis->TauIn = opac.taumin; 00461 00462 /* transitions down to 2s are special since 'alpha' (2s-2p) has 00463 * small optical depth, so use 3 to 2p instead */ 00464 f = Transitions[ipH_LIKE][nelem][ipH3s][ipH2p].Emis->TauTot/ 00465 Transitions[ipH_LIKE][nelem][ipH3s][ipH2p].Emis->opacity* BalmerTauOn; 00466 00467 ipLo = ipH2s; 00468 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ ) 00469 { 00470 if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].ipCont <= 0 ) 00471 continue; 00472 00473 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot = opac.taumin + 00474 f* Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->opacity; 00475 ASSERT(Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot > 0.); 00476 } 00477 00478 /* this is for rest of lines, scale from opacity */ 00479 for( ipLo=ipH2p; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ ) 00480 { 00481 long ipISO = ipH_LIKE, ipNS, ipNPlusOneP; 00482 00483 /* scale everything with same factor we use for n+1 P -> nS */ 00484 ipNS = iso.QuantumNumbers2Index[ipH_LIKE][nelem][ N_(ipLo) ][0][2]; 00485 if( ( N_(ipLo) + 1 ) <= 00486 (iso.n_HighestResolved_max[ipH_LIKE][nelem] + iso.nCollapsed_max[ipH_LIKE][nelem] ) ) 00487 { 00488 ipNPlusOneP = iso.QuantumNumbers2Index[ipH_LIKE][nelem][ N_(ipLo)+1 ][1][2]; 00489 } 00490 else 00491 { 00492 ipNPlusOneP = ipNS + 1; 00493 } 00494 00495 #if 1 /* old way */ 00496 ipNS = ipLo; 00497 ipNPlusOneP = ipNS + 1; 00498 #endif 00499 00500 if( Transitions[ipH_LIKE][nelem][ipNPlusOneP][ipNS].ipCont <= 0 ) 00501 { 00502 f = SMALLFLOAT; 00503 } 00504 else 00505 { 00506 f = Transitions[ipH_LIKE][nelem][ipNPlusOneP][ipNS].Emis->TauTot/ 00507 Transitions[ipH_LIKE][nelem][ipNPlusOneP][ipNS].Emis->opacity* 00508 BalmerTauOn; 00509 } 00510 00511 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ ) 00512 { 00513 if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].ipCont <= 0 ) 00514 continue; 00515 00516 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot = opac.taumin + 00517 f* Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->opacity; 00518 ASSERT(Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot > 0.); 00519 } 00520 } 00521 00522 /* this loop is over all possible H lines, do some final cleanup */ 00523 for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ ) 00524 { 00525 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ ) 00526 { 00527 if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].ipCont <= 0 ) 00528 continue; 00529 00530 /* TauCon is line optical depth to inner face used for continuum pumping rate, 00531 * not equal to TauIn in case of static sphere since TauCon will never 00532 * count the far side line optical depth */ 00533 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauCon = 00534 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn; 00535 00536 /* make sure inward optical depth is not larger than half total */ 00537 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn = 00538 MIN2( Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn , 00539 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot/2.f ); 00540 ASSERT(Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn > 0.f); 00541 00542 /* this is fraction of line that goes inward, not known until second iteration*/ 00543 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->FracInwd = 0.5; 00544 } 00545 } 00546 } 00547 } 00548 00549 /* initialize all he-like optical depths */ 00550 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ ) 00551 { 00552 if( dense.lgElmtOn[nelem] ) 00553 { 00554 for( ipLo=0; ipLo < (iso.numLevels_max[ipHE_LIKE][nelem] - 1); ipLo++ ) 00555 { 00556 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ ) 00557 { 00558 /* set all inward optical depths to taumin, regardless of abundance 00559 * this is a very small number, 1e-20 */ 00560 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->TauIn = opac.taumin; 00561 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->TauTot = 2.f*opac.taumin; 00562 } 00563 } 00564 } 00565 } 00566 00567 /* now do helium like sequence if case b */ 00568 if( opac.lgCaseB ) 00569 { 00570 for( nelem=1; nelem < LIMELM; nelem++ ) 00571 { 00572 if( dense.lgElmtOn[nelem] ) 00573 { 00574 double Aprev; 00575 realnum ratio; 00576 /* La may be case B, tlamin set to 1e9 by default with case b command */ 00577 Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn = opac.tlamin; 00578 00579 Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauCon = 00580 Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn; 00581 00582 Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauTot = 00583 2.f*Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn; 00584 00585 ratio = opac.tlamin/Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->opacity; 00586 00587 /* this will be the trans prob of the previous lyman line, will use this to 00588 * find the next one up in the series */ 00589 Aprev = Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->Aul; 00590 00591 i = ipHe2p1P+1; 00592 /* >>chng 02 jan 05, remove explicit list of lyman lines, use As to guess 00593 * which are which - this will work for any number of levels */ 00594 for( i=ipHe2p1P+1; i < iso.numLevels_max[ipHE_LIKE][nelem]; i++ ) 00595 { 00596 if( Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].ipCont <= 0 ) 00597 continue; 00598 00599 /* >>chng 02 mar 15, add test for collapsed levels - As are very much 00600 * smaller at boundary between collapsed/resolved so this test is needed*/ 00601 if( Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->Aul> Aprev/10. || 00602 StatesElem[ipHE_LIKE][nelem][i].S < 0 ) 00603 { 00604 Aprev = Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->Aul; 00605 00606 Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn = 00607 ratio*Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->opacity; 00608 /* reset line optical depth to continuum source */ 00609 Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauCon = 00610 Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn; 00611 Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauTot = 00612 2.f*Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn; 00613 } 00614 } 00615 } 00616 } 00617 } 00618 00619 /* begin sanity check, total greater than 1/0.9 time inward */ 00620 lgHit = false; 00621 for( nelem=0; nelem < LIMELM; nelem++ ) 00622 { 00623 if( dense.lgElmtOn[nelem] ) 00624 { 00625 for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ ) 00626 { 00627 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ ) 00628 { 00629 if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].ipCont <= 0 ) 00630 continue; 00631 00632 if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot < 00633 0.9f*Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn ) 00634 { 00635 fprintf(ioQQQ, 00636 "RT_tau_init insanity for h line, Z=%li lo=%li hi=%li tot=%g in=%g \n", 00637 nelem , ipLo, ipHi , Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot , 00638 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn ); 00639 lgHit = true; 00640 } 00641 } 00642 } 00643 } 00644 } 00645 if( lgHit ) 00646 { 00647 fprintf( ioQQQ," stopping due to insanity in RT_tau_init\n"); 00648 ShowMe(); 00649 cdEXIT(EXIT_FAILURE); 00650 } 00651 /*end sanity check */ 00652 00653 /* fix offset for effective column density optical depth */ 00654 rt.tauxry = opac.TauAbsGeo[0][rt.ipxry-1]; 00655 00656 /* this is flag detected by dest prob routines to see whether ots rates are 00657 * oscaillating - damp them out if so */ 00658 conv.lgOscilOTS = false; 00659 00660 /* now that optical depths have been incremented, do escape prob, this 00661 * is located here instead on in cloudy.c (where it belongs) because 00662 * information generated by RT_line_all is needed for the following printout */ 00663 00664 RT_line_all(true , false ); 00665 00666 /* rest of routine is printout in case of trace */ 00667 if( trace.lgTrace ) 00668 { 00669 /* default is h-like */ 00670 ipISO = ipH_LIKE; 00671 if( trace.lgIsoTraceFull[ipHE_LIKE] ) 00672 ipISO = ipHE_LIKE; 00673 00674 if( trace.lgIsoTraceFull[ipISO] ) 00675 { 00676 fprintf( ioQQQ, "\n\n up Transitions[ipISO][hi][lo].TauTot array as set in RT_tau_init ipZTrace (nelem)= %ld\n", 00677 trace.ipIsoTrace[ipH_LIKE] ); 00678 for( ipHi=2; ipHi < iso.numLevels_max[ipISO][trace.ipIsoTrace[ipISO]]; ipHi++ ) 00679 { 00680 fprintf( ioQQQ, " %3ld", ipHi ); 00681 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00682 { 00683 if( Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul <= iso.SmallA ) 00684 continue; 00685 00686 fprintf( ioQQQ,PrintEfmt("%9.2e", 00687 Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->TauTot )); 00688 } 00689 fprintf( ioQQQ, "\n" ); 00690 } 00691 00692 fprintf( ioQQQ, "\n\n TauIn array\n" ); 00693 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][trace.ipIsoTrace[ipISO]]; ipHi++ ) 00694 { 00695 fprintf( ioQQQ, " %3ld", ipHi ); 00696 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00697 { 00698 if( Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul <= iso.SmallA ) 00699 continue; 00700 00701 fprintf( ioQQQ,PrintEfmt("%9.2e", 00702 Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->TauIn )); 00703 } 00704 fprintf( ioQQQ, "\n" ); 00705 } 00706 00707 fprintf( ioQQQ, "\n\n Aul As array\n" ); 00708 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][trace.ipIsoTrace[ipISO]]; ipHi++ ) 00709 { 00710 fprintf( ioQQQ, " %3ld", ipHi ); 00711 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00712 { 00713 if( Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul <= iso.SmallA ) 00714 continue; 00715 00716 fprintf( ioQQQ,PrintEfmt("%9.2e", 00717 Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul) ); 00718 } 00719 fprintf( ioQQQ, "\n" ); 00720 } 00721 00722 fprintf( ioQQQ, "\n\n Aul*esc array\n" ); 00723 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][trace.ipIsoTrace[ipISO]]; ipHi++ ) 00724 { 00725 fprintf( ioQQQ, " %3ld", ipHi ); 00726 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00727 { 00728 if( Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul <= iso.SmallA ) 00729 continue; 00730 00731 fprintf( ioQQQ,PrintEfmt("%9.2e", 00732 Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul* 00733 (Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Pdest + 00734 Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Pelec_esc + 00735 Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Pesc) )); 00736 } 00737 fprintf( ioQQQ, "\n" ); 00738 } 00739 00740 fprintf( ioQQQ, "\n\n H opac array\n" ); 00741 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][trace.ipIsoTrace[ipISO]]; ipHi++ ) 00742 { 00743 fprintf( ioQQQ, " %3ld", ipHi ); 00744 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00745 { 00746 if( Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul <= iso.SmallA ) 00747 continue; 00748 00749 fprintf( ioQQQ,PrintEfmt("%9.2e", 00750 Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->opacity )); 00751 } 00752 fprintf( ioQQQ, "\n" ); 00753 } 00754 } 00755 00756 else 00757 { 00758 fprintf( ioQQQ, " RT_tau_init called.\n" ); 00759 } 00760 } 00761 00762 ASSERT( Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauIn> 0. ); 00763 ASSERT( Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->PopOpc>= 0. ); 00764 return; 00765 }