cloudy  trunk
 All Data Structures Namespaces 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-2022 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 "iso.h"
6 #include "dense.h"
7 #include "phycon.h"
8 #include "rt.h"
9 
10 STATIC realnum strkar( long nLo, long nHi );
11 
12 void RT_stark(void)
13 {
14  long int ipLo,
15  ipHi;
16 
17  double aa , ah,
18  stark,
19  strkla;
20 
21  DEBUG_ENTRY( "RT_stark()" );
22 
23  /* only evaluate one time per zone */
24  static long int nZoneEval=-1;
25  if( nzone==nZoneEval )
26  return;
27  nZoneEval = nzone;
28 
29  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
30  {
31  /* loop over all iso-electronic sequences */
32  for( long nelem=ipISO; nelem<LIMELM; ++nelem )
33  {
34  if( nelem >= 2 && !dense.lgElmtOn[nelem] )
35  continue;
36 
37  t_iso_sp* sp = &iso_sp[ipISO][nelem];
38 
39  if( !rt.lgStarkON )
40  {
41  for( ipHi=0; ipHi < sp->numLevels_max; ipHi++ )
42  {
43  for( ipLo=0; ipLo < sp->numLevels_max; ipLo++ )
44  {
45  sp->ex[ipHi][ipLo].pestrk = 0.;
46  sp->ex[ipHi][ipLo].pestrk_up = 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 *= powpq( (double)(nelem+1), -9, 2 );
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  long ipHi = iso_ctrl.nLyaLevel[ipISO];
72  /* Lyman lines always have outer optical depths */
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)max(0.,sp->trans(ipHi,0).Emis().TauIn());
77  aa = powpq( aa, -3, 4 );
78  sp->ex[ipHi][0].pestrk_up = strkla/2.*MAX2(1.,aa);
79 
84  sp->ex[ipHi][0].pestrk_up =
85  MIN2(.01,sp->ex[ipHi][0].pestrk_up);
86  sp->ex[ipHi][0].pestrk_up = 0.;
87  sp->ex[ipHi][0].pestrk = sp->ex[ipHi][0].pestrk_up *
88  sp->trans(ipHi,0).Emis().Aul();
89 
90  /* >>chng 06 aug 28, from numLevels_max to _local. */
91  for( ipHi=3; ipHi < sp->numLevels_local; ipHi++ )
92  {
93  if( sp->trans(ipHi,0).ipCont() <= 0 )
94  continue;
95 
96  sp->ex[ipHi][0].pestrk_up = stark / 2. *
97  strkar( sp->st[0].n(), sp->st[ipHi].n() ) *
98  powpq(MAX2(1.,sp->trans(ipHi,0).Emis().TauIn()),-3,4);
99 
100  sp->ex[ipHi][0].pestrk_up = MIN2(.01,sp->ex[ipHi][0].pestrk_up);
101  sp->ex[ipHi][0].pestrk = sp->trans(ipHi,0).Emis().Aul()*
102  sp->ex[ipHi][0].pestrk_up;
103  }
104 
105  /* zero out rates above iso.numLevels_local */
106  for( ipHi=sp->numLevels_local; ipHi < sp->numLevels_max; ipHi++ )
107  {
108  sp->ex[ipHi][0].pestrk_up = 0.;
109  sp->ex[ipHi][0].pestrk = 0.;
110  }
111 
112  /* all other lines */
113  for( ipLo=ipH2s; ipLo < (sp->numLevels_local - 1); ipLo++ )
114  {
115  for( ipHi=ipLo + 1; ipHi < sp->numLevels_local; ipHi++ )
116  {
117  if( sp->trans(ipHi,ipLo).ipCont() <= 0 )
118  continue;
119 
120  aa = stark *
121  strkar( sp->st[ipLo].n(), sp->st[ipHi].n() ) *
122  powpq(MAX2(1.,sp->trans(ipHi,ipLo).Emis().TauIn()),-3,4);
123  sp->ex[ipHi][ipLo].pestrk_up =
124  (realnum)MIN2(.01,aa);
125 
126  sp->ex[ipHi][ipLo].pestrk = sp->trans(ipHi,ipLo).Emis().Aul()*
127  sp->ex[ipHi][ipLo].pestrk_up;
128  }
129  }
130 
131  /* zero out rates above iso.numLevels_local */
132  for( ipLo=(sp->numLevels_local - 1); ipLo<(sp->numLevels_max - 1); ipLo++ )
133  {
134  for( ipHi=ipLo + 1; ipHi < sp->numLevels_max; ipHi++ )
135  {
136  sp->ex[ipHi][ipLo].pestrk_up = 0.;
137  sp->ex[ipHi][ipLo].pestrk = 0.;
138  }
139  }
140  }
141  }
142 
143  return;
144 }
145 
146 STATIC realnum strkar( long nLo, long nHi )
147 {
148  return (realnum)pow((realnum)( nLo * nHi ),(realnum)1.2f);
149 }
150 
qList st
Definition: iso.h:482
double te03
Definition: phycon.h:58
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
const int NISO
Definition: cddefines.h:311
t_phycon phycon
Definition: phycon.cpp:6
long int nzone
Definition: cddefines.cpp:14
#define MIN2(a, b)
Definition: cddefines.h:803
t_dense dense
Definition: global.cpp:15
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
double te01
Definition: phycon.h:58
#define STATIC
Definition: cddefines.h:118
EmissionList::reference Emis() const
Definition: transition.h:447
STATIC realnum strkar(long nLo, long nHi)
Definition: rt_stark.cpp:146
long & ipCont() const
Definition: transition.h:489
float realnum
Definition: cddefines.h:124
long max(int a, long b)
Definition: cddefines.h:817
void RT_stark(void)
Definition: rt_stark.cpp:12
bool lgElmtOn[LIMELM]
Definition: dense.h:160
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
Definition: iso.h:470
int nLyaLevel[NISO]
Definition: iso.h:397
multi_arr< extra_tr, 2 > ex
Definition: iso.h:478
bool lgStarkON
Definition: rt.h:211
const int ipH2s
Definition: iso.h:30
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:308
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
double powpq(double x, int p, int q)
Definition: service.cpp:630
double te10
Definition: phycon.h:58
double eden
Definition: dense.h:201
#define MAX2(a, b)
Definition: cddefines.h:824
long int numLevels_max
Definition: iso.h:524
double sqrte
Definition: phycon.h:58
realnum & Aul() const
Definition: emission.h:690
long int numLevels_local
Definition: iso.h:529
realnum & TauIn() const
Definition: emission.h:470
t_rt rt
Definition: rt.cpp:5