root/branches/newmole/source/atmdat_adfa.cpp

Revision 2542, 21.5 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/*phfit derive photoionization cross sections for first 30 elements */
4#include "cddefines.h"
5#include "physconst.h"
6#include "atmdat.h"
7#include "iso.h"
8
9/** constructor: read in all the ADfA data files */
10t_ADfA::t_ADfA()
11{
12        DEBUG_ENTRY( "t_ADfA()" );
13
14        /* this option, use the new atmdat_rad_rec recombination rates */
15        version = PHFIT_UNDEF;
16
17        double help[9];
18        const long VERSION_MAGIC = 20061204L;
19
20        const static char chFile[] = "phfit.dat";
21
22        FILE *io = open_data( chFile, "r" );
23
24        bool lgErr = false;
25        long i=-1, j=-1, k=-1, n;
26
27        lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
28        if( lgErr || i != VERSION_MAGIC )
29        {
30                fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile, i );
31                fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
32                cdEXIT(EXIT_FAILURE);
33        }
34
35        for( n=0; n < 7; n++ )
36                lgErr = lgErr || ( fscanf( io, "%ld", &L[n] ) != 1 );
37        for( n=0; n < 30; n++ )
38                lgErr = lgErr || ( fscanf( io, "%ld", &NINN[n] ) != 1 );
39        for( n=0; n < 30; n++ )
40                lgErr = lgErr || ( fscanf( io, "%ld", &NTOT[n] ) != 1 );
41        while( true )
42        {
43                lgErr = lgErr || ( fscanf( io, "%ld %ld %ld", &i, &j, &k ) != 3 );
44                if( i == -1 && j == -1 && k == -1 )
45                        break;
46                lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le", &help[0], &help[1],
47                                           &help[2], &help[3], &help[4], &help[5] ) != 6 );
48                for( int l=0; l < 6; ++l )
49                        PH1[i][j][k][l] = (realnum)help[l];
50        }
51        while( true )
52        {
53                lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
54                if( i == -1 && j == -1 )
55                        break;
56                lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le", &help[0], &help[1],
57                                           &help[2], &help[3], &help[4], &help[5], &help[6] ) != 7 );
58                for( int l=0; l < 7; ++l )
59                        PH2[i][j][l]  = (realnum)help[l];
60        }
61        fclose( io );
62
63        ASSERT( !lgErr );
64
65        const static char chFile2[] = "hpfit.dat";
66
67        io = open_data( chFile2, "r" );
68
69        lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
70        if( lgErr || i != VERSION_MAGIC )
71        {
72                fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile2, i );
73                fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
74                cdEXIT(EXIT_FAILURE);
75        }
76
77        for( i=0; i < NHYDRO_MAX_LEVEL; i++ )
78        {
79                lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le", &help[0], &help[1],
80                                           &help[2], &help[3], &help[4] ) != 5 );
81                for( int l=0; l < 5; ++l )
82                        PHH[i][l] = (realnum)help[l];
83        }
84
85        fclose( io );
86
87        ASSERT( !lgErr );
88
89        const static char chFile3[] = "rec_lines.dat";
90
91        io = open_data( chFile3, "r" );
92
93        lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
94        if( lgErr || i != VERSION_MAGIC )
95        {
96                fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile3, i );
97                fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
98                cdEXIT(EXIT_FAILURE);
99        }
100
101        for( i=0; i < 110; i++ )
102        {
103                lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2],
104                                           &help[3], &help[4], &help[5], &help[6], &help[7] ) != 8 );
105                for( int l=0; l < 8; ++l )
106                        P[l][i] = (realnum)help[l];
107        }
108
109
110        for( i=0; i < 405; i++ )
111        {
112                lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2],
113                                           &help[3], &help[4], &help[5], &help[6], &help[7], &help[8] ) != 9 );
114                for( int l=0; l < 9; ++l )
115                        ST[l][i] = (realnum)help[l];
116        }
117
118        fclose( io );
119
120        ASSERT( !lgErr );
121
122        const static char chFile4[] = "rad_rec.dat";
123
124        io = open_data( chFile4, "r" );
125
126        lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
127        if( lgErr || i != VERSION_MAGIC )
128        {
129                fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile4, i );
130                fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
131                cdEXIT(EXIT_FAILURE);
132        }
133
134        while( true )
135        {
136                lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
137                if( i == -1 && j == -1 )
138                        break;
139                lgErr = lgErr || ( fscanf( io, "%le %le", &help[0], &help[1] ) != 2 );
140                for( int l=0; l < 2; ++l )
141                        rrec[i][j][l] = (realnum)help[l];
142        }
143        while( true )
144        {
145                lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
146                if( i == -1 && j == -1 )
147                        break;
148                lgErr = lgErr || ( fscanf( io, "%le %le %le %le", &help[0], &help[1],
149                                           &help[2], &help[3] ) != 4 );
150                for( int l=0; l < 4; ++l )
151                        rnew[i][j][l] = (realnum)help[l];
152        }
153        while( true )
154        {
155                lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
156                if( i == -1 )
157                        break;
158                lgErr = lgErr || ( fscanf( io, "%le %le %le", &help[0], &help[1], &help[2] ) != 3 );
159                for( int l=0; l < 3; ++l )
160                        fe[i][l] = (realnum)help[l];
161        }
162
163        fclose( io );
164
165        ASSERT( !lgErr );
166
167        const static char chFile5[] = "h_rad_rec.dat";
168
169        io = open_data( chFile5, "r" );
170
171        lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
172        if( lgErr || i != VERSION_MAGIC )
173        {
174                fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile5, i );
175                fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
176                cdEXIT(EXIT_FAILURE);
177        }
178
179        for( i=0; i < NHYDRO_MAX_LEVEL; i++ )
180        {
181                lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2],
182                                           &help[3], &help[4], &help[5], &help[6], &help[7], &help[8] ) != 9 );
183                for( int l=0; l < 9; ++l )
184                        HRF[i][l] = (realnum)help[l];
185        }
186
187        fclose( io );
188
189        ASSERT( !lgErr );
190
191        const static char chFile6[] = "h_phot_cs.dat";
192
193        io = open_data( chFile6, "r" );
194
195        lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
196        if( lgErr || i != VERSION_MAGIC )
197        {
198                fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile6, i );
199                fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
200                cdEXIT(EXIT_FAILURE);
201        }
202
203        for( i=0; i < NHYDRO_MAX_LEVEL; i++ )
204        {
205                lgErr = lgErr || ( fscanf( io, "%le", &help[0] ) != 1 );
206                STH[i] = (realnum)help[0];
207        }
208
209        fclose( io );
210
211        ASSERT( !lgErr );
212
213        const static char chFile7[] = "coll_ion.dat";
214
215        io = open_data( chFile7, "r" );
216
217        lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
218        if( lgErr || i != VERSION_MAGIC )
219        {
220                fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile7, i );
221                fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
222                cdEXIT(EXIT_FAILURE);
223        }
224
225        while( true )
226        {
227                lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
228                if( i == -1 && j == -1 )
229                        break;
230                lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le", &CF[i][j][0], &CF[i][j][1],
231                                           &CF[i][j][2], &CF[i][j][3], &CF[i][j][4] ) != 5 );
232        }
233
234        fclose( io );
235
236        ASSERT( !lgErr );
237
238        /*refer HI      cs      Anderson, H., Ballance, C.P., Badnell, N.R.,
239         *refercon      & Summers, H.P  2000, J Phys B, 33, 1255 */
240        const static char chFile8[] = "h_coll_str.dat";
241
242        io = open_data( chFile8, "r" );
243
244        lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
245        if( lgErr || i != VERSION_MAGIC )
246        {
247                fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile8, i );
248                fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
249                cdEXIT(EXIT_FAILURE);
250        }
251
252        while( true )
253        {
254                lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
255                if( i == -1 && j == -1 )
256                        break;
257                lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le", &HCS[i-2][j-1][0], &HCS[i-2][j-1][1],
258                                           &HCS[i-2][j-1][2], &HCS[i-2][j-1][3], &HCS[i-2][j-1][4], &HCS[i-2][j-1][5],
259                                           &HCS[i-2][j-1][6], &HCS[i-2][j-1][7] ) != 8 );
260        }
261
262        fclose( io );
263
264        ASSERT( !lgErr );
265}
266
267double t_ADfA::phfit(long int nz,
268                     long int ne,
269                     long int is,
270                     double e)
271{
272        long int nint,
273          nout;
274        double a,
275          b,
276          crs,
277          einn,
278          p1,
279          q,
280          x,
281          y,
282          z;
283
284        DEBUG_ENTRY( "phfit()" );
285
286        /*** Version 3. October 8, 1996.
287         *** Written by D. A. Verner, verner@pa.uky.edu
288         *** Inner-shell ionization energies of some low-ionized species are slightly
289         *** improved to fit smoothly the experimental inner-shell ionization energies
290         *** of neutral atoms.
291         ******************************************************************************
292         *** This subroutine calculates partial photoionization cross sections
293         *** for all ionization stages of all atoms from H to Zn (Z=30) by use of
294         *** the following fit parameters:
295         *** Outer shells of the Opacity Project (OP) elements:
296         *** >>refer    all     photocs Verner, Ferland, Korista, Yakovlev, 1996, ApJ, in press.
297         *** Inner shells of all elements, and outer shells of the non-OP elements:
298         ***  Verner and Yakovlev, 1995, A&AS, 109, 125
299         *** Input parameters:  nz - atomic number from 1 to 30 (integer)
300         ***          ne - number of electrons from 1 to iz (integer)
301         ***          is - shell number (integer)
302         ***          e - photon energy, eV
303         ***          version - enum, PHFIT96 (default): calculates
304         ***                 new cross sections, PHFIT95: calculates
305         ***                 only old Hartree-Slater cross sections
306         *** Output parameter:  photoionization cross section, Mb
307         *** Shell numbers:
308         *** 1 - 1s, 2 - 2s, 3 - 2p, 4 - 3s, 5 - 3p, 6 - 3d, 7 - 4s.
309         *** If a species in the ground state has no electrons on the given shell,
310         *** the subroutine returns 0.
311         ****************************************************************************** */
312
313        crs = 0.0;
314        if( nz < 1 || nz > 30 )
315        {
316                return(crs);
317        }
318
319        if( ne < 1 || ne > nz )
320        {
321                return(crs);
322        }
323
324        nout = NTOT[ne-1];
325        if( nz == ne && nz > 18 )
326                nout = 7;
327        if( nz == (ne + 1) && ((((nz == 20 || nz == 21) || nz == 22) ||
328          nz == 25) || nz == 26) )
329                nout = 7;
330        if( is > nout )
331        {
332                return(crs);
333        }
334
335        if( (is == 6 && (nz == 20 || nz == 19)) && ne >= 19 )
336        {
337                return(crs);
338        }
339
340        ASSERT( is >= 1 && is <= 7 );
341
342        if( e < PH1[is-1][ne-1][nz-1][0] )
343        {
344                return(crs);
345        }
346
347        nint = NINN[ne-1];
348        if( ((nz == 15 || nz == 17) || nz == 19) || (nz > 20 && nz != 26) )
349        {
350                einn = 0.0;
351        }
352        else
353        {
354                if( ne < 3 )
355                {
356                        einn = 1.0e30;
357                }
358                else
359                {
360                        einn = PH1[nint-1][ne-1][nz-1][0];
361                }
362        }
363
364        if( (is <= nint || e >= einn) || version == PHFIT95 )
365        {
366                p1 = -PH1[is-1][ne-1][nz-1][4];
367                y = e/PH1[is-1][ne-1][nz-1][1];
368                q = -0.5*p1 - L[is-1] - 5.5;
369                a = PH1[is-1][ne-1][nz-1][2]*(POW2(y - 1.0) +
370                        POW2(PH1[is-1][ne-1][nz-1][5]));
371                b = sqrt(y/PH1[is-1][ne-1][nz-1][3]) + 1.0;
372                crs = a*pow(y,q)*pow(b,p1);
373        }
374        else
375        {
376                if( (is < nout && is > nint) && e < einn )
377                {
378                        return(crs);
379                }
380                p1 = -PH2[ne-1][nz-1][3];
381                q = -0.5*p1 - 5.5;
382                x = e/PH2[ne-1][nz-1][0] - PH2[ne-1][nz-1][5];
383                z = sqrt(x*x+POW2(PH2[ne-1][nz-1][6]));
384                a = PH2[ne-1][nz-1][1]*(POW2(x - 1.0) +
385                        POW2(PH2[ne-1][nz-1][4]));
386                b = 1.0 + sqrt(z/PH2[ne-1][nz-1][2]);
387                crs = a*pow(z,q)*pow(b,p1);
388        }
389        return(crs);
390}
391
392/* same as hpfit, but energy is relative to threshold */
393double t_ADfA::hpfit_rel(long int iz,
394                         long int n,
395                         double e)
396{
397        long m;
398        double crs , ex , eth;
399
400        if( n == 0 )
401        {
402                m = 1;
403        }
404        else
405        {
406                if( n == 1 )
407                {
408                        m = 2;
409                }
410                else
411                {
412                        m = n;
413                }
414        }
415
416        eth = ph1(0,0,iz-1,0)/POW2((double)m);
417        ex = MAX2(1. , e/eth );
418
419        crs = hpfit( iz , n , ex );
420        ASSERT( crs > 0. );
421
422        return crs;
423}
424
425double t_ADfA::hpfit(long int iz,
426                     long int n,
427                     double e)
428{
429        long int l,
430          m;
431        double cs,
432          eth,
433          ex,
434          q,
435          x;
436
437        DEBUG_ENTRY( "hpfit()" );
438
439        /*state specific photoionization cross sections for model hydrogen atom
440         * Version 1, September 23, 1997
441         ******************************************************************************
442         *** This subroutine calculates state-specific photoionization cross sections
443         *** for hydrogen and hydrogen-like ions.
444         *** Input parameters:  iz - atomic number
445         ***          n  - shell number, from 0 to 400:
446         ***                                    0 - 1s
447         ***                                    1 - 2s
448         ***                                    2 - 2p
449         ***                                    3 - 3
450         ***                                    ......
451         ***          e  - photon energy, eV
452         *** return value - cross section, cm^(-2)     
453         *******************************************************************************/
454
455        cs = 0.0;
456
457        ASSERT( iz > 0 && e>0. );
458
459        if( n >= NHYDRO_MAX_LEVEL )
460        {
461                fprintf( ioQQQ, " hpfit called with too large n, =%li\n" , n );
462                cdEXIT(EXIT_FAILURE);
463        }
464
465        l = 0;
466        if( n == 2 )
467        {
468                l = 1;
469        }
470        q = 3.5 + l - 0.5*PHH[n][1];
471
472        if( n == 0 )
473        {
474                m = 1;
475        }
476        else
477        {
478                if( n == 1 )
479                {
480                        m = 2;
481                }
482                else
483                {
484                        m = n;
485                }
486        }
487
488        eth = ph1(0,0,iz-1,0)/POW2((double)m);
489        ex = MAX2(1. , e/eth );
490
491        /* Don't just force to be at least one...make sure e/eth is close to one or greater.    */
492        ASSERT( e/eth > 0.95 );
493
494        if( ex < 1.0 )
495        {
496                return(0.);
497        }
498
499        x = ex/PHH[n][0];
500        cs = (PHH[n][4]*pow(1.0 + ((double)PHH[n][2]/x),(double)PHH[n][3])/
501          pow(x,q)/pow(1.0 + sqrt(x),(double)PHH[n][1])*8.79737e-17/
502          POW2((double)iz));
503        return(cs);
504}
505
506void t_ADfA::rec_lines(double t,
507                       realnum r[][471])
508{
509        long int i,
510          j,
511          ipj1,
512          ipj2;
513
514        double a,
515          a1,
516          dr[4][405],
517          p1,
518          p2,
519          p3,
520          p4,
521          p5,
522          p6,
523          rr[4][110],
524          te,
525          x,
526          z;
527
528        static long jd[6]={143,145,157,360,376,379};
529
530        static long ip[38]={7,9,12,13,14,16,18,19,20,21,22,44,45,49,50,
531          52,53,54,55,56,57,58,59,60,66,67,78,83,84,87,88,95,96,97,100,
532          101,103,104};
533
534        static long id[38]={7,3,10,27,23,49,51,64,38,47,60,103,101,112,
535          120,114,143,145,157,152,169,183,200,163,225,223,237,232,235,
536          249,247,300,276,278,376,360,379,384};
537
538        DEBUG_ENTRY( "rec_lines()" );
539
540        /*effective recombination coefficients for lines of C, N, O, by D. Verner
541         * Version 2, April 30, 1997
542         ******************************************************************************
543         *** This subroutine calculates effective recombination coefficients
544         *** for 110 permitted recombination lines of C, N, O (Pequignot, Petitjean,
545         *** & Boisson, 1991, A&A, 251, 680) and 405 permitted dielectronic
546         *** recombination lines (Nussbaumer & Storey, 1984, A&AS, 56, 293)
547         *** Input parameter:   t  - temperature, K
548         *** Output parameters: r(i,j), i=1,471
549         ***          r(i,1) - atomic number
550         ***          r(i,2) - number of electrons
551         ***          r(i,3) - wavelength, angstrom
552         ***          r(i,4) - rate coefficient, cm^3 s^(-1)
553         ****************************************************************************** */
554
555        for( i=0; i < 110; i++ )
556        {
557                rr[0][i] = P[0][i];
558                rr[1][i] = P[1][i];
559                rr[2][i] = P[2][i];
560                z = P[0][i] - P[1][i] + 1.0;
561                te = 1.0e-04*t/z/z;
562                p1 = P[3][i];
563                p2 = P[4][i];
564                p3 = P[5][i];
565                p4 = P[6][i];
566                if( te < 0.004 )
567                {
568                        a1 = p1*pow(0.004,p2)/(1.0 + p3*pow(0.004,p4));
569                        a = a1/sqrt(te/0.004);
570                }
571                else
572                {
573                        if( te > 2.0 )
574                        {
575                                a1 = p1*pow(2.0,p2)/(1.0 + p3*pow(2.0,p4));
576                                a = a1/pow(te/2.0,1.5);
577                        }
578                        else
579                        {
580                                a = p1*pow(te,p2)/(1.0 + p3*pow(te,p4));
581                        }
582                }
583                rr[3][i] = 1.0e-13*z*a*P[7][i];
584        }
585
586        for( i=0; i < 405; i++ )
587        {
588                dr[0][i] = ST[0][i];
589                dr[1][i] = ST[1][i];
590                dr[2][i] = ST[2][i];
591                te = 1.0e-04*t;
592                p1 = ST[3][i];
593                p2 = ST[4][i];
594                p3 = ST[5][i];
595                p4 = ST[6][i];
596                p5 = ST[7][i];
597                p6 = ST[8][i];
598                if( te < p6 )
599                {
600                        x = p5*(1.0/te - 1.0/p6);
601                        if( x > 80.0 )
602                        {
603                                a = 0.0;
604                        }
605                        else
606                        {
607                                a1 = (p1/p6 + p2 + p3*p6 + p4*p6*p6)/pow(p6,1.5)/exp(p5/
608                                  p6);
609                                a = a1/exp(x);
610                        }
611                }
612                else
613                {
614                        if( te > 6.0 )
615                        {
616                                a1 = (p1/6.0 + p2 + p3*6.0 + p4*36.0)/pow(6.0,1.5)/
617                                  exp(p5/6.0);
618                                a = a1/pow(te/6.0,1.5);
619                        }
620                        else
621                        {
622                                a = (p1/te + p2 + p3*te + p4*te*te)/pow(te,1.5)/exp(p5/
623                                  te);
624                        }
625                }
626                dr[3][i] = 1.0e-12*a;
627        }
628
629        for( i=0; i < 6; i++ )
630        {
631                ipj1 = jd[i];
632                ipj2 = ipj1 + 1;
633                dr[3][ipj1-1] += dr[3][ipj2-1];
634                dr[0][ipj2-1] = 0.0;
635        }
636
637        for( i=0; i < 38; i++ )
638        {
639                ipj1 = ip[i];
640                ipj2 = id[i];
641                rr[3][ipj1-1] += dr[3][ipj2-1];
642                dr[0][ipj2-1] = 0.0;
643        }
644
645        for( i=0; i < 110; i++ )
646        {
647                r[0][i] = (realnum)rr[0][i];
648                r[1][i] = (realnum)rr[1][i];
649                r[2][i] = (realnum)rr[2][i];
650                r[3][i] = (realnum)rr[3][i];
651        }
652
653        j = 110;
654        for( i=0; i < 405; i++ )
655        {
656                if( dr[0][i] > 1.0 )
657                {
658                        j += 1;
659                        r[0][j-1] = (realnum)dr[0][i];
660                        r[1][j-1] = (realnum)dr[1][i];
661                        r[2][j-1] = (realnum)dr[2][i];
662                        r[3][j-1] = (realnum)dr[3][i];
663                }
664        }
665        return;
666}
667
668double t_ADfA::rad_rec(long int iz,
669                       long int in,
670                       double t)
671{
672        /*
673         *** Version 4. June 29, 1999.
674         *** Written by D. A. Verner, verner@pa.uky.edu
675         ******************************************************************************
676         *** This subroutine calculates rates of radiative recombination for all ions
677         *** of all elements from H through Zn by use of the following fits:
678         *** H-like, He-like, Li-like, Na-like -
679         *** >>refer    all     reccoef Verner & Ferland, 1996, ApJS, 103, 467
680         *** Other ions of C, N, O, Ne - Pequignot et al. 1991, A&A, 251, 680,
681         ***    refitted by Verner & Ferland formula to ensure correct asymptotes
682         *** Fe XVII-XXIII -
683         *** >>refer    Fe17-23 recom   Arnaud & Raymond, 1992, ApJ, 398, 394
684         *** Fe I-XV - refitted by Verner & Ferland formula to ensure correct asymptotes
685         *** Other ions of Mg, Si, S, Ar, Ca, Fe, Ni -
686         ***                      -
687         *** >>refer    all     recom   Shull & Van Steenberg, 1982, ApJS, 48, 95
688         *** Other ions of Na, Al -
689         *** >>refer    Na, Al  recom   Landini & Monsignori Fossi, 1990, A&AS, 82, 229
690         *** Other ions of F, P, Cl, K, Ti, Cr, Mn, Co (excluding Ti I-II, Cr I-IV,
691         *** Mn I-V, Co I)        -
692         *** >>refer    many    recom   Landini & Monsignori Fossi, 1991, A&AS, 91, 183
693         *** All other species    - interpolations of the power-law fits
694         *** Input parameters:  iz - atomic number
695         ***                    in - number of electrons from 1 to iz
696         ***                    t  - temperature, K
697         *** return result:  - rate coefficient, cm^3 s^(-1)
698         ******************************************************************************
699         */
700        double tt;
701        double rate;
702
703        DEBUG_ENTRY( "rad_rec()" );
704
705        rate = 0.0;
706        if( iz < 1 || iz > 30 )
707        {
708                fprintf( ioQQQ, " rad_rec called with insane atomic number, =%4ld\n",
709                  iz );
710                cdEXIT(EXIT_FAILURE);
711        }
712        if( in < 1 || in > iz )
713        {
714                fprintf( ioQQQ, " rad_rec called with insane number elec =%4ld\n",
715                  in );
716                cdEXIT(EXIT_FAILURE);
717        }
718        if( (((in <= 3 || in == 11) || (iz > 5 && iz < 9)) || iz == 10) ||
719          (iz == 26 && in > 11) )
720        {
721                tt = sqrt(t/rnew[in-1][iz-1][2]);
722                rate =
723                  rnew[in-1][iz-1][0]/(tt*pow(tt + 1.0,1.0 - rnew[in-1][iz-1][1])*
724                  pow(1.0 + sqrt(t/rnew[in-1][iz-1][3]),1.0 + rnew[in-1][iz-1][1]));
725        }
726        else
727        {
728                tt = t*1.0e-04;
729                if( iz == 26 && in <= 13 )
730                {
731                        rate = fe[in-1][0]/pow(tt,fe[in-1][1] +
732                          fe[in-1][2]*log10(tt));
733                }
734                else
735                {
736                        rate = rrec[in-1][iz-1][0]/pow(tt,(double)rrec[in-1][iz-1][1]);
737                }
738        }
739
740        return rate;
741}
742
743double t_ADfA::H_rad_rec(long int iz,
744                         long int n,
745                         double t)
746{
747        /*
748         * Version 4, October 9, 1997
749         ******************************************************************************
750         *** This subroutine calculates state-specific recombination rates
751         *** for hydrogen and hydrogen-like ions.
752         *** Input parameters:  iz - atomic number
753         ***          n  - shell number, from 0 to 400:
754         ***                                    0 - 1s
755         ***                                    1 - 2s
756         ***                                    2 - 2p
757         ***                                    3 - 3
758         ***                                    ......
759         ***          t  - temperature, K
760         *** Output parameter:  r  - rate coefficient, cm^3 s^(-1)
761         *** If n is negative, the subroutine returns the total recombination
762         *** rate coefficient
763         ******************************************************************************
764         */
765        double rate,
766          TeScaled,
767          x,
768          x1,
769          x2;
770
771        DEBUG_ENTRY( "H_rad_rec()" );
772
773        rate = 0.0;
774
775        /* iz is charge, must be 1 or greater */
776        ASSERT( iz > 0 );
777
778        /* n is level number, must be less than dim or hydro vectors */
779        ASSERT( n < NHYDRO_MAX_LEVEL );
780
781        TeScaled = t/POW2((double)iz);
782
783        if( n < 0 )
784        {
785                x1 = sqrt(TeScaled/3.148);
786                x2 = sqrt(TeScaled/7.036e05);
787                rate = 7.982e-11/x1/pow(1.0 + x1,0.252)/pow(1.0 + x2,1.748);
788        }
789        else
790        {
791                x = log10(TeScaled);
792                rate = (HRF[n][0] + HRF[n][2]*x + HRF[n][4]*
793                  x*x + HRF[n][6]*powi(x,3) + HRF[n][8]*powi(x,4))/
794                  (1.0 + HRF[n][1]*x + HRF[n][3]*x*x + HRF[n][5]*
795                  powi(x,3) + HRF[n][7]*powi(x,4));
796                rate = pow(10.0,rate)/TeScaled;
797        }
798        rate *= iz;
799
800        return rate;
801}
802
803/*coll_ion D Verner's routine to compute collisional ionization rate coefficients,
804 * returns collisional ionization rate coefficient cm^3 s^-1*/
805double t_ADfA::coll_ion(
806        /* atomic number, 1 for hydrogen */
807        long int iz,
808        /* stage of ionization, 1 for atom */
809        long int in,
810        /* temperature */
811        double t)
812{
813        double rate, te, u;
814
815        DEBUG_ENTRY( "coll_ion()" );
816        /*D Verner's routine to compute collisional ionization rate coefficients
817         * Version 3, April 21, 1997
818         * Cu (Z=29) and Zn (Z=30) are added (fits from Ni, correct thresholds).
819         ******************************************************************************
820         *** This subroutine calculates rates of direct collisional ionization
821         *** for all ionization stages of all elements from H to Ni (Z=28)
822         *** by use of the fits from
823         *>>refer       all     collion Voronov, G. S., 1997, At. Data Nucl. Data Tables, 65, 1
824         *** Input parameters:  iz - atomic number on pphysical scale, H is 1
825         ***          in - number of electrons from 1 to iz
826         ***          t  - temperature, K
827         *** Output parameter:  c  - rate coefficient, cm^3 s^(-1)
828         ****************************************************************************** */
829        rate = 0.0;
830
831        if( iz < 1 || iz > 30 )
832        {
833                /* return zero rate is atomic number outside range of code */
834                return( 0. );
835        }
836
837        if( in < 1 || in > iz )
838        {
839                /* return zero rate is ion stage is impossible */
840                return( 0. );
841        }
842
843        te = t*EVRYD/TE1RYD;
844        u = CF[in-1][iz-1][0]/te;
845        if( u > 80.0 )
846        {
847                return( 0. );
848        }
849
850        rate = (CF[in-1][iz-1][2]*(1.0 + CF[in-1][iz-1][1]*
851          sqrt(u))/(CF[in-1][iz-1][3] + u)*pow(u,CF[in-1][iz-1][4])*
852          exp(-u));
853
854        return(rate);
855}
856
857realnum t_ADfA::h_coll_str( long ipLo, long ipHi, long ipTe )
858{
859        realnum rate;
860
861        DEBUG_ENTRY( "h_coll_str()" );
862
863        ASSERT( ipLo < ipHi );
864
865#       if !defined NDEBUG
866        long ipISO = ipH_LIKE;
867        long nelem = ipHYDROGEN;
868#       endif
869        ASSERT( N_(ipLo) < N_(ipHi) );
870        ASSERT( N_(ipHi) <= 5 );
871
872        rate = (realnum)HCS[ipHi-1][ipLo][ipTe];
873
874        return rate;
875}
Note: See TracBrowser for help on using the browser.