00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #include "cddefines.h"
00016 #include "phycon.h"
00017 #include "elementnames.h"
00018 #include "atmdat.h"
00019 #include "iso.h"
00020 #include "ionbal.h"
00021 #include "dense.h"
00022 #include "taulines.h"
00023
00024 static const int MAX_FIT_PAR_DR = 9;
00025 static double ***DRFitParPart1;
00026 static double ***DRFitParPart2;
00027 static int **nDRFitPar;
00028
00029 static const int MAX_FIT_PAR_RR = 6;
00030 static double ***RRFitPar;
00031
00032
00033 static bool **lgDRBadnellDefined ,
00034 **lgDRBadnellDefinedPart2,
00035 **lgRRBadnellDefined,
00036 **lgDR_BadWeb_exist;
00037 static bool lgMustMallocRec=true;
00038 static double RecNoise[LIMELM],
00039 DR_Badnell_rate_coef_mean_ion[LIMELM];
00040
00041 static char chDRDataSource[LIMELM][LIMELM][10];
00042 static char chRRDataSource[LIMELM][LIMELM][10];
00043
00044
00045
00046
00047
00048 #if defined(PRINT_DR) || defined(PRINT_RR)
00049 static const char FILE_NAME_OUT[] = "array.out";
00050 #endif
00051
00052
00053
00054
00055
00056
00057
00058
00059 STATIC double CollisSuppres(
00060
00061 long int atomic_number,
00062
00063 long int ionic_charge,
00064
00065 double eden,
00066
00067 double T )
00068 {
00069
00070
00071 const double mu = 0.372;
00072 const double w = 4.969;
00073 const double x_a0 = 0.608;
00074
00075
00076
00077
00078 const double T_0 = 1.e5;
00079 const double q_0 = 3.;
00080
00081
00082
00083
00084
00085 const double c = 10.0;
00086
00087 double s, snew, x_a, E_c;
00088 long int iso_sequence;
00089
00090 eden = log10(eden);
00091
00092 iso_sequence = atomic_number - ionic_charge;
00093 ASSERT( iso_sequence>=0 );
00094
00095 x_a = x_a0 + log10( pow( (ionic_charge/q_0), 7. ) * sqrt( T/T_0 ) );
00096
00097
00098 if( eden >= x_a )
00099 {
00100 s = ( mu/( 1. + pow((eden-x_a)/w, 2.) ) +
00101 (1. - mu) * exp( -LN_TWO * pow((eden-x_a)/w, 2.) ) );
00102 }
00103 else
00104 {
00105 s = 1.;
00106 }
00107
00108
00109 realnum ionchar = ionic_charge/10.;
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120 if( iso_sequence == 3 )
00121 {
00122 E_c = 1.96274 + ionchar*(20.30014 + ionchar*(-0.97103 + ionchar*( 0.85453 + ionchar*( 0.13547 + 0.02401*ionchar))));
00123 }
00124 else if( iso_sequence == 4 )
00125 {
00126 E_c = 5.78908 + ionchar*(34.08270 + ionchar*( 1.51729 + ionchar*(-1.21227 + ionchar*( 0.77559 - 0.00410*ionchar))));
00127 }
00128 else if( iso_sequence == 7 )
00129 {
00130 E_c = 11.37092 + ionchar*(36.22053 + ionchar*( 7.08448 + ionchar*(-5.16840 + ionchar*( 2.45056 - 0.16961*ionchar))));
00131 }
00132 else if( iso_sequence == 11 )
00133 {
00134 E_c = 2.24809 + ionchar*(22.27768 + ionchar*(-1.12285 + ionchar*( 0.90267 + ionchar*(-0.03860 + 0.01468*ionchar))));
00135 }
00136 else if( iso_sequence == 12 )
00137 {
00138 E_c = 2.74508 + ionchar*(19.18623 + ionchar*(-0.54317 + ionchar*( 0.78685 + ionchar*(-0.04249 + 0.01357*ionchar))));
00139 }
00140 else if( iso_sequence == 15 )
00141 {
00142 E_c = 1.42762 + ionchar*( 3.90778 + ionchar*( 0.73119 + ionchar*(-1.91404 + ionchar*( 1.05059 - 0.08992*ionchar))));
00143 }
00144 else if( iso_sequence == 1 || iso_sequence == 2 || iso_sequence == 10 )
00145 {
00146
00147
00148 E_c = 1.e10;
00149 }
00150 else
00151 {
00152
00153 E_c = 0.0;
00154
00155
00156 }
00157
00158
00159
00160
00161
00162 snew = 1. + (s-1.)*exp(-(E_c*EVDEGK)/(c*T));
00163
00164 ASSERT( snew >=0. && snew <= 1. );
00165 return snew;
00166 }
00167
00181 STATIC double Badnell_DR_rate_eval(
00182
00183 int nAtomicNumberCScale,
00184
00185 int n_core_e_before_recomb )
00186 {
00187
00188 double RateCoefficient, sum;
00189 int i;
00190
00191 DEBUG_ENTRY( "Badnell_DR_rate_eval()" );
00192 ASSERT( nAtomicNumberCScale>=0 && nAtomicNumberCScale<LIMELM );
00193
00194 if( nAtomicNumberCScale==ipIRON && n_core_e_before_recomb>=12 &&
00195 n_core_e_before_recomb<=18 )
00196 {
00197
00198
00199
00200
00201
00202
00203
00204 double cFe_q[7][8] =
00205 {
00206 {5.636e-4, 7.390e-3, 3.635e-2, 1.693e-1, 3.315e-2, 2.288e-1, 7.316e-2, 0.},
00207 {1.090e-3, 7.801e-3, 1.132e-2, 4.740e-2, 1.990e-1, 3.379e-2, 1.140e-1, 1.250e-1},
00208 {3.266e-3, 7.637e-3, 1.005e-2, 2.527e-2, 6.389e-2, 1.564e-1, 0., 0.},
00209 {1.074e-3, 6.080e-3, 1.887e-2, 2.540e-2, 7.580e-2, 2.773e-1, 0., 0.},
00210 {9.073e-4, 3.777e-3, 1.027e-2, 3.321e-2, 8.529e-2, 2.778e-1, 0., 0.},
00211 {5.335e-4, 1.827e-3, 4.851e-3, 2.710e-2, 8.226e-2, 3.147e-1, 0., 0.},
00212 {7.421e-4, 2.526e-3, 4.605e-3, 1.489e-2, 5.891e-2, 2.318e-1, 0., 0.}
00213 };
00214
00215
00216 double EFe_q[7][8] =
00217 {
00218 {3.628e3, 2.432e4, 1.226e5, 4.351e5, 1.411e6, 6.589e6, 1.030e7, 0},
00219 {1.246e3, 1.063e4, 4.719e4, 1.952e5, 5.637e5, 2.248e6, 7.202e6, 3.999e9},
00220 {1.242e3, 1.001e4, 4.466e4, 1.497e5, 3.919e5, 6.853e5, 0. , 0.},
00221 {1.387e3, 1.048e4, 3.955e4, 1.461e5, 4.010e5, 7.208e5, 0. , 0.},
00222 {1.525e3, 1.071e4, 4.033e4, 1.564e5, 4.196e5, 7.580e5, 0. , 0.},
00223 {2.032e3, 1.018e4, 4.638e4, 1.698e5, 4.499e5, 7.880e5, 0. , 0.},
00224 {3.468e3, 1.353e4, 3.690e4, 1.957e5, 4.630e5, 8.202e5, 0. , 0.}
00225 };
00226
00227 long int nion = n_core_e_before_recomb - 12;
00228 ASSERT( nion>=0 && nion <=6 );
00229
00230 sum = 0;
00231 i = 0;
00232
00233 for(i=0; i<8; ++i )
00234 {
00235 sum += (cFe_q[nion][i] * sexp( EFe_q[nion][i]/phycon.te));
00236 }
00237
00238
00239 RateCoefficient = sum / phycon.te32;
00240 strcpy(chDRDataSource[nAtomicNumberCScale][nAtomicNumberCScale-n_core_e_before_recomb] ,
00241 "Bad06D");
00242
00243 return RateCoefficient;
00244 }
00245
00246
00247 else if( nAtomicNumberCScale < n_core_e_before_recomb )
00248 {
00249 RateCoefficient = -2;
00250 }
00251
00252 else if( nAtomicNumberCScale >= LIMELM )
00253 {
00254 RateCoefficient = -2;
00255 }
00256
00257 else if( !lgDRBadnellDefined[nAtomicNumberCScale][n_core_e_before_recomb] )
00258 {
00259 RateCoefficient = -1;
00260 }
00261 else if( lgDRBadnellDefined[nAtomicNumberCScale][n_core_e_before_recomb] )
00262 {
00263
00264 sum = 0;
00265 i = 0;
00266
00267 for(i=0; i<nDRFitPar[nAtomicNumberCScale][n_core_e_before_recomb]; ++i )
00268 {
00269 sum += (DRFitParPart1[nAtomicNumberCScale][n_core_e_before_recomb][i] *
00270 sexp( DRFitParPart2[nAtomicNumberCScale][n_core_e_before_recomb][i]/phycon.te));
00271 }
00272
00273 strcpy(chDRDataSource[nAtomicNumberCScale][nAtomicNumberCScale-n_core_e_before_recomb] ,
00274 "BadWeb");
00275
00276
00277 RateCoefficient = sum / phycon.te32;
00278 }
00279
00280 else
00281 {
00282 RateCoefficient = -99;
00283 }
00284
00285 ASSERT( RateCoefficient < 1e-6 );
00286
00287 return RateCoefficient;
00288 }
00289
00294 STATIC double Badnell_RR_rate_eval(
00295
00296 int nAtomicNumberCScale,
00297
00298 int n_core_e_before_recomb )
00299 {
00300 double RateCoefficient;
00301 double B, D, F;
00302
00303 DEBUG_ENTRY( "Badnell_RR_rate_eval()" );
00304
00305 ASSERT( nAtomicNumberCScale>=0 && nAtomicNumberCScale<LIMELM );
00306
00307 if( nAtomicNumberCScale==ipIRON &&
00308 n_core_e_before_recomb>=12 && n_core_e_before_recomb<=18 )
00309 {
00310
00311
00312
00313
00314
00315
00316
00317 double parFeq[7][6] ={
00318 {1.179e-9 , 0.7096, 4.508e2, 3.393e7, 0.0154, 3.977e6},
00319 {1.050e-9 , 0.6939, 4.568e2, 3.987e7, 0.0066, 5.451e5},
00320 {9.832e-10, 0.7146, 3.597e2, 3.808e7, 0.0045, 3.952e5},
00321 {8.303e-10, 0.7156, 3.531e2, 3.554e7, 0.0132, 2.951e5},
00322 {1.052e-9 , 0.7370, 1.639e2, 2.924e7, 0.0224, 4.291e5},
00323 {1.338e-9 , 0.7495, 7.242e1, 2.453e7, 0.0404, 4.199e5},
00324 {1.263e-9 , 0.7532, 5.209e1, 2.169e7, 0.0421, 2.917e5}
00325 };
00326
00327 double temp;
00328
00329 long int nion = n_core_e_before_recomb - 12;
00330 ASSERT( nion>=0 && nion <=6 );
00331
00332 temp = -parFeq[nion][5]/phycon.te;
00333 B = parFeq[nion][1] + parFeq[nion][4]*exp(temp);
00334 D = sqrt(phycon.te/parFeq[nion][2]);
00335 F = sqrt(phycon.te/parFeq[nion][3]);
00336 RateCoefficient = parFeq[nion][0]/(D*pow((1.+D),(1.-B))*pow((1.+F),(1.+B)));
00337 strcpy(chRRDataSource[nAtomicNumberCScale][nAtomicNumberCScale-n_core_e_before_recomb] ,"Bad06");
00338
00339 return RateCoefficient;
00340 }
00341
00342
00343 else if( nAtomicNumberCScale < n_core_e_before_recomb )
00344 {
00345 RateCoefficient = -2;
00346 }
00347
00348 else if( nAtomicNumberCScale >= LIMELM )
00349 {
00350 RateCoefficient = -2;
00351 }
00352
00353 else if( !lgRRBadnellDefined[nAtomicNumberCScale][n_core_e_before_recomb] )
00354 {
00355 RateCoefficient = -1;
00356 }
00357
00358 else if( lgRRBadnellDefined[nAtomicNumberCScale][n_core_e_before_recomb] )
00359 {
00360
00361
00362
00363 double temp;
00364
00365 temp = -RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][5]/phycon.te;
00366 B = RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][1] +
00367 RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][4]*exp(temp);
00368 D = sqrt(phycon.te/RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][2]);
00369 F = sqrt(phycon.te/RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][3]);
00370 RateCoefficient = RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][0]/(D*pow((1.+D),(1.-B))*pow((1.+F),(1.+B)));
00371 strcpy(chRRDataSource[nAtomicNumberCScale][nAtomicNumberCScale-n_core_e_before_recomb] ,"Bad06");
00372 }
00373
00374
00375 else
00376 RateCoefficient = -99;
00377
00378 return RateCoefficient;
00379 }
00380
00381
00382
00383
00384
00385 void Badnell_rec_init( void )
00386 {
00387
00388 double par_C[MAX_FIT_PAR_DR];
00389 double par_E[MAX_FIT_PAR_DR];
00390 char chLine[INPUT_LINE_LENGTH];
00391 int NuclearCharge=-1, NumberElectrons=-1;
00392 int count, number;
00393 double temp_par[MAX_FIT_PAR_RR];
00394 int M_state, W_state;
00395
00396 const int NBLOCK = 2;
00397 int data_begin_line[NBLOCK];
00398 int length_of_line;
00399 FILE *ioDATA;
00400 const char* chFilename;
00401 int yr, mo, dy;
00402 char *chs;
00403
00404 const int BIGGEST_INDEX_TO_USE = 103;
00405
00406
00407 long TheirIndexToOurIndex[BIGGEST_INDEX_TO_USE];
00408 char string[120];
00409 double value;
00410 bool lgEOL;
00411 long int i1;
00412 long INDX=0,INDP=0,N=0,S=0,L=0,J=0,maxINDX=0,loopindex=0,max_N_of_data=-1;
00413 bool lgFlag = true;
00414
00415 static int nCalled = 0;
00416
00417 const char* cdDATAFILE[] =
00418 {
00419
00420 "",
00421 "UTA/nrb00_h_he1ic12.dat",
00422 "UTA/nrb00_h_li2ic12.dat",
00423 "UTA/nrb00_h_be3ic12.dat",
00424 "UTA/nrb00_h_b4ic12.dat",
00425 "UTA/nrb00_h_c5ic12.dat",
00426 "UTA/nrb00_h_n6ic12.dat",
00427 "UTA/nrb00_h_o7ic12.dat",
00428 "UTA/nrb00_h_f8ic12.dat",
00429 "UTA/nrb00_h_ne9ic12.dat",
00430 "UTA/nrb00_h_na10ic12.dat",
00431 "UTA/nrb00_h_mg11ic12.dat",
00432 "UTA/nrb00_h_al12ic12.dat",
00433 "UTA/nrb00_h_si13ic12.dat",
00434 "UTA/nrb00_h_p14ic12.dat",
00435 "UTA/nrb00_h_s15ic12.dat",
00436 "UTA/nrb00_h_cl16ic12.dat",
00437 "UTA/nrb00_h_ar17ic12.dat",
00438 "UTA/nrb00_h_k18ic12.dat",
00439 "UTA/nrb00_h_ca19ic12.dat",
00440 "UTA/nrb00_h_sc20ic12.dat",
00441 "UTA/nrb00_h_ti21ic12.dat",
00442 "UTA/nrb00_h_v22ic12.dat",
00443 "UTA/nrb00_h_cr23ic12.dat",
00444 "UTA/nrb00_h_mn24ic12.dat",
00445 "UTA/nrb00_h_fe25ic12.dat",
00446 "UTA/nrb00_h_co26ic12.dat",
00447 "UTA/nrb00_h_ni27ic12.dat",
00448 "UTA/nrb00_h_cu28ic12.dat",
00449 "UTA/nrb00_h_zn29ic12.dat"
00450 };
00451
00452
00453 DEBUG_ENTRY( "Badnell_rec_init()" );
00454
00455
00456 if( nCalled > 0 )
00457 {
00458 return;
00459 }
00460 ++nCalled;
00461
00462 # if defined(PRINT_DR) || defined(PRINT_RR)
00463 FILE *ofp = open_data( FILE_NAME_OUT, "w", AS_LOCAL_ONLY );
00464 # endif
00465
00466 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00467 {
00468 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
00469 {
00470 if( nelem < 2 || dense.lgElmtOn[nelem] )
00471 {
00472 for( long ipHi=0; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
00473 {
00474 for( long k=0; k<NUM_DR_TEMPS; ++k )
00475 iso_sp[ipISO][nelem].fb[ipHi].DielecRecombVsTemp[k] = 0.;
00476 }
00477 }
00478 }
00479 }
00480
00481
00482
00483 for( long nelem=ipHELIUM; nelem<LIMELM; nelem++)
00484 {
00485 if( nelem < 2 || dense.lgElmtOn[nelem] )
00486 {
00487 ioDATA= open_data( cdDATAFILE[nelem], "r" );
00488
00489 lgFlag = true;
00490 ASSERT(ioDATA);
00491
00492 for( long i=0; i<BIGGEST_INDEX_TO_USE; i++ )
00493 TheirIndexToOurIndex[i] = -1;
00494
00495
00496 while(lgFlag)
00497 {
00498 if(read_whole_line(string,sizeof(string),ioDATA)!=NULL)
00499 {
00500 if( nMatch("INDX INDP ",string) )
00501 {
00502
00503 if( read_whole_line( string , (int)sizeof(string) , ioDATA ) == NULL )
00504 {
00505 fprintf( ioQQQ, " Badnell data file appears to be corrupted.\n");
00506 cdEXIT(EXIT_FAILURE);
00507 }
00508
00509
00510 while( read_whole_line(string, (int)sizeof(string), ioDATA) != NULL )
00511 {
00512 if( strcmp(string,"\n")==0 )
00513 {
00514 lgFlag = false;
00515 break;
00516 }
00517
00518 i1=3;
00519 INDX=(long)FFmtRead(string,&i1,sizeof(string),&lgEOL);
00520 if( INDX >= BIGGEST_INDEX_TO_USE )
00521 {
00522 INDX--;
00523 lgFlag = false;
00524 break;
00525 }
00526
00527 ASSERT( INDX < BIGGEST_INDEX_TO_USE );
00528
00529 INDP=(long)FFmtRead(string,&i1,sizeof(string),&lgEOL);
00530 ASSERT( INDP >= 1 );
00531
00532 if(INDP==1)
00533 {
00534 if( (i1=nMatch("1S1 ",string)) > 0 )
00535 {
00536 i1 += 4;
00537 N=(long)FFmtRead(string,&i1,sizeof(string),&lgEOL);
00538 ASSERT( N>=1 );
00539 }
00540 else
00541 {
00542 TotalInsanity();
00543 }
00544
00545 if( (i1=nMatch(" (",string)) > 0 )
00546 {
00547 i1 += 6;
00548 S=(long)FFmtRead(string,&i1,sizeof(string),&lgEOL);
00549
00550 ASSERT( S==1 || S==3 );
00551 }
00552 else
00553 {
00554 TotalInsanity();
00555 }
00556
00557
00558 i1++;
00559 L=(long)FFmtRead(string,&i1,sizeof(string),&lgEOL);
00560 ASSERT( L >= 0 && L < N );
00561
00562
00563 i1 += 2;
00564 J=(long)FFmtRead(string,&i1,sizeof(string),&lgEOL);
00565 ASSERT( J <= ( L + (int)((S+1)/2) ) &&
00566 J >= ( L - (int)((S+1)/2) ) && J >= 0 );
00567
00568
00569 if( N<= iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max + iso_sp[ipHE_LIKE][nelem].nCollapsed_max )
00570 TheirIndexToOurIndex[INDX] = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[N][L][S];
00571 else
00572 {
00573
00574
00575 INDX--;
00576 lgFlag = false;
00577 break;
00578 }
00579
00580
00581 if( N==2 && L==1 && S==3 )
00582 {
00583 if( J==0 )
00584 TheirIndexToOurIndex[INDX] = 3;
00585 else if( J==1 )
00586 TheirIndexToOurIndex[INDX] = 4;
00587 else
00588 {
00589 ASSERT( J==2 );
00590 ASSERT( TheirIndexToOurIndex[INDX] == 5 );
00591 }
00592 }
00593 max_N_of_data = MAX2( max_N_of_data, N );
00594 }
00595 else
00596 {
00597
00598 lgFlag = false;
00599 }
00600 }
00601 }
00602 }
00603 else
00604 {
00605
00606 lgFlag = false;
00607 }
00608 }
00609
00610 maxINDX =INDX;
00611 ASSERT( maxINDX > 0 );
00612 ASSERT( maxINDX < BIGGEST_INDEX_TO_USE );
00613
00614 INDX = 0;
00615 lgFlag = true;
00616 while(lgFlag)
00617 {
00618 if(read_whole_line(string,sizeof(string),ioDATA)!=NULL)
00619 {
00620
00621 if( nMatch("INDX TE= ",string) )
00622 {
00623 lgFlag = false;
00624
00625
00626 if( read_whole_line( string , (int)sizeof(string) , ioDATA ) == NULL )
00627 {
00628 fprintf( ioQQQ, " Badnell data file appears to be corrupted.\n");
00629 cdEXIT(EXIT_FAILURE);
00630 }
00631
00632
00633 while( read_whole_line(string, (int)sizeof(string), ioDATA) != NULL )
00634 {
00635
00636 if( nMatch("PRTF",string) || INDX >= maxINDX || INDX<0 )
00637 break;
00638
00639 i1=3;
00640 INDX=(long)FFmtRead(string,&i1,sizeof(string),&lgEOL);
00641 if( INDX>maxINDX )
00642 break;
00643
00644 freeBound *fb;
00645
00646 if( TheirIndexToOurIndex[INDX] < iso_sp[ipHE_LIKE][nelem].numLevels_max &&
00647 TheirIndexToOurIndex[INDX] > 0 )
00648 fb = &iso_sp[ipHE_LIKE][nelem].fb[TheirIndexToOurIndex[INDX]];
00649 else
00650 continue;
00651
00652 for(loopindex=0;loopindex<10;loopindex++)
00653 {
00654 value=(double)FFmtRead(string,&i1,sizeof(string),&lgEOL);
00655 fb->DielecRecombVsTemp[loopindex] += value;
00656 }
00657
00658
00659 if( read_whole_line( string , (int)sizeof(string) , ioDATA ) == NULL )
00660 {
00661 fprintf( ioQQQ, " Badnell data file appears to be corrupted.\n");
00662 cdEXIT(EXIT_FAILURE);
00663 }
00664
00665
00666 i1 = 13;
00667 for(loopindex=10;loopindex<19;loopindex++)
00668 {
00669 value=(double)FFmtRead(string,&i1,sizeof(string),&lgEOL);
00670 fb->DielecRecombVsTemp[loopindex] += value;
00671 }
00672 }
00673 }
00674 }
00675 else
00676 lgFlag = false;
00677 }
00678 fclose(ioDATA);
00679 ASSERT( maxINDX > 0 );
00680 ASSERT( maxINDX < BIGGEST_INDEX_TO_USE );
00681 ASSERT( max_N_of_data > 0 );
00682
00683 if( max_N_of_data < iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max + iso_sp[ipHE_LIKE][nelem].nCollapsed_max )
00684 {
00685 long indexOfMaxN;
00686 L = -1;
00687 S = -1;
00688
00689
00690 for( long i=TheirIndexToOurIndex[maxINDX]+1;
00691 i<iso_sp[ipHE_LIKE][nelem].numLevels_max-iso_sp[ipHE_LIKE][nelem].nCollapsed_max; i++ )
00692 {
00693 long ipISO = ipHE_LIKE;
00694 L = L_(i);
00695 S = S_(i);
00696
00697 if( L > 4 )
00698 continue;
00699
00700 indexOfMaxN = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[max_N_of_data][L][S];
00701 for(loopindex=0;loopindex<19;loopindex++)
00702 {
00703 iso_sp[ipHE_LIKE][nelem].fb[i].DielecRecombVsTemp[loopindex] =
00704 iso_sp[ipHE_LIKE][nelem].fb[indexOfMaxN].DielecRecombVsTemp[loopindex] *
00705 pow3( (double)max_N_of_data/(double)iso_sp[ipHE_LIKE][nelem].st[i].n());
00706 }
00707 }
00708
00709
00710 indexOfMaxN =
00711 iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max][1][1];
00712
00713
00714 for( long i=iso_sp[ipHE_LIKE][nelem].numLevels_max-iso_sp[ipHE_LIKE][nelem].nCollapsed_max;
00715 i<iso_sp[ipHE_LIKE][nelem].numLevels_max; i++ )
00716 {
00717 for(loopindex=0;loopindex<19;loopindex++)
00718 {
00719 iso_sp[ipHE_LIKE][nelem].fb[i].DielecRecombVsTemp[loopindex] =
00720 iso_sp[ipHE_LIKE][nelem].fb[indexOfMaxN].DielecRecombVsTemp[loopindex] *
00721 pow3( (double)iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max/
00722 (double)iso_sp[ipHE_LIKE][nelem].st[i].n());
00723 }
00724 }
00725 }
00726 }
00727 }
00728
00729 for( long i=0; i<NBLOCK; ++i )
00730 {
00731
00732 data_begin_line[i] = INT_MIN;
00733 }
00734
00735 chFilename = "badnell_dr.dat";
00736 ioDATA = open_data( chFilename, "r" );
00737
00738 count = 0;
00739 number = 0;
00740
00741
00742
00743 while( read_whole_line(chLine, (int)sizeof(chLine), ioDATA) != NULL )
00744 {
00745 count++;
00746
00747 if( chLine[2]=='Z' )
00748 {
00749
00750
00751 data_begin_line[number] = count;
00752 ASSERT( number < NBLOCK );
00753 number++;
00754 }
00755 }
00756
00757
00758 if( lgMustMallocRec )
00759 {
00760 nDRFitPar = (int**)MALLOC( LIMELM*sizeof( int*) );
00761 lgDRBadnellDefined = (bool **)MALLOC( LIMELM*sizeof(bool*) );
00762 lgDR_BadWeb_exist = (bool **)MALLOC( LIMELM*sizeof(bool*) );
00763 lgDRBadnellDefinedPart2 = (bool **)MALLOC( LIMELM*sizeof(bool*) );
00764 lgRRBadnellDefined = (bool **)MALLOC( LIMELM*sizeof(bool*) );
00765
00766 DRFitParPart1 = (double ***)MALLOC( LIMELM*sizeof(double**) );
00767 DRFitParPart2 = (double ***)MALLOC( LIMELM*sizeof(double**) );
00768 RRFitPar = (double ***)MALLOC( LIMELM*sizeof(double**) );
00769 }
00770
00771 for( long nelem=0; nelem<LIMELM; nelem++ )
00772 {
00773 if( lgMustMallocRec )
00774 {
00775 nDRFitPar[nelem] = (int*)MALLOC( (nelem+1)*sizeof( int) );
00776 lgDR_BadWeb_exist[nelem] = (bool *)MALLOC( (nelem+1)*sizeof(bool) );
00777 lgDRBadnellDefined[nelem] = (bool *)MALLOC( (nelem+1)*sizeof(bool) );
00778 lgDRBadnellDefinedPart2[nelem] = (bool *)MALLOC( (nelem+1)*sizeof(bool) );
00779 lgRRBadnellDefined[nelem] = (bool *)MALLOC( (nelem+1)*sizeof(bool) );
00780
00781 DRFitParPart1[nelem] = (double **)MALLOC( (nelem+1)*sizeof(double*) );
00782 DRFitParPart2[nelem] = (double **)MALLOC( (nelem+1)*sizeof(double*) );
00783 RRFitPar[nelem] = (double **)MALLOC( (nelem+1)*sizeof(double*) );
00784 }
00785 for( long ion=0; ion<nelem+1; ++ion )
00786 {
00787 if( lgMustMallocRec )
00788 {
00789 DRFitParPart1[nelem][ion] = (double *)MALLOC( MAX_FIT_PAR_DR*sizeof(double) );
00790 DRFitParPart2[nelem][ion] = (double *)MALLOC( MAX_FIT_PAR_DR*sizeof(double) );
00791 RRFitPar[nelem][ion] = (double *)MALLOC( MAX_FIT_PAR_RR*sizeof(double) );
00792 }
00793 lgDRBadnellDefined[nelem][ion] = false;
00794 lgDRBadnellDefinedPart2[nelem][ion] = false;
00795 lgRRBadnellDefined[nelem][ion] = false;
00796
00797
00798 for( long k=0; k<MAX_FIT_PAR_DR; k++ )
00799 {
00800 DRFitParPart1[nelem][ion][k] = 0;
00801 DRFitParPart2[nelem][ion][k] = 0;
00802 }
00803 for( long k=0; k<MAX_FIT_PAR_RR; k++ )
00804 {
00805 RRFitPar[nelem][ion][k] = 0;
00806 }
00807 }
00808 }
00809 lgMustMallocRec = false;
00810
00811 count = 0;
00812
00813 fseek(ioDATA, 0, SEEK_SET);
00814
00815 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00816 {
00817 fprintf( ioQQQ, " DISASTER PROBLEM Badnell_rec_init could not read first line of badnell_dr.dat.\n");
00818 cdEXIT(EXIT_FAILURE);
00819 }
00820 count++;
00821
00822
00823 if( (chs = strchr_s(chLine, ')'))==NULL )
00824 {
00825
00826 fprintf( ioQQQ, " DISASTER PROBLEM Badnell_rec_init data file incorrect format.\n");
00827 cdEXIT(EXIT_FAILURE);
00828 }
00829
00830 ++chs;
00831 sscanf(chs, "%4i%2i%2i",&yr, &mo, &dy);
00832
00833 int dr_yr = 2012, dr_mo = 6, dr_dy = 28;
00834 if((yr != dr_yr) || (mo != dr_mo) || (dy != dr_dy))
00835 {
00836 fprintf(ioQQQ,
00837 "DISASTER PROBLEM Badnell_rec_init The version of %s I found (%i %i %i) is not the current version (%i %i %i).\n",
00838 chFilename, yr, mo, dy, dr_yr, dr_mo, dr_dy);
00839 fprintf(ioQQQ," The first line of the file is the following\n %s\n", chLine );
00840 cdEXIT(EXIT_FAILURE);
00841 }
00842
00843 while( read_whole_line(chLine, (int)sizeof(chLine), ioDATA) != NULL )
00844 {
00845 count++;
00846 length_of_line = (int)strlen(chLine);
00847
00848
00849 if( count > data_begin_line[0] && count < data_begin_line[1] && length_of_line >3 )
00850 {
00851
00852 for( long i=0; i<MAX_FIT_PAR_DR; i++ )
00853 {
00854 par_C[i] = 0;
00855 }
00856 sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf%lf%lf%lf",
00857 &NuclearCharge, &NumberElectrons, &M_state, &W_state, &par_C[0], &par_C[1], &par_C[2],
00858 &par_C[3], &par_C[4], &par_C[5], &par_C[6], &par_C[7], &par_C[8]);
00859
00860
00861 long int NuclearChargeM1 = NuclearCharge-1;
00862
00863 if(M_state == 1 && NuclearChargeM1 < LIMELM )
00864 {
00865
00866 ASSERT( NumberElectrons < LIMELM );
00867 ASSERT( NuclearChargeM1 < LIMELM );
00868 lgDRBadnellDefined[NuclearChargeM1][NumberElectrons] = true;
00869
00870
00871 nDRFitPar[NuclearChargeM1][NumberElectrons] = 9;
00872 for( long i=8; i>=0; i-- )
00873 {
00874 if( par_C[i] == 0 )
00875 --nDRFitPar[NuclearChargeM1][NumberElectrons];
00876 else
00877 break;
00878 }
00879
00880
00881 for( long i=0; i<9; i++ )
00882 DRFitParPart1[NuclearChargeM1][NumberElectrons][i] = par_C[i];
00883 }
00884 }
00885 }
00886
00887
00888 fseek(ioDATA, 0, SEEK_SET);
00889 count = 0;
00890 while( read_whole_line(chLine, (int)sizeof(chLine), ioDATA) != NULL )
00891 {
00892 count++;
00893 length_of_line = (int)strlen(chLine);
00894 if( count > data_begin_line[1] && length_of_line > 3 )
00895 {
00896
00897
00898 for( long i=0; i<MAX_FIT_PAR_DR; i++ )
00899 {
00900 par_E[i] = 0;
00901 }
00902 sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf%lf%lf%lf",
00903 &NuclearCharge, &NumberElectrons, &M_state, &W_state, &par_E[0], &par_E[1], &par_E[2],
00904 &par_E[3], &par_E[4], &par_E[5], &par_E[6], &par_E[7], &par_E[8]);
00905
00906 long int NuclearChargeM1 = NuclearCharge-1;
00907
00908 if(M_state == 1 && NuclearChargeM1<LIMELM)
00909 {
00910 ASSERT( NumberElectrons < LIMELM );
00911 ASSERT( NuclearChargeM1 < LIMELM );
00912 lgDRBadnellDefinedPart2[NuclearChargeM1][NumberElectrons] = true;
00913
00914
00915 nDRFitPar[NuclearChargeM1][NumberElectrons] = 9;
00916 for( long i=8; i>=0; i-- )
00917 {
00918 if( par_E[i] == 0 )
00919 --nDRFitPar[NuclearChargeM1][NumberElectrons];
00920 else
00921 break;
00922 }
00923
00924
00925 for( long i=0; i<nDRFitPar[NuclearChargeM1][NumberElectrons]; i++ )
00926 DRFitParPart2[NuclearChargeM1][NumberElectrons][i] = par_E[i];
00927 }
00928 }
00929 }
00930
00931 fclose( ioDATA );
00932
00933
00934 # ifdef PRINT_DR
00935 for( long nelem=0; nelem<LIMELM; nelem++ )
00936 {
00937 for( int ion=0; ion<nelem+1;++ion )
00938 {
00939 if( lgDRBadnellDefined[nelem][ion] )
00940 {
00941 fprintf(ofp, "%i %i %e %e %e %e %e %e %e %e %e\n",
00942 nelem, ion, DRFitParPart1[nelem][ion][0],
00943 DRFitParPart1[nelem][ion][1], DRFitParPart1[nelem][ion][2],
00944 DRFitParPart1[nelem][ion][3], DRFitParPart1[nelem][ion][4],
00945 DRFitParPart1[nelem][ion][5], DRFitParPart1[nelem][ion][6],
00946 DRFitParPart1[nelem][ion][7], DRFitParPart1[nelem][ion][8]);
00947 }
00948 }
00949 }
00950 for( long nelem=0; nelem<LIMELM; nelem++ )
00951 {
00952 for( int ion=0; ion<nelem+1; ion++ )
00953 {
00954 if( lgDRBadnellDefinedPart2[nelem][ion] )
00955 {
00956 fprintf(ofp, "%i %i %e %e %e %e %e %e %e %e %e\n",
00957 nelem, ion, DRFitParPart2[nelem][ion][0],
00958 DRFitParPart2[nelem][ion][1], DRFitParPart2[nelem][ion][2],
00959 DRFitParPart2[nelem][ion][3], DRFitParPart2[nelem][ion][4],
00960 DRFitParPart2[nelem][ion][5], DRFitParPart2[nelem][ion][6],
00961 DRFitParPart2[nelem][ion][7], DRFitParPart2[nelem][ion][8]);
00962 }
00963 }
00964 }
00965 fclose(ofp);
00966 # endif
00967
00968
00969
00970 bool lgDRBadnellBothDefined = true;
00971 for( int nelem=0; nelem<LIMELM; nelem++ )
00972 {
00973 for( int ion=0; ion<nelem+1; ion++ )
00974 {
00975
00976
00977 if( lgDRBadnellDefined[nelem][ion] != lgDRBadnellDefinedPart2[nelem][ion] )
00978 {
00979 fprintf( ioQQQ, "DR %i, RR %i: %c %c\n", nelem, ion,
00980 TorF(lgDRBadnellDefined[nelem][ion]),
00981 TorF(lgDRBadnellDefinedPart2[nelem][ion]));
00982 fprintf( ioQQQ, "PROBLEM ion_recomb_Badnell first and second half of Badnell DR not consistent.\n");
00983 lgDRBadnellBothDefined = false;
00984 }
00985 }
00986 }
00987
00988 if( !lgDRBadnellBothDefined )
00989 {
00990
00991 fprintf(ioQQQ,
00992 "DISASTER PROBLEM The DR data files are corrupted - part 1 and 2 do not agree.\n");
00993 fprintf(ioQQQ," Start again with a fresh copy of the data directory\n" );
00994 cdEXIT(EXIT_FAILURE);
00995 }
00996
00997
00998 chFilename = "badnell_rr.dat";
00999 ioDATA = open_data( chFilename, "r" );
01000
01001
01002 {
01003 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01004 {
01005 fprintf( ioQQQ, " DISASTER PROBLEM Badnell_rec_init could not read first line of badnell_rr.dat.\n");
01006 cdEXIT(EXIT_FAILURE);
01007 }
01008
01009 if( (chs = strchr_s(chLine, ')'))==NULL )
01010 {
01011
01012 fprintf( ioQQQ, " DISASTER PROBLEM Badnell_rec_init data file incorrect format.\n");
01013 cdEXIT(EXIT_FAILURE);
01014 }
01015 ++chs;
01016 sscanf(chs, "%4i%2i%2i", &yr, &mo, &dy);
01017 int rr_yr = 2011, rr_mo = 4, rr_dy = 12;
01018 if((yr != rr_yr)||(mo != rr_mo)||(dy != rr_dy))
01019 {
01020 fprintf(ioQQQ,"DISASTER PROBLEM The version of %s I found (%i %i %i) is not the current version (%i %i %i).\n",
01021 chFilename, yr, mo, dy, rr_yr, rr_mo, rr_dy);
01022 fprintf(ioQQQ," The line was as follows:\n %s\n", chLine );
01023 cdEXIT(EXIT_FAILURE);
01024 }
01025 }
01026
01027 while( read_whole_line(chLine, (int)sizeof(chLine), ioDATA) != NULL )
01028 {
01029
01030 for( long i=0; i<MAX_FIT_PAR_RR; i++ )
01031 {
01032 temp_par[i] = 0;
01033 }
01034 if(chLine[0] != '#')
01035 {
01036 sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf",
01037 &NuclearCharge, &NumberElectrons, &M_state, &W_state, &temp_par[0], &temp_par[1],
01038 &temp_par[2], &temp_par[3], &temp_par[4], &temp_par[5]);
01039 long NuclearChargeM1 = NuclearCharge-1;
01040
01041 if(M_state == 1 && NuclearChargeM1<LIMELM)
01042 {
01043 ASSERT( NuclearChargeM1 < LIMELM );
01044 ASSERT( NumberElectrons <= LIMELM );
01045
01046 lgRRBadnellDefined[NuclearChargeM1][NumberElectrons] = true;
01047
01048 for( long i=0; i<MAX_FIT_PAR_RR; i++ )
01049 RRFitPar[NuclearChargeM1][NumberElectrons][i] = temp_par[i];
01050 }
01051 }
01052 }
01053
01054
01055 # ifdef PRINT_RR
01056 count = 0;
01057 for( long nelem=0; nelem<LIMELM; nelem++ )
01058 {
01059 for( long ion=0; ion<nelem+1; ion++ )
01060 {
01061 if( lgRRBadnellDefined[nelem][ion] )
01062 {
01063 fprintf(ofp, "%li %li %e %e %e %e %e %e\n",
01064 nelem, ion, RRFitPar[nelem][ion][0],
01065 RRFitPar[nelem][ion][1], RRFitPar[nelem][ion][2],
01066 RRFitPar[nelem][ion][3],
01067 RRFitPar[nelem][ion][4], RRFitPar[nelem][ion][5]);
01068 count++;
01069 }
01070 }
01071 }
01072 fprintf(ofp, "total lines are %i ", count);
01073
01074 fclose(ofp);
01075 # endif
01076
01077 fclose(ioDATA);
01078
01079 {
01080 enum {DEBUG_LOC=false};
01081 if( DEBUG_LOC )
01082 {
01083 for( int nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
01084 {
01085 fprintf(ioQQQ,"\nDEBUG rr rec\t%i",nelem);
01086 for( int ion=0; ion<=nelem; ++ion )
01087 {
01088 fprintf(ioQQQ,"\t%.2e", Badnell_RR_rate_eval(nelem+1 , nelem-ion ) );
01089 }
01090 fprintf(ioQQQ,"\n");
01091 fprintf(ioQQQ,"DEBUG dr rec\t%i",nelem);
01092 for( int ion=0; ion<=nelem; ++ion )
01093 {
01094 fprintf(ioQQQ,"\t%.2e", Badnell_DR_rate_eval(nelem+1 , nelem-ion ) );
01095 }
01096 fprintf(ioQQQ,"\n");
01097 }
01098 cdEXIT(EXIT_FAILURE);
01099 }
01100 }
01101
01102
01103
01104 if( ionbal.guess_noise !=0. )
01105 {
01106 for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
01107
01108
01109 RecNoise[nelem] = pow(10., RandGauss( 0. , ionbal.guess_noise ) );
01110 }
01111 else
01112 {
01113 for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
01114 RecNoise[nelem] = 1.;
01115 }
01116
01117
01118 for( long nelem=0; nelem<LIMELM; ++nelem )
01119 DR_Badnell_rate_coef_mean_ion[nelem] = 0.;
01120
01121 return;
01122 }
01123
01124
01125 void ion_recom_calculate( void )
01126 {
01127 static double TeUsed = -1 , EdenUsed = -1.;
01128
01129 DEBUG_ENTRY( "ion_recom_calculate()" );
01130
01131
01132 if( fp_equal(phycon.te,TeUsed) && fp_equal( dense.eden , EdenUsed ))
01133 return;
01134
01135
01136
01137
01138 TeUsed = phycon.te;
01139 EdenUsed = dense.eden;
01140
01141 for( long nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
01142 {
01143
01144 for( long ion=0; ion < nelem+1; ++ion )
01145 {
01146 long int n_bnd_elec_before_recom ,
01147 n_bnd_elec_after_recom;
01148
01149 n_bnd_elec_before_recom = nelem-ion;
01150 n_bnd_elec_after_recom = nelem-ion+1;
01151
01152
01153 ionbal.DR_Badnell_rate_coef[nelem][ion] = -1.;
01154 ionbal.RR_rate_coef_used[nelem][ion] = 0.;
01155 strcpy( chDRDataSource[nelem][ion] , "none" );
01156 strcpy( chRRDataSource[nelem][ion] , "none" );
01157
01158
01159 if( (ionbal.DR_Badnell_rate_coef[nelem][ion] =
01160 Badnell_DR_rate_eval(
01161
01162 nelem,
01163
01164
01165 n_bnd_elec_before_recom )) >= 0. )
01166 {
01167 lgDR_BadWeb_exist[nelem][ion] = true;
01168 }
01169 else
01170 {
01171
01172 lgDR_BadWeb_exist[nelem][ion] = false;
01173 }
01174
01175
01176
01177 ionbal.RR_Verner_rate_coef[nelem][ion] =
01178 t_ADfA::Inst().rad_rec(
01179
01180 nelem+1 ,
01181
01182 n_bnd_elec_after_recom ,
01183 phycon.te );
01184
01185
01186 if( (ionbal.RR_Badnell_rate_coef[nelem][ion] = Badnell_RR_rate_eval(
01187
01188 nelem,
01189
01190 n_bnd_elec_before_recom )) >= 0. )
01191 {
01192 ionbal.RR_rate_coef_used[nelem][ion] = ionbal.RR_Badnell_rate_coef[nelem][ion];
01193 }
01194 else
01195 {
01196 strcpy( chRRDataSource[nelem][ion] , "Verner" );
01197 ionbal.RR_rate_coef_used[nelem][ion] = ionbal.RR_Verner_rate_coef[nelem][ion];
01198 }
01199 }
01200
01201 ionbal.DR_Badnell_rate_coef[nelem][nelem] = 0.;
01202 strcpy(chDRDataSource[nelem][nelem] , "NA");
01203 }
01204
01205
01206 double Fe_Gu_c[9][6] = {
01207 { 2.50507e-11, 5.60226e-11, 1.85001e-10, 3.57495e-9, 1.66321e-7, 0. },
01208 { 9.19610e-11, 2.92460e-10, 1.02120e-9, 1.14852e-8, 3.25418e-7, 0. },
01209 { 9.02625e-11, 6.22962e-10, 5.77545e-9, 1.78847e-8, 3.40610e-7, 0. },
01210 { 9.04286e-12, 9.68148e-10, 4.83636e-9, 2.48159e-8, 3.96815e-7, 0. },
01211 { 6.77873e-10, 1.47252e-9, 5.31202e-9, 2.54793e-8, 3.47407e-7, 0. },
01212 { 1.29742e-9, 4.10172e-9, 1.23605e-8, 2.33615e-8, 2.97261e-7, 0. },
01213 { 8.78027e-10, 2.31680e-9, 3.49333e-9, 1.16927e-8, 8.18537e-8, 1.54740e-7 },
01214 { 2.23178e-10, 1.87313e-9, 2.86171e-9, 1.38575e-8, 1.17803e-7, 1.06251e-7 },
01215 { 2.17263e-10, 7.35929e-10, 2.81276e-9, 1.32411e-8, 1.15761e-7, 4.80389e-8 }
01216 },
01217
01218 Fe_Gu_E[9][6] = {
01219 { 8.30501e-2, 8.52897e-1, 3.40225e0, 2.23053e1, 6.80367e1, 0. },
01220 { 1.44392e-1, 9.23999e-1, 5.45498e0, 2.04301e1, 7.06112e1, 0. },
01221 { 5.79132e-2, 1.27852e0, 3.22439e0, 1.79602e1, 6.96277e1, 0. },
01222 { 1.02421e-1, 1.79393e0, 4.83226e0, 1.91117e1, 6.80858e1, 0. },
01223 { 1.24630e-1, 6.86045e-1, 3.09611e0, 1.44023e1, 6.42820e1, 0. },
01224 { 1.34459e-1, 6.63028e-1, 2.61753e0, 1.30392e1, 6.10222e1, 0. },
01225 { 7.79748e-2, 5.35522e-1, 1.88407e0, 8.38459e0, 3.38613e1, 7.89706e1 },
01226 { 8.83019e-2, 6.12756e-1, 2.36035e0, 9.61736e0, 3.64467e1, 8.72406e1 },
01227 { 1.51322e-1, 5.63155e-1, 2.57013e0, 9.08166e0, 3.69528e1, 1.08067e2 }
01228 };
01229
01230
01231 double te_eV32 = sqrt(pow3(phycon.te_eV));
01232
01233
01234
01235
01236
01237 for( long ion=0; ion<9; ion++ )
01238 {
01239
01240 if( ionbal.DR_Badnell_rate_coef[ipIRON][ion+5]<0. )
01241 {
01242 double fitSum = 0;
01243 for( long i=0; i<6; i++ )
01244 {
01245 fitSum += Fe_Gu_c[ion][i] * sexp( Fe_Gu_E[ion][i]/phycon.te_eV );
01246 }
01247 strcpy(chDRDataSource[ipIRON][ion+5] , "GuPC");
01248 lgDR_BadWeb_exist[ipIRON][ion+5] = true;
01249 ionbal.DR_Badnell_rate_coef[ipIRON][ion+5] = fitSum / te_eV32;
01250 }
01251 }
01252
01253
01254
01255 double BadnelDR_RateSave[LIMELM] =
01256 {
01257 3.78e-13, 1.70e-12, 8.14e-12, 1.60e-11, 2.38e-11,
01258 6.42e-11, 5.97e-11, 1.47e-10, 1.11e-10, 3.26e-10,
01259 1.88e-10, 2.06e-10, 4.14e-10, 3.97e-10, 2.07e-10,
01260 2.46e-10, 3.38e-10, 3.15e-10, 9.70e-11, 6.49e-11,
01261 6.93e-10, 3.70e-10, 3.29e-11, 4.96e-11, 5.03e-11,
01262 2.91e-12, 4.62e-14, 0.00e+00, 0.00e+00, 0.00e+00
01263 };
01264 for( long nelem=0; nelem < LIMELM; ++nelem )
01265 {
01266 DR_Badnell_rate_coef_mean_ion[nelem] =
01267 BadnelDR_RateSave[nelem] * RecNoise[nelem] *
01268
01269 ionbal.DR_mean_scale[nelem];
01270 }
01271
01272
01273
01274 for( long ion=0; ion < ipIRON+1; ++ion )
01275 {
01276 if( ionbal.DR_Badnell_rate_coef[ipIRON][ion] < 0. )
01277 {
01278 ionbal.DR_Badnell_rate_coef[ipIRON][ion] =
01279 DR_Badnell_rate_coef_mean_ion[ion]
01280 + atmdat_dielrec_fe(ion+1, phycon.te );
01281 strcpy(chDRDataSource[ipIRON][ion] , "mean+");
01282 }
01283 }
01284
01285 for( long nelem=0; nelem < LIMELM; ++nelem )
01286 {
01287 for( long ion=0; ion < nelem+1; ++ion )
01288 if( ionbal.DR_Badnell_rate_coef[nelem][ion] < 0. )
01289 {
01290 strcpy(chDRDataSource[nelem][ion] , "mean");
01291 ionbal.DR_Badnell_rate_coef[nelem][ion] = DR_Badnell_rate_coef_mean_ion[ion];
01292 }
01293 }
01294
01295
01296 for( long nelem=ipLITHIUM; nelem < LIMELM; ++nelem )
01297 {
01298 for( long ion=0; ion < nelem-1; ++ion )
01299 {
01300
01301
01302
01303
01304
01305 ionbal.DR_Badnell_rate_coef[nelem][ion] *= CollisSuppres(
01306
01307
01308 nelem+1,
01309
01310 ion+1,
01311
01312 dense.eden,
01313
01314 phycon.te );
01315
01316 ASSERT(ionbal.DR_Badnell_rate_coef[nelem][ion] >= 0);
01317 ASSERT(ionbal.RR_rate_coef_used[nelem][ion] >= 0);
01318 }
01319 }
01320
01321
01322 if( ionbal.lgRecom_Badnell_print )
01323 {
01324
01325 fprintf(ioQQQ,"\n\n RR recombination data sources \n" );
01326
01327 for( long loop=0;loop<30;loop+=10)
01328 {
01329 fprintf(ioQQQ,"\n\n ");
01330 for(long ion=loop; ion<loop+10; ++ion )
01331 {
01332 fprintf(ioQQQ,"&%7li",ion);
01333 }
01334 fprintf(ioQQQ,"\\\\\n" );
01335 for( long nelem=loop; nelem<LIMELM; ++nelem )
01336 {
01337 fprintf(ioQQQ,"%2li %5s ",nelem+1 , elementnames.chElementNameShort[nelem] );
01338 long limit = MIN2(nelem+1,loop+10);
01339 for( long ion=loop; ion<limit; ++ion )
01340 {
01341 fprintf(ioQQQ,"&%7s",chRRDataSource[nelem][ion] );
01342 }
01343 for( long ion=limit; ion<loop+10; ++ion )
01344 {
01345 fprintf(ioQQQ,"&%7s",chRRDataSource[nelem][ion] );
01346 }
01347 fprintf(ioQQQ,"\\\\\n" );
01348 }
01349 }
01350 fprintf(ioQQQ,"\nData sources\n");
01351 fprintf(ioQQQ,"Bad06: Badnell, N., 2006, ApJ, 167, 334B\n");
01352 fprintf(ioQQQ,"Verner: Verner & Ferland, 1996, ApJS, 103, 467\n");
01353
01354 fprintf(ioQQQ,"\n\n DR recombination data sources \n" );
01355
01356 for( long loop=0;loop<30;loop+=10)
01357 {
01358 fprintf(ioQQQ,"\n\n ");
01359 for(long ion=loop; ion<loop+10; ++ion )
01360 {
01361 fprintf(ioQQQ,"&%7li",ion);
01362 }
01363 fprintf(ioQQQ,"\\\\\n" );
01364 for( long nelem=loop; nelem<LIMELM; ++nelem )
01365 {
01366 fprintf(ioQQQ,"%2li %5s ",
01367 nelem+1 , elementnames.chElementNameShort[nelem] );
01368 long limit = MIN2(nelem+1,loop+10);
01369 for( long ion=loop; ion<limit; ++ion )
01370 {
01371 fprintf(ioQQQ,"&%7s",chDRDataSource[nelem][ion] );
01372 }
01373 for( long ion=limit; ion<loop+10; ++ion )
01374 {
01375 fprintf(ioQQQ,"&%7s",chDRDataSource[nelem][ion] );
01376 }
01377 fprintf(ioQQQ,"\\\\\n" );
01378 }
01379 }
01380 fprintf(ioQQQ,"\nData sources\nBadWeb: Badnell web site http://amdpp.phys.strath.ac.uk/tamoc/DR/\n");
01381 fprintf(ioQQQ,"Bad06D: Badnell, N., 2006, ApJ, 651, L73\n");
01382 fprintf(ioQQQ,"GuPC: Gu, M. private communication\n");
01383
01384 fprintf(ioQQQ,"\n\nDEBUG Badnell recombination RR, then DR, T=%.3e\n", phycon.te );
01385 for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
01386 {
01387 fprintf(ioQQQ,"nelem=%li %s, RR then DR\n",
01388 nelem , elementnames.chElementNameShort[nelem] );
01389 for( long ion=0; ion<nelem+1; ++ion )
01390 {
01391 fprintf(ioQQQ," %.2e", ionbal.RR_rate_coef_used[nelem][ion] );
01392 }
01393 fprintf(ioQQQ,"\n" );
01394 for( long ion=0; ion<nelem+1; ++ion )
01395 {
01396 fprintf(ioQQQ," %.2e", ionbal.DR_Badnell_rate_coef[nelem][ion] );
01397 }
01398 fprintf(ioQQQ,"\n\n" );
01399 }
01400
01401 fprintf(ioQQQ,"mean DR recombination ion mean \n" );
01402 for( long ion=0; ion<LIMELM; ++ion )
01403 {
01404 fprintf(ioQQQ," %2li %.2e \n",
01405 ion ,
01406 DR_Badnell_rate_coef_mean_ion[ion] );
01407 }
01408
01409 fprintf( ioQQQ, "\n\nCollisSuppres finds following dielectronic"
01410 " recom suppression factors, eden=%10.3e\n", dense.eden );
01411 fprintf( ioQQQ, "nelem ion fac \n" );
01412 for( long nelem=0; nelem<LIMELM; ++nelem )
01413 {
01414 for( long ion=0; ion < nelem+1; ion++ )
01415 {
01416 fprintf( ioQQQ, "%3ld %4ld %10.3e\n", nelem+1 , ion+1,
01417 CollisSuppres(
01418
01419
01420 nelem+1,
01421
01422 ion+1,
01423
01424 dense.eden,
01425
01426 phycon.te )
01427 );
01428
01429 }
01430 fprintf( ioQQQ, "\n");
01431 }
01432
01433 cdEXIT( EXIT_SUCCESS );
01434 }
01435 return;
01436 }