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 /*esc_CRDwing_1side fundamental escape probability radiative transfer routine, for complete redistribution */ 00004 /*esc_PRD_1side fundamental escape probability radiative transfer routine for incomplete redistribution */ 00005 /*RTesc_lya escape prob for hydrogen atom Lya, using Hummer and Kunasz results, 00006 * called by hydropesc */ 00007 /*esc_PRD escape probability radiative transfer for incomplete redistribution */ 00008 /*esc_CRDwing escape probability for CRD with wings */ 00009 /*esc_CRDcore escape probability for CRD with no wings */ 00010 /*esca0k2 derive Hummer's K2 escape probability for Doppler core only */ 00011 /*RT_DestProb returns line destruction probability due to continuum opacity */ 00012 /*RT_LineWidth determine half width of any line with known optical depths */ 00013 #include "cddefines.h" 00014 #include "physconst.h" 00015 #define SCALE 2. 00016 #include "dense.h" 00017 #include "conv.h" 00018 #include "rfield.h" 00019 #include "opacity.h" 00020 #include "lines_service.h" 00021 #include "taulines.h" 00022 #include "doppvel.h" 00023 #include "pressure.h" 00024 #include "wind.h" 00025 #include "rt.h" 00026 /* 00027 *variables used to pass continuum optical depth and temp to 00028 *routine that integrates over continuum 00029 */ 00030 static double chnukt_ContTkt, chnukt_ctau; 00031 00032 /*escmase escape probability for negative (masing) optical depths,*/ 00033 STATIC double escmase(double tau); 00034 /*RTesc_lya_1side fit Hummer and Kunasz escape probability for hydrogen atom Lya */ 00035 STATIC void RTesc_lya_1side(double taume, 00036 double beta, 00037 realnum *esc, 00038 realnum *dest, 00039 /* position of line in frequency array on c scale */ 00040 long ipLine ); 00041 00042 double esc_PRD_1side(double tau, 00043 double a) 00044 { 00045 double atau, 00046 b, 00047 escinc_v; 00048 00049 DEBUG_ENTRY( "esc_PRD_1side()" ); 00050 ASSERT( a>0.0 ); 00051 00052 /* this is one of the three fundamental escape probability routines 00053 * the three are esc_CRDwing_1side, esc_PRD_1side, and RTesc_lya 00054 * it computes esc prob for incomplete redistribution 00055 * */ 00056 # if 0 00057 if( strcmp(rt.chEscFunSubord,"SIMP") == 0 ) 00058 { 00059 /* this set with "escape simple" command, used for debugging */ 00060 escinc_v = 1./(1. + tau); 00061 return( escinc_v ); 00062 } 00063 # endif 00064 00065 if( tau < 0. ) 00066 { 00067 /* line mased */ 00068 escinc_v = escmase(tau); 00069 } 00070 else if( tau < 10. ) 00071 { 00072 /* linear part of doppler core */ 00073 escinc_v = 1./(1. + 1.6*tau); 00074 } 00075 else 00076 { 00077 /* first find coefficient b(tau) */ 00078 atau = a*tau; 00079 if( atau > 1. ) 00080 { 00081 b = 1.6 + (3.*pow(2.*a,-0.12))/(1. + atau); 00082 } 00083 else 00084 { 00085 b = 1.6 + (3.*pow(2.*a,-0.12))/(1. + 1./sqrt(atau)); 00086 } 00087 b = MIN2(6.,b); 00088 00089 escinc_v = 1./(1. + b*tau); 00090 } 00091 return( escinc_v ); 00092 } 00093 00094 /*esc_CRDwing_1side fundamental escape probability radiative transfer routine, for complete redistribution */ 00095 double esc_CRDwing_1side(double tau, 00096 double a ) 00097 { 00098 double esccom_v; 00099 00100 DEBUG_ENTRY( "esc_CRDwing_1side()" ); 00101 00102 /* this is one of the three fundamental escape probability routines 00103 * the three are esc_CRDwing_1side, esc_PRD_1side, and RTesc_lya 00104 * it computes esc prob for complete redistribution with wings 00105 * computes escape prob for complete redistribution in one direction 00106 * */ 00107 00108 /* this is the only case that this routine computes, 00109 * and is the usual case for subordinate lines, 00110 * complete redistribution with damping wings */ 00111 esccom_v = esca0k2(tau); 00112 if( tau > 1e3 ) 00113 { 00114 esccom_v += 0.333*sqrt(a/(SQRTPI*tau)); 00115 } 00116 return( esccom_v ); 00117 } 00118 00119 /*RTesc_lya escape prob for hydrogen atom Lya, using 00120 >>refer La escp Hummer, D.G., & Kunasz, P.B., 1980, ApJ, 236, 609 00121 * called by hydropesc, return value is escape probability */ 00122 double RTesc_lya( 00123 /* the inward escape probability */ 00124 double *esin, 00125 /* the destruction probility */ 00126 double *dest, 00127 /* abundance of the species */ 00128 double abund, 00129 /* element number, 0 for H */ 00130 long int nelem) 00131 { 00132 double beta, 00133 conopc, 00134 escla_v; 00135 realnum dstin, 00136 dstout; 00137 00138 DEBUG_ENTRY( "RTesc_lya()" ); 00139 00140 /* 00141 * this is one of the three fundamental escape probability functions 00142 * the three are esc_CRDwing_1side, esc_PRD_1side, and RTesc_lya 00143 * evaluate esc prob for LA 00144 * optical depth in outer direction always defined 00145 */ 00146 00147 /* check charge */ 00148 ASSERT( nelem >= 0 && nelem < LIMELM ); 00149 00150 if( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot - 00151 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn < 0. ) 00152 { 00153 /* this is the case if we overrun the optical depth scale 00154 * just leave things as they are */ 00155 escla_v = Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pesc; 00156 rt.fracin = Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->FracInwd; 00157 *esin = rt.fracin; 00158 *dest = Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pdest; 00159 return( escla_v ); 00160 } 00161 00162 /* incomplete redistribution */ 00163 conopc = opac.opacity_abs[Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].ipCont-1]; 00164 if( abund > 0. ) 00165 { 00166 /* the continuous opacity is positive, we have a valid soln */ 00167 beta = conopc/(abund/SQRTPI*Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity/ 00168 DoppVel.doppler[nelem] + conopc); 00169 } 00170 else 00171 { 00172 /* abundance is zero, set miniumum dest prob */ 00173 beta = 1e-10; 00174 } 00175 00176 /* find rt.wayin, the escape prob in inward direction */ 00177 RTesc_lya_1side( 00178 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn, 00179 beta, 00180 &rt.wayin, 00181 &dstin , 00182 /* position of line in energy array on C scale */ 00183 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].ipCont-1); 00184 00185 ASSERT( (rt.wayin <= 1.) && (rt.wayin >= 0.) && (dstin <= 1.) && (dstin >= 0.) ); 00186 00187 /* find rt.wayout, the escape prob in outward direction */ 00188 RTesc_lya_1side(MAX2(Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot/100., 00189 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot-Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn), 00190 beta, 00191 &rt.wayout, 00192 &dstout, 00193 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].ipCont-1); 00194 00195 ASSERT( (rt.wayout <= 1.) && (rt.wayout >= 0.) && (dstout <= 1.) && (dstout >= 0.) ); 00196 00197 /* esc prob is mean of in and out */ 00198 escla_v = (rt.wayin + rt.wayout)/2.; 00199 /* the inward escaping part of the line */ 00200 *esin = rt.wayin; 00201 00202 /* dest prob is mean of in and out */ 00203 *dest = (dstin + dstout)/2.f; 00204 /* >>chng 02 oct 02, sum of escape and dest prob must be less then unity, 00205 * for very thin models this forces dest prob to go to zero, 00206 * rather than the value of DEST0, caught by Jon Slavin */ 00207 *dest = (realnum)MIN2( *dest , 1.-escla_v ); 00208 /* but dest prob can't be negative */ 00209 *dest = (realnum)MAX2(0., *dest ); 00210 00211 /* fraction of line emitted in inward direction */ 00212 rt.fracin = rt.wayin/(rt.wayin + rt.wayout); 00213 ASSERT( escla_v >=0. && *dest>=0. && *esin>=0. ); 00214 return( escla_v ); 00215 } 00216 00217 /*esc_PRD escape probability radiative transfer for incomplete redistribution */ 00218 double esc_PRD(double tau, 00219 double tau_out, 00220 double damp ) 00221 { 00222 double escgrd_v, 00223 tt; 00224 00225 DEBUG_ENTRY( "esc_PRD()" ); 00226 00227 /* find escape prob for incomp redis, average of two 1-sided probs*/ 00228 00229 if( iteration > 1 ) 00230 { 00231 /* outward optical depth if defined */ 00232 tt = tau_out - tau; 00233 /* help convergence by not letting tau go to zero at back edge of 00234 * when there was a bad guess for the total optical depth 00235 * note that this test is seldom hit since RTMakeStat does check 00236 * for overrun */ 00237 if( tt < 0. ) 00238 { 00239 tt = tau/SCALE; 00240 } 00241 00242 rt.wayin = (realnum)esc_PRD_1side(tau,damp); 00243 rt.wayout = (realnum)esc_PRD_1side(tt,damp); 00244 rt.fracin = rt.wayin/(rt.wayin + rt.wayout); 00245 escgrd_v = 0.5*(rt.wayin + rt.wayout); 00246 } 00247 else 00248 { 00249 /* outward optical depth not defined, dont estimate fraction out */ 00250 rt.fracin = 0.; 00251 rt.wayout = 1.; 00252 escgrd_v = esc_PRD_1side(tau,damp); 00253 rt.wayin = (realnum)escgrd_v; 00254 } 00255 00256 ASSERT( escgrd_v > 0. ); 00257 return( escgrd_v ); 00258 } 00259 00260 /*esc_CRDwing escape probability radiative transfer for CRDS in core only */ 00261 double esc_CRDwing(double tau_in, 00262 double tau_out, 00263 double damp) 00264 { 00265 double escgrd_v, 00266 tt; 00267 00268 DEBUG_ENTRY( "esc_CRDwing()" ); 00269 00270 /* find escape prob for CRD with damping wings, average of two 1-sided probs*/ 00271 00272 /* crd with wings */ 00273 if( iteration > 1 ) 00274 { 00275 /* outward optical depth if defined */ 00276 /* >>chng 03 jun 07, add test for masers here */ 00277 if( tau_out <0 || tau_in < 0. ) 00278 { 00279 /* we have a maser, use smallest optical depth to damp it out */ 00280 tt = MIN2( tau_out , tau_in ); 00281 tau_in = tt; 00282 } 00283 else 00284 { 00285 tt = tau_out - tau_in; 00286 /* help convergence by not letting tau_in go to zero at back edge of 00287 * when there was a bad guess for the total optical depth 00288 * note that this test is seldom hit since RTMakeStat does check 00289 * for overrun */ 00290 if( tt < 0. ) 00291 { 00292 tt = tau_in/SCALE; 00293 } 00294 } 00295 00296 rt.wayin = (realnum)esc_CRDwing_1side(tau_in,damp); 00297 rt.wayout = (realnum)esc_CRDwing_1side(tt,damp); 00298 rt.fracin = rt.wayin/(rt.wayin + rt.wayout); 00299 escgrd_v = 0.5*(rt.wayin + rt.wayout); 00300 } 00301 else 00302 { 00303 /* outward optical depth not defined, dont estimate fraction out */ 00304 rt.fracin = 0.; 00305 rt.wayout = 1.; 00306 escgrd_v = esc_CRDwing_1side(tau_in,damp); 00307 rt.wayin = (realnum)escgrd_v; 00308 } 00309 00310 ASSERT( escgrd_v > 0. ); 00311 return( escgrd_v ); 00312 } 00313 00314 /*esc_CRDwing escape probability radiative transfer for incomplete redistribution */ 00315 double esc_CRDcore(double tau_in, 00316 double tau_out) 00317 { 00318 double escgrd_v, 00319 tt; 00320 00321 DEBUG_ENTRY( "esc_CRDcore()" ); 00322 00323 /* find escape prob for CRD with damping wings, average of two 1-sided probs*/ 00324 00325 /* crd with wings */ 00326 if( iteration > 1 ) 00327 { 00328 /* outward optical depth if defined */ 00329 /* >>chng 03 jun 07, add test for masers here */ 00330 if( tau_out <0 || tau_in < 0. ) 00331 { 00332 /* we have a maser, use smallest optical depth to damp it out */ 00333 tt = MIN2( tau_out , tau_in ); 00334 tau_in = tt; 00335 } 00336 else 00337 { 00338 tt = tau_out - tau_in; 00339 /* help convergence by not letting tau_in go to zero at back edge of 00340 * when there was a bad guess for the total optical depth 00341 * note that this test is seldom hit since RTMakeStat does check 00342 * for overrun */ 00343 if( tt < 0. ) 00344 { 00345 tt = tau_in/SCALE; 00346 } 00347 } 00348 00349 rt.wayin = (realnum)esca0k2(tau_in); 00350 rt.wayout = (realnum)esca0k2(tt); 00351 rt.fracin = rt.wayin/(rt.wayin + rt.wayout); 00352 escgrd_v = 0.5*(rt.wayin + rt.wayout); 00353 } 00354 else 00355 { 00356 /* outward optical depth not defined, dont estimate fraction out */ 00357 rt.fracin = 0.; 00358 rt.wayout = 1.; 00359 escgrd_v = esca0k2(tau_in); 00360 rt.wayin = (realnum)escgrd_v; 00361 } 00362 00363 ASSERT( escgrd_v > 0. ); 00364 return( escgrd_v ); 00365 } 00366 00367 /*esca0k2 derive Hummer's K2 escape probability for Doppler core only */ 00368 double esca0k2(double taume) 00369 { 00370 double arg, 00371 esca0k2_v, 00372 suma, 00373 sumb, 00374 sumc, 00375 sumd, 00376 tau; 00377 static double a[5]={1.00,-0.1117897,-0.1249099917,-9.136358767e-3, 00378 -3.370280896e-4}; 00379 static double b[6]={1.00,0.1566124168,9.013261660e-3,1.908481163e-4, 00380 -1.547417750e-7,-6.657439727e-9}; 00381 static double c[5]={1.000,19.15049608,100.7986843,129.5307533,-31.43372468}; 00382 static double d[6]={1.00,19.68910391,110.2576321,169.4911399,-16.69969409, 00383 -36.664480000}; 00384 00385 DEBUG_ENTRY( "esca0k2()" ); 00386 00387 /* compute Hummer's K2 escape probability function for a=0 00388 * using approx from 00389 * >>refer line escp Hummer, D.G., xxxx, JQRST, 26, 187. 00390 * 00391 * convert to David's opacity */ 00392 tau = taume*SQRTPI; 00393 00394 if( tau < 0. ) 00395 { 00396 /* the line mased */ 00397 esca0k2_v = escmase(taume); 00398 00399 } 00400 else if( tau < 0.01 ) 00401 { 00402 esca0k2_v = 1. - 2.*tau; 00403 00404 } 00405 else if( tau <= 11. ) 00406 { 00407 suma = a[0] + tau*(a[1] + tau*(a[2] + tau*(a[3] + a[4]*tau))); 00408 sumb = b[0] + tau*(b[1] + tau*(b[2] + tau*(b[3] + tau*(b[4] + 00409 b[5]*tau)))); 00410 esca0k2_v = tau/2.5066283*log(tau/SQRTPI) + suma/sumb; 00411 00412 } 00413 else 00414 { 00415 /* large optical depth limit */ 00416 arg = 1./log(tau/SQRTPI); 00417 sumc = c[0] + arg*(c[1] + arg*(c[2] + arg*(c[3] + c[4]*arg))); 00418 sumd = d[0] + arg*(d[1] + arg*(d[2] + arg*(d[3] + arg*(d[4] + 00419 d[5]*arg)))); 00420 esca0k2_v = (sumc/sumd)/(2.*tau*sqrt(log(tau/SQRTPI))); 00421 } 00422 return( esca0k2_v ); 00423 } 00424 00425 /*escmase escape probability for negative (masing) optical depths */ 00426 STATIC void FindNeg( void ) 00427 { 00428 long int i; 00429 00430 DEBUG_ENTRY( "FindNeg()" ); 00431 00432 /* do the level 1 lines */ 00433 for( i=1; i <= nLevel1; i++ ) 00434 { 00435 /* check if a line was a strong maser */ 00436 if( TauLines[i].Emis->TauIn < -1. ) 00437 DumpLine(&TauLines[i]); 00438 } 00439 00440 /* Generic atoms & molecules from databases 00441 * added by Humeshkar Nemala*/ 00442 for( i = 0; i <linesAdded2; i++) 00443 { 00444 /* check if a line was a strong maser */ 00445 if(atmolEmis[i].TauIn < -1. ) 00446 DumpLine(atmolEmis[i].tran); 00447 } 00448 00449 /* now do the level 2 lines */ 00450 for( i=0; i < nWindLine; i++ ) 00451 { 00452 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 00453 { 00454 /* check if a line was a strong maser */ 00455 if( TauLine2[i].Emis->TauIn < -1. ) 00456 DumpLine(&TauLine2[i]); 00457 } 00458 } 00459 00460 /* now do the hyperfine structure lines */ 00461 for( i=0; i < nHFLines; i++ ) 00462 { 00463 /* check if a line was a strong maser */ 00464 if( HFLines[i].Emis->TauIn < -1. ) 00465 DumpLine(&HFLines[i]); 00466 } 00467 00468 /* now do the co carbon monoxide lines */ 00469 for( i=0; i < nCORotate; i++ ) 00470 { 00471 /* check if a line was a strong maser */ 00472 if( C12O16Rotate[i].Emis->TauIn < -1. ) 00473 DumpLine(&C12O16Rotate[i]); 00474 /* check if a line was a strong maser */ 00475 if( C13O16Rotate[i].Emis->TauIn < -1. ) 00476 DumpLine(&C13O16Rotate[i]); 00477 } 00478 return; 00479 } 00480 00481 STATIC double escmase(double tau) 00482 { 00483 double escmase_v; 00484 00485 DEBUG_ENTRY( "escmase()" ); 00486 00487 /* this is the only routine that computes maser escape probabilities */ 00488 ASSERT( tau <= 0. ); 00489 00490 if( tau > -0.1 ) 00491 { 00492 escmase_v = 1. - tau*(0.5 + tau/6.); 00493 } 00494 else if( tau > -30. ) 00495 { 00496 escmase_v = (1. - exp(-tau))/tau; 00497 } 00498 else 00499 { 00500 fprintf( ioQQQ, " DISASTER escmase called with 2big tau%10.2e\n", 00501 tau ); 00502 fprintf( ioQQQ, " This is zone number%4ld\n", nzone ); 00503 FindNeg(); 00504 ShowMe(); 00505 cdEXIT(EXIT_FAILURE); 00506 } 00507 00508 ASSERT( escmase_v >= 1. ); 00509 return( escmase_v ); 00510 } 00511 00512 /*escConE2 one of the forms of the continuum escape probability */ 00513 /*cone2 generate e2 function needed for continuum transport */ 00514 STATIC double cone2(double t); 00515 00516 double escConE2(double x) 00517 { 00518 double conesc_v; 00519 00520 DEBUG_ENTRY( "escConE2()" ); 00521 00522 conesc_v = exp(-chnukt_ContTkt*(x-1.))/ 00523 x*cone2(chnukt_ctau/x/x/x); 00524 return( conesc_v ); 00525 } 00526 00527 /*cone2 generate second exponential integral e2 function needed for continuum transport */ 00528 STATIC double cone2(double t) 00529 { 00530 double cone2_v, 00531 remain, 00532 tln; 00533 00534 DEBUG_ENTRY( "cone2()" ); 00535 00536 if( t < 80. ) 00537 { 00538 tln = exp(-t); 00539 } 00540 else 00541 { 00542 cone2_v = 0.; 00543 return( cone2_v ); 00544 } 00545 00546 /* fit of second exponential integral; 00547 * T is optical depth, and TLN is EXP(-t) 00548 * */ 00549 if( t < 0.3 ) 00550 { 00551 remain = (1.998069357 + t*(66.4037741 + t*107.2041376))/(1. + 00552 t*(37.4009646 + t*105.0388805)); 00553 00554 } 00555 else if( t < 20. ) 00556 { 00557 remain = (1.823707708 + t*2.395042899)/(1. + t*(2.488885899 - 00558 t*0.00430538)); 00559 00560 } 00561 else 00562 { 00563 remain = 1.; 00564 } 00565 00566 cone2_v = remain*tln/(2. + t); 00567 return( cone2_v ); 00568 } 00569 00570 /* a continuum escape probability */ 00571 STATIC double conrec(double x) 00572 { 00573 double conrec_v; 00574 00575 DEBUG_ENTRY( "conrec()" ); 00576 00577 conrec_v = exp(-chnukt_ContTkt*(x-1.))/x; 00578 return( conrec_v ); 00579 } 00580 00581 /*esccon continuum escape probability */ 00582 double esccon(double tau, 00583 double hnukt) 00584 { 00585 double dinc, 00586 escpcn_v, 00587 sumesc, 00588 sumrec; 00589 00590 DEBUG_ENTRY( "esccon()" ); 00591 00592 /* computes continuum escape probabilities */ 00593 if( tau < 0.01 ) 00594 { 00595 escpcn_v = 1.; 00596 return( escpcn_v ); 00597 } 00598 00599 else if( hnukt > 1. && tau > 100. ) 00600 { 00601 escpcn_v = 1e-20; 00602 return( escpcn_v ); 00603 } 00604 00605 chnukt_ContTkt = hnukt; 00606 chnukt_ctau = tau; 00607 00608 dinc = 10./hnukt; 00609 sumrec = qg32(1.,1.+dinc,conrec); 00610 sumesc = qg32(1.,1.+dinc,escConE2); 00611 00612 if( sumrec > 0. ) 00613 { 00614 escpcn_v = sumesc/sumrec; 00615 } 00616 else 00617 { 00618 escpcn_v = 0.; 00619 } 00620 return( escpcn_v ); 00621 } 00622 00623 /*RTesc_lya_1side fit Hummer and Kunasz escape probability for hydrogen atom Lya */ 00624 STATIC void RTesc_lya_1side(double taume, 00625 double beta, 00626 realnum *esc, 00627 realnum *dest, 00628 /* position of line in frequency array on c scale */ 00629 long ipLine ) 00630 { 00631 double esc0, 00632 fac, 00633 fac1, 00634 fac2, 00635 tau, 00636 taucon, 00637 taulog; 00638 00639 /* DEST0 is the smallest destruction probability to return 00640 * in high metallicity models, in rt.h 00641 const double DEST0=1e-8;*/ 00642 00643 DEBUG_ENTRY( "RTesc_lya_1side()" ); 00644 00645 /* fits to numerical results of Hummer and Kunasz Ap.J. 80 */ 00646 tau = taume*SQRTPI; 00647 00648 /* this is the real escape probability */ 00649 esc0 = 1./((0.6451 + tau)*(0.47 + 1.08/(1. + 7.3e-6*tau))); 00650 00651 esc0 = MAX2(0.,esc0); 00652 esc0 = MIN2(1.,esc0); 00653 00654 if( tau > 0. ) 00655 { 00656 taulog = log10(MIN2(1e8,tau)); 00657 } 00658 else 00659 { 00660 /* the line mased 00661 *>>chng 06 sep 08, kill xLaMase 00662 hydro.xLaMase = MIN2(hydro.xLaMase,(realnum)tau);*/ 00663 taulog = 0.; 00664 *dest = 0.; 00665 *esc = (realnum)esc0; 00666 } 00667 00668 if( beta > 0. ) 00669 { 00670 taucon = MIN2(2.,beta*tau); 00671 00672 if( taucon > 1e-3 ) 00673 { 00674 fac1 = -1.25 + 0.475*taulog; 00675 fac2 = -0.485 + 0.1615*taulog; 00676 fac = -fac1*taucon + fac2*POW2(taucon); 00677 fac = pow(10.,fac); 00678 fac = MIN2(1.,fac); 00679 } 00680 else 00681 { 00682 fac = 1.; 00683 } 00684 00685 *esc = (realnum)(esc0*fac); 00686 /* MIN puts cat at 50 */ 00687 *dest = (realnum)(beta/(0.30972 - MIN2(.28972,0.03541667*taulog))); 00688 } 00689 00690 else 00691 { 00692 *dest = 0.; 00693 *esc = (realnum)esc0; 00694 } 00695 00696 *dest = MIN2(*dest,1.f-*esc); 00697 *dest = MAX2(0.f,*dest); 00698 00699 /* >>chng 99 apr 12, limit destruction prob in case where gas dominated by scattering. 00700 * in this case scattering is much more likely than absorption on this event */ 00701 *dest = (realnum)( (1. - opac.albedo[ipLine]) * *dest + opac.albedo[ipLine]*DEST0); 00702 /* this is for debugging H Lya */ 00703 { 00704 /*@-redef@*/ 00705 enum {BUG=false}; 00706 /*@+redef@*/ 00707 if( BUG ) 00708 { 00709 fprintf(ioQQQ,"scatdest tau %.2e beta%.2e 1-al%.2e al%.2e dest%.2e \n", 00710 taume, 00711 beta, 00712 (1. - opac.albedo[ipLine]), 00713 opac.albedo[ipLine] , 00714 *dest 00715 ); 00716 } 00717 } 00718 return; 00719 } 00720 00721 /*RT_DestProb returns line destruction probability due to continuum opacity */ 00722 double RT_DestProb( 00723 /* abundance of species */ 00724 double abund, 00725 /* its line absorption cross section */ 00726 double crsec, 00727 /* pointer to energy within continuum array, to get background opacity, 00728 * this is on the f not c scale */ 00729 long int ipanu, 00730 /* line width */ 00731 double widl, 00732 /* escape probability */ 00733 double escp, 00734 /* type of redistribution function */ 00735 int nCore) 00736 { 00737 /* this will be the value we shall return */ 00738 double eovrlp_v; 00739 00740 double conopc, 00741 beta; 00742 00743 /* DEST0 is the smallest destruction probability to return 00744 * in high metallicity models 00745 * this was set to 1e-8 until 99nov18, 00746 * in cooling flow model the actual Lya ots dest prob was 1e-16, 00747 * and this lower limit of 1e-8 caused energy balance problems, 00748 * since dest prob was 8 orders of magnitude too great. 00749 * >>chng 99 nov 18, to 1e-20, but beware that comments indicate that 00750 * this will cause problems with high metallicity clouds(?) */ 00751 /* >>chng 00 jun 04, to 0 since very feeble ionization clouds, with almost zero opacity, 00752 * this was a LARGE number */ 00753 /*const double DEST0=1e-20; 00754 const double DEST0=0.;*/ 00755 00756 DEBUG_ENTRY( "RT_DestProb()" ); 00757 00758 /* computes "escape probability" due to continuum destruction of 00759 * 00760 * if esc prob gt 1 then line is masing - return small number for dest prob */ 00761 /* >>>chng 99 apr 10, return min dest when scattering greater than abs */ 00762 /* no idea of opacity whatsoever, on very first soln for this model */ 00763 /* >>chng 05 mar 20, add test on line being above upper bound of frequency 00764 * do not want to evaluate opacity in this case since off scale */ 00765 if( escp >= 1.0 || !conv.nTotalIoniz || ipanu >= rfield.nflux ) 00766 { 00767 eovrlp_v = 0.; 00768 return( eovrlp_v ); 00769 } 00770 00771 /* find continuum opacity */ 00772 conopc = opac.opacity_abs[ipanu-1]; 00773 00774 ASSERT( crsec > 0. ); 00775 00776 /* may be no population, cannot use above test since return 0 not DEST0 */ 00777 if( abund <= 0. || conopc <= 0. ) 00778 { 00779 /* do not set this to DEST0 since energy not then conserved */ 00780 eovrlp_v = 0.; 00781 return( eovrlp_v ); 00782 } 00783 00784 /* fac of 1.7 convert to Hummer convention for line opacity */ 00785 beta = conopc/(abund*SQRTPI*crsec/widl + conopc); 00786 /* >>chng 04 may 10, rm * 1-pesc) 00787 beta = MIN2(beta,(1.-escp)); */ 00788 00789 if( nCore == ipDEST_INCOM ) 00790 { 00791 /* fits to 00792 * >>>refer la esc Hummer and Kunasz 1980 Ap.J. 236,609. 00793 * the max value of 1e-3 is so that we do not go too far 00794 * beyond what Hummer and Kunasz did, discussed in 00795 * >>refer rt esc proc Ferland, G.J., 1999, ApJ, 512, 247 */ 00798 eovrlp_v = MIN2(1e-3,8.5*beta); 00799 } 00800 else if( nCore == ipDEST_K2 ) 00801 { 00802 /* Doppler core only; a=0., Hummer 68 00803 eovrlp_v = RT_DestHummer(beta);*/ 00804 eovrlp_v = MIN2(1e-3,8.5*beta); 00805 } 00806 else if( nCore == ipDEST_SIMPL ) 00807 { 00808 /* this for debugging only 00809 eovrlp_v = 8.5*beta;*/ 00810 /* >>chng 04 may 13, use same min function */ 00811 eovrlp_v = MIN2(1e-3,8.5*beta); 00812 } 00813 else 00814 { 00815 fprintf( ioQQQ, " chCore of %i not understood by RT_DestProb.\n", 00816 nCore ); 00817 cdEXIT(EXIT_FAILURE); 00818 } 00819 00820 /* renorm to unity */ 00821 eovrlp_v /= 1. + eovrlp_v; 00822 00823 /* multiply by 1-escape prob, since no destruction when optically thin */ 00824 eovrlp_v *= 1. - escp; 00825 00826 /*check results in bounds */ 00827 ASSERT( eovrlp_v >= 0. ); 00828 ASSERT( eovrlp_v <= 1. ); 00829 00830 { 00831 /* debugging code for Lya problems */ 00832 /*@-redef@*/ 00833 enum {DEBUG_LOC=false}; 00834 /*@+redef@*/ 00835 if( DEBUG_LOC ) 00836 { 00837 if( rfield.anu[ipanu-1]>0.73 && rfield.anu[ipanu-1]<0.76 && 00838 fp_equal( abund, Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->PopOpc ) ) 00839 { 00840 fprintf(ioQQQ,"%li RT_DestProb\t%g\n", 00841 nzone, eovrlp_v ); 00842 } 00843 } 00844 } 00845 00846 /* >>chng 04 may 10, rm min */ 00847 /* this hack removed since no fundamental reason for it to be here, 00848 * this should be added to scattering escape, if included at all */ 00849 # if 0 00850 /* >>chng 99 apr 12, limit destruction prob in case where gas dominated by scattering. 00851 * in this case scattering is much more likely than absorption on this event 00852 eovrlp_v = (1. - opac.albedo[ipanu-1]) * eovrlp_v + 00853 opac.albedo[ipanu-1]*DEST0; */ 00854 /* >>chng 01 aug 11, add factor of 3 for increase in mean free path, and min on 0 */ 00855 /*eovrlp_v = MAX2(DEST0,1. - 3.*opac.albedo[ipanu-1]) * eovrlp_v + 00856 opac.albedo[ipanu-1]*DEST0;*/ 00857 eovrlp_v = POW2(1. - opac.albedo[ipanu-1]) * eovrlp_v + 00858 opac.albedo[ipanu-1]*0.; 00859 # endif 00860 00861 return( eovrlp_v ); 00862 } 00863 00864 /*RT_LineWidth compute line width (cm/sec), using optical depth array information 00865 * this is where effects of wind are done */ 00866 double RT_LineWidth(const transition * t) 00867 { 00868 double RT_LineWidth_v, 00869 aa, 00870 atau, 00871 b, 00872 r, 00873 vth; 00874 realnum tau; 00875 00876 DEBUG_ENTRY( "RT_LineWidth()" ); 00877 00878 /* uses line width from 00879 * >>refer esc prob Bonilha et al. Ap.J. (1979) 233 649 00880 * return value is half velocity width*(1-ESC PROB) [cm s-1] 00881 * this assumes incomplete redistribution, damp.tau^1/3 width */ 00882 00883 /* thermal width */ 00884 vth = DoppVel.doppler[ t->Hi->nelem-1 ]; 00885 00886 /* optical depth in outer direction is defined 00887 * on second and later iterations. 00888 * smaller of inner and outer optical depths is chosen for esc prob */ 00889 if( iteration > 1 ) 00890 { 00891 /* optical depth to shielded face */ 00892 realnum tauout = t->Emis->TauTot - t->Emis->TauIn; 00893 00894 /* >>chng 99 apr 22 use smaller optical depth */ 00895 tau = MIN2( t->Emis->TauIn , tauout ); 00896 } 00897 else 00898 { 00899 tau = t->Emis->TauIn; 00900 } 00901 /* do not evaluate line width if quite optically thin - will be dominated 00902 * by noise in this case */ 00903 if( tau <1e-3 ) 00904 return 0; 00905 00906 double Pesc = esc_PRD_1side( tau , t->Emis->damp); 00907 00908 /* max optical depth is thermalization length */ 00909 realnum therm = (realnum)(5.3e16/dense.eden); 00910 if( tau > therm ) 00911 { 00912 pressure.lgPradDen = true; 00913 tau = therm; 00914 } 00915 00916 /* >>chng 01 jun 23, use wind vel instead of rt since rt deleted */ 00917 /* >>chng 04 may 13, use thermal for subsonic cases */ 00922 if( wind.windv <= 0. ) 00923 { 00924 /* static geometry */ 00925 /* esc prob has noise if smaller than FLT_EPSILON, or is masing */ 00926 if( (tau-opac.taumin)/100. < FLT_EPSILON ) 00927 { 00928 RT_LineWidth_v = 0.; 00929 } 00930 else if( tau <= 20. ) 00931 { 00932 atau = -6.907755; 00933 if( tau > 1e-3 ) 00934 atau = log(tau); 00935 aa = 4.8 + 5.2*tau + (4.*tau - 1.)*atau; 00936 b = 6.5*tau - atau; 00937 RT_LineWidth_v = vth*0.8862*aa/b*(1. - MIN2( 1., 00938 ( Pesc + t->Emis->Pelec_esc + t->Emis->Pdest) ) ); 00939 /* small number roundoff can dominate this process */ 00940 if( Pesc < 10. * FLT_EPSILON ) 00941 RT_LineWidth_v = 0.; 00942 } 00943 else 00944 { 00945 ASSERT( t->Emis->damp*tau >= 0.); 00946 atau = log(MAX2(0.0001,tau)); 00947 aa = 1. + 2.*atau/pow(1. + 0.3*t->Emis->damp*tau,0.6667) + pow(6.5* 00948 t->Emis->damp*tau,0.333); 00949 b = 1.6 + 1.5/(1. + 0.20*t->Emis->damp*tau); 00950 RT_LineWidth_v = vth*0.8862*aa/b*(1. - MIN2( 1. , 00951 (Pesc+ t->Emis->Pelec_esc + t->Emis->Pdest)) ); 00952 } 00953 00954 /* we want full width, not half width */ 00955 RT_LineWidth_v *= 2.; 00956 00957 } 00958 else 00959 { 00960 /* wind */ 00961 r = t->Emis->damp*tau/PI; 00962 if( r <= 1. ) 00963 { 00964 RT_LineWidth_v = vth*sqrt(log(MAX2(1.,tau))*.2821); 00965 } 00966 else 00967 { 00968 RT_LineWidth_v = 2.*fabs(wind.windv0); 00969 if( r*vth <= RT_LineWidth_v ) 00970 { 00971 RT_LineWidth_v = vth*r*log(RT_LineWidth_v/(r*vth)); 00972 } 00973 } 00974 } 00975 00976 ASSERT( RT_LineWidth_v >= 0. ); 00977 return( RT_LineWidth_v ); 00978 } 00979 00980 /*RT_DestHummer evaluate Hummer's betaF(beta) function */ 00981 double RT_DestHummer(double beta) /* beta is ratio of continuum to mean line opacity, 00982 * returns dest prob = beta F(beta) */ 00983 { 00984 double fhummr_v, 00985 x; 00986 00987 DEBUG_ENTRY( "RT_DestHummer()" ); 00988 00989 /* evaluates Hummer's F(beta) function for case where damping 00990 * constant is zero, are returns beta*F(beta) 00991 * fit to Table 1, page 80, of Hummer MNRAS 138, 73-108. 00992 * beta is ratio of continuum to line opacity; FUMMER is 00993 * product of his F() times beta; the total destruction prob 00994 * this beta is Hummer's normalization of the Voigt function */ 00995 00996 ASSERT( beta >= 0.);/* non-positive is unphysical */ 00997 if( beta <= 0. ) 00998 { 00999 fhummr_v = 0.; 01000 } 01001 else 01002 { 01003 x = log10(beta); 01004 if( x < -5.5 ) 01005 { 01006 fhummr_v = 3.8363 - 0.56329*x; 01007 } 01008 else if( x < -3.5 ) 01009 { 01010 fhummr_v = 2.79153 - 0.75325*x; 01011 } 01012 else if( x < -2. ) 01013 { 01014 fhummr_v = 1.8446 - 1.0238*x; 01015 } 01016 else 01017 { 01018 fhummr_v = 0.72500 - 1.5836*x; 01019 } 01020 fhummr_v *= beta; 01021 } 01022 return( fhummr_v ); 01023 }