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