cloudy trunk
Loading...
Searching...
No Matches
rt_line_all.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_all do escape and destruction probabilities for all lines in code. */
4#include "cddefines.h"
5#include "taulines.h"
6#include "atomfeii.h"
7#include "dense.h"
8#include "conv.h"
9#include "atoms.h"
10#include "rfield.h"
11#include "wind.h"
12#include "iso.h"
13#include "h2.h"
14#include "opacity.h"
15#include "trace.h"
16#include "lines_service.h"
17#include "atmdat.h"
18#include "hydrogenic.h"
19#include "rt.h"
20#include "cosmology.h"
21#include "physconst.h"
22#include "doppvel.h"
23#include "mole.h"
24
25/*RT_line_all do escape and destruction probabilities for all lines in code. */
26void RT_line_all( void )
27{
28 long int i,
29 ion,
30 ipISO,
31 nelem;
32 long ipHi , ipLo;
33
34 /* this flag says whether to update the line escape probabilities */
35 bool lgPescUpdate = conv.lgFirstSweepThisZone || conv.lgIonStageTrimed;
36
37 DEBUG_ENTRY( "RT_line_all()" );
38
39 if( trace.lgTrace )
40 fprintf( ioQQQ, " RT_line_all called\n" );
41
42 /* rfield.lgDoLineTrans - skip line transfer if requested with no line transfer
43 * conv.nPres2Ioniz is zero during search phase and on first call in this zone */
44 if( !rfield.lgDoLineTrans && conv.nPres2Ioniz )
45 {
46 return;
47 }
48
49 /* this array is huge and takes significant time to zero out or update,
50 * only do so when needed, */
51 if( conv.lgLastSweepThisZone )
52 {
53 /* zero out fine opacity array */
54 memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine*sizeof(realnum) );
55
56 if( 0 && cosmology.lgDo )
57 {
58 realnum dVel = (realnum)SPEEDLIGHT * ( 1.f - 1.f/POW2(1.f+cosmology.redshift_start-cosmology.redshift_current) );
59 rfield.ipFineConVelShift = -(long int)( dVel/ rfield.fine_opac_velocity_width + 0.5 );
60 }
61 else
62 {
63 /* this is offset within fine opacity array caused by Doppler shift
64 * between first zone, the standard of rest, and current position
65 * in case of acceleration away from star in outflowing wind
66 * dWind is positive, this means that the rest frame original
67 * velocity is red shifted to lower energy */
68 realnum dWind = wind.windv - wind.windv0;
69
70 /* will add ipVelShift to index of original energy so redshift has
71 * to be negative */
72 rfield.ipFineConVelShift = -(long int)(dWind / rfield.fine_opac_velocity_width+0.5);
73 }
74 }
75
76# if 0
77 /* this code is a copy of the code within line_one that does cloaking
78 * within this zone. it can be used to see how a particular line
79 * is being treated. */
80 {
81#include "doppvel.h"
82 double dTau , aa;
83 ipISO = 0; nelem = 0;ipLo = 0;
84 ipHi = 23;
85 dTau = iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc() *
86 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().opacity() /
87 GetDopplerWidth(dense.AtomicWeight[nelem]) + opac.opacity_abs[iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont()-1];
88 dTau *= radius.drad_x_fillfac;
89 aa = log(1. + dTau ) / SDIV(dTau);
90 fprintf(ioQQQ,"DEBUG dTau\t%.2f\t%.5e\t%.5e\t%.5e\n",fnzone,dTau,
91 radius.drad_x_fillfac,
92 aa);
93 }
94# endif
95
96 /* find Stark escape probabilities for hydrogen itself */
97 if( lgPescUpdate )
98 RT_stark();
99
100 /*if( lgUpdateFineOpac )
101 fprintf(ioQQQ,"debuggg\tlgUpdateFineOpac in rt_line_all\n");*/
102 for( ipISO=ipH_LIKE; ipISO < NISO; ++ipISO )
103 {
104 /* loop over all iso-electronic sequences */
105 for( nelem=ipISO; nelem < LIMELM; ++nelem )
106 {
107 /* parent ion stage, for H is 1, for He is 1 for He-like and
108 * 2 for H-like */
109 ion = nelem+1-ipISO;
110
111 /* element turned off */
112 if( !dense.lgElmtOn[nelem] )
113 continue;
114 /* need we consider this ion? */
115 if( ion <= dense.IonHigh[nelem] )
116 {
117 /* loop over all lines */
118 for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ++ipHi )
119 {
120 for( ipLo=0; ipLo < ipHi; ++ipLo )
121 {
122 /* negative ipCont means this is not a real line, so do not
123 * transfer it */
124 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() < 0 )
125 continue;
126
127 /* generate escape prob, pumping rate, destruction prob,
128 * inward outward fracs */
129 fixit(); // should this use pestrk_up or pestrk?
130 RT_line_one( iso_sp[ipISO][nelem].trans(ipHi,ipLo),
131 true,(realnum)iso_sp[ipISO][nelem].ex[ipHi][ipLo].pestrk_up,
132 GetDopplerWidth(dense.AtomicWeight[nelem]));
133
134 /* set true to print pump rates*/
135 enum {DEBUG_LOC=false};
136 if( DEBUG_LOC && nelem==1&& ipLo==0 /*&& iteration==2*/ )
137 {
138 fprintf(ioQQQ,"DEBUG pdest\t%3li\t%.2f\t%.3e\n",
139 ipHi ,
140 fnzone,
141 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pdest());
142 }
143 }
144 }
145
146 /* this is option to not do destruction probabilities
147 * for case b no pdest option */
148 if( opac.lgCaseB_no_pdest )
149 {
150 ipLo = 0;
151 /* okay to let this go to numLevels_max. */
152 for( ipHi=ipLo+1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ++ipHi )
153 {
154 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() <= 0 )
155 continue;
156
157 /* hose the previously computed destruction probability */
158 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pdest() = SMALLFLOAT;
159 }
160 }
161
162 ipLo = 0;
163 /* these are the extra Lyman lines for the iso sequences */
164 /*for( ipHi=2; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )*/
165 /* only update if significant abundance and need to update fine opac */
166 if( dense.xIonDense[nelem][ion] > 1e-30 && ( conv.lgFirstSweepThisZone || conv.lgLastSweepThisZone ) )
167 {
168 for( ipHi=iso_sp[ipISO][nelem].st[iso_sp[ipISO][nelem].numLevels_local-1].n()+1; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )
169 {
170 TransitionList::iterator tr = ExtraLymanLines[ipISO][nelem].begin()+ipExtraLymanLines[ipISO][nelem][ipHi];
171 /* we just want the population of the ground state */
172 (*tr).Emis().PopOpc() = iso_sp[ipISO][nelem].st[0].Pop();
173 (*(*tr).Lo()).Pop() =
174 iso_sp[ipISO][nelem].st[ipLo].Pop();
175
176 /* actually do the work */
177 RT_line_one( *tr, true, 0.f,
178 GetDopplerWidth(dense.AtomicWeight[nelem]));
179 }
180 }
181
182 }/* if nelem if ion <=dense.IonHigh */
183 }/* loop over nelem */
184 }/* loop over ipISO */
185
186 /* note that pesc and dest are updated no matter what escprob logic we
187 * specify, and that these are not updated if we have overrun the
188 * optical depth scale. only update here in that case */
189 if( lgTauGood( iso_sp[ipH_LIKE][ipHYDROGEN].trans(iso_ctrl.nLyaLevel[ipH_LIKE],0) ) )
190 {
191 /* add on destruction of hydrogen Lya by FeII
192 * now add in FeII deexcitation via overlap,
193 * but only as loss, not photoionization, source
194 * dstfe2lya is Ly-alpha deexcit by FeII overlap - conversion into FeII em */
195 /* find FeII overlap destruction rate,
196 * this does NOTHING when large FeII atom is turned on,
197 * in this case evaluation is done in call to FeIILevelPops */
199 /*fprintf(ioQQQ,"DEBUG fe2 %.2e %.2e\n", hydro.dstfe2lya ,
200 hydro.HLineWidth);*/
201
202 /* >>chng 00 jan 06, let dest be large than 1 to desaturate the atom */
203 /* >>chng 01 apr 01, add test for tout >= 0.,
204 * must not add to Pdest when it has not been refreshed here */
205 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pdest() += hydro.dstfe2lya;
206 }
207
208 /* is continuum pumping of H Lyman lines included? yes, but turned off
209 * with atom h-like Lyman pumping off command */
210 if( !hydro.lgLymanPumping )
211 {
212 ipISO = ipH_LIKE;
213 nelem = ipHYDROGEN;
214 ipLo = 0;
215 for( ipHi=ipLo+1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ++ipHi )
216 {
217 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().pump() = 0.;
218 }
219 }
220
221 {
222 /* following should be set true to print ots contributors for he-like lines*/
223 enum {DEBUG_LOC=false};
224 if( DEBUG_LOC && nzone>433 /*&& iteration==2*/ )
225 {
226 /* option to dump a line */
228 }
229 }
230
231 /* level 1 lines */
232 for( i=1; i <= nLevel1; i++ )
233 {
234 RT_line_one( TauLines[i], true,0.f, GetDopplerWidth(dense.AtomicWeight[(*TauLines[i].Hi()).nelem()-1]) );
235 }
236
237 {
238 // 09 mar 28, this transition, between the highest two levels in the
239 // very simple version of the model Fe II atom, is troublesome due
240 // to Lya pumping, which can invert the highest two level populations.
241 // The electron scattering escape probability shows noise as a result.
242 // The sim nova_photos will throw an assert if this is removed
243 // This model is unphysical since highest two levels are close together
244 // and Lya is populating the highest one. The problem does not occur
245 // when the large FeII atom is turned on - that has realistic level
246 // energies.
247 // evaluate electron scattering escape probability one time per zone
248 static realnum P_elec_esc_ipTFe56;
249 static long nSave;
250 if( nzone <= 1 || nzone!=nSave )
251 {
252 P_elec_esc_ipTFe56 = TauLines[ipTFe56].Emis().Pelec_esc();
253 nSave = nzone;
254 }
255 else
256 TauLines[ipTFe56].Emis().Pelec_esc() = P_elec_esc_ipTFe56;
257 }
258
259 /* external database lines */
260 for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
261 {
262 if( dBaseSpecies[ipSpecies].lgActive )
263 {
264 realnum DopplerWidth = GetDopplerWidth( dBaseSpecies[ipSpecies].fmolweight );
265 for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin();
266 tr != dBaseTrans[ipSpecies].end(); ++tr)
267 {
268 int ipHi = (*tr).ipHi();
269 if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local)
270 continue;
271 if( (*tr).ipCont() > 0 )
272 {
273 RT_line_one( *tr, true, 0.f, DopplerWidth );
274 }
275 }
276 }
277 }
278
279 /* this is a major time sink for this routine - only evaluate on last
280 * sweep when fine opacities are updated since only effect of UTAs is
281 * to pump inner shell lines and add to total opacity */
282 if( conv.lgFirstSweepThisZone || conv.lgLastSweepThisZone )
283 {
284 for( i=0; i < nUTA; i++ )
285 {
286 /* these are not defined in cooling routines so must be set here */
287 UTALines[i].Emis().PopOpc() = dense.xIonDense[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1];
288 (*UTALines[i].Lo()).Pop() = dense.xIonDense[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1];
289 (*UTALines[i].Hi()).Pop() = 0.;
290 RT_line_one( UTALines[i], true,0.f,
291 GetDopplerWidth(dense.AtomicWeight[(*UTALines[i].Hi()).nelem()-1]) );
292 }
293 }
294
295 /* do satellite lines for iso sequences gt H-like
296 * H-like ions only have one electron, no satellite lines. */
297 for( ipISO=ipHE_LIKE; ipISO<NISO; ++ipISO )
298 {
299 /* loop over all iso-electronic sequences */
300 for( nelem=ipISO; nelem<LIMELM; ++nelem )
301 {
302 if( dense.lgElmtOn[nelem] && iso_ctrl.lgDielRecom[ipISO] )
303 {
304 for( ipLo=0; ipLo < iso_sp[ipISO][nelem].numLevels_local; ++ipLo )
305 {
306 RT_line_one( SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][ipLo]], true,0.f,
307 GetDopplerWidth(dense.AtomicWeight[nelem]) );
308 }
309 }
310 }
311 }
312
313 /* the large H2 molecule */
314 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
315 (*diatom)->H2_RTMake( );
316
317 return;
318}
long ipTFe56
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
double fnzone
Definition cddefines.cpp:15
const int LIMELM
Definition cddefines.h:258
#define POW2
Definition cddefines.h:929
const int NISO
Definition cddefines.h:261
float realnum
Definition cddefines.h:103
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void fixit(void)
Definition service.cpp:991
static t_fe2ovr_la & Inst()
Definition cddefines.h:175
TransitionProxy::iterator iterator
Definition transition.h:280
void atoms_fe2ovr(void)
t_conv conv
Definition conv.cpp:5
t_cosmology cosmology
Definition cosmology.cpp:11
const realnum SMALLFLOAT
Definition cpu.h:191
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_hydro hydro
Definition hydrogenic.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipH1s
Definition iso.h:27
const int ipHE_LIKE
Definition iso.h:63
const int ipH2p
Definition iso.h:29
const int ipH_LIKE
Definition iso.h:62
t_opac opac
Definition opacity.cpp:5
UNUSED const double SPEEDLIGHT
Definition physconst.h:100
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
void RT_line_one(const TransitionProxy &t, bool lgShield_this_zone, realnum pestrk, realnum DopplerWidth)
void RT_stark(void)
Definition rt_stark.cpp:14
void RT_line_all(void)
static double * ex
Definition species2.cpp:28
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)
vector< TransitionList > dBaseTrans
Definition taulines.cpp:17
multi_arr< int, 3 > ipSatelliteLines
Definition taulines.cpp:37
long int nUTA
Definition taulines.cpp:26
long int nLevel1
Definition taulines.cpp:28
multi_arr< int, 3 > ipExtraLymanLines
Definition taulines.cpp:24
TransitionList TauLines("TauLines", &AnonStates)
species * dBaseSpecies
Definition taulines.cpp:14
t_trace trace
Definition trace.cpp:5
void DumpLine(const TransitionProxy &t)
bool lgTauGood(const TransitionProxy &t)
Definition transition.h:570
Wind wind
Definition wind.cpp:5