cloudy  trunk
 All Data Structures 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-2008 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  collider;
23 
24  DEBUG_ENTRY( "PrtHydroTrace1a()" );
25 
26  /* >>chng 05 aug 17, must use real electron density for collisional ionization of H
27  * since in Leiden v4 pdr with its artificial high temperature coll ion can be important
28  * H on H is homonuclear and scaling laws for other elements does not apply
29  * next two were multiplied by dense.EdenHCorr and changed to dense.eden for H,
30  * EdenHCorr for rest */
31  if( nelem==ipHYDROGEN )
32  {
33  /* special version for H0 onto H0 */
34  collider = dense.EdenHontoHCorr;
35  }
36  else
37  {
38  collider = dense.EdenHCorr;
39  }
40 
41  /* identify how atom is ionized for full trace */
42  if( iso.xIonSimple[ipH_LIKE][nelem] > 0. )
43  {
44  /* fraction of ionization due to photoionization */
45  phtfrc = iso.gamnc[ipH_LIKE][nelem][ipH1s]/((dense.eden*(iso.RadRec_effec[ipH_LIKE][nelem] +
46  ionbal.CotaRate[nelem]) )*
47  iso.xIonSimple[ipH_LIKE][nelem]);
48 
49  /* fraction of ionization due to collisional ionization
50  * >>chng 05 aug 17 from EdenHCorr to collider */
51  colfrc = (iso.ColIoniz[ipH_LIKE][nelem][ipH1s]*collider )/
52  ((dense.eden*(iso.RadRec_effec[ipH_LIKE][nelem] +
53  ionbal.CotaRate[0]) )*
54  iso.xIonSimple[ipH_LIKE][nelem]);
55 
56  /* fraction of ionization due to secondary ionization */
57  secfrc = secondaries.csupra[nelem][nelem]/((dense.eden*(iso.RadRec_effec[ipH_LIKE][nelem] +
58  ionbal.CotaRate[0]) )*
59  iso.xIonSimple[ipH_LIKE][nelem]);
60  }
61  else
62  {
63  phtfrc = 0.;
64  colfrc = 0.;
65  secfrc = 0.;
66  }
67 
68  fprintf( ioQQQ, " HydroLevel Z=%2ld called, simple II/I=",nelem);
70  fprintf( ioQQQ," PhotFrc:");
71  PrintE93( ioQQQ,phtfrc);
72  fprintf(ioQQQ," ColFrc:");
73  PrintE93( ioQQQ,colfrc);
74  fprintf( ioQQQ," SecFrc");
75  PrintE93( ioQQQ, secfrc);
76  fprintf( ioQQQ," Te:");
78  fprintf( ioQQQ," eden:");
80  fprintf( ioQQQ,"\n");
81  return;
82 }
83 
84 /*PrtHydroTrace1 print trace info for hydrogen-like species */
85 STATIC void PrtHydroTrace1(long nelem )
86 {
87  long int ipHi , ipLo , i;
88 
89  DEBUG_ENTRY( "PrtHydroTrace1()" );
90 
91  fprintf( ioQQQ,
92  " HydroLevel%3ld finds arrays, with optical depths defined? %li induced 2ph=%12.3e\n",
93  nelem, iteration, Transitions[ipH_LIKE][nelem][ipH2s][ipH1s].Emis->pump );
94  /* 06 aug 28, from numLevels_max to _local. */
95  for( ipHi=ipH2s; ipHi < iso.numLevels_local[ipH_LIKE][nelem]; ipHi++ )
96  {
97  fprintf( ioQQQ, "up:%2ld", ipHi );
98  fprintf( ioQQQ, "lo" );
99  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
100  {
101  fprintf( ioQQQ, "%9ld", ipLo );
102  }
103  fprintf( ioQQQ, "\n" );
104 
105  fprintf( ioQQQ, "%3ld", ipHi );
106  fprintf( ioQQQ, " A*esc" );
107  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
108  {
109  fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Aul*
110  Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Pesc ));
111  }
112  fprintf( ioQQQ, "\n" );
113 
114  fprintf( ioQQQ, "%3ld", ipHi );
115  fprintf( ioQQQ, " A*ees" );
116  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
117  {
118  fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Aul*
119  Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Pelec_esc ));
120  }
121  fprintf( ioQQQ, "\n" );
122 
123  fprintf( ioQQQ, "%3ld", ipHi );
124  fprintf( ioQQQ, " tauin" );
125  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
126  {
127  fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn ));
128  }
129  fprintf( ioQQQ, "\n" );
130 
131  fprintf( ioQQQ, "%3ld", ipHi );
132  fprintf( ioQQQ, " t tot" );
133  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
134  {
135  fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot ));
136  }
137  fprintf( ioQQQ, "\n" );
138 
139  fprintf( ioQQQ, "%3ld", ipHi );
140  fprintf( ioQQQ, " Esc " );
141  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
142  {
143  fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Pesc ));
144  }
145  fprintf( ioQQQ, "\n" );
146 
147  fprintf( ioQQQ, "%3ld", ipHi );
148  fprintf( ioQQQ, " Eesc " );
149  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
150  {
151  fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Pelec_esc ));
152  }
153  fprintf( ioQQQ, "\n" );
154 
155  fprintf( ioQQQ, "%3ld", ipHi );
156  fprintf( ioQQQ, " Dest " );
157  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
158  {
159  fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Pdest) );
160  }
161  fprintf( ioQQQ, "\n" );
162 
163  fprintf( ioQQQ, "%3ld", ipHi );
164  fprintf( ioQQQ, " A*dst" );
165  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
166  {
167  fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Aul*
168  Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Pdest ));
169  }
170  fprintf( ioQQQ, "\n" );
171 
172  fprintf( ioQQQ, "%3ld", ipHi );
173  fprintf( ioQQQ, " StrkE" );
174  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
175  {
176  fprintf( ioQQQ,PrintEfmt("%9.2e", iso.pestrk[ipH_LIKE][nelem][ipLo][ipHi] ));
177  }
178  fprintf( ioQQQ, "\n" );
179 
180  fprintf( ioQQQ, "%3ld", ipHi );
181  fprintf( ioQQQ, " B(ul)" );
182  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
183  {
184  fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->pump*
185  StatesElem[ipH_LIKE][nelem][ipLo].g/StatesElem[ipH_LIKE][nelem][ipHi].g ));
186  }
187  fprintf( ioQQQ, "\n" );
188 
189  fprintf( ioQQQ, "%3ld", ipHi );
190  fprintf( ioQQQ, " tcont" );
191  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
192  {
193  fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauCon ));
194  }
195  fprintf( ioQQQ, "\n" );
196 
197  fprintf( ioQQQ, "%3ld", ipHi );
198  fprintf( ioQQQ, " C(ul)" );
199  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
200  {
201  fprintf( ioQQQ,PrintEfmt("%9.2e",
202  Transitions[ipH_LIKE][nelem][ipHi][ipLo].Coll.ColUL*dense.eden ));
203  }
204  fprintf( ioQQQ, "\n" );
205 
206  if( ipHi == 2 )
207  {
208  fprintf( ioQQQ, " FeIIo");
209  fprintf( ioQQQ,PrintEfmt("%9.2e",
210  hydro.dstfe2lya* Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Aul ));
211  fprintf( ioQQQ, "\n");
212  }
213  }
214 
215  fprintf( ioQQQ, " " );
216  /* 06 aug 28, from numLevels_max to _local. */
217  for( i=1; i < iso.numLevels_local[ipH_LIKE][nelem]; i++ )
218  {
219  fprintf( ioQQQ, "%9ld", i );
220  }
221  fprintf( ioQQQ, "\n" );
222  return;
223 }
224 
225 /*HydroLevel calls iso_level to solve for ionization balance
226  * and level populations of model hydrogen atom */
227 void HydroLevel(long int nelem)
228 {
229  long int i;
230  double collider;
231 
232  int ipISO = ipH_LIKE;
233 
234  DEBUG_ENTRY( "HydroLevel()" );
235 
236  /* check that we were called with valid charge */
237  ASSERT( nelem >= 0);
238  ASSERT( nelem < LIMELM );
239 
240  /* >>chng 05 aug 17, must use real electron density for collisional ionization of H
241  * since in Leiden v4 pdr with its artificial high temperature coll ion can be important
242  * H on H is homonuclear and scaling laws for other elements does not apply
243  * next two were multiplied by dense.EdenHCorr and changed to dense.eden for H,
244  * EdenHCorr for rest */
245  if( ipISO == ipH_LIKE && nelem==ipHYDROGEN )
246  {
247  /* special version for H0 onto H0 */
248  collider = dense.EdenHontoHCorr;
249  }
250  else
251  {
252  collider = dense.EdenHCorr;
253  }
254 
255  /* option to print some rates */
256  if( (trace.lgTrace && trace.lgIsoTraceFull[ipISO]) && (nelem == trace.ipIsoTrace[ipISO]) )
257  PrtHydroTrace1(nelem);
258 
259  if( trace.lgHBug && trace.lgTrace )
260  PrtHydroTrace1a(nelem);
261 
262  /* this is main trace h-like printout */
263  if( (trace.lgIsoTraceFull[ipISO] && trace.lgTrace) && (nelem == trace.ipIsoTrace[ipISO]) )
264  {
265  fprintf( ioQQQ, " HLEV HGAMNC" );
266  PrintE93( ioQQQ, iso.gamnc[ipISO][nelem][ipH1s] );
267  /* 06 aug 28, from numLevels_max to _local. */
268  for( i=ipH2s; i < iso.numLevels_local[ipISO][nelem]; i++ )
269  {
270  fprintf(ioQQQ,PrintEfmt("%9.2e", iso.gamnc[ipISO][nelem][i] ));
271  }
272  fprintf( ioQQQ, "\n" );
273 
274  fprintf( ioQQQ, " HLEV TOTCAP" );
275  /* 06 aug 28, from numLevels_max to _local. */
276  for( i=1; i < iso.numLevels_local[ipISO][nelem]; i++ )
277  {
278  fprintf(ioQQQ,PrintEfmt("%9.2e", iso.RateCont2Level[ipISO][nelem][i]/dense.eden ));
279  }
280  fprintf( ioQQQ," tot");
281  fprintf( ioQQQ,PrintEfmt("%10.2e", ionbal.RateRecomTot[nelem][nelem-ipISO]/dense.eden ) );
282  fprintf( ioQQQ, "\n" );
283 
284  fprintf( ioQQQ, " HLEV IND Rc" );
285  /* 06 aug 28, from numLevels_max to _local. */
286  for( i=ipH1s; i < iso.numLevels_local[ipISO][nelem]; i++ )
287  {
288  fprintf(ioQQQ,PrintEfmt("%9.2e", iso.RecomInducRate[ipISO][nelem][i]*iso.PopLTE[ipISO][nelem][i] ));
289  }
290  fprintf( ioQQQ, "\n" );
291 
292  /* incuded recombination rate coefficients */
293  fprintf( ioQQQ, " IND Rc LTE " );
294  /* 06 aug 28, from numLevels_max to _local. */
295  for( i=ipH1s; i < iso.numLevels_local[ipISO][nelem]; i++ )
296  {
297  fprintf(ioQQQ,PrintEfmt("%9.2e",
298  iso.gamnc[ipISO][nelem][i]*iso.PopLTE[ipISO][nelem][i] ));
299  }
300  fprintf( ioQQQ, "\n" );
301 
302  /* LTE level populations */
303  fprintf( ioQQQ, " HLEV HLTE" );
304  /* 06 aug 28, from numLevels_max to _local. */
305  for( i=ipH1s; i < iso.numLevels_local[ipISO][nelem]; i++ )
306  {
307  fprintf(ioQQQ,PrintEfmt("%9.2e", iso.PopLTE[ipISO][nelem][i] ));
308  }
309  fprintf( ioQQQ, "\n" );
310 
311  /* fraction of total ionization due to collisions from given level */
312  fprintf( ioQQQ, " HLEVfr cion" );
313  /* 06 aug 28, from numLevels_max to _local. */
314  for( i=ipH1s; i < iso.numLevels_local[ipISO][nelem]; i++ )
315  {
316  fprintf(ioQQQ,PrintEfmt("%9.2e",
317  iso.ColIoniz[ipISO][nelem][i]*
318  StatesElem[ipISO][nelem][i].Pop*collider/MAX2(1e-30,iso.RateLevel2Cont[ipISO][nelem][i]) ) );
319  }
320  fprintf( ioQQQ, "\n" );
321 
322  /* fraction of total ionization due to photoionization from given level */
323  if( ionbal.RateRecomTot[nelem][nelem]> 0. )
324  {
325  fprintf( ioQQQ, " HLEVfrPhIon" );
326  /* 06 aug 28, from numLevels_max to _local. */
327  for( i=ipH1s; i < iso.numLevels_local[ipISO][nelem]; i++ )
328  {
329  fprintf(ioQQQ,PrintEfmt("%9.2e",
330  iso.gamnc[ipISO][nelem][i]*StatesElem[ipISO][ipHYDROGEN][i].Pop/MAX2(1e-30,iso.RateLevel2Cont[ipISO][nelem][i]) ) );
331  }
332  fprintf( ioQQQ, "\n" );
333  }
334 
335  fprintf( ioQQQ, " HLEV HN" );
336  /* 06 aug 28, from numLevels_max to _local. */
337  for( i=ipH1s; i < iso.numLevels_local[ipISO][nelem]; i++ )
338  {
339  fprintf(ioQQQ,PrintEfmt("%9.2e", StatesElem[ipISO][nelem][i].Pop ));
340  }
341  fprintf( ioQQQ, "\n" );
342 
343  fprintf( ioQQQ, " HLEV b(n)" );
344  /* 06 aug 28, from numLevels_max to _local. */
345  for( i=ipH1s; i < iso.numLevels_local[ipISO][nelem]; i++ )
346  {
347  fprintf(ioQQQ,PrintEfmt("%9.2e", iso.DepartCoef[ipISO][nelem][i] ));
348  }
349  fprintf( ioQQQ, "\n" );
350 
351  fprintf( ioQQQ, " HLEV X12tot");
352  fprintf(ioQQQ,PrintEfmt("%9.2e", secondaries.x12tot ));
353  fprintf( ioQQQ," Grn dest:");
354  fprintf(ioQQQ,PrintEfmt("%9.2e",
355  ionbal.RateIonizTot[nelem][nelem] ));
356  fprintf(ioQQQ, "\n");
357  }
358 
359 
360  /* >>chng 05 mar 24,
361  * renormalize the populations and emission of H atom to agree with chemistry
362  * these were not being kept parallel with chemistry, and caused large changes in O+
363  * abundance when finally done */
364  fixit(); /* this probably does not need to be called here. */
365  if( nelem == ipHYDROGEN )
366  HydroRenorm();
367 
368  if( trace.lgTrace )
369  {
370  /* iso.RecomTotal[nelem] is gross rec coef, computed here while filling in matrix
371  * elements, all physical processes included.
372  * RadRec_effec is total effective radiative only */
373  fprintf( ioQQQ, " HydroLevel%3ld return %s te=",
374  nelem,
375  iso.chTypeAtomUsed[ipISO][nelem] );
377  fprintf( ioQQQ," ion/atom=%.4e",iso.pop_ion_ov_neut[ipISO][nelem]);
378 
379  fprintf( ioQQQ," simple=%.4e",iso.xIonSimple[ipISO][nelem]);
380 
381  fprintf( ioQQQ," b1=%.2e",iso.DepartCoef[ipISO][nelem][ipH1s]);
382 
383  fprintf( ioQQQ," ion rate=%.4e",ionbal.RateIonizTot[nelem][nelem-ipISO]);
384 
385  fprintf( ioQQQ," TotRec=%.4e",ionbal.RateRecomTot[nelem][nelem-ipISO]);
386 
387  fprintf( ioQQQ," RadRec=%.4e",iso.RadRec_effec[ipISO][nelem]);
388  fprintf( ioQQQ, "\n");
389  }
390  return;
391 }

Generated for cloudy by doxygen 1.8.4