00001 /* This file is part of Cloudy and is copyright (C)1978-2007 by Gary J. Ferland 00002 * For conditions of distribution and use see copyright notice in license.txt */ 00003 /*RT_line_static do line radiative transfer, 00004 * evaluates escape and destruction probability, 00005 * called by HydroPEsc, RT_line_all */ 00006 /* ND wind distinction is not here - that is done where optical depths are incremented 00007 * with line opacity - rt_line_one_tauinc and in line width for rad pressure, in 00008 * rt escprob */ 00009 #include "cddefines.h" 00010 #include "rfield.h" 00011 #include "doppvel.h" 00012 #include "dense.h" 00013 #include "opacity.h" 00014 #include "lines_service.h" 00015 #include "conv.h" 00016 #include "radius.h" 00017 #include "rt.h" 00018 00019 /*RT_line_static do line radiative transfer, 00020 * called by RT_line_all */ 00021 static void RT_line_static( 00022 /* the em line we will work on */ 00023 EmLine * t , 00024 /* when this is true do both esc and dest probs, 00025 * when falst, only dest probs */ 00026 bool lgDoEsc , 00027 /* this is option to not include line self shielding across this zone. 00028 * this can cause pump to depend on zone thickness, and leads to unstable 00029 * feedback in some models with the large H2 molecule, due to Solomon 00030 * process depending on zone thickness and level populations. */ 00031 bool lgShield_this_zone ); 00032 00033 /*RT_line_one do rt for emission line structure - calls RT_line_static or RT_line_wind */ 00034 void RT_line_one( 00035 /* the em line we will work on */ 00036 EmLine *t , 00037 /* when this is true do both esc and dest probs, 00038 * when falst, only dest probs */ 00039 bool lgDoEsc , 00040 /* if true then we want to update the fine opacity scale */ 00041 bool lgDoFine_opac_update , 00042 /* this is option to not include line self shielding across this zone. 00043 * this can cause pump to depend on zone thickness, and leads to unstable 00044 * feedback in some models with the large H2 molecule, due to Solomon 00045 * process depending on zone thickness and level populations. */ 00046 bool lgShield_this_zone ) 00047 { 00048 static long int nTau[100], 00049 n; 00050 00051 DEBUG_ENTRY( "RT_line_one()" ); 00052 00053 /*>>chng 06 apr 05, effects of simply returning upon call with no pop */ 00054 /* skip line transfer if requested with 'no line transfer' command, but never skip Lya */ 00055 if( (t->PopLo<=SMALLFLOAT && conv.nTotalIoniz) || 00056 (!rfield.lgDoLineTrans && (t->iRedisFun != ipLY_A) ) ) 00057 { 00058 DEBUG_EXIT( "RT_line_one()" ); 00059 return; 00060 } 00061 00062 { 00063 /* option to keep track of population values during calls, 00064 * print out data to make histogram */ 00065 /*@-redef@*/ 00066 enum {DEBUG_LOC=false}; 00067 /*@+redef@*/ 00068 if( DEBUG_LOC ) 00069 { 00070 if( nzone==0 ) 00071 { 00072 for(n=0; n<100; ++n ) 00073 nTau[n] = 0; 00074 } 00075 if( t->PopLo<=SMALLFLOAT ) 00076 n = 0; 00077 else 00078 n = (long)log10( t->PopLo )+37; 00079 n = MIN2( n , 99 ); 00080 n = MAX2( n , 0 ); 00081 /*if( n==37 && nzone > 5 ) 00082 { 00083 fprintf(ioQQQ,"HIT\n"); 00084 }*/ 00085 ++nTau[n]; 00086 if( nzone > 183 ) 00087 { 00088 for(n=0; n<100; ++n ) 00089 fprintf(ioQQQ,"%li\t%li\n", n , nTau[n] ); 00090 cdEXIT(1); 00091 } 00092 } 00093 } 00094 00095 /* always call single routine */ 00096 RT_line_static( 00097 /* the line we will work on */ 00098 t , 00099 /* when this is true do both esc and dest probs, 00100 * when falst, only dest probs */ 00101 lgDoEsc , 00102 lgShield_this_zone); 00103 00104 /* >>chng 03 feb 12, add this array */ 00105 /* define fine opacity fine grid fine mesh */ 00106 /* rfield.lgOpacityFine flag set flast with no fine opacities command */ 00107 if( lgDoFine_opac_update && t->ipFine >= 0 && t->PopOpc > 0. && 00108 t->ipFine<rfield.nfine && rfield.lgOpacityFine ) 00109 { 00110 long int i , ip_lo , ip_hi; 00111 00112 /* number of fine opacity cells corresponding to one doppler width for current 00113 * temperature, velocity field, and nuclear mass, 00114 * rfield.fine_opac_velocity_width is width per cell, cm/s */ 00115 float cells_wide_1x = DoppVel.doppler[t->nelem-1]/rfield.fine_opac_velocity_width; 00116 00117 /* line center opacity - want float since will add to fine opacity array, 00118 * which is float */ 00119 float opac_line = (float)t->PopOpc * t->opacity / DoppVel.doppler[t->nelem-1]; 00120 00121 /* this is effective optical depth to this point. Do not work with line if this product 00122 * is substantially less than unity */ 00123 double dTauEffec = opac_line*radius.depth_x_fillfac; 00124 00125 /* >>chng 04 mar 29, use dTauEffec to choose whether line is important */ 00126 /*if( t->TauIn < 0.0001 )*/ 00127 /*if( dTauEffec > 0.1 )*/ 00128 /* >>chng 05 feb 27, keep all opacities */ 00129 if( dTauEffec > SMALLFLOAT ) 00130 { 00131 /* >>chng 04 mar 29, rewrite to use symmetry between 00132 * two halves of line profile */ 00133 /* width of optically thick line, do 4x with exponential core, 00134 * must be at least one cell */ 00135 long int nCells_core = (long)(cells_wide_1x*4.+1.5f); 00136 float opacity; 00137 00138 /* core is symmetric - make sure both upper and lower bounds 00139 * are within continuum energy bin */ 00140 if( t->ipFine - nCells_core < 1 ) 00141 nCells_core = t->ipFine - 1; 00142 if( t->ipFine + nCells_core > rfield.nfine ) 00143 nCells_core = t->ipFine -rfield.nfine - 1; 00144 00145 /* want core to be at least one cell wide */ 00146 nCells_core = MAX2( 1 , nCells_core ); 00147 00148 /* line center itself, must not double count here */ 00149 rfield.fine_opac_zone[t->ipFine] += opac_line; 00150 for( i=1; i<nCells_core; ++i ) 00151 { 00152 /* square of distance from line center in units of doppler width */ 00153 double xsq = POW2(i/cells_wide_1x); 00154 00155 /* the line opacity at that point, includes doppler core and damping wings */ 00156 opacity = opac_line*(float)( 00157 exp( -xsq ) + t->damp/MAX2(1., (1.772f*xsq)) ); 00158 00159 rfield.fine_opac_zone[t->ipFine+i] += opacity; 00160 rfield.fine_opac_zone[t->ipFine-i] += opacity; 00161 } 00162 00163 /* include damping wings if optical depth is large */ 00164 /* >>chng 04 mar 29, rewrite to use symmetry in two halves of line */ 00165 if( dTauEffec*t->damp/9. > 0.1 ) 00166 { 00167 /* do damping wings if very optically thick 00168 * this is number of cells to extend the damping wings - cells_wide is one dop width 00169 * 56.419 is 1 / (0.01 * sqrt pi) ) */ 00170 /* tests with th85orion, stopping at half -h2 point, showed that 0.01 was 00171 * needed for outer edge, given the defin of dTauEffec */ 00172 float x = (float)sqrt( dTauEffec * t->damp *56.419 ) * cells_wide_1x; 00173 long int nCells_damp; 00174 /* >>chng 04 mar 24, add test on size of x, which can exceed 00175 * limits of long in extreme optical depths */ 00176 if( x<LONG_MAX ) 00177 { 00178 nCells_damp = (long)x; 00179 00180 if( t->ipFine-nCells_damp < 1 ) 00181 nCells_damp = t->ipFine-1; 00182 00183 if( t->ipFine+nCells_damp > rfield.nfine ) 00184 nCells_damp = rfield.nfine - t->ipFine-1; 00185 } 00186 else 00187 { 00188 /* x was too big, just set to extreme range, which is 00189 * number of cells to nearest boundary */ 00190 nCells_damp = MIN2( rfield.nfine-t->ipFine , t->ipFine )-1; 00191 } 00192 00193 ASSERT( nCells_core>0 ); 00194 for( i=nCells_core; i<nCells_damp; ++i ) 00195 { 00196 /* distance from line center in doppler units */ 00197 float xsq = POW2(i/cells_wide_1x); 00198 /* term in denominator is x^2 */ 00199 opacity = opac_line*t->damp/(1.772f*xsq ); 00200 rfield.fine_opac_zone[t->ipFine+i] += opacity; 00201 rfield.fine_opac_zone[t->ipFine-i] += opacity; 00202 } 00203 } 00204 } 00205 00206 else if( dTauEffec> 0.001 ) 00207 { 00208 /* optically thin line - use square topped line profile */ 00209 ip_lo = t->ipFine - (long)(cells_wide_1x+0.5f); 00210 ip_hi = t->ipFine + (long)(cells_wide_1x+0.5f); 00211 ip_lo = MAX2(1, ip_lo ); 00212 /* make sure at least one cell wide */ 00213 ip_hi = MAX2( ip_hi , t->ipFine+1); 00214 ip_hi = MIN2( rfield.nfine , ip_hi ); 00215 for( i=ip_lo; i<ip_hi; ++i ) 00216 { 00217 rfield.fine_opac_zone[i] += opac_line; 00218 } 00219 } 00220 } 00221 00222 DEBUG_EXIT( "RT_line_one()" ); 00223 00224 return; 00225 } 00226 00227 /*RT_line_static do line radiative transfer for static geometry */ 00228 static void RT_line_static( 00229 /* the em line we will work on */ 00230 EmLine * t , 00231 /* when this is true do both esc and dest probs, 00232 * when false, only dest probs */ 00233 bool lgDoEsc , 00234 /* this is option to not include line self shielding across this zone. 00235 * this can cause pump to depend on zone thickness, and leads to unstable 00236 * feedback in some models with the large H2 molecule, due to Solomon 00237 * process depending on zone thickness and level populations. */ 00238 bool lgShield_this_zone ) 00239 { 00240 bool lgGoodTau; 00241 long int nelem; 00242 int nRedis = -1; 00243 00244 DEBUG_ENTRY( "RT_line_static()" ); 00245 00246 /* this is the routine that computes much of the radiative transfer 00247 * physics for lines for static case. */ 00248 00249 /*ASSERT( t->IonStg < LIMELM+4 );*/ 00250 /*ASSERT( t->nelem-1 < LIMELM+4 );*/ 00251 00252 /* >>chng 01 aug 18, do not evaluate if no population */ 00253 /* molecules are evaluated here, these must be valid */ 00254 /* >>chng 01 oct 25, add test for lgDoEsc, this is to make sure 00255 * that rates are evaluated at start of a new iteration */ 00256 if( dense.xIonDense[ t->nelem-1][t->IonStg-1] == 0. && !lgDoEsc ) 00257 { 00258 /* zero population, return after setting everything with side effects */ 00259 t->Pesc = 1.f; 00260 00261 /* inward escaping fraction */ 00262 t->FracInwd = 0.5; 00263 00264 /* pumping rate */ 00265 t->pump = 0.; 00266 00267 /* destruction probability */ 00268 t->Pdest = 0.; 00269 t->Pelec_esc = 0.; 00270 00271 /* damping constant */ 00272 t->damp = 0.; 00273 00274 DEBUG_EXIT( "RT_line_static()" ); 00275 00276 return; 00277 } 00278 /* check that we don't have a runaway maser */ 00279 else if( t->TauIn < -30. ) 00280 { 00281 fprintf( ioQQQ, "PROBLEM RT_line_static called with large negative optical depth, zone %.2f, setting lgAbort true.\n", 00282 fnzone ); 00283 DumpLine(t); 00284 /* >>>chng 00 apr 23, return busted true instead of exit here, so that code will 00285 * not hang under MPI */ 00286 lgAbort = true; 00287 00288 DEBUG_EXIT( "RT_line_static()" ); 00289 return; 00290 } 00291 00292 /* this checks if we have overrun the optical depth scale 00293 * possible for line to mase, reason for sec test 00294 * >>chng 97 jan 08, on iterations where the first iteration had NO species, 00295 * but it does exist on later iterations, added test for existence of ion */ 00296 /* >>chng 03 may 14, go to function visible to whole code */ 00297 /* test on lgTauOutOn is to say tau ok if this is first iteration 00298 * and outward optical depths turned off */ 00299 if( !opac.lgTauOutOn || lgTauGood( t ) ) 00300 { 00301 lgGoodTau = true; 00302 } 00303 else 00304 { 00305 /* do not do anything if we have overrun the scale, */ 00306 lgGoodTau = false; 00307 } 00308 00309 /* atomic number of this element on the C scale*/ 00310 nelem = t->nelem-1; 00311 00312 /* damping constant for the line, save it */ 00313 t->damp = t->dampXvel / DoppVel.doppler[t->nelem-1]; 00314 /*if( t->damp <0. ) 00315 { 00316 fprintf(ioQQQ,"DEBUG hit it \n"); 00317 }*/ 00318 ASSERT( t->damp >= 0. ); 00319 00320 /* static solution - which type of line will determine 00321 * which redistribution function */ 00322 /* iRedisFun == 1 - alpha resonance line, partial redistribution, 00323 * ipPRD == 1 */ 00324 if( t->iRedisFun == ipPRD ) 00325 { 00326 /* incomplete redistribution with wings */ 00327 if( lgDoEsc ) 00328 { 00329 if( lgGoodTau ) 00330 { 00331 t->Pesc = (float)esc_PRD( t->TauIn , t->TauTot , t->damp ); 00332 00333 /* inward escaping fraction */ 00334 t->FracInwd = rt.fracin; 00335 } 00336 00337 # if 0 00338 /* update pumping even if over run optical depth scale */ 00339 if( rfield.lgInducProcess ) 00340 { 00341 /* the inward-looking escape probability, TauCon is optical depth to 00342 * the continuum source, including all overlapping lines*/ 00343 rt.wayin = (float)esc_PRD_1side(t->TauCon,t->damp); 00344 } 00345 # endif 00346 } 00347 nRedis = ipDEST_INCOM; 00348 } 00349 00350 /* complete redistribution without wings - t->ipLnRedis is ipCRD == -1 */ 00351 else if( t->iRedisFun == ipCRD ) 00352 { 00353 if( lgDoEsc ) 00354 { 00355 if( lgGoodTau ) 00356 { 00357 /* >>chng 01 mar -6, escsub will call any of several esc prob routines, 00358 * depending of how core is set. We always want core-only for this option, 00359 * so call esca0k2(tau) directly */ 00360 t->Pesc = (float)esc_CRDcore( t->TauIn , t->TauTot ); 00361 00362 /* inward escaping fraction */ 00363 t->FracInwd = rt.fracin; 00364 } 00365 00366 # if 0 00367 /* continuum pumping rate 00368 * the "no induced" command causes continuum pumping to be set to 0 */ 00369 if( rfield.lgInducProcess ) 00370 { 00371 /* wayin = esc_CRDwing_1side( t(ipLnTauCon),damp ) 00372 * this is to make it more like C90 was 00373 * complete redistribution with doppler core only */ 00374 rt.wayin = (float)esca0k2(t->TauCon); 00375 } 00376 # endif 00377 } 00378 nRedis = ipDEST_K2; 00379 } 00380 00381 /* chng 01 feb 27, add CRD with wings, = 2 */ 00382 else if( t->iRedisFun == ipCRDW ) 00383 { 00384 /* complete redistribution with damping wings */ 00385 if( lgDoEsc ) 00386 { 00387 if( lgGoodTau ) 00388 { 00389 t->Pesc = (float)esc_CRDwing( t->TauIn , t->TauTot , t->damp ); 00390 00391 /* inward escaping fraction */ 00392 t->FracInwd = rt.fracin; 00393 } 00394 # if 0 00395 /* continuum pumping rate 00396 * the "no induced" command causes continuum pumping to be set to 0 */ 00397 if( rfield.lgInducProcess ) 00398 { 00399 /* wayin = esc_CRDwing_1side( t(ipLnTauCon),damp ) 00400 * this is to make it more like C90 was 00401 * complete redistribution with doppler core only */ 00402 /* >>chng 01 may 29, had been esc_PRD_1side, as per above comment, 00403 * but this should include damping wings. */ 00404 /*rt.wayin = (float)esc_PRD_1side(t->TauCon,damp);*/ 00405 rt.wayin = (float)esc_CRDwing_1side(t->TauCon,t->damp); 00406 } 00407 # endif 00408 } 00409 00410 nRedis = ipDEST_K2; 00411 } 00412 00413 /* chng 03 may 14, special Lya case*/ 00414 else if( t->iRedisFun == ipLY_A ) 00415 { 00416 /* incomplete redistribution with wings, for special case of Lya 00417 * uses fits to Hummer & Kunasz numerical results 00418 * this routine is different because escape and dest probs 00419 * are evaluated together, so no test of lgDoEsc */ 00420 if( lgGoodTau ) 00421 { 00422 double dest , esin; 00423 00424 /* this will always evaluate escape prob, no matter what lgDoEsc is. 00425 * The destruction prob comes back as dest */ 00426 t->Pesc = (float)RTesc_lya(&esin , &dest , t->PopOpc , nelem); 00427 00428 /* this is current destruction rate */ 00429 t->Pdest = (float)dest; 00430 00431 /* this is fraction of line which is inward escaping */ 00432 t->FracInwd = rt.fracin; 00433 } 00434 # if 0 00435 /* update pumping even if over run optical depth scale */ 00436 if( rfield.lgInducProcess ) 00437 { 00438 /* the inward-looking escape probability */ 00439 rt.wayin = (float)esc_PRD_1side(t->TauCon,t->damp); 00440 } 00441 # endif 00442 } 00443 else 00444 { 00445 fprintf( ioQQQ, " RT_line_static called with redistribution function zero\n" ); 00446 ShowMe(); 00447 puts( "[Stop in RT_line_static]" ); 00448 cdEXIT(EXIT_FAILURE); 00449 } 00450 00451 /* pumping by incident and diffuse continua */ 00452 /* >>chng 03 may 15, this block had been repeated all four times above, move down 00453 * here to use common code, only rt.wayin is defined above */ 00454 /* >>chng 03 may 17, do not evaluate when drad is zero, before first 00455 * zone thickness is established */ 00456 /* >>chng 06 nov 29, reordered if clauses, no point in doing bigger loop if "no induced" is on. */ 00457 /* option to kill induced processes */ 00458 if( !rfield.lgInducProcess ) 00459 { 00460 t->pump = 0.; 00461 } 00462 else if( (lgDoEsc || (t->iRedisFun == ipLY_A) ) ) 00463 { 00464 float dTau; 00465 double pumpsave = t->pump; 00466 float shield_continuum; 00467 00468 /* >>chng 04 apr 18, break out line continuum shielding into this function */ 00469 shield_continuum = RT_continuum_shield_fcn( t ); 00470 00471 /* continuum upward pumping rate, A gu/gl abs prob occnum 00472 * the "no induced" command causes continuum pumping to be set to 0 */ 00473 /*t->pump = t->Aul * t->gHi / t->gLo * rt.wayin **/ 00474 /* >>chng 05 feb 16, add outward diffuse emission */ 00475 /*t->pump = t->Aul * t->gHi / t->gLo * shield_continuum * 00476 rfield.OccNumbIncidCont[t->ipCont-1];*/ 00477 t->pump = t->Aul * t->gHi / t->gLo * shield_continuum *( 00478 rfield.OccNumbIncidCont[t->ipCont-1] + rfield.OccNumbContEmitOut[t->ipCont-1] ); 00479 00480 /* >>chng 99 dec 02, add correction for line optical depth across zone */ 00481 /* correction for attenuation within line across zone */ 00482 /* >>chng 99 mar 18, from 0 to 1e-8 in following to avoid underflow, 00483 * this caused an overestimate of pumping rates */ 00484 /* >>chng 00 dec 19, nova.in oscillated since correction factor was 00485 * very large across zones that were already very optically thick, dr changed, 00486 * causing feedback between ionization and thickness 00487 * the intention of this correction was to only apply when on the 00488 * verge of optically thick, 00489 * change lower bound on product t->dTau * radius.drad_x_fillfac to 1e-3 to avoid 00490 * constantly getting unity, also check on relative change optical depth/optical depth*/ 00491 /* this can drive oscillations in very neutral gas, as blr92.in, blr89.in 00492 * or nova.in, when higher lyman lines have produce 00493 * t->dTau * radius.drad_x_fillfac about unity, where small changes in 00494 * pump can change ionization, dTau, and hence pump */ 00495 /*if( t->dTau * radius.drad_x_fillfac > 1e-8 )*/ 00496 /* >>chng 02 jun 14, nova.in again - upper limit had been 10 - but 00497 * there were major jumps in heating rate at this point - increase 10 -> 1e10 00498 * not clear why there should be an upper limit at all, other than to save time? */ 00499 /* update dTau since affects pump rate */ 00500 /* >>chng 03 jun 18, include option to not shield across this zone */ 00501 /* >>chng 03 jul 19, add continuum opacity */ 00502 dTau = (float)((t->PopOpc * t->opacity / DoppVel.doppler[nelem] + opac.opacity_abs[t->ipCont-1])* 00503 radius.drad_x_fillfac_mean); 00504 /* >>chng 04 mar 04, do not want to include maser action across this zone */ 00505 dTau = MIN2(1.f, dTau); 00506 /*if( 0 && lgShield_this_zone && dTau * radius.drad_x_fillfac > 1e-3 && dTau * radius.drad_x_fillfac < 1e10 )*/ 00507 /* >>chng 04 mar 04, use mean of past few zone thicknesses rather than 00508 * current zone thickness since this can induce wild oscillations, see 00509 * nova.in test case */ 00510 if( lgShield_this_zone && dTau > 1e-3 && dTau < 1e10 ) 00511 { 00512 /* during very first search dTau is -1e38*/ 00513 /* >>chng 04 mar 04, use mean, see above comment */ 00514 /*t->pump *= log(1. + dTau * radius.drad_x_fillfac ) / (dTau * radius.drad_x_fillfac);*/ 00515 t->pump *= log(1. + dTau ) / dTau; 00516 } 00517 00518 /* >>chng 00 Oct 03, add pumping by local diffuse continuum, 00519 * rfield.DiffPumpOn is unity unless process disabled by setting to 0 in parsedont.c 00520 * with no DIFFuse LINE PUMPing command */ 00521 if( dTau*rfield.DiffPumpOn > SMALLFLOAT ) 00522 { 00523 /* >>chng 02 mar 14, added brackets around scale*rfield.OccNumbDiffCont so 00524 * that smallest and largest number are multiplied first to prevent overflow, PvH */ 00525 double scale; 00526 if( opac.lgTauOutOn ) 00527 { 00528 if( lgGoodTau ) 00529 { 00530 /* >>chng 03 jul 19, add check for tau = 1 in continuum */ 00531 /* scale for tau = 1 in the continuum */ 00532 double scale_con = 1. / SDIV(opac.opacity_abs[t->ipCont-1]+opac.opacity_sct[t->ipCont-1]); 00533 /* inward looking depth, dTau is opacity in line, cm^-1, 00534 * so multiplying by this scale is same as dividing by opacity in 00535 * source function */ 00536 double scale_in = MIN3( 1./ dTau , radius.depth , scale_con); 00537 double scale_out = MIN3( 1./ dTau , radius.Depth2Go , scale_con ); 00538 t->pump += t->Aul * t->gHi / t->gLo * 0.5 *((scale_in+scale_out)* 00539 rfield.OccNumbDiffCont[t->ipCont-1]); 00540 } 00541 } 00542 else 00543 { 00544 /* >>chng 03 jul 19, add check for tau = 1 in continuum */ 00545 /* scale for tau = 1 in the continuum */ 00546 double scale_con = 1. / SDIV(opac.opacity_abs[t->ipCont-1]+opac.opacity_sct[t->ipCont-1]); 00547 /* first iteration, only inward looking depths really known */ 00548 scale = 1./ dTau; 00549 scale = MIN3( scale , radius.depth , scale_con ); 00550 t->pump += t->Aul * t->gHi / t->gLo * (scale* 00551 rfield.OccNumbDiffCont[t->ipCont-1]); 00552 } 00553 } 00554 /* this stop oscillations from developing in very deep blr models 00555 * like blr89 and blr98, but must not be applied too early since 00556 * line optical depths change greatly in first few zones */ 00557 if( nzone>50 ) 00558 t->pump = (pumpsave + t->pump) / 2.; 00559 } 00560 00561 00562 /* escape following scattering off an electron */ 00563 /* >>chng 02 may 08, put elec scat here, had been only for h-like species, 00564 * in hydropesc */ 00565 /* this is turned off with the no scattering escape command */ 00566 if( rt.lgElecScatEscape ) 00567 { 00568 /* in the EmLine structure PopOpc is the pop relative to the ion, 00569 * but has been converted to a physical population before this routine 00570 * was called. No need to convert here */ 00571 double opac_line = t->PopOpc * t->opacity/DoppVel.doppler[nelem]; 00572 00573 if( opac_line > SMALLFLOAT ) 00574 { 00575 /* the old escape probability */ 00576 double eesc_save = t->Pelec_esc; 00577 /* the opacity in electron scattering */ 00578 double es = dense.eden*6.65e-25; 00579 /* >>chng 04 may 12, use total scattering 00580 double es = opac.opacity_sct[t->ipCont-1]; */ 00581 /* the optical depth in electron scattering*/ 00582 double taues = es * radius.depth; 00583 /* this is equation 5 of 00584 *>>refer line desp Netzer, H., Elitzur, M., & Ferland, G. J. 1985, ApJ, 299, 752*/ 00585 /* >>chng 03 may 07, add correction for total optical depth to electrons */ 00586 double opacity_ratio = es/(es+opac_line) * (taues/(1.+taues)); 00587 /* >>chng 04 may 12, just do line to total opac 00588 double opacity_ratio = es/(es+opac_line);*/ 00589 t->Pelec_esc = (float)opacity_ratio * MAX2(0.f,1.f-t->Pesc-t->Pdest); 00590 /*KLUDGE >>chng 03 feb 02 00591 * take mean of old and new if electron scattering escape */ 00592 if( !conv.lgSearch ) 00593 t->Pelec_esc = t->Pelec_esc*0.75f + (float)eesc_save*0.25f; 00594 } 00595 else 00596 t->Pelec_esc = 0.; 00597 } 00598 00599 /* >>>chng 99 dec 18, did not have the test for lgGoodTau, and so clobbered 00600 * dest prob when overrun */ 00601 /* only do this if not Lya special case, since dest prob already done */ 00602 if( lgGoodTau && t->iRedisFun!=ipLY_A) 00603 { 00604 double DestSave = t->Pdest; 00605 t->Pdest = (float)RT_DestProb( 00606 /* the abundance of the species */ 00607 t->PopOpc, 00608 /* line center opacity in funny units (needs a vel) */ 00609 t->opacity, 00610 /* array index for line in continuum array, 00611 * on f not c scale */ 00612 t->ipCont, 00613 /* line width in velocity units */ 00614 DoppVel.doppler[nelem], 00615 /* escape probability */ 00616 t->Pesc, 00617 /* redistribution function */ 00618 nRedis); 00619 00620 /* >>chng 03 may 10, this had been 3/4 new + 1/4 old, but had to be combined 00621 ^ with half and half when this routine was called, to converge very high Z high U 00622 * quasar models. bring all these factors together here */ 00623 /*t->Pdest = 0.75f*t->Pdest + 0.25f * DestSave;*/ 00624 # define OLDFAC (0.625) 00625 if( nzone>1 ) 00626 t->Pdest = (float)((1.-OLDFAC)*t->Pdest + OLDFAC * DestSave); 00627 } 00628 00629 DEBUG_EXIT( "RT_line_static()" ); 00630 00631 return; 00632 } 00633 00634 #if 0 00635 /*RT_line_wind do line radiative transfer for wind geometry */ 00636 static void RT_line_wind(EmLine * t , bool lgDoEsc); 00637 /*RT_line_wind do line radiative transfer for wind geometry */ 00638 static void RT_line_wind(EmLine * t , bool lgDoEsc) 00639 { 00640 float velocity; 00641 00642 DEBUG_ENTRY( "RT_line_wind()" ); 00643 00644 /* 00645 * this is the routine that computes much of the radiative transfer 00646 * physics for lines in wind case. Lines diveded into two types, 00647 * high quality (level1) and lower quality (atom_level2) g-bar lines 00648 */ 00649 00650 /* >>chng 01 aug 18, do not evaluation if no population */ 00651 /* >>chng 01 oct 25, add test for lgDoEsc, this is to make sure 00652 * that rates are evaluated at start of a new iteration */ 00653 if( dense.xIonDense[ t->nelem-1][t->IonStg-1] == 0. && !lgDoEsc ) 00654 { 00655 /* zero population, return after setting everything with side effects */ 00656 t->Pesc = 1.f; 00657 00658 /* inward escaping fraction */ 00659 t->FracInwd = 0.5; 00660 00661 /* pumping rate */ 00662 t->pump = 0.; 00663 00664 /* destruction probability */ 00665 t->Pdest = 0.; 00666 t->Pelec_esc = 0.; 00667 00668 DEBUG_EXIT( "RT_line_wind()" ); 00669 return; 00670 } 00671 00672 /*confirm pops and cs ok */ 00673 ASSERT( t->ColOvTot <= 1. && t->ColOvTot >= 0. ); 00674 00675 /* t->nelem is atomic number of species */ 00676 ASSERT( t->nelem >= 1 && t->ipCont > 0 ); 00677 00678 /* line width, thermal + turbulence */ 00679 velocity = DoppVel.doppler[t->nelem-1]; 00680 00681 t->damp = 0.; 00682 00683 /* use only one-sided escape probabilities 00684 * get escape probability, this gets fractions too */ 00685 if( lgDoEsc ) 00686 { 00687 t->Pesc = esc_CRDwing_1side(t->TauIn,t->damp); 00688 00689 /* inward fraction */ 00690 t->FracInwd = 0.5; 00691 00692 /* continuum pumping */ 00693 if( rfield.lgInducProcess ) 00694 { 00695 /* >>chng 05 feb 16, include diffuse outward emission */ 00696 /*t->pump = t->Aul*t->gHi/t->gLo*t->Pesc* 00697 rfield.OccNumbIncidCont[t->ipCont-1];*/ 00698 t->pump = t->Aul*t->gHi/t->gLo*t->Pesc*( 00699 rfield.OccNumbIncidCont[t->ipCont-1] + rfield.OccNumbContEmitOut[t->ipCont-1]); 00700 } 00701 else 00702 { 00703 t->pump = 0.; 00704 } 00705 } 00706 00707 /* get destruction probability */ 00708 t->Pdest = 00709 RT_DestProb( 00710 t->PopOpc, 00711 /* its line absorption cross section */ 00712 t->opacity, 00713 /* index for line in continuum array, on f not c scale */ 00714 t->ipCont, 00715 velocity, 00716 t->Pesc, 00717 ipDEST_K2); 00720 DEBUG_EXIT( "RT_line_wind()" ); 00721 00722 return; 00723 } 00724 #endif 00725 00726