cloudy  trunk
 All Data Structures Namespaces 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-2022 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 "conv.h"
8 #include "atmdat.h"
9 #include "secondaries.h"
10 #include "iso.h"
11 #include "phycon.h"
12 #include "rfield.h"
13 #include "trace.h"
14 #include "freebound.h"
15 #include "dense.h"
16 /*lint -e778 constant expression evaluates to zero - in HMRATE macro */
17 /*iso_charge_transfer_update update rate of ct ionization and recombination for H atoms */
18 void iso_charge_transfer_update(long nelem1)
19 {
20  DEBUG_ENTRY( "iso_charge_transfer_update()" );
21 
22  if( nelem1 >= t_atmdat::NCX )
23  return;
24 
25  atmdat.CharExcIonTotal[nelem1] = 0.;
26  atmdat.CharExcRecTotal[nelem1] = 0.;
27 
28  // find total charge transfer rate coefficients
29  // charge transfer reactions of lighter species
30  for ( long nelem=ipHYDROGEN; nelem < nelem1; ++nelem)
31  {
32  atmdat.CharExcIonTotal[nelem1] += atmdat.CharExcIonOf[nelem][nelem1][0]*dense.xIonDense[nelem][1];
33  long ipISO=nelem;
34  atmdat.CharExcRecTotal[nelem1] += atmdat.CharExcRecTo[nelem][nelem1][0]*iso_sp[ipISO][nelem].st[0].Pop();
35  }
36 
37  // and of heavier species
38  for( long nelem=nelem1+1; nelem<LIMELM; ++nelem)
39  {
40  for( long ion=0; ion<=nelem; ++ion )
41  {
42  // charge transfer ionization of nelem1, recombination for other species, cm-3 s-1
43  atmdat.CharExcIonTotal[nelem1] += atmdat.CharExcRecTo[nelem1][nelem][ion]*dense.xIonDense[nelem][ion+1];
44  // charge transfer recombination of nelem1, ionization for other species, cm-3 s-1
45  atmdat.CharExcRecTotal[nelem1] += atmdat.CharExcIonOf[nelem1][nelem][ion]*dense.xIonDense[nelem][ion];
46  }
47  }
48 
49  return;
50 }
51 
52 /*iso_ionize_recombine find state specific creation and destruction rates for iso sequences */
54  /* iso sequence, hydrogen or helium for now */
55  long ipISO ,
56  /* the chemical element, 0 for hydrogen */
57  long int nelem )
58 {
59  long int level;
60  double Recom3Body;
61 
62  DEBUG_ENTRY( "iso_ionize_recombine()" );
63 
64  ASSERT( ipISO >= 0 && ipISO < NISO );
65  ASSERT( nelem >= 0 && nelem < LIMELM );
66 
67  /* find recombination and ionization elements, will use to get simple estimate
68  * of ionization ratio below */
69  /* >>chng 06 jul 20, level should go to numLevels_local instead of numLevels_max */
70 
71  fixit("must apply error to iso.gamnc");
72 
73  for( level=ipH1s; level< iso_sp[ipISO][nelem].numLevels_local; ++level)
74  {
75  // This term is intended to capture LTE in DR-dominated plasmas.
76  // It is scaled by Occ num / Occ num at STE to ensure correct limits w.r.t. ambient radiation.
77  double DR_reverse = 0.;
78  if( ipISO > ipH_LIKE && iso_sp[ipISO][nelem].fb[level].PopLTE > SMALLFLOAT )
79  {
80  long indexIP = iso_sp[ipISO][nelem].fb[level].ipIsoLevNIonCon-1;
81  double OccNum = rfield.OccNumbIncidCont[indexIP] + rfield.OccNumbContEmitOut[indexIP];
82  double OccNumBB = 1./(exp( iso_sp[ipISO][nelem].fb[level].xIsoLevNIonRyd / phycon.te_ryd ) - 1. );
83  double RelOccNum = MIN2( 1., OccNum / OccNumBB );
84  //fprintf( ioQQQ, "DEBUGGG DR_reverse %li\t%li\t%e\t%e\t%e\t%e\t%e\t%e\n", level, indexIP, RelOccNum, OccNum, OccNumBB,
85  // rfield.anu(indexIP), iso_sp[ipISO][nelem].fb[level].xIsoLevNIonRyd, rfield.anu(indexIP+1) );
86  DR_reverse = iso_sp[ipISO][nelem].fb[level].DielecRecomb * RelOccNum /
87  iso_sp[ipISO][nelem].fb[level].PopLTE;
88  }
89 
90  /* all process moving level to continuum, units s-1 */
91  iso_sp[ipISO][nelem].fb[level].RateLevel2Cont = iso_sp[ipISO][nelem].fb[level].gamnc +
92  iso_sp[ipISO][nelem].fb[level].ColIoniz* dense.EdenHCorr +
93  DR_reverse +
94  secondaries.csupra[nelem][nelem-ipISO]*iso_ctrl.lgColl_ionize[ipISO];
95 
96  /* all processes from continuum to level n, units s-1 */
97  iso_sp[ipISO][nelem].fb[level].RateCont2Level = (
98  /* radiative recombination */
99  iso_sp[ipISO][nelem].fb[level].RadRecomb[ipRecRad]*
100  iso_sp[ipISO][nelem].fb[level].RadRecomb[ipRecNetEsc] +
101 
102  /* dielectronic recombination */
103  iso_sp[ipISO][nelem].fb[level].DielecRecomb +
104 
105  /* induced recombination */
106  iso_sp[ipISO][nelem].fb[level].RecomInducRate*iso_sp[ipISO][nelem].fb[level].PopLTE +
107 
108  /* collisional or three body recombination */
109  /* PopLTE(level,nelem) is only LTE pop when mult by n_e n_H */
110  iso_sp[ipISO][nelem].fb[level].ColIoniz*dense.EdenHCorr*iso_sp[ipISO][nelem].fb[level].PopLTE
111  ) * dense.eden;
112 
113  ASSERT( iso_sp[ipISO][nelem].fb[level].RadRecomb[ipRecRad] >= 0. );
114  ASSERT( iso_sp[ipISO][nelem].fb[level].RadRecomb[ipRecNetEsc] >= 0. );
115  ASSERT( iso_sp[ipISO][nelem].fb[level].DielecRecomb >= 0. );
116  ASSERT( iso_sp[ipISO][nelem].fb[level].RecomInducRate >= 0. );
117  ASSERT( iso_sp[ipISO][nelem].fb[level].PopLTE >= 0. );
118  ASSERT( iso_sp[ipISO][nelem].fb[level].ColIoniz >= 0. );
119  ASSERT( iso_sp[ipISO][nelem].fb[level].RateCont2Level >= 0. );
120 
121  if( iso_ctrl.lgRandErrGen[ipISO] )
122  {
123  iso_sp[ipISO][nelem].fb[level].RateCont2Level *=
124  iso_sp[ipISO][nelem].ex[ iso_sp[ipISO][nelem].numLevels_max ][level].ErrorFactor[IPRAD];
125  }
126  }
127 
128  /* now create sums of recombination and ionization rates for ISO species */
129  ionbal.RateRecomIso[nelem][ipISO] = 0.;
130  ionbal.RR_rate_coef_used[nelem][nelem-ipISO] = 0.;
131  Recom3Body = 0.;
132  /* >>chng 06 jul 20, level should go to numLevels_local instead of numLevels_max */
133  for( level=0; level< iso_sp[ipISO][nelem].numLevels_local; ++level)
134  {
135 
136  /* units of ionbal.RateRecomTot are s-1,
137  * equivalent ionization term is done after level populations are known */
138  ionbal.RateRecomIso[nelem][ipISO] += iso_sp[ipISO][nelem].fb[level].RateCont2Level;
139 
140  /* just the radiative recombination rate coef, cm3 s-1 */
141  ionbal.RR_rate_coef_used[nelem][nelem-ipISO] += iso_sp[ipISO][nelem].fb[level].RadRecomb[ipRecRad]*
142  iso_sp[ipISO][nelem].fb[level].RadRecomb[ipRecNetEsc];
143 
144  /* >>chng 05 jul 11, from > to >=, some very opt thick sims did block escape to zero */
145  ASSERT( ionbal.RR_rate_coef_used[nelem][nelem-ipISO]>= 0. );
146 
147  /* this is three-body recombination rate coef by itself -
148  * need factor of eden to become rate */
149  Recom3Body += iso_sp[ipISO][nelem].fb[level].ColIoniz*dense.EdenHCorr*iso_sp[ipISO][nelem].fb[level].PopLTE;
150  }
151 
152  /* fraction of total recombinations due to three body - when most are due to this
153  * small changes in temperature can produce large changes in recombination coefficient,
154  * and so in ionization */
155  iso_sp[ipISO][nelem].RecomCollisFrac = Recom3Body* dense.eden / ionbal.RateRecomIso[nelem][ipISO];
156 
157  /* very first pass through here rate RateIoniz not yet evaluated */
158  if( conv.nTotalIoniz==0 )
159  ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] = iso_sp[ipISO][nelem].fb[0].RateLevel2Cont;
160 
161  /* get simple estimate of atom to ion ratio */
162  if( ionbal.RateRecomIso[nelem][ipISO] > 0. )
163  {
164  iso_sp[ipISO][nelem].xIonSimple = ionbal.RateIonizTot(nelem,nelem-ipISO)/ionbal.RateRecomIso[nelem][ipISO];
165  }
166  else
167  {
168  iso_sp[ipISO][nelem].xIonSimple = 0.;
169  }
170 
171  if( trace.lgTrace && (trace.lgHBug||trace.lgHeBug) )
172  {
173  fprintf( ioQQQ, " iso_ionize_recombine iso=%2ld Z=%2ld Level2Cont[0] %10.2e RateRecomTot %10.2e xIonSimple %10.2e\n",
174  ipISO, nelem, iso_sp[ipISO][nelem].fb[0].RateLevel2Cont, ionbal.RateRecomIso[nelem][ipISO], iso_sp[ipISO][nelem].xIonSimple );
175  }
176 
177  return;
178 }
179 /*lint +e778 constant expression evaluates to zero - in HMRATE macro */
t_atmdat atmdat
Definition: atmdat.cpp:6
qList st
Definition: iso.h:482
double EdenHCorr
Definition: dense.h:227
bool lgHeBug
Definition: trace.h:79
bool lgColl_ionize[NISO]
Definition: iso.h:365
const realnum SMALLFLOAT
Definition: cpu.h:246
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
const int NISO
Definition: cddefines.h:311
t_conv conv
Definition: conv.cpp:5
realnum * OccNumbContEmitOut
Definition: rfield.h:62
t_phycon phycon
Definition: phycon.cpp:6
const int ipRecNetEsc
Definition: cddefines.h:331
FILE * ioQQQ
Definition: cddefines.cpp:7
bool lgRandErrGen[NISO]
Definition: iso.h:430
vector< freeBound > fb
Definition: iso.h:481
#define MIN2(a, b)
Definition: cddefines.h:803
#define IPRAD
Definition: iso.h:88
t_dense dense
Definition: global.cpp:15
double ** RateRecomIso
Definition: ionbal.h:197
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
double CharExcRecTotal[NCX]
Definition: atmdat.h:310
double CharExcIonTotal[NCX]
Definition: atmdat.h:310
t_trace trace
Definition: trace.cpp:5
t_ionbal ionbal
Definition: ionbal.cpp:8
const int ipH1s
Definition: iso.h:29
bool lgTrace
Definition: trace.h:12
void iso_charge_transfer_update(long nelem)
t_rfield rfield
Definition: rfield.cpp:9
const int ipRecRad
Definition: cddefines.h:333
realnum * OccNumbIncidCont
Definition: rfield.h:117
double RecomCollisFrac
Definition: iso.h:551
long int nTotalIoniz
Definition: conv.h:159
multi_arr< extra_tr, 2 > ex
Definition: iso.h:478
#define ASSERT(exp)
Definition: cddefines.h:613
bool lgHBug
Definition: trace.h:82
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:308
double *** RateIoniz
Definition: ionbal.h:180
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
double te_ryd
Definition: phycon.h:27
double eden
Definition: dense.h:201
double RateIonizTot(long nelem, long ion) const
Definition: ionbal.cpp:224
double xIonSimple
Definition: iso.h:496
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
realnum ** csupra
Definition: secondaries.h:33
long int numLevels_max
Definition: iso.h:524
#define fixit(a)
Definition: cddefines.h:417
t_secondaries secondaries
Definition: secondaries.cpp:5
double CharExcIonOf[NCX][LIMELM][LIMELM+1]
Definition: atmdat.h:300
const int ipHYDROGEN
Definition: cddefines.h:349
double CharExcRecTo[NCX][LIMELM][LIMELM+1]
Definition: atmdat.h:300
void iso_ionize_recombine(long ipISO, long nelem)
long int numLevels_local
Definition: iso.h:529
double ** RR_rate_coef_used
Definition: ionbal.h:210