cloudy  trunk
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 */
17 STATIC 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:");
60  PrintE93( ioQQQ,phycon.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 */
68 STATIC 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 */
209 void 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 }
t_trace::lgIsoTraceFull
bool lgIsoTraceFull[NISO]
Definition: trace.h:88
t_dense::eden
double eden
Definition: dense.h:190
t_dense::EdenHCorr
double EdenHCorr
Definition: dense.h:216
colliders
ColliderList colliders
Definition: collision.cpp:57
dense
t_dense dense
Definition: dense.cpp:24
secondaries.h
t_ionbal::CotaRate
realnum CotaRate[LIMELM]
Definition: ionbal.h:242
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
PrintE93
void PrintE93(FILE *, double)
Definition: service.cpp:838
t_hydro::dstfe2lya
realnum dstfe2lya
Definition: hydrogenic.h:60
STATIC
#define STATIC
Definition: cddefines.h:97
t_ionbal::RateRecomTot
double ** RateRecomTot
Definition: ionbal.h:198
t_iso_sp::RadRec_effec
double RadRec_effec
Definition: iso.h:517
t_iso_sp::numLevels_local
long int numLevels_local
Definition: iso.h:498
phycon
t_phycon phycon
Definition: phycon.cpp:6
trace.h
ex
static double * ex
Definition: species2.cpp:28
t_iso_sp::xIonSimple
double xIonSimple
Definition: iso.h:465
t_secondaries::x12tot
realnum x12tot
Definition: secondaries.h:53
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
iso.h
PrtHydroTrace1
STATIC void PrtHydroTrace1(long nelem)
Definition: hydrolevel.cpp:68
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:408
dense.h
trace
t_trace trace
Definition: trace.cpp:5
cddefines.h
ipH2s
const int ipH2s
Definition: iso.h:28
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:452
MAX2
#define MAX2
Definition: cddefines.h:782
ionbal
t_ionbal ionbal
Definition: ionbal.cpp:5
LIMELM
const int LIMELM
Definition: cddefines.h:258
PrtHydroTrace1a
STATIC void PrtHydroTrace1a(long nelem)
Definition: hydrolevel.cpp:17
HydroLevel
void HydroLevel(long int nelem)
Definition: hydrolevel.cpp:209
hydro
t_hydro hydro
Definition: hydrogenic.cpp:5
iteration
long int iteration
Definition: cddefines.cpp:16
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
hydrogenic.h
ipH2p
const int ipH2p
Definition: iso.h:29
t_secondaries::csupra
realnum ** csupra
Definition: secondaries.h:21
ionbal.h
secondaries
t_secondaries secondaries
Definition: secondaries.cpp:5
t_trace::ipIsoTrace
long int ipIsoTrace[NISO]
Definition: trace.h:91
t_ionbal::RateIonizTot
double RateIonizTot(long nelem, long ion)
Definition: ionbal.h:254
PrintEfmt
#define PrintEfmt(F, V)
Definition: cddefines.h:1472
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
taulines.h
phycon.h
t_trace::lgHBug
bool lgHBug
Definition: trace.h:85
ipH1s
const int ipH1s
Definition: iso.h:27
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:613
t_phycon::te
double te
Definition: phycon.h:11
t_ionbal::RateIoniz
double *** RateIoniz
Definition: ionbal.h:184
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
g
static double * g
Definition: species2.cpp:28
Generated by   doxygen 1.8.17