cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
rt_escprob.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*esc_CRDwing_1side fundamental escape probability radiative transfer routine, for complete redistribution */
4 /*esc_PRD_1side fundamental escape probability radiative transfer routine for incomplete redistribution */
5 /*RTesc_lya escape prob for hydrogen atom Lya, using Hummer and Kunasz results,
6  * called by hydropesc */
7 /*esc_PRD escape probability radiative transfer for incomplete redistribution */
8 /*esc_CRDwing escape probability for CRD with wings */
9 /*esc_CRDcore escape probability for CRD with no wings */
10 /*esca0k2 derive Hummer's K2 escape probability for Doppler core only */
11 /*RT_DestProb returns line destruction probability due to continuum opacity */
12 /*RT_LineWidth determine half width of any line with known optical depths */
13 #include "cddefines.h"
14 #include "physconst.h"
15 #define SCALE 2.
16 #include "dense.h"
17 #include "conv.h"
18 #include "rfield.h"
19 #include "opacity.h"
20 #include "lines_service.h"
21 #include "taulines.h"
22 #include "doppvel.h"
23 #include "pressure.h"
24 #include "wind.h"
25 #include "rt.h"
26 /*
27  *variables used to pass continuum optical depth and temp to
28  *routine that integrates over continuum
29  */
30 static double chnukt_ContTkt, chnukt_ctau;
31 
32 /*escmase escape probability for negative (masing) optical depths,*/
33 STATIC double escmase(double tau);
34 /*RTesc_lya_1side fit Hummer and Kunasz escape probability for hydrogen atom Lya */
35 STATIC void RTesc_lya_1side(double taume,
36  double beta,
37  realnum *esc,
38  realnum *dest,
39  /* position of line in frequency array on c scale */
40  long ipLine );
41 
42 double esc_PRD_1side(double tau,
43  double a)
44 {
45  double atau,
46  b,
47  escinc_v;
48 
49  DEBUG_ENTRY( "esc_PRD_1side()" );
50  ASSERT( a>0.0 );
51 
52  /* this is one of the three fundamental escape probability routines
53  * the three are esc_CRDwing_1side, esc_PRD_1side, and RTesc_lya
54  * it computes esc prob for incomplete redistribution
55  * */
56 # if 0
57  if( strcmp(rt.chEscFunSubord,"SIMP") == 0 )
58  {
59  /* this set with "escape simple" command, used for debugging */
60  escinc_v = 1./(1. + tau);
61  return( escinc_v );
62  }
63 # endif
64 
65  if( tau < 0. )
66  {
67  /* line mased */
68  escinc_v = escmase(tau);
69  }
70  else if( tau < 10. )
71  {
72  /* linear part of doppler core */
73  escinc_v = 1./(1. + 1.6*tau);
74  }
75  else
76  {
77  /* first find coefficient b(tau) */
78  atau = a*tau;
79  if( atau > 1. )
80  {
81  b = 1.6 + (3.*pow(2.*a,-0.12))/(1. + atau);
82  }
83  else
84  {
85  b = 1.6 + (3.*pow(2.*a,-0.12))/(1. + 1./sqrt(atau));
86  }
87  b = MIN2(6.,b);
88 
89  escinc_v = 1./(1. + b*tau);
90  }
91  return( escinc_v );
92 }
93 
94 /*esc_CRDwing_1side fundamental escape probability radiative transfer routine, for complete redistribution */
95 double esc_CRDwing_1side(double tau,
96  double a )
97 {
98  double esccom_v;
99 
100  DEBUG_ENTRY( "esc_CRDwing_1side()" );
101 
102  /* this is one of the three fundamental escape probability routines
103  * the three are esc_CRDwing_1side, esc_PRD_1side, and RTesc_lya
104  * it computes esc prob for complete redistribution with wings
105  * computes escape prob for complete redistribution in one direction
106  * */
107 
108  /* this is the only case that this routine computes,
109  * and is the usual case for subordinate lines,
110  * complete redistribution with damping wings */
111  esccom_v = esca0k2(tau);
112  if( tau > 1e3 )
113  {
114  esccom_v += 0.333*sqrt(a/(SQRTPI*tau));
115  }
116  return( esccom_v );
117 }
118 
119 /*RTesc_lya escape prob for hydrogen atom Lya, using
120  >>refer La escp Hummer, D.G., & Kunasz, P.B., 1980, ApJ, 236, 609
121  * called by hydropesc, return value is escape probability */
122 double RTesc_lya(
123  /* the inward escape probability */
124  double *esin,
125  /* the destruction probility */
126  double *dest,
127  /* abundance of the species */
128  double abund,
129  /* element number, 0 for H */
130  long int nelem)
131 {
132  double beta,
133  conopc,
134  escla_v;
135  realnum dstin,
136  dstout;
137 
138  DEBUG_ENTRY( "RTesc_lya()" );
139 
140  /*
141  * this is one of the three fundamental escape probability functions
142  * the three are esc_CRDwing_1side, esc_PRD_1side, and RTesc_lya
143  * evaluate esc prob for LA
144  * optical depth in outer direction always defined
145  */
146 
147  /* check charge */
148  ASSERT( nelem >= 0 && nelem < LIMELM );
149 
150  if( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot -
151  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn < 0. )
152  {
153  /* this is the case if we overrun the optical depth scale
154  * just leave things as they are */
155  escla_v = Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pesc;
156  rt.fracin = Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->FracInwd;
157  *esin = rt.fracin;
158  *dest = Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pdest;
159  return( escla_v );
160  }
161 
162  /* incomplete redistribution */
163  conopc = opac.opacity_abs[Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].ipCont-1];
164  if( abund > 0. )
165  {
166  /* the continuous opacity is positive, we have a valid soln */
167  beta = conopc/(abund/SQRTPI*Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity/
168  DoppVel.doppler[nelem] + conopc);
169  }
170  else
171  {
172  /* abundance is zero, set miniumum dest prob */
173  beta = 1e-10;
174  }
175 
176  /* find rt.wayin, the escape prob in inward direction */
178  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn,
179  beta,
180  &rt.wayin,
181  &dstin ,
182  /* position of line in energy array on C scale */
183  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].ipCont-1);
184 
185  ASSERT( (rt.wayin <= 1.) && (rt.wayin >= 0.) && (dstin <= 1.) && (dstin >= 0.) );
186 
187  /* find rt.wayout, the escape prob in outward direction */
188  RTesc_lya_1side(MAX2(Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot/100.,
189  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot-Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn),
190  beta,
191  &rt.wayout,
192  &dstout,
193  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].ipCont-1);
194 
195  ASSERT( (rt.wayout <= 1.) && (rt.wayout >= 0.) && (dstout <= 1.) && (dstout >= 0.) );
196 
197  /* esc prob is mean of in and out */
198  escla_v = (rt.wayin + rt.wayout)/2.;
199  /* the inward escaping part of the line */
200  *esin = rt.wayin;
201 
202  /* dest prob is mean of in and out */
203  *dest = (dstin + dstout)/2.f;
204  /* >>chng 02 oct 02, sum of escape and dest prob must be less then unity,
205  * for very thin models this forces dest prob to go to zero,
206  * rather than the value of DEST0, caught by Jon Slavin */
207  *dest = (realnum)MIN2( *dest , 1.-escla_v );
208  /* but dest prob can't be negative */
209  *dest = (realnum)MAX2(0., *dest );
210 
211  /* fraction of line emitted in inward direction */
212  rt.fracin = rt.wayin/(rt.wayin + rt.wayout);
213  ASSERT( escla_v >=0. && *dest>=0. && *esin>=0. );
214  return( escla_v );
215 }
216 
217 /*esc_PRD escape probability radiative transfer for incomplete redistribution */
218 double esc_PRD(double tau,
219  double tau_out,
220  double damp )
221 {
222  double escgrd_v,
223  tt;
224 
225  DEBUG_ENTRY( "esc_PRD()" );
226 
227  /* find escape prob for incomp redis, average of two 1-sided probs*/
228 
229  if( iteration > 1 )
230  {
231  /* outward optical depth if defined */
232  tt = tau_out - tau;
233  /* help convergence by not letting tau go to zero at back edge of
234  * when there was a bad guess for the total optical depth
235  * note that this test is seldom hit since RTMakeStat does check
236  * for overrun */
237  if( tt < 0. )
238  {
239  tt = tau/SCALE;
240  }
241 
242  rt.wayin = (realnum)esc_PRD_1side(tau,damp);
243  rt.wayout = (realnum)esc_PRD_1side(tt,damp);
244  rt.fracin = rt.wayin/(rt.wayin + rt.wayout);
245  escgrd_v = 0.5*(rt.wayin + rt.wayout);
246  }
247  else
248  {
249  /* outward optical depth not defined, dont estimate fraction out */
250  rt.fracin = 0.;
251  rt.wayout = 1.;
252  escgrd_v = esc_PRD_1side(tau,damp);
253  rt.wayin = (realnum)escgrd_v;
254  }
255 
256  ASSERT( escgrd_v > 0. );
257  return( escgrd_v );
258 }
259 
260 /*esc_CRDwing escape probability radiative transfer for CRDS in core only */
261 double esc_CRDwing(double tau_in,
262  double tau_out,
263  double damp)
264 {
265  double escgrd_v,
266  tt;
267 
268  DEBUG_ENTRY( "esc_CRDwing()" );
269 
270  /* find escape prob for CRD with damping wings, average of two 1-sided probs*/
271 
272  /* crd with wings */
273  if( iteration > 1 )
274  {
275  /* outward optical depth if defined */
276  /* >>chng 03 jun 07, add test for masers here */
277  if( tau_out <0 || tau_in < 0. )
278  {
279  /* we have a maser, use smallest optical depth to damp it out */
280  tt = MIN2( tau_out , tau_in );
281  tau_in = tt;
282  }
283  else
284  {
285  tt = tau_out - tau_in;
286  /* help convergence by not letting tau_in go to zero at back edge of
287  * when there was a bad guess for the total optical depth
288  * note that this test is seldom hit since RTMakeStat does check
289  * for overrun */
290  if( tt < 0. )
291  {
292  tt = tau_in/SCALE;
293  }
294  }
295 
296  rt.wayin = (realnum)esc_CRDwing_1side(tau_in,damp);
297  rt.wayout = (realnum)esc_CRDwing_1side(tt,damp);
298  rt.fracin = rt.wayin/(rt.wayin + rt.wayout);
299  escgrd_v = 0.5*(rt.wayin + rt.wayout);
300  }
301  else
302  {
303  /* outward optical depth not defined, dont estimate fraction out */
304  rt.fracin = 0.;
305  rt.wayout = 1.;
306  escgrd_v = esc_CRDwing_1side(tau_in,damp);
307  rt.wayin = (realnum)escgrd_v;
308  }
309 
310  ASSERT( escgrd_v > 0. );
311  return( escgrd_v );
312 }
313 
314 /*esc_CRDwing escape probability radiative transfer for incomplete redistribution */
315 double esc_CRDcore(double tau_in,
316  double tau_out)
317 {
318  double escgrd_v,
319  tt;
320 
321  DEBUG_ENTRY( "esc_CRDcore()" );
322 
323  /* find escape prob for CRD with damping wings, average of two 1-sided probs*/
324 
325  /* crd with wings */
326  if( iteration > 1 )
327  {
328  /* outward optical depth if defined */
329  /* >>chng 03 jun 07, add test for masers here */
330  if( tau_out <0 || tau_in < 0. )
331  {
332  /* we have a maser, use smallest optical depth to damp it out */
333  tt = MIN2( tau_out , tau_in );
334  tau_in = tt;
335  }
336  else
337  {
338  tt = tau_out - tau_in;
339  /* help convergence by not letting tau_in go to zero at back edge of
340  * when there was a bad guess for the total optical depth
341  * note that this test is seldom hit since RTMakeStat does check
342  * for overrun */
343  if( tt < 0. )
344  {
345  tt = tau_in/SCALE;
346  }
347  }
348 
349  rt.wayin = (realnum)esca0k2(tau_in);
350  rt.wayout = (realnum)esca0k2(tt);
351  rt.fracin = rt.wayin/(rt.wayin + rt.wayout);
352  escgrd_v = 0.5*(rt.wayin + rt.wayout);
353  }
354  else
355  {
356  /* outward optical depth not defined, dont estimate fraction out */
357  rt.fracin = 0.;
358  rt.wayout = 1.;
359  escgrd_v = esca0k2(tau_in);
360  rt.wayin = (realnum)escgrd_v;
361  }
362 
363  ASSERT( escgrd_v > 0. );
364  return( escgrd_v );
365 }
366 
367 /*esca0k2 derive Hummer's K2 escape probability for Doppler core only */
368 double esca0k2(double taume)
369 {
370  double arg,
371  esca0k2_v,
372  suma,
373  sumb,
374  sumc,
375  sumd,
376  tau;
377  static double a[5]={1.00,-0.1117897,-0.1249099917,-9.136358767e-3,
378  -3.370280896e-4};
379  static double b[6]={1.00,0.1566124168,9.013261660e-3,1.908481163e-4,
380  -1.547417750e-7,-6.657439727e-9};
381  static double c[5]={1.000,19.15049608,100.7986843,129.5307533,-31.43372468};
382  static double d[6]={1.00,19.68910391,110.2576321,169.4911399,-16.69969409,
383  -36.664480000};
384 
385  DEBUG_ENTRY( "esca0k2()" );
386 
387  /* compute Hummer's K2 escape probability function for a=0
388  * using approx from
389  * >>refer line escp Hummer, D.G., xxxx, JQRST, 26, 187.
390  *
391  * convert to David's opacity */
392  tau = taume*SQRTPI;
393 
394  if( tau < 0. )
395  {
396  /* the line mased */
397  esca0k2_v = escmase(taume);
398 
399  }
400  else if( tau < 0.01 )
401  {
402  esca0k2_v = 1. - 2.*tau;
403 
404  }
405  else if( tau <= 11. )
406  {
407  suma = a[0] + tau*(a[1] + tau*(a[2] + tau*(a[3] + a[4]*tau)));
408  sumb = b[0] + tau*(b[1] + tau*(b[2] + tau*(b[3] + tau*(b[4] +
409  b[5]*tau))));
410  esca0k2_v = tau/2.5066283*log(tau/SQRTPI) + suma/sumb;
411 
412  }
413  else
414  {
415  /* large optical depth limit */
416  arg = 1./log(tau/SQRTPI);
417  sumc = c[0] + arg*(c[1] + arg*(c[2] + arg*(c[3] + c[4]*arg)));
418  sumd = d[0] + arg*(d[1] + arg*(d[2] + arg*(d[3] + arg*(d[4] +
419  d[5]*arg))));
420  esca0k2_v = (sumc/sumd)/(2.*tau*sqrt(log(tau/SQRTPI)));
421  }
422  return( esca0k2_v );
423 }
424 
425 /*escmase escape probability for negative (masing) optical depths */
426 STATIC void FindNeg( void )
427 {
428  long int i;
429 
430  DEBUG_ENTRY( "FindNeg()" );
431 
432  /* do the level 1 lines */
433  for( i=1; i <= nLevel1; i++ )
434  {
435  /* check if a line was a strong maser */
436  if( TauLines[i].Emis->TauIn < -1. )
437  DumpLine(&TauLines[i]);
438  }
439 
440  /* Generic atoms & molecules from databases
441  * added by Humeshkar Nemala*/
442  for( i = 0; i <linesAdded2; i++)
443  {
444  /* check if a line was a strong maser */
445  if(atmolEmis[i].TauIn < -1. )
446  DumpLine(atmolEmis[i].tran);
447  }
448 
449  /* now do the level 2 lines */
450  for( i=0; i < nWindLine; i++ )
451  {
452  if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
453  {
454  /* check if a line was a strong maser */
455  if( TauLine2[i].Emis->TauIn < -1. )
456  DumpLine(&TauLine2[i]);
457  }
458  }
459 
460  /* now do the hyperfine structure lines */
461  for( i=0; i < nHFLines; i++ )
462  {
463  /* check if a line was a strong maser */
464  if( HFLines[i].Emis->TauIn < -1. )
465  DumpLine(&HFLines[i]);
466  }
467 
468  /* now do the co carbon monoxide lines */
469  for( i=0; i < nCORotate; i++ )
470  {
471  /* check if a line was a strong maser */
472  if( C12O16Rotate[i].Emis->TauIn < -1. )
473  DumpLine(&C12O16Rotate[i]);
474  /* check if a line was a strong maser */
475  if( C13O16Rotate[i].Emis->TauIn < -1. )
476  DumpLine(&C13O16Rotate[i]);
477  }
478  return;
479 }
480 
481 STATIC double escmase(double tau)
482 {
483  double escmase_v;
484 
485  DEBUG_ENTRY( "escmase()" );
486 
487  /* this is the only routine that computes maser escape probabilities */
488  ASSERT( tau <= 0. );
489 
490  if( tau > -0.1 )
491  {
492  escmase_v = 1. - tau*(0.5 + tau/6.);
493  }
494  else if( tau > -30. )
495  {
496  escmase_v = (1. - exp(-tau))/tau;
497  }
498  else
499  {
500  fprintf( ioQQQ, " DISASTER escmase called with 2big tau%10.2e\n",
501  tau );
502  fprintf( ioQQQ, " This is zone number%4ld\n", nzone );
503  FindNeg();
504  ShowMe();
505  cdEXIT(EXIT_FAILURE);
506  }
507 
508  ASSERT( escmase_v >= 1. );
509  return( escmase_v );
510 }
511 
512 /*escConE2 one of the forms of the continuum escape probability */
513 /*cone2 generate e2 function needed for continuum transport */
514 STATIC double cone2(double t);
515 
516 double escConE2(double x)
517 {
518  double conesc_v;
519 
520  DEBUG_ENTRY( "escConE2()" );
521 
522  conesc_v = exp(-chnukt_ContTkt*(x-1.))/
523  x*cone2(chnukt_ctau/x/x/x);
524  return( conesc_v );
525 }
526 
527 /*cone2 generate second exponential integral e2 function needed for continuum transport */
528 STATIC double cone2(double t)
529 {
530  double cone2_v,
531  remain,
532  tln;
533 
534  DEBUG_ENTRY( "cone2()" );
535 
536  if( t < 80. )
537  {
538  tln = exp(-t);
539  }
540  else
541  {
542  cone2_v = 0.;
543  return( cone2_v );
544  }
545 
546  /* fit of second exponential integral;
547  * T is optical depth, and TLN is EXP(-t)
548  * */
549  if( t < 0.3 )
550  {
551  remain = (1.998069357 + t*(66.4037741 + t*107.2041376))/(1. +
552  t*(37.4009646 + t*105.0388805));
553 
554  }
555  else if( t < 20. )
556  {
557  remain = (1.823707708 + t*2.395042899)/(1. + t*(2.488885899 -
558  t*0.00430538));
559 
560  }
561  else
562  {
563  remain = 1.;
564  }
565 
566  cone2_v = remain*tln/(2. + t);
567  return( cone2_v );
568 }
569 
570 /* a continuum escape probability */
571 STATIC double conrec(double x)
572 {
573  double conrec_v;
574 
575  DEBUG_ENTRY( "conrec()" );
576 
577  conrec_v = exp(-chnukt_ContTkt*(x-1.))/x;
578  return( conrec_v );
579 }
580 
581 /*esccon continuum escape probability */
582 double esccon(double tau,
583  double hnukt)
584 {
585  double dinc,
586  escpcn_v,
587  sumesc,
588  sumrec;
589 
590  DEBUG_ENTRY( "esccon()" );
591 
592  /* computes continuum escape probabilities */
593  if( tau < 0.01 )
594  {
595  escpcn_v = 1.;
596  return( escpcn_v );
597  }
598 
599  else if( hnukt > 1. && tau > 100. )
600  {
601  escpcn_v = 1e-20;
602  return( escpcn_v );
603  }
604 
605  chnukt_ContTkt = hnukt;
606  chnukt_ctau = tau;
607 
608  dinc = 10./hnukt;
609  sumrec = qg32(1.,1.+dinc,conrec);
610  sumesc = qg32(1.,1.+dinc,escConE2);
611 
612  if( sumrec > 0. )
613  {
614  escpcn_v = sumesc/sumrec;
615  }
616  else
617  {
618  escpcn_v = 0.;
619  }
620  return( escpcn_v );
621 }
622 
623 /*RTesc_lya_1side fit Hummer and Kunasz escape probability for hydrogen atom Lya */
624 STATIC void RTesc_lya_1side(double taume,
625  double beta,
626  realnum *esc,
627  realnum *dest,
628  /* position of line in frequency array on c scale */
629  long ipLine )
630 {
631  double esc0,
632  fac,
633  fac1,
634  fac2,
635  tau,
636  taucon,
637  taulog;
638 
639  /* DEST0 is the smallest destruction probability to return
640  * in high metallicity models, in rt.h
641  const double DEST0=1e-8;*/
642 
643  DEBUG_ENTRY( "RTesc_lya_1side()" );
644 
645  /* fits to numerical results of Hummer and Kunasz Ap.J. 80 */
646  tau = taume*SQRTPI;
647 
648  /* this is the real escape probability */
649  esc0 = 1./((0.6451 + tau)*(0.47 + 1.08/(1. + 7.3e-6*tau)));
650 
651  esc0 = MAX2(0.,esc0);
652  esc0 = MIN2(1.,esc0);
653 
654  if( tau > 0. )
655  {
656  taulog = log10(MIN2(1e8,tau));
657  }
658  else
659  {
660  /* the line mased
661  *>>chng 06 sep 08, kill xLaMase
662  hydro.xLaMase = MIN2(hydro.xLaMase,(realnum)tau);*/
663  taulog = 0.;
664  *dest = 0.;
665  *esc = (realnum)esc0;
666  }
667 
668  if( beta > 0. )
669  {
670  taucon = MIN2(2.,beta*tau);
671 
672  if( taucon > 1e-3 )
673  {
674  fac1 = -1.25 + 0.475*taulog;
675  fac2 = -0.485 + 0.1615*taulog;
676  fac = -fac1*taucon + fac2*POW2(taucon);
677  fac = pow(10.,fac);
678  fac = MIN2(1.,fac);
679  }
680  else
681  {
682  fac = 1.;
683  }
684 
685  *esc = (realnum)(esc0*fac);
686  /* MIN puts cat at 50 */
687  *dest = (realnum)(beta/(0.30972 - MIN2(.28972,0.03541667*taulog)));
688  }
689 
690  else
691  {
692  *dest = 0.;
693  *esc = (realnum)esc0;
694  }
695 
696  *dest = MIN2(*dest,1.f-*esc);
697  *dest = MAX2(0.f,*dest);
698 
699  /* >>chng 99 apr 12, limit destruction prob in case where gas dominated by scattering.
700  * in this case scattering is much more likely than absorption on this event */
701  *dest = (realnum)( (1. - opac.albedo[ipLine]) * *dest + opac.albedo[ipLine]*DEST0);
702  /* this is for debugging H Lya */
703  {
704  /*@-redef@*/
705  enum {BUG=false};
706  /*@+redef@*/
707  if( BUG )
708  {
709  fprintf(ioQQQ,"scatdest tau %.2e beta%.2e 1-al%.2e al%.2e dest%.2e \n",
710  taume,
711  beta,
712  (1. - opac.albedo[ipLine]),
713  opac.albedo[ipLine] ,
714  *dest
715  );
716  }
717  }
718  return;
719 }
720 
721 /*RT_DestProb returns line destruction probability due to continuum opacity */
722 double RT_DestProb(
723  /* abundance of species */
724  double abund,
725  /* its line absorption cross section */
726  double crsec,
727  /* pointer to energy within continuum array, to get background opacity,
728  * this is on the f not c scale */
729  long int ipanu,
730  /* line width */
731  double widl,
732  /* escape probability */
733  double escp,
734  /* type of redistribution function */
735  int nCore)
736 {
737  /* this will be the value we shall return */
738  double eovrlp_v;
739 
740  double conopc,
741  beta;
742 
743  /* DEST0 is the smallest destruction probability to return
744  * in high metallicity models
745  * this was set to 1e-8 until 99nov18,
746  * in cooling flow model the actual Lya ots dest prob was 1e-16,
747  * and this lower limit of 1e-8 caused energy balance problems,
748  * since dest prob was 8 orders of magnitude too great.
749  * >>chng 99 nov 18, to 1e-20, but beware that comments indicate that
750  * this will cause problems with high metallicity clouds(?) */
751  /* >>chng 00 jun 04, to 0 since very feeble ionization clouds, with almost zero opacity,
752  * this was a LARGE number */
753  /*const double DEST0=1e-20;
754  const double DEST0=0.;*/
755 
756  DEBUG_ENTRY( "RT_DestProb()" );
757 
758  /* computes "escape probability" due to continuum destruction of
759  *
760  * if esc prob gt 1 then line is masing - return small number for dest prob */
761  /* >>>chng 99 apr 10, return min dest when scattering greater than abs */
762  /* no idea of opacity whatsoever, on very first soln for this model */
763  /* >>chng 05 mar 20, add test on line being above upper bound of frequency
764  * do not want to evaluate opacity in this case since off scale */
765  if( escp >= 1.0 || !conv.nTotalIoniz || ipanu >= rfield.nflux )
766  {
767  eovrlp_v = 0.;
768  return( eovrlp_v );
769  }
770 
771  /* find continuum opacity */
772  conopc = opac.opacity_abs[ipanu-1];
773 
774  ASSERT( crsec > 0. );
775 
776  /* may be no population, cannot use above test since return 0 not DEST0 */
777  if( abund <= 0. || conopc <= 0. )
778  {
779  /* do not set this to DEST0 since energy not then conserved */
780  eovrlp_v = 0.;
781  return( eovrlp_v );
782  }
783 
784  /* fac of 1.7 convert to Hummer convention for line opacity */
785  beta = conopc/(abund*SQRTPI*crsec/widl + conopc);
786  /* >>chng 04 may 10, rm * 1-pesc)
787  beta = MIN2(beta,(1.-escp)); */
788 
789  if( nCore == ipDEST_INCOM )
790  {
791  /* fits to
792  * >>>refer la esc Hummer and Kunasz 1980 Ap.J. 236,609.
793  * the max value of 1e-3 is so that we do not go too far
794  * beyond what Hummer and Kunasz did, discussed in
795  * >>refer rt esc proc Ferland, G.J., 1999, ApJ, 512, 247 */
798  eovrlp_v = MIN2(1e-3,8.5*beta);
799  }
800  else if( nCore == ipDEST_K2 )
801  {
802  /* Doppler core only; a=0., Hummer 68
803  eovrlp_v = RT_DestHummer(beta);*/
804  eovrlp_v = MIN2(1e-3,8.5*beta);
805  }
806  else if( nCore == ipDEST_SIMPL )
807  {
808  /* this for debugging only
809  eovrlp_v = 8.5*beta;*/
810  /* >>chng 04 may 13, use same min function */
811  eovrlp_v = MIN2(1e-3,8.5*beta);
812  }
813  else
814  {
815  fprintf( ioQQQ, " chCore of %i not understood by RT_DestProb.\n",
816  nCore );
817  cdEXIT(EXIT_FAILURE);
818  }
819 
820  /* renorm to unity */
821  eovrlp_v /= 1. + eovrlp_v;
822 
823  /* multiply by 1-escape prob, since no destruction when optically thin */
824  eovrlp_v *= 1. - escp;
825 
826  /*check results in bounds */
827  ASSERT( eovrlp_v >= 0. );
828  ASSERT( eovrlp_v <= 1. );
829 
830  {
831  /* debugging code for Lya problems */
832  /*@-redef@*/
833  enum {DEBUG_LOC=false};
834  /*@+redef@*/
835  if( DEBUG_LOC )
836  {
837  if( rfield.anu[ipanu-1]>0.73 && rfield.anu[ipanu-1]<0.76 &&
838  fp_equal( abund, Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->PopOpc ) )
839  {
840  fprintf(ioQQQ,"%li RT_DestProb\t%g\n",
841  nzone, eovrlp_v );
842  }
843  }
844  }
845 
846  /* >>chng 04 may 10, rm min */
847  /* this hack removed since no fundamental reason for it to be here,
848  * this should be added to scattering escape, if included at all */
849 # if 0
850  /* >>chng 99 apr 12, limit destruction prob in case where gas dominated by scattering.
851  * in this case scattering is much more likely than absorption on this event
852  eovrlp_v = (1. - opac.albedo[ipanu-1]) * eovrlp_v +
853  opac.albedo[ipanu-1]*DEST0; */
854  /* >>chng 01 aug 11, add factor of 3 for increase in mean free path, and min on 0 */
855  /*eovrlp_v = MAX2(DEST0,1. - 3.*opac.albedo[ipanu-1]) * eovrlp_v +
856  opac.albedo[ipanu-1]*DEST0;*/
857  eovrlp_v = POW2(1. - opac.albedo[ipanu-1]) * eovrlp_v +
858  opac.albedo[ipanu-1]*0.;
859 # endif
860 
861  return( eovrlp_v );
862 }
863 
864 /*RT_LineWidth compute line width (cm/sec), using optical depth array information
865  * this is where effects of wind are done */
866 double RT_LineWidth(const transition * t)
867 {
868  double RT_LineWidth_v,
869  aa,
870  atau,
871  b,
872  r,
873  vth;
874  realnum tau;
875 
876  DEBUG_ENTRY( "RT_LineWidth()" );
877 
878  /* uses line width from
879  * >>refer esc prob Bonilha et al. Ap.J. (1979) 233 649
880  * return value is half velocity width*(1-ESC PROB) [cm s-1]
881  * this assumes incomplete redistribution, damp.tau^1/3 width */
882 
883  /* thermal width */
884  vth = DoppVel.doppler[ t->Hi->nelem-1 ];
885 
886  /* optical depth in outer direction is defined
887  * on second and later iterations.
888  * smaller of inner and outer optical depths is chosen for esc prob */
889  if( iteration > 1 )
890  {
891  /* optical depth to shielded face */
892  realnum tauout = t->Emis->TauTot - t->Emis->TauIn;
893 
894  /* >>chng 99 apr 22 use smaller optical depth */
895  tau = MIN2( t->Emis->TauIn , tauout );
896  }
897  else
898  {
899  tau = t->Emis->TauIn;
900  }
901  /* do not evaluate line width if quite optically thin - will be dominated
902  * by noise in this case */
903  if( tau <1e-3 )
904  return 0;
905 
906  double Pesc = esc_PRD_1side( tau , t->Emis->damp);
907 
908  /* max optical depth is thermalization length */
909  realnum therm = (realnum)(5.3e16/dense.eden);
910  if( tau > therm )
911  {
912  pressure.lgPradDen = true;
913  tau = therm;
914  }
915 
916  /* >>chng 01 jun 23, use wind vel instead of rt since rt deleted */
917  /* >>chng 04 may 13, use thermal for subsonic cases */
922  if( wind.windv <= 0. )
923  {
924  /* static geometry */
925  /* esc prob has noise if smaller than FLT_EPSILON, or is masing */
926  if( (tau-opac.taumin)/100. < FLT_EPSILON )
927  {
928  RT_LineWidth_v = 0.;
929  }
930  else if( tau <= 20. )
931  {
932  atau = -6.907755;
933  if( tau > 1e-3 )
934  atau = log(tau);
935  aa = 4.8 + 5.2*tau + (4.*tau - 1.)*atau;
936  b = 6.5*tau - atau;
937  RT_LineWidth_v = vth*0.8862*aa/b*(1. - MIN2( 1.,
938  ( Pesc + t->Emis->Pelec_esc + t->Emis->Pdest) ) );
939  /* small number roundoff can dominate this process */
940  if( Pesc < 10. * FLT_EPSILON )
941  RT_LineWidth_v = 0.;
942  }
943  else
944  {
945  ASSERT( t->Emis->damp*tau >= 0.);
946  atau = log(MAX2(0.0001,tau));
947  aa = 1. + 2.*atau/pow(1. + 0.3*t->Emis->damp*tau,0.6667) + pow(6.5*
948  t->Emis->damp*tau,0.333);
949  b = 1.6 + 1.5/(1. + 0.20*t->Emis->damp*tau);
950  RT_LineWidth_v = vth*0.8862*aa/b*(1. - MIN2( 1. ,
951  (Pesc+ t->Emis->Pelec_esc + t->Emis->Pdest)) );
952  }
953 
954  /* we want full width, not half width */
955  RT_LineWidth_v *= 2.;
956 
957  }
958  else
959  {
960  /* wind */
961  r = t->Emis->damp*tau/PI;
962  if( r <= 1. )
963  {
964  RT_LineWidth_v = vth*sqrt(log(MAX2(1.,tau))*.2821);
965  }
966  else
967  {
968  RT_LineWidth_v = 2.*fabs(wind.windv0);
969  if( r*vth <= RT_LineWidth_v )
970  {
971  RT_LineWidth_v = vth*r*log(RT_LineWidth_v/(r*vth));
972  }
973  }
974  }
975 
976  ASSERT( RT_LineWidth_v >= 0. );
977  return( RT_LineWidth_v );
978 }
979 
980 /*RT_DestHummer evaluate Hummer's betaF(beta) function */
981 double RT_DestHummer(double beta) /* beta is ratio of continuum to mean line opacity,
982  * returns dest prob = beta F(beta) */
983 {
984  double fhummr_v,
985  x;
986 
987  DEBUG_ENTRY( "RT_DestHummer()" );
988 
989  /* evaluates Hummer's F(beta) function for case where damping
990  * constant is zero, are returns beta*F(beta)
991  * fit to Table 1, page 80, of Hummer MNRAS 138, 73-108.
992  * beta is ratio of continuum to line opacity; FUMMER is
993  * product of his F() times beta; the total destruction prob
994  * this beta is Hummer's normalization of the Voigt function */
995 
996  ASSERT( beta >= 0.);/* non-positive is unphysical */
997  if( beta <= 0. )
998  {
999  fhummr_v = 0.;
1000  }
1001  else
1002  {
1003  x = log10(beta);
1004  if( x < -5.5 )
1005  {
1006  fhummr_v = 3.8363 - 0.56329*x;
1007  }
1008  else if( x < -3.5 )
1009  {
1010  fhummr_v = 2.79153 - 0.75325*x;
1011  }
1012  else if( x < -2. )
1013  {
1014  fhummr_v = 1.8446 - 1.0238*x;
1015  }
1016  else
1017  {
1018  fhummr_v = 0.72500 - 1.5836*x;
1019  }
1020  fhummr_v *= beta;
1021  }
1022  return( fhummr_v );
1023 }

Generated for cloudy by doxygen 1.8.4