cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
iso_ionize_recombine.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 /*iso_ionize_recombine find state specific creation and destruction rates for iso sequences */
4 /*ChargeTransferUpdate update rate of ct ionization and recombination for H atoms */
5 #include "cddefines.h"
6 #include "ionbal.h"
7 #include "atmdat.h"
8 #include "dense.h"
9 #include "hmi.h"
10 #include "mole.h"
11 #include "secondaries.h"
12 #include "iso.h"
13 #include "trace.h"
14 /*lint -e778 constant expression evaluates to zero - in HMRATE macro */
15 /*ChargeTransferUpdate update rate of ct ionization and recombination for H atoms */
17 {
18  long int ion;
19  long int nelem;
20  /* find total rate for charge transfer ionization of hydrogen,
21  * recombination for other element is ionization of hydrogen */
23  for( nelem=ipHELIUM; nelem<LIMELM; ++nelem)
24  {
25  /* ion on the C scale, 0 is atom, so goes up to nelem+1,
26  * for helium nelem=1, ions go from 1 to 2 */
27  for( ion=1; ion<=nelem+1; ++ion )
28  {
29  double one;
30  /* we intentionally skip CT with O+ since this is in hmole */
31  /* >>chng 05 aug 17, add logic on !ionbal.lgHO_ct_chem
32  * this is flag that is true when H O charge transfer is done in
33  * chemistry, not here. set false with set HO charge transfer ionization
34  * command, so we do it here - introduced by NA for the Leiden hacks */
35  if( (nelem == ipOXYGEN) && (ion == 1) && ionbal.lgHO_ct_chem )
36  continue;
37  /* charge transfer ionization of H, recombination for other species */
38  one = atmdat.HCharExcRecTo[nelem][ion-1]*dense.xIonDense[nelem][ion]; //units s-1
39  /* this, when multiplied by dense.xIonDense[ipHYDROGEN][0], is the charge transfer
40  * ionization RATE of neutral hydrogen, units cm-3 s-1 */
41  atmdat.HCharExcIonTotal += one;
42  }
43  }
44 
45  /* >>chng 01 may 07, add this in */
46  /* charge transfer recombination of hydrogen,
47  * which is ionization of the heavy element */
49  for( nelem=ipHELIUM; nelem<LIMELM; ++nelem)
50  {
51  /* this is ion on the abundances scale, 1 is atom, so goes up to nelem+1,
52  * for helium nelem=1, ion must go up to 3 */
53  for( ion=0; ion<=nelem; ++ion )
54  {
55  /* option to skip Oo => O+ */
56  /* >>chng 05 aug 17, add logic on !ionbal.lgHO_ct_chem
57  * this is flag that is true when H O charge transfer is done in
58  * chemistry, not here. set false with set HO charge transfer ionization
59  * command, so we do it here - introduced by NA for the Leiden hacks */
60  if( (nelem == ipOXYGEN) && (ion == 0) && ionbal.lgHO_ct_chem )
61  continue;
62  /* charge transfer ionization of H, recombination for other species */
63  /* this, when multiplied by dense.xIonDense[ipHYDROGEN][1], is the charge transfer
64  * recombination RATE of neutral hydrogen, units cm-3 s-1 */
65  atmdat.HCharExcRecTotal += atmdat.HCharExcIonOf[nelem][ion]*dense.xIonDense[nelem][ion];
66  }
67  }
68 
69  /* >>logic checked 02 may 02, think it's right */
70  /* add on charge transfer ionization of helium,
71  * recombination for other element is ionization of helium */
73  /* loop up from Li, the next heavier element */
74  for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem)
75  {
76  double hold_one = 0.;
77  /* check that array bounds not exceeded */
78  /* ion on the C scale, 0 is atom, so goes up to nelem+1,
79  * for helium nelem=1, ions go from 1 to 2 */
80  for( ion=1; ion<=nelem+1; ++ion )
81  {
82  /* charge transfer ionization of He, recombination for other species */
83  hold_one = atmdat.HeCharExcRecTo[nelem][ion-1]*dense.xIonDense[nelem][ion];
84  atmdat.HeCharExcIonTotal += hold_one;
85  }
86  }
87 
88  /* >>chng 04 jul 02,
89  * add on charge transfer reactions of He-H */
91 
92  /* charge transfer recombination of helium,
93  * which is ionization of the heavy element */
95  for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem)
96  {
97  /* this is ion on the physics scale, 1 is atom, so goes up to nelem+1,
98  * for helium nelem=1, ion must go up to 3 */
99  for( ion=0; ion<=nelem; ++ion )
100  {
101  /* charge transfer recombination of He, ionization for other species */
103  atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[nelem][ion];
104  }
105  }
106  /* >>chng 04 jul 02
107  * Add on charge transfer reactions of He+ +H0 -> He0 + H+ */
109 
110  return;
111 }
112 
113 /*iso_ionize_recombine find state specific creation and destruction rates for iso sequences */
115  /* iso sequence, hydrogen or helium for now */
116  long ipISO ,
117  /* the chemical element, 0 for hydrogen */
118  long int nelem )
119 {
120  long int level;
121  double Recom3Body;
122 
123  DEBUG_ENTRY( "iso_ionize_recombine()" );
124 
125  ASSERT( ipISO >= 0 && ipISO < NISO );
126  ASSERT( nelem >= 0 && nelem < LIMELM );
127 
128  /* get total charge transfer ionization rate if this is hydrogen itself,
129  * routine also does CT for helium but all models must have hydrogen, so
130  * he is done at this stage */
131  if( ipISO==ipH_LIKE && nelem == ipHYDROGEN )
133 
134  /* find recombination and ionization elements, will use to get simple estimate
135  * of ionization ratio below */
136  /* >>chng 06 jul 20, level should go to numLevels_local instead of numLevels_max */
137 
138  fixit(); /* must apply error to iso.gamnc */
139 
140  for( level=ipH1s; level< iso.numLevels_local[ipISO][nelem]; ++level)
141  {
142  /* all process moving level to continuum, units s-1 */
143  iso.RateLevel2Cont[ipISO][nelem][level] = iso.gamnc[ipISO][nelem][level] +
144  iso.ColIoniz[ipISO][nelem][level]* dense.EdenHCorr +
145  secondaries.csupra[nelem][nelem-ipISO]*iso.lgColl_ionize[ipISO];
146 
147  /* all processes from continuum to level n, units s-1 */
148  iso.RateCont2Level[ipISO][nelem][level] = (
149  /* radiative recombination */
150  iso.RadRecomb[ipISO][nelem][level][ipRecRad]*
151  iso.RadRecomb[ipISO][nelem][level][ipRecNetEsc] +
152 
153  /* dielectronic recombination */
154  iso.DielecRecomb[ipISO][nelem][level] +
155 
156  /* induced recombination */
157  iso.RecomInducRate[ipISO][nelem][level]*iso.PopLTE[ipISO][nelem][level] +
158 
159  /* collisional or three body recombination */
160  /* PopLTE(level,nelem) is only LTE pop when mult by n_e n_H */
161  iso.ColIoniz[ipISO][nelem][level]*dense.EdenHCorr*iso.PopLTE[ipISO][nelem][level]
162  ) * dense.eden;
163  }
164 
165  /* now charge transfer - all into/from ground, two cases, H and not H */
166  if( nelem==ipHYDROGEN )
167  {
168  /* charge transfer, hydrogen onto everything else */
169  /* charge exchange ionization, these come out of every level */
170  for( level=0; level< iso.numLevels_local[ipISO][nelem]; ++level)
171  {
172  iso.RateLevel2Cont[ipISO][nelem][level] += atmdat.HCharExcIonTotal;
173  }
174  /* charge exchange recombination */
175  iso.RateCont2Level[ipISO][nelem][0] += atmdat.HCharExcRecTotal;
176  }
177  else if( nelem==ipHELIUM && ipISO==ipHE_LIKE )
178  {
179  /* add CO here, only for he itself,
180  * destruction of CO but recombination for he */
181 
182  /* The following terms destroy He+ and thereby need to be included
183  * here. Also included are the contributions
184  * to the recombination due to hmole_step.*/
185 
186  fixit(); // should these chemistry rates be accounted for in mole.source and mole.sink?
187  double COsource = CO_source_rate("He+");
188  for( level=0; level< iso.numLevels_local[ipISO][nelem]; ++level)
189  {
190  iso.RateLevel2Cont[ipISO][nelem][level] += COsource;
191  fixit(); //need reverse rates for hmi rates just below
192  }
193  iso.RateCont2Level[ipISO][nelem][0] += CO_sink_rate("He+") +
195 
196  /* this is ioniz of He due to ct with all other gas constituents */
197  for( level=0; level< iso.numLevels_local[ipISO][nelem]; ++level)
198  {
199  iso.RateLevel2Cont[ipISO][nelem][level] += atmdat.HeCharExcIonTotal;
200  }
201  /* this is recom of He due to ct with all other gas constituents */
202  iso.RateCont2Level[ipISO][nelem][0] += atmdat.HeCharExcRecTotal;
203  }
204  else
205  {
206  for( level=0; level< iso.numLevels_local[ipISO][nelem]; ++level)
207  {
208  iso.RateLevel2Cont[ipISO][nelem][level] +=
209  atmdat.HeCharExcIonOf[nelem][nelem-ipISO]*dense.xIonDense[ipHELIUM][1] +
210  atmdat.HCharExcIonOf[nelem][nelem-ipISO]*dense.xIonDense[ipHYDROGEN][1];
211  }
212 
213  iso.RateCont2Level[ipISO][nelem][0] +=
214  atmdat.HeCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHELIUM][0] +
215  atmdat.HCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHYDROGEN][0];
216  }
217 
218 
219  /* now create sums of recombination and ionization rates for ISO species */
220  ionbal.RateRecomTot[nelem][nelem-ipISO] = 0.;
221  ionbal.RR_rate_coef_used[nelem][nelem-ipISO] = 0.;
222  Recom3Body = 0.;
223  /* >>chng 06 jul 20, level should go to numLevels_local instead of numLevels_max */
224  for( level=0; level< iso.numLevels_local[ipISO][nelem]; ++level)
225  {
226 
227  /* units of ionbal.RateRecomTot are s-1,
228  * equivalent ionization term is done after level populations are known */
229  ionbal.RateRecomTot[nelem][nelem-ipISO] += iso.RateCont2Level[ipISO][nelem][level];
230 
231  /* just the radiative recombination rate coef, cm3 s-1 */
232  ionbal.RR_rate_coef_used[nelem][nelem-ipISO] += iso.RadRecomb[ipISO][nelem][level][ipRecRad]*
233  iso.RadRecomb[ipISO][nelem][level][ipRecNetEsc];
234 
235  /* >>chng 05 jul 11, from > to >=, some very opt thick sims did block escape to zero */
236  ASSERT( ionbal.RR_rate_coef_used[nelem][nelem-ipISO]>= 0. );
237 
238  /* this is three-body recombination rate coef by itself -
239  * need factor of eden to become rate */
240  Recom3Body += iso.ColIoniz[ipISO][nelem][level]*dense.EdenHCorr*iso.PopLTE[ipISO][nelem][level];
241  }
242 
243  /* fraction of total recombs due to three body - when most are due to this
244  * small changes in temperature can produce large changes in rec coef,
245  * and so in ionization */
246  iso.RecomCollisFrac[ipISO][nelem] = Recom3Body* dense.eden / ionbal.RateRecomTot[nelem][nelem-ipISO];
247 
248  /* get simple estimate of level of ionization */
249  if( ionbal.RateRecomTot[nelem][nelem-ipISO] > 0. )
250  {
251  iso.xIonSimple[ipISO][nelem] = iso.RateLevel2Cont[ipISO][nelem][ipH1s]/ionbal.RateRecomTot[nelem][nelem-ipISO];
252  }
253  else
254  {
255  iso.xIonSimple[ipISO][nelem] = 0.;
256  }
257 
258  if( trace.lgTrace )
259  {
260  fprintf( ioQQQ, " iso_ionize_recombine iso=%2ld Z=%2ld Level2Cont[0] %10.2e RateRecomTot %10.2e xIonSimple %10.2e\n",
261  ipISO, nelem, iso.RateLevel2Cont[ipISO][nelem][0], ionbal.RateRecomTot[nelem][nelem-ipISO], iso.xIonSimple[ipISO][nelem] );
262  }
263 
264  return;
265 }
266 /*lint +e778 constant expression evaluates to zero - in HMRATE macro */

Generated for cloudy by doxygen 1.8.3.1