00001
00002
00003
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "atmdat.h"
00007 #include "iso.h"
00008 #include "taulines.h"
00009
00011 t_ADfA::t_ADfA()
00012 {
00013 DEBUG_ENTRY( "t_ADfA()" );
00014
00015
00016 version = PHFIT_UNDEF;
00017
00018 double help[9];
00019 const long VERSION_MAGIC = 20061204L;
00020
00021 static const char chFile[] = "phfit.dat";
00022
00023 FILE *io = open_data( chFile, "r" );
00024
00025 bool lgErr = false;
00026 long i=-1, j=-1, k=-1, n;
00027
00028 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00029 if( lgErr || i != VERSION_MAGIC )
00030 {
00031 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile, i );
00032 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00033 cdEXIT(EXIT_FAILURE);
00034 }
00035
00036 for( n=0; n < 7; n++ )
00037 lgErr = lgErr || ( fscanf( io, "%ld", &L[n] ) != 1 );
00038 for( n=0; n < 30; n++ )
00039 lgErr = lgErr || ( fscanf( io, "%ld", &NINN[n] ) != 1 );
00040 for( n=0; n < 30; n++ )
00041 lgErr = lgErr || ( fscanf( io, "%ld", &NTOT[n] ) != 1 );
00042 while( true )
00043 {
00044 lgErr = lgErr || ( fscanf( io, "%ld %ld %ld", &i, &j, &k ) != 3 );
00045 if( i == -1 && j == -1 && k == -1 )
00046 break;
00047 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le", &help[0], &help[1],
00048 &help[2], &help[3], &help[4], &help[5] ) != 6 );
00049 for( int l=0; l < 6; ++l )
00050 PH1[i][j][k][l] = (realnum)help[l];
00051 }
00052 while( true )
00053 {
00054 lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
00055 if( i == -1 && j == -1 )
00056 break;
00057 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le", &help[0], &help[1],
00058 &help[2], &help[3], &help[4], &help[5], &help[6] ) != 7 );
00059 for( int l=0; l < 7; ++l )
00060 PH2[i][j][l] = (realnum)help[l];
00061 }
00062 fclose( io );
00063
00064 ASSERT( !lgErr );
00065
00066 static const char chFile2[] = "hpfit.dat";
00067
00068 io = open_data( chFile2, "r" );
00069
00070 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00071 if( lgErr || i != VERSION_MAGIC )
00072 {
00073 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile2, i );
00074 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00075 cdEXIT(EXIT_FAILURE);
00076 }
00077
00078 for( i=0; i < NHYDRO_MAX_LEVEL; i++ )
00079 {
00080 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le", &help[0], &help[1],
00081 &help[2], &help[3], &help[4] ) != 5 );
00082 for( int l=0; l < 5; ++l )
00083 PHH[i][l] = (realnum)help[l];
00084 }
00085
00086 fclose( io );
00087
00088 ASSERT( !lgErr );
00089
00090 static const char chFile3[] = "rec_lines.dat";
00091
00092 io = open_data( chFile3, "r" );
00093
00094 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00095 if( lgErr || i != VERSION_MAGIC )
00096 {
00097 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile3, i );
00098 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00099 cdEXIT(EXIT_FAILURE);
00100 }
00101
00102 for( i=0; i < 110; i++ )
00103 {
00104 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2],
00105 &help[3], &help[4], &help[5], &help[6], &help[7] ) != 8 );
00106 for( int l=0; l < 8; ++l )
00107 P[l][i] = (realnum)help[l];
00108 }
00109
00110
00111 for( i=0; i < 405; i++ )
00112 {
00113 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2],
00114 &help[3], &help[4], &help[5], &help[6], &help[7], &help[8] ) != 9 );
00115 for( int l=0; l < 9; ++l )
00116 ST[l][i] = (realnum)help[l];
00117 }
00118
00119 fclose( io );
00120
00121 ASSERT( !lgErr );
00122
00123 static const char chFile4[] = "rad_rec.dat";
00124
00125 io = open_data( chFile4, "r" );
00126
00127 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00128 if( lgErr || i != VERSION_MAGIC )
00129 {
00130 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile4, i );
00131 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00132 cdEXIT(EXIT_FAILURE);
00133 }
00134
00135 while( true )
00136 {
00137 lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
00138 if( i == -1 && j == -1 )
00139 break;
00140 lgErr = lgErr || ( fscanf( io, "%le %le", &help[0], &help[1] ) != 2 );
00141 for( int l=0; l < 2; ++l )
00142 rrec[i][j][l] = (realnum)help[l];
00143 }
00144 while( true )
00145 {
00146 lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
00147 if( i == -1 && j == -1 )
00148 break;
00149 lgErr = lgErr || ( fscanf( io, "%le %le %le %le", &help[0], &help[1],
00150 &help[2], &help[3] ) != 4 );
00151 for( int l=0; l < 4; ++l )
00152 rnew[i][j][l] = (realnum)help[l];
00153 }
00154 while( true )
00155 {
00156 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00157 if( i == -1 )
00158 break;
00159 lgErr = lgErr || ( fscanf( io, "%le %le %le", &help[0], &help[1], &help[2] ) != 3 );
00160 for( int l=0; l < 3; ++l )
00161 fe[i][l] = (realnum)help[l];
00162 }
00163
00164 fclose( io );
00165
00166 ASSERT( !lgErr );
00167
00168 static const char chFile5[] = "h_rad_rec.dat";
00169
00170 io = open_data( chFile5, "r" );
00171
00172 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00173 if( lgErr || i != VERSION_MAGIC )
00174 {
00175 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile5, i );
00176 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00177 cdEXIT(EXIT_FAILURE);
00178 }
00179
00180 for( i=0; i < NHYDRO_MAX_LEVEL; i++ )
00181 {
00182 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2],
00183 &help[3], &help[4], &help[5], &help[6], &help[7], &help[8] ) != 9 );
00184 for( int l=0; l < 9; ++l )
00185 HRF[i][l] = (realnum)help[l];
00186 }
00187
00188 fclose( io );
00189
00190 ASSERT( !lgErr );
00191
00192 static const char chFile6[] = "h_phot_cs.dat";
00193
00194 io = open_data( chFile6, "r" );
00195
00196 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00197 if( lgErr || i != VERSION_MAGIC )
00198 {
00199 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile6, i );
00200 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00201 cdEXIT(EXIT_FAILURE);
00202 }
00203
00204 for( i=0; i < NHYDRO_MAX_LEVEL; i++ )
00205 {
00206 lgErr = lgErr || ( fscanf( io, "%le", &help[0] ) != 1 );
00207 STH[i] = (realnum)help[0];
00208 }
00209
00210 fclose( io );
00211
00212 ASSERT( !lgErr );
00213
00214 static const char chFile7[] = "coll_ion.dat";
00215
00216 io = open_data( chFile7, "r" );
00217
00218 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00219 if( lgErr || i != VERSION_MAGIC )
00220 {
00221 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile7, i );
00222 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00223 cdEXIT(EXIT_FAILURE);
00224 }
00225
00226 while( true )
00227 {
00228 lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
00229 if( i == -1 && j == -1 )
00230 break;
00231 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le", &CF[i][j][0], &CF[i][j][1],
00232 &CF[i][j][2], &CF[i][j][3], &CF[i][j][4] ) != 5 );
00233 }
00234
00235 fclose( io );
00236
00237 ASSERT( !lgErr );
00238
00239
00240
00241 static const char chFile8[] = "h_coll_str.dat";
00242
00243 io = open_data( chFile8, "r" );
00244
00245 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00246 if( lgErr || i != VERSION_MAGIC )
00247 {
00248 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile8, i );
00249 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00250 cdEXIT(EXIT_FAILURE);
00251 }
00252
00253 while( true )
00254 {
00255 lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
00256 if( i == -1 && j == -1 )
00257 break;
00258 lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le", &HCS[i-2][j-1][0], &HCS[i-2][j-1][1],
00259 &HCS[i-2][j-1][2], &HCS[i-2][j-1][3], &HCS[i-2][j-1][4], &HCS[i-2][j-1][5],
00260 &HCS[i-2][j-1][6], &HCS[i-2][j-1][7] ) != 8 );
00261 }
00262
00263 fclose( io );
00264
00265 ASSERT( !lgErr );
00266 }
00267
00268 double t_ADfA::phfit(long int nz,
00269 long int ne,
00270 long int is,
00271 double e)
00272 {
00273 long int nint,
00274 nout;
00275 double a,
00276 b,
00277 crs,
00278 einn,
00279 p1,
00280 q,
00281 x,
00282 y,
00283 z;
00284
00285 DEBUG_ENTRY( "phfit()" );
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314 crs = 0.0;
00315 if( nz < 1 || nz > 30 )
00316 {
00317 return(crs);
00318 }
00319
00320 if( ne < 1 || ne > nz )
00321 {
00322 return(crs);
00323 }
00324
00325 nout = NTOT[ne-1];
00326 if( nz == ne && nz > 18 )
00327 nout = 7;
00328 if( nz == (ne + 1) && ((((nz == 20 || nz == 21) || nz == 22) ||
00329 nz == 25) || nz == 26) )
00330 nout = 7;
00331 if( is > nout )
00332 {
00333 return(crs);
00334 }
00335
00336 if( (is == 6 && (nz == 20 || nz == 19)) && ne >= 19 )
00337 {
00338 return(crs);
00339 }
00340
00341 ASSERT( is >= 1 && is <= 7 );
00342
00343 if( e < PH1[is-1][ne-1][nz-1][0] )
00344 {
00345 return(crs);
00346 }
00347
00348 nint = NINN[ne-1];
00349 if( ((nz == 15 || nz == 17) || nz == 19) || (nz > 20 && nz != 26) )
00350 {
00351 einn = 0.0;
00352 }
00353 else
00354 {
00355 if( ne < 3 )
00356 {
00357 einn = 1.0e30;
00358 }
00359 else
00360 {
00361 einn = PH1[nint-1][ne-1][nz-1][0];
00362 }
00363 }
00364
00365 if( (is <= nint || e >= einn) || version == PHFIT95 )
00366 {
00367 p1 = -PH1[is-1][ne-1][nz-1][4];
00368 y = e/PH1[is-1][ne-1][nz-1][1];
00369 q = -0.5*p1 - L[is-1] - 5.5;
00370 a = PH1[is-1][ne-1][nz-1][2]*(POW2(y - 1.0) +
00371 POW2(PH1[is-1][ne-1][nz-1][5]));
00372 b = sqrt(y/PH1[is-1][ne-1][nz-1][3]) + 1.0;
00373 crs = a*pow(y,q)*pow(b,p1);
00374 }
00375 else
00376 {
00377 if( (is < nout && is > nint) && e < einn )
00378 {
00379 return(crs);
00380 }
00381 p1 = -PH2[ne-1][nz-1][3];
00382 q = -0.5*p1 - 5.5;
00383 x = e/PH2[ne-1][nz-1][0] - PH2[ne-1][nz-1][5];
00384 z = sqrt(x*x+POW2(PH2[ne-1][nz-1][6]));
00385 a = PH2[ne-1][nz-1][1]*(POW2(x - 1.0) +
00386 POW2(PH2[ne-1][nz-1][4]));
00387 b = 1.0 + sqrt(z/PH2[ne-1][nz-1][2]);
00388 crs = a*pow(z,q)*pow(b,p1);
00389 }
00390 return(crs);
00391 }
00392
00393 double t_ADfA::hpfit(long int iz,
00394 long int n,
00395 double e)
00396 {
00397 long int l,
00398 m;
00399 double cs,
00400 eth,
00401 ex,
00402 q,
00403 x;
00404
00405 DEBUG_ENTRY( "hpfit()" );
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423 cs = 0.0;
00424
00425 ASSERT( iz > 0 && e>0. );
00426
00427 if( n >= NHYDRO_MAX_LEVEL )
00428 {
00429 fprintf( ioQQQ, " hpfit called with too large n, =%li\n" , n );
00430 cdEXIT(EXIT_FAILURE);
00431 }
00432
00433 l = 0;
00434 if( n == 2 )
00435 {
00436 l = 1;
00437 }
00438 q = 3.5 + l - 0.5*PHH[n][1];
00439
00440 if( n == 0 )
00441 {
00442 m = 1;
00443 }
00444 else
00445 {
00446 if( n == 1 )
00447 {
00448 m = 2;
00449 }
00450 else
00451 {
00452 m = n;
00453 }
00454 }
00455
00456 eth = ph1(0,0,iz-1,0)/POW2((double)m);
00457 ex = MAX2(1. , e/eth );
00458
00459
00460 ASSERT( e/eth > 0.95 );
00461
00462 if( ex < 1.0 )
00463 {
00464 return(0.);
00465 }
00466
00467 x = ex/PHH[n][0];
00468 cs = (PHH[n][4]*pow(1.0 + ((double)PHH[n][2]/x),(double)PHH[n][3])/
00469 pow(x,q)/pow(1.0 + sqrt(x),(double)PHH[n][1])*8.79737e-17/
00470 POW2((double)iz));
00471 return(cs);
00472 }
00473
00474 void t_ADfA::rec_lines(double t,
00475 realnum r[][471])
00476 {
00477 long int i,
00478 j,
00479 ipj1,
00480 ipj2;
00481
00482 double a,
00483 a1,
00484 dr[4][405],
00485 p1,
00486 p2,
00487 p3,
00488 p4,
00489 p5,
00490 p6,
00491 rr[4][110],
00492 te,
00493 x,
00494 z;
00495
00496 static long jd[6]={143,145,157,360,376,379};
00497
00498 static long ip[38]={7,9,12,13,14,16,18,19,20,21,22,44,45,49,50,
00499 52,53,54,55,56,57,58,59,60,66,67,78,83,84,87,88,95,96,97,100,
00500 101,103,104};
00501
00502 static long id[38]={7,3,10,27,23,49,51,64,38,47,60,103,101,112,
00503 120,114,143,145,157,152,169,183,200,163,225,223,237,232,235,
00504 249,247,300,276,278,376,360,379,384};
00505
00506 DEBUG_ENTRY( "rec_lines()" );
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523 for( i=0; i < 110; i++ )
00524 {
00525 rr[0][i] = P[0][i];
00526 rr[1][i] = P[1][i];
00527 rr[2][i] = P[2][i];
00528 z = P[0][i] - P[1][i] + 1.0;
00529 te = 1.0e-04*t/z/z;
00530 p1 = P[3][i];
00531 p2 = P[4][i];
00532 p3 = P[5][i];
00533 p4 = P[6][i];
00534 if( te < 0.004 )
00535 {
00536 a1 = p1*pow(0.004,p2)/(1.0 + p3*pow(0.004,p4));
00537 a = a1/sqrt(te/0.004);
00538 }
00539 else
00540 {
00541 if( te > 2.0 )
00542 {
00543 a1 = p1*pow(2.0,p2)/(1.0 + p3*pow(2.0,p4));
00544 a = a1/pow(te/2.0,1.5);
00545 }
00546 else
00547 {
00548 a = p1*pow(te,p2)/(1.0 + p3*pow(te,p4));
00549 }
00550 }
00551 rr[3][i] = 1.0e-13*z*a*P[7][i];
00552 }
00553
00554 for( i=0; i < 405; i++ )
00555 {
00556 dr[0][i] = ST[0][i];
00557 dr[1][i] = ST[1][i];
00558 dr[2][i] = ST[2][i];
00559 te = 1.0e-04*t;
00560 p1 = ST[3][i];
00561 p2 = ST[4][i];
00562 p3 = ST[5][i];
00563 p4 = ST[6][i];
00564 p5 = ST[7][i];
00565 p6 = ST[8][i];
00566 if( te < p6 )
00567 {
00568 x = p5*(1.0/te - 1.0/p6);
00569 if( x > 80.0 )
00570 {
00571 a = 0.0;
00572 }
00573 else
00574 {
00575 a1 = (p1/p6 + p2 + p3*p6 + p4*p6*p6)/pow(p6,1.5)/exp(p5/
00576 p6);
00577 a = a1/exp(x);
00578 }
00579 }
00580 else
00581 {
00582 if( te > 6.0 )
00583 {
00584 a1 = (p1/6.0 + p2 + p3*6.0 + p4*36.0)/pow(6.0,1.5)/
00585 exp(p5/6.0);
00586 a = a1/pow(te/6.0,1.5);
00587 }
00588 else
00589 {
00590 a = (p1/te + p2 + p3*te + p4*te*te)/pow(te,1.5)/exp(p5/
00591 te);
00592 }
00593 }
00594 dr[3][i] = 1.0e-12*a;
00595 }
00596
00597 for( i=0; i < 6; i++ )
00598 {
00599 ipj1 = jd[i];
00600 ipj2 = ipj1 + 1;
00601 dr[3][ipj1-1] += dr[3][ipj2-1];
00602 dr[0][ipj2-1] = 0.0;
00603 }
00604
00605 for( i=0; i < 38; i++ )
00606 {
00607 ipj1 = ip[i];
00608 ipj2 = id[i];
00609 rr[3][ipj1-1] += dr[3][ipj2-1];
00610 dr[0][ipj2-1] = 0.0;
00611 }
00612
00613 for( i=0; i < 110; i++ )
00614 {
00615 r[0][i] = (realnum)rr[0][i];
00616 r[1][i] = (realnum)rr[1][i];
00617 r[2][i] = (realnum)rr[2][i];
00618 r[3][i] = (realnum)rr[3][i];
00619 }
00620
00621 j = 110;
00622 for( i=0; i < 405; i++ )
00623 {
00624 if( dr[0][i] > 1.0 )
00625 {
00626 j += 1;
00627 r[0][j-1] = (realnum)dr[0][i];
00628 r[1][j-1] = (realnum)dr[1][i];
00629 r[2][j-1] = (realnum)dr[2][i];
00630 r[3][j-1] = (realnum)dr[3][i];
00631 }
00632 }
00633 return;
00634 }
00635
00636 double t_ADfA::rad_rec(long int iz,
00637 long int in,
00638 double t)
00639 {
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668 double tt;
00669 double rate;
00670
00671 DEBUG_ENTRY( "rad_rec()" );
00672
00673 rate = 0.0;
00674 if( iz < 1 || iz > 30 )
00675 {
00676 fprintf( ioQQQ, " rad_rec called with insane atomic number, =%4ld\n",
00677 iz );
00678 cdEXIT(EXIT_FAILURE);
00679 }
00680 if( in < 1 || in > iz )
00681 {
00682 fprintf( ioQQQ, " rad_rec called with insane number elec =%4ld\n",
00683 in );
00684 cdEXIT(EXIT_FAILURE);
00685 }
00686 if( (((in <= 3 || in == 11) || (iz > 5 && iz < 9)) || iz == 10) ||
00687 (iz == 26 && in > 11) )
00688 {
00689 tt = sqrt(t/rnew[in-1][iz-1][2]);
00690 rate =
00691 rnew[in-1][iz-1][0]/(tt*pow(tt + 1.0,1.0 - rnew[in-1][iz-1][1])*
00692 pow(1.0 + sqrt(t/rnew[in-1][iz-1][3]),1.0 + rnew[in-1][iz-1][1]));
00693 }
00694 else
00695 {
00696 tt = t*1.0e-04;
00697 if( iz == 26 && in <= 13 )
00698 {
00699 rate = fe[in-1][0]/pow(tt,fe[in-1][1] +
00700 fe[in-1][2]*log10(tt));
00701 }
00702 else
00703 {
00704 rate = rrec[in-1][iz-1][0]/pow(tt,(double)rrec[in-1][iz-1][1]);
00705 }
00706 }
00707
00708 return rate;
00709 }
00710
00711 double t_ADfA::H_rad_rec(long int iz,
00712 long int n,
00713 double t)
00714 {
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733 double rate,
00734 TeScaled,
00735 x,
00736 x1,
00737 x2;
00738
00739 DEBUG_ENTRY( "H_rad_rec()" );
00740
00741 rate = 0.0;
00742
00743
00744 ASSERT( iz > 0 );
00745
00746
00747 ASSERT( n < NHYDRO_MAX_LEVEL );
00748
00749 TeScaled = t/POW2((double)iz);
00750
00751 if( n < 0 )
00752 {
00753 x1 = sqrt(TeScaled/3.148);
00754 x2 = sqrt(TeScaled/7.036e05);
00755 rate = 7.982e-11/x1/pow(1.0 + x1,0.252)/pow(1.0 + x2,1.748);
00756 }
00757 else
00758 {
00759 x = log10(TeScaled);
00760 rate = (HRF[n][0] + HRF[n][2]*x + HRF[n][4]*
00761 x*x + HRF[n][6]*powi(x,3) + HRF[n][8]*powi(x,4))/
00762 (1.0 + HRF[n][1]*x + HRF[n][3]*x*x + HRF[n][5]*
00763 powi(x,3) + HRF[n][7]*powi(x,4));
00764 rate = pow(10.0,rate)/TeScaled;
00765 }
00766 rate *= iz;
00767
00768 return rate;
00769 }
00770
00771
00772
00773 double t_ADfA::coll_ion(
00774
00775 long int iz,
00776
00777 long int in,
00778
00779 double t)
00780 {
00781 double rate, te, u;
00782
00783 DEBUG_ENTRY( "coll_ion()" );
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797 rate = 0.0;
00798
00799 if( iz < 1 || iz > 30 )
00800 {
00801
00802 return( 0. );
00803 }
00804
00805 if( in < 1 || in > iz )
00806 {
00807
00808 return( 0. );
00809 }
00810
00811 te = t*EVRYD/TE1RYD;
00812 u = CF[in-1][iz-1][0]/te;
00813 if( u > 80.0 )
00814 {
00815 return( 0. );
00816 }
00817
00818 rate = (CF[in-1][iz-1][2]*(1.0 + CF[in-1][iz-1][1]*
00819 sqrt(u))/(CF[in-1][iz-1][3] + u)*pow(u,CF[in-1][iz-1][4])*
00820 exp(-u));
00821
00822 return(rate);
00823 }
00824
00825 realnum t_ADfA::h_coll_str( long ipLo, long ipHi, long ipTe )
00826 {
00827 realnum rate;
00828
00829 DEBUG_ENTRY( "h_coll_str()" );
00830
00831 ASSERT( ipLo < ipHi );
00832
00833 # if !defined NDEBUG
00834 long ipISO = ipH_LIKE;
00835 long nelem = ipHYDROGEN;
00836 # endif
00837 ASSERT( N_(ipLo) < N_(ipHi) );
00838 ASSERT( N_(ipHi) <= 5 );
00839
00840 rate = (realnum)HCS[ipHi-1][ipLo][ipTe];
00841
00842 return rate;
00843 }