cloudy trunk
Loading...
Searching...
No Matches
rt_tau_inc.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_tau_inc increment optical depths once per zone, called after radius_increment */
4#include "cddefines.h"
5#include "taulines.h"
6#include "iso.h"
7#include "rfield.h"
8#include "trace.h"
9#include "dense.h"
10#include "hyperfine.h"
11#include "wind.h"
12#include "prt.h"
13#include "conv.h"
14#include "h2.h"
15#include "mole.h"
16#include "hmi.h"
17#include "opacity.h"
18#include "cooling.h"
19#include "thermal.h"
20#include "radius.h"
21#include "atomfeii.h"
22#include "rt.h"
23#include "doppvel.h"
24#include "mole.h"
25
26/*RT_tau_inc increment optical depths once per zone, called after radius_increment */
27void RT_tau_inc(void)
28{
29
30 long int i,
31 ipHi,
32 ipLo;
33
34 DEBUG_ENTRY( "RT_tau_inc()" );
35
36 if( trace.lgTrace )
37 {
38 fprintf( ioQQQ, " RT_tau_inc called.\n" );
39 }
40
41 /* call RT_line_all one last time in this zone, to get fine opacities defined */
42 ASSERT( !conv.lgFirstSweepThisZone );
43 conv.lgLastSweepThisZone = true;
44 RT_line_all( );
45
46 /* rfield.lgOpacityFine flag set false with no fine opacities command
47 * tests show that always evaluating this changes fast run of
48 * pn_paris from 26.7 sec to 35.1 sec
49 * but must always update fine opacities since used for transmission */
50 if( rfield.lgOpacityFine )
51 {
52 /* increment the fine opacity array */
53 for( long i=0; i<rfield.nfine; ++i )
54 {
55 realnum tauzone = rfield.fine_opac_zone[i]*(realnum)radius.drad_x_fillfac;
56 rfield.fine_opt_depth[i] += tauzone;
57 }
58 rfield.trans_coef_total_stale = true;
59 }
60
61 /* this may have updated some escape/destruction rates - force update
62 * to all cooling lines */
63 CoolEvaluate( &thermal.ctot );
64
65 if( nzone <=1 )
66 {
67 opac.telec = (realnum)(radius.drad_x_fillfac*dense.eden*6.65e-25);
68 opac.thmin = (realnum)(radius.drad_x_fillfac*findspecieslocal("H-")->den*3.9e-17*
69 (1. - rfield.ContBoltz[hmi.iphmin-1]/ hmi.hmidep));
70 }
71 else
72 {
73 opac.telec += (realnum)(radius.drad_x_fillfac*dense.eden*6.65e-25);
74 opac.thmin += (realnum)(radius.drad_x_fillfac*findspecieslocal("H-")->den*3.9e-17*
75 (1. - rfield.ContBoltz[hmi.iphmin-1]/ hmi.hmidep));
76 }
77
78 /* prevent maser runaway */
79 rt.dTauMase = 0;
80 rt.mas_species = 0;
81 rt.mas_ion = 0;
82 rt.mas_hi = 0;
83 rt.mas_lo = 0;
84
85 /* all lines in iso sequences */
86 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
87 {
88 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
89 {
90 /* this is the parent ion, for HI lines, is 1,
91 * for element He is 1 for He-like (HeI) and 2 for H-like (HeII) */
92 int ion = nelem+1-ipISO;
93 /* do not evaluate in case where trivial parent ion */
94 if( ion <=dense.IonHigh[nelem] && dense.xIonDense[nelem][ion] > dense.density_low_limit )
95 {
96 if( iso_ctrl.lgDielRecom[ipISO] )
97 {
98 // SatelliteLines are indexed by lower level
99 for( ipLo=0; ipLo < iso_sp[ipISO][nelem].numLevels_local; ipLo++ )
100 {
101 RT_line_one_tauinc(SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][ipLo]], ipISO, nelem, -1, ipLo,
102 GetDopplerWidth(dense.AtomicWeight[nelem]) );
103 }
104 }
105
106 for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
107 {
108 for( ipLo=0; ipLo < ipHi; ipLo++ )
109 {
110 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() <= 0 )
111 continue;
112
113 /* actually do the work */
114 RT_line_one_tauinc(iso_sp[ipISO][nelem].trans(ipHi,ipLo), ipISO, nelem, ipHi, ipLo,
115 GetDopplerWidth(dense.AtomicWeight[nelem]) );
116 }
117 }
118 ipLo = 0;
119 /* these are the extra Lyman lines */
120 for( ipHi=iso_sp[ipISO][nelem].st[iso_sp[ipISO][nelem].numLevels_local-1].n()+1; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )
121 {
122 TransitionList::iterator tr = ExtraLymanLines[ipISO][nelem].begin()+ipExtraLymanLines[ipISO][nelem][ipHi];
123 (*tr).Emis().PopOpc() = iso_sp[ipISO][nelem].st[0].Pop();
124
125 /* actually do the work */
126 RT_line_one_tauinc(*tr, -1 ,ipISO, nelem, ipHi,
127 GetDopplerWidth(dense.AtomicWeight[nelem]) );
128 }
129 }
130 }
131 }
132
133 /* increment optical depths for all heavy element lines
134 * same routine does wind and static,
135 * does not start from 0 since first line is dummy */
136 for( i=1; i <= nLevel1; i++ )
137 {
138 RT_line_one_tauinc(TauLines[i], -2, -2, -2, i, GetDopplerWidth(dense.AtomicWeight[(*TauLines[i].Hi()).nelem()-1]) );
139 }
140
141 /* all lines in cooling with g-bar */
142 for( i=0; i < nWindLine; i++ )
143 {
144 /* do not include H-like or He-like in the level two lines since
145 * these are already counted in iso sequences */
146 if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
147 {
148 RT_line_one_tauinc(TauLine2[i], -3, -3, -3, i, GetDopplerWidth(dense.AtomicWeight[(*TauLine2[i].Hi()).nelem()-1]) );
149 }
150 }
151
152 /* the block of inner shell lines */
153 for( i=0; i < nUTA; i++ )
154 {
155 /* populations have not been set */
156 UTALines[i].Emis().PopOpc() = dense.xIonDense[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1];
157 (*UTALines[i].Lo()).Pop() = dense.xIonDense[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1];
158 (*UTALines[i].Hi()).Pop() = 0.;
159 RT_line_one_tauinc(UTALines[i], -4 , -4 , -4 , i, GetDopplerWidth(dense.AtomicWeight[(*UTALines[i].Hi()).nelem()-1]) );
160 }
161
162 /* all hyper fine structure lines */
163 for( i=0; i < nHFLines; i++ )
164 {
165 /* remember current gas-phase abundances */
166 realnum save = dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1];
167
168 /* bail if no abundance */
169 if( save<=0. ) continue;
170
171 /* set gas-phase abundance to total times isotope ratio */
172 dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1] *= hyperfine.HFLabundance[i];
173
174 RT_line_one_tauinc(HFLines[i] , -5 , -5 , -5 , i, GetDopplerWidth(dense.AtomicWeight[(*HFLines[i].Hi()).nelem()-1]) );
175
176 /* put the correct gas-phase abundance back in the array */
177 dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1] = save;
178 }
179
180 /* do large FeII atom if this is enabled */
182
183 /* increment optical depth for the H2 molecule */
184 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
185 (*diatom)->H2_RT_tau_inc();
186
187 /* database Lines*/
188 for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
189 {
190 if( dBaseSpecies[ipSpecies].lgActive )
191 {
192 realnum DopplerWidth = GetDopplerWidth( dBaseSpecies[ipSpecies].fmolweight );
193 for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin();
194 tr != dBaseTrans[ipSpecies].end(); ++tr)
195 {
196 int ipHi = (*tr).ipHi();
197 if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local || (*tr).ipCont() <= 0)
198 continue;
199 int ipLo = (*tr).ipLo();
200
201 RT_line_one_tauinc( *tr, -10, ipSpecies, ipHi, ipLo, DopplerWidth );
202 }
203 }
204 }
205
206 /* following is for static atmosphere */
207 if( wind.lgStatic() )
208 {
209 /* iron fe feii fe2 - overlapping feii lines */
211 }
212
213 if( trace.lgTrace && trace.lgOptcBug )
214 {
215 fprintf( ioQQQ, " RT_tau_inc updated optical depths:\n" );
216 prtmet();
217 }
218
219 if( trace.lgTrace )
220 fprintf( ioQQQ, " RT_tau_inc returns.\n" );
221
222 return;
223}
void FeII_RT_TauInc(void)
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
const int LIMELM
Definition cddefines.h:258
const int NISO
Definition cddefines.h:261
float realnum
Definition cddefines.h:103
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long nWindLine
Definition cdinit.cpp:19
static t_fe2ovr_la & Inst()
Definition cddefines.h:175
TransitionProxy::iterator iterator
Definition transition.h:280
t_conv conv
Definition conv.cpp:5
void CoolEvaluate(double *tot)
Definition cool_eval.cpp:45
t_dense dense
Definition dense.cpp:24
realnum GetDopplerWidth(realnum massAMU)
vector< diatomics * > diatoms
Definition h2.cpp:8
vector< diatomics * >::iterator diatom_iter
Definition h2.h:13
t_hmi hmi
Definition hmi.cpp:5
t_hyperfine hyperfine
Definition hyperfine.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipH_LIKE
Definition iso.h:62
molezone * findspecieslocal(const char buf[])
t_opac opac
Definition opacity.cpp:5
void prtmet(void)
Definition prt_met.cpp:15
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
t_rt rt
Definition rt.cpp:5
void RT_line_all(void)
void RT_line_one_tauinc(const TransitionProxy &t, long int mas_species, long int mas_ion, long int mas_hi, long int mas_lo, realnum DopplerWidth)
void RT_tau_inc(void)
t_save save
Definition save.cpp:5
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
t_thermal thermal
Definition thermal.cpp:5
t_trace trace
Definition trace.cpp:5
Wind wind
Definition wind.cpp:5