cloudy: /home66/gary/public_html/cloudy/c13_branch/source/atmdat

00001 /* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
00002  * others.  For conditions of distribution and use see copyright notice in license.txt */
00003 /*phfit derive photoionization cross sections for first 30 elements */
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "atmdat.h"
00007 #include "iso.h"
00008 #include "taulines.h"
00009 #include "thirdparty.h"
00010 
00012 t_ADfA::t_ADfA()
00013 {
00014         DEBUG_ENTRY( "t_ADfA()" );
00015 
00016         /* this option, use the new atmdat_rad_rec recombination rates */
00017         version = PHFIT_UNDEF;
00018 
00019         double help[9];
00020         const long VERSION_MAGIC = 20061204L;
00021 
00022         static const char chFile[] = "phfit.dat";
00023 
00024         FILE *io = open_data( chFile, "r" );
00025 
00026         bool lgErr = false;
00027         long i=-1, j=-1, k=-1, n;
00028 
00029         lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00030         if( lgErr || i != VERSION_MAGIC )
00031         {
00032                 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile, i );
00033                 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00034                 cdEXIT(EXIT_FAILURE);
00035         }
00036 
00037         for( n=0; n < 7; n++ )
00038                 lgErr = lgErr || ( fscanf( io, "%ld", &L[n] ) != 1 );
00039         for( n=0; n < 30; n++ )
00040                 lgErr = lgErr || ( fscanf( io, "%ld", &NINN[n] ) != 1 );
00041         for( n=0; n < 30; n++ )
00042                 lgErr = lgErr || ( fscanf( io, "%ld", &NTOT[n] ) != 1 );
00043         while( true )
00044         {
00045                 lgErr = lgErr || ( fscanf( io, "%ld %ld %ld", &i, &j, &k ) != 3 );
00046                 if( i == -1 && j == -1 && k == -1 )
00047                         break;
00048                 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le", &help[0], &help[1],
00049                                            &help[2], &help[3], &help[4], &help[5] ) != 6 );
00050                 for( int l=0; l < 6; ++l )
00051                         PH1[i][j][k][l] = (realnum)help[l];
00052         }
00053         while( true )
00054         {
00055                 lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
00056                 if( i == -1 && j == -1 )
00057                         break;
00058                 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le", &help[0], &help[1],
00059                                            &help[2], &help[3], &help[4], &help[5], &help[6] ) != 7 );
00060                 for( int l=0; l < 7; ++l )
00061                         PH2[i][j][l]  = (realnum)help[l];
00062         }
00063         fclose( io );
00064 
00065         ASSERT( !lgErr );
00066 
00067         static const char chFile2[] = "hpfit.dat";
00068 
00069         io = open_data( chFile2, "r" );
00070 
00071         lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00072         if( lgErr || i != VERSION_MAGIC )
00073         {
00074                 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile2, i );
00075                 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00076                 cdEXIT(EXIT_FAILURE);
00077         }
00078 
00079         for( i=0; i < NHYDRO_MAX_LEVEL; i++ )
00080         {
00081                 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le", &help[0], &help[1],
00082                                            &help[2], &help[3], &help[4] ) != 5 );
00083                 for( int l=0; l < 5; ++l )
00084                         PHH[i][l] = (realnum)help[l];
00085         }
00086 
00087         fclose( io );
00088 
00089         ASSERT( !lgErr );
00090 
00091         static const char chFile3[] = "rec_lines.dat";
00092 
00093         io = open_data( chFile3, "r" );
00094 
00095         lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00096         if( lgErr || i != VERSION_MAGIC )
00097         {
00098                 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile3, i );
00099                 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00100                 cdEXIT(EXIT_FAILURE);
00101         }
00102 
00103         for( i=0; i < 110; i++ )
00104         {
00105                 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2],
00106                                            &help[3], &help[4], &help[5], &help[6], &help[7] ) != 8 );
00107                 for( int l=0; l < 8; ++l )
00108                         P[l][i] = (realnum)help[l];
00109         }
00110 
00111 
00112         for( i=0; i < 405; i++ )
00113         {
00114                 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2],
00115                                            &help[3], &help[4], &help[5], &help[6], &help[7], &help[8] ) != 9 );
00116                 for( int l=0; l < 9; ++l )
00117                         ST[l][i] = (realnum)help[l];
00118         }
00119 
00120         fclose( io );
00121 
00122         ASSERT( !lgErr );
00123 
00124         static const char chFile4[] = "rad_rec.dat";
00125 
00126         io = open_data( chFile4, "r" );
00127 
00128         lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00129         if( lgErr || i != VERSION_MAGIC )
00130         {
00131                 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile4, i );
00132                 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00133                 cdEXIT(EXIT_FAILURE);
00134         }
00135 
00136         while( true )
00137         {
00138                 lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
00139                 if( i == -1 && j == -1 )
00140                         break;
00141                 lgErr = lgErr || ( fscanf( io, "%le %le", &help[0], &help[1] ) != 2 );
00142                 for( int l=0; l < 2; ++l )
00143                         rrec[i][j][l] = (realnum)help[l];
00144         }
00145         while( true )
00146         {
00147                 lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
00148                 if( i == -1 && j == -1 )
00149                         break;
00150                 lgErr = lgErr || ( fscanf( io, "%le %le %le %le", &help[0], &help[1],
00151                                            &help[2], &help[3] ) != 4 );
00152                 for( int l=0; l < 4; ++l )
00153                         rnew[i][j][l] = (realnum)help[l];
00154         }
00155         while( true )
00156         {
00157                 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00158                 if( i == -1 )
00159                         break;
00160                 lgErr = lgErr || ( fscanf( io, "%le %le %le", &help[0], &help[1], &help[2] ) != 3 );
00161                 for( int l=0; l < 3; ++l )
00162                         fe[i][l] = (realnum)help[l];
00163         }
00164 
00165         fclose( io );
00166 
00167         ASSERT( !lgErr );
00168 
00169         static const char chFile5[] = "h_rad_rec.dat";
00170 
00171         io = open_data( chFile5, "r" );
00172 
00173         lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00174         if( lgErr || i != VERSION_MAGIC )
00175         {
00176                 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile5, i );
00177                 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00178                 cdEXIT(EXIT_FAILURE);
00179         }
00180 
00181         for( i=0; i < NHYDRO_MAX_LEVEL; i++ )
00182         {
00183                 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2],
00184                                            &help[3], &help[4], &help[5], &help[6], &help[7], &help[8] ) != 9 );
00185                 for( int l=0; l < 9; ++l )
00186                         HRF[i][l] = (realnum)help[l];
00187         }
00188 
00189         fclose( io );
00190 
00191         ASSERT( !lgErr );
00192 
00193         static const char chFile6[] = "h_phot_cs.dat";
00194 
00195         io = open_data( chFile6, "r" );
00196 
00197         lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00198         if( lgErr || i != VERSION_MAGIC )
00199         {
00200                 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile6, i );
00201                 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00202                 cdEXIT(EXIT_FAILURE);
00203         }
00204 
00205         for( i=0; i < NHYDRO_MAX_LEVEL; i++ )
00206         {
00207                 lgErr = lgErr || ( fscanf( io, "%le", &help[0] ) != 1 );
00208                 STH[i] = (realnum)help[0];
00209         }
00210 
00211         fclose( io );
00212 
00213         ASSERT( !lgErr );
00214 
00215         static const char chFile7[] = "coll_ion.dat";
00216 
00217         io = open_data( chFile7, "r" );
00218 
00219         lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00220         if( lgErr || i != VERSION_MAGIC )
00221         {
00222                 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile7, i );
00223                 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00224                 cdEXIT(EXIT_FAILURE);
00225         }
00226 
00227         while( true )
00228         {
00229                 lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
00230                 if( i == -1 && j == -1 )
00231                         break;
00232                 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le", &CF[i][j][0], &CF[i][j][1],
00233                                            &CF[i][j][2], &CF[i][j][3], &CF[i][j][4] ) != 5 );
00234         }
00235 
00236         fclose( io );
00237 
00238         ASSERT( !lgErr );
00239 
00240         /*refer HI      cs      Anderson, H., Ballance, C.P., Badnell, N.R., 
00241          *refercon      & Summers, H.P  2000, J Phys B, 33, 1255 */
00242         static const char chFile8[] = "h_coll_str.dat";
00243 
00244         io = open_data( chFile8, "r" );
00245 
00246         lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00247         if( lgErr || i != VERSION_MAGIC )
00248         {
00249                 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile8, i );
00250                 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00251                 cdEXIT(EXIT_FAILURE);
00252         }
00253 
00254         while( true )
00255         {
00256                 lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
00257                 if( i == -1 && j == -1 )
00258                         break;
00259                 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le", &HCS[i-2][j-1][0], &HCS[i-2][j-1][1],
00260                                            &HCS[i-2][j-1][2], &HCS[i-2][j-1][3], &HCS[i-2][j-1][4], &HCS[i-2][j-1][5],
00261                                            &HCS[i-2][j-1][6], &HCS[i-2][j-1][7] ) != 8 );
00262         }
00263 
00264         fclose( io );
00265 
00266         ASSERT( !lgErr );
00267 }
00268 
00269 double t_ADfA::phfit(long int nz, 
00270                      long int ne,
00271                      long int is, 
00272                      double e)
00273 {
00274         long int nint, 
00275           nout;
00276         double a, 
00277           b, 
00278           crs, 
00279           einn, 
00280           p1, 
00281           q, 
00282           x, 
00283           y, 
00284           z;
00285 
00286         DEBUG_ENTRY( "phfit()" );
00287 
00288         /*** Version 3. October 8, 1996.
00289          *** Written by D. A. Verner, verner@pa.uky.edu
00290          *** Inner-shell ionization energies of some low-ionized species are slightly
00291          *** improved to fit smoothly the experimental inner-shell ionization energies 
00292          *** of neutral atoms.
00293          ******************************************************************************
00294          *** This subroutine calculates partial photoionization cross sections
00295          *** for all ionization stages of all atoms from H to Zn (Z=30) by use of
00296          *** the following fit parameters:
00297          *** Outer shells of the Opacity Project (OP) elements:
00298          *** >>refer    all     photo_cs        Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487.
00299          *** Inner shells of all elements, and outer shells of the non-OP elements:
00300          ***  Verner and Yakovlev, 1995, A&AS, 109, 125
00301          *** Input parameters:  nz - atomic number from 1 to 30 (integer) 
00302          ***          ne - number of electrons from 1 to iz (integer)
00303          ***          is - shell number (integer)
00304          ***          e - photon energy, eV 
00305          ***          version - enum, PHFIT96 (default): calculates 
00306          ***                 new cross sections, PHFIT95: calculates
00307          ***                 only old Hartree-Slater cross sections
00308          *** Output parameter:  photoionization cross section, Mb
00309          *** Shell numbers:
00310          *** 1 - 1s, 2 - 2s, 3 - 2p, 4 - 3s, 5 - 3p, 6 - 3d, 7 - 4s. 
00311          *** If a species in the ground state has no electrons on the given shell,
00312          *** the subroutine returns 0.
00313          ****************************************************************************** */
00314 
00315         crs = 0.0;
00316         if( nz < 1 || nz > 30 )
00317         { 
00318                 return(crs);
00319         }
00320 
00321         if( ne < 1 || ne > nz )
00322         { 
00323                 return(crs);
00324         }
00325 
00326         nout = NTOT[ne-1];
00327         if( nz == ne && nz > 18 )
00328                 nout = 7;
00329         if( nz == (ne + 1) && ((((nz == 20 || nz == 21) || nz == 22) || 
00330           nz == 25) || nz == 26) )
00331                 nout = 7;
00332         if( is > nout )
00333         { 
00334                 return(crs);
00335         }
00336 
00337         if( (is == 6 && (nz == 20 || nz == 19)) && ne >= 19 )
00338         { 
00339                 return(crs);
00340         }
00341 
00342         ASSERT( is >= 1 && is <= 7 );
00343 
00344         if( e < PH1[is-1][ne-1][nz-1][0] )
00345         { 
00346                 return(crs);
00347         }
00348 
00349         nint = NINN[ne-1];
00350         if( ((nz == 15 || nz == 17) || nz == 19) || (nz > 20 && nz != 26) )
00351         {
00352                 einn = 0.0;
00353         }
00354         else
00355         {
00356                 if( ne < 3 )
00357                 {
00358                         einn = 1.0e30;
00359                 }
00360                 else
00361                 {
00362                         einn = PH1[nint-1][ne-1][nz-1][0];
00363                 }
00364         }
00365 
00366         if( (is <= nint || e >= einn) || version == PHFIT95 )
00367         {
00368                 p1 = -PH1[is-1][ne-1][nz-1][4];
00369                 y = e/PH1[is-1][ne-1][nz-1][1];
00370                 q = -0.5*p1 - L[is-1] - 5.5;
00371                 a = PH1[is-1][ne-1][nz-1][2]*(POW2(y - 1.0) + 
00372                         POW2(PH1[is-1][ne-1][nz-1][5]));
00373                 b = sqrt(y/PH1[is-1][ne-1][nz-1][3]) + 1.0;
00374                 crs = a*pow(y,q)*pow(b,p1);
00375         }
00376         else
00377         {
00378                 if( (is < nout && is > nint) && e < einn )
00379                 { 
00380                         return(crs);
00381                 }
00382                 p1 = -PH2[ne-1][nz-1][3];
00383                 q = -0.5*p1 - 5.5;
00384                 x = e/PH2[ne-1][nz-1][0] - PH2[ne-1][nz-1][5];
00385                 z = sqrt(x*x+POW2(PH2[ne-1][nz-1][6]));
00386                 a = PH2[ne-1][nz-1][1]*(POW2(x - 1.0) + 
00387                         POW2(PH2[ne-1][nz-1][4]));
00388                 b = 1.0 + sqrt(z/PH2[ne-1][nz-1][2]);
00389                 crs = a*pow(z,q)*pow(b,p1);
00390         }
00391         return(crs);
00392 }
00393 
00394 double t_ADfA::hpfit(long int iz,
00395                      long int n,
00396                      double e)
00397 {
00398         long int l, 
00399           m;
00400         double cs,
00401           eth, 
00402           ex, 
00403           q, 
00404           x;
00405 
00406         DEBUG_ENTRY( "hpfit()" );
00407 
00408         /*state specific photoionization cross sections for model hydrogen atom
00409          * Version 1, September 23, 1997
00410          ******************************************************************************
00411          *** This subroutine calculates state-specific photoionization cross sections
00412          *** for hydrogen and hydrogen-like ions.
00413          *** Input parameters:  iz - atomic number 
00414          ***          n  - shell number, from 0 to 400:
00415          ***                                    0 - 1s
00416          ***                                    1 - 2s
00417          ***                                    2 - 2p
00418          ***                                    3 - 3 
00419          ***                                    ......
00420          ***          e  - photon energy, eV
00421          *** return value - cross section, cm^(-2)     
00422          *******************************************************************************/
00423 
00424         cs = 0.0;
00425 
00426         ASSERT( iz > 0 && e>0. );
00427 
00428         if( n >= NHYDRO_MAX_LEVEL )
00429         { 
00430                 fprintf( ioQQQ, " hpfit called with too large n, =%li\n" , n );
00431                 cdEXIT(EXIT_FAILURE);
00432         }
00433 
00434         l = 0;
00435         if( n == 2 )
00436         {
00437                 l = 1;
00438         }
00439         q = 3.5 + l - 0.5*PHH[n][1];
00440 
00441         if( n == 0 )
00442         {
00443                 m = 1;
00444         }
00445         else
00446         {
00447                 if( n == 1 )
00448                 {
00449                         m = 2;
00450                 }
00451                 else
00452                 {
00453                         m = n;
00454                 }
00455         }
00456 
00457         eth = ph1(0,0,iz-1,0)/POW2((double)m);
00458         ex = MAX2(1. , e/eth );
00459 
00460         /* Don't just force to be at least one...make sure e/eth is close to one or greater.    */
00461         ASSERT( e/eth > 0.95 );
00462 
00463         if( ex < 1.0 )
00464         { 
00465                 return(0.);
00466         }
00467 
00468         x = ex/PHH[n][0];
00469         cs = (PHH[n][4]*pow(1.0 + ((double)PHH[n][2]/x),(double)PHH[n][3])/
00470           pow(x,q)/pow(1.0 + sqrt(x),(double)PHH[n][1])*8.79737e-17/
00471           POW2((double)iz));
00472         return(cs);
00473 }
00474 
00475 void t_ADfA::rec_lines(double t, 
00476                        realnum r[][471])
00477 {
00478         long int i, 
00479           j, 
00480           ipj1, 
00481           ipj2;
00482 
00483         double a, 
00484           a1, 
00485           dr[4][405], 
00486           p1, 
00487           p2, 
00488           p3, 
00489           p4, 
00490           p5, 
00491           p6, 
00492           rr[4][110], 
00493           te, 
00494           x, 
00495           z;
00496 
00497         static long jd[6]={143,145,157,360,376,379};
00498 
00499         static long ip[38]={7,9,12,13,14,16,18,19,20,21,22,44,45,49,50,
00500           52,53,54,55,56,57,58,59,60,66,67,78,83,84,87,88,95,96,97,100,
00501           101,103,104};
00502 
00503         static long id[38]={7,3,10,27,23,49,51,64,38,47,60,103,101,112,
00504           120,114,143,145,157,152,169,183,200,163,225,223,237,232,235,
00505           249,247,300,276,278,376,360,379,384};
00506 
00507         DEBUG_ENTRY( "rec_lines()" );
00508 
00509         /*effective recombination coefficients for lines of C, N, O, by D. Verner
00510          * Version 2, April 30, 1997
00511          ******************************************************************************
00512          *** This subroutine calculates effective recombination coefficients
00513          *** for 110 permitted recombination lines of C, N, O (Pequignot, Petitjean,
00514          *** & Boisson, 1991, A&A, 251, 680) and 405 permitted dielectronic 
00515          *** recombination lines (Nussbaumer & Storey, 1984, A&AS, 56, 293)
00516          *** Input parameter:   t  - temperature, K
00517          *** Output parameters: r(i,j), i=1,471
00518          ***          r(i,1) - atomic number
00519          ***          r(i,2) - number of electrons
00520          ***          r(i,3) - wavelength, angstrom
00521          ***          r(i,4) - rate coefficient, cm^3 s^(-1)
00522          ****************************************************************************** */
00523 
00524         for( i=0; i < 110; i++ )
00525         {
00526                 rr[0][i] = P[0][i];
00527                 rr[1][i] = P[1][i];
00528                 rr[2][i] = P[2][i];
00529                 z = P[0][i] - P[1][i] + 1.0;
00530                 te = 1.0e-04*t/z/z;
00531                 p1 = P[3][i];
00532                 p2 = P[4][i];
00533                 p3 = P[5][i];
00534                 p4 = P[6][i];
00535                 if( te < 0.004 )
00536                 {
00537                         a1 = p1*pow(0.004,p2)/(1.0 + p3*pow(0.004,p4));
00538                         a = a1/sqrt(te/0.004);
00539                 }
00540                 else
00541                 {
00542                         if( te > 2.0 )
00543                         {
00544                                 a1 = p1*pow(2.0,p2)/(1.0 + p3*pow(2.0,p4));
00545                                 a = a1/pow(te/2.0,1.5);
00546                         }
00547                         else
00548                         {
00549                                 a = p1*pow(te,p2)/(1.0 + p3*pow(te,p4));
00550                         }
00551                 }
00552                 rr[3][i] = 1.0e-13*z*a*P[7][i];
00553         }
00554 
00555         for( i=0; i < 405; i++ )
00556         {
00557                 dr[0][i] = ST[0][i];
00558                 dr[1][i] = ST[1][i];
00559                 dr[2][i] = ST[2][i];
00560                 te = 1.0e-04*t;
00561                 p1 = ST[3][i];
00562                 p2 = ST[4][i];
00563                 p3 = ST[5][i];
00564                 p4 = ST[6][i];
00565                 p5 = ST[7][i];
00566                 p6 = ST[8][i];
00567                 if( te < p6 )
00568                 {
00569                         x = p5*(1.0/te - 1.0/p6);
00570                         if( x > 80.0 )
00571                         {
00572                                 a = 0.0;
00573                         }
00574                         else
00575                         {
00576                                 a1 = (p1/p6 + p2 + p3*p6 + p4*p6*p6)/pow(p6,1.5)/exp(p5/
00577                                   p6);
00578                                 a = a1/exp(x);
00579                         }
00580                 }
00581                 else
00582                 {
00583                         if( te > 6.0 )
00584                         {
00585                                 a1 = (p1/6.0 + p2 + p3*6.0 + p4*36.0)/pow(6.0,1.5)/
00586                                   exp(p5/6.0);
00587                                 a = a1/pow(te/6.0,1.5);
00588                         }
00589                         else
00590                         {
00591                                 a = (p1/te + p2 + p3*te + p4*te*te)/pow(te,1.5)/exp(p5/
00592                                   te);
00593                         }
00594                 }
00595                 dr[3][i] = 1.0e-12*a;
00596         }
00597 
00598         for( i=0; i < 6; i++ )
00599         {
00600                 ipj1 = jd[i];
00601                 ipj2 = ipj1 + 1;
00602                 dr[3][ipj1-1] += dr[3][ipj2-1];
00603                 dr[0][ipj2-1] = 0.0;
00604         }
00605 
00606         for( i=0; i < 38; i++ )
00607         {
00608                 ipj1 = ip[i];
00609                 ipj2 = id[i];
00610                 rr[3][ipj1-1] += dr[3][ipj2-1];
00611                 dr[0][ipj2-1] = 0.0;
00612         }
00613 
00614         for( i=0; i < 110; i++ )
00615         {
00616                 r[0][i] = (realnum)rr[0][i];
00617                 r[1][i] = (realnum)rr[1][i];
00618                 r[2][i] = (realnum)rr[2][i];
00619                 r[3][i] = (realnum)rr[3][i];
00620         }
00621 
00622         j = 110;
00623         for( i=0; i < 405; i++ )
00624         {
00625                 if( dr[0][i] > 1.0 )
00626                 {
00627                         j += 1;
00628                         r[0][j-1] = (realnum)dr[0][i];
00629                         r[1][j-1] = (realnum)dr[1][i];
00630                         r[2][j-1] = (realnum)dr[2][i];
00631                         r[3][j-1] = (realnum)dr[3][i];
00632                 }
00633         }
00634         return;
00635 }
00636 
00637 double t_ADfA::rad_rec(long int iz, 
00638                        long int in, 
00639                        double t)
00640 {
00641         /*
00642          *** Version 4. June 29, 1999.
00643          *** Written by D. A. Verner, verner@pa.uky.edu 
00644          ******************************************************************************
00645          *** This subroutine calculates rates of radiative recombination for all ions
00646          *** of all elements from H through Zn by use of the following fits:
00647          *** H-like, He-like, Li-like, Na-like - 
00648          *** >>refer    all     rec_coef        Verner, D. A., & Ferland, G. J. 1996, ApJS, 103, 467
00649          *** Other ions of C, N, O, Ne - Pequignot et al. 1991, A&A, 251, 680,
00650          ***    refitted by Verner & Ferland formula to ensure correct asymptotes
00651          *** Fe XVII-XXIII - 
00652          *** >>refer    Fe17-23 recom   Arnaud, M., & Raymond, J. 1992, ApJ, 398, 394
00653          *** Fe I-XV - refitted by Verner & Ferland formula to ensure correct asymptotes
00654          *** Other ions of Mg, Si, S, Ar, Ca, Fe, Ni - 
00655          ***                      -
00656          *** >>refer    all     recom   Shull, M. J., & Van Steenberg, M. 1982, ApJS, 48, 95
00657          *** Other ions of Na, Al - 
00658          *** >>refer    Na, Al  recom   Landini, M., & Monsignori-Fossi, B.C. 1990, A&AS, 82, 229
00659          *** Other ions of F, P, Cl, K, Ti, Cr, Mn, Co (excluding Ti I-II, Cr I-IV,
00660          *** Mn I-V, Co I)        - 
00661          *** >>refer    many    recom   Landini, M., & Monsignori-Fossi, B.C. 1991, A&AS, 91, 183
00662          *** All other species    - interpolations of the power-law fits
00663          *** Input parameters:  iz - atomic number 
00664          ***                    in - number of electrons from 1 to iz 
00665          ***                    t  - temperature, K
00666          *** return result:  - rate coefficient, cm^3 s^(-1)
00667          ******************************************************************************
00668          */
00669         double tt;
00670         double rate;
00671 
00672         DEBUG_ENTRY( "rad_rec()" );
00673 
00674         rate = 0.0;
00675         if( iz < 1 || iz > 30 )
00676         {
00677                 fprintf( ioQQQ, " rad_rec called with insane atomic number, =%4ld\n", 
00678                   iz );
00679                 cdEXIT(EXIT_FAILURE);
00680         }
00681         if( in < 1 || in > iz )
00682         {
00683                 fprintf( ioQQQ, " rad_rec called with insane number elec =%4ld\n", 
00684                   in );
00685                 cdEXIT(EXIT_FAILURE);
00686         }
00687         if( (((in <= 3 || in == 11) || (iz > 5 && iz < 9)) || iz == 10) || 
00688           (iz == 26 && in > 11) )
00689         {
00690                 tt = sqrt(t/rnew[in-1][iz-1][2]);
00691                 rate = 
00692                   rnew[in-1][iz-1][0]/(tt*pow(tt + 1.0,1.0 - rnew[in-1][iz-1][1])*
00693                   pow(1.0 + sqrt(t/rnew[in-1][iz-1][3]),1.0 + rnew[in-1][iz-1][1]));
00694         }
00695         else
00696         {
00697                 tt = t*1.0e-04;
00698                 if( iz == 26 && in <= 13 )
00699                 {
00700                         rate = fe[in-1][0]/pow(tt,fe[in-1][1] + 
00701                           fe[in-1][2]*log10(tt));
00702                 }
00703                 else
00704                 {
00705                         rate = rrec[in-1][iz-1][0]/pow(tt,(double)rrec[in-1][iz-1][1]);
00706                 }
00707         }
00708 
00709         return rate;
00710 }
00711 
00712 double t_ADfA::H_rad_rec(long int iz,
00713                          long int n,
00714                          double t)
00715 {
00716         /*
00717          * Version 4, October 9, 1997
00718          ******************************************************************************
00719          *** This subroutine calculates state-specific recombination rates 
00720          *** for hydrogen and hydrogen-like ions.
00721          *** Input parameters:  iz - atomic number 
00722          ***          n  - shell number, from 0 to 400:
00723          ***                                    0 - 1s
00724          ***                                    1 - 2s
00725          ***                                    2 - 2p
00726          ***                                    3 - 3 
00727          ***                                    ......
00728          ***          t  - temperature, K
00729          *** Output parameter:  r  - rate coefficient, cm^3 s^(-1)
00730          *** If n is negative, the subroutine returns the total recombination 
00731          *** rate coefficient
00732          ******************************************************************************
00733          */
00734         double rate,
00735           TeScaled, 
00736           x, 
00737           x1, 
00738           x2;
00739 
00740         DEBUG_ENTRY( "H_rad_rec()" );
00741 
00742         rate = 0.0;
00743 
00744         /* iz is charge, must be 1 or greater */
00745         ASSERT( iz > 0 );
00746 
00747         /* n is level number, must be less than dim or hydro vectors */
00748         ASSERT( n < NHYDRO_MAX_LEVEL );
00749 
00750         TeScaled = t/POW2((double)iz);
00751 
00752         if( n < 0 )
00753         {
00754                 x1 = sqrt(TeScaled/3.148);
00755                 x2 = sqrt(TeScaled/7.036e05);
00756                 rate = 7.982e-11/x1/pow(1.0 + x1,0.252)/pow(1.0 + x2,1.748);
00757         }
00758         else
00759         {
00760                 x = log10(TeScaled);
00761                 rate = (HRF[n][0] + HRF[n][2]*x + HRF[n][4]*
00762                   x*x + HRF[n][6]*powi(x,3) + HRF[n][8]*powi(x,4))/
00763                   (1.0 + HRF[n][1]*x + HRF[n][3]*x*x + HRF[n][5]*
00764                   powi(x,3) + HRF[n][7]*powi(x,4));
00765                 rate = pow(10.0,rate)/TeScaled;
00766         }
00767         rate *= iz;
00768 
00769         return rate;
00770 }
00771 
00772 /*coll_ion D Verner's routine to compute collisional ionization rate coefficients,
00773  * returns collisional ionization rate coefficient cm^3 s^-1*/
00774 double t_ADfA::coll_ion(
00775         /* atomic number, 1 for hydrogen */
00776         long int iz, 
00777         /* number of bound electrons before ionization*/
00778         long int in, 
00779         /* temperature */
00780         double t)
00781 {
00782         double rate, te, u;
00783 
00784         DEBUG_ENTRY( "coll_ion()" );
00785         /*D Verner's routine to compute collisional ionization rate coefficients
00786          * Version 3, April 21, 1997
00787          * Cu (Z=29) and Zn (Z=30) are added (fits from Ni, correct thresholds).
00788          ******************************************************************************
00789          *** This subroutine calculates rates of direct collisional ionization 
00790          *** for all ionization stages of all elements from H to Ni (Z=28)
00791          *** by use of the fits from
00792          *>>refer       all     coll_ion        Voronov, G. S. 1997, At. Data Nucl. Data Tables, 65, 1
00793          *** Input parameters:  iz - atomic number on pphysical scale, H is 1
00794          ***          in - number of electrons from 1 to iz 
00795          ***          t  - temperature, K
00796          *** Output parameter:  c  - rate coefficient, cm^3 s^(-1)
00797          ****************************************************************************** */
00798         rate = 0.0;
00799 
00800         te = t*EVRYD/TE1RYD;
00801         u = CF[in-1][iz-1][0]/te;
00802         if( u > 80.0 )
00803         { 
00804                 return( 0. );
00805         }
00806 
00807         rate = (CF[in-1][iz-1][2]*(1.0 + CF[in-1][iz-1][1]*
00808           sqrt(u))/(CF[in-1][iz-1][3] + u)*pow(u,CF[in-1][iz-1][4])*
00809           exp(-u));
00810 
00811         return(rate);
00812 }
00813 
00814 double t_ADfA::coll_ion_wrapper(
00815         /* (atomic number - 1), 0 for hydrogen */
00816         long int z,
00817         /* stage of ionization, 0 for atom */
00818         long int n,
00819         /* temperature K */
00820         double t)
00821 {
00822         double rate = 0.0;
00823 
00824         DEBUG_ENTRY( "coll_ion_wrapper()" );
00825 
00826         if( z < 0 || z > LIMELM-1 )
00827         {
00828                 /* return zero rate is atomic number outside range of code */
00829                 return( 0. );
00830         }
00831 
00832         if( n < 0 || n > z )
00833         {
00834                 /* return zero rate is ion stage is impossible */
00835                 return( 0. );
00836         }
00837 
00838         /*Get the rate coefficients from either Dima or Hybrid.
00839          * The enum variable CIRCData is set to a default in init_defaults_preparse
00840          * It can be changed parse_set via the input command set collisional ionization XXX */
00841         if( atmdat.CIRCData == t_atmdat::DIMA)
00842         {
00843                 /* collisional ionization by thermal electrons
00844                 * >>chng 97 mar 19, to Dima's new routine using
00845                 * >>refer       all     coll_ion        Voronov, G. S. 1997, At. Data Nucl. Data Tables, 65, 1
00846                 *coll_ion(atomic number, number of electrons, temperature)*/
00847                 rate = coll_ion( z+1, z+1-n, t );
00848         }
00849         else if( atmdat.CIRCData == t_atmdat::HYBRID)
00850         {
00851                 // Get the rate coefficient by scaling Dima to Dere.
00852                 rate = coll_ion_hybrid( z, n, t );
00853         }
00854         else
00855         {
00856                 //CIRCData is an enum so it must be equal to one of the enum values.
00857                 TotalInsanity();
00858         }
00859 
00860         ASSERT(rate >= 0.0);
00861         return(rate);
00862 }
00863 
00864 /*coll_ion D Verner's routine to compute collisional ionization rate coefficients,
00865  * returns collisional ionization rate coefficient cm^3 s^-1*/
00866 double t_ADfA::coll_ion_hybrid(
00867         /* (atomic number - 1)  0 for hydrogen */
00868         long int nelem,
00869         /* stage of ionization, 0 for atom */
00870         long int ion,
00871         /* temperature */
00872         double t)
00873 {
00874 
00875         DEBUG_ENTRY( "coll_ion_hybrid()" );
00876 
00880         static const double DereRatio[30][30]= 
00881                 {{0.9063,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
00882                 {1.0389,1.0686,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
00883                 {0.6466,1.0772,0.963,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
00884                 {0.9715,0.7079,0.973,0.9899,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
00885                 {1.0781,1.3336,1.0032,1.1102,0.9874,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
00886                 {1.0499,0.913,1.0377,1.299,1.2874,1.0656,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
00887                 {1.0421,1.1966,1.0842,0.9619,1.0583,1.155,1.0648,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
00888                 {1.041,1.1181,1.0164,1.0271,0.9789,0.6829,1.1805,1.0672,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
00889                 {0.2636,0.658,1.0431,1.0749,0.894,1.059,0.9888,1.1426,1.0633,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
00890                 {0.8089,1.1395,1.1041,1.0449,1.1365,1.089,1.0147,0.9135,1.336,1.0429,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
00891                 {0.9545,1.1302,1.1105,1.2067,1.2569,1.079,0.8866,0.9924,0.9933,1.0488,1.0602,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
00892                 {0.3793,0.9857,1.1152,1.2826,1.3006,1.2416,1.0893,0.93,0.9953,0.9992,1.3479,1.0622,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
00893                 {0.7458,0.9975,1.0251,0.9553,1.5023,1.2136,1.2566,1.2417,1.2381,1.3755,1.1113,1.1629,1.0589,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
00894                 {0.7328,0.8798,0.4492,0.8221,1.7874,1.5418,1.5777,1.3036,1.124,1.0667,1.0009,1.002,1.1284,1.0673,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
00895                 {0.7075,0.9805,0.9451,0.5937,1.2061,1.1772,1.3818,1.4073,1.2643,1.1095,0.9903,1.3716,1.0058,1.0966,1.0517,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
00896                 {1.3572,0.8925,0.8119,0.8991,0.6633,1.4519,1.2073,1.4058,1.4519,1.2726,1.1428,0.9818,1.3643,1.0074,1.1836,1.1005,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
00897                 {0.5412,0.8428,0.9237,0.819,0.9151,0.7909,1.7424,1.2653,1.4513,1.4871,1.2633,1.1373,0.9969,0.9919,1.0091,1.2077,1.105,0,0,0,0,0,0,0,0,0,0,0,0,0},
00898                 {0.9242,0.8644,0.9752,0.8562,0.6998,0.6273,0.8686,2.0326,1.0364,1.4687,1.4812,1.2701,1.1376,0.9932,1.0001,1.0084,1.2036,1.1074,0,0,0,0,0,0,0,0,0,0,0,0},
00899                 {0.6381,0.9604,1.3556,0.9364,0.9678,1.0604,1.3067,1.0434,2.0722,0.979,1.5318,1.4968,1.2713,1.1427,1.0036,0.9953,1.03,1.2071,1.1083,0,0,0,0,0,0,0,0,0,0,0},
00900                 {0.7652,1.1668,1.0422,0.8705,0.9193,0.9982,1.077,1.3875,1.1544,2.4653,1.3188,1.5375,1.5307,1.2776,1.1427,1.0117,0.9872,1.0118,1.1899,1.1119,0,0,0,0,0,0,0,0,0,0},
00901                 {1.0951,0.5891,0.9222,1.2903,1.287,1.0563,1.0444,1.1405,1.4841,1.2583,2.7347,0.9844,1.034,1.0412,1.1062,1.2211,1.6728,1.6347,1.6554,1.8095,1.9561,0,0,0,0,0,0,0,0,0},
00902                 {0.9789,0.7389,1.3107,0.9274,0.9921,0.9812,0.8971,0.9936,1.0079,1.0377,1.2073,2.7198,0.9995,1.037,1.0433,1.1093,1.2323,1.6673,1.6231,1.6666,1.8089,1.7302,0,0,0,0,0,0,0,0},
00903                 {0.6169,0.4618,1.3566,3.3812,1.1951,1.1746,1.0133,0.9195,1.0548,1.0608,1.1016,1.285,3.1081,0.9986,1.0094,1.0379,1.005,0.9932,1.0651,0.8787,0.952,0.9862,1.0083,0,0,0,0,0,0,0},
00904                 {1.0603,0.262,0.9237,2.4323,1.8345,1.2197,1.0924,1.0205,0.9379,1.0946,1.1001,1.1741,1.3608,3.152,0.9692,0.8866,1.0556,1.1127,1.2289,1.6828,1.6298,1.6469,1.8082,1.0605,0,0,0,0,0,0},
00905                 {0.9061,0.7287,0.9676,1.9425,1.5785,1.9028,1.3827,1.1136,1.0368,0.9516,1.1389,1.1594,1.1642,1.4161,2.8162,0.9744,0.9246,1.0644,1.1118,1.2332,1.6892,1.6311,1.6347,1.8073,1.0722,0,0,0,0,0},
00906                 {0.9904,1.0568,1.824,1.6578,1.5033,0.9704,0.8933,0.8579,1.084,1.0308,0.9637,1.1885,1.2179,1.2979,1.4846,3.0344,0.9879,0.8927,1.0534,1.1233,1.2242,1.6794,1.6255,1.6487,1.8141,1.7629,0,0,0,0},
00907                 {0.9667,1.2622,0.959,2.8444,1.304,1.5632,1.709,2.298,1.427,1.0885,1.0285,0.9767,1.2237,1.2632,1.3585,1.5495,3.1385,0.9959,0.9327,1.0904,1.0357,1.0125,1.0027,0.9788,0.8975,1.0539,1.048,0,0,0},
00908                 {1.0004,0.8721,1.3073,0.9564,2.7856,1.6197,1.5516,2.229,2.1494,1.4916,1.0871,1.0359,0.9768,1.2269,1.3265,1.4169,1.597,3.4077,0.9979,0.8869,1.0635,1.1063,1.2267,1.6914,1.6401,1.6216,1.0598,1.0363,0,0},
00909                 {0.728,1.3939,0.4822,0.626,0.5463,2.2491,1.064,1.1998,1.3083,1.6545,1.1149,0.8826,0.8566,0.8196,1.1173,1.17,1.2445,1.394,2.9832,0.8715,0.7703,0.929,0.968,1.0828,1.4973,1.4513,1.4476,0.9621,0.9353,0},
00910                         {1.0766,0.6459,0.7688,0.2358,0.3709,0.3476,1.6754,0.8288,0.9184,1.0686,1.4546,0.9515,0.7272,0.7048,0.6911,0.9722,1.0308,1.0988,1.1959,2.6404,0.7644,0.6768,0.8181,0.8566,0.9689,1.3344,1.3031,1.3345,0.8676,0.8478}};
00911         /* Get the hybrid coeffs from the Dima rate coefficient times the average Dere to Dima ratio.
00912          */
00913         ASSERT( nelem>=0 && nelem<LIMELM && ion>=0 && ion<=nelem );
00914         double rate = coll_ion(nelem+1,nelem+1-ion,t) * DereRatio[nelem][ion];
00915         ASSERT( rate >=0. );
00916         return(rate);
00917 }
00918 realnum t_ADfA::h_coll_str( long ipLo, long ipHi, long ipTe )
00919 {
00920         realnum rate;
00921 
00922         DEBUG_ENTRY( "h_coll_str()" );
00923 
00924         ASSERT( ipLo < ipHi );
00925 
00926 #       if !defined NDEBUG
00927         long ipISO = ipH_LIKE;
00928         long nelem = ipHYDROGEN;
00929 #       endif
00930         ASSERT( N_(ipLo) < N_(ipHi) );
00931         ASSERT( N_(ipHi) <= 5 );
00932 
00933         rate = (realnum)HCS[ipHi-1][ipLo][ipTe];
00934 
00935         return rate;
00936 }
/home66/gary/public_html/cloudy/c13_branch/source/atmdat_adfa.cpp
