00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #include "cddefines.h"
00016 #include "phycon.h"
00017 #include "path.h"
00018 #include "elementnames.h"
00019 #include "atmdat.h"
00020 #include "ionbal.h"
00021 #include "dense.h"
00022
00023 static const int MAX_FIT_PAR_DR = 9;
00024
00025 static const int ATOMIC_NUM_BIG = 55;
00026
00027 static const int NELECTRONS = 12;
00028 static double C[ATOMIC_NUM_BIG][NELECTRONS][MAX_FIT_PAR_DR];
00029 static double E[ATOMIC_NUM_BIG][NELECTRONS][MAX_FIT_PAR_DR];
00030 static int defn_C[ATOMIC_NUM_BIG][NELECTRONS];
00031 static int defn_E[ATOMIC_NUM_BIG][NELECTRONS];
00032 static int nn_counter[ATOMIC_NUM_BIG][NELECTRONS];
00033
00034 static const int MAX_FIT_PAR_RR = 6;
00035 static double par[ATOMIC_NUM_BIG][NELECTRONS][MAX_FIT_PAR_RR];
00036 static int defn[ATOMIC_NUM_BIG][NELECTRONS];
00037
00038
00039
00040
00041 #if defined(PRINT_DR) || defined(PRINT_RR)
00042 static const char FILE_NAME_OUT[] = "array.out";
00043 #endif
00044
00045
00046 void ion_recom_calculate( void )
00047 {
00048 static double TeUsed = -1;
00049 long int ion ,
00050 nelem ,
00051 i;
00052
00053 DEBUG_ENTRY( "ion_recom_calculate()" );
00054
00055 if( fabs(phycon.te/TeUsed - 1. ) > 0.001 )
00056 {
00057
00058 long int nsumrec[LIMELM];
00059 for( ion=0; ion<LIMELM; ++ion )
00060 {
00061 ionbal.DR_Badnell_rate_coef_mean_ion[ion] = 0.;
00062 nsumrec[ion] = 0;
00063 }
00064 TeUsed = phycon.te;
00065
00066
00067
00068 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00069 {
00070 for( ion=0; ion<nelem+1; ++ion )
00071 {
00072 long int n_bnd_elec_before_recom ,
00073 n_bnd_elec_after_recom;
00074
00075 n_bnd_elec_before_recom = nelem-ion;
00076 n_bnd_elec_after_recom = nelem-ion+1;
00077 # define DR2SMALL 1e-15
00078
00079 if( (ionbal.DR_Badnell_rate_coef[nelem][ion] = Badnell_DR_rate_eval(
00080
00081 nelem+1,
00082
00083 n_bnd_elec_before_recom )) >= 0. )
00084 {
00085 ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = true;
00086 if( ionbal.DR_Badnell_rate_coef[nelem][ion] > DR2SMALL)
00087 {
00088
00089 ++nsumrec[ion];
00090 ionbal.DR_Badnell_rate_coef_mean_ion[ion] +=
00091 log10(ionbal.DR_Badnell_rate_coef[nelem][ion]);
00092 }
00093 }
00094 else
00095 {
00096 ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = false;
00097 ionbal.DR_Badnell_rate_coef[nelem][ion] = 0.;
00098 }
00099
00100 if( (ionbal.RR_Badnell_rate_coef[nelem][ion] = Badnell_RR_rate_eval(
00101
00102 nelem+1,
00103
00104 n_bnd_elec_before_recom )) >= 0. )
00105 {
00106 ionbal.lgRR_Badnell_rate_coef_exist[nelem][ion] = true;
00107 }
00108 else
00109 {
00110 ionbal.lgRR_Badnell_rate_coef_exist[nelem][ion] = false;
00111 ionbal.RR_Badnell_rate_coef[nelem][ion] = 0.;
00112 }
00113
00114
00115 ionbal.RR_Verner_rate_coef[nelem][ion] =
00116 t_ADfA::Inst().rad_rec(
00117
00118 nelem+1 ,
00119
00120 n_bnd_elec_after_recom ,
00121 phycon.te );
00122 }
00123 }
00124
00125
00126
00127
00128
00129 double fitSum = 0;
00130
00131 double fe13_c[10]={ 3.55e-4, 2.40e-3, 7.83e-3, 1.10e-2, 3.30e-2, 1.45e-1, 8.50e-2, 2.59e-2, 8.93e-3, 9.80e-3 },
00132 fe13_E[10]={ 2.19e-2, 1.79e-1, 7.53e-1, 2.21e0, 9.57e0, 3.09e1, 6.37e1, 2.19e2, 1.50e3, 7.86e3 };
00133
00134 double Fe_Gu_c[9][6] = {
00135 { 2.50507e-11, 5.60226e-11, 1.85001e-10, 3.57495e-9, 1.66321e-7, 0. },
00136 { 9.19610e-11, 2.92460e-10, 1.02120e-9, 1.14852e-8, 3.25418e-7, 0. },
00137 { 9.02625e-11, 6.22962e-10, 5.77545e-9, 1.78847e-8, 3.40610e-7, 0. },
00138 { 9.04286e-12, 9.68148e-10, 4.83636e-9, 2.48159e-8, 3.96815e-7, 0. },
00139 { 6.77873e-10, 1.47252e-9, 5.31202e-9, 2.54793e-8, 3.47407e-7, 0. },
00140 { 1.29742e-9, 4.10172e-9, 1.23605e-8, 2.33615e-8, 2.97261e-7, 0. },
00141 { 8.78027e-10, 2.31680e-9, 3.49333e-9, 1.16927e-8, 8.18537e-8, 1.54740e-7 },
00142 { 2.23178e-10, 1.87313e-9, 2.86171e-9, 1.38575e-8, 1.17803e-7, 1.06251e-7 },
00143 { 2.17263e-10, 7.35929e-10, 2.81276e-9, 1.32411e-8, 1.15761e-7, 4.80389e-8 }
00144 },
00145
00146 Fe_Gu_E[9][6] = {
00147 { 8.30501e-2, 8.52897e-1, 3.40225e0, 2.23053e1, 6.80367e1, 0. },
00148 { 1.44392e-1, 9.23999e-1, 5.45498e0, 2.04301e1, 7.06112e1, 0. },
00149 { 5.79132e-2, 1.27852e0, 3.22439e0, 1.79602e1, 6.96277e1, 0. },
00150 { 1.02421e-1, 1.79393e0, 4.83226e0, 1.91117e1, 6.80858e1, 0. },
00151 { 1.24630e-1, 6.86045e-1, 3.09611e0, 1.44023e1, 6.42820e1, 0. },
00152 { 1.34459e-1, 6.63028e-1, 2.61753e0, 1.30392e1, 6.10222e1, 0. },
00153 { 7.79748e-2, 5.35522e-1, 1.88407e0, 8.38459e0, 3.38613e1, 7.89706e1 },
00154 { 8.83019e-2, 6.12756e-1, 2.36035e0, 9.61736e0, 3.64467e1, 8.72406e1 },
00155 { 1.51322e-1, 5.63155e-1, 2.57013e0, 9.08166e0, 3.69528e1, 1.08067e2 }
00156 };
00157
00158 double s_c[5][2] = {
00159 {0.1565e-3, 0.1617e-2},
00160 {0.3026e-3, 0.2076e-1},
00161 {0.3177e-1, 0.6309e-3},
00162 {0.2464e-1, 0.5553e-3},
00163 {0.1924e-1, 0.5111e-3}
00164 },
00165
00166 s_E[5][2] = {
00167 {0.1157e6, 0.1868e6},
00168 {0.9662e5, 0.1998e6},
00169 {0.1928e6, 0.6126e5},
00170 {0.1806e6, 0.3142e5},
00171 {0.1519e6, 0.4978e5}
00172 };
00173
00174
00175 double te_eV32 = pow( phycon.te_eV, 1.5);
00176
00177
00178 nelem = ipIRON;
00179 ion = 12;
00180 if( !ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] )
00181 {
00182 for( i=0; i<10; i++ )
00183 {
00184 fitSum += fe13_c[i] * sexp( fe13_E[i]/phycon.te_eV );
00185 }
00186
00187
00188 ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = true;
00189 ionbal.DR_Badnell_rate_coef[nelem][ion] = fitSum / phycon.te32;
00190 ++nsumrec[ion];
00191 ionbal.DR_Badnell_rate_coef_mean_ion[ion] +=
00192 log10(ionbal.DR_Badnell_rate_coef[nelem][ion]);
00193 }
00194
00195
00196
00197
00198
00199 for( ion=0; ion<9; ion++ )
00200 {
00201
00202 if( !ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion+5] )
00203 {
00204 fitSum = 0;
00205 for( i=0; i<6; i++ )
00206 {
00207 fitSum += Fe_Gu_c[ion][i] * sexp( Fe_Gu_E[ion][i]/phycon.te_eV );
00208 }
00209 ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion+5] = true;
00210 ionbal.DR_Badnell_rate_coef[nelem][ion+5] = fitSum / te_eV32;
00211 if( ionbal.DR_Badnell_rate_coef[nelem][ion+5] > DR2SMALL )
00212 {
00213 ++nsumrec[ion+5];
00214 ionbal.DR_Badnell_rate_coef_mean_ion[ion+5] +=
00215 log10(ionbal.DR_Badnell_rate_coef[nelem][ion+5]);
00216 }
00217 }
00218 }
00219
00220
00221
00222 for( ion=0; ion<LIMELM; ++ion )
00223 {
00224 if( nsumrec[ion] > 0 )
00225 ionbal.DR_Badnell_rate_coef_mean_ion[ion] =
00226 pow(10., ionbal.DR_Badnell_rate_coef_mean_ion[ion]/nsumrec[ion]);
00227
00228 }
00229
00230
00231 nelem = ipSULPHUR;
00232
00233 for( ion=0; ion<5; ion++ )
00234 {
00235
00236
00237
00238 if( !ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] )
00239 {
00240
00241
00242 fitSum = 0;
00243 for( i=0; i<2; i++ )
00244 {
00245 fitSum += s_c[ion][i] * sexp( s_E[ion][i]/phycon.te );
00246 }
00247 ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = true;
00248 ionbal.DR_Badnell_rate_coef[nelem][ion] = fitSum / phycon.te32;
00249
00250
00251
00252
00253
00254
00255 if( ionbal.nDR_S_guess==0 )
00256 {
00257
00258 ionbal.DR_Badnell_rate_coef[nelem][ion] =
00259 MAX2( ionbal.DR_Badnell_rate_coef[nelem][ion] ,
00260 ionbal.DR_Badnell_rate_coef_mean_ion[ion]);
00261 }
00262 else if( ionbal.nDR_S_guess==1 )
00263 {
00264
00265 ionbal.DR_Badnell_rate_coef[nelem][ion] =
00266 ionbal.DR_Badnell_rate_coef[nelem][ion];
00267 }
00268 else if( ionbal.nDR_S_guess==2 )
00269 {
00270
00271 ionbal.DR_Badnell_rate_coef[nelem][ion] =
00272 ionbal.DR_Badnell_rate_coef[ipOXYGEN][ion]*
00273 ionbal.DR_S_scale[ion];
00274 }
00275 else
00276 TotalInsanity();
00277 }
00278 }
00279
00280
00281 if( ionbal.lgRecom_Badnell_print )
00282 {
00283 fprintf(ioQQQ,"DEBUG Badnell recombination RR, then DR, T=%.3e\n", phycon.te );
00284 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00285 {
00286 fprintf(ioQQQ,"nelem=%li %s, RR then DR\n",
00287 nelem , elementnames.chElementNameShort[nelem] );
00288 for( ion=0; ion<nelem+1; ++ion )
00289 {
00290 if( ionbal.lgRR_Badnell_rate_coef_exist[nelem][ion] )
00291 {
00292 fprintf(ioQQQ," %.2e", ionbal.RR_Badnell_rate_coef[nelem][ion] );
00293 }
00294 else
00295 {
00296 fprintf(ioQQQ," %.2e", -1. );
00297 }
00298 }
00299 fprintf(ioQQQ,"\n" );
00300 for( ion=0; ion<nelem+1; ++ion )
00301 {
00302 if( ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] )
00303 {
00304 fprintf(ioQQQ," %.2e", ionbal.DR_Badnell_rate_coef[nelem][ion] );
00305 }
00306 else
00307 {
00308 fprintf(ioQQQ," %.2e", -1. );
00309 }
00310 }
00311 fprintf(ioQQQ,"\n\n" );
00312 }
00313
00314 fprintf(ioQQQ,"mean DR recombination ion mean stddev stddev/mean \n" );
00315 for( ion=0; ion<LIMELM; ++ion )
00316 {
00317 double stddev;
00318 stddev = 0.;
00319 for( nelem=ion; nelem<LIMELM; ++nelem )
00320 {
00321 stddev += POW2(
00322 ionbal.DR_Badnell_rate_coef[nelem][ion]-
00323 ionbal.DR_Badnell_rate_coef_mean_ion[ion] );
00324 }
00325 stddev = sqrt( stddev / MAX2( 1 , nsumrec[ion]-1 ) );
00326 fprintf(ioQQQ," %2li %.2e %.2e %.2e\n",
00327 ion ,
00328 ionbal.DR_Badnell_rate_coef_mean_ion[ion] ,
00329 stddev ,
00330 stddev / SDIV(ionbal.DR_Badnell_rate_coef_mean_ion[ion]) );
00331 }
00332 cdEXIT( EXIT_SUCCESS );
00333 }
00334 }
00335
00336 DEBUG_EXIT( "ion_recom_calculate()" );
00337
00338 return;
00339 }
00340
00341
00342
00343
00344 void Badnell_rec_init( void )
00345 {
00346
00347 double par_C[MAX_FIT_PAR_DR];
00348 double par_E[MAX_FIT_PAR_DR];
00349 char chLine[INPUT_LINE_LENGTH];
00350 int z_value, n_value;
00351 int i, j, k;
00352 int count, number;
00353 double temp_par[MAX_FIT_PAR_RR];
00354 int M_state, W_state;
00355
00356
00357 const int NBLOCK = 2;
00358 int data_begin_line[NBLOCK];
00359 int length_of_line;
00360 FILE *ioDATA;
00361 char chFilename[FILENAME_PATH_LENGTH_2];
00362 int yr, mo, dy;
00363 char *chs;
00364
00365 DEBUG_ENTRY( "Badnell_rec_init()" );
00366
00367 # if defined(PRINT_DR) || defined(PRINT_RR)
00368 FILE *ofp;
00369 ofp= fopen(FILE_NAME_OUT, "w");
00370 if( ofp == NULL )
00371 {
00372 printf("Can't open %s\n", FILE_NAME_OUT);
00373 exit(1);
00374 }
00375 # endif
00376
00377 for( i=0; i<NBLOCK; ++i )
00378 {
00379
00380 data_begin_line[i] = INT_MIN;
00381 }
00382
00383
00384
00385 if( lgDataPathSet )
00386 {
00387
00388 strcpy( chFilename , chDataPath );
00389 strcat( chFilename , "badnell_dr.dat" );
00390 }
00391 else
00392 {
00393
00394 strcpy( chFilename , "badnell_dr.dat" );
00395 }
00396
00397 if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
00398 {
00399 printf("Badnell_rec_init cannot open %s\n", chFilename);
00400
00401 path_not_set();
00402 cdEXIT(EXIT_FAILURE);
00403 }
00404
00405 count = 0;
00406 number = 0;
00407
00408
00409
00410 while( fgets(chLine, (int)sizeof(chLine), ioDATA) != NULL )
00411 {
00412 count++;
00413
00414 if( chLine[2]=='Z' )
00415 {
00416
00417
00418 data_begin_line[number] = count;
00419 ASSERT( number < NBLOCK );
00420 number++;
00421 }
00422 }
00423
00424
00425 for( i=0; i<ATOMIC_NUM_BIG; i++ )
00426 {
00427 for( j=0; j<NELECTRONS; j++ )
00428 {
00429 defn_C[i][j] = 0;
00430 defn_E[i][j] = 0;
00431 }
00432 }
00433
00434
00435 for( i=0; i<ATOMIC_NUM_BIG; i++ )
00436 {
00437 for( j=0; j<NELECTRONS; j++ )
00438 {
00439 for( k=0; k<MAX_FIT_PAR_DR; k++ )
00440 {
00441 C[i][j][k] = 0;
00442 E[i][j][k] = 0;
00443 }
00444 }
00445 }
00446
00447 count = 0;
00448
00449 fseek(ioDATA, 0, SEEK_SET);
00450
00451 fgets(chLine, (int)sizeof(chLine), ioDATA);
00452 count++;
00453 {
00454 int dr_yr = 2006, dr_mo = 3, dr_dy = 11;
00455 chs = strchr(chLine, ')');
00456 ++chs;
00457 sscanf(chs, "%4i%2i%2i",&yr, &mo, &dy);
00458 if((yr != dr_yr) || (mo != dr_mo) || (dy != dr_dy))
00459 {
00460 printf("The version of %s I found (%i %i %i) is not the current version (%i %i %i).\n",
00461 chFilename, yr, mo, dy, dr_yr, dr_mo, dr_dy);
00462 cdEXIT(EXIT_FAILURE);
00463 }
00464 }
00465
00466 while( fgets(chLine, (int)sizeof(chLine), ioDATA) != NULL )
00467 {
00468 count++;
00469 length_of_line = (int)strlen(chLine);
00470
00471
00472 if( count > data_begin_line[0] && count < data_begin_line[1] && length_of_line >3 )
00473 {
00474
00475 for( i=0; i<MAX_FIT_PAR_DR; i++ )
00476 {
00477 par_C[i] = 0;
00478 }
00479 sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf%lf%lf%lf",
00480 &z_value, &n_value, &M_state, &W_state, &par_C[0], &par_C[1], &par_C[2],
00481 &par_C[3], &par_C[4], &par_C[5], &par_C[6], &par_C[7], &par_C[8]);
00482
00483 if(M_state == 1)
00484 {
00485
00486 defn_C[z_value][n_value] = 1;
00487
00488
00489 nn_counter[z_value][n_value] = 8;
00490 while( par_C[nn_counter[z_value][n_value]] == 0 )
00491 {
00492 --nn_counter[z_value][n_value];
00493 if( nn_counter[z_value][n_value] < 0 )
00494 break;
00495 }
00496 ++nn_counter[z_value][n_value];
00497
00498
00499 for( i=0; i<nn_counter[z_value][n_value]; i++ )
00500 C[z_value][n_value][i] = par_C[i];
00501 }
00502 }
00503 }
00504
00505
00506 fseek(ioDATA, 0, SEEK_SET);
00507 count = 0;
00508 while( fgets(chLine, (int)sizeof(chLine), ioDATA) != NULL )
00509 {
00510 count++;
00511 length_of_line = (int)strlen(chLine);
00512 if( count > data_begin_line[1] && length_of_line > 3 )
00513 {
00514
00515
00516 for( i=0; i<MAX_FIT_PAR_DR; i++ )
00517 {
00518 par_E[i] = 0;
00519 }
00520 sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf%lf%lf%lf",
00521 &z_value, &n_value, &M_state, &W_state, &par_E[0], &par_E[1], &par_E[2],
00522 &par_E[3], &par_E[4], &par_E[5], &par_E[6], &par_E[7], &par_E[8]);
00523
00524 if(M_state == 1)
00525 {
00526
00527 nn_counter[z_value][n_value] = 8;
00528 while( par_E[nn_counter[z_value][n_value]] == 0 )
00529 {
00530 --nn_counter[z_value][n_value];
00531
00532 if( nn_counter[z_value][n_value] < 0 )
00533 break;
00534 }
00535 ++nn_counter[z_value][n_value];
00536
00537
00538 defn_E[z_value][n_value] = 1;
00539
00540
00541 for( i=0; i<nn_counter[z_value][n_value]; i++ )
00542 E[z_value][n_value][i] = par_E[i];
00543 }
00544 }
00545 }
00546
00547
00548 # ifdef PRINT_DR
00549 for( i=0; i<ATOMIC_NUM_BIG; i++ )
00550 {
00551 for( j=0; j<NELECTRONS; j++ )
00552 {
00553 if( defn_C[i][j] == 1 )
00554 {
00555 fprintf(ofp, "%i %i %e %e %e %e %e %e %e %e %e\n",
00556 i, j, C[i][j][0], C[i][j][1], C[i][j][2], C[i][j][3], C[i][j][4],
00557 C[i][j][5], C[i][j][6], C[i][j][7], C[i][j][8]);
00558 }
00559 }
00560 }
00561 for( i=0; i<ATOMIC_NUM_BIG; i++ )
00562 {
00563 for( j=0; j<NELECTRONS; j++ )
00564 {
00565 if( defn_E[i][j] == 1 )
00566 {
00567 fprintf(ofp, "%i %i %e %e %e %e %e %e %e %e %e\n",
00568 i, j, E[i][j][0], E[i][j][1], E[i][j][2], E[i][j][3], E[i][j][4],
00569 E[i][j][5], E[i][j][6], E[i][j][7], E[i][j][8]);
00570 }
00571 }
00572 }
00573 fclose(ofp);
00574 # endif
00575
00576
00577 for( i=0; i<ATOMIC_NUM_BIG; i++ )
00578 {
00579 for( j=0; j<NELECTRONS; j++ )
00580 {
00581 if( defn_C[i][j] != defn_E[i][j] )
00582 {
00583 printf("C%i, E%i: %i %i\n", i, j, defn_C[i][j], defn_E[i][j]);
00584 printf("Critical error in defn\n");
00585 }
00586 }
00587 }
00588
00589
00590
00591 DEBUG_EXIT( "Badnell_rec_init()" );
00592
00593
00594
00595
00596 if( lgDataPathSet )
00597 {
00598
00599 strcpy( chFilename , chDataPath );
00600 strcat( chFilename , "badnell_rr.dat" );
00601 }
00602 else
00603 {
00604
00605 strcpy( chFilename , "badnell_rr.dat" );
00606 }
00607
00608 if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
00609 {
00610 printf("Badnell_rec_init cannot open %s\n", chFilename);
00611 exit(1);
00612 }
00613
00614
00615 for( i=0; i<ATOMIC_NUM_BIG; i++ )
00616 {
00617 for( j=0; j<NELECTRONS; j++ )
00618 {
00619 defn[i][j] = 0;
00620 }
00621 }
00622
00623
00624 for( i=0; i<ATOMIC_NUM_BIG; i++ )
00625 {
00626 for( j=0; j<NELECTRONS; j++ )
00627 {
00628 for( k=0; k<MAX_FIT_PAR_RR; k++ )
00629 {
00630 par[i][j][k] = 0;
00631 }
00632 }
00633 }
00634
00635
00636 {
00637 int rr_yr = 2005, rr_mo = 12, rr_dy = 27;
00638 fgets(chLine, (int)sizeof(chLine), ioDATA);
00639 chs = strchr(chLine, ')');
00640 ++chs;
00641 sscanf(chs, "%4i%2i%2i", &yr, &mo, &dy);
00642 if((yr != rr_yr)||(mo != rr_mo)||(dy != rr_dy))
00643 {
00644 printf("The version of %s I found (%i %i %i) is not the current version (%i %i %i).\n",
00645 chFilename, yr, mo, dy, rr_yr, rr_mo, rr_dy);
00646 cdEXIT(EXIT_FAILURE);
00647 }
00648 }
00649
00650 while( fgets(chLine, (int)sizeof(chLine), ioDATA) != NULL )
00651 {
00652
00653
00654 for( i=0; i<MAX_FIT_PAR_RR; i++ )
00655 {
00656 temp_par[i] = 0;
00657 }
00658 if(chLine[0] != '#')
00659 {
00660 sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf",
00661 &z_value, &n_value, &M_state, &W_state, &temp_par[0], &temp_par[1],
00662 &temp_par[2], &temp_par[3], &temp_par[4], &temp_par[5]);
00663
00664 if(M_state == 1)
00665 {
00666
00667 defn[z_value][n_value] = 1;
00668
00669 for( i=0; i<MAX_FIT_PAR_RR; i++ )
00670 par[z_value][n_value][i] = temp_par[i];
00671 }
00672 }
00673 }
00674
00675
00676 # ifdef PRINT_RR
00677 count = 0;
00678 for( i=0; i<ATOMIC_NUM_BIG; i++ )
00679 {
00680 for( j=0; j<NELECTRONS; j++ )
00681 {
00682 if( defn[i][j] == 1 )
00683 {
00684 fprintf(ofp, "%i %i %e %e %e %e %e %e\n",
00685 i, j, par[i][j][0], par[i][j][1], par[i][j][2], par[i][j][3],
00686 par[i][j][4], par[i][j][5]);
00687 count++;
00688 }
00689 }
00690 }
00691 fprintf(ofp, "total lines are %i ", count);
00692
00693 fclose(ofp);
00694 # endif
00695
00696 fclose(ioDATA);
00697
00698 {
00699
00700 enum {DEBUG_LOC=false};
00701
00702 if( DEBUG_LOC )
00703 {
00704 long int nelem;
00705 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00706 {
00707 int ion;
00708 fprintf(ioQQQ,"\nDEBUG rr rec\t%li",nelem);
00709 for( ion=0; ion<=nelem; ++ion )
00710 {
00711 fprintf(ioQQQ,"\t%.2e", Badnell_RR_rate_eval(nelem+1 , nelem-ion ) );
00712 }
00713 fprintf(ioQQQ,"\n");
00714 fprintf(ioQQQ,"DEBUG dr rec\t%li",nelem);
00715 for( ion=0; ion<=nelem; ++ion )
00716 {
00717 fprintf(ioQQQ,"\t%.2e", Badnell_DR_rate_eval(nelem+1 , nelem-ion ) );
00718 }
00719 fprintf(ioQQQ,"\n");
00720 }
00721 cdEXIT(EXIT_FAILURE);
00722 }
00723 }
00724
00725 return;
00726 }
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736 double Badnell_DR_rate_eval(
00737
00738 int n_atomic_number,
00739
00740 int n_core_e_before_recomb )
00741 {
00742
00743 double k, sum;
00744 int i;
00745
00746 DEBUG_ENTRY( "Badnell_DR_rate_eval()" );
00747
00748 # if 0
00749 if( n_atomic_number == 26 && n_core_e_before_recomb==13 )
00750 {
00751
00752
00753
00754
00755
00756 double cFe13[8] =
00757 {1.090e-3, 7.801e-3, 1.132e-2, 4.740e-2, 1.990e-1, 3.379e-2, 1.140e-1, 1.250e-1};
00758
00759 double EFe13[8] =
00760 {1.246e3, 1.063e4, 4.719e4, 1.952e5, 5.637e5, 2.248e6, 7.202e6, 3.999e9};
00761
00762 sum = 0;
00763 i = 0;
00764
00765 for(i=0; i<8; ++i )
00766 {
00767 sum += (cFe13[i] * sexp( EFe13[i]/phycon.te));
00768 }
00769
00770
00771 k = sum / phycon.te32;
00772
00773 DEBUG_EXIT( "Badnell_DR_rate_eval()" );
00774
00775 return k;
00776
00777 }
00778 # endif
00779 if( n_atomic_number==26 && n_core_e_before_recomb>=12 && n_core_e_before_recomb<=18 )
00780 {
00781
00782
00783
00784
00785
00786 double cFe_q[7][8] =
00787 {
00788 {5.636e-4, 7.390e-3, 3.635e-2, 1.693e-1, 3.315e-2, 2.288e-1, 7.316e-2, 0.},
00789 {1.090e-3, 7.801e-3, 1.132e-2, 4.740e-2, 1.990e-1, 3.379e-2, 1.140e-1, 1.250e-1},
00790 {3.266e-3, 7.637e-3, 1.005e-2, 2.527e-2, 6.389e-2, 1.564e-1, 0., 0.},
00791 {1.074e-3, 6.080e-3, 1.887e-2, 2.540e-2, 7.580e-2, 2.773e-1, 0., 0.},
00792 {9.073e-4, 3.777e-3, 1.027e-2, 3.321e-2, 8.529e-2, 2.778e-1, 0., 0.},
00793 {5.335e-4, 1.827e-3, 4.851e-3, 2.710e-2, 8.226e-2, 3.147e-1, 0., 0.},
00794 {7.421e-4, 2.526e-3, 4.605e-3, 1.489e-2, 5.891e-2, 2.318e-1, 0., 0.}
00795 };
00796
00797
00798 double EFe_q[7][8] =
00799 {
00800 {3.628e3, 2.432e4, 1.226e5, 4.351e5, 1.411e6, 6.589e6, 1.030e7, 0},
00801 {1.246e3, 1.063e4, 4.719e4, 1.952e5, 5.637e5, 2.248e6, 7.202e6, 3.999e9},
00802 {1.242e3, 1.001e4, 4.466e4, 1.497e5, 3.919e5, 6.853e5, 0. , 0.},
00803 {1.387e3, 1.048e4, 3.955e4, 1.461e5, 4.010e5, 7.208e5, 0. , 0.},
00804 {1.525e3, 1.071e4, 4.033e4, 1.564e5, 4.196e5, 7.580e5, 0. , 0.},
00805 {2.032e3, 1.018e4, 4.638e4, 1.698e5, 4.499e5, 7.880e5, 0. , 0.},
00806 {3.468e3, 1.353e4, 3.690e4, 1.957e5, 4.630e5, 8.202e5, 0. , 0.}
00807 };
00808
00809 long int nion = n_core_e_before_recomb - 12;
00810 ASSERT( nion>=0 && nion <=6 );
00811
00812 sum = 0;
00813 i = 0;
00814
00815 for(i=0; i<8; ++i )
00816 {
00817 sum += (cFe_q[nion][i] * sexp( EFe_q[nion][i]/phycon.te));
00818 }
00819
00820
00821 k = sum / phycon.te32;
00822
00823 DEBUG_EXIT( "Badnell_DR_rate_eval()" );
00824
00825 return k;
00826
00827
00828 }
00829
00830 else if( n_atomic_number <= n_core_e_before_recomb )
00831 {
00832 k = -2;
00833 }
00834
00835 else if( n_atomic_number < 2 || n_atomic_number > 30 ||
00836 n_core_e_before_recomb < 1 || n_core_e_before_recomb > 11 )
00837 {
00838 k = -2;
00839 }
00840
00841 else if( defn_C[n_atomic_number][n_core_e_before_recomb] == 0 )
00842 {
00843 k = -1;
00844 }
00845 else if( defn_C[n_atomic_number][n_core_e_before_recomb] == 1 )
00846 {
00847 sum = 0;
00848 i = 0;
00849
00850 for(i=0; i<nn_counter[n_atomic_number][n_core_e_before_recomb]; ++i )
00851 {
00852 sum += (C[n_atomic_number][n_core_e_before_recomb][i] *
00853 sexp( E[n_atomic_number][n_core_e_before_recomb][i]/phycon.te));
00854 }
00855
00856
00857 k = sum / phycon.te32;
00858 }
00859
00860 else
00861 {
00862 k = -99;
00863 }
00864
00865 DEBUG_EXIT( "Badnell_DR_rate_eval()" );
00866
00867 return k;
00868
00869 }
00870
00871 double Badnell_RR_rate_eval(
00872
00873 int n_atomic_number,
00874
00875 int n_core_e_before_recomb )
00876 {
00877 double k;
00878 double B, D, F;
00879
00880 DEBUG_ENTRY( "Badnell_DR_rate_eval()" );
00881
00882 # if 0
00883 if( n_atomic_number == 26 && n_core_e_before_recomb==13 )
00884 {
00885
00886
00887
00888 double parFe13[6] =
00889 {4.321e-10, 0.6091, 2.255e03, 4.962e07, 0.0356, 1.006e05};
00890
00891
00892 double temp;
00893
00894 temp = -parFe13[5]/phycon.te;
00895 B = parFe13[1] + parFe13[4]*exp(temp);
00896 D = sqrt(phycon.te/parFe13[2]);
00897 F = sqrt(phycon.te/parFe13[3]);
00898 k = parFe13[0]/(D*pow((1.+D),(1.-B))*pow((1.+F),(1.+B)));
00899
00900 DEBUG_EXIT( "Badnell_RR_rate_eval()" );
00901
00902 return k;
00903 }
00904 # endif
00905 if( n_atomic_number==26 && n_core_e_before_recomb>=12 && n_core_e_before_recomb<=18 )
00906 {
00907
00908 double parFeq[7][6] ={
00909 {1.179e-9 , 0.7096, 4.508e2, 3.393e7, 0.0154, 3.977e6},
00910 {1.050e-9 , 0.6939, 4.568e2, 3.987e7, 0.0066, 5.451e5},
00911 {9.832e-10, 0.7146, 3.597e2, 3.808e7, 0.0045, 3.952e5},
00912 {8.303e-10, 0.7156, 3.531e2, 3.554e7, 0.0132, 2.951e5},
00913 {1.052e-9 , 0.7370, 1.639e2, 2.924e7, 0.0224, 4.291e5},
00914 {1.338e-9 , 0.7495, 7.242e1, 2.453e7, 0.0404, 4.199e5},
00915 {1.263e-9 , 0.7532, 5.209e1, 2.169e7, 0.0421, 2.917e5}
00916 };
00917
00918 double temp;
00919
00920 long int nion = n_core_e_before_recomb - 12;
00921 ASSERT( nion>=0 && nion <=6 );
00922
00923 temp = -parFeq[nion][5]/phycon.te;
00924 B = parFeq[nion][1] + parFeq[nion][4]*exp(temp);
00925 D = sqrt(phycon.te/parFeq[nion][2]);
00926 F = sqrt(phycon.te/parFeq[nion][3]);
00927 k = parFeq[nion][0]/(D*pow((1.+D),(1.-B))*pow((1.+F),(1.+B)));
00928
00929 DEBUG_EXIT( "Badnell_RR_rate_eval()" );
00930
00931 return k;
00932 }
00933
00934
00935
00936 else if( n_atomic_number <= n_core_e_before_recomb )
00937 {
00938 k = -2;
00939 }
00940
00941 else if( n_atomic_number < 1 || n_atomic_number > 30 ||
00942 n_core_e_before_recomb < 0 || n_core_e_before_recomb > 11 )
00943 {
00944 k = -2;
00945 }
00946
00947 else if( defn[n_atomic_number][n_core_e_before_recomb] == 0 )
00948 {
00949 k = -1;
00950 }
00951
00952 else if( defn[n_atomic_number][n_core_e_before_recomb] == 1 )
00953 {
00954
00955
00956
00957 double temp;
00958
00959 temp = -par[n_atomic_number][n_core_e_before_recomb][5]/phycon.te;
00960 B = par[n_atomic_number][n_core_e_before_recomb][1] + par[n_atomic_number][n_core_e_before_recomb][4]*exp(temp);
00961 D = sqrt(phycon.te/par[n_atomic_number][n_core_e_before_recomb][2]);
00962 F = sqrt(phycon.te/par[n_atomic_number][n_core_e_before_recomb][3]);
00963 k = par[n_atomic_number][n_core_e_before_recomb][0]/(D*pow((1.+D),(1.-B))*pow((1.+F),(1.+B)));
00964
00965 }
00966
00967
00968 else
00969 {
00970 k = -99;
00971 }
00972
00973 DEBUG_EXIT( "Badnell_DR_rate_eval()" );
00974 return k;
00975 }