cloudy trunk
Loading...
Searching...
No Matches
rt_line_driving.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_line_driving derive radiative acceleration due to line absorption of incident continuum,
4 * return value is line radiative acceleration */
5#include "cddefines.h"
6#include "physconst.h"
7#include "iso.h"
8#include "dense.h"
9#include "taulines.h"
10#include "h2.h"
11#include "atomfeii.h"
12#include "rt.h"
13
14/*RT_line_driving derive radiative acceleration due to line absorption of incident continuum,
15 * return value is line radiative acceleration */
16double RT_line_driving(void)
17{
18 long int i,
19 ipHi,
20 nelem,
21 ipLo,
22 ipISO;
23
24 double AllHeavy,
25 AllRest,
26 OneLine,
27 fe2drive,
28 forlin_v,
29 h2drive,
30 accel_iso[NISO];
31
32 /* following used for debugging */
33 /* double
34 RestMax,
35 HeavMax,
36 hydromax;
37 long int
38 ipRestMax,
39 ihmax; */
40
41 DEBUG_ENTRY( "RT_line_driving()" );
42
43 /* this function finds the total rate the gas absorbs energy
44 * this result is divided by the calling routine to find the
45 * momentum absorbed by the gas, and eventually the radiative acceleration
46 *
47 * the energy absorbed by the line is
48 * Abundance * energy * A *(g_up/g_lo) * occnum * escape prob
49 * where occnum is the photon occupation number, and the g's are
50 * the ratios of statistical weights */
51
52 /* total energy absorbed in this zone, per cubic cm
53 * do hydrogen first */
54
55 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
56 {
57 accel_iso[ipISO] = 0;
58 for( nelem=ipISO; nelem < LIMELM; nelem++ )
59 {
60 if( (dense.IonHigh[nelem] >= nelem + 1-ipISO) )
61 {
62 for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
63 {
64 /* do not put in highest level since its not real */
65 for( ipLo=0; ipLo < ipHi - 1; ipLo++ )
66 {
67 /* do not include bogus lines */
68 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() > 0 )
69 {
70 OneLine = iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().pump()*
71 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyErg()*
72 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc();
73
74 accel_iso[ipISO] += OneLine;
75 }
76 }
77 }
78
79 if( iso_ctrl.lgDielRecom[ipISO] )
80 {
81 // SatelliteLines are indexed by lower level, summed over satellite levels
82 for( ipLo=0; ipLo < iso_sp[ipISO][nelem].numLevels_local; ipLo++ )
83 {
84 /* do not include bogus lines */
85 TransitionList::iterator tr = SatelliteLines[ipISO][nelem].begin()+ipSatelliteLines[ipISO][nelem][ipLo];
86 if((*tr).ipCont() > 0 )
87 {
88 OneLine = (*tr).Emis().pump()*
89 (*tr).EnergyErg()*
90 (*tr).Emis().PopOpc();
91
92 accel_iso[ipISO] += OneLine;
93 }
94 }
95 }
96
97 for( ipHi=iso_sp[ipISO][nelem].st[iso_sp[ipISO][nelem].numLevels_local-1].n()+1; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )
98 {
99 /* do not include bogus lines */
100 TransitionList::iterator tr = ExtraLymanLines[ipISO][nelem].begin()+ipExtraLymanLines[ipISO][nelem][ipHi];
101 if( (*tr).ipCont() > 0 )
102 {
103 OneLine = (*tr).Emis().pump()*
104 (*tr).EnergyErg()*
105 (*tr).Emis().PopOpc();
106
107 accel_iso[ipISO] += OneLine;
108 }
109
110 }
111 }
112 }
113 }
114
115 /* all heavy element lines in calculation of cooling
116 * these are the level 1 lines */
117 AllHeavy = 0.;
118 for( i=1; i <= nLevel1; i++ )
119 {
120 OneLine =
121 TauLines[i].Emis().pump()*
122 TauLines[i].EnergyErg()*
123 TauLines[i].Emis().PopOpc();
124 AllHeavy += OneLine;
125 }
126
127 /* all heavy element lines treated with g-bar
128 * these are the level 2 lines, f should be ok */
129 AllRest = 0.;
130 for( i=0; i < nWindLine; i++ )
131 {
132 OneLine =
133 TauLine2[i].Emis().pump()*
134 TauLine2[i].EnergyErg()*
135 TauLine2[i].Emis().PopOpc();
136 AllRest += OneLine;
137 }
138 for( i=0; i < nUTA; i++ )
139 {
140 OneLine =
141 UTALines[i].Emis().pump()*
142 UTALines[i].EnergyErg()*
143 UTALines[i].Emis().PopOpc();
144 AllRest += OneLine;
145 }
146 for( i=0; i < nHFLines; i++ )
147 {
148 OneLine =
149 HFLines[i].Emis().pump()*
150 HFLines[i].EnergyErg()*
151 HFLines[i].Emis().PopOpc();
152 AllRest += OneLine;
153 }
154
155 for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
156 {
157 if( dBaseSpecies[ipSpecies].lgActive )
158 {
159 for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin();
160 tr != dBaseTrans[ipSpecies].end(); ++tr)
161 {
162 int ipHi = (*tr).ipHi();
163 if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local || (*tr).ipCont() <= 0)
164 continue;
165 OneLine = (*tr).EnergyErg()*(*tr).Emis().pump()*(*tr).Emis().PopOpc();
166 AllRest += OneLine;
167 }
168 }
169 }
170
171 /* the H2 molecule */
172 h2drive = 0.;
173 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
174 h2drive += (*diatom)->H2_Accel();
175
176 /* The large model FeII atom */
177 fe2drive = 0.;
178 FeIIAccel(&fe2drive);
179
180 forlin_v = AllHeavy + accel_iso[ipH_LIKE] + accel_iso[ipHE_LIKE] +
181 fe2drive + h2drive + AllRest;
182
183 /*fprintf(ioQQQ," wind te %e %e %e %e %e\n",
184 AllHeavy , HydroAccel , fe2drive , he1l , AllRest );*/
185 return( forlin_v );
186}
void FeIIAccel(double *fe2drive)
const int LIMELM
Definition cddefines.h:258
const int NISO
Definition cddefines.h:261
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long nWindLine
Definition cdinit.cpp:19
TransitionProxy::iterator iterator
Definition transition.h:280
t_dense dense
Definition dense.cpp:24
vector< diatomics * > diatoms
Definition h2.cpp:8
vector< diatomics * >::iterator diatom_iter
Definition h2.h:13
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipHE_LIKE
Definition iso.h:63
const int ipH_LIKE
Definition iso.h:62
double RT_line_driving(void)
long int nSpecies
Definition taulines.cpp:21
vector< vector< TransitionList > > SatelliteLines
Definition taulines.cpp:38
vector< vector< TransitionList > > ExtraLymanLines
Definition taulines.cpp:25
TransitionList UTALines("UTALines", &AnonStates)
TransitionList TauLine2("TauLine2", &AnonStates)
vector< TransitionList > dBaseTrans
Definition taulines.cpp:17
TransitionList HFLines("HFLines", &AnonStates)
multi_arr< int, 3 > ipSatelliteLines
Definition taulines.cpp:37
long int nUTA
Definition taulines.cpp:26
long int nLevel1
Definition taulines.cpp:28
long int nHFLines
Definition taulines.cpp:31
multi_arr< int, 3 > ipExtraLymanLines
Definition taulines.cpp:24
TransitionList TauLines("TauLines", &AnonStates)
species * dBaseSpecies
Definition taulines.cpp:14