cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
rt_stark.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_stark compute stark broadening escape probabilities using Puetter formalism */
4 #include "cddefines.h"
5 #include "taulines.h"
6 #include "iso.h"
7 #include "dense.h"
8 #include "hydrogenic.h"
9 #include "phycon.h"
10 #include "rt.h"
11 
12 void RT_stark(void)
13 {
14  long int ipLo,
15  ipHi,
16  nelem,
17  ipISO;
18 
19  double aa , ah,
20  stark,
21  strkla;
22 
23  DEBUG_ENTRY( "RT_stark()" );
24 
25  /* only evaluate one time per zone */
26  static long int nZoneEval=-1;
27  if( nzone==nZoneEval )
28  return;
29  nZoneEval = nzone;
30 
31  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
32  {
33  /* loop over all iso-electronic sequences */
34  for( nelem=ipISO; nelem<LIMELM; ++nelem )
35  {
36  if( nelem >= 2 && !dense.lgElmtOn[nelem] )
37  continue;
38 
39  if( !rt.lgStarkON || dense.eden < 1e8 )
40  {
41  for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
42  {
43  for( ipLo=0; ipLo < iso.numLevels_max[ipISO][nelem]; ipLo++ )
44  {
45  iso.pestrk[ipISO][nelem][ipHi][ipLo] = 0.;
46  iso.pestrk[ipISO][nelem][ipLo][ipHi] = 0.;
47  }
48  }
49  continue;
50  }
51 
52  /* evaluate Stark escape probability from
53  * >>ref Puetter Ap.J. 251, 446. */
54 
55  /* coefficients for Stark broadening escape probability
56  * to be Puetters AH, equation 9b, needs factor of (Z^-4.5 * (nu*nl)^3 * xl) */
57  ah = 6.9e-6*1000./1e12/(phycon.sqrte*phycon.te10*phycon.te10*
59 
60  /* include Z factor */
61  ah *= pow( (double)(nelem+1), -4.5 );
62 
63  /* coefficient for all lines except Ly alpha */
64  /* equation 10b, except missing tau^-0.6 */
65  stark = 0.264*pow(ah,0.4);
66 
67  /* coefficient for Ly alpha */
68  /* first few factors resemble equation 13c...what about the rest? */
69  strkla = 0.538*ah*4.*9.875*(phycon.sqrte/phycon.te10/phycon.te03);
70 
71  /* Lyman lines always have outer optical depths */
72  /*ASSERT( Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauIn> 0. );*/
73  /* >>chng 02 mar 31, put in max, crashed on some first iteration
74  * with negative optical depths,
75  * NB did not understand why neg optical depths started */
76  aa = (realnum)SDIV(Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn);
77  aa = pow( aa ,-0.75);
78  iso.pestrk[ipISO][nelem][0][iso.nLyaLevel[ipISO]] = strkla/2.*MAX2(1.,aa);
79 
84  iso.pestrk[ipISO][nelem][0][iso.nLyaLevel[ipISO]] = MIN2(.01,iso.pestrk[ipISO][nelem][0][iso.nLyaLevel[ipISO]]);
85  iso.pestrk[ipISO][nelem][0][iso.nLyaLevel[ipISO]] = 0.;
86  iso.pestrk[ipISO][nelem][iso.nLyaLevel[ipISO]][0] =
87  iso.pestrk[ipISO][nelem][0][iso.nLyaLevel[ipISO]]*Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Aul;
88 
89 
90  /* >>chng 06 aug 28, from numLevels_max to _local. */
91  for( ipHi=3; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ )
92  {
93  if( Transitions[ipISO][nelem][ipHi][ipH1s].ipCont <= 0 )
94  continue;
95 
96  iso.pestrk[ipISO][nelem][0][ipHi] = stark*iso.strkar[ipISO][nelem][0][ipHi]/2.*pow(MAX2(1.,
97  Transitions[ipISO][nelem][ipHi][ipH1s].Emis->TauIn),-0.75);
98 
99  iso.pestrk[ipISO][nelem][0][ipHi] = MIN2(.01,iso.pestrk[ipISO][nelem][0][ipHi]);
100  iso.pestrk[ipISO][nelem][ipHi][0] = Transitions[ipISO][nelem][ipHi][ipH1s].Emis->Aul*
101  iso.pestrk[ipISO][nelem][0][ipHi];
102  }
103 
104  /* zero out rates above iso.numLevels_local */
105  for( ipHi=iso.numLevels_local[ipISO][nelem]; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
106  {
107  iso.pestrk[ipISO][nelem][0][ipHi] = 0.;
108  iso.pestrk[ipISO][nelem][ipHi][0] = 0.;
109  }
110 
111  /* all other lines */
112  for( ipLo=ipH2s; ipLo < (iso.numLevels_local[ipISO][nelem] - 1); ipLo++ )
113  {
114  for( ipHi=ipLo + 1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ )
115  {
116  if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 )
117  continue;
118 
119  aa = stark*iso.strkar[ipISO][nelem][ipLo][ipHi]*
120  pow(MAX2(1.,Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn),-0.75);
121  iso.pestrk[ipISO][nelem][ipLo][ipHi] =
122  (realnum)MIN2(.01,aa);
123 
124  iso.pestrk[ipISO][nelem][ipHi][ipLo] = Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul*
125  iso.pestrk[ipISO][nelem][ipLo][ipHi];
126  }
127  }
128 
129  /* zero out rates above iso.numLevels_local */
130  for( ipLo=(iso.numLevels_local[ipISO][nelem] - 1); ipLo<(iso.numLevels_max[ipISO][nelem] - 1); ipLo++ )
131  {
132  for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
133  {
134  iso.pestrk[ipISO][nelem][ipLo][ipHi] = 0.;
135  iso.pestrk[ipISO][nelem][ipHi][ipLo] = 0.;
136  }
137  }
138  }
139  }
140 
141  return;
142 }

Generated for cloudy by doxygen 1.8.4