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