cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
hydrot2low.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 /*HydroT2Low called to do hydrogenic level populations when temp too low for matrix,
4  * forces total of level populations to add up to iso.xIonSimple[ipISO][nelem],
5  * so results of this routine will always agree with that value */
6 #include "cddefines.h"
7 #include "taulines.h"
8 #include "iso.h"
9 #include "secondaries.h"
10 #include "ionbal.h"
11 #include "dense.h"
12 #include "trace.h"
13 #include "phycon.h"
14 #include "hydrogenic.h"
15 
16 void HydroT2Low( long int ipISO, long int nelem )
17 {
18  long int iup,
19  level,
20  low;
21  double entries,
22  exits,
23  sum,
24  renorm,
25  collider;
26 
27  DEBUG_ENTRY( "HydroT2Low()" );
28 
29  /* check that we were called with valid charge */
30  ASSERT( nelem >= 0);
31  ASSERT( nelem < LIMELM );
32  /* >>chng 05 aug 17, use eden as collider for H and corrected eden for heavier
33  * nuclei - in hot PDR H is collisionally ionized, but not by other H0 since
34  * collision is homonuclear */
35  if( nelem==ipHYDROGEN )
36  {
37  /* special version for H0 onto H0 */
38  collider = dense.EdenHontoHCorr;
39  }
40  else
41  {
42  collider = dense.EdenHCorr;
43  }
44 
45  /* do radiative and downward collisional cascades */
46  /* 06 aug 28, from numLevels_max to _local. */
47  for( level=iso.numLevels_local[ipISO][nelem]-1; level >= 0; level-- )
48  {
49  /* captures into this level from the continuum,
50  * radiative rec plus three body rec, units s-1 */
51  entries = iso.RateCont2Level[ipISO][nelem][level];
52 
53  /* sum over all levels between current level and top of atom */
54  /* 06 aug 28, from numLevels_max to _local. */
55  for( iup=level + 1; iup < iso.numLevels_local[ipISO][nelem]; iup++ )
56  {
57  /* radiative decays from upper level to here, s-1 */
58  double RadDecay;
59 
60  RadDecay = MAX2( iso.SmallA,
61  Transitions[ipISO][nelem][iup][level].Emis->Aul*
62  (Transitions[ipISO][nelem][iup][level].Emis->Pesc +
63  Transitions[ipISO][nelem][iup][level].Emis->Pelec_esc +
64  Transitions[ipISO][nelem][iup][level].Emis->Pdest) );
65 
66  /* rate upper level (rel to continuum) decays to here, units s-1 */
67  entries += StatesElem[ipISO][nelem][iup].Pop*(RadDecay +
68  Transitions[ipISO][nelem][iup][level].Coll.ColUL*collider);
69  }
70 
71  /* now find rate this levels decays to all lower levels and continuum
72  * continuum ionization rate, s-1, first */
73  exits = iso.gamnc[ipISO][nelem][level];
74  for( low=ipH1s; low <= (level - 1); low++ )
75  {
76  /* radiative decays to all lower levels - s-1 */
77  double RadDecay;
78 
79  RadDecay = MAX2( iso.SmallA,
80  Transitions[ipISO][nelem][level][low].Emis->Aul*
81  (Transitions[ipISO][nelem][level][low].Emis->Pesc +
82  Transitions[ipISO][nelem][level][low].Emis->Pelec_esc +
83  Transitions[ipISO][nelem][level][low].Emis->Pdest) );
84 
85  /* add rad and collisional decays */
86  exits += RadDecay +
87  Transitions[ipISO][nelem][level][low].Coll.ColUL*collider;
88  }
89  if( exits > 1e-25 )
90  {
91  StatesElem[ipISO][nelem][level].Pop = entries/exits;
92 
93  }
94  else
95  {
96  StatesElem[ipISO][nelem][level].Pop = 0.;
97  }
98  }
99 
100  /* ===================================================================
101  *
102  * at this point all destruction processes have been established
103  *
104  * ==================================================================== */
105 
106  /* do simple sums for atom to ion, then fill in for rest */
107  StatesElem[ipISO][nelem][ipH1s].Pop =
108  ionbal.RateRecomTot[nelem][nelem-ipISO]/
109  iso.RateLevel2Cont[ipISO][nelem][ipH1s];
110 
111  if( ipISO==ipH_LIKE )
112  {
113 
114  StatesElem[ipISO][nelem][ipH2s].Pop =
115  (iso.RadRec_caseB[ipISO][nelem]*0.33*dense.eden +
116  StatesElem[ipISO][nelem][ipH1s].Pop*secondaries.Hx12[ipISO][nelem][ipH2s] +
117  StatesElem[ipISO][nelem][ipH2p].Pop*
118  Transitions[ipISO][nelem][ipH2p][ipH2s].Coll.ColUL*collider)/
119  (Transitions[ipISO][nelem][ipH2s][ipH1s].Emis->Aul +
120  Transitions[ipISO][nelem][ipH2s][ipH1s].Coll.ColUL*collider +
121  Transitions[ipISO][nelem][ipH2p][ipH2s].Coll.ColUL*collider*6. );
122 
123  StatesElem[ipISO][nelem][ipH2p].Pop =
124  (iso.RadRec_caseB[ipISO][nelem]*0.67*dense.eden +
125  StatesElem[ipISO][nelem][ipH1s].Pop*secondaries.Hx12[ipISO][nelem][ipH2p] +
126  StatesElem[ipISO][nelem][ipH2s].Pop*
127  Transitions[ipISO][nelem][ipH2p][ipH2s].Coll.ColUL*collider*6.)/
128  (Transitions[ipISO][nelem][ipH2p][ipH1s].Emis->Aul*
129  (Transitions[ipISO][nelem][ipH2p][ipH1s].Emis->Pesc +
130  Transitions[ipISO][nelem][ipH2p][ipH1s].Emis->Pelec_esc +
131  Transitions[ipISO][nelem][ipH2p][ipH1s].Emis->Pdest) +
132  Transitions[ipISO][nelem][ipH2p][ipH1s].Coll.ColUL * collider +
133  Transitions[ipISO][nelem][ipH2p][ipH2s].Coll.ColUL * collider);
134  }
135 
136 
137  /* now renormalize to simple ionization ratio calculated in hydrolevel */
138  sum = 0.;
139  /* 06 aug 28, from numLevels_max to _local. */
140  for( level=0; level < iso.numLevels_local[ipISO][nelem]; level++ )
141  {
142  sum += StatesElem[ipISO][nelem][level].Pop;
143  }
144 
145  /* rescale so that inverse of sum is iso.xIonSimple[ipISO][nelem] */
146  renorm = 1. / SDIV( iso.xIonSimple[ipISO][nelem]) / SDIV( sum);
147 
148  for( level=0; level < iso.numLevels_local[ipISO][nelem]; level++ )
149  {
150  StatesElem[ipISO][nelem][level].Pop *= renorm;
151  ASSERT( StatesElem[ipISO][nelem][level].Pop<BIGFLOAT );
152  if( iso.PopLTE[ipISO][nelem][level] > 1e-25 )
153  {
154  iso.DepartCoef[ipISO][nelem][level] =
155  StatesElem[ipISO][nelem][level].Pop/
156  (iso.PopLTE[ipISO][nelem][level]*dense.eden);
157  }
158  else
159  {
160  iso.DepartCoef[ipISO][nelem][level] = 0.;
161  }
162  }
163 
164  if( trace.lgHBug && trace.lgTrace )
165  {
166  fprintf( ioQQQ, " LOW TE,=%10.3e HN(1)=%10.3e rec=%10.3e Hgamnc(1s)=%10.3e\n",
167  phycon.te, StatesElem[ipISO][nelem][ipH1s].Pop, ionbal.RateRecomTot[nelem][nelem-ipISO],
168  iso.gamnc[ipISO][nelem][ipH1s] );
169  }
170 
171  if( trace.lgTrace )
172  {
173  fprintf( ioQQQ,
174  " HydroT2Low return, LOW TE used, HII/HI=%.3e simple=%.3e IonRate=%.3e RecCo=%.3e\n",
175  iso.pop_ion_ov_neut[ipISO][nelem],
176  iso.xIonSimple[ipISO][nelem],
177  iso.RateLevel2Cont[ipISO][nelem][ipH1s],
178  ionbal.RateRecomTot[nelem][nelem-ipISO] );
179  }
180  return;
181 }

Generated for cloudy by doxygen 1.8.1.1