cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
rt_continuum_shield_fcn.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_continuum_shield_fcn computing continuum shielding due to single line */
4 /*conpmp local continuum pumping rate radiative transfer for all lines */
5 /*conpmp local continuum pumping rate radiative transfer for all lines */
6 #include "cddefines.h"
7 #include "rt.h"
8 /*conpmp local continuum pumping rate radiative transfer for all lines */
9 STATIC double conpmp(transition * t);
10 STATIC inline double FITTED( double t );
11 static double PumpDamp , PumpTau;
12 
13 /*rt_continuum_shield_fcn computing continuum shielding due to single line */
15 {
16  double value;
17 
18  DEBUG_ENTRY( "rt_continuum_shield_fcn()" );
19 
20  value = -1.;
21 
23  {
24  /* set continuum shielding pesc - shieding based on escape probability */
25  if( t->Emis->iRedisFun == ipPRD )
26  {
27  value = esc_PRD_1side(t->Emis->TauCon,t->Emis->damp);
28  }
29  else if( t->Emis->iRedisFun == ipCRD )
30  {
31  value = esca0k2(t->Emis->TauCon);
32  }
33  else if( t->Emis->iRedisFun == ipCRDW )
34  {
35  value = esc_CRDwing_1side(t->Emis->TauCon,t->Emis->damp);
36  }
37  else if( t->Emis->iRedisFun == ipLY_A )
38  {
39  value = esc_PRD_1side(t->Emis->TauCon,t->Emis->damp);
40  }
41  else
42  TotalInsanity();
43  }
45  {
46  /* set continuum shielding Federman - this is the default */
47  double core, wings;
48 
49  /* these expressions implement the appendix of
50  * >>refer line shielding Federman, S.R., Glassgold, A.E., &
51  * >>refercon Kwan, J. 1979, ApJ, 227, 466 */
52  /* doppler core - equation A8 */
53  if( t->Emis->TauCon < 2. )
54  {
55  core = sexp( t->Emis->TauCon * 0.66666 );
56  }
57  else if( t->Emis->TauCon < 10. )
58  {
59  core = 0.638 * pow(t->Emis->TauCon,(realnum)-1.25f );
60  }
61  else if( t->Emis->TauCon < 100. )
62  {
63  core = 0.505 * pow(t->Emis->TauCon,(realnum)-1.15f );
64  }
65  else
66  {
67  core = 0.344 * pow(t->Emis->TauCon,(realnum)-1.0667f );
68  }
69 
70  /* do we add damping wings? */
71  wings = 0.;
72  if( t->Emis->TauCon > 1.f && t->Emis->damp>0. )
73  {
74  /* equation A6 */
75  double t1 = 3.02*pow(t->Emis->damp*1e3,-0.064 );
76  double u1 = sqrt(t->Emis->TauCon*t->Emis->damp )/SDIV(t1);
77  wings = t->Emis->damp/SDIV(t1)/sqrt( 0.78540 + POW2(u1) );
78  /* add very large optical depth tail to converge this with respect
79  * to escape probabilities - if this function falls off more slowly
80  * than escape probability then upper level will become overpopulated.
81  * original paper was not intended for this regime */
82  if( t->Emis->TauCon>1e7 )
83  wings *= pow( t->Emis->TauCon/1e7,-1.1 );
84  }
85  value = core + wings;
86  /* some x-ray lines have vastly large damping constants, greater than 1.
87  * in these cases the previous wings value does not work - approximation
88  * is for small damping constant - do not let pump efficiency exceed unity
89  * in this case */
90  if( t->Emis->TauCon>0. )
91  value = MIN2(1., value );
92  }
94  {
95  /* set continuum shielding ferland */
96  value = conpmp( t );
97  }
98  else if( rt.nLineContShield == 0 )
99  {
100  /* set continuum shielding none */
101  value = 1.;
102  }
103  else
104  {
105  TotalInsanity();
106  }
107 
108  /* the returned pump shield function must be greater than zero,
109  * and less than 1 if a maser did not occur */
110  ASSERT( value>=0 && (value<=1.||t->Emis->TauCon<0.) );
111 
112  return value;
113 }
114 
115 /*opfun routine used to get continuum pumping of lines
116  * used in conpmp in call to qg32 */
117 STATIC double opfun(double x)
118 {
119  double opfun_v,
120  v;
121 
122  DEBUG_ENTRY( "opfun()" );
123 
124  v = vfun(PumpDamp,x);
125  opfun_v = sexp(PumpTau*v)*v;
126  return( opfun_v );
127 }
128 
129 static const double BREAK = 3.;
130 /* fit to results for tau less than 10 */
131 // #define FITTED(t) ((0.98925439 + 0.084594094*(t))/(1. + (t)*(0.64794212 + (t)*0.44743976)))
132 STATIC inline double FITTED( double t )
133 {
134  return (0.98925439 + 0.084594094*t)/(1. + t*(0.64794212 + t*0.44743976));
135 }
136 
137 /*conpmp local continuum pumping rate radiative transfer for all lines */
139 {
140  double a0,
141  conpmp_v,
142  tau,
143  yinc1,
144  yinc2;
145 
146  DEBUG_ENTRY( "conpmp()" );
147 
148  /* tau used will be optical depth in center of next zone
149  * >>chng 96 july 6, had been ipLnTauIn, did not work when sphere static set */
150  tau = t->Emis->TauCon;
151  /* compute pumping probability */
152  if( tau <= 10. )
153  {
154  /* for tau<10 a does not matter, and one fit for all */
155  conpmp_v = FITTED(tau);
156  }
157  else if( tau > 1e6 )
158  {
159  /* this far in winds line opacity well below electron scattering
160  * so ignore for this problem */
161  conpmp_v = 0.;
162  }
163  else
164  {
165  /* following two are passed on to later subs */
166  PumpDamp = t->Emis->damp;
167  PumpTau = tau;
168  a0 = 0.886227*(1. + PumpDamp);
169  yinc1 = qg32(0.,BREAK,opfun);
170  yinc2 = qg32(BREAK,100.,opfun);
171  conpmp_v = (yinc1 + yinc2)/a0;
172  }
173 
174  /* EscProb is escape probability, will not allow conpmp to be greater than it
175  * on second iteration with thick lines, pump prob=1 and esc=0.5
176  * conpmp = MIN( conpmp , t->t(ipLnEscP) )
177  * */
178  return( conpmp_v );
179 }

Generated for cloudy by doxygen 1.8.3.1