cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
rt_line_one.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 /*RT_line_static do line radiative transfer,
4  * evaluates escape and destruction probability,
5  * called by HydroPEsc, RT_line_all */
6 /* ND wind distinction is not here - that is done where optical depths are incremented
7  * with line opacity - rt_line_one_tauinc and in line width for rad pressure, in
8  * rt escprob */
9 #include "cddefines.h"
10 #include "rfield.h"
11 #include "doppvel.h"
12 #include "dense.h"
13 #include "opacity.h"
14 #include "lines_service.h"
15 #include "conv.h"
16 #include "radius.h"
17 #include "rt.h"
18 #include "physconst.h"
19 
20 /*RT_line_static do line radiative transfer,
21  * called by RT_line_all */
23  /* the em line we will work on */
24  transition * t ,
25  /* when this is true do both esc and dest probs,
26  * when falst, only dest probs */
27  bool lgDoEsc ,
28  /* this is option to not include line self shielding across this zone.
29  * this can cause pump to depend on zone thickness, and leads to unstable
30  * feedback in some models with the large H2 molecule, due to Solomon
31  * process depending on zone thickness and level populations. */
32  bool lgShield_this_zone );
33 
34 /*RT_line_one do rt for emission line structure - calls RT_line_static or RT_line_wind */
36  /* the em line we will work on */
37  transition *t ,
38  /* when this is true do both esc and dest probs,
39  * when falst, only dest probs */
40  bool lgDoEsc ,
41  /* if true then we want to update the fine opacity scale */
42  bool lgDoFine_opac_update ,
43  /* this is option to not include line self shielding across this zone.
44  * this can cause pump to depend on zone thickness, and leads to unstable
45  * feedback in some models with the large H2 molecule, due to Solomon
46  * process depending on zone thickness and level populations. */
47  bool lgShield_this_zone )
48 {
49  static long int nTau[100],
50  n;
51 
52  DEBUG_ENTRY( "RT_line_one()" );
53 
54  /*>>chng 06 apr 05, effects of simply returning upon call with no pop */
55  /* skip line transfer if requested with 'no line transfer' command, but never skip Lya */
56  if( (t->Lo->Pop<=SMALLFLOAT && conv.nTotalIoniz) ||
57  (!rfield.lgDoLineTrans && (t->Emis->iRedisFun != ipLY_A) ) )
58  {
59  return;
60  }
61 
62  {
63  /* option to keep track of population values during calls,
64  * print out data to make histogram */
65  /*@-redef@*/
66  enum {DEBUG_LOC=false};
67  /*@+redef@*/
68  if( DEBUG_LOC )
69  {
70  if( nzone==0 )
71  {
72  for(n=0; n<100; ++n )
73  nTau[n] = 0;
74  }
75  if( t->Lo->Pop<=SMALLFLOAT )
76  n = 0;
77  else
78  n = (long)log10( t->Lo->Pop )+37;
79  n = MIN2( n , 99 );
80  n = MAX2( n , 0 );
81  /*if( n==37 && nzone > 5 )
82  {
83  fprintf(ioQQQ,"HIT\n");
84  }*/
85  ++nTau[n];
86  if( nzone > 183 )
87  {
88  for(n=0; n<100; ++n )
89  fprintf(ioQQQ,"%li\t%li\n", n , nTau[n] );
90  cdEXIT(EXIT_SUCCESS);
91  }
92  }
93  }
94 
95  /* always call single routine */
97  /* the line we will work on */
98  t ,
99  /* when this is true do both esc and dest probs,
100  * when falst, only dest probs */
101  lgDoEsc ,
102  lgShield_this_zone);
103 
104  /* this is line center frequency, including bulk motion of gas */
105  long int ipLineCenter = t->Emis->ipFine + rfield.ipFineConVelShift;
106 
107  /* define fine opacity fine grid fine mesh */
108  /* rfield.lgOpacityFine flag set flast with no fine opacities command */
109  if( lgDoFine_opac_update && ipLineCenter >= 0 && t->Emis->PopOpc > 0. &&
110  ipLineCenter<rfield.nfine && rfield.lgOpacityFine )
111  {
112  long int i;
113 
114  /* number of fine opacity cells corresponding to one doppler width for current
115  * temperature, velocity field, and nuclear mass,
116  * rfield.fine_opac_velocity_width is width per cell, cm/s */
118 
119  /* line center opacity - want realnum since will add to fine opacity array,
120  * which is realnum */
121  realnum opac_line = (realnum)t->Emis->PopOpc * t->Emis->opacity / DoppVel.doppler[t->Hi->nelem-1];
122 
123  /* this is effective optical depth to this point. Do not work with line if this product
124  * is substantially less than unity */
125  double dTauEffec = opac_line*radius.depth_x_fillfac;
126 
127  /* this is line center for current velocity, which is measured
128  * relative to illuminated face */
129 
130  /* use dTauEffec to choose whether line is important */
131 
132  if( dTauEffec > SMALLFLOAT )
133  {
134  /* width of optically thick line, do 4x with exponential core,
135  * must be at least one cell, but profile is symmetric */
136  long int nCells_core = (long)(cells_wide_1x*4.f+1.5f);
137  realnum opacity;
138 
139  /* core is symmetric - make sure both upper and lower bounds
140  * are within continuum energy bin */
141  if( ipLineCenter - nCells_core < 1 )
142  nCells_core = ipLineCenter - 1;
143  if( ipLineCenter + nCells_core > rfield.nfine )
144  nCells_core = ipLineCenter -rfield.nfine - 1;
145 
146  /* want core to be at least one cell wide */
147  nCells_core = MAX2( 1 , nCells_core );
148 
149  /* line center itself, must not double count here */
150  rfield.fine_opac_zone[ipLineCenter] += opac_line;
151  for( i=1; i<nCells_core; ++i )
152  {
153  /* square of distance from line center in units of doppler width */
154  realnum xsq = POW2(i/cells_wide_1x);
155 
156  /* the line opacity at that point, includes doppler core and damping wings */
157  opacity = opac_line*( sexp(xsq) + t->Emis->damp / MAX2(1.f,(1.772f*xsq)) );
158 
159  rfield.fine_opac_zone[ipLineCenter+i] += opacity;
160  rfield.fine_opac_zone[ipLineCenter-i] += opacity;
161  }
162 
163  /* include damping wings if optical depth is large */
164  /* >>chng 04 mar 29, rewrite to use symmetry in two halves of line */
165  if( dTauEffec*t->Emis->damp/9. > 0.1 )
166  {
167  /* do damping wings if very optically thick
168  * this is number of cells to extend the damping wings - cells_wide is one dop width
169  * 56.419 is 1 / (0.01 * sqrt pi) ) */
170  /* tests with th85orion, stopping at half -h2 point, showed that 0.01 was
171  * needed for outer edge, given the definition of dTauEffec */
172  realnum x = (realnum)sqrt( dTauEffec * t->Emis->damp *56.419 ) * cells_wide_1x;
173  long int nCells_damp;
174  /* >>chng 04 mar 24, add test on size of x, which can exceed
175  * limits of long in extreme optical depths */
176  if( x<LONG_MAX )
177  {
178  nCells_damp = (long)x;
179 
180  if( ipLineCenter-nCells_damp < 1 )
181  nCells_damp = ipLineCenter-1;
182 
183  if( ipLineCenter+nCells_damp > rfield.nfine )
184  nCells_damp = rfield.nfine - ipLineCenter-1;
185  }
186  else
187  {
188  /* x was too big, just set to extreme range, which is
189  * number of cells to nearest boundary */
190  nCells_damp = MIN2( rfield.nfine-ipLineCenter , ipLineCenter )-1;
191  }
192 
193  ASSERT( nCells_core>0 );
194  for( i=nCells_core; i<nCells_damp; ++i )
195  {
196  /* distance from line center in doppler units */
197  realnum xsq = POW2(i/cells_wide_1x);
198  /* term in denominator is x^2 */
199  opacity = opac_line*t->Emis->damp/(1.772f*xsq );
200  rfield.fine_opac_zone[ipLineCenter+i] += opacity;
201  rfield.fine_opac_zone[ipLineCenter-i] += opacity;
202  }
203  }
204  }
205  }
206  return;
207 }
208 
209 /*RT_line_static do line radiative transfer for static geometry */
211  /* the em line we will work on */
212  transition * t ,
213  /* when this is true do both esc and dest probs,
214  * when false, only dest probs */
215  bool lgDoEsc ,
216  /* this is option to not include line self shielding across this zone.
217  * this can cause pump to depend on zone thickness, and leads to unstable
218  * feedback in some models with the large H2 molecule, due to Solomon
219  * process depending on zone thickness and level populations. */
220  bool lgShield_this_zone )
221 {
222  bool lgGoodTau;
223  long int nelem;
224  int nRedis = -1;
225 
226  DEBUG_ENTRY( "RT_line_static()" );
227 
228  /* this is the routine that computes much of the radiative transfer
229  * physics for lines for static case. */
230 
231  /*ASSERT( t->IonStg < LIMELM+4 );*/
232  /*ASSERT( t->nelem-1 < LIMELM+4 );*/
233 
234  /* >>chng 01 aug 18, do not evaluate if no population */
235  /* molecules are evaluated here, these must be valid */
236  /* >>chng 01 oct 25, add test for lgDoEsc, this is to make sure
237  * that rates are evaluated at start of a new iteration */
238  if( dense.xIonDense[ t->Hi->nelem-1][t->Hi->IonStg-1] == 0. && !lgDoEsc )
239  {
240  /* zero population, return after setting everything with side effects */
241  t->Emis->Pesc = 1.f;
242 
243  /* inward escaping fraction */
244  t->Emis->FracInwd = 0.5;
245 
246  /* pumping rate */
247  t->Emis->pump = 0.;
248 
249  /* destruction probability */
250  t->Emis->Pdest = 0.;
251  t->Emis->Pelec_esc = 0.;
252 
253  return;
254  }
255  /* check that we don't have a runaway maser */
256  else if( t->Emis->TauIn < -30. )
257  {
258  fprintf( ioQQQ, "PROBLEM RT_line_static called with large negative optical depth, zone %.2f, setting lgAbort true.\n",
259  fnzone );
260  DumpLine(t);
261  /* >>>chng 00 apr 23, return busted true instead of exit here, so that code will
262  * not hang under MPI */
263  lgAbort = true;
264  return;
265  }
266 
267  /* this checks if we have overrun the optical depth scale,
268  * in which case the inward optical depth is greater than the
269  * previous iteration's total optical depth.
270  * We do not reevaluate escape probabilities if the optical depth
271  * scale has been overrun due to huge bogus change in solution
272  * that would result */
273  if( iteration == 1 || lgTauGood( t ) )
274  {
275  lgGoodTau = true;
276  }
277  else
278  {
279  /* do not reevaluate escape probabilities if we have overrun
280  * the optical depth scale */
281  lgGoodTau = false;
282  }
283 
284  /* atomic number of this element on the C scale*/
285  nelem = t->Hi->nelem-1;
286 
287  /* damping constant for the line, save it */
288  t->Emis->damp = t->Emis->dampXvel / DoppVel.doppler[t->Hi->nelem-1];
289  /*if( t->damp <0. )
290  {
291  fprintf(ioQQQ,"DEBUG hit it \n");
292  }*/
293  ASSERT( t->Emis->damp > 0. ||
294  (opac.lgCaseB && t->Emis->damp >= 0.) );
295 
296  /* static solution - which type of line will determine
297  * which redistribution function */
298  /* iRedisFun == 1 - alpha resonance line, partial redistribution,
299  * ipPRD == 1 */
300  if( t->Emis->iRedisFun == ipPRD )
301  {
302  /* incomplete redistribution with wings */
303  if( lgDoEsc )
304  {
305  if( lgGoodTau )
306  {
307  t->Emis->Pesc = (realnum)esc_PRD( t->Emis->TauIn , t->Emis->TauTot , t->Emis->damp );
308 
309  /* inward escaping fraction */
310  t->Emis->FracInwd = rt.fracin;
311  }
312  }
313  nRedis = ipDEST_INCOM;
314  }
315 
316  /* complete redistribution without wings - t->ipLnRedis is ipCRD == -1 */
317  else if( t->Emis->iRedisFun == ipCRD )
318  {
319  if( lgDoEsc )
320  {
321  if( lgGoodTau )
322  {
323  /* >>chng 01 mar -6, escsub will call any of several esc prob routines,
324  * depending of how core is set. We always want core-only for this option,
325  * so call esca0k2(tau) directly */
326  t->Emis->Pesc = (realnum)esc_CRDcore( t->Emis->TauIn , t->Emis->TauTot );
327 
328  /* inward escaping fraction */
329  t->Emis->FracInwd = rt.fracin;
330  }
331  }
332  nRedis = ipDEST_K2;
333  }
334 
335  /* chng 01 feb 27, add CRD with wings, = 2 */
336  else if( t->Emis->iRedisFun == ipCRDW )
337  {
338  /* complete redistribution with damping wings */
339  if( lgDoEsc )
340  {
341  if( lgGoodTau )
342  {
343  t->Emis->Pesc = (realnum)esc_CRDwing( t->Emis->TauIn , t->Emis->TauTot , t->Emis->damp );
344 
345  /* inward escaping fraction */
346  t->Emis->FracInwd = rt.fracin;
347  }
348  }
349  nRedis = ipDEST_K2;
350  }
351 
352  /* chng 03 may 14, special Lya case*/
353  else if( t->Emis->iRedisFun == ipLY_A )
354  {
355  /* incomplete redistribution with wings, for special case of Lya
356  * uses fits to Hummer & Kunasz numerical results
357  * this routine is different because escape and dest probs
358  * are evaluated together, so no test of lgDoEsc */
359  if( lgGoodTau )
360  {
361  double dest , esin;
362 
363  /* this will always evaluate escape prob, no matter what lgDoEsc is.
364  * The destruction prob comes back as dest */
365  t->Emis->Pesc = (realnum)RTesc_lya(&esin , &dest , t->Emis->PopOpc , nelem);
366 
367  /* this is current destruction rate */
368  t->Emis->Pdest = (realnum)dest;
369 
370  /* this is fraction of line which is inward escaping */
371  t->Emis->FracInwd = rt.fracin;
372  }
373  }
374  else
375  {
376  fprintf( ioQQQ, " RT_line_static called with redistribution function zero\n" );
377  ShowMe();
378  cdEXIT(EXIT_FAILURE);
379  }
380 
381  /* pumping by incident and diffuse continua */
382  /* >>chng 03 may 15, this block had been repeated all four times above, move down
383  * here to use common code, only rt.wayin is defined above */
384  /* >>chng 03 may 17, do not evaluate when drad is zero, before first
385  * zone thickness is established */
386  /* >>chng 06 nov 29, reordered if clauses, no point in doing bigger loop if "no induced" is on. */
387  /* option to kill induced processes */
388  if( !rfield.lgInducProcess )
389  {
390  t->Emis->pump = 0.;
391  }
392  else if( (lgDoEsc || (t->Emis->iRedisFun == ipLY_A) ) )
393  {
394  double dTau;
395  double pumpsave = t->Emis->pump;
396  double shield_continuum;
397 
398  /* >>chng 04 apr 18, break out line continuum shielding into this function */
399  shield_continuum = RT_continuum_shield_fcn( t );
400 
401  /* continuum upward pumping rate, A gu/gl abs prob occnum
402  * the "no induced" command causes continuum pumping to be set to 0 */
403  /*t->pump = t->Aul * t->gHi / t->gLo * rt.wayin **/
404  /* >>chng 05 feb 16, add outward diffuse emission */
405  /*t->pump = t->Aul * t->gHi / t->gLo * shield_continuum *
406  rfield.OccNumbIncidCont[t->ipCont-1];*/
407  t->Emis->pump = t->Emis->Aul * t->Hi->g / t->Lo->g * shield_continuum *(
409 
410  /* >>chng 99 dec 02, add correction for line optical depth across zone */
411  /* correction for attenuation within line across zone */
412  /* >>chng 99 mar 18, from 0 to 1e-8 in following to avoid underflow,
413  * this caused an overestimate of pumping rates */
414  /* >>chng 00 dec 19, nova.in oscillated since correction factor was
415  * very large across zones that were already very optically thick, dr changed,
416  * causing feedback between ionization and thickness
417  * the intention of this correction was to only apply when on the
418  * verge of optically thick,
419  * change lower bound on product t->dTau * radius.drad_x_fillfac to 1e-3 to avoid
420  * constantly getting unity, also check on relative change optical depth/optical depth*/
421  /* this can drive oscillations in very neutral gas, as blr92.in, blr89.in
422  * or nova.in, when higher lyman lines have produce
423  * t->dTau * radius.drad_x_fillfac about unity, where small changes in
424  * pump can change ionization, dTau, and hence pump */
425  /*if( t->dTau * radius.drad_x_fillfac > 1e-8 )*/
426  /* >>chng 02 jun 14, nova.in again - upper limit had been 10 - but
427  * there were major jumps in heating rate at this point - increase 10 -> 1e10
428  * not clear why there should be an upper limit at all, other than to save time? */
429  /* update dTau since affects pump rate */
430  /* >>chng 03 jun 18, include option to not shield across this zone */
431  /* >>chng 03 jul 19, add continuum opacity */
432  dTau = (t->Emis->PopOpc * t->Emis->opacity / DoppVel.doppler[nelem] + opac.opacity_abs[t->ipCont-1])*
434  /* >>chng 04 mar 04, do not want to include maser action across this zone */
435  dTau = MIN2(1., dTau);
436  /*if( 0 && lgShield_this_zone && dTau * radius.drad_x_fillfac > 1e-3 && dTau * radius.drad_x_fillfac < 1e10 )*/
437  /* >>chng 04 mar 04, use mean of past few zone thicknesses rather than
438  * current zone thickness since this can induce wild oscillations, see
439  * nova.in test case */
440  if( lgShield_this_zone && dTau > 1e-3 && dTau < 1e10 )
441  {
442  /* during very first search dTau is -1e38*/
443  /* >>chng 04 mar 04, use mean, see above comment */
444  /*t->pump *= log(1. + dTau * radius.drad_x_fillfac ) / (dTau * radius.drad_x_fillfac);*/
445  t->Emis->pump *= log(1. + dTau ) / dTau;
446  }
447 
448  /* >>chng 00 Oct 03, add pumping by local diffuse continuum,
449  * rfield.DiffPumpOn is unity unless process disabled by setting to 0 in parsedont.c
450  * with no DIFFuse LINE PUMPing command */
451  if( dTau*rfield.DiffPumpOn > SMALLFLOAT )
452  {
453  /* >>chng 02 mar 14, added brackets around scale*rfield.OccNumbDiffCont so
454  * that smallest and largest number are multiplied first to prevent overflow, PvH */
455  double scale;
456  if( iteration > 1 )
457  {
458  if( lgGoodTau )
459  {
460  /* >>chng 03 jul 19, add check for tau = 1 in continuum */
461  /* scale for tau = 1 in the continuum */
462  double scale_con = 1. / SDIV(opac.opacity_abs[t->ipCont-1]+opac.opacity_sct[t->ipCont-1]);
463  /* inward looking depth, dTau is opacity in line, cm^-1,
464  * so multiplying by this scale is same as dividing by opacity in
465  * source function */
466  double scale_in = MIN3( 1./ dTau , radius.depth , scale_con);
467  double scale_out = MIN3( 1./ dTau , radius.Depth2Go , scale_con );
468  t->Emis->pump += t->Emis->Aul * t->Hi->g / t->Lo->g * 0.5 *((scale_in+scale_out)*
470  }
471  }
472  else
473  {
474  /* >>chng 03 jul 19, add check for tau = 1 in continuum */
475  /* scale for tau = 1 in the continuum */
476  double scale_con = 1. / SDIV(opac.opacity_abs[t->ipCont-1]+opac.opacity_sct[t->ipCont-1]);
477  /* first iteration, only inward looking depths really known */
478  scale = 1./ dTau;
479  scale = MIN3( scale , radius.depth , scale_con );
480  t->Emis->pump += t->Emis->Aul * t->Hi->g / t->Lo->g * (scale*
482  }
483  }
484  /* this stop oscillations from developing in very deep blr models
485  * like blr89 and blr98, but must not be applied too early since
486  * line optical depths change greatly in first few zones */
487  if( nzone>50 )
488  t->Emis->pump = (pumpsave + t->Emis->pump) / 2.;
489  }
490 
491 
492  /* escape following scattering off an electron */
493  /* this is turned off with the no scattering escape command */
494  if( !rt.lgElecScatEscape )
495  {
496  t->Emis->Pelec_esc = 0.;
497  }
498  else
499  {
500  /* in the transition structure PopOpc is the pop relative to the ion,
501  * but has been converted to a physical population before this routine
502  * was called. No need to convert here */
503  double opac_line = t->Emis->PopOpc * t->Emis->opacity/DoppVel.doppler[nelem];
504 
505  if( opac_line > SMALLFLOAT )
506  {
507  /* the old escape probability */
508  double eesc_save = t->Emis->Pelec_esc;
509  /* the opacity in electron scattering */
510  double es = dense.eden*6.65e-25;
511  fixit(); // should we use total opacity here instead of opac_line?
512  /* this is equation 5 of
513  *>>refer line desp Netzer, H., Elitzur, M., & Ferland, G. J. 1985, ApJ, 299, 752*/
514  double opacity_ratio = es/(es+opac_line) * (1.-t->Emis->Pesc);
515  /* keep total probability not more than 1. */
516  t->Emis->Pelec_esc = (realnum)opacity_ratio * MAX2(0.f,1.f-t->Emis->Pesc-t->Emis->Pdest);
517  /*KLUDGE >>chng 03 feb 02
518  * take mean of old and new if electron scattering escape */
519  if( !conv.lgSearch )
520  t->Emis->Pelec_esc = t->Emis->Pelec_esc*0.75f + (realnum)eesc_save*0.25f;
521  }
522  else
523  t->Emis->Pelec_esc = 0.;
524  }
525 
526  /* >>>chng 99 dec 18, did not have the test for lgGoodTau, and so clobbered
527  * dest prob when overrun */
528  /* only do this if not Lya special case, since dest prob already done */
529  if( lgGoodTau && t->Emis->iRedisFun!=ipLY_A && t->Emis->opacity > 0. )
530  {
531  double DestSave = t->Emis->Pdest;
532  t->Emis->Pdest = (realnum)RT_DestProb(
533  /* the abundance of the species */
534  t->Emis->PopOpc,
535  /* line center opacity in funny units (needs a vel) */
536  t->Emis->opacity,
537  /* array index for line in continuum array,
538  * on f not c scale */
539  t->ipCont,
540  /* line width in velocity units */
541  DoppVel.doppler[nelem],
542  /* escape probability */
543  t->Emis->Pesc,
544  /* redistribution function */
545  nRedis);
546 
547  /* >>chng 03 may 10, this had been 3/4 new + 1/4 old, but had to be combined
548  ^ with half and half when this routine was called, to converge very high Z high U
549  * quasar models. bring all these factors together here */
550  /*t->Pdest = 0.75f*t->Pdest + 0.25f * DestSave;*/
551 # define OLDFAC (0.625)
552  if( nzone>1 )
553  t->Emis->Pdest = (realnum)((1.-OLDFAC)*t->Emis->Pdest + OLDFAC * DestSave);
554  }
555 
556  /* zero escape probability and pumping term if energy is below plasma frequency */
557  if( t->EnergyErg / EN1RYD < rfield.plsfrq )
558  {
559  t->Emis->Pesc = 0.;
560  t->Emis->Pdest = 0.;
561  t->Emis->pump = 0.;
562  }
563 
564  return;
565 }
566 
567 #if 0
568 /*RT_line_wind do line radiative transfer for wind geometry */
569 STATIC void RT_line_wind(transition * t , bool lgDoEsc);
570 /*RT_line_wind do line radiative transfer for wind geometry */
571 STATIC void RT_line_wind(transition * t , bool lgDoEsc)
572 {
573  realnum velocity;
574 
575  DEBUG_ENTRY( "RT_line_wind()" );
576 
577  /*
578  * this is the routine that computes much of the radiative transfer
579  * physics for lines in wind case. Lines diveded into two types,
580  * high quality (level1) and lower quality (atom_level2) g-bar lines
581  */
582 
583  /* >>chng 01 aug 18, do not evaluation if no population */
584  /* >>chng 01 oct 25, add test for lgDoEsc, this is to make sure
585  * that rates are evaluated at start of a new iteration */
586  if( dense.xIonDense[ t->nelem-1][t->IonStg-1] == 0. && !lgDoEsc )
587  {
588  /* zero population, return after setting everything with side effects */
589  t->Pesc = 1.f;
590 
591  /* inward escaping fraction */
592  t->FracInwd = 0.5;
593 
594  /* pumping rate */
595  t->pump = 0.;
596 
597  /* destruction probability */
598  t->Pdest = 0.;
599  t->Pelec_esc = 0.;
600  return;
601  }
602 
603  /*confirm pops and cs ok */
604  ASSERT( t->ColOvTot <= 1. && t->ColOvTot >= 0. );
605 
606  /* t->nelem is atomic number of species */
607  ASSERT( t->nelem >= 1 && t->ipCont > 0 );
608 
609  /* line width, thermal + turbulence */
610  velocity = DoppVel.doppler[t->nelem-1];
611 
612  t->damp = 0.;
613 
614  /* use only one-sided escape probabilities
615  * get escape probability, this gets fractions too */
616  if( lgDoEsc )
617  {
618  t->Pesc = esc_CRDwing_1side(t->TauIn,t->damp);
619 
620  /* inward fraction */
621  t->FracInwd = 0.5;
622 
623  /* continuum pumping */
624  if( rfield.lgInducProcess )
625  {
626  /* >>chng 05 feb 16, include diffuse outward emission */
627  /*t->pump = t->Aul*t->gHi/t->gLo*t->Pesc*
628  rfield.OccNumbIncidCont[t->ipCont-1];*/
629  t->pump = t->Aul*t->gHi/t->gLo*t->Pesc*(
631  }
632  else
633  {
634  t->pump = 0.;
635  }
636  }
637 
638  /* get destruction probability */
639  t->Pdest =
640  RT_DestProb(
641  t->PopOpc,
642  /* its line absorption cross section */
643  t->opacity,
644  /* index for line in continuum array, on f not c scale */
645  t->ipCont,
646  velocity,
647  t->Pesc,
648  ipDEST_K2);
651  return;
652 }
653 #endif
654 

Generated for cloudy by doxygen 1.8.4