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