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

Generated for cloudy by doxygen 1.8.1.1