root/branches/newmole/source/atmdat_char_tran.cpp

Revision 2542, 57.6 kB (checked in by rjrw, 3 weeks ago)

Merged from trunk r2448:2541

  • Property svn:eol-style set to native
Line 
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+ */
24STATIC double HCTIon(long int ion,
25  long int nelem);
26
27/*HCTRecom H charge transfer recombination, H + A+ => H+ + A */
28STATIC double HCTRecom(long int ion,
29  long int nelem);
30
31/* the translated block data */
32STATIC 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  */
36static double CTIonData[LIMELM][4][8];
37static double CTRecombData[LIMELM][4][7];
38/* this will be flag for whether or not charge transfer data
39 * have been initialized */
40static bool lgCTDataDefined = false;
41
42/*ChargTranEval update charge trans rate coefficients if temperature has changed */
43void ChargTranEval(
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 || !fp_equal(phycon.te,TeUsed) )
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(!mole.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(!mole.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 */
282                        /** \todo       1       update ct to Kimura et al. (1996) */
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));
291                        /** \todo       1       above rate not intended for very low temperatures - find ref for low-T rate,
292                         * probably is 1e-9 like above */
293                }
294                else
295                {
296                        /** \todo       0       these should be values at 1e5 K */
297                        atmdat.HCharExcIonOf[ipMAGNESIUM][0] = 0.;
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 */
317                        /** \todo       2       not currently used - include as deexcitation process */
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 */
377                atmdat.HeCharExcRecTo[ipCARBON][2] = 4.6e-19*phycon.tesqrd;
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*/
393                atmdat.HeCharExcIonOf[ipCARBON][1] =
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 */
418                atmdat.HeCharExcIonOf[ipNITROGEN][1] =
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 */
425                atmdat.HeCharExcRecTo[ipOXYGEN][1] = 3.2e-14*phycon.te/phycon.te05;
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 */
436                atmdat.HeCharExcIonOf[ipOXYGEN][0] =
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 */
447                atmdat.HeCharExcRecTo[ipNEON][2] = 1e-16*phycon.sqrte;
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 */
459                atmdat.HeCharExcRecTo[ipMAGNESIUM][3] = 1.4e-10*phycon.te30;
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 */
468                atmdat.HeCharExcRecTo[ipSILICON][2] += phycon.sqrte*phycon.te10*phycon.te10*
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 */
473                atmdat.HeCharExcRecTo[ipSILICON][3] = 2.54e-11*phycon.sqrte/phycon.te03/
474                  phycon.te01/phycon.te01;
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 */
483                atmdat.HeCharExcIonOf[ipSILICON][1] =
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 */
488                atmdat.HeCharExcIonOf[ipSILICON][2] =
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 */
495                atmdat.HeCharExcRecTo[ipSULPHUR][2] = phycon.sqrte*1.1e-11;
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 */
505                atmdat.HeCharExcIonOf[ipSULPHUR][1] =
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 */
510                atmdat.HeCharExcIonOf[ipSULPHUR][2] =
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(!mole.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
541                        atmdat.HeCharExcRecTo[ipCARBON][0] = 0;
542                        atmdat.HeCharExcIonOf[ipOXYGEN][0] = 0;
543                        atmdat.HeCharExcRecTo[ipOXYGEN][0] = 0;
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 *================================================================================*/
570double 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 */
583        ASSERT( lgCTDataDefined );
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]*
608                                StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop*dense.xIonDense[ipHYDROGEN][1]*
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 */
629                atmdat.HCharHeatMax = MAX2(atmdat.HCharHeatMax,
630                        SumCTHeat_v/thermal.htot );
631
632                atmdat.HCharCoolMax = MAX2(atmdat.HCharCoolMax,
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]*
660                                                StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop*dense.xIonDense[ipHYDROGEN][1]*
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]*
669                                                        StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop*dense.xIonDense[ipHYDROGEN][1]*
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 *================================================================================*/
692STATIC 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                // make sure that the Boltzmann factor always uses the real temperature, even
751                // outside the range of validity, so that the CT rate properly goes to zero
752                // when the the gas temperature goes to zero.
753                HCTIon_v = CTIonData[nelem][ion][0]*1e-9*(pow(tused,CTIonData[nelem][ion][1]))*
754                  (1. + CTIonData[nelem][ion][2]*exp(CTIonData[nelem][ion][3]*tused))*
755                  exp(-CTIonData[nelem][ion][6]*1.e4/phycon.te);
756        }
757        return( HCTIon_v );
758}
759
760/*================================================================================*
761 *================================================================================*/
762STATIC double HCTRecom(
763        /* ion is stage of ionization on C scale, 0 for rec to atom */
764        long int ion,   
765        /* nelem is atomic number of element on C scale, 1 = he up to LIMELM */
766        /* HCTRecom(1,5) would be C+2 + H => C+ + H+ */
767        long int nelem)
768{
769        double HCTRecom_v,
770          tused;
771
772        DEBUG_ENTRY( "HCTRecom()" );
773        /*
774         * H charge transfer recombination using Jim Kingdon's block ctdata.for
775         */
776
777        /* set up the rate coefficients if this is first call */
778        if( !lgCTDataDefined )
779        {
780                if( trace.lgTrace )
781                {
782                        fprintf( ioQQQ,"       HCTIon doing 1-time init of charge transfer data\n");
783                }
784                lgCTDataDefined = true;
785                MakeHCTData();
786        }
787
788        /* this is check that data have been set up properly, will
789         * fail if arrays are not initialized properly */
790        ASSERT( CTRecombData[1][0][0] > 0. );
791
792        /* use Dalgarno estimate for highly ionized species, number reset with
793         * set charge transfer command */
794        if( ion > 3 )
795        {
796                /* >>chng 96 nov 25, added this option, default is 1.92e-9
797                 * Dalgarno's charge transfer */
798                HCTRecom_v = atmdat.HCTAlex*((double)ion+1.);
799                return( HCTRecom_v );
800        }
801
802        /* check that ion stage within bound for this atom */
803        ASSERT( ion >= 0 && ion <= nelem );
804
805        /* now check the element is valid, this is routine HCTIon */
806        ASSERT( nelem > 0 && nelem < LIMELM );
807
808        tused = MAX2((double)phycon.te,CTRecombData[nelem][ion][4]);
809        tused = MIN2(tused,CTRecombData[nelem][ion][5]);
810        tused *= 1e-4;
811
812        if( tused == 0. )
813        {
814                HCTRecom_v = 0.;
815                return( HCTRecom_v );
816        }
817
818        /* the interpolation equation */
819        HCTRecom_v = CTRecombData[nelem][ion][0]*1e-9*(pow(tused,CTRecombData[nelem][ion][1]))*
820          (1. + CTRecombData[nelem][ion][2]*sexp(-CTRecombData[nelem][ion][3]*tused));
821
822        /* in sexp negative sign not typo - there are negative signs already
823         * in coefficient, and sexp has implicit negative sign */
824        return( HCTRecom_v );
825}
826
827/*================================================================================*
828 *================================================================================*/
829/*block data with Jim Kingdon's charge transfer data */
830/* >>refer      H       CT      Kingdon, J, B., & Ferland, G.J. 1996, ApJS, 106, 205 */
831/*
832 * first dimension is atomic number of atom, 0 for H
833 * second dimension is ionization stage,
834 * 1=+0 for parent, etc
835 * third dimension is atomic number of atom
836 * second dimension is ionization stage,
837 * 1=+1 for parent, etc
838 */
839
840/* digital form of the fits to the charge transfer
841 * ionization rate coefficients
842 *
843 * Note: First parameter is in units of 1e-9!
844 * Note: Seventh parameter is in units of 1e4 K */
845
846/* digital form of the fits to the charge transfer
847 * recombination rate coefficients (total)
848 *
849 * Note: First parameter is in units of 1e-9!
850 * recombination
851 */
852
853/* holds data for charge transfer fits */
854STATIC void MakeHC