cloudy trunk
Loading...
Searching...
No Matches
rt_stark.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 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
12STATIC realnum strkar( long nLo, long nHi );
13
14void RT_stark(void)
15{
16 long int ipLo,
17 ipHi;
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( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
32 {
33 /* loop over all iso-electronic sequences */
34 for( long nelem=ipISO; nelem<LIMELM; ++nelem )
35 {
36 if( nelem >= 2 && !dense.lgElmtOn[nelem] )
37 continue;
38
39 t_iso_sp* sp = &iso_sp[ipISO][nelem];
40
41 if( !rt.lgStarkON )
42 {
43 for( ipHi=0; ipHi < sp->numLevels_max; ipHi++ )
44 {
45 for( ipLo=0; ipLo < sp->numLevels_max; ipLo++ )
46 {
47 sp->ex[ipHi][ipLo].pestrk = 0.;
48 sp->ex[ipHi][ipLo].pestrk_up = 0.;
49 }
50 }
51 continue;
52 }
53
54 /* evaluate Stark escape probability from
55 * >>ref Puetter Ap.J. 251, 446. */
56
57 /* coefficients for Stark broadening escape probability
58 * to be Puetters AH, equation 9b, needs factor of (Z^-4.5 * (nu*nl)^3 * xl) */
59 ah = 6.9e-6*1000./1e12/(phycon.sqrte*phycon.te10*phycon.te10*
60 phycon.te03*phycon.te01*phycon.te01)*dense.eden;
61
62 /* include Z factor */
63 ah *= pow( (double)(nelem+1), -4.5 );
64
65 /* coefficient for all lines except Ly alpha */
66 /* equation 10b, except missing tau^-0.6 */
67 stark = 0.264*pow(ah,0.4);
68
69 /* coefficient for Ly alpha */
70 /* first few factors resemble equation 13c...what about the rest? */
71 strkla = 0.538*ah*4.*9.875*(phycon.sqrte/phycon.te10/phycon.te03);
72
73 long ipHi = iso_ctrl.nLyaLevel[ipISO];
74 /* Lyman lines always have outer optical depths */
75 /* >>chng 02 mar 31, put in max, crashed on some first iteration
76 * with negative optical depths,
77 * NB did not understand why neg optical depths started */
78 aa = (realnum)SDIV(sp->trans(ipHi,0).Emis().TauIn());
79 aa = pow( aa ,-0.75);
80 sp->ex[ipHi][0].pestrk_up = strkla/2.*MAX2(1.,aa);
81
86 sp->ex[ipHi][0].pestrk_up =
87 MIN2(.01,sp->ex[ipHi][0].pestrk_up);
88 sp->ex[ipHi][0].pestrk_up = 0.;
89 sp->ex[ipHi][0].pestrk = sp->ex[ipHi][0].pestrk_up *
90 sp->trans(ipHi,0).Emis().Aul();
91
92 /* >>chng 06 aug 28, from numLevels_max to _local. */
93 for( ipHi=3; ipHi < sp->numLevels_local; ipHi++ )
94 {
95 if( sp->trans(ipHi,0).ipCont() <= 0 )
96 continue;
97
98 sp->ex[ipHi][0].pestrk_up = stark / 2. *
99 strkar( sp->st[0].n(), sp->st[ipHi].n() ) *
100 pow(MAX2(1.,sp->trans(ipHi,0).Emis().TauIn()),-0.75);
101
102 sp->ex[ipHi][0].pestrk_up = MIN2(.01,sp->ex[ipHi][0].pestrk_up);
103 sp->ex[ipHi][0].pestrk = sp->trans(ipHi,0).Emis().Aul()*
104 sp->ex[ipHi][0].pestrk_up;
105 }
106
107 /* zero out rates above iso.numLevels_local */
108 for( ipHi=sp->numLevels_local; ipHi < sp->numLevels_max; ipHi++ )
109 {
110 sp->ex[ipHi][0].pestrk_up = 0.;
111 sp->ex[ipHi][0].pestrk = 0.;
112 }
113
114 /* all other lines */
115 for( ipLo=ipH2s; ipLo < (sp->numLevels_local - 1); ipLo++ )
116 {
117 for( ipHi=ipLo + 1; ipHi < sp->numLevels_local; ipHi++ )
118 {
119 if( sp->trans(ipHi,ipLo).ipCont() <= 0 )
120 continue;
121
122 aa = stark *
123 strkar( sp->st[ipLo].n(), sp->st[ipHi].n() ) *
124 pow(MAX2(1.,sp->trans(ipHi,ipLo).Emis().TauIn()),-0.75);
125 sp->ex[ipHi][ipLo].pestrk_up =
126 (realnum)MIN2(.01,aa);
127
128 sp->ex[ipHi][ipLo].pestrk = sp->trans(ipHi,ipLo).Emis().Aul()*
129 sp->ex[ipHi][ipLo].pestrk_up;
130 }
131 }
132
133 /* zero out rates above iso.numLevels_local */
134 for( ipLo=(sp->numLevels_local - 1); ipLo<(sp->numLevels_max - 1); ipLo++ )
135 {
136 for( ipHi=ipLo + 1; ipHi < sp->numLevels_max; ipHi++ )
137 {
138 sp->ex[ipHi][ipLo].pestrk_up = 0.;
139 sp->ex[ipHi][ipLo].pestrk = 0.;
140 }
141 }
142 }
143 }
144
145 return;
146}
147
148STATIC realnum strkar( long nLo, long nHi )
149{
150 return (realnum)pow((realnum)( nLo * nHi ),(realnum)1.2f);
151}
152
long int nzone
Definition cddefines.cpp:14
#define STATIC
Definition cddefines.h:97
#define MIN2
Definition cddefines.h:761
const int LIMELM
Definition cddefines.h:258
const int NISO
Definition cddefines.h:261
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
sys_float SDIV(sys_float x)
Definition cddefines.h:952
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
realnum & Aul() const
Definition emission.h:613
realnum & TauIn() const
Definition emission.h:423
long & ipCont() const
Definition transition.h:450
EmissionList::reference Emis() const
Definition transition.h:408
long int numLevels_max
Definition iso.h:493
multi_arr< extra_tr, 2 > ex
Definition iso.h:449
TransitionProxy trans(const long ipHi, const long ipLo)
Definition iso.h:444
long int numLevels_local
Definition iso.h:498
qList st
Definition iso.h:453
t_dense dense
Definition dense.cpp:24
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipH2s
Definition iso.h:28
const int ipH_LIKE
Definition iso.h:62
t_phycon phycon
Definition phycon.cpp:6
t_rt rt
Definition rt.cpp:5
void RT_stark(void)
Definition rt_stark.cpp:14
STATIC realnum strkar(long nLo, long nHi)
Definition rt_stark.cpp:148