cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 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 00010 t_ADfA::t_ADfA() 00011 { 00012 DEBUG_ENTRY( "t_ADfA()" ); 00013 00014 /* this option, use the new atmdat_rad_rec recombination rates */ 00015 version = PHFIT_UNDEF; 00016 00017 double help[9]; 00018 const long VERSION_MAGIC = 20061204L; 00019 00020 const static char chFile[] = "phfit.dat"; 00021 00022 FILE *io = open_data( chFile, "r" ); 00023 00024 bool lgErr = false; 00025 long i=-1, j=-1, k=-1, n; 00026 00027 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 ); 00028 if( lgErr || i != VERSION_MAGIC ) 00029 { 00030 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile, i ); 00031 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC ); 00032 cdEXIT(EXIT_FAILURE); 00033 } 00034 00035 for( n=0; n < 7; n++ ) 00036 lgErr = lgErr || ( fscanf( io, "%ld", &L[n] ) != 1 ); 00037 for( n=0; n < 30; n++ ) 00038 lgErr = lgErr || ( fscanf( io, "%ld", &NINN[n] ) != 1 ); 00039 for( n=0; n < 30; n++ ) 00040 lgErr = lgErr || ( fscanf( io, "%ld", &NTOT[n] ) != 1 ); 00041 while( true ) 00042 { 00043 lgErr = lgErr || ( fscanf( io, "%ld %ld %ld", &i, &j, &k ) != 3 ); 00044 if( i == -1 && j == -1 && k == -1 ) 00045 break; 00046 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le", &help[0], &help[1], 00047 &help[2], &help[3], &help[4], &help[5] ) != 6 ); 00048 for( int l=0; l < 6; ++l ) 00049 PH1[i][j][k][l] = (realnum)help[l]; 00050 } 00051 while( true ) 00052 { 00053 lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 ); 00054 if( i == -1 && j == -1 ) 00055 break; 00056 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le", &help[0], &help[1], 00057 &help[2], &help[3], &help[4], &help[5], &help[6] ) != 7 ); 00058 for( int l=0; l < 7; ++l ) 00059 PH2[i][j][l] = (realnum)help[l]; 00060 } 00061 fclose( io ); 00062 00063 ASSERT( !lgErr ); 00064 00065 const static char chFile2[] = "hpfit.dat"; 00066 00067 io = open_data( chFile2, "r" ); 00068 00069 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 ); 00070 if( lgErr || i != VERSION_MAGIC ) 00071 { 00072 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile2, i ); 00073 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC ); 00074 cdEXIT(EXIT_FAILURE); 00075 } 00076 00077 for( i=0; i < NHYDRO_MAX_LEVEL; i++ ) 00078 { 00079 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le", &help[0], &help[1], 00080 &help[2], &help[3], &help[4] ) != 5 ); 00081 for( int l=0; l < 5; ++l ) 00082 PHH[i][l] = (realnum)help[l]; 00083 } 00084 00085 fclose( io ); 00086 00087 ASSERT( !lgErr ); 00088 00089 const static char chFile3[] = "rec_lines.dat"; 00090 00091 io = open_data( chFile3, "r" ); 00092 00093 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 ); 00094 if( lgErr || i != VERSION_MAGIC ) 00095 { 00096 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile3, i ); 00097 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC ); 00098 cdEXIT(EXIT_FAILURE); 00099 } 00100 00101 for( i=0; i < 110; i++ ) 00102 { 00103 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2], 00104 &help[3], &help[4], &help[5], &help[6], &help[7] ) != 8 ); 00105 for( int l=0; l < 8; ++l ) 00106 P[l][i] = (realnum)help[l]; 00107 } 00108 00109 00110 for( i=0; i < 405; i++ ) 00111 { 00112 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2], 00113 &help[3], &help[4], &help[5], &help[6], &help[7], &help[8] ) != 9 ); 00114 for( int l=0; l < 9; ++l ) 00115 ST[l][i] = (realnum)help[l]; 00116 } 00117 00118 fclose( io ); 00119 00120 ASSERT( !lgErr ); 00121 00122 const static char chFile4[] = "rad_rec.dat"; 00123 00124 io = open_data( chFile4, "r" ); 00125 00126 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 ); 00127 if( lgErr || i != VERSION_MAGIC ) 00128 { 00129 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile4, i ); 00130 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC ); 00131 cdEXIT(EXIT_FAILURE); 00132 } 00133 00134 while( true ) 00135 { 00136 lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 ); 00137 if( i == -1 && j == -1 ) 00138 break; 00139 lgErr = lgErr || ( fscanf( io, "%le %le", &help[0], &help[1] ) != 2 ); 00140 for( int l=0; l < 2; ++l ) 00141 rrec[i][j][l] = (realnum)help[l]; 00142 } 00143 while( true ) 00144 { 00145 lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 ); 00146 if( i == -1 && j == -1 ) 00147 break; 00148 lgErr = lgErr || ( fscanf( io, "%le %le %le %le", &help[0], &help[1], 00149 &help[2], &help[3] ) != 4 ); 00150 for( int l=0; l < 4; ++l ) 00151 rnew[i][j][l] = (realnum)help[l]; 00152 } 00153 while( true ) 00154 { 00155 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 ); 00156 if( i == -1 ) 00157 break; 00158 lgErr = lgErr || ( fscanf( io, "%le %le %le", &help[0], &help[1], &help[2] ) != 3 ); 00159 for( int l=0; l < 3; ++l ) 00160 fe[i][l] = (realnum)help[l]; 00161 } 00162 00163 fclose( io ); 00164 00165 ASSERT( !lgErr ); 00166 00167 const static char chFile5[] = "h_rad_rec.dat"; 00168 00169 io = open_data( chFile5, "r" ); 00170 00171 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 ); 00172 if( lgErr || i != VERSION_MAGIC ) 00173 { 00174 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile5, i ); 00175 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC ); 00176 cdEXIT(EXIT_FAILURE); 00177 } 00178 00179 for( i=0; i < NHYDRO_MAX_LEVEL; i++ ) 00180 { 00181 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2], 00182 &help[3], &help[4], &help[5], &help[6], &help[7], &help[8] ) != 9 ); 00183 for( int l=0; l < 9; ++l ) 00184 HRF[i][l] = (realnum)help[l]; 00185 } 00186 00187 fclose( io ); 00188 00189 ASSERT( !lgErr ); 00190 00191 const static char chFile6[] = "h_phot_cs.dat"; 00192 00193 io = open_data( chFile6, "r" ); 00194 00195 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 ); 00196 if( lgErr || i != VERSION_MAGIC ) 00197 { 00198 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile6, i ); 00199 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC ); 00200 cdEXIT(EXIT_FAILURE); 00201 } 00202 00203 for( i=0; i < NHYDRO_MAX_LEVEL; i++ ) 00204 { 00205 lgErr = lgErr || ( fscanf( io, "%le", &help[0] ) != 1 ); 00206 STH[i] = (realnum)help[0]; 00207 } 00208 00209 fclose( io ); 00210 00211 ASSERT( !lgErr ); 00212 00213 const static char chFile7[] = "coll_ion.dat"; 00214 00215 io = open_data( chFile7, "r" ); 00216 00217 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 ); 00218 if( lgErr || i != VERSION_MAGIC ) 00219 { 00220 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile7, i ); 00221 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC ); 00222 cdEXIT(EXIT_FAILURE); 00223 } 00224 00225 while( true ) 00226 { 00227 lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 ); 00228 if( i == -1 && j == -1 ) 00229 break; 00230 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le", &CF[i][j][0], &CF[i][j][1], 00231 &CF[i][j][2], &CF[i][j][3], &CF[i][j][4] ) != 5 ); 00232 } 00233 00234 fclose( io ); 00235 00236 ASSERT( !lgErr ); 00237 00238 const static char chFile8[] = "h_coll_str.dat"; 00239 00240 io = open_data( chFile8, "r" ); 00241 00242 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 ); 00243 if( lgErr || i != VERSION_MAGIC ) 00244 { 00245 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile8, i ); 00246 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC ); 00247 cdEXIT(EXIT_FAILURE); 00248 } 00249 00250 while( true ) 00251 { 00252 lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 ); 00253 if( i == -1 && j == -1 ) 00254 break; 00255 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le", &HCS[i-2][j-1][0], &HCS[i-2][j-1][1], 00256 &HCS[i-2][j-1][2], &HCS[i-2][j-1][3], &HCS[i-2][j-1][4], &HCS[i-2][j-1][5], 00257 &HCS[i-2][j-1][6], &HCS[i-2][j-1][7] ) != 8 ); 00258 } 00259 00260 fclose( io ); 00261 00262 ASSERT( !lgErr ); 00263 } 00264 00265 double t_ADfA::phfit(long int nz, 00266 long int ne, 00267 long int is, 00268 double e) 00269 { 00270 long int nint, 00271 nout; 00272 double a, 00273 b, 00274 crs, 00275 einn, 00276 p1, 00277 q, 00278 x, 00279 y, 00280 z; 00281 00282 DEBUG_ENTRY( "phfit()" ); 00283 00284 /*** Version 3. October 8, 1996. 00285 *** Written by D. A. Verner, verner@pa.uky.edu 00286 *** Inner-shell ionization energies of some low-ionized species are slightly 00287 *** improved to fit smoothly the experimental inner-shell ionization energies 00288 *** of neutral atoms. 00289 ****************************************************************************** 00290 *** This subroutine calculates partial photoionization cross sections 00291 *** for all ionization stages of all atoms from H to Zn (Z=30) by use of 00292 *** the following fit parameters: 00293 *** Outer shells of the Opacity Project (OP) elements: 00294 *** >>refer all photocs Verner, Ferland, Korista, Yakovlev, 1996, ApJ, in press. 00295 *** Inner shells of all elements, and outer shells of the non-OP elements: 00296 *** Verner and Yakovlev, 1995, A&AS, 109, 125 00297 *** Input parameters: nz - atomic number from 1 to 30 (integer) 00298 *** ne - number of electrons from 1 to iz (integer) 00299 *** is - shell number (integer) 00300 *** e - photon energy, eV 00301 *** version - enum, PHFIT96 (default): calculates 00302 *** new cross sections, PHFIT95: calculates 00303 *** only old Hartree-Slater cross sections 00304 *** Output parameter: photoionization cross section, Mb 00305 *** Shell numbers: 00306 *** 1 - 1s, 2 - 2s, 3 - 2p, 4 - 3s, 5 - 3p, 6 - 3d, 7 - 4s. 00307 *** If a species in the ground state has no electrons on the given shell, 00308 *** the subroutine returns 0. 00309 ****************************************************************************** */ 00310 00311 crs = 0.0; 00312 if( nz < 1 || nz > 30 ) 00313 { 00314 return(crs); 00315 } 00316 00317 if( ne < 1 || ne > nz ) 00318 { 00319 return(crs); 00320 } 00321 00322 nout = NTOT[ne-1]; 00323 if( nz == ne && nz > 18 ) 00324 nout = 7; 00325 if( nz == (ne + 1) && ((((nz == 20 || nz == 21) || nz == 22) || 00326 nz == 25) || nz == 26) ) 00327 nout = 7; 00328 if( is > nout ) 00329 { 00330 return(crs); 00331 } 00332 00333 if( (is == 6 && (nz == 20 || nz == 19)) && ne >= 19 ) 00334 { 00335 return(crs); 00336 } 00337 00338 ASSERT( is >= 1 && is <= 7 ); 00339 00340 if( e < PH1[is-1][ne-1][nz-1][0] ) 00341 { 00342 return(crs); 00343 } 00344 00345 nint = NINN[ne-1]; 00346 if( ((nz == 15 || nz == 17) || nz == 19) || (nz > 20 && nz != 26) ) 00347 { 00348 einn = 0.0; 00349 } 00350 else 00351 { 00352 if( ne < 3 ) 00353 { 00354 einn = 1.0e30; 00355 } 00356 else 00357 { 00358 einn = PH1[nint-1][ne-1][nz-1][0]; 00359 } 00360 } 00361 00362 if( (is <= nint || e >= einn) || version == PHFIT95 ) 00363 { 00364 p1 = -PH1[is-1][ne-1][nz-1][4]; 00365 y = e/PH1[is-1][ne-1][nz-1][1]; 00366 q = -0.5*p1 - L[is-1] - 5.5; 00367 a = PH1[is-1][ne-1][nz-1][2]*(POW2(y - 1.0) + 00368 POW2(PH1[is-1][ne-1][nz-1][5])); 00369 b = sqrt(y/PH1[is-1][ne-1][nz-1][3]) + 1.0; 00370 crs = a*pow(y,q)*pow(b,p1); 00371 } 00372 else 00373 { 00374 if( (is < nout && is > nint) && e < einn ) 00375 { 00376 return(crs); 00377 } 00378 p1 = -PH2[ne-1][nz-1][3]; 00379 q = -0.5*p1 - 5.5; 00380 x = e/PH2[ne-1][nz-1][0] - PH2[ne-1][nz-1][5]; 00381 z = sqrt(x*x+POW2(PH2[ne-1][nz-1][6])); 00382 a = PH2[ne-1][nz-1][1]*(POW2(x - 1.0) + 00383 POW2(PH2[ne-1][nz-1][4])); 00384 b = 1.0 + sqrt(z/PH2[ne-1][nz-1][2]); 00385 crs = a*pow(z,q)*pow(b,p1); 00386 } 00387 return(crs); 00388 } 00389 00390 /* same as hpfit, but energy is relative to threshold */ 00391 double t_ADfA::hpfit_rel(long int iz, 00392 long int n, 00393 double e) 00394 { 00395 long m; 00396 double crs , ex , eth; 00397 00398 if( n == 0 ) 00399 { 00400 m = 1; 00401 } 00402 else 00403 { 00404 if( n == 1 ) 00405 { 00406 m = 2; 00407 } 00408 else 00409 { 00410 m = n; 00411 } 00412 } 00413 00414 eth = ph1(0,0,iz-1,0)/POW2((double)m); 00415 ex = MAX2(1. , e/eth ); 00416 00417 crs = hpfit( iz , n , ex ); 00418 ASSERT( crs > 0. ); 00419 00420 return crs; 00421 } 00422 00423 double t_ADfA::hpfit(long int iz, 00424 long int n, 00425 double e) 00426 { 00427 long int l, 00428 m; 00429 double cs, 00430 eth, 00431 ex, 00432 q, 00433 x; 00434 00435 DEBUG_ENTRY( "hpfit()" ); 00436 00437 /*state specific photoionization cross sections for model hydrogen atom 00438 * Version 1, September 23, 1997 00439 ****************************************************************************** 00440 *** This subroutine calculates state-specific photoionization cross sections 00441 *** for hydrogen and hydrogen-like ions. 00442 *** Input parameters: iz - atomic number 00443 *** n - shell number, from 0 to 400: 00444 *** 0 - 1s 00445 *** 1 - 2s 00446 *** 2 - 2p 00447 *** 3 - 3 00448 *** ...... 00449 *** e - photon energy, eV 00450 *** return value - cross section, cm^(-2) 00451 *******************************************************************************/ 00452 00453 cs = 0.0; 00454 00455 ASSERT( iz > 0 && e>0. ); 00456 00457 if( n >= NHYDRO_MAX_LEVEL ) 00458 { 00459 fprintf( ioQQQ, " hpfit called with too large n, =%li\n" , n ); 00460 cdEXIT(EXIT_FAILURE); 00461 } 00462 00463 l = 0; 00464 if( n == 2 ) 00465 { 00466 l = 1; 00467 } 00468 q = 3.5 + l - 0.5*PHH[n][1]; 00469 00470 if( n == 0 ) 00471 { 00472 m = 1; 00473 } 00474 else 00475 { 00476 if( n == 1 ) 00477 { 00478 m = 2; 00479 } 00480 else 00481 { 00482 m = n; 00483 } 00484 } 00485 00486 eth = ph1(0,0,iz-1,0)/POW2((double)m); 00487 ex = MAX2(1. , e/eth ); 00488 00489 /* Don't just force to be at least one...make sure e/eth is close to one or greater. */ 00490 ASSERT( e/eth > 0.95 ); 00491 00492 if( ex < 1.0 ) 00493 { 00494 return(0.); 00495 } 00496 00497 x = ex/PHH[n][0]; 00498 cs = (PHH[n][4]*pow(1.0 + ((double)PHH[n][2]/x),(double)PHH[n][3])/ 00499 pow(x,q)/pow(1.0 + sqrt(x),(double)PHH[n][1])*8.79737e-17/ 00500 POW2((double)iz)); 00501 return(cs); 00502 } 00503 00504 void t_ADfA::rec_lines(double t, 00505 realnum r[][471]) 00506 { 00507 long int i, 00508 j, 00509 ipj1, 00510 ipj2; 00511 00512 double a, 00513 a1, 00514 dr[4][405], 00515 p1, 00516 p2, 00517 p3, 00518 p4, 00519 p5, 00520 p6, 00521 rr[4][110], 00522 te, 00523 x, 00524 z; 00525 00526 static long jd[6]={143,145,157,360,376,379}; 00527 00528 static long ip[38]={7,9,12,13,14,16,18,19,20,21,22,44,45,49,50, 00529 52,53,54,55,56,57,58,59,60,66,67,78,83,84,87,88,95,96,97,100, 00530 101,103,104}; 00531 00532 static long id[38]={7,3,10,27,23,49,51,64,38,47,60,103,101,112, 00533 120,114,143,145,157,152,169,183,200,163,225,223,237,232,235, 00534 249,247,300,276,278,376,360,379,384}; 00535 00536 DEBUG_ENTRY( "rec_lines()" ); 00537 00538 /*effective recombination coefficients for lines of C, N, O, by D. Verner 00539 * Version 2, April 30, 1997 00540 ****************************************************************************** 00541 *** This subroutine calculates effective recombination coefficients 00542 *** for 110 permitted recombination lines of C, N, O (Pequignot, Petitjean, 00543 *** & Boisson, 1991, A&A, 251, 680) and 405 permitted dielectronic 00544 *** recombination lines (Nussbaumer & Storey, 1984, A&AS, 56, 293) 00545 *** Input parameter: t - temperature, K 00546 *** Output parameters: r(i,j), i=1,471 00547 *** r(i,1) - atomic number 00548 *** r(i,2) - number of electrons 00549 *** r(i,3) - wavelength, angstrom 00550 *** r(i,4) - rate coefficient, cm^3 s^(-1) 00551 ****************************************************************************** */ 00552 00553 for( i=0; i < 110; i++ ) 00554 { 00555 rr[0][i] = P[0][i]; 00556 rr[1][i] = P[1][i]; 00557 rr[2][i] = P[2][i]; 00558 z = P[0][i] - P[1][i] + 1.0; 00559 te = 1.0e-04*t/z/z; 00560 p1 = P[3][i]; 00561 p2 = P[4][i]; 00562 p3 = P[5][i]; 00563 p4 = P[6][i]; 00564 if( te < 0.004 ) 00565 { 00566 a1 = p1*pow(0.004,p2)/(1.0 + p3*pow(0.004,p4)); 00567 a = a1/sqrt(te/0.004); 00568 } 00569 else 00570 { 00571 if( te > 2.0 ) 00572 { 00573 a1 = p1*pow(2.0,p2)/(1.0 + p3*pow(2.0,p4)); 00574 a = a1/pow(te/2.0,1.5); 00575 } 00576 else 00577 { 00578 a = p1*pow(te,p2)/(1.0 + p3*pow(te,p4)); 00579 } 00580 } 00581 rr[3][i] = 1.0e-13*z*a*P[7][i]; 00582 } 00583 00584 for( i=0; i < 405; i++ ) 00585 { 00586 dr[0][i] = ST[0][i]; 00587 dr[1][i] = ST[1][i]; 00588 dr[2][i] = ST[2][i]; 00589 te = 1.0e-04*t; 00590 p1 = ST[3][i]; 00591 p2 = ST[4][i]; 00592 p3 = ST[5][i]; 00593 p4 = ST[6][i]; 00594 p5 = ST[7][i]; 00595 p6 = ST[8][i]; 00596 if( te < p6 ) 00597 { 00598 x = p5*(1.0/te - 1.0/p6); 00599 if( x > 80.0 ) 00600 { 00601 a = 0.0; 00602 } 00603 else 00604 { 00605 a1 = (p1/p6 + p2 + p3*p6 + p4*p6*p6)/pow(p6,1.5)/exp(p5/ 00606 p6); 00607 a = a1/exp(x); 00608 } 00609 } 00610 else 00611 { 00612 if( te > 6.0 ) 00613 { 00614 a1 = (p1/6.0 + p2 + p3*6.0 + p4*36.0)/pow(6.0,1.5)/ 00615 exp(p5/6.0); 00616 a = a1/pow(te/6.0,1.5); 00617 } 00618 else 00619 { 00620 a = (p1/te + p2 + p3*te + p4*te*te)/pow(te,1.5)/exp(p5/ 00621 te); 00622 } 00623 } 00624 dr[3][i] = 1.0e-12*a; 00625 } 00626 00627 for( i=0; i < 6; i++ ) 00628 { 00629 ipj1 = jd[i]; 00630 ipj2 = ipj1 + 1; 00631 dr[3][ipj1-1] += dr[3][ipj2-1]; 00632 dr[0][ipj2-1] = 0.0; 00633 } 00634 00635 for( i=0; i < 38; i++ ) 00636 { 00637 ipj1 = ip[i]; 00638 ipj2 = id[i]; 00639 rr[3][ipj1-1] += dr[3][ipj2-1]; 00640 dr[0][ipj2-1] = 0.0; 00641 } 00642 00643 for( i=0; i < 110; i++ ) 00644 { 00645 r[0][i] = (realnum)rr[0][i]; 00646 r[1][i] = (realnum)rr[1][i]; 00647 r[2][i] = (realnum)rr[2][i]; 00648 r[3][i] = (realnum)rr[3][i]; 00649 } 00650 00651 j = 110; 00652 for( i=0; i < 405; i++ ) 00653 { 00654 if( dr[0][i] > 1.0 ) 00655 { 00656 j += 1; 00657 r[0][j-1] = (realnum)dr[0][i]; 00658 r[1][j-1] = (realnum)dr[1][i]; 00659 r[2][j-1] = (realnum)dr[2][i]; 00660 r[3][j-1] = (realnum)dr[3][i]; 00661 } 00662 } 00663 return; 00664 } 00665 00666 double t_ADfA::rad_rec(long int iz, 00667 long int in, 00668 double t) 00669 { 00670 /* 00671 *** Version 4. June 29, 1999. 00672 *** Written by D. A. Verner, verner@pa.uky.edu 00673 ****************************************************************************** 00674 *** This subroutine calculates rates of radiative recombination for all ions 00675 *** of all elements from H through Zn by use of the following fits: 00676 *** H-like, He-like, Li-like, Na-like - 00677 *** >>refer all reccoef Verner & Ferland, 1996, ApJS, 103, 467 00678 *** Other ions of C, N, O, Ne - Pequignot et al. 1991, A&A, 251, 680, 00679 *** refitted by Verner & Ferland formula to ensure correct asymptotes 00680 *** Fe XVII-XXIII - 00681 *** >>refer Fe17-23 recom Arnaud & Raymond, 1992, ApJ, 398, 394 00682 *** Fe I-XV - refitted by Verner & Ferland formula to ensure correct asymptotes 00683 *** Other ions of Mg, Si, S, Ar, Ca, Fe, Ni - 00684 *** - 00685 *** >>refer all recom Shull & Van Steenberg, 1982, ApJS, 48, 95 00686 *** Other ions of Na, Al - 00687 *** >>refer Na, Al recom Landini & Monsignori Fossi, 1990, A&AS, 82, 229 00688 *** Other ions of F, P, Cl, K, Ti, Cr, Mn, Co (excluding Ti I-II, Cr I-IV, 00689 *** Mn I-V, Co I) - 00690 *** >>refer many recom Landini & Monsignori Fossi, 1991, A&AS, 91, 183 00691 *** All other species - interpolations of the power-law fits 00692 *** Input parameters: iz - atomic number 00693 *** in - number of electrons from 1 to iz 00694 *** t - temperature, K 00695 *** return result: - rate coefficient, cm^3 s^(-1) 00696 ****************************************************************************** 00697 */ 00698 double tt; 00699 double rate; 00700 00701 DEBUG_ENTRY( "rad_rec()" ); 00702 00703 rate = 0.0; 00704 if( iz < 1 || iz > 30 ) 00705 { 00706 fprintf( ioQQQ, " rad_rec called with insane atomic number, =%4ld\n", 00707 iz ); 00708 cdEXIT(EXIT_FAILURE); 00709 } 00710 if( in < 1 || in > iz ) 00711 { 00712 fprintf( ioQQQ, " rad_rec called with insane number elec =%4ld\n", 00713 in ); 00714 cdEXIT(EXIT_FAILURE); 00715 } 00716 if( (((in <= 3 || in == 11) || (iz > 5 && iz < 9)) || iz == 10) || 00717 (iz == 26 && in > 11) ) 00718 { 00719 tt = sqrt(t/rnew[in-1][iz-1][2]); 00720 rate = 00721 rnew[in-1][iz-1][0]/(tt*pow(tt + 1.0,1.0 - rnew[in-1][iz-1][1])* 00722 pow(1.0 + sqrt(t/rnew[in-1][iz-1][3]),1.0 + rnew[in-1][iz-1][1])); 00723 } 00724 else 00725 { 00726 tt = t*1.0e-04; 00727 if( iz == 26 && in <= 13 ) 00728 { 00729 rate = fe[in-1][0]/pow(tt,fe[in-1][1] + 00730 fe[in-1][2]*log10(tt)); 00731 } 00732 else 00733 { 00734 rate = rrec[in-1][iz-1][0]/pow(tt,(double)rrec[in-1][iz-1][1]); 00735 } 00736 } 00737 00738 return rate; 00739 } 00740 00741 double t_ADfA::H_rad_rec(long int iz, 00742 long int n, 00743 double t) 00744 { 00745 /* 00746 * Version 4, October 9, 1997 00747 ****************************************************************************** 00748 *** This subroutine calculates state-specific recombination rates 00749 *** for hydrogen and hydrogen-like ions. 00750 *** Input parameters: iz - atomic number 00751 *** n - shell number, from 0 to 400: 00752 *** 0 - 1s 00753 *** 1 - 2s 00754 *** 2 - 2p 00755 *** 3 - 3 00756 *** ...... 00757 *** t - temperature, K 00758 *** Output parameter: r - rate coefficient, cm^3 s^(-1) 00759 *** If n is negative, the subroutine returns the total recombination 00760 *** rate coefficient 00761 ****************************************************************************** 00762 */ 00763 double rate, 00764 TeScaled, 00765 x, 00766 x1, 00767 x2; 00768 00769 DEBUG_ENTRY( "H_rad_rec()" ); 00770 00771 rate = 0.0; 00772 00773 /* iz is charge, must be 1 or greater */ 00774 ASSERT( iz > 0 ); 00775 00776 /* n is level number, must be less than dim or hydro vectors */ 00777 ASSERT( n < NHYDRO_MAX_LEVEL ); 00778 00779 TeScaled = t/POW2((double)iz); 00780 00781 if( n < 0 ) 00782 { 00783 x1 = sqrt(TeScaled/3.148); 00784 x2 = sqrt(TeScaled/7.036e05); 00785 rate = 7.982e-11/x1/pow(1.0 + x1,0.252)/pow(1.0 + x2,1.748); 00786 } 00787 else 00788 { 00789 x = log10(TeScaled); 00790 rate = (HRF[n][0] + HRF[n][2]*x + HRF[n][4]* 00791 x*x + HRF[n][6]*powi(x,3) + HRF[n][8]*powi(x,4))/ 00792 (1.0 + HRF[n][1]*x + HRF[n][3]*x*x + HRF[n][5]* 00793 powi(x,3) + HRF[n][7]*powi(x,4)); 00794 rate = pow(10.0,rate)/TeScaled; 00795 } 00796 rate *= iz; 00797 00798 return rate; 00799 } 00800 00801 /*coll_ion D Verner's routine to compute collisional ionization rate coefficients, 00802 * returns collisional ionization rate coefficient cm^3 s^-1*/ 00803 double t_ADfA::coll_ion( 00804 /* atomic number, 1 for hydrogen */ 00805 long int iz, 00806 /* stage of ionization, 1 for atom */ 00807 long int in, 00808 /* temperature */ 00809 double t) 00810 { 00811 double rate, te, u; 00812 00813 DEBUG_ENTRY( "coll_ion()" ); 00814 /*D Verner's routine to compute collisional ionization rate coefficients 00815 * Version 3, April 21, 1997 00816 * Cu (Z=29) and Zn (Z=30) are added (fits from Ni, correct thresholds). 00817 ****************************************************************************** 00818 *** This subroutine calculates rates of direct collisional ionization 00819 *** for all ionization stages of all elements from H to Ni (Z=28) 00820 *** by use of the fits from 00821 *>>refer all collion Voronov, G. S., 1997, At. Data Nucl. Data Tables, 65, 1 00822 *** Input parameters: iz - atomic number on pphysical scale, H is 1 00823 *** in - number of electrons from 1 to iz 00824 *** t - temperature, K 00825 *** Output parameter: c - rate coefficient, cm^3 s^(-1) 00826 ****************************************************************************** */ 00827 rate = 0.0; 00828 00829 if( iz < 1 || iz > 30 ) 00830 { 00831 /* return zero rate is atomic number outside range of code */ 00832 return( 0. ); 00833 } 00834 00835 if( in < 1 || in > iz ) 00836 { 00837 /* return zero rate is ion stage is impossible */ 00838 return( 0. ); 00839 } 00840 00841 te = t*EVRYD/TE1RYD; 00842 u = CF[in-1][iz-1][0]/te; 00843 if( u > 80.0 ) 00844 { 00845 return( 0. ); 00846 } 00847 00848 rate = (CF[in-1][iz-1][2]*(1.0 + CF[in-1][iz-1][1]* 00849 sqrt(u))/(CF[in-1][iz-1][3] + u)*pow(u,CF[in-1][iz-1][4])* 00850 exp(-u)); 00851 00852 return(rate); 00853 } 00854 00855 realnum t_ADfA::h_coll_str( long ipLo, long ipHi, long ipTe ) 00856 { 00857 realnum rate; 00858 00859 DEBUG_ENTRY( "h_coll_str()" ); 00860 00861 ASSERT( ipLo < ipHi ); 00862 00863 #ifndef NDEBUG 00864 long ipISO = ipH_LIKE; 00865 long nelem = ipHYDROGEN; 00866 #endif 00867 ASSERT( N_(ipLo) < N_(ipHi) ); 00868 ASSERT( N_(ipHi) <= 5 ); 00869 00870 rate = (realnum)HCS[ipHi-1][ipLo][ipTe]; 00871 00872 return rate; 00873 }