cloudy  trunk
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+ */
25 STATIC double HCTIon(long int ion,
26  long int nelem);
27 
28 /*HCTRecom H charge transfer recombination, H + A+ => H+ + A */
29 STATIC double HCTRecom(long int ion,
30  long int nelem);
31 
32 /* the translated block data */
33 STATIC 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 */
37 static double CTIonData[LIMELM][4][8];
38 static double CTRecombData[LIMELM][4][7];
39 /* this will be flag for whether or not charge transfer data
40  * have been initialized */
41 static bool lgCTDataDefined = false;
42 
43 /*ChargTranEval update charge trans rate coefficients if temperature has changed */
44 void 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 */
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 */
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 */
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 */
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 */
211 
212  /* >>refer P CT Pequignot, D., & Aldrovandi, S. M. V. 1986, A&A, 161, 169-176 */
213  /* P0 + H+ => P+ + H0 */
215 
216  /* >>refer Cl CT Pequignot, D., & Aldrovandi, S. M. V. 1986, A&A, 161, 169-176 */
217  /* Cl0 + H+ => Cl+ + H0 */
219 
220  /* >>refer Ti CT Pequignot, D., & Aldrovandi, S. M. V. 1986, A&A, 161, 169-176 */
221  /* Ti0 + H+ => Cl+ + H0 */
223 
224  /* >>refer Mn CT Pequignot, D., & Aldrovandi, S. M. V. 1986, A&A, 161, 169-176 */
225  /* Mn0 + H+ => Mn+ + H0 */
227 
228  /* >>refer Ni CT Pequignot, D., & Aldrovandi, S. M. V. 1986, A&A, 161, 169-176 */
229  /* Ni0 + H+ => Ni+ + H0 */
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  {
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 */
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*/
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 */
415 
416  /* ionization */
417  /* N+1 + He+ => N+2 + He */
418  /* >>refer N1 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
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 */
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 */
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 */
449  /* Ne+4 + He => Ne+3 + He+ */
450  /* >>refer Ne4 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
502 
503  /* ionization */
504  /* S+1 + He+ => S+2 + He */
505  /* >>refer S1 CT Butler, S. E., & Dalgarno, A. 1980b, ApJ, 241, 838 */
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 */
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 */
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 
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 
541 
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  *================================================================================*/
569 double ChargTranSumHeat(void)
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]*
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()*
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 *
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*
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()*
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 *
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 */
676  SumCTHeat_v/thermal.htot );
677 
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  *================================================================================*/
736 STATIC double HCTIon(
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 */
898 STATIC void MakeHCTData(void)
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 */
1729 void 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 }
thermal.h
t_atmdat::CharExcIonOf
double CharExcIonOf[NCX][LIMELM][LIMELM+1]
Definition: atmdat.h:152
ipOXYGEN
const int ipOXYGEN
Definition: cddefines.h:312
lgCTDataDefined
static bool lgCTDataDefined
Definition: atmdat_char_tran.cpp:41
dense
t_dense dense
Definition: dense.cpp:24
ipMANGANESE
const int ipMANGANESE
Definition: cddefines.h:329
elementnames.h
t_phycon::tesqrd
double tesqrd
Definition: phycon.h:26
t_atmdat::HCharCoolMax
double HCharCoolMax
Definition: atmdat.h:155
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
TempChange
void TempChange(double TempNew, bool lgForceUpdate)
Definition: temp_change.cpp:51
realnum
float realnum
Definition: cddefines.h:103
t_elementnames::chIonStage
char chIonStage[LIMELM+1][CHARS_ION_STAGE]
Definition: elementnames.h:29
conv.h
STATIC
#define STATIC
Definition: cddefines.h:97
ipLITHIUM
const int ipLITHIUM
Definition: cddefines.h:307
ipCARBON
const int ipCARBON
Definition: cddefines.h:310
mole.h
abund.h
ipMAGNESIUM
const int ipMAGNESIUM
Definition: cddefines.h:316
ipIRON
const int ipIRON
Definition: cddefines.h:330
t_phycon::te03
double te03
Definition: phycon.h:59
phycon
t_phycon phycon
Definition: phycon.cpp:6
t_atmdat::HCharHeatOn
double HCharHeatOn
Definition: atmdat.h:156
trace.h
ipNEON
const int ipNEON
Definition: cddefines.h:314
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:952
ChargTranEval
void ChargTranEval(void)
Definition: atmdat_char_tran.cpp:44
HCTRecom
STATIC double HCTRecom(long int ion, long int nelem)
Definition: atmdat_char_tran.cpp:806
Heavy
t_Heavy Heavy
Definition: heavy.cpp:5
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
iso.h
ChargTranPun
void ChargTranPun(FILE *ipPnunit, char *chSave)
Definition: atmdat_char_tran.cpp:1729
ipNITROGEN
const int ipNITROGEN
Definition: cddefines.h:311
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
HMRATE
#define HMRATE(a, b, c)
Definition: cddefines.h:1046
atmdat.h
MIN2
#define MIN2
Definition: cddefines.h:761
t_phycon::te_eV
double te_eV
Definition: phycon.h:14
sexp
sys_float sexp(sys_float x)
Definition: service.cpp:914
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
dense.h
ipSILICON
const int ipSILICON
Definition: cddefines.h:318
t_phycon::te01
double te01
Definition: phycon.h:61
trace
t_trace trace
Definition: trace.cpp:5
t_Heavy::Valence_IP_Ryd
double Valence_IP_Ryd[LIMELM][LIMELM]
Definition: heavy.h:24
cddefines.h
HCTIon
STATIC double HCTIon(long int ion, long int nelem)
Definition: atmdat_char_tran.cpp:736
ipSULPHUR
const int ipSULPHUR
Definition: cddefines.h:320
EN1EV
const UNUSED double EN1EV
Definition: physconst.h:192
e2
double e2(double t)
Definition: service.cpp:299
thermal
t_thermal thermal
Definition: thermal.cpp:5
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
heavy.h
t_phycon::alnte
double alnte
Definition: phycon.h:85
ipALUMINIUM
const int ipALUMINIUM
Definition: cddefines.h:317
ipCHLORINE
const int ipCHLORINE
Definition: cddefines.h:321
MAX2
#define MAX2
Definition: cddefines.h:782
t_phycon::te30
double te30
Definition: phycon.h:53
LIMELM
const int LIMELM
Definition: cddefines.h:258
ipPHOSPHORUS
const int ipPHOSPHORUS
Definition: cddefines.h:319
t_mole_global::lgLeidenHack
bool lgLeidenHack
Definition: mole.h:286
t_dense::IonLow
long int IonLow[LIMELM+1]
Definition: dense.h:119
ipPOTASSIUM
const int ipPOTASSIUM
Definition: cddefines.h:323
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
t_iso_sp::st
qList st
Definition: iso.h:453
t_atmdat::HCharExcRecTo_N0_2D
double HCharExcRecTo_N0_2D
Definition: atmdat.h:159
t_dense::IonHigh
long int IonHigh[LIMELM+1]
Definition: dense.h:120
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
MakeHCTData
STATIC void MakeHCTData(void)
Definition: atmdat_char_tran.cpp:898
abund
t_abund abund
Definition: abund.cpp:5
t_thermal::htot
double htot
Definition: thermal.h:149
ChargTranSumHeat
double ChargTranSumHeat(void)
Definition: atmdat_char_tran.cpp:569
CTRecombData
static double CTRecombData[LIMELM][4][7]
Definition: atmdat_char_tran.cpp:38
t_conv::nTotalIoniz
long int nTotalIoniz
Definition: conv.h:166
t_phycon::te10
double te10
Definition: phycon.h:55
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:812
physconst.h
t_atmdat::HCTAlex
double HCTAlex
Definition: atmdat.h:173
t_abund::lgAGN
bool lgAGN[LIMELM]
Definition: abund.h:44
conv
t_conv conv
Definition: conv.cpp:5
FRAC
#define FRAC
fixit
void fixit(void)
Definition: service.cpp:991
t_phycon::te20
double te20
Definition: phycon.h:54
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
taulines.h
ipSODIUM
const int ipSODIUM
Definition: cddefines.h:315
phycon.h
atmdat
t_atmdat atmdat
Definition: atmdat.cpp:6
ipH1s
const int ipH1s
Definition: iso.h:27
t_atmdat::lgCTOn
bool lgCTOn
Definition: atmdat.h:177
t_phycon::sqrte
double sqrte
Definition: phycon.h:48
t_phycon::te05
double te05
Definition: phycon.h:57
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
t_atmdat::CharExcRecTo
double CharExcRecTo[NCX][LIMELM][LIMELM+1]
Definition: atmdat.h:153
ipARGON
const int ipARGON
Definition: cddefines.h:322
CTIonData
static double CTIonData[LIMELM][4][8]
Definition: atmdat_char_tran.cpp:37
t_phycon::te
double te
Definition: phycon.h:11
t_atmdat::HCharHeatMax
double HCharHeatMax
Definition: atmdat.h:154
mole_global
t_mole_global mole_global
Definition: mole.cpp:6
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
ipNICKEL
const int ipNICKEL
Definition: cddefines.h:332
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
ipTITANIUM
const int ipTITANIUM
Definition: cddefines.h:326