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_line_all do escape and destruction probabilities for all lines in code. 00004 * Called with false arg by ionize, to only redo destruction probabilities, 00005 * and with true by cloudy to do both escape and destruction */ 00006 #include "cddefines.h" 00007 #include "taulines.h" 00008 #include "atomfeii.h" 00009 #include "dense.h" 00010 #include "conv.h" 00011 #include "atoms.h" 00012 #include "rfield.h" 00013 #include "wind.h" 00014 #include "iso.h" 00015 #include "h2.h" 00016 #include "opacity.h" 00017 #include "trace.h" 00018 #include "lines_service.h" 00019 #include "atmdat.h" 00020 #include "hydrogenic.h" 00021 #include "rt.h" 00022 00023 void RT_line_all( 00024 /* this is true if we want to do both escape and destruction probabilities, 00025 * and false if only destruction probabilities are needed */ 00026 bool lgDoEsc , 00027 /* flag saying whether to update fine opacities */ 00028 bool lgUpdateFineOpac ) 00029 { 00030 long int i, 00031 ion, 00032 ipISO, 00033 nelem; 00034 long ipHi , ipLo; 00035 double factor, 00036 coloi; 00037 double SaveLyaPesc[NISO][LIMELM] , 00038 SaveLyaPdest[NISO][LIMELM]; 00039 /* this flag says whether to update the line RT */ 00040 bool lgRT_update = lgUpdateFineOpac || lgDoEsc || conv.lgIonStageTrimed; 00041 00042 DEBUG_ENTRY( "RT_line_all()" ); 00043 00044 if( trace.lgTrace ) 00045 fprintf( ioQQQ, " RT_line_all called\n" ); 00046 00047 /* skip line transfer if requested, and not first call in this zone */ 00048 if( !rfield.lgDoLineTrans && conv.nPres2Ioniz ) 00049 { 00050 return; 00051 } 00052 00053 /* this array is huge and takes significant time to zero out or update, 00054 * only do so when needed, */ 00055 rfield.lgFine_opac_update = lgUpdateFineOpac; 00056 if( rfield.lgFine_opac_update ) 00057 { 00058 /* zero out fine opacity array */ 00059 memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine*sizeof(realnum) ); 00060 00061 /* this is offset within fine opacity array caused by Doppler shift 00062 * between first zone, the standard of rest, and current position 00063 * in case of acceleration away from star in outflowing wind 00064 * dWind is positive, this means that the rest frame original 00065 * velocity is red shifted to lower energy */ 00066 realnum dWind = wind.windv - wind.windv0; 00067 /* will add ipVelShift to index of original energy so redshift has 00068 * to be negative */ 00069 rfield.ipFineConVelShift = -(long int)(dWind / rfield.fine_opac_velocity_width+0.5); 00070 } 00071 00072 # if 0 00073 /* this code is a copy of the code within line_one that does cloaking 00074 * within this zone. it can be used to see how a particular line 00075 * is being treated. */ 00076 { 00077 #include "doppvel.h" 00078 double dTau , aa; 00079 ipISO = 0; nelem = 0;ipLo = 0; 00080 ipHi = 23; 00081 dTau = Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc * 00082 Transitions[ipISO][nelem][ipHi][ipLo].Emis->opacity / 00083 DoppVel.doppler[nelem] + opac.opacity_abs[Transitions[ipISO][nelem][ipHi][ipLo].ipCont-1]; 00084 dTau *= radius.drad_x_fillfac_mean; 00085 aa = log(1. + dTau ) / SDIV(dTau); 00086 fprintf(ioQQQ,"DEBUG dTau\t%.2f\t%.5e\t%.5e\t%.5e\n",fnzone,dTau, 00087 radius.drad_x_fillfac_mean, 00088 aa); 00089 } 00090 # endif 00091 00092 /* find Stark escape probabilities for hydrogen itself */ 00093 RT_stark(); 00094 00095 /*if( lgUpdateFineOpac ) 00096 fprintf(ioQQQ,"debuggg\tlgUpdateFineOpac in rt_line_all\n");*/ 00097 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00098 { 00099 /* loop over all iso-electronic sequences */ 00100 for( nelem=ipISO; nelem<LIMELM; ++nelem ) 00101 { 00102 /* parent ion stage, for H is 1, for He is 1 for He-like and 00103 * 2 for H-like */ 00104 ion = nelem+1-ipISO; 00105 00106 /* element turned off */ 00107 if( !dense.lgElmtOn[nelem] ) 00108 continue; 00109 /* need we consider this ion? */ 00110 if( ion <=dense.IonHigh[nelem] ) 00111 { 00112 /* save escape and destruction probs for each Lya since 00113 * these are special cases, and quantities are not damped 00114 * in routine that evaluates each, test is to avoid 00115 * evaluating iso sequence that does not exist for an element, 00116 * as in He-like hydrogen */ 00117 if( ipISO<=nelem ) 00118 { 00119 SaveLyaPesc[ipISO][nelem] = Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc; 00120 SaveLyaPdest[ipISO][nelem] = Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pdest; 00121 } 00122 00123 /* convert pops to per unit vol rather than per ion */ 00124 if( dense.xIonDense[nelem][ion] > 1e-30 ) 00125 { 00126 factor = dense.xIonDense[nelem][ion]; 00127 } 00128 else 00129 { 00130 /* case where almost no parent ion - this will make 00131 * very large line opacity, so background dest small */ 00132 factor = 1.; 00133 } 00134 00135 /* loop over all lines */ 00136 for( ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ++ipHi ) 00137 { 00138 for( ipLo=0; ipLo < ipHi; ++ipLo ) 00139 { 00140 /* negative ipCont means this is not a real line, so do not 00141 * transfer it */ 00142 if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont < 0 ) 00143 continue; 00144 00145 double SavePopOpc = Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc; 00146 /* must temporarily make ipLnPopOpc physical */ 00147 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc *= factor; 00148 /*if( ipISO==0 && nelem==0 && ipHi==12 && ipLo==11 ) 00149 fprintf(ioQQQ,"DEBUG rt loop %li %li %li %li %.3e %.3e %.3e\n", 00150 ipISO, nelem, ipHi , ipLo , 00151 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc , 00152 Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn, 00153 Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauTot);*/ 00154 00155 /* generate escape prob, pumping rate, destruction prob, 00156 * inward outward fracs */ 00157 RT_line_one(&Transitions[ipISO][nelem][ipHi][ipLo] , lgDoEsc , lgUpdateFineOpac,true); 00158 00159 00160 { 00161 /* set true to print pump rates*/ 00162 enum {DEBUG_LOC=false}; 00163 if( DEBUG_LOC && nelem==1&& ipLo==0 /*&& iteration==2*/ ) 00164 { 00165 fprintf(ioQQQ,"DEBUG pdest\t%3li\t%.2f\t%.3e\n", 00166 ipHi , 00167 fnzone, 00168 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pdest); 00169 } 00170 } 00171 00172 /* go back to original units so that final correction ok */ 00173 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc = SavePopOpc; 00174 } 00175 } 00176 00177 if( ipISO > 0 && lgDoEsc ) 00178 { 00179 /* do not let too much Lya escape outward since so important 00180 * H Lya is special case done above */ 00181 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc *= 00182 opac.ExpmTau[Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].ipCont-1]; 00183 } 00184 00185 /* now update two photon induced rates */ 00186 atmdat_2phot_rate(ipISO , nelem); 00187 00188 /* this is option to not do destruction probabilities 00189 * for case b no pdest option */ 00190 if( opac.lgCaseB_no_pdest ) 00191 { 00192 ipLo = 0; 00193 /* okay to let this go to numLevels_max. */ 00194 for( ipHi=ipLo+1; ipHi < iso.numLevels_max[ipISO][nelem]; ++ipHi ) 00195 { 00196 if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 ) 00197 continue; 00198 00199 /* hose the previously computed destruction probability */ 00200 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pdest = SMALLFLOAT; 00201 } 00202 } 00203 00204 ipLo = 0; 00205 /* these are the extra lyman lines for the iso sequences */ 00206 /*for( ipHi=2; ipHi < iso.nLyman[ipISO]; ipHi++ )*/ 00207 /* only update if significant abundance and need to update fine opac */ 00208 if( dense.xIonDense[nelem][ion] > 1e-30 && lgUpdateFineOpac ) 00209 { 00210 /* we just want the population of the ground state */ 00211 factor = StatesElem[ipISO][nelem][0].Pop * dense.xIonDense[nelem][ion]; 00212 for( ipHi=StatesElem[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1].n; ipHi < iso.nLyman[ipISO]; ipHi++ ) 00213 { 00214 /* must make ipLnPopOpc physical */ 00215 ExtraLymanLines[ipISO][nelem][ipHi].Emis->PopOpc = factor; 00216 ExtraLymanLines[ipISO][nelem][ipHi].Lo->Pop = 00217 StatesElem[ipISO][nelem][ipLo].Pop; 00218 00219 /* actually do the work */ 00220 RT_line_one(&ExtraLymanLines[ipISO][nelem][ipHi] , lgDoEsc , lgUpdateFineOpac,true); 00221 00222 /* reset to funny population of ground state */ 00223 ExtraLymanLines[ipISO][nelem][ipHi].Emis->PopOpc = 00224 StatesElem[ipISO][nelem][0].Pop; 00225 } 00226 } 00227 00228 }/* if nelem if ion <=dense.IonHigh */ 00229 }/* loop over nelem */ 00230 }/* loop over ipISO */ 00231 00232 00233 /* add on Stark escape probabilities for H itself */ 00234 for( ipISO=0; ipISO<NISO; ipISO++ ) 00235 { 00236 for( nelem=ipISO; nelem<LIMELM; nelem++ ) 00237 { 00238 if( nelem < 2 || dense.lgElmtOn[nelem] ) 00239 { 00240 /* note that pesc and dest are updated no matter what escprob logic we 00241 * specify, and that these are not updated if we have overrun the 00242 * optical depth scale. only update here in that case */ 00243 if( lgTauGood( &Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0] ) ) 00244 { 00245 if( lgDoEsc ) 00246 { 00247 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00248 { 00249 for( ipLo=0; ipLo < ipHi; ipLo++ ) 00250 { 00251 if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 ) 00252 continue; 00253 00254 /* >>chng 03 jun 07, do not clobber esp prob when line is masing - 00255 * this had effect of preventing total escape prob from getting larger than 1 */ 00256 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc<1. ) 00257 { 00258 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc = MIN2(1.f, 00259 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc+ 00260 (realnum)iso.pestrk[ipISO][nelem][ipLo][ipHi]); 00261 } 00262 } 00263 } 00264 } 00265 /* always add on stark broadening term to Lya*/ 00266 else if( Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc < 1. ) 00267 { 00268 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc = MIN2(1.f, 00269 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc+ 00270 (realnum)iso.pestrk[ipISO][nelem][0][iso.nLyaLevel[ipISO]]); 00271 } 00272 } 00273 } 00274 } 00275 } 00276 00277 /* note that pesc and dest are updated no matter what escprob logic we 00278 * specify, and that these are not updated if we have overrun the 00279 * optical depth scale. only update here in that case */ 00280 if( lgTauGood( &Transitions[ipH_LIKE][ipHYDROGEN][iso.nLyaLevel[ipH_LIKE]][0] ) ) 00281 { 00282 /* do the 8446 problem */ 00283 atom_oi_calc(&coloi); 00284 00285 Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc = (realnum)(atoms.pmph31/ 00286 Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul); 00287 00288 if( trace.lgTrace && trace.lgIsoTraceFull[ipH_LIKE] ) 00289 { 00290 fprintf( ioQQQ, " RT_line_all calls P8446 who found pmph31=%10.2e\n", 00291 atoms.pmph31 ); 00292 } 00293 00294 /* add on destruction of hydrogen Lya by FeII 00295 * now add in FeII deexcitation via overlap, 00296 * but only as loss, not photoionization, source 00297 * dstfe2lya is Ly-alpha deexcit by FeII overlap - conversion into FeII em */ 00298 /* find FeII overlap destruction rate, 00299 * this does NOTHING when large FeII atom is turned on, 00300 * in this case evaluation is done in call to FeIILevelPops */ 00301 t_fe2ovr_la::Inst().atoms_fe2ovr(); 00302 /*fprintf(ioQQQ,"DEBUG fe2 %.2e %.2e\n", hydro.dstfe2lya , 00303 hydro.HLineWidth);*/ 00304 00305 /* >>chng 00 jan 06, let dest be large than 1 to desaturate the atom */ 00306 /* >>chng 01 apr 01, add test for tout >= 0., 00307 * must not add to Pdest when it has not been refreshed here */ 00308 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest += hydro.dstfe2lya; 00309 } 00310 00311 /* take mean of old and new to damp oscillations in Lya (only) */ 00312 if( nzone > 1 ) 00313 { 00314 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00315 { 00316 /* loop over all iso-electronic sequences */ 00317 for( nelem=ipISO; nelem<LIMELM; ++nelem ) 00318 { 00319 ion = nelem+1-ipISO; 00320 if( ion <=dense.IonHigh[nelem] && ipISO<=nelem ) 00321 { 00322 /* now take mean of old and new escape dest probs for Lya */ 00323 /* only ipLY_A redist did not already take average */ 00324 if( Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->iRedisFun==ipLY_A && 00325 lgTauGood( &Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0] ) ) 00326 { 00327 /* the Lya routine is different in that escape and dest 00328 * probs are always updated, unless we have overrun the 00329 * optical depth scale */ 00330 /*lint -e644 SaveLyaPdest not init */ 00331 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pdest = 00332 ((realnum)SaveLyaPdest[ipISO][nelem] + 00333 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pdest) / 2.f; 00334 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc = 00335 ((realnum)SaveLyaPesc[ipISO][nelem] + 00336 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc) / 2.f; 00337 /*lint +e644 SaveLyaPdest not init */ 00338 } 00339 } 00340 } 00341 } 00342 } 00343 /*5986 fprintf(ioQQQ,"DEBUG he rt_all\t%li\t%li\t%.4e\t%.4e\n", 00344 iteration, nzone, 00345 Transitions[1][1][10][5].Emis->TauIn, 00346 Transitions[1][1][10][5].Emis->Pesc);*/ 00347 00348 /* is continuum pumping of H Lyman lines included? yes, but turned off 00349 * with atom h-like Lyman pumping off command */ 00350 if( !hydro.lgLymanPumping ) 00351 { 00352 ipISO = ipH_LIKE; 00353 nelem = ipHYDROGEN; 00354 ipLo = 0; 00355 for( ipHi=ipLo+1; ipHi < iso.numLevels_max[ipISO][nelem]; ++ipHi ) 00356 { 00357 /* negative ipCont means this is not a real line, so do not 00358 * transfer it */ 00359 Transitions[ipISO][nelem][ipHi][ipLo].Emis->pump = 0.; 00360 } 00361 } 00362 00363 { 00364 /* following should be set true to print ots contributors for he-like lines*/ 00365 /*@-redef@*/ 00366 enum {DEBUG_LOC=false}; 00367 /*@+redef@*/ 00368 if( DEBUG_LOC && nzone>433 /*&& iteration==2*/ ) 00369 { 00370 /* option to dump a line */ 00371 DumpLine(&Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][0] ); 00372 } 00373 } 00374 00375 /* level 1 lines */ 00376 for( i=1; i <= nLevel1; i++ ) 00377 { 00378 RT_line_one(&TauLines[i] , lgDoEsc , lgUpdateFineOpac,true); 00379 } 00380 /* co carbon monoxide */ 00381 for( i=0; i < nCORotate; i++ ) 00382 { 00383 RT_line_one(&C12O16Rotate[i] , lgDoEsc , lgUpdateFineOpac,true); 00384 RT_line_one(&C13O16Rotate[i] , lgDoEsc , lgUpdateFineOpac,true); 00385 } 00386 for( i=0; i < nHFLines; i++ ) 00387 { 00388 RT_line_one(&HFLines[i] , lgDoEsc , lgUpdateFineOpac,true); 00389 } 00390 /* Database Atoms & Molecules*/ 00391 for( i=0; i < linesAdded2; i++) 00392 { 00393 RT_line_one(atmolEmis[i].tran , lgDoEsc , lgUpdateFineOpac,true); 00394 } 00395 /* this is a major time sink for this routine */ 00396 if( lgRT_update ) 00397 { 00398 for( i=0; i < nUTA; i++ ) 00399 { 00400 /* Aul = -1 means this line is superceded by better calculations in 00401 * different data set */ 00402 if( UTALines[i].Emis->Aul > 0. ) 00403 { 00404 /* these are not defined in cooling routines so must be set here */ 00405 UTALines[i].Emis->PopOpc = dense.xIonDense[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1]; 00406 UTALines[i].Lo->Pop = dense.xIonDense[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1]; 00407 UTALines[i].Hi->Pop = 0.; 00408 RT_line_one(&UTALines[i] , lgDoEsc , lgUpdateFineOpac,true); 00409 } 00410 } 00411 } 00412 00413 /* h-like ions only have one electron, no satellite lines. */ 00414 for( ipISO=ipHE_LIKE; ipISO<NISO; ++ipISO ) 00415 { 00416 /* loop over all iso-electronic sequences */ 00417 for( nelem=ipISO; nelem<LIMELM; ++nelem ) 00418 { 00419 if( dense.lgElmtOn[nelem] ) 00420 { 00421 for( ipLo=1; ipLo < (iso.numLevels_local[ipISO][nelem]-1); ++ipLo ) 00422 { 00423 RT_line_one(&SatelliteLines[ipISO][nelem][ipLo], lgDoEsc, lgUpdateFineOpac, true); 00424 } 00425 } 00426 } 00427 } 00428 00429 /* >>chng 03 aug 22, must call RT_line_one every single time, since that routine 00430 * establishes the fine opacities for the level 2 lines */ 00431 /* >>chng 01 aug 11, the level 2 lines were only updated when lgDoEsc was true, 00432 * this meant that dest probs were not updated but once, causing very high-Z models 00433 * to develop numerical oscialltions. removed if, now evaluate level 2 always */ 00434 /* only update their dest/esc probs one time per zone */ 00435 /* lgLevel2_OTS_Imp set true in dimacool if ots rates were significant */ 00436 /* level 2 heavy element lines in cooling with only g-bar,*/ 00437 for( i=0; i < nWindLine; i++ ) 00438 { 00439 /* >>chng 02 sug 11, do not include ions done in iso-sequences */ 00440 /*if( TauLine2[i].IonStg <= TauLine2[i].nelem-1 )*/ 00441 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 00442 { 00443 RT_line_one(&TauLine2[i] , lgDoEsc , lgUpdateFineOpac,true); 00444 } 00445 } 00446 00447 /* the large H2 molecule */ 00448 H2_RTMake( lgDoEsc , lgUpdateFineOpac); 00449 00450 /* The large model FeII atom */ 00451 FeII_RT_Make( lgDoEsc , lgUpdateFineOpac); 00452 return; 00453 }