00001
00002
00003
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
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
00241
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
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
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
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
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
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
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
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
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
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
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
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
00745 ASSERT( iz > 0 );
00746
00747
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
00773
00774 double t_ADfA::coll_ion(
00775
00776 long int iz,
00777
00778 long int in,
00779
00780 double t)
00781 {
00782 double rate, te, u;
00783
00784 DEBUG_ENTRY( "coll_ion()" );
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
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
00816 long int z,
00817
00818 long int n,
00819
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
00829 return( 0. );
00830 }
00831
00832 if( n < 0 || n > z )
00833 {
00834
00835 return( 0. );
00836 }
00837
00838
00839
00840
00841 if( atmdat.CIRCData == t_atmdat::DIMA)
00842 {
00843
00844
00845
00846
00847 rate = coll_ion( z+1, z+1-n, t );
00848 }
00849 else if( atmdat.CIRCData == t_atmdat::HYBRID)
00850 {
00851
00852 rate = coll_ion_hybrid( z, n, t );
00853 }
00854 else
00855 {
00856
00857 TotalInsanity();
00858 }
00859
00860 ASSERT(rate >= 0.0);
00861 return(rate);
00862 }
00863
00864
00865
00866 double t_ADfA::coll_ion_hybrid(
00867
00868 long int nelem,
00869
00870 long int ion,
00871
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
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 }