cloudy trunk
Loading...
Searching...
No Matches
atmdat_char_tran.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/*ChargTranEval fill in the CharExcIonOf[ipHYDROGEN] and Rec arrays with Kingdon's fitted CT with H */
4/*ChargTranSumHeat sum net heating due to charge transfer, called in HeatSum */
5/*HCTIon H charge transfer ionization*/
6/*HCTRecom H charge transfer recombination */
7/*ChargTranPun, save charge transfer coefficient */
8/*MakeHCTData holds data for charge transfer fits */
9#include "cddefines.h"
10#include "phycon.h"
11#include "physconst.h"
12#include "abund.h"
13#include "dense.h"
14#include "iso.h"
15#include "thermal.h"
16#include "mole.h"
17#include "elementnames.h"
18#include "heavy.h"
19#include "trace.h"
20#include "conv.h"
21#include "atmdat.h"
22#include "taulines.h"
23
24/*HCTion H charge transfer ionization, H+ + A => H + A+ */
25STATIC double HCTIon(long int ion,
26 long int nelem);
27
28/*HCTRecom H charge transfer recombination, H + A+ => H+ + A */
29STATIC double HCTRecom(long int ion,
30 long int nelem);
31
32/* the translated block data */
33STATIC void MakeHCTData(void);
34
35/* the structures for storing the charge transfer data, these are filled in
36 * at the end of this file, in what used to be a block data */
37static double CTIonData[LIMELM][4][8];
38static double CTRecombData[LIMELM][4][7];
39/* this will be flag for whether or not charge transfer data
40 * have been initialized */
41static bool lgCTDataDefined = false;
42
43/*ChargTranEval update charge trans rate coefficients if temperature has changed */
44void ChargTranEval( void )
45{
46 long int ion,
47 nelem;
48 double a, b, c, a_op, b_op, c_op, d_op, e_op, f_op, a_o,
49 b_o, c_o, d_o, e_o, f_o, g_o;
50
51 static double TeUsed = -1.;
52
53 DEBUG_ENTRY( "ChargTranEval()" );
54
55 /* first is to force reevaluation on very first call */
56 if( !conv.nTotalIoniz || !fp_equal(phycon.te,TeUsed) )
57 {
58 /* refresh hydrogen charge transfer arrays */
59 /* >>chng 01 apr 25, lower limit had been 2, lithium, changed to 1, helium */
60 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
61 {
62 for( ion=0; ion <= nelem; ion++ )
63 {
64 /* >>chng 01 apr 28, add factors to turn off ct,
65 * had previously been handled downstream */
66
67 /* CharExcIonOf[ipHYDROGEN][nelem][ion]*hii is ion => ion+1 for nelem */
68 /* charge transfer ionization O^0 + H+ -> O+ + H0
69 * is CharExcIonOf[ipHYDROGEN][ipOXYGEN][0]*dense.xIonDense[ipHYDROGEN][1]
70 * charge transfer recombination of atomic hydrogen is
71 * CharExcIonOf[ipHYDROGEN][ipOXYGEN][0]*dense.xIonDense[ipOXYGEN][0] */
72 atmdat.CharExcIonOf[ipHYDROGEN][nelem][ion] = HCTIon(ion,nelem);
73
74 /* CharExcRecTo[ipHYDROGEN][nelem][ion]*hi is ion+1 => ion of nelem */
75 /* charge transfer recombination O+ + H0 -> O^0 + H+ is
76 * CharExcRecTo[ipHYDROGEN][ipOXYGEN][0]*dense.xIonDense[ipHYDROGEN][0]
77 * charge transfer ionization of atomic hydrogen is
78 * CharExcRecTo[ipHYDROGEN][ipOXYGEN][0]*dense.xIonDense[ipOXYGEN][1] */
79 atmdat.CharExcRecTo[ipHYDROGEN][nelem][ion] = HCTRecom(ion,nelem);
80 }
81
82 /* zero out helium charge transfer arrays */
83 for( ion=0; ion < LIMELM; ion++ )
84 {
85 atmdat.CharExcIonOf[ipHELIUM][nelem][ion] = 0.;
86 atmdat.CharExcRecTo[ipHELIUM][nelem][ion] = 0.;
87 }
88 }
89
90 /* The above included only the radiative charge transfer from
91 * Stancil et al 1998. must explicitly add on the ct fitted by Kingdon & Ferland,
92 * The process H0 + He+ -> He0 + H+ */
93 atmdat.CharExcRecTo[ipHYDROGEN][ipHELIUM][0] += 7.47e-15*pow(phycon.te/1e4,2.06)*
94 (1.+9.93*sexp(3.89e-4*phycon.te) );
95
96 /* >>chng 04 jun 21 -- NPA. Put in the charge transfer rate for:
97 He+ + H => He + H+ as defined in the UMIST database. This is only
98 used if the "set UMIST rates" command is used */
99 if(mole_global.lgLeidenHack)
100 {
101 atmdat.CharExcRecTo[ipHYDROGEN][ipHELIUM][0] = 4.85e-15*pow(phycon.te/300, 0.18);
102 }
103 /* >>chng 04 jun 21 -- NPA. update the charge transfer rates between hydrogen
104 and oxygen to:
105 >>refer O CT Stancil, et al. 1999, A&AS, 140, 225-234
106 Instead of using the UMIST rate, the program TableCurve was used
107 to generate a fit to the data listed in Tables 2, 3, and 4 of the
108 aforementioned reference. The following fitting equations agree
109 very well with the published data. */
110
111 /* At or below 10K, just use the value the formula's below give
112 at 10K.*/
113 /* do both O+ -> O and O -> O+ for low T limit */
114 if(phycon.te <= 10. )
115 {
116 atmdat.CharExcRecTo[ipHYDROGEN][ipOXYGEN][0] = 3.744e-10;
117 atmdat.CharExcIonOf[ipHYDROGEN][ipOXYGEN][0] = 4.749e-20;
118 }
119
120 /* this does O+ -> O for all higher temperatures */
121 if( phycon.te > 10.)
122 {
123 a_op = 2.3344302e-10;
124 b_op = 2.3651505e-10;
125 c_op = -1.3146803e-10;
126 d_op = 2.9979994e-11;
127 e_op = -2.8577012e-12;
128 f_op = 1.1963502e-13;
129
130 double lnte = phycon.alnte;
131
132 /* equation rank 53 of TableCurve */
133 atmdat.CharExcRecTo[ipHYDROGEN][ipOXYGEN][0] =
134 ((((f_op*lnte + e_op)*lnte + d_op)*lnte + c_op)*lnte + b_op)*lnte + a_op;
135 }
136
137 /* now do O -> O+
138 * The next two equations were chosen to match up at 200K, so that there
139 *are no discontinuities */
140 if((phycon.te > 10.) && (phycon.te <= 190.))
141 {
142 a = -21.134531;
143 b = -242.06831;
144 c = 84.761441;
145
146 /* equation rank 2 of TableCurve */
147 atmdat.CharExcIonOf[ipHYDROGEN][ipOXYGEN][0] = exp((a + b/SDIV(phycon.te) + c/SDIV(phycon.tesqrd)));
148 }
149
150 else if((phycon.te > 190.) && (phycon.te <= 200.))
151 {
152
153 /* We are using two fitting formula's for this rate, in order to assure no
154 sudden "jumps" in the rate, the rate between 190-200K is made to
155 increase linearly. The formula below gets the same answer as the equation
156 above at 190K, and gets the same answer as the the formula below this one
157 at 200K*/
158 atmdat.CharExcIonOf[ipHYDROGEN][ipOXYGEN][0] = 2.18733e-12*(phycon.te-190) + 1.85823e-10;
159 }
160
161 else if(phycon.te > 200.)
162 {
163
164 a_o = -7.6767404e-14;
165 b_o = -3.7282001e-13;
166 c_o = -1.488594e-12;
167 d_o = -3.6606214e-12;
168 e_o = 2.0699463e-12;
169 f_o = -2.6139493e-13;
170 g_o = 1.1580844e-14;
171
172 double lnte = phycon.alnte;
173
174 /* equation rank 120 of TableCurve */
175 atmdat.CharExcIonOf[ipHYDROGEN][ipOXYGEN][0] =
176 (((((g_o*lnte + f_o)*lnte + e_o)*lnte + d_o)*lnte + c_o)*lnte + b_o)*lnte + a_o;
177 }
178
179 /* use UMIST rates for charge transfer if UMIST command is used - disagree
180 * by 20% at 5000K and diverge at high temperature */
181 if(mole_global.lgLeidenHack)
182 {
183 atmdat.CharExcIonOf[ipHYDROGEN][ipOXYGEN][0] = HMRATE(7.31e-10,0.23,225.9);
184 atmdat.CharExcRecTo[ipHYDROGEN][ipOXYGEN][0] = HMRATE(5.66e-10,0.36,-8.60);
185 }
186
187 /* >>chng 01 may 07, following had all been +=, ok if above was zero.
188 * changed to = and added HCTOn */
189 /* >>chng 01 jan 30, add following block of CT reactions */
190 /* ionization, added as per Phillip Stancil OK, number comes from
191 * >>refer Fe CT Tielens, A. G. G. M., & Hollenbach, D. 1985a, ApJ, 294, 722-746 */
192 /* >>refer Fe CT Prasad, S. S., & Huntress, W. T. 1980, ApJS, 43, 1-35 */
193 /* the actual rate comes from the following paper: */
194 /* >>refer Fe CT Pequignot, D., & Aldrovandi, S. M. V. 1986, A&A, 161, 169-176 */
195 /* Fe0 + H+ => Fe+ + H0 */
196 /*>>chng 05 sep 15, GS, old rate had problem in predicting observed Fe I column density along HD185418.
197 *>> refer Private communication with Stancil, data taken from ORNL web site,
198 * "There is a well known problem with the Fe charge transfer rate coefficients: i.e., there are no accurate calculations nor or there
199 any experiments. For Fe + H+ -> Fe+ + H, I estimated for Gary a few years ago the value of 5.4e-9. So mid way between the two values
200 you are using. I have some notes on it in my office, but not with me. See: http://cfadc.phy.ornl.gov/astro/ps/data/home.html
201 value changed from 3e-9 to 5.4e-9 */
202
203 atmdat.CharExcIonOf[ipHYDROGEN][ipIRON][0] = 5.4e-9;
204 /*>>chng 06 sep 20 - following sets removes Fe ionization ct to be similar to Mg */
205 /*atmdat.CharExcIonOf[ipHYDROGEN][ipIRON][0] = 0.;broken();rm this line */
206
207 /* all remaining entries are from Pequignot & Aldrovandi*/
208 /* >>refer Al CT Pequignot, D., & Aldrovandi, S. M. V. 1986, A&A, 161, 169-176 */
209 /* Al0 + H+ => Al+ + H0 */
210 atmdat.CharExcIonOf[ipHYDROGEN][ipALUMINIUM][0] = 3e-9;
211
212 /* >>refer P CT Pequignot, D., & Aldrovandi, S. M. V. 1986, A&A, 161, 169-176 */
213 /* P0 + H+ => P+ + H0 */
214 atmdat.CharExcIonOf[ipHYDROGEN][ipPHOSPHORUS][0] = 1e-9;
215
216 /* >>refer Cl CT Pequignot, D., & Aldrovandi, S. M. V. 1986, A&A, 161, 169-176 */
217 /* Cl0 + H+ => Cl+ + H0 */
218 atmdat.CharExcIonOf[ipHYDROGEN][ipCHLORINE][0] = 1e-9;
219
220 /* >>refer Ti CT Pequignot, D., & Aldrovandi, S. M. V. 1986, A&A, 161, 169-176 */
221 /* Ti0 + H+ => Cl+ + H0 */
222 atmdat.CharExcIonOf[ipHYDROGEN][ipTITANIUM][0] = 3e-9;
223
224 /* >>refer Mn CT Pequignot, D., & Aldrovandi, S. M. V. 1986, A&A, 161, 169-176 */
225 /* Mn0 + H+ => Mn+ + H0 */
226 atmdat.CharExcIonOf[ipHYDROGEN][ipMANGANESE][0] = 3e-9;
227
228 /* >>refer Ni CT Pequignot, D., & Aldrovandi, S. M. V. 1986, A&A, 161, 169-176 */
229 /* Ni0 + H+ => Ni+ + H0 */
230 atmdat.CharExcIonOf[ipHYDROGEN][ipNICKEL][0] = 3e-9;
231
232 /* >>chng 01 feb 02, add following CT reaction from */
233 /* >>refer Na0 CT Dutta, C. M., Nordlander, P., Kimura, M., & Dalgarno, A. 2001, Phys. Rev. A, 63, 022709 */
234 /* this is roughly their number around 500K - they do not give explicit values, rather
235 * a small figure. Previous calculations were 5 orders of mag smaller at this temp.
236 * ND this deposits electron into n=2 */
237 /* Na0 + H+ => Na+ + H0(n=2) */
238 atmdat.CharExcIonOf[ipHYDROGEN][ipSODIUM][0] = 7e-12;
239
240 /* >>chng 05 sep 15,GS, add following CT reaction from */
241 /* >>refer Na0 CT Watanabe, A., Dutta, C. M., Nordlander, P., et al. 2002, Phys. Rev. A, 66, 044701 */
242 /* this is roughly their number around 50K - they do not give explicit values, rather
243 * a small figure. this deposits electron into n=1
244 * Na0 + H+ => Na+ + H0(n=1)
245 * add to previous rate which was for population of n=2 */
246 atmdat.CharExcIonOf[ipHYDROGEN][ipSODIUM][0] += 0.7e-12;
247
248 /* >>chng 05 sep 15, GS, add following CT reaction from
249 * >>refer K0 CT Watanabe, A., Dutta, C. M., Nordlander, P., et al. 2002, Phys. Rev. A, 66, 044701
250 * this is roughly their number around 50K - they do not give explicit values, rather
251 * a small figure.
252 * K0 + H+ => K+ + H0(n=1) */
253 atmdat.CharExcIonOf[ipHYDROGEN][ipPOTASSIUM][0] = 1.25e-12;
254
255 /* >>chng 05 sep 15, GS, add following CT reaction from
256 * >>refer S0 CT ORNL data base for charge transfer
257 * This rate is valid for 1e3 to 1e4. Due to the small value, I did not put any limit on temp.
258 * Earlier, other reactions also assume constant value
259 * S0 + H+ => H + S+ */
260 atmdat.CharExcIonOf[ipHYDROGEN][ipSULPHUR][0] = 1.e-14;
261
262 if( phycon.te < 1e5 )
263 {
264
265 /* >>chng 05 sep 15, GS, add following CT reaction from
266 * >>refer Mg0 CT ORNL data base for charge transfer,
267 * this rate is valid for temp 5e3 to 3e4, The rate goes down very fast in low temp. So I did not put a lower cut of for temp
268 * Mg0 + H+ => H + Mg+ */
269 atmdat.CharExcIonOf[ipHYDROGEN][ipMAGNESIUM][0] = 9.76e-12*pow((phycon.te/1e4),3.14)*(1. + 55.54*sexp(1.12*phycon.te/1e4));
270 /*>>chng 06 jul 20, do not allow this to fall below UMIST rate - above fit not intended for
271 * very low temperatures */
272 /*>>chng 06 aug 01, UMIST is bogus - email exchange with Phillip Stancil, late July 2006 */
273 /*atmdat.CharExcIonOf[ipHYDROGEN][ipMAGNESIUM][0] = MAX2( 1.1e-9 , atmdat.CharExcIonOf[ipHYDROGEN][ipMAGNESIUM][0]);*/
274 /*>>chng 06 sep 20 - following sets Mg ionization ct to Fe */
275 /*atmdat.CharExcIonOf[ipHYDROGEN][ipMAGNESIUM][0] = 5.4e-9;broken(); rm this line */
276
277 /* >>chng 05 sep 15, GS, add following CT reaction from
278 * >>refer Si0 CT ORNL data base for charge transfer
279 * this rate is valid for temp 1e3 to 2e5, The rate goes down very fast in low temp. So I did not put a lower cut of for temp
280 * Si0 + H+ => H + Si+ */
281 atmdat.CharExcIonOf[ipHYDROGEN][ipSILICON][0] = 0.92e-12*pow((phycon.te/1e4),1.15)*(1. + 0.80*sexp(0.24*phycon.te/1e4));
282 /*>>chng 06 jul 20, do not allow this to fall below UMIST rate - above fit not intended for
283 * very low temperatures */
285 /*>>chng 06 aug 01, UMIST rate is bogus as per Phillip Stancil emails of late July 2006 */
286 /*atmdat.CharExcIonOf[ipHYDROGEN][ipSILICON][0] = MAX2( 9.9e-10 , atmdat.CharExcIonOf[ipHYDROGEN][ipSILICON][0]);*/
287
288 /* >>chng 05 sep 15, GS, add following CT reaction from
289 * >>refer Li0 CT ORNL data base for charge transfer
290 * this rate is valid for temp 1e2 to 1e4, The rate goes down very fast in low temp. So I did not put a lower cut of for temp
291 * Li0 + H+ => H + Li+ */
292 atmdat.CharExcIonOf[ipHYDROGEN][ipLITHIUM][0] = 2.84e-12*pow((phycon.te/1e4),1.99)*(1. + 375.54*sexp(54.07*phycon.te/1e4));
295 }
296 else
297 {
299 atmdat.CharExcIonOf[ipHYDROGEN][ipMAGNESIUM][0] = 0.;
300 atmdat.CharExcIonOf[ipHYDROGEN][ipSILICON][0] = 0.;
301 atmdat.CharExcIonOf[ipHYDROGEN][ipLITHIUM][0] = 0.;
302 }
303
304 {
305 /*>>chng 06 jul 07, Terry Yun add these charge transfer reactions */
306 /*>>refer N0 CT Lin, C. Y., Stancil, P. C., Gu, J. P., et al. 2005, Phys. Rev. A, 71, 062708
307 * and combined with data from
308 *>>refer N0 CT Butler, S. E., & Dalgarno, A. 1979, ApJ, 234, 765 */
309
310 /* natural log of te */
311 double tefac = phycon.te * phycon.alnte;
312
313 /* N(4S) + H+ -> N+(3P) + H */
314 /* >>chng 06 jul 10, add exp for endoergic reaction */
315 double ct_from_n0grhp_to_npgrh0 = (1.64e-16*phycon.te - 8.76e-17*tefac + 2.41e-20*phycon.tesqrd + 9.83e-13*phycon.alnte )*
316 sexp( 10985./phycon.te ) * atmdat.lgCTOn;
317
318 /* N(2D) + H+ -> N+(3P) + H */
320 /*double ct_from_n0exhp_to_npgrh0 = 1.51e-15*phycon.te -1.61e-16*tefac + 7.74e-21*phycon.tesqrd + 1.34e-16*phycon.alnte;*/
321
322 /* N+(3P) + H -> N(4S) + H+ endoergic */
323 double ct_from_npgrh0_to_n0grhp = (1.56e-15*phycon.te - 1.79e-16*tefac + 1.15e-20*phycon.tesqrd + 1.08e-13*phycon.alnte) * atmdat.lgCTOn;
324
325 /* N+(3P) + H0 -> N(2D) + H+ */
326 /* >>chng 06 jul 10, add exp for endoergic reaction */
327 atmdat.HCharExcRecTo_N0_2D = (6.83e-16*phycon.te - 7.40e-17*tefac + 3.73e-21*phycon.tesqrd + 1.75e-15*phycon.alnte)*
328 sexp( 16680./phycon.te ) * atmdat.lgCTOn;
329
330 /* these rates are from the ground state into all possible states of the
331 * species that is produced */
332 atmdat.CharExcIonOf[ipHYDROGEN][ipNITROGEN][0] = ct_from_n0grhp_to_npgrh0;
333 atmdat.CharExcRecTo[ipHYDROGEN][ipNITROGEN][0] = ct_from_npgrh0_to_n0grhp + atmdat.HCharExcRecTo_N0_2D;
334 }
335
336 /*>>chng 06 aug 01, update O++ and N++ -- H0 CT recombination
337 *>>refer O3 CT Barragan, P., Errea, L. F., Mendez, L., et al. 2006, ApJ, 636, 544 */
338 /* O+2 + H -> O+ + H+ */
339 if( phycon.te <= 1500. )
340 {
341 atmdat.CharExcRecTo[ipHYDROGEN][ipOXYGEN][1] = 0.5337e-9*pow( (phycon.te/100.) ,-0.076);
342 }
343 else
344 {
345 atmdat.CharExcRecTo[ipHYDROGEN][ipOXYGEN][1] = 0.4344e-9 +
346 0.6340e-9*pow( log10(phycon.te/1500.) ,2.06 );
347 }
348
349 /* N+2 + H -> N+ + H+ */
350 if( phycon.te <= 1500. )
351 {
352 atmdat.CharExcRecTo[ipHYDROGEN][ipNITROGEN][1] = 0.8692e-9*pow( (phycon.te/1500.) ,0.17);
353 }
354 else if( phycon.te <= 20000. )
355 {
356 atmdat.CharExcRecTo[ipHYDROGEN][ipNITROGEN][1] = 0.9703e-9*pow( (phycon.te/10000.) ,0.058);
357 }
358 else
359 {
360 atmdat.CharExcRecTo[ipHYDROGEN][ipNITROGEN][1] = 1.0101e-9 +
361 1.4589e-9*pow( log10(phycon.te/20000.) ,2.06 );
362 }
363
364 /* ===================== helium charge transfer ====================*/
365
366 /* atmdat.CharExcIonOf[ipHELIUM] is ionization, */
367 /* [0] is Atom^0 + He+ => Atom+1 + He0
368 * [n] is Atom^+n + He+ => Atom^+n-1 + He0 */
369
370 /* atmdat.CharExcRecTo[ipHELIUM] is recombination */
371 /* [0] is Atom^+1 + He0 => Atom^0 + He^+
372 * [n] is Atom^+n+1 + He0 => Atom^+n + He^+ */
373
374 /* Carbon */
375 /* recombination */
376 /* C+3 + He => C+2 + He+ */
377 /* >>refer C3 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
378 atmdat.CharExcRecTo[ipHELIUM][ipCARBON][2] = 4.6e-19*phycon.tesqrd;
379
380 /* C+4 + He => C+3 + He+ */
381 /* >>refer C4 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
382 atmdat.CharExcRecTo[ipHELIUM][ipCARBON][3] = 1e-14;
383
384 /* ionization */
385 /* C0 + He+ => C+ + He0 */
386 /* unknown reference, from older version of the code */
387 /*atmdat.CharExcIonOf[ipHELIUM][ipCARBON][0] = 4.17e-17*(phycon.te/phycon.te03);*/
388
389 /* >>chng 04 jun 21 -- update this rate to that given in the UMIST database - NPA */
390 atmdat.CharExcIonOf[ipHELIUM][ipCARBON][0] = 6.3e-15*pow((phycon.te/300),0.75);
391
392 /* C+1 + He+ => C+2 + He */
393 /* >>refer C1 CT Butler, S. E., Heil, T. G., & Dalgarno, A. 1980, ApJ, 241, 442*/
394 atmdat.CharExcIonOf[ipHELIUM][ipCARBON][1] =
395 5e-20*phycon.tesqrd*sexp(0.07e-4*phycon.te)*sexp(6.29/phycon.te_eV);
396
397 /* nitrogen */
398 /* recombination */
399 /* N+2 => N+ Butler and Dalgarno 1980B
400 * ct with update
401 * >>refer N2 CT Sun, Y., Sadeghpour, H.R., Kirby, K., et al. CfA preprint 4208
402 * this agrees with exp
403 * >>refer N2 CT Fang, Z., & Kwong, V. H. S. 1997, ApJ, 474, 529 */
404 atmdat.CharExcRecTo[ipHELIUM][ipNITROGEN][1] = 0.8e-10;
405
406 /* N+3 => N+2 */
407 /* >>refer N3 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
408 atmdat.CharExcRecTo[ipHELIUM][ipNITROGEN][2] = 1.5e-10;
409
410 /* ce rate from quantal calc of ce,
411 * >>refer N4 CT Feickert, C. A., Blint, R. J., Surratt, G. T., & Watson, W.D. 1984, ApJ, 286, 371,
412 * >>refer N4 CT Rittby, M., Elander, N., Brandas, E., & Barany, A. 1984, J. Phys. B, 17, L677.
413 * CR = 1.E-9 + 8E-12 * TE10 * SQRTE */
414 atmdat.CharExcRecTo[ipHELIUM][ipNITROGEN][3] = 2e-9;
415
416 /* ionization */
417 /* N+1 + He+ => N+2 + He */
418 /* >>refer N1 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
419 atmdat.CharExcIonOf[ipHELIUM][ipNITROGEN][1] =
420 3.7e-20*phycon.tesqrd*sexp(0.063e-4*phycon.te)*sexp(1.44/phycon.te_eV);
421
422 /* oxygen */
423 /* recombination */
424 /* O+2 + He => O+1 + He+ */
425 /* >>refer O2 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
426 atmdat.CharExcRecTo[ipHELIUM][ipOXYGEN][1] = 3.2e-14*phycon.te/phycon.te05;
427 /* O+3 + He => O+2 + He+ */
428 /* >>refer O3 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
429 atmdat.CharExcRecTo[ipHELIUM][ipOXYGEN][2] = 1e-9;
430 /* O+4 + He => O+3 + He+ */
431 /* >>refer O4 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
432 atmdat.CharExcRecTo[ipHELIUM][ipOXYGEN][3] = 6e-10;
433
434 /* ionization */
435 /* O0 + He+ => O+ + He0 */
436 /* >>refer O0 CT Zhao, L. B., Stancil, P. C., Gu, J. P., et al. 2004, ApJ, 615, 1063 */
437 atmdat.CharExcIonOf[ipHELIUM][ipOXYGEN][0] =
438 4.991E-15 * pow( phycon.te / 1e4, 0.3794 )* sexp( phycon.te/1.121E6 ) +
439 2.780E-15 * pow( phycon.te / 1e4, -0.2163 )* exp( -1. * MIN2(1e7, phycon.te)/(-8.158E5) );
440
441 /* neon */
442 /* recombination */
443 /* Ne+2 + He => Ne+1 + He+ */
444 /* >>refer Ne2 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
445 atmdat.CharExcRecTo[ipHELIUM][ipNEON][1] = 1e-14;
446 /* Ne+3 + He => Ne+2 + He+ */
447 /* >>refer Ne3 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
448 atmdat.CharExcRecTo[ipHELIUM][ipNEON][2] = 1e-16*phycon.sqrte;
449 /* Ne+4 + He => Ne+3 + He+ */
450 /* >>refer Ne4 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
451 atmdat.CharExcRecTo[ipHELIUM][ipNEON][3] = 1.7e-11*phycon.sqrte;
452
453 /* magnesium */
454 /* recombination */
455 /* Mg+3 + Heo => Mg+2 */
456 /* >>refer Mg3 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
457 atmdat.CharExcRecTo[ipHELIUM][ipMAGNESIUM][2] = 7.5e-10;
458 /* Mg+4 + Heo => Mg+3 */
459 /* >>refer Mg4 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
460 atmdat.CharExcRecTo[ipHELIUM][ipMAGNESIUM][3] = 1.4e-10*phycon.te30;
461
462
463 /* silicon */
464 /* recombination */
465 /* Si+3 +He => Si+2 + He+ */
466 /* >>refer Si3 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838
467 * scale by 1.3 to bring into agreement with
468 * >>refer Si3 CT Fang, Z., & Kwong, V. H. S. 1997, ApJ, 483, 527 */
469 atmdat.CharExcRecTo[ipHELIUM][ipSILICON][2] += phycon.sqrte*phycon.te10*phycon.te10*
470 1.3*1.5e-12;
471
472 /* Si+4 + Heo => Si+3
473 * >>refer Si4 CT Opradolce, L., McCarrol, R., & Valiron, P. 1985, A&A, 148, 229 */
474 atmdat.CharExcRecTo[ipHELIUM][ipSILICON][3] = 2.54e-11*phycon.sqrte/phycon.te03/
475 phycon.te01/phycon.te01;
476
477 /* ionization */
478 /* Si0 + He+ => Si+ + He0 */
479 /* >>refer Si0 CT Prasad, S. S., & Huntress, W. T. 1980, ApJS, 43, 1-35 */
480 atmdat.CharExcIonOf[ipHELIUM][ipSILICON][0] = 3.3e-9;
481
482 /* Si+1 + He+ => Si+2 + He */
483 /* >>refer Si1 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
484 atmdat.CharExcIonOf[ipHELIUM][ipSILICON][1] =
485 1.5e-11*phycon.te20*phycon.te05*sexp(6.91/phycon.te_eV);
486
487 /* Si+2 + He+ => Si+3 + He */
488 /* >>refer Si2 CT Gargaud, M., McCarroll, R., & Valiron, P. 1982, A&AS, 45, 603 */
489 atmdat.CharExcIonOf[ipHELIUM][ipSILICON][2] =
490 1.15e-11*phycon.sqrte*sexp(8.88/phycon.te_eV);
491
492 /* sulphur */
493 /* recombination */
494 /* S+3 + Heo => S+2 */
495 /* >>refer S3 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
496 atmdat.CharExcRecTo[ipHELIUM][ipSULPHUR][2] = phycon.sqrte*1.1e-11;
497
498 /* S+4 + Heo => S+3 */
499 /* >>refer S4 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
500 /* >>chng 04 jul 01, from [ipSULPHUR][2] to [3] - bug */
501 atmdat.CharExcRecTo[ipHELIUM][ipSULPHUR][3] = 4.8e-14*phycon.te30;
502
503 /* ionization */
504 /* S+1 + He+ => S+2 + He */
505 /* >>refer S1 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
506 atmdat.CharExcIonOf[ipHELIUM][ipSULPHUR][1] =
507 4.4e-16*phycon.te*phycon.te20*sexp(0.036e-4*phycon.te)*sexp(9.2/phycon.te_eV);
508
509 /* S+2 + He+ => S+3 + He */
510 /* >>refer S2 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
511 atmdat.CharExcIonOf[ipHELIUM][ipSULPHUR][2] =
512 5.5e-18*phycon.te*phycon.sqrte*phycon.te10*sexp(0.046e-4*phycon.te)*sexp(10.5/phycon.te_eV);
513
514 /* Argon */
515 /* recombination */
516 /* >>refer Ar2 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
517 atmdat.CharExcRecTo[ipHELIUM][ipARGON][1] = 1.3e-10;
518
519 /* >>refer Ar3 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
520 atmdat.CharExcRecTo[ipHELIUM][ipARGON][2] = 1.e-14;
521
522 /* >>refer Ar4 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
523 atmdat.CharExcRecTo[ipHELIUM][ipARGON][3] = 1.6e-8/phycon.te30;
524
525 /* ionization */
526 /* Ar+1 + He+ => Ar+2 + He0 */
527 /* >>refer Ar1 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
528 atmdat.CharExcIonOf[ipHELIUM][ipARGON][1] = 1.1e-10;
529
530 TeUsed = phycon.te;
531
532 if(mole_global.lgLeidenHack)
533 {
534 /* Set all charge transfer rates equal to zero that do not appear
535 in the UMIST database. This if statement is only performed
536 if the "set UMIST rates" command is used */
537
538 atmdat.CharExcIonOf[ipHYDROGEN][ipHELIUM][0] = 0;
539 atmdat.CharExcIonOf[ipHYDROGEN][ipCARBON][0] = 0;
540 atmdat.CharExcRecTo[ipHYDROGEN][ipCARBON][0] = 0;
541
542 atmdat.CharExcRecTo[ipHELIUM][ipCARBON][0] = 0;
543 atmdat.CharExcIonOf[ipHELIUM][ipOXYGEN][0] = 0;
544 atmdat.CharExcRecTo[ipHELIUM][ipOXYGEN][0] = 0;
545 }
546
547
548 /* this is set false with the no charge transfer command */
549 if( !atmdat.lgCTOn )
550 {
551 for( nelem=0; nelem< LIMELM; ++nelem )
552 {
553 for( ion=0; ion<LIMELM; ++ion )
554 {
555 atmdat.CharExcIonOf[ipHYDROGEN][nelem][ion] = 0.;
556 atmdat.CharExcRecTo[ipHYDROGEN][nelem][ion] = 0.;
557 atmdat.CharExcIonOf[ipHELIUM][nelem][ion] = 0.;
558 atmdat.CharExcRecTo[ipHELIUM][nelem][ion] = 0.;
559 }
560 }
561 }
562 }
563
564 return;
565}
566
567/*================================================================================*
568 *================================================================================*/
570{
571 long int ion,
572 nelem;
573 double SumCTHeat_v;
574
575 DEBUG_ENTRY( "ChargTranSumHeat()" );
576
577 /* second dimension is ionization stage,
578 * 1=+1 for parent, etc
579 * third dimension is atomic weight of atom */
580
581 /* make sure data are initialized */
583
584 SumCTHeat_v = 0.;
585 /* >>chng 01 apr 25, lower limit had been 0 should have been 1 (helium) */
586 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
587 {
588 /* >>chng >>01 apr 25, loops had been to LIMELM, which may have done no harm
589 * since extra array elements were set to zero, but is incorrect since the physical
590 * limit is the number of stages of ionization */
591 int limit = MIN2(4, nelem);
592 /* this first group of lower stages of ionization have exact rate coefficients */
593 for( ion=0; ion < limit; ion++ )
594 {
595 /* CTIonData[nelem][ion][7] and CTRecombData[nelem][ion][6] are the energy deficits in eV,
596 * atmdat.CharExcIonOf[ipHYDROGEN][nelem][ion] and atmdat.CharExcIonOf[ipHYDROGEN][nelem][ion]
597 * save the rate coefficients
598 * this is sum of heat exchange in eV s^-1 cm^-3 */
599 SumCTHeat_v +=
600
601 /* heating due to ionization of heavy element, recombination of hydrogen */
602 atmdat.CharExcIonOf[ipHYDROGEN][nelem][ion]*CTIonData[nelem][ion][7]*
603 dense.xIonDense[ipHYDROGEN][1]*
604 dense.xIonDense[nelem][ion] +
605
606 /* heating due to recombination of heavy element, ionization of hydrogen */
607 atmdat.CharExcRecTo[ipHYDROGEN][nelem][ion]*CTRecombData[nelem][ion][6]*
608 //iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop()*
609 dense.xIonDense[ipHYDROGEN][0]*
610 dense.xIonDense[nelem][ion+1];
611 }
612
613 /* >>chng >>01 apr 25, following loop had been to LIMELM, change to nelem */
614 /* following do not have exact energies, so use 2.86*(Z-1) */
615 for( ion=4; ion < nelem; ion++ )
616 {
617 SumCTHeat_v +=
618 atmdat.CharExcRecTo[ipHYDROGEN][nelem][ion]* 2.86*(double)ion *
619 dense.xIonDense[ipHYDROGEN][0]*
620 dense.xIonDense[nelem][ion+1];
621 }
622 }
623
624#if 0
625 /* charge exchange with helium */
626 for( nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
627 {
628 /* >>chng >>01 apr 25, loops had been to LIMELM, which may have done no harm
629 * since extra array elements were set to zero, but is incorrect since the physical
630 * limit is the number of stages of ionization */
631 int limit = MIN2(4, nelem);
632 /* this first group of lower stages of ionization have exact rate coefficients */
633 for( ion=0; ion < limit; ion++ )
634 {
635 fixit(); // fill in these barriers
636 double barrier_rec_eV = CTIonData[nelem][ion][7];
637 double barrier_ion_eV = CTRecombData[nelem][ion][6];
638
639 /* atmdat.CharExcIonOf[ipHELIUM][nelem][ion] and atmdat.CharExcIonOf[ipHELIUM][nelem][ion]
640 * save the rate coefficients
641 * this is sum of heat exchange in eV s^-1 cm^-3 */
642 SumCTHeat_v +=
643
644 /* heating due to ionization of heavy element, recombination of helium */
645 atmdat.CharExcIonOf[ipHELIUM][nelem][ion]*barrier_rec_eV*
646 dense.xIonDense[ipHELIUM][1]*
647 dense.xIonDense[nelem][ion] +
648
649 /* heating due to recombination of heavy element, ionization of helium */
650 atmdat.CharExcRecTo[ipHELIUM][nelem][ion]*barrier_ion_eV*
651 //iso_sp[ipHE_LIKE][ipHELIUM].st[ipH1s].Pop()*
652 dense.xIonDense[ipHELIUM][0]*
653 dense.xIonDense[nelem][ion+1];
654 }
655
656 /* >>chng >>01 apr 25, following loop had been to LIMELM, change to nelem */
657 /* following do not have exact energies, so use 2.86*(Z-1) */
658 for( ion=4; ion < nelem; ion++ )
659 {
660 SumCTHeat_v +=
661 atmdat.CharExcRecTo[ipHELIUM][nelem][ion]* 2.86*(double)ion *
662 dense.xIonDense[ipHELIUM][0]*
663 dense.xIonDense[nelem][ion+1];
664 }
665 }
666#endif
667
668 /* convert from eV to ergs, HCharHeatOn usually 1, set to 0 with no CTHeat,
669 * EN1EV is ergs in 1 eV, 1.602176e-012*/
670 SumCTHeat_v *= EN1EV * atmdat.HCharHeatOn;
671
672 if( thermal.htot > 1e-35 )
673 {
674 /* remember largest fractions of heating and cooling for comment */
675 atmdat.HCharHeatMax = MAX2(atmdat.HCharHeatMax,
676 SumCTHeat_v/thermal.htot );
677
678 atmdat.HCharCoolMax = MAX2(atmdat.HCharCoolMax,
679 -SumCTHeat_v/thermal.htot);
680 }
681
682 /* debug code to print out the contributors to total CT heating */
683 {
684 enum {DEBUG_LOC=false};
685 if( DEBUG_LOC)
686 {
687# define FRAC 0.1
688 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
689 {
690 /* >>chng >>01 apr 25, loops had been to LIMELM, which may have done no harm
691 * since extra array elements were set to zero, but is incorrect since the physical
692 * limit is the number of stages of ionization */
693 int limit = MIN2(4, nelem);
694 /* this first group of lower stages of ionization have exact rate coefficients */
695 for( ion=0; ion < limit; ion++ )
696 {
697 if(
698 /* heating due to ionization of heavy element, recombination of hydrogen */
699 (atmdat.CharExcIonOf[ipHYDROGEN][nelem][ion]*CTIonData[nelem][ion][7]*
700 (double)dense.xIonDense[ipHYDROGEN][1]*(double)dense.xIonDense[nelem][ion] +
701
702 /* heating due to recombination of heavy element, ionization of hydrogen */
703 atmdat.CharExcRecTo[ipHYDROGEN][nelem][ion]*CTRecombData[nelem][ion][6]*
704 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop()*
705 (double)dense.xIonDense[nelem][ion+1])/SumCTHeat_v> FRAC )
706
707 fprintf(ioQQQ,"DEBUG ct %li %li %.3f\n", nelem, ion,
708 (atmdat.CharExcIonOf[ipHYDROGEN][nelem][ion]*CTIonData[nelem][ion][7]*
709 (double)dense.xIonDense[ipHYDROGEN][1]*(double)dense.xIonDense[nelem][ion] +
710
711 /* heating due to recombination of heavy element, ionization of hydrogen */
712 atmdat.CharExcRecTo[ipHYDROGEN][nelem][ion]*CTRecombData[nelem][ion][6]*
713 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop()*
714 (double)dense.xIonDense[nelem][ion+1]) );
715 }
716
717 for( ion=4; ion < nelem; ion++ )
718 {
719 if(
720 (atmdat.CharExcRecTo[ipHYDROGEN][nelem][ion]* 2.86*(double)ion *
721 (double)dense.xIonDense[ipHYDROGEN][0]*(double)dense.xIonDense[nelem][ion+1])/SumCTHeat_v> FRAC )
722 fprintf(ioQQQ,"DEBUG ct %li %li %.3f\n", nelem, ion,
723 (atmdat.CharExcRecTo[ipHYDROGEN][nelem][ion]* 2.86*(double)ion *
724 (double)dense.xIonDense[ipHYDROGEN][0]*(double)dense.xIonDense[nelem][ion+1]) );
725 }
726 }
727# undef FRAC
728 fprintf(ioQQQ,"DEBUT ct tot %.e3\n", SumCTHeat_v / thermal.htot );
729 }
730 }
731 return( SumCTHeat_v );
732}
733
734/*================================================================================*
735 *================================================================================*/
737 /* ion is stage of ionization on C scale, 0 for atom */
738 long int ion,
739 /* nelem is atomic number of element on C scale, 1 to 29 */
740 /* HCTIon(1,5) is C+ + H+ => C++ + H */
741 long int nelem)
742{
743 double HCTIon_v,
744 tused;
745
746 DEBUG_ENTRY( "HCTIon()" );
747 /* H charge transfer ionization, using Jim Kingdon's ctdata.for */
748
749 /* set up the rate coefficients if this is first call */
750 if( !lgCTDataDefined )
751 {
752 if( trace.lgTrace )
753 {
754 fprintf( ioQQQ," HCTIon doing 1-time init of charge transfer data\n");
755 }
756 lgCTDataDefined = true;
757 MakeHCTData();
758 }
759
760 /* check that data have been linked in,
761 * error here would mean that data below had not been loaded */
762 ASSERT( CTIonData[2][0][0] > 0. );
763
764 /* return zero for highly ionized species */
765 if( ion >= 3 )
766 {
767 HCTIon_v = 0.;
768 return( HCTIon_v );
769 }
770
771 /*begin sanity checks */
772 /* check that ionization stage is ok for this element*/
773 ASSERT( ion >= 0);
774 ASSERT( ion <= nelem );
775
776 /* now check the element is valid, this is routine HCTIon */
777 ASSERT( nelem > 0 );
778 ASSERT( nelem < LIMELM );
779
780 /* may be no entry, first coefficient is zero in this case */
781 if( CTIonData[nelem][ion][0] <= 0. )
782 {
783 HCTIon_v = 0.;
784
785 }
786
787 else
788 {
789 /* Make sure te is between temp. boundaries; set constant outside of range */
790 tused = MAX2((double)phycon.te,CTIonData[nelem][ion][4]);
791 tused = MIN2(tused,CTIonData[nelem][ion][5]);
792 tused *= 1e-4;
793
794 // make sure that the Boltzmann factor always uses the real temperature, even
795 // outside the range of validity, so that the CT rate properly goes to zero
796 // when the the gas temperature goes to zero.
797 HCTIon_v = CTIonData[nelem][ion][0]*1e-9*(pow(tused,CTIonData[nelem][ion][1]))*
798 (1. + CTIonData[nelem][ion][2]*exp(CTIonData[nelem][ion][3]*tused))*
799 exp(-CTIonData[nelem][ion][6]*1.e4/phycon.te);
800 }
801 return( HCTIon_v );
802}
803
804/*================================================================================*
805 *================================================================================*/
807 /* ion is stage of ionization on C scale, 0 for rec to atom */
808 long int ion,
809 /* nelem is atomic number of element on C scale, 1 = he up to LIMELM */
810 /* HCTRecom(1,5) would be C+2 + H => C+ + H+ */
811 long int nelem)
812{
813 double HCTRecom_v,
814 tused;
815
816 DEBUG_ENTRY( "HCTRecom()" );
817 /*
818 * H charge transfer recombination using Jim Kingdon's block ctdata.for
819 */
820
821 /* set up the rate coefficients if this is first call */
822 if( !lgCTDataDefined )
823 {
824 if( trace.lgTrace )
825 {
826 fprintf( ioQQQ," HCTIon doing 1-time init of charge transfer data\n");
827 }
828 lgCTDataDefined = true;
829 MakeHCTData();
830 }
831
832 /* this is check that data have been set up properly, will
833 * fail if arrays are not initialized properly */
834 ASSERT( CTRecombData[1][0][0] > 0. );
835
836 /* use Dalgarno estimate for highly ionized species, number reset with
837 * set charge transfer command */
838 if( ion > 3 )
839 {
840 /* >>chng 96 nov 25, added this option, default is 1.92e-9
841 * Dalgarno's charge transfer */
842 HCTRecom_v = atmdat.HCTAlex*((double)ion+1.);
843 return( HCTRecom_v );
844 }
845
846 /* check that ion stage within bound for this atom */
847 ASSERT( ion >= 0 && ion <= nelem );
848
849 /* now check the element is valid, this is routine HCTIon */
850 ASSERT( nelem > 0 && nelem < LIMELM );
851
852 tused = MAX2((double)phycon.te,CTRecombData[nelem][ion][4]);
853 tused = MIN2(tused,CTRecombData[nelem][ion][5]);
854 tused *= 1e-4;
855
856 if( tused == 0. )
857 {
858 HCTRecom_v = 0.;
859 return( HCTRecom_v );
860 }
861
862 /* the interpolation equation */
863 HCTRecom_v = CTRecombData[nelem][ion][0]*1e-9*(pow(tused,CTRecombData[nelem][ion][1]))*
864 (1. + CTRecombData[nelem][ion][2]*sexp(-CTRecombData[nelem][ion][3]*tused));
865
866 /* in sexp negative sign not typo - there are negative signs already
867 * in coefficient, and sexp has implicit negative sign */
868 return( HCTRecom_v );
869}
870
871/*================================================================================*
872 *================================================================================*/
873/*block data with Jim Kingdon's charge transfer data */
874/* >>refer H CT Kingdon, J. B., & Ferland, G. J. 1996, ApJS, 106, 205 */
875/*
876 * first dimension is atomic number of atom, 0 for H
877 * second dimension is ionization stage,
878 * 1=+0 for parent, etc
879 * third dimension is atomic number of atom
880 * second dimension is ionization stage,
881 * 1=+1 for parent, etc
882 */
883
884/* digital form of the fits to the charge transfer
885 * ionization rate coefficients
886 *
887 * Note: First parameter is in units of 1e-9!
888 * Note: Seventh parameter is in units of 1e4 K */
889
890/* digital form of the fits to the charge transfer
891 * recombination rate coefficients (total)
892 *
893 * Note: First parameter is in units of 1e-9!
894 * recombination
895 */
896
897/* holds data for charge transfer fits */
899{
900 long int i,
901 j,
902 nelem,
903 _r;
904
905 DEBUG_ENTRY( "MakeHCTData()" );
906
907 /* >>chng 01 apr 24, zero out this block, as per PvH comments that
908 * translated block data's do not fully initialize arrays */
909 /* first zero out entire arrays, since some may not have charge transfer data */
910 for( nelem=0; nelem<LIMELM; ++nelem )
911 {
912 for( i=0; i<4; ++i )
913 {
914 for( j=0; j<7; ++j )
915 {
916 CTIonData[nelem][i][j] = 0.;
917 CTRecombData[nelem][i][j] = 0.;
918 }
919 CTIonData[nelem][i][7] = 0.;
920 }
921 }
922
923 /*
924 * following are coefficients for charge transfer ionization,
925 * H+ + A => H + A+
926 */
927 /* Lithium +0 */
928 { static double _itmp0[] = {2.84e-3 , 1.99 , 375.54 , -54.07 , 1e2 , 1e4 , 0.,
929 -10.};
930
931 for( i=1, _r = 0; i <= 8; i++ )
932 {
933 CTIonData[2][0][i-1] = _itmp0[_r++];
934 }
935 }
936
937 /* C+0 ionization */
938 { static double _itmp1[] = {1.07e-6 , 3.15 , 176.43 , -4.29 , 1e3 , 1e5 , 0. ,2.34};
939 for( i=1, _r = 0; i <= 8; i++ )
940 {
941 CTIonData[5][0][i-1] = _itmp1[_r++];
942 }
943 }
944 { static double _itmp2[] = {4.55e-3,-0.29,-0.92,-8.38,1e2,5e4,
945 1.086,-0.94};
946 for( i=1, _r = 0; i <= 8; i++ )
947 {
948 CTIonData[6][0][i-1] = _itmp2[_r++];
949 }
950 }
951 /* oxygen */
952 { static double _itmp3[] = {7.40e-2,0.47,24.37,-0.74,1e1,1e4,
953 0.023,-0.02};
954 for( i=1, _r = 0; i <= 8; i++ )
955 {
956 CTIonData[7][0][i-1] = _itmp3[_r++];
957 }
958 }
959 { static double _itmp4[] = {3.34e-6,9.31,2632.31,-3.04,1e3,
960 2e4,0.0,-1.74};
961 for( i=1, _r = 0; i <= 8; i++ )
962 {
963 CTIonData[10][0][i-1] = _itmp4[_r++];
964 }
965 }
966 { static double _itmp5[] = {9.76e-3,3.14,55.54,-1.12,5e3,3e4,
967 0.0,1.52};
968 for( i=1, _r = 0; i <= 8; i++ )
969 {
970 CTIonData[11][0][i-1] = _itmp5[_r++];
971 }
972 }
973 { static double _itmp6[] = {7.60e-5,0.00,-1.97,-4.32,1e4,3e5,
974 1.670,-1.44};
975 for( i=1, _r = 0; i <= 8; i++ )
976 {
977 CTIonData[11][1][i-1] = _itmp6[_r++];
978 }
979 }
980 { static double _itmp7[] = {0.92,1.15,0.80,-0.24,1e3,2e5,0.0,
981 0.12};
982 for( i=1, _r = 0; i <= 8; i++ )
983 {
984 CTIonData[13][0][i-1] = _itmp7[_r++];
985 }
986 }
987 /* Si+1 ionization */
988 { static double _itmp8[] = {2.26 , 7.36e-2 , -0.43 , -0.11 , 2e3 , 1e5 , 3.031
989 ,-2.72};
990 for( i=1, _r = 0; i <= 8; i++ )
991 {
992 CTIonData[13][1][i-1] = _itmp8[_r++];
993 }
994 }
995 { static double _itmp9[] = {1.00e-5,0.00,0.00,0.00,1e3,1e4,
996 0.0,-3.24};
997 for( i=1, _r = 0; i <= 8; i++ )
998 {
999 CTIonData[15][0][i-1] = _itmp9[_r++];
1000 }
1001 }
1002 { static double _itmp10[] = {4.39,0.61,-0.89,-3.56,1e3,3e4,
1003 3.349,-2.89};
1004 for( i=1, _r = 0; i <= 8; i++ )
1005 {
1006 CTIonData[23][1][i-1] = _itmp10[_r++];
1007 }
1008 }
1009 { static double _itmp11[] = {2.83e-1,6.80e-3,6.44e-2,-9.70,
1010 1e3,3e4,2.368,-2.04};
1011 for( i=1, _r = 0; i <= 8; i++ )
1012 {
1013 CTIonData[24][1][i-1] = _itmp11[_r++];
1014 }
1015 }
1016 { static double _itmp12[] = {2.10,7.72e-2,-0.41,-7.31,1e4,1e5,
1017 3.005,-2.56};
1018 for( i=1, _r = 0; i <= 8; i++ )
1019 {
1020 CTIonData[25][1][i-1] = _itmp12[_r++];
1021 }
1022 }
1023 { static double _itmp13[] = {1.20e-2,3.49,24.41,-1.26,1e3,3e4,
1024 4.044,-3.49};
1025 for( i=1, _r = 0; i <= 8; i++ )
1026 {
1027 CTIonData[26][1][i-1] = _itmp13[_r++];
1028 }
1029 }
1030 /* CT recombination, A+n + H => A+n-1 + H+ */
1031 /* >>chng 01 may 03, first coefficient multiplied by 0.25, as per comment in
1032 * >>refer Li CT Stancil, P. C., & Zygelman, B. 1996, ApJ, 472, 102
1033 * which corrected the error in
1034 * >>refer He CT Zygelman, B., Dalgarno, A., Kimura, M., & Lane, N. F. 1989, Phys. Rev. A, 40, 2340
1035 * this was used in the original Kingdon & Ferland paper so no correction required
1036 * >>chng 04 apr 27, He was in error above as well, factor of 4, noted in
1037 * >>refer He CT Stancil, P. C., Lepp, S., & Dalgarno, A. 1998, ApJ, 509, 1
1038 */
1039 { static double _itmp14[] = {/*7.47e-6*/1.87e-6,2.06,9.93,-3.89,6e3,1e5,
1040 10.99};
1041 for( i=1, _r = 0; i <= 7; i++ )
1042 {
1043 CTRecombData[1][0][i-1] = _itmp14[_r++];
1044 }
1045 }
1046 { static double _itmp15[] = {1.00e-5,0.,0.,0.,1e3,1e7,-40.81};
1047 for( i=1, _r = 0; i <= 7; i++ )
1048 {
1049 CTRecombData[1][1][i-1] = _itmp15[_r++];
1050 }
1051 }
1052 for( i=1; i <= 7; i++ )
1053 {
1054 CTRecombData[2][0][i-1] = 0.f;
1055 }
1056 { static double _itmp16[] = {1.26,0.96,3.02,-0.65,1e3,3e4,3.02};
1057 for( i=1, _r = 0; i <= 7; i++ )
1058 {
1059 CTRecombData[2][1][i-1] = _itmp16[_r++];
1060 }
1061 }
1062 { static double _itmp17[] = {1.00e-5,0.,0.,0.,2e3,5e4,-108.83};
1063 for( i=1, _r = 0; i <= 7; i++ )
1064 {
1065 CTRecombData[2][2][i-1] = _itmp17[_r++];
1066 }
1067 }
1068 for( i=1; i <= 7; i++ )
1069 {
1070 CTRecombData[3][0][i-1] = 0.f;
1071 }
1072 { static double _itmp18[] = {1.00e-5,0.,0.,0.,2e3,5e4,-4.61};
1073 for( i=1, _r = 0; i <= 7; i++ )
1074 {
1075 CTRecombData[3][1][i-1] = _itmp18[_r++];
1076 }
1077 }
1078 { static double _itmp19[] = {1.00e-5,0.,0.,0.,2e3,5e4,-140.26};
1079 for( i=1, _r = 0; i <= 7; i++ )
1080 {
1081 CTRecombData[3][2][i-1] = _itmp19[_r++];
1082 }
1083 }
1084 { static double _itmp20[] = {5.17,0.82,-0.69,-1.12,2e3,5e4,
1085 10.59};
1086 for( i=1, _r = 0; i <= 7; i++ )
1087 {
1088 CTRecombData[3][3][i-1] = _itmp20[_r++];
1089 }
1090 }
1091 for( i=1; i <= 7; i++ )
1092 {
1093 CTRecombData[4][0][i-1] = 0.f;
1094 }
1095 { static double _itmp21[] = {2.00e-2,0.,0.,0.,1e3,1e9,2.46};
1096 for( i=1, _r = 0; i <= 7; i++ )
1097 {
1098 CTRecombData[4][1][i-1] = _itmp21[_r++];
1099 }
1100 }
1101 { static double _itmp22[] = {1.00e-5,0.,0.,0.,2e3,5e4,-24.33};
1102 for( i=1, _r = 0; i <= 7; i++ )
1103 {
1104 CTRecombData[4][2][i-1] = _itmp22[_r++];
1105 }
1106 }
1107 /* B+4 recombinatino */
1108 { static double _itmp23[] = {2.74 , 0.93 , -0.61 , -1.13 , 2e3 , 5e4 ,
1109 11.};
1110 for( i=1, _r = 0; i <= 7; i++ )
1111 {
1112 CTRecombData[4][3][i-1] = _itmp23[_r++];
1113 }
1114 }
1115 /* C+1 recombinatino */
1116 { static double _itmp24[] = {4.88e-7 , 3.25 , -1.12 , -0.21 , 5.5e3 , 1e5 ,
1117 -2.34};
1118 for( i=1, _r = 0; i <= 7; i++ )
1119 {
1120 CTRecombData[5][0][i-1] = _itmp24[_r++];
1121 }
1122 }
1123 { static double _itmp25[] = {1.67e-4,2.79,304.72,-4.07,5e3,
1124 5e4,4.01};
1125 for( i=1, _r = 0; i <= 7; i++ )
1126 {
1127 CTRecombData[5][1][i-1] = _itmp25[_r++];
1128 }
1129 }
1130 { static double _itmp26[] = {3.25,0.21,0.19,-3.29,1e3,1e5,5.73};
1131 for( i=1, _r = 0; i <= 7; i++ )
1132 {
1133 CTRecombData[5][2][i-1] = _itmp26[_r++];
1134 }
1135 }
1136 { static double _itmp27[] = {332.46,-0.11,-9.95e-1,-1.58e-3,
1137 1e1,1e5,11.30};
1138 for( i=1, _r = 0; i <= 7; i++ )
1139 {
1140 CTRecombData[5][3][i-1] = _itmp27[_r++];
1141 }
1142 }
1143 { static double _itmp28[] = {1.01e-3,-0.29,-0.92,-8.38,1e2,
1144 5e4,0.94};
1145 for( i=1, _r = 0; i <= 7; i++ )
1146 {
1147 CTRecombData[6][0][i-1] = _itmp28[_r++];
1148 }
1149 }
1150 { static double _itmp29[] = {3.05e-1,0.60,2.65,-0.93,1e3,1e5,
1151 4.56};
1152 for( i=1, _r = 0; i <= 7; i++ )
1153 {
1154 CTRecombData[6][1][i-1] = _itmp29[_r++];
1155 }
1156 }
1157 { static double _itmp30[] = {4.54,0.57,-0.65,-0.89,1e1,1e5,
1158 6.40};
1159 for( i=1, _r = 0; i <= 7; i++ )
1160 {
1161 CTRecombData[6][2][i-1] = _itmp30[_r++];
1162 }
1163 }
1164 /* N+4 recombination */
1165 { static double _itmp31[] = { 2.95 , 0.55 , -0.39 , -1.07 , 1e3 , 1e6 ,
1166 11.};
1167 for( i=1, _r = 0; i <= 7; i++ )
1168 {
1169 CTRecombData[6][3][i-1] = _itmp31[_r++];
1170 }
1171 }
1172 { static double _itmp32[] = {1.04,3.15e-2,-0.61,-9.73,1e1,1e4,
1173 0.02};
1174 for( i=1, _r = 0; i <= 7; i++ )
1175 {
1176 CTRecombData[7][0][i-1] = _itmp32[_r++];
1177 }
1178 }
1179 { static double _itmp33[] = {1.04,0.27,2.02,-5.92,1e2,1e5,6.65};
1180 for( i=1, _r = 0; i <= 7; i++ )
1181 {
1182 CTRecombData[7][1][i-1] = _itmp33[_r++];
1183 }
1184 }
1185 { static double _itmp34[] = {3.98,0.26,0.56,-2.62,1e3,5e4,5.};
1186 for( i=1, _r = 0; i <= 7; i++ )
1187 {
1188 CTRecombData[7][2][i-1] = _itmp34[_r++];
1189 }
1190 }
1191 { static double _itmp35[] = {2.52e-1,0.63,2.08,-4.16,1e3,3e4,
1192 8.47};
1193 for( i=1, _r = 0; i <= 7; i++ )
1194 {
1195 CTRecombData[7][3][i-1] = _itmp35[_r++];
1196 }
1197 }
1198 for( i=1; i <= 7; i++ )
1199 {
1200 CTRecombData[8][0][i-1] = 0.f;
1201 }
1202 { static double _itmp36[] = {1.00e-5,0.,0.,0.,2e3,5e4,-21.37};
1203 for( i=1, _r = 0; i <= 7; i++ )
1204 {
1205 CTRecombData[8][1][i-1] = _itmp36[_r++];
1206 }
1207 }
1208 { static double _itmp37[] = {9.86,0.29,-0.21,-1.15,2e3,5e4,
1209 5.6};
1210 for( i=1, _r = 0; i <= 7; i++ )
1211 {
1212 CTRecombData[8][2][i-1] = _itmp37[_r++];
1213 }
1214 }
1215 { static double _itmp38[] = {7.15e-1,1.21,-0.70,-0.85,2e3,5e4,
1216 11.8};
1217 for( i=1, _r = 0; i <= 7; i++ )
1218 {
1219 CTRecombData[8][3][i-1] = _itmp38[_r++];
1220 }
1221 }
1222 for( i=1; i <= 7; i++ )
1223 {
1224 CTRecombData[9][0][i-1] = 0.f;
1225 }
1226 { static double _itmp39[] = {1.00e-5,0.,0.,0.,5e3,5e4,-27.36};
1227 for( i=1, _r = 0; i <= 7; i++ )
1228 {
1229 CTRecombData[9][1][i-1] = _itmp39[_r++];
1230 }
1231 }
1232 { static double _itmp40[] = {14.73,4.52e-2,-0.84,-0.31,5e3,
1233 5e4,5.82};
1234 for( i=1, _r = 0; i <= 7; i++ )
1235 {
1236 CTRecombData[9][2][i-1] = _itmp40[_r++];
1237 }
1238 }
1239 { static double _itmp41[] = {6.47,0.54,3.59,-5.22,1e3,3e4,8.60};
1240 for( i=1, _r = 0; i <= 7; i++ )
1241 {
1242 CTRecombData[9][3][i-1] = _itmp41[_r++];
1243 }
1244 }
1245 for( i=1; i <= 7; i++ )
1246 {
1247 CTRecombData[10][0][i-1] = 0.f;
1248 }
1249 { static double _itmp42[] = {1.00e-5,0.,0.,0.,2e3,5e4,-33.68};
1250 for( i=1, _r = 0; i <= 7; i++ )
1251 {
1252 CTRecombData[10][1][i-1] = _itmp42[_r++];
1253 }
1254 }
1255 { static double _itmp43[] = {1.33,1.15,1.20,-0.32,2e3,5e4,6.25};
1256 for( i=1, _r = 0; i <= 7; i++ )
1257 {
1258 CTRecombData[10][2][i-1] = _itmp43[_r++];
1259 }
1260 }
1261 { static double _itmp44[] = {1.01e-1,1.34,10.05,-6.41,2e3,5e4,
1262 11.};
1263 for( i=1, _r = 0; i <= 7; i++ )
1264 {
1265 CTRecombData[10][3][i-1] = _itmp44[_r++];
1266 }
1267 }
1268 for( i=1; i <= 7; i++ )
1269 {
1270 CTRecombData[11][0][i-1] = 0.f;
1271 }
1272 { static double _itmp45[] = {8.58e-5,2.49e-3,2.93e-2,-4.33,
1273 1e3,3e4,1.44};
1274 for( i=1, _r = 0; i <= 7; i++ )
1275 {
1276 CTRecombData[11][1][i-1] = _itmp45[_r++];
1277 }
1278 }
1279 { static double _itmp46[] = {6.49,0.53,2.82,-7.63,1e3,3e4,5.73};
1280 for( i=1, _r = 0; i <= 7; i++ )
1281 {
1282 CTRecombData[11][2][i-1] = _itmp46[_r++];
1283 }
1284 }
1285 { static double _itmp47[] = {6.36,0.55,3.86,-5.19,1e3,3e4,8.60};
1286 for( i=1, _r = 0; i <= 7; i++ )
1287 {
1288 CTRecombData[11][3][i-1] = _itmp47[_r++];
1289 }
1290 }
1291 for( i=1; i <= 7; i++ )
1292 {
1293 CTRecombData[12][0][i-1] = 0.f;
1294 }
1295 { static double _itmp48[] = {1.00e-5,0.,0.,0.,1e3,3e4,-5.23};
1296 for( i=1, _r = 0; i <= 7; i++ )
1297 {
1298 CTRecombData[12][1][i-1] = _itmp48[_r++];
1299 }
1300 }
1301 { static double _itmp49[] = {7.11e-5,4.12,1.72e4,-22.24,1e3,
1302 3e4,8.17};
1303 for( i=1, _r = 0; i <= 7; i++ )
1304 {
1305 CTRecombData[12][2][i-1] = _itmp49[_r++];
1306 }
1307 }
1308 { static double _itmp50[] = {7.52e-1,0.77,6.24,-5.67,1e3,3e4,
1309 8.};
1310 for( i=1, _r = 0; i <= 7; i++ )
1311 {
1312 CTRecombData[12][3][i-1] = _itmp50[_r++];
1313 }
1314 }
1315 for( i=1; i <= 7; i++ )
1316 {
1317 CTRecombData[13][0][i-1] = 0.f;
1318 }
1319 /* Si+2 recombination */
1320 { static double _itmp51[] = {6.77 , 7.36e-2 , -0.43 , -0.11 , 5e2 , 1e5 ,
1321 2.72};
1322 for( i=1, _r = 0; i <= 7; i++ )
1323 {
1324 CTRecombData[13][1][i-1] = _itmp51[_r++];
1325 }
1326 }
1327 { static double _itmp52[] = {4.90e-1,-8.74e-2,-0.36,-0.79,1e3,
1328 3e4,4.23};
1329 for( i=1, _r = 0; i <= 7; i++ )
1330 {
1331 CTRecombData[13][2][i-1] = _itmp52[_r++];
1332 }
1333 }
1334 { static double _itmp53[] = {7.58,0.37,1.06,-4.09,1e3,5e4,7.49};
1335 for( i=1, _r = 0; i <= 7; i++ )
1336 {
1337 CTRecombData[13][3][i-1] = _itmp53[_r++];
1338 }
1339 }
1340 for( i=1; i <= 7; i++ )
1341 {
1342 CTRecombData[14][0][i-1] = 0.f;
1343 }
1344 { static double _itmp54[] = {1.74e-4,3.84,36.06,-0.97,1e3,3e4,
1345 3.45};
1346 for( i=1, _r = 0; i <= 7; i++ )
1347 {
1348 CTRecombData[14][1][i-1] = _itmp54[_r++];
1349 }
1350 }
1351 { static double _itmp55[] = {9.46e-2,-5.58e-2,0.77,-6.43,1e3,
1352 3e4,7.29};
1353 for( i=1, _r = 0; i <= 7; i++ )
1354 {
1355 CTRecombData[14][2][i-1] = _itmp55[_r++];
1356 }
1357 }
1358 { static double _itmp56[] = {5.37,0.47,2.21,-8.52,1e3,3e4,9.71};
1359 for( i=1, _r = 0; i <= 7; i++ )
1360 {
1361 CTRecombData[14][3][i-1] = _itmp56[_r++];
1362 }
1363 }
1364 { static double _itmp57[] = {3.82e-7,11.10,2.57e4,-8.22,1e3,
1365 1e4,-3.24};
1366 for( i=1, _r = 0; i <= 7; i++ )
1367 {
1368 CTRecombData[15][0][i-1] = _itmp57[_r++];
1369 }
1370 }
1371 { static double _itmp58[] = {1.00e-5,0.,0.,0.,1e3,3e4,-9.73};
1372 for( i=1, _r = 0; i <= 7; i++ )
1373 {
1374 CTRecombData[15][1][i-1] = _itmp58[_r++];
1375 }
1376 }
1377 { static double _itmp59[] = {2.29,4.02e-2,1.59,-6.06,1e3,3e4,
1378 5.73};
1379 for( i=1, _r = 0; i <= 7; i++ )
1380 {
1381 CTRecombData[15][2][i-1] = _itmp59[_r++];
1382 }
1383 }
1384 { static double _itmp60[] = {6.44,0.13,2.69,-5.69,1e3,3e4,8.60};
1385 for( i=1, _r = 0; i <= 7; i++ )
1386 {
1387 CTRecombData[15][3][i-1] = _itmp60[_r++];
1388 }
1389 }
1390 for( i=1; i <= 7; i++ )
1391 {
1392 CTRecombData[16][0][i-1] = 0.f;
1393 }
1394 { static double _itmp61[] = {1.00e-5,0.,0.,0.,1e3,3e4,-10.21};
1395 for( i=1, _r = 0; i <= 7; i++ )
1396 {
1397 CTRecombData[16][1][i-1] = _itmp61[_r++];
1398 }
1399 }
1400 { static double _itmp62[] = {1.88,0.32,1.77,-5.70,1e3,3e4,8.};
1401 for( i=1, _r = 0; i <= 7; i++ )
1402 {
1403 CTRecombData[16][2][i-1] = _itmp62[_r++];
1404 }
1405 }
1406 { static double _itmp63[] = {7.27,0.29,1.04,-10.14,1e3,3e4,
1407 9.};
1408 for( i=1, _r = 0; i <= 7; i++ )
1409 {
1410 CTRecombData[16][3][i-1] = _itmp63[_r++];
1411 }
1412 }
1413 for( i=1; i <= 7; i++ )
1414 {
1415 CTRecombData[17][0][i-1] = 0.f;
1416 }
1417 { static double _itmp64[] = {1.00e-5,0.,0.,0.,1e3,3e4,-14.03};
1418 for( i=1, _r = 0; i <= 7; i++ )
1419 {
1420 CTRecombData[17][1][i-1] = _itmp64[_r++];
1421 }
1422 }
1423 { static double _itmp65[] = {4.57,0.27,-0.18,-1.57,1e3,3e4,
1424 5.73};
1425 for( i=1, _r = 0; i <= 7; i++ )
1426 {
1427 CTRecombData[17][2][i-1] = _itmp65[_r++];
1428 }
1429 }
1430 { static double _itmp66[] = {6.37,0.85,10.21,-6.22,1e3,3e4,
1431 8.60};
1432 for( i=1, _r = 0; i <= 7; i++ )
1433 {
1434 CTRecombData[17][3][i-1] = _itmp66[_r++];
1435 }
1436 }
1437 for( i=1; i <= 7; i++ )
1438 {
1439 CTRecombData[18][0][i-1] = 0.f;
1440 }
1441 { static double _itmp67[] = {1.00e-5,0.,0.,0.,1e3,3e4,-18.02};
1442 for( i=1, _r = 0; i <= 7; i++ )
1443 {
1444 CTRecombData[18][1][i-1] = _itmp67[_r++];
1445 }
1446 }
1447 { static double _itmp68[] = {4.76,0.44,-0.56,-0.88,1e3,3e4,
1448 6.};
1449 for( i=1, _r = 0; i <= 7; i++ )
1450 {
1451 CTRecombData[18][2][i-1] = _itmp68[_r++];
1452 }
1453 }
1454 { static double _itmp69[] = {1.00e-5,0.,0.,0.,1e3,3e4,-47.3};
1455 for( i=1, _r = 0; i <= 7; i++ )
1456 {
1457 CTRecombData[18][3][i-1] = _itmp69[_r++];
1458 }
1459 }
1460 for( i=1; i <= 7; i++ )
1461 {
1462 CTRecombData[19][0][i-1] = 0.f;
1463 }
1464 { static double _itmp70[] = {0.,0.,0.,0.,1e1,1e9,0.};
1465 for( i=1, _r = 0; i <= 7; i++ )
1466 {
1467 CTRecombData[19][1][i-1] = _itmp70[_r++];
1468 }
1469 }
1470 { static double _itmp71[] = {3.17e-2,2.12,12.06,-0.40,1e3,3e4,
1471 6.6};
1472 for( i=1, _r = 0; i <= 7; i++ )
1473 {
1474 CTRecombData[19][2][i-1] = _itmp71[_r++];
1475 }
1476 }
1477 { static double _itmp72[] = {2.68,0.69,-0.68,-4.47,1e3,3e4,
1478 9.9};
1479 for( i=1, _r = 0; i <= 7; i++ )
1480 {
1481 CTRecombData[19][3][i-1] = _itmp72[_r++];
1482 }
1483 }
1484 for( i=1; i <= 7; i++ )
1485 {
1486 CTRecombData[20][0][i-1] = 0.f;
1487 }
1488 { static double _itmp73[] = {0.,0.,0.,0.,1e1,1e9,0.};
1489 for( i=1, _r = 0; i <= 7; i++ )
1490 {
1491 CTRecombData[20][1][i-1] = _itmp73[_r++];
1492 }
1493 }
1494 { static double _itmp74[] = {7.22e-3,2.34,411.50,-13.24,1e3,
1495 3e4,3.5};
1496 for( i=1, _r = 0; i <= 7; i++ )
1497 {
1498 CTRecombData[20][2][i-1] = _itmp74[_r++];
1499 }
1500 }
1501 { static double _itmp75[] = {1.20e-1,1.48,4.00,-9.33,1e3,3e4,
1502 10.61};
1503 for( i=1, _r = 0; i <= 7; i++ )
1504 {
1505 CTRecombData[20][3][i-1] = _itmp75[_r++];
1506 }
1507 }
1508 for( i=1; i <= 7; i++ )
1509 {
1510 CTRecombData[21][0][i-1] = 0.f;
1511 }
1512 { static double _itmp76[] = {0.,0.,0.,0.,1e1,1e9,0.};
1513 for( i=1, _r = 0; i <= 7; i++ )
1514 {
1515 CTRecombData[21][1][i-1] = _itmp76[_r++];
1516 }
1517 }
1518 { static double _itmp77[] = {6.34e-1,6.87e-3,0.18,-8.04,1e3,
1519 3e4,4.3};
1520 for( i=1, _r = 0; i <= 7; i++ )
1521 {
1522 CTRecombData[21][2][i-1] = _itmp77[_r++];
1523 }
1524 }
1525 { static double _itmp78[] = {4.37e-3,1.25,40.02,-8.05,1e3,3e4,
1526 5.3};
1527 for( i=1, _r = 0; i <= 7; i++ )
1528 {
1529 CTRecombData[21][3][i-1] = _itmp78[_r++];
1530 }
1531 }
1532 for( i=1; i <= 7; i++ )
1533 {
1534 CTRecombData[22][0][i-1] = 0.f;
1535 }
1536 { static double _itmp79[] = {1.00e-5,0.,0.,0.,1e3,3e4,-1.05};
1537 for( i=1, _r = 0; i <= 7; i++ )
1538 {
1539 CTRecombData[22][1][i-1] = _itmp79[_r++];
1540 }
1541 }
1542 { static double _itmp80[] = {5.12,-2.18e-2,-0.24,-0.83,1e3,
1543 3e4,4.7};
1544 for( i=1, _r = 0; i <= 7; i++ )
1545 {
1546 CTRecombData[22][2][i-1] = _itmp80[_r++];
1547 }
1548 }
1549 { static double _itmp81[] = {1.96e-1,-8.53e-3,0.28,-6.46,1e3,
1550 3e4,6.2};
1551 for( i=1, _r = 0; i <= 7; i++ )
1552 {
1553 CTRecombData[22][3][i-1] = _itmp81[_r++];
1554 }
1555 }
1556 for( i=1; i <= 7; i++ )
1557 {
1558 CTRecombData[23][0][i-1] = 0.f;
1559 }
1560 { static double _itmp82[] = {5.27e-1,0.61,-0.89,-3.56,1e3,3e4,
1561 2.89};
1562 for( i=1, _r = 0; i <= 7; i++ )
1563 {
1564 CTRecombData[23][1][i-1] = _itmp82[_r++];
1565 }
1566 }
1567 { static double _itmp83[] = {10.90,0.24,0.26,-11.94,1e3,3e4,
1568 5.4};
1569 for( i=1, _r = 0; i <= 7; i++ )
1570 {
1571 CTRecombData[23][2][i-1] = _itmp83[_r++];
1572 }
1573 }
1574 { static double _itmp84[] = {1.18,0.20,0.77,-7.09,1e3,3e4,6.6};
1575 for( i=1, _r = 0; i <= 7; i++ )
1576 {
1577 CTRecombData[23][3][i-1] = _itmp84[_r++];
1578 }
1579 }
1580 for( i=1; i <= 7; i++ )
1581 {
1582 CTRecombData[24][0][i-1] = 0.f;
1583 }
1584 { static double _itmp85[] = {1.65e-1,6.80e-3,6.44e-2,-9.70,
1585 1e3,3e4,2.04};
1586 for( i=1, _r = 0; i <= 7; i++ )
1587 {
1588 CTRecombData[24][1][i-1] = _itmp85[_r++];
1589 }
1590 }
1591 { static double _itmp86[] = {14.20,0.34,-0.41,-1.19,1e3,3e4,
1592 6.};
1593 for( i=1, _r = 0; i <= 7; i++ )
1594 {
1595 CTRecombData[24][2][i-1] = _itmp86[_r++];
1596 }
1597 }
1598 { static double _itmp87[] = {4.43e-1,0.91,10.76,-7.49,1e3,3e4,
1599 7.};
1600 for( i=1, _r = 0; i <= 7; i++ )
1601 {
1602 CTRecombData[24][3][i-1] = _itmp87[_r++];
1603 }
1604 }
1605 for( i=1; i <= 7; i++ )
1606 {
1607 CTRecombData[25][0][i-1] = 0.f;
1608 }
1609 { static double _itmp88[] = {1.26,7.72e-2,-0.41,-7.31,1e3,1e5,
1610 2.56};
1611 for( i=1, _r = 0; i <= 7; i++ )
1612 {
1613 CTRecombData[25][1][i-1] = _itmp88[_r++];
1614 }
1615 }
1616 { static double _itmp89[] = {3.42,0.51,-2.06,-8.99,1e3,1e5,
1617 6.3};
1618 for( i=1, _r = 0; i <= 7; i++ )
1619 {
1620 CTRecombData[25][2][i-1] = _itmp89[_r++];
1621 }
1622 }
1623 { static double _itmp90[] = {14.60,3.57e-2,-0.92,-0.37,1e3,
1624 3e4,10.};
1625 for( i=1, _r = 0; i <= 7; i++ )
1626 {
1627 CTRecombData[25][3][i-1] = _itmp90[_r++];
1628 }
1629 }
1630 for( i=1; i <= 7; i++ )
1631 {
1632 CTRecombData[26][0][i-1] = 0.f;
1633 }
1634 { static double _itmp91[] = {5.30,0.24,-0.91,-0.47,1e3,3e4,
1635 2.9};
1636 for( i=1, _r = 0; i <= 7; i++ )
1637 {
1638 CTRecombData[26][1][i-1] = _itmp91[_r++];
1639 }
1640 }
1641 { static double _itmp92[] = {3.26,0.87,2.85,-9.23,1e3,3e4,6.};
1642 for( i=1, _r = 0; i <= 7; i++ )
1643 {
1644 CTRecombData[26][2][i-1] = _itmp92[_r++];
1645 }
1646 }
1647 { static double _itmp93[] = {1.03,0.58,-0.89,-0.66,1e3,3e4,
1648 10.51};
1649 for( i=1, _r = 0; i <= 7; i++ )
1650 {
1651 CTRecombData[26][3][i-1] = _itmp93[_r++];
1652 }
1653 }
1654 for( i=1; i <= 7; i++ )
1655 {
1656 CTRecombData[27][0][i-1] = 0.f;
1657 }
1658 { static double _itmp94[] = {1.05,1.28,6.54,-1.81,1e3,1e5,3.0};
1659 for( i=1, _r = 0; i <= 7; i++ )
1660 {
1661 CTRecombData[27][1][i-1] = _itmp94[_r++];
1662 }
1663 }
1664 { static double _itmp95[] = {9.73,0.35,0.90,-5.33,1e3,3e4,5.2};
1665 for( i=1, _r = 0; i <= 7; i++ )
1666 {
1667 CTRecombData[27][2][i-1] = _itmp95[_r++];
1668 }
1669 }
1670 { static double _itmp96[] = {6.14,0.25,-0.91,-0.42,1e3,3e4,
1671 10.};
1672 for( i=1, _r = 0; i <= 7; i++ )
1673 {
1674 CTRecombData[27][3][i-1] = _itmp96[_r++];
1675 }
1676 }
1677 for( i=1; i <= 7; i++ )
1678 {
1679 CTRecombData[28][0][i-1] = 0.f;
1680 }
1681 { static double _itmp97[] = {1.47e-3,3.51,23.91,-0.93,1e3,3e4,
1682 3.44};
1683 for( i=1, _r = 0; i <= 7; i++ )
1684 {
1685 CTRecombData[28][1][i-1] = _itmp97[_r++];
1686 }
1687 }
1688 { static double _itmp98[] = {9.26,0.37,0.40,-10.73,1e3,3e4,
1689 5.6};
1690 for( i=1, _r = 0; i <= 7; i++ )
1691 {
1692 CTRecombData[28][2][i-1] = _itmp98[_r++];
1693 }
1694 }
1695 { static double _itmp99[] = {11.59,0.20,0.80,-6.62,1e3,3e4,
1696 9.};
1697 for( i=1, _r = 0; i <= 7; i++ )
1698 {
1699 CTRecombData[28][3][i-1] = _itmp99[_r++];
1700 }
1701 }
1702 for( i=1; i <= 7; i++ )
1703 {
1704 CTRecombData[29][0][i-1] = 0.f;
1705 }
1706 { static double _itmp100[] = {1.00e-5,0.,0.,0.,1e3,3e4,-4.37};
1707 for( i=1, _r = 0; i <= 7; i++ )
1708 {
1709 CTRecombData[29][1][i-1] = _itmp100[_r++];
1710 }
1711 }
1712 { static double _itmp101[] = {6.96e-4,4.24,26.06,-1.24,1e3,
1713 3e4,7.8};
1714 for( i=1, _r = 0; i <= 7; i++ )
1715 {
1716 CTRecombData[29][2][i-1] = _itmp101[_r++];
1717 }
1718 }
1719 { static double _itmp102[] = {1.33e-2,1.56,-0.92,-1.20,1e3,
1720 3e4,11.73};
1721 for( i=1, _r = 0; i <= 7; i++ )
1722 {
1723 CTRecombData[29][3][i-1] = _itmp102[_r++];
1724 }
1725 }
1726}
1727
1728/* ChargTranPun, save charge transfer coefficient */
1729void ChargTranPun( FILE* ipPnunit , char* chSave )
1730{
1731 long int j, jj;
1732 /* restore temp when done with this routine */
1733 double TempSave = phycon.te;
1734
1735 DEBUG_ENTRY( "ChargTranPun()" );
1736
1737 /* this is the usual charge transfer option */
1738 if( strcmp( chSave,"CHAR") == 0 )
1739 {
1740 /* charge exchange rate coefficients, entered with
1741 * save charge transfer command. Queries Jim Kingdon's
1742 * CT fits and routines to get H - He and higher CT rates
1743 * rates are evaluated at the current temperature, which can
1744 * be specified with the constant temperature command */
1745 /* first group is charge transfer recombination,
1746 * the process A+x + H => A+x-1 + H+ */
1747 fprintf( ipPnunit, "#element\tion\n");
1748 for( j=1; j < LIMELM; j++ )
1749 {
1750 /* first number is atomic number, so add 1 to get onto physical scale */
1751 /* >>chng 00 may 09, caught by Jon Slavin */
1752 fprintf( ipPnunit, "%s\t", elementnames.chElementSym[j] );
1753 /*fprintf( ipPnunit, "%3ld\t", j+1 );*/
1754
1755 for( jj=0; jj < j; jj++ )
1756 {
1757 fprintf( ipPnunit, "%.2e\t",
1758 HCTRecom(jj,j) );
1759 }
1760 fprintf( ipPnunit, "\n" );
1761 }
1762
1763 /* second group is charge transfer ionization,
1764 * the process A+x + H+ => A+x+1 + H */
1765 fprintf( ipPnunit, "\n#ionization rates, atomic number\n");
1766 for( j=1; j < LIMELM; j++ )
1767 {
1768 fprintf( ipPnunit, "%s\t", elementnames.chElementSym[j] );
1769 for( jj=0; jj < j; jj++ )
1770 {
1771 fprintf( ipPnunit, "%.2e\t",
1772 HCTIon(jj,j) );
1773 }
1774 fprintf( ipPnunit, "\n" );
1775 }
1776 }
1777
1778 /* this is the charge transfer to produce output for AGN3,
1779 * invoked with the save charge transfer AGN command */
1780 else if( strcmp( chSave,"CHAG") == 0 )
1781 {
1782 /* this is boundary between two tables */
1783 double BreakEnergy = 100./13.0;
1784 realnum teinit = 5000.;
1785 realnum tefinal = 20000.;
1786 realnum te1;
1787 long int nelem, ion;
1788
1789 te1 = teinit;
1790 fprintf(ipPnunit,"H ioniz\n X+i\\Te");
1791 while( te1 <= tefinal )
1792 {
1793 fprintf(ipPnunit,"\t%.0f K",te1);
1794 te1 *= 2.;
1795 }
1796 fprintf(ipPnunit,"\n");
1797
1798 /* make sure rates already evaluated at least one time */
1799 ChargTranEval();
1800
1801 /* loop over all elements, H charge transfer ionization */
1802 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
1803 {
1804 /* this list of elements included in the AGN tables is defined in zeroabun.c */
1805 if( abund.lgAGN[nelem] )
1806 {
1807 for( ion=0; ion<=nelem; ++ion )
1808 {
1809 /* skip high ionization CT */
1810 if( Heavy.Valence_IP_Ryd[nelem][ion] > BreakEnergy )
1811 break;
1812 /* most of these are actually zero */
1813 if( atmdat.CharExcIonOf[ipHYDROGEN][nelem][ion] == 0 )
1814 continue;
1815
1816 /* print chemical symbol */
1817 fprintf(ipPnunit,"%s",
1818 elementnames.chElementSym[nelem]);
1819 /* now ionization stage */
1820 if( ion==0 )
1821 {
1822 fprintf(ipPnunit,"0 ");
1823 }
1824 else if( ion==1 )
1825 {
1826 fprintf(ipPnunit,"+ ");
1827 }
1828 else
1829 {
1830 fprintf(ipPnunit,"+%li",ion);
1831 }
1832
1833 /* fully define the new temperature */
1834 TempChange(teinit , false);
1835
1836 while( phycon.te <= tefinal )
1837 {
1838 dense.IonLow[nelem] = 0;
1839 dense.IonHigh[nelem] = nelem+1;
1840 ChargTranEval();
1841
1842 fprintf(ipPnunit,"\t%.2e",atmdat.CharExcIonOf[ipHYDROGEN][nelem][ion]);
1843 TempChange(phycon.te *2.f , false);
1844 }
1845 fprintf(ipPnunit,"\n");
1846 }
1847 fprintf(ipPnunit,"\n");
1848 }
1849 }
1850
1851 te1 = teinit;
1852 fprintf(ipPnunit,"H recom\n X+i\\Te");
1853 while( te1 <= tefinal )
1854 {
1855 fprintf(ipPnunit,"\t%.0f K",te1);
1856 te1 *= 2.;
1857 }
1858 fprintf(ipPnunit,"\n");
1859
1860 /* loop over all elements, H charge transfer recombination */
1861 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
1862 {
1863 /* this list of elements included in the AGN tables is defined in zeroabun.c */
1864 if( abund.lgAGN[nelem] )
1865 {
1866 for( ion=0; ion<=nelem; ++ion )
1867 {
1868 /* skip high ionization CT */
1869 if( Heavy.Valence_IP_Ryd[nelem][ion] > BreakEnergy )
1870 break;
1871 /* most of these are actually zero */
1872 if( atmdat.CharExcRecTo[ipHYDROGEN][nelem][ion] == 0 )
1873 continue;
1874
1875 /* print chemical symbol */
1876 fprintf(ipPnunit,"%s",
1877 elementnames.chElementSym[nelem]);
1878 /* now ionization stage */
1879 if( ion==0 )
1880 {
1881 fprintf(ipPnunit,"0 ");
1882 }
1883 else if( ion==1 )
1884 {
1885 fprintf(ipPnunit,"+ ");
1886 }
1887 else
1888 {
1889 fprintf(ipPnunit,"+%li",ion);
1890 }
1891
1892 /* fully define the new temperature */
1893 TempChange(teinit , false);
1894 while( phycon.te <= tefinal )
1895 {
1896 dense.IonLow[nelem] = 0;
1897 dense.IonHigh[nelem] = nelem+1;
1898 ChargTranEval();
1899
1900 fprintf(ipPnunit,"\t%.2e",atmdat.CharExcRecTo[ipHYDROGEN][nelem][ion]);
1901 TempChange(phycon.te *2.f , false);
1902 }
1903 fprintf(ipPnunit,"\n");
1904 }
1905 fprintf(ipPnunit,"\n");
1906 }
1907 }
1908
1909# if 0
1910 te1 = teinit;
1911 fprintf(ipPnunit,"He recom\n Elem\\Te");
1912 while( te1 <= tefinal )
1913 {
1914 fprintf(ipPnunit,"\t%.0f",te1);
1915 te1 *= 2.;
1916 }
1917 fprintf(ipPnunit,"\n");
1918
1919 /* loop over all elements, H charge transfer recombination */
1920 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
1921 {
1922 /* this list of elements included in the AGN tables is defined in zeroabun.c */
1923 if( abund.lgAGN[nelem] )
1924 {
1925 for( ion=0; ion<=nelem; ++ion )
1926 {
1927 /* most of these are actually zero */
1928 if( atmdat.CharExcRecTo[ipHELIUM][nelem][ion] == 0 )
1929 continue;
1930 fprintf(ipPnunit,"%.2s%.2s",
1931 elementnames.chElementSym[nelem],
1932 elementnames.chIonStage[ion]);
1933
1934 /* fully define the new temperature */
1935 TempChange(teinit , false);
1936 while( phycon.te <= tefinal )
1937 {
1938 dense.IonLow[nelem] = 0;
1939 dense.IonHigh[nelem] = nelem+1;
1940 ChargTranEval();
1941
1942 fprintf(ipPnunit,"\t%.2e",atmdat.CharExcRecTo[ipHELIUM][nelem][ion]);
1943 TempChange(phycon.te *2.fprintf , false);
1944 }
1945 fprintf(ipPnunit,"\n");
1946 }
1947 fprintf(ipPnunit,"\n");
1948 }
1949 }
1950
1951
1952 te1 = teinit;
1953 fprintf(ipPnunit,"He ioniz\n Elem\\Te");
1954 while( te1 <= tefinal )
1955 {
1956 fprintf(ipPnunit,"\t%.0f",te1);
1957 te1 *= 2.;
1958 }
1959 fprintf(ipPnunit,"\n");
1960
1961 /* loop over all elements, H charge transfer recombination */
1962 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
1963 {
1964 /* this list of elements included in the AGN tables is defined in zeroabun.c */
1965 if( abund.lgAGN[nelem] )
1966 {
1967 for( ion=0; ion<=nelem; ++ion )
1968 {
1969 /* most of these are actually zero */
1970 if( atmdat.CharExcIonOf[ipHELIUM][nelem][ion] == 0 )
1971 continue;
1972 fprintf(ipPnunit,"%.2s%.2s",
1973 elementnames.chElementSym[nelem],
1974 elementnames.chIonStage[ion]);
1975
1976 /* fully define the new temperature */
1977 TempChange(teinit , false);
1978 while( phycon.te <= tefinal )
1979 {
1980 dense.IonLow[nelem] = 0;
1981 dense.IonHigh[nelem] = nelem+1;
1982 ChargTranEval();
1983
1984 fprintf(ipPnunit,"\t%.2e",atmdat.CharExcIonOf[ipHELIUM][nelem][ion]);
1985 TempChange(phycon.te*2.f , true);
1986 }
1987 fprintf(ipPnunit,"\n");
1988 }
1989 fprintf(ipPnunit,"\n");
1990 }
1991 }
1992# endif
1993 }
1994 else
1995 {
1996 fprintf( ioQQQ, " save charge keyword insane\n" );
1998 }
1999
2000 TempChange(TempSave , false);
2001 return;
2002}
t_abund abund
Definition abund.cpp:5
t_atmdat atmdat
Definition atmdat.cpp:6
double ChargTranSumHeat(void)
void ChargTranPun(FILE *ipPnunit, char *chSave)
STATIC double HCTIon(long int ion, long int nelem)
#define FRAC
void ChargTranEval(void)
static double CTIonData[LIMELM][4][8]
static bool lgCTDataDefined
static double CTRecombData[LIMELM][4][7]
STATIC double HCTRecom(long int ion, long int nelem)
STATIC void MakeHCTData(void)
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
sys_float sexp(sys_float x)
Definition service.cpp:914
const int ipSILICON
Definition cddefines.h:318
#define STATIC
Definition cddefines.h:97
const int ipNITROGEN
Definition cddefines.h:311
const int ipIRON
Definition cddefines.h:330
const int ipCARBON
Definition cddefines.h:310
const int ipALUMINIUM
Definition cddefines.h:317
#define MIN2
Definition cddefines.h:761
const int ipOXYGEN
Definition cddefines.h:312
const int LIMELM
Definition cddefines.h:258
const int ipPHOSPHORUS
Definition cddefines.h:319
const int ipLITHIUM
Definition cddefines.h:307
const int ipMAGNESIUM
Definition cddefines.h:316
#define EXIT_FAILURE
Definition cddefines.h:140
const int ipSODIUM
Definition cddefines.h:315
const int ipTITANIUM
Definition cddefines.h:326
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
const int ipNICKEL
Definition cddefines.h:332
const int ipSULPHUR
Definition cddefines.h:320
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
double e2(double t)
Definition service.cpp:299
const int ipPOTASSIUM
Definition cddefines.h:323
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
const int ipNEON
Definition cddefines.h:314
const int ipCHLORINE
Definition cddefines.h:321
const int ipMANGANESE
Definition cddefines.h:329
const int ipARGON
Definition cddefines.h:322
void fixit(void)
Definition service.cpp:991
#define HMRATE(a, b, c)
Definition cddefines.h:1046
t_conv conv
Definition conv.cpp:5
t_dense dense
Definition dense.cpp:24
t_elementnames elementnames
t_Heavy Heavy
Definition heavy.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH1s
Definition iso.h:27
const int ipH_LIKE
Definition iso.h:62
t_mole_global mole_global
Definition mole.cpp:6
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double EN1EV
Definition physconst.h:192
void TempChange(double TempNew, bool lgForceUpdate)
t_thermal thermal
Definition thermal.cpp:5
t_trace trace
Definition trace.cpp:5