cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
hydrolevel.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 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 "iso.h"
8 #include "secondaries.h"
9 #include "trace.h"
10 #include "phycon.h"
11 #include "ionbal.h"
12 #include "hydrogenic.h"
13 #include "freebound.h"
14 #include "two_photon.h"
15 #include "dense.h"
16 
17 /*PrtHydroTrace1a print trace info for hydrogen-like species */
18 STATIC void PrtHydroTrace1a( long ipISO, long nelem )
19 {
20  double colfrc,
21  phtfrc,
22  secfrc;
23 
24  DEBUG_ENTRY( "PrtHydroTrace1a()" );
25 
26  /* identify how atom is ionized for full trace */
27  if( iso_sp[ipISO][nelem].xIonSimple > 0. )
28  {
29  /* fraction of ionization due to photoionization */
30  phtfrc = iso_sp[ipISO][nelem].fb[ipH1s].gamnc/((dense.eden*(iso_sp[ipISO][nelem].RadRec_effec +
31  ionbal.CotaRate[nelem]) )*
32  iso_sp[ipISO][nelem].xIonSimple);
33 
34  /* fraction of ionization due to collisional ionization */
35  colfrc = (iso_sp[ipISO][nelem].fb[ipH1s].ColIoniz*dense.EdenHCorr )/
36  ((dense.eden*(iso_sp[ipISO][nelem].RadRec_effec +
37  ionbal.CotaRate[0]) )*
38  iso_sp[ipISO][nelem].xIonSimple);
39 
40  /* fraction of ionization due to secondary ionization */
41  secfrc = secondaries.csupra[nelem][nelem]/((dense.eden*(iso_sp[ipISO][nelem].RadRec_effec +
42  ionbal.CotaRate[0]) )*
43  iso_sp[ipISO][nelem].xIonSimple);
44  }
45  else
46  {
47  phtfrc = 0.;
48  colfrc = 0.;
49  secfrc = 0.;
50  }
51 
52  fprintf( ioQQQ, " HydroLevel Z=%2ld called, simple II/I=",nelem);
53  PrintE93( ioQQQ, iso_sp[ipISO][nelem].xIonSimple);
54  fprintf( ioQQQ," PhotFrc:");
55  PrintE93( ioQQQ,phtfrc);
56  fprintf(ioQQQ," ColFrc:");
57  PrintE93( ioQQQ,colfrc);
58  fprintf( ioQQQ," SecFrc");
59  PrintE93( ioQQQ, secfrc);
60  fprintf( ioQQQ," Te:");
61  PrintE93( ioQQQ,phycon.te);
62  fprintf( ioQQQ," eden:");
63  PrintE93( ioQQQ,dense.eden);
64  fprintf( ioQQQ,"\n");
65  return;
66 }
67 
68 /*PrtHydroTrace1 print trace info for hydrogen-like species */
69 STATIC void PrtHydroTrace1( long ipISO, long nelem )
70 {
71  long int ipHi , ipLo , i;
72 
73  DEBUG_ENTRY( "PrtHydroTrace1()" );
74 
75  fprintf( ioQQQ,
76  " HydroLevel%3ld finds arrays, with optical depths defined? %li induced 2ph=%12.3e\n",
77  nelem, iteration, iso_sp[ipISO][nelem].TwoNu[0].induc_up );
78  /* 06 aug 28, from numLevels_max to _local. */
79  for( ipHi=ipH2s; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
80  {
81  fprintf( ioQQQ, "up:%2ld", ipHi );
82  fprintf( ioQQQ, "lo" );
83  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
84  {
85  fprintf( ioQQQ, "%9ld", ipLo );
86  }
87  fprintf( ioQQQ, "\n" );
88 
89  fprintf( ioQQQ, "%3ld", ipHi );
90  fprintf( ioQQQ, " A*esc" );
91  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
92  {
93  fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul()*
94  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pesc() ));
95  }
96  fprintf( ioQQQ, "\n" );
97 
98  fprintf( ioQQQ, "%3ld", ipHi );
99  fprintf( ioQQQ, " A*ees" );
100  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
101  {
102  fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul()*
103  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pelec_esc() ));
104  }
105  fprintf( ioQQQ, "\n" );
106 
107  fprintf( ioQQQ, "%3ld", ipHi );
108  fprintf( ioQQQ, " tauin" );
109  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
110  {
111  fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn() ));
112  }
113  fprintf( ioQQQ, "\n" );
114 
115  fprintf( ioQQQ, "%3ld", ipHi );
116  fprintf( ioQQQ, " t tot" );
117  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
118  {
119  fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauTot() ));
120  }
121  fprintf( ioQQQ, "\n" );
122 
123  fprintf( ioQQQ, "%3ld", ipHi );
124  fprintf( ioQQQ, " Esc " );
125  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
126  {
127  fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pesc() ));
128  }
129  fprintf( ioQQQ, "\n" );
130 
131  fprintf( ioQQQ, "%3ld", ipHi );
132  fprintf( ioQQQ, " Eesc " );
133  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
134  {
135  fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pelec_esc() ));
136  }
137  fprintf( ioQQQ, "\n" );
138 
139  fprintf( ioQQQ, "%3ld", ipHi );
140  fprintf( ioQQQ, " Dest " );
141  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
142  {
143  fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pdest()) );
144  }
145  fprintf( ioQQQ, "\n" );
146 
147  fprintf( ioQQQ, "%3ld", ipHi );
148  fprintf( ioQQQ, " A*dst" );
149  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
150  {
151  fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul()*
152  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pdest() ));
153  }
154  fprintf( ioQQQ, "\n" );
155 
156  fprintf( ioQQQ, "%3ld", ipHi );
157  fprintf( ioQQQ, " StrkE" );
158  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
159  {
160  fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].ex[ipHi][ipLo].pestrk_up ));
161  }
162  fprintf( ioQQQ, "\n" );
163 
164  fprintf( ioQQQ, "%3ld", ipHi );
165  fprintf( ioQQQ, " B(ul)" );
166  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
167  {
168  fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().pump()*
169  iso_sp[ipISO][nelem].st[ipLo].g()/iso_sp[ipISO][nelem].st[ipHi].g() ));
170  }
171  fprintf( ioQQQ, "\n" );
172 
173  fprintf( ioQQQ, "%3ld", ipHi );
174  fprintf( ioQQQ, " tcont" );
175  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
176  {
177  fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauCon() ));
178  }
179  fprintf( ioQQQ, "\n" );
180 
181  fprintf( ioQQQ, "%3ld", ipHi );
182  fprintf( ioQQQ, " C(ul)" );
183  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
184  {
185  fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].trans(ipHi,ipLo).Coll().ColUL( colliders ) ));
186  }
187  fprintf( ioQQQ, "\n" );
188 
189  if( ipHi == 2 && ipISO==ipH_LIKE && nelem==ipHYDROGEN )
190  {
191  fprintf( ioQQQ, " FeIIo");
192  fprintf( ioQQQ,PrintEfmt("%9.2e",
193  hydro.dstfe2lya* iso_sp[ipISO][nelem].trans(ipH2p,ipH1s).Emis().Aul() ));
194  fprintf( ioQQQ, "\n");
195  }
196  }
197 
198  fprintf( ioQQQ, " " );
199  /* 06 aug 28, from numLevels_max to _local. */
200  for( i=1; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
201  {
202  fprintf( ioQQQ, "%9ld", i );
203  }
204  fprintf( ioQQQ, "\n" );
205  return;
206 }
207 
208 /*HydroLevel calls iso_level to solve for ionization balance
209  * and level populations of model hydrogen atom */
210 void HydroLevel(long ipISO, long nelem)
211 {
212  long int i;
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(ipISO, nelem);
223 
224  if( trace.lgHBug && trace.lgTrace )
225  PrtHydroTrace1a(ipISO, 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].st[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].st[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_secondaries::x12tot
realnum x12tot
Definition: secondaries.h:65
PrintE93
void PrintE93(FILE *, double)
Definition: service.cpp:1019
t_dense::EdenHCorr
double EdenHCorr
Definition: dense.h:227
two_photon.h
t_trace::lgIsoTraceFull
bool lgIsoTraceFull[NISO]
Definition: trace.h:85
PrintEfmt
#define PrintEfmt(F, V)
Definition: cddefines.h:1364
phycon
t_phycon phycon
Definition: phycon.cpp:6
t_trace::ipIsoTrace
long int ipIsoTrace[NISO]
Definition: trace.h:88
PrtHydroTrace1a
STATIC void PrtHydroTrace1a(long ipISO, long nelem)
Definition: hydrolevel.cpp:18
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:481
dense
t_dense dense
Definition: global.cpp:15
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
t_ionbal::RateRecomTot
double ** RateRecomTot
Definition: ionbal.h:191
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
iteration
long int iteration
Definition: cddefines.cpp:16
trace
t_trace trace
Definition: trace.cpp:5
ionbal
t_ionbal ionbal
Definition: ionbal.cpp:8
colliders
ColliderList colliders(dense)
ipH1s
const int ipH1s
Definition: iso.h:29
STATIC
#define STATIC
Definition: cddefines.h:118
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:447
hydro
t_hydro hydro
Definition: hydrogenic.cpp:5
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
PrtHydroTrace1
STATIC void PrtHydroTrace1(long ipISO, long nelem)
Definition: hydrolevel.cpp:69
ipH2p
const int ipH2p
Definition: iso.h:31
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:617
ipH2s
const int ipH2s
Definition: iso.h:30
t_iso_sp::RadRec_effec
double RadRec_effec
Definition: iso.h:548
t_trace::lgHBug
bool lgHBug
Definition: trace.h:82
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:64
LIMELM
const int LIMELM
Definition: cddefines.h:307
t_ionbal::RateIoniz
double *** RateIoniz
Definition: ionbal.h:177
t_hydro::dstfe2lya
realnum dstfe2lya
Definition: hydrogenic.h:96
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
t_dense::eden
double eden
Definition: dense.h:201
MAX2
#define MAX2(a, b)
Definition: cddefines.h:828
t_ionbal::RateIonizTot
double RateIonizTot(long nelem, long ion) const
Definition: ionbal.cpp:223
t_iso_sp::xIonSimple
double xIonSimple
Definition: iso.h:496
fprintf
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
t_secondaries::csupra
realnum ** csupra
Definition: secondaries.h:33
HydroLevel
void HydroLevel(long ipISO, long int ipZ)
secondaries
t_secondaries secondaries
Definition: secondaries.cpp:5
t_phycon::te
double te
Definition: phycon.h:21
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:348
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:668
t_iso_sp::numLevels_local
long int numLevels_local
Definition: iso.h:529
t_ionbal::CotaRate
realnum CotaRate[LIMELM]
Definition: ionbal.h:236
Generated by   doxygen 1.8.5