cloudy trunk
Loading...
Searching...
No Matches
hydrolevel.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/*HydroLevel calls iso_level to solve for ionization balance
4 * and level populations of model hydrogen atom */
5/*PrtHydroTrace1 print trace info for hydrogen-like species */
6#include "cddefines.h"
7#include "taulines.h"
8#include "iso.h"
9#include "dense.h"
10#include "secondaries.h"
11#include "trace.h"
12#include "phycon.h"
13#include "ionbal.h"
14#include "hydrogenic.h"
15
16/*PrtHydroTrace1a print trace info for hydrogen-like species */
17STATIC void PrtHydroTrace1a(long nelem )
18{
19 double colfrc,
20 phtfrc,
21 secfrc;
22
23 DEBUG_ENTRY( "PrtHydroTrace1a()" );
24
25 /* identify how atom is ionized for full trace */
26 if( iso_sp[ipH_LIKE][nelem].xIonSimple > 0. )
27 {
28 /* fraction of ionization due to photoionization */
29 phtfrc = iso_sp[ipH_LIKE][nelem].fb[ipH1s].gamnc/((dense.eden*(iso_sp[ipH_LIKE][nelem].RadRec_effec +
30 ionbal.CotaRate[nelem]) )*
31 iso_sp[ipH_LIKE][nelem].xIonSimple);
32
33 /* fraction of ionization due to collisional ionization */
34 colfrc = (iso_sp[ipH_LIKE][nelem].fb[ipH1s].ColIoniz*dense.EdenHCorr )/
35 ((dense.eden*(iso_sp[ipH_LIKE][nelem].RadRec_effec +
36 ionbal.CotaRate[0]) )*
37 iso_sp[ipH_LIKE][nelem].xIonSimple);
38
39 /* fraction of ionization due to secondary ionization */
40 secfrc = secondaries.csupra[nelem][nelem]/((dense.eden*(iso_sp[ipH_LIKE][nelem].RadRec_effec +
41 ionbal.CotaRate[0]) )*
42 iso_sp[ipH_LIKE][nelem].xIonSimple);
43 }
44 else
45 {
46 phtfrc = 0.;
47 colfrc = 0.;
48 secfrc = 0.;
49 }
50
51 fprintf( ioQQQ, " HydroLevel Z=%2ld called, simple II/I=",nelem);
52 PrintE93( ioQQQ, iso_sp[ipH_LIKE][nelem].xIonSimple);
53 fprintf( ioQQQ," PhotFrc:");
54 PrintE93( ioQQQ,phtfrc);
55 fprintf(ioQQQ," ColFrc:");
56 PrintE93( ioQQQ,colfrc);
57 fprintf( ioQQQ," SecFrc");
58 PrintE93( ioQQQ, secfrc);
59 fprintf( ioQQQ," Te:");
61 fprintf( ioQQQ," eden:");
62 PrintE93( ioQQQ,dense.eden);
63 fprintf( ioQQQ,"\n");
64 return;
65}
66
67/*PrtHydroTrace1 print trace info for hydrogen-like species */
68STATIC void PrtHydroTrace1(long nelem )
69{
70 long int ipHi , ipLo , i;
71
72 DEBUG_ENTRY( "PrtHydroTrace1()" );
73
74 fprintf( ioQQQ,
75 " HydroLevel%3ld finds arrays, with optical depths defined? %li induced 2ph=%12.3e\n",
76 nelem, iteration, iso_sp[ipH_LIKE][nelem].TwoNu[0].induc_up );
77 /* 06 aug 28, from numLevels_max to _local. */
78 for( ipHi=ipH2s; ipHi < iso_sp[ipH_LIKE][nelem].numLevels_local; ipHi++ )
79 {
80 fprintf( ioQQQ, "up:%2ld", ipHi );
81 fprintf( ioQQQ, "lo" );
82 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
83 {
84 fprintf( ioQQQ, "%9ld", ipLo );
85 }
86 fprintf( ioQQQ, "\n" );
87
88 fprintf( ioQQQ, "%3ld", ipHi );
89 fprintf( ioQQQ, " A*esc" );
90 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
91 {
92 fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul()*
93 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pesc() ));
94 }
95 fprintf( ioQQQ, "\n" );
96
97 fprintf( ioQQQ, "%3ld", ipHi );
98 fprintf( ioQQQ, " A*ees" );
99 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
100 {
101 fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul()*
102 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pelec_esc() ));
103 }
104 fprintf( ioQQQ, "\n" );
105
106 fprintf( ioQQQ, "%3ld", ipHi );
107 fprintf( ioQQQ, " tauin" );
108 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
109 {
110 fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().TauIn() ));
111 }
112 fprintf( ioQQQ, "\n" );
113
114 fprintf( ioQQQ, "%3ld", ipHi );
115 fprintf( ioQQQ, " t tot" );
116 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
117 {
118 fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().TauTot() ));
119 }
120 fprintf( ioQQQ, "\n" );
121
122 fprintf( ioQQQ, "%3ld", ipHi );
123 fprintf( ioQQQ, " Esc " );
124 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
125 {
126 fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pesc() ));
127 }
128 fprintf( ioQQQ, "\n" );
129
130 fprintf( ioQQQ, "%3ld", ipHi );
131 fprintf( ioQQQ, " Eesc " );
132 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
133 {
134 fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pelec_esc() ));
135 }
136 fprintf( ioQQQ, "\n" );
137
138 fprintf( ioQQQ, "%3ld", ipHi );
139 fprintf( ioQQQ, " Dest " );
140 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
141 {
142 fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pdest()) );
143 }
144 fprintf( ioQQQ, "\n" );
145
146 fprintf( ioQQQ, "%3ld", ipHi );
147 fprintf( ioQQQ, " A*dst" );
148 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
149 {
150 fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul()*
151 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pdest() ));
152 }
153 fprintf( ioQQQ, "\n" );
154
155 fprintf( ioQQQ, "%3ld", ipHi );
156 fprintf( ioQQQ, " StrkE" );
157 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
158 {
159 fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].ex[ipHi][ipLo].pestrk_up ));
160 }
161 fprintf( ioQQQ, "\n" );
162
163 fprintf( ioQQQ, "%3ld", ipHi );
164 fprintf( ioQQQ, " B(ul)" );
165 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
166 {
167 fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().pump()*
168 iso_sp[ipH_LIKE][nelem].st[ipLo].g()/iso_sp[ipH_LIKE][nelem].st[ipHi].g() ));
169 }
170 fprintf( ioQQQ, "\n" );
171
172 fprintf( ioQQQ, "%3ld", ipHi );
173 fprintf( ioQQQ, " tcont" );
174 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
175 {
176 fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().TauCon() ));
177 }
178 fprintf( ioQQQ, "\n" );
179
180 fprintf( ioQQQ, "%3ld", ipHi );
181 fprintf( ioQQQ, " C(ul)" );
182 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
183 {
184 fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Coll().ColUL( colliders ) ));
185 }
186 fprintf( ioQQQ, "\n" );
187
188 if( ipHi == 2 )
189 {
190 fprintf( ioQQQ, " FeIIo");
191 fprintf( ioQQQ,PrintEfmt("%9.2e",
192 hydro.dstfe2lya* iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Aul() ));
193 fprintf( ioQQQ, "\n");
194 }
195 }
196
197 fprintf( ioQQQ, " " );
198 /* 06 aug 28, from numLevels_max to _local. */
199 for( i=1; i < iso_sp[ipH_LIKE][nelem].numLevels_local; i++ )
200 {
201 fprintf( ioQQQ, "%9ld", i );
202 }
203 fprintf( ioQQQ, "\n" );
204 return;
205}
206
207/*HydroLevel calls iso_level to solve for ionization balance
208 * and level populations of model hydrogen atom */
209void HydroLevel(long int nelem)
210{
211 long int i;
212 int ipISO = ipH_LIKE;
213
214 DEBUG_ENTRY( "HydroLevel()" );
215
216 /* check that we were called with valid charge */
217 ASSERT( nelem >= 0);
218 ASSERT( nelem < LIMELM );
219
220 /* option to print some rates */
221 if( (trace.lgTrace && trace.lgIsoTraceFull[ipISO]) && (nelem == trace.ipIsoTrace[ipISO]) )
222 PrtHydroTrace1(nelem);
223
224 if( trace.lgHBug && trace.lgTrace )
225 PrtHydroTrace1a(nelem);
226
227 /* this is main trace h-like printout */
228 if( (trace.lgIsoTraceFull[ipISO] && trace.lgTrace) && (nelem == trace.ipIsoTrace[ipISO]) )
229 {
230 fprintf( ioQQQ, " HLEV HGAMNC" );
231 PrintE93( ioQQQ, iso_sp[ipISO][nelem].fb[ipH1s].gamnc );
232 /* 06 aug 28, from numLevels_max to _local. */
233 for( i=ipH2s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
234 {
235 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].fb[i].gamnc ));
236 }
237 fprintf( ioQQQ, "\n" );
238
239 fprintf( ioQQQ, " HLEV TOTCAP" );
240 /* 06 aug 28, from numLevels_max to _local. */
241 for( i=1; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
242 {
243 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].fb[i].RateCont2Level/dense.eden ));
244 }
245 fprintf( ioQQQ," tot");
246 fprintf( ioQQQ,PrintEfmt("%10.2e", ionbal.RateRecomTot[nelem][nelem-ipISO]/dense.eden ) );
247 fprintf( ioQQQ, "\n" );
248
249 fprintf( ioQQQ, " HLEV IND Rc" );
250 /* 06 aug 28, from numLevels_max to _local. */
251 for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
252 {
253 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].fb[i].RecomInducRate*iso_sp[ipISO][nelem].fb[i].PopLTE ));
254 }
255 fprintf( ioQQQ, "\n" );
256
257 /* incuded recombination rate coefficients */
258 fprintf( ioQQQ, " IND Rc LTE " );
259 /* 06 aug 28, from numLevels_max to _local. */
260 for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
261 {
262 fprintf(ioQQQ,PrintEfmt("%9.2e",
263 iso_sp[ipISO][nelem].fb[i].gamnc*iso_sp[ipISO][nelem].fb[i].PopLTE ));
264 }
265 fprintf( ioQQQ, "\n" );
266
267 /* LTE level populations */
268 fprintf( ioQQQ, " HLEV HLTE" );
269 /* 06 aug 28, from numLevels_max to _local. */
270 for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
271 {
272 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].fb[i].PopLTE ));
273 }
274 fprintf( ioQQQ, "\n" );
275
276 /* fraction of total ionization due to collisions from given level */
277 fprintf( ioQQQ, " HLEVfr cion" );
278 /* 06 aug 28, from numLevels_max to _local. */
279 for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
280 {
281 fprintf(ioQQQ,PrintEfmt("%9.2e",
282 iso_sp[ipISO][nelem].fb[i].ColIoniz*dense.EdenHCorr/MAX2(1e-30,iso_sp[ipISO][nelem].fb[i].RateLevel2Cont) ) );
283 }
284 fprintf( ioQQQ, "\n" );
285
286 /* fraction of total ionization due to photoionization from given level */
287 if( ionbal.RateRecomTot[nelem][nelem]> 0. )
288 {
289 fprintf( ioQQQ, " HLEVfrPhIon" );
290 /* 06 aug 28, from numLevels_max to _local. */
291 for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
292 {
293 fprintf(ioQQQ,PrintEfmt("%9.2e",
294 iso_sp[ipISO][nelem].fb[i].gamnc/MAX2(1e-30,iso_sp[ipISO][nelem].fb[i].RateLevel2Cont) ) );
295 }
296 fprintf( ioQQQ, "\n" );
297 }
298
299 fprintf( ioQQQ, " HLEV HN" );
300 /* 06 aug 28, from numLevels_max to _local. */
301 for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
302 {
303 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].st[i].Pop() ));
304 }
305 fprintf( ioQQQ, "\n" );
306
307 fprintf( ioQQQ, " HLEV b(n)" );
308 /* 06 aug 28, from numLevels_max to _local. */
309 for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
310 {
311 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].fb[i].DepartCoef ));
312 }
313 fprintf( ioQQQ, "\n" );
314
315 fprintf( ioQQQ, " HLEV X12tot");
316 fprintf(ioQQQ,PrintEfmt("%9.2e", secondaries.x12tot ));
317 fprintf( ioQQQ," Grn dest:");
318 fprintf(ioQQQ,PrintEfmt("%9.2e",
319 ionbal.RateIoniz[nelem][nelem][nelem+1] ));
320 fprintf(ioQQQ, "\n");
321 }
322
323 if( trace.lgTrace )
324 {
325 /* iso.RecomTotal[nelem] is gross rec coef, computed here while filling in matrix
326 * elements, all physical processes included.
327 * RadRec_effec is total effective radiative only */
328 fprintf( ioQQQ, " HydroLevel Z:%2ld return %s te=",
329 nelem,
330 iso_sp[ipISO][nelem].chTypeAtomUsed );
331 PrintE93( ioQQQ,phycon.te);
332 fprintf( ioQQQ," density=%.4e", dense.xIonDense[nelem][nelem-ipISO] );
333
334 fprintf( ioQQQ," simple=%.4e",iso_sp[ipISO][nelem].xIonSimple);
335
336 fprintf( ioQQQ," b1=%.2e",iso_sp[ipISO][nelem].fb[ipH1s].DepartCoef);
337
338 fprintf( ioQQQ," ion rate=%.4e",ionbal.RateIonizTot(nelem,nelem-ipISO) );
339
340 fprintf( ioQQQ," TotRec=%.4e",ionbal.RateRecomTot[nelem][nelem-ipISO]);
341
342 fprintf( ioQQQ," RadRec=%.4e",iso_sp[ipISO][nelem].RadRec_effec);
343 fprintf( ioQQQ, "\n");
344 }
345 return;
346}
FILE * ioQQQ
Definition cddefines.cpp:7
long int iteration
Definition cddefines.cpp:16
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define PrintEfmt(F, V)
Definition cddefines.h:1472
void PrintE93(FILE *, double)
Definition service.cpp:838
const int LIMELM
Definition cddefines.h:258
#define MAX2
Definition cddefines.h:782
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
ColliderList colliders
Definition collision.cpp:57
t_dense dense
Definition dense.cpp:24
t_hydro hydro
Definition hydrogenic.cpp:5
void HydroLevel(long int nelem)
STATIC void PrtHydroTrace1(long nelem)
STATIC void PrtHydroTrace1a(long nelem)
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH1s
Definition iso.h:27
const int ipH2p
Definition iso.h:29
const int ipH2s
Definition iso.h:28
const int ipH_LIKE
Definition iso.h:62
t_phycon phycon
Definition phycon.cpp:6
t_secondaries secondaries
static double * ex
Definition species2.cpp:28
static double * g
Definition species2.cpp:28
t_trace trace
Definition trace.cpp:5