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 static long int *nsumrec;
00032
00033
00034 static bool **lgDRBadnellDefined ,
00035 **lgDRBadnellDefinedPart2,
00036 **lgRRBadnellDefined;
00037 static bool lgMustMallocRec=true;
00038
00039
00040
00041
00042
00043 #if defined(PRINT_DR) || defined(PRINT_RR)
00044 static const char FILE_NAME_OUT[] = "array.out";
00045 #endif
00046
00060 STATIC double Badnell_DR_rate_eval(
00061
00062 int nAtomicNumberCScale,
00063
00064 int n_core_e_before_recomb )
00065 {
00066
00067 double RateCoefficient, sum;
00068 int i;
00069
00070 DEBUG_ENTRY( "Badnell_DR_rate_eval()" );
00071 ASSERT( nAtomicNumberCScale>=0 && nAtomicNumberCScale<LIMELM );
00072
00073 if( nAtomicNumberCScale==ipIRON && n_core_e_before_recomb>=12 &&
00074 n_core_e_before_recomb<=18 )
00075 {
00076
00077
00078
00079
00080
00081
00082
00083 double cFe_q[7][8] =
00084 {
00085 {5.636e-4, 7.390e-3, 3.635e-2, 1.693e-1, 3.315e-2, 2.288e-1, 7.316e-2, 0.},
00086 {1.090e-3, 7.801e-3, 1.132e-2, 4.740e-2, 1.990e-1, 3.379e-2, 1.140e-1, 1.250e-1},
00087 {3.266e-3, 7.637e-3, 1.005e-2, 2.527e-2, 6.389e-2, 1.564e-1, 0., 0.},
00088 {1.074e-3, 6.080e-3, 1.887e-2, 2.540e-2, 7.580e-2, 2.773e-1, 0., 0.},
00089 {9.073e-4, 3.777e-3, 1.027e-2, 3.321e-2, 8.529e-2, 2.778e-1, 0., 0.},
00090 {5.335e-4, 1.827e-3, 4.851e-3, 2.710e-2, 8.226e-2, 3.147e-1, 0., 0.},
00091 {7.421e-4, 2.526e-3, 4.605e-3, 1.489e-2, 5.891e-2, 2.318e-1, 0., 0.}
00092 };
00093
00094
00095 double EFe_q[7][8] =
00096 {
00097 {3.628e3, 2.432e4, 1.226e5, 4.351e5, 1.411e6, 6.589e6, 1.030e7, 0},
00098 {1.246e3, 1.063e4, 4.719e4, 1.952e5, 5.637e5, 2.248e6, 7.202e6, 3.999e9},
00099 {1.242e3, 1.001e4, 4.466e4, 1.497e5, 3.919e5, 6.853e5, 0. , 0.},
00100 {1.387e3, 1.048e4, 3.955e4, 1.461e5, 4.010e5, 7.208e5, 0. , 0.},
00101 {1.525e3, 1.071e4, 4.033e4, 1.564e5, 4.196e5, 7.580e5, 0. , 0.},
00102 {2.032e3, 1.018e4, 4.638e4, 1.698e5, 4.499e5, 7.880e5, 0. , 0.},
00103 {3.468e3, 1.353e4, 3.690e4, 1.957e5, 4.630e5, 8.202e5, 0. , 0.}
00104 };
00105
00106 long int nion = n_core_e_before_recomb - 12;
00107 ASSERT( nion>=0 && nion <=6 );
00108
00109 sum = 0;
00110 i = 0;
00111
00112 for(i=0; i<8; ++i )
00113 {
00114 sum += (cFe_q[nion][i] * sexp( EFe_q[nion][i]/phycon.te));
00115 }
00116
00117
00118 RateCoefficient = sum / phycon.te32;
00119
00120 return RateCoefficient;
00121 }
00122
00123
00124 else if( nAtomicNumberCScale < n_core_e_before_recomb )
00125 {
00126 RateCoefficient = -2;
00127 }
00128
00129 else if( nAtomicNumberCScale >= LIMELM )
00130 {
00131 RateCoefficient = -2;
00132 }
00133
00134 else if( !lgDRBadnellDefined[nAtomicNumberCScale][n_core_e_before_recomb] )
00135 {
00136 RateCoefficient = -1;
00137 }
00138 else if( lgDRBadnellDefined[nAtomicNumberCScale][n_core_e_before_recomb] )
00139 {
00140
00141 sum = 0;
00142 i = 0;
00143
00144 for(i=0; i<nDRFitPar[nAtomicNumberCScale][n_core_e_before_recomb]; ++i )
00145 {
00146 sum += (DRFitParPart1[nAtomicNumberCScale][n_core_e_before_recomb][i] *
00147 sexp( DRFitParPart2[nAtomicNumberCScale][n_core_e_before_recomb][i]/phycon.te));
00148 }
00149
00150
00151 RateCoefficient = sum / phycon.te32;
00152 }
00153
00154 else
00155 {
00156 RateCoefficient = -99;
00157 }
00158
00159 return RateCoefficient;
00160 }
00161
00166 STATIC double Badnell_RR_rate_eval(
00167
00168 int nAtomicNumberCScale,
00169
00170 int n_core_e_before_recomb )
00171 {
00172 double RateCoefficient;
00173 double B, D, F;
00174
00175 DEBUG_ENTRY( "Badnell_RR_rate_eval()" );
00176
00177 ASSERT( nAtomicNumberCScale>=0 && nAtomicNumberCScale<LIMELM );
00178
00179 if( nAtomicNumberCScale==ipIRON &&
00180 n_core_e_before_recomb>=12 && n_core_e_before_recomb<=18 )
00181 {
00182
00183
00184
00185
00186
00187
00188
00189 double parFeq[7][6] ={
00190 {1.179e-9 , 0.7096, 4.508e2, 3.393e7, 0.0154, 3.977e6},
00191 {1.050e-9 , 0.6939, 4.568e2, 3.987e7, 0.0066, 5.451e5},
00192 {9.832e-10, 0.7146, 3.597e2, 3.808e7, 0.0045, 3.952e5},
00193 {8.303e-10, 0.7156, 3.531e2, 3.554e7, 0.0132, 2.951e5},
00194 {1.052e-9 , 0.7370, 1.639e2, 2.924e7, 0.0224, 4.291e5},
00195 {1.338e-9 , 0.7495, 7.242e1, 2.453e7, 0.0404, 4.199e5},
00196 {1.263e-9 , 0.7532, 5.209e1, 2.169e7, 0.0421, 2.917e5}
00197 };
00198
00199 double temp;
00200
00201 long int nion = n_core_e_before_recomb - 12;
00202 ASSERT( nion>=0 && nion <=6 );
00203
00204 temp = -parFeq[nion][5]/phycon.te;
00205 B = parFeq[nion][1] + parFeq[nion][4]*exp(temp);
00206 D = sqrt(phycon.te/parFeq[nion][2]);
00207 F = sqrt(phycon.te/parFeq[nion][3]);
00208 RateCoefficient = parFeq[nion][0]/(D*pow((1.+D),(1.-B))*pow((1.+F),(1.+B)));
00209
00210 return RateCoefficient;
00211 }
00212
00213
00214 else if( nAtomicNumberCScale < n_core_e_before_recomb )
00215 {
00216 RateCoefficient = -2;
00217 }
00218
00219 else if( nAtomicNumberCScale >= LIMELM )
00220 {
00221 RateCoefficient = -2;
00222 }
00223
00224 else if( !lgRRBadnellDefined[nAtomicNumberCScale][n_core_e_before_recomb] )
00225 {
00226 RateCoefficient = -1;
00227 }
00228
00229 else if( lgRRBadnellDefined[nAtomicNumberCScale][n_core_e_before_recomb] )
00230 {
00231
00232
00233
00234 double temp;
00235
00236 temp = -RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][5]/phycon.te;
00237 B = RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][1] +
00238 RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][4]*exp(temp);
00239 D = sqrt(phycon.te/RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][2]);
00240 F = sqrt(phycon.te/RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][3]);
00241 RateCoefficient = RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][0]/(D*pow((1.+D),(1.-B))*pow((1.+F),(1.+B)));
00242 }
00243
00244
00245 else
00246 RateCoefficient = -99;
00247
00248 return RateCoefficient;
00249 }
00250
00251
00252
00253
00254
00255 void Badnell_rec_init( void )
00256 {
00257
00258 double par_C[MAX_FIT_PAR_DR];
00259 double par_E[MAX_FIT_PAR_DR];
00260 char chLine[INPUT_LINE_LENGTH];
00261 int NuclearCharge=-1, NumberElectrons=-1;
00262 int i, k;
00263 int count, number;
00264 double temp_par[MAX_FIT_PAR_RR];
00265 int M_state, W_state,
00266 nelem;
00267
00268 const int NBLOCK = 2;
00269 int data_begin_line[NBLOCK];
00270 int length_of_line;
00271 FILE *ioDATA;
00272 const char* chFilename;
00273 int yr, mo, dy;
00274 char *chs;
00275
00276 const int BIGGEST_INDEX_TO_USE = 103;
00277
00278
00279 long TheirIndexToOurIndex[BIGGEST_INDEX_TO_USE];
00280 char string[120];
00281 double value;
00282 bool lgEOL;
00283 long int i1, ipISO = ipHE_LIKE;
00284 long INDX=0,INDP=0,N=0,S=0,L=0,J=0,maxINDX=0,loopindex=0,max_N_of_data=-1;
00285 bool lgFlag = true;
00286
00287 static int nCalled = 0;
00288
00289 const char* cdDATAFILE[] =
00290 {
00291
00292 "",
00293 "nrb00#h_he1ic12.dat",
00294 "nrb00#h_li2ic12.dat",
00295 "nrb00#h_be3ic12.dat",
00296 "nrb00#h_b4ic12.dat",
00297 "nrb00#h_c5ic12.dat",
00298 "nrb00#h_n6ic12.dat",
00299 "nrb00#h_o7ic12.dat",
00300 "nrb00#h_f8ic12.dat",
00301 "nrb00#h_ne9ic12.dat",
00302 "nrb00#h_na10ic12.dat",
00303 "nrb00#h_mg11ic12.dat",
00304 "nrb00#h_al12ic12.dat",
00305 "nrb00#h_si13ic12.dat",
00306 "nrb00#h_p14ic12.dat",
00307 "nrb00#h_s15ic12.dat",
00308 "nrb00#h_cl16ic12.dat",
00309 "nrb00#h_ar17ic12.dat",
00310 "nrb00#h_k18ic12.dat",
00311 "nrb00#h_ca19ic12.dat",
00312 "nrb00#h_sc20ic12.dat",
00313 "nrb00#h_ti21ic12.dat",
00314 "nrb00#h_v22ic12.dat",
00315 "nrb00#h_cr23ic12.dat",
00316 "nrb00#h_mn24ic12.dat",
00317 "nrb00#h_fe25ic12.dat",
00318 "nrb00#h_co26ic12.dat",
00319 "nrb00#h_ni27ic12.dat",
00320 "nrb00#h_cu28ic12.dat",
00321 "nrb00#h_zn29ic12.dat"
00322 };
00323
00324
00325 DEBUG_ENTRY( "Badnell_rec_init()" );
00326
00327
00328 if( nCalled > 0 )
00329 {
00330 return;
00331 }
00332 ++nCalled;
00333
00334 # if defined(PRINT_DR) || defined(PRINT_RR)
00335 FILE *ofp = open_data( FILE_NAME_OUT, "w", AS_LOCAL_ONLY );
00336 # endif
00337
00338
00339
00340 for(nelem=ipHELIUM;nelem<LIMELM;nelem++)
00341 {
00342 if( nelem < 2 || dense.lgElmtOn[nelem] )
00343 {
00344 ioDATA= open_data( cdDATAFILE[nelem], "r" );
00345
00346 lgFlag = true;
00347 ASSERT(ioDATA);
00348
00349 for( i=0; i<BIGGEST_INDEX_TO_USE; i++ )
00350 TheirIndexToOurIndex[i] = -1;
00351
00352
00353 while(lgFlag)
00354 {
00355 if(read_whole_line(string,sizeof(string),ioDATA)!=NULL)
00356 {
00357 if( nMatch("INDX INDP ",string) )
00358 {
00359
00360 if( read_whole_line( string , (int)sizeof(string) , ioDATA ) == NULL )
00361 {
00362 fprintf( ioQQQ, " Badnell data file appears to be corrupted.\n");
00363 cdEXIT(EXIT_FAILURE);
00364 }
00365
00366
00367 while( read_whole_line(string, (int)sizeof(string), ioDATA) != NULL )
00368 {
00369 if( strcmp(string,"\n")==0 )
00370 {
00371 lgFlag = false;
00372 break;
00373 }
00374
00375 i1=3;
00376 INDX=(long)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
00377 if( INDX >= BIGGEST_INDEX_TO_USE )
00378 {
00379 INDX--;
00380 lgFlag = false;
00381 break;
00382 }
00383
00384 ASSERT( INDX < BIGGEST_INDEX_TO_USE );
00385
00386 INDP=(long)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
00387 ASSERT( INDP >= 1 );
00388
00389 if(INDP==1)
00390 {
00391 if( (i1=nMatch("1S1 ",string)) > 0 )
00392 {
00393 i1 += 4;
00394 N=(long)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
00395 ASSERT( N>=1 );
00396 }
00397 else
00398 {
00399 TotalInsanity();
00400 }
00401
00402 if( (i1=nMatch(" (",string)) > 0 )
00403 {
00404 i1 += 6;
00405 S=(long)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
00406
00407 ASSERT( S==1 || S==3 );
00408 }
00409 else
00410 {
00411 TotalInsanity();
00412 }
00413
00414
00415 i1++;
00416 L=(long)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
00417 ASSERT( L >= 0 && L < N );
00418
00419
00420 i1 += 2;
00421 J=(long)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
00422 ASSERT( J <= ( L + (int)((S+1)/2) ) &&
00423 J >= ( L - (int)((S+1)/2) ) && J >= 0 );
00424
00425
00426 if( N<= iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] )
00427 TheirIndexToOurIndex[INDX] = iso.QuantumNumbers2Index[ipHE_LIKE][nelem][N][L][S];
00428 else
00429 {
00430
00431
00432 INDX--;
00433 lgFlag = false;
00434 break;
00435 }
00436
00437
00438 if( N==2 && L==1 && S==3 )
00439 {
00440 if( J==0 )
00441 TheirIndexToOurIndex[INDX] = 3;
00442 else if( J==1 )
00443 TheirIndexToOurIndex[INDX] = 4;
00444 else
00445 {
00446 ASSERT( J==2 );
00447 ASSERT( TheirIndexToOurIndex[INDX] == 5 );
00448 }
00449 }
00450 max_N_of_data = MAX2( max_N_of_data, N );
00451 }
00452 else
00453 {
00454
00455 lgFlag = false;
00456 }
00457 }
00458 }
00459 }
00460 else
00461 {
00462
00463 lgFlag = false;
00464 }
00465 }
00466
00467 maxINDX =INDX;
00468 ASSERT( maxINDX > 0 );
00469 ASSERT( maxINDX < BIGGEST_INDEX_TO_USE );
00470
00471 INDX = 0;
00472 lgFlag = true;
00473 while(lgFlag)
00474 {
00475 if(read_whole_line(string,sizeof(string),ioDATA)!=NULL)
00476 {
00477
00478 if( nMatch("INDX TE= ",string) )
00479 {
00480 lgFlag = false;
00481
00482
00483 if( read_whole_line( string , (int)sizeof(string) , ioDATA ) == NULL )
00484 {
00485 fprintf( ioQQQ, " Badnell data file appears to be corrupted.\n");
00486 cdEXIT(EXIT_FAILURE);
00487 }
00488
00489
00490 while( read_whole_line(string, (int)sizeof(string), ioDATA) != NULL )
00491 {
00492
00493 if( nMatch("PRTF",string) || INDX >= maxINDX || INDX<0 )
00494 break;
00495
00496 i1=3;
00497 INDX=(long)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
00498 if( INDX>maxINDX )
00499 break;
00500
00501 for(loopindex=0;loopindex<10;loopindex++)
00502 {
00503 value=(double)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
00504 if( TheirIndexToOurIndex[INDX] < iso.numLevels_max[ipHE_LIKE][nelem] &&
00505 TheirIndexToOurIndex[INDX] > 0 )
00506 {
00507 iso.DielecRecombVsTemp[ipHE_LIKE][nelem][TheirIndexToOurIndex[INDX]][loopindex] += value;
00508 }
00509 }
00510
00511
00512 if( read_whole_line( string , (int)sizeof(string) , ioDATA ) == NULL )
00513 {
00514 fprintf( ioQQQ, " Badnell data file appears to be corrupted.\n");
00515 cdEXIT(EXIT_FAILURE);
00516 }
00517
00518
00519 i1 = 13;
00520 for(loopindex=10;loopindex<19;loopindex++)
00521 {
00522 value=(double)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
00523 if( TheirIndexToOurIndex[INDX] < iso.numLevels_max[ipHE_LIKE][nelem] &&
00524 TheirIndexToOurIndex[INDX] > 0 )
00525 {
00526 iso.DielecRecombVsTemp[ipHE_LIKE][nelem][TheirIndexToOurIndex[INDX]][loopindex] += value;
00527 }
00528 }
00529 }
00530 }
00531 }
00532 else
00533 lgFlag = false;
00534 }
00535 fclose(ioDATA);
00536 ASSERT( maxINDX > 0 );
00537 ASSERT( maxINDX < BIGGEST_INDEX_TO_USE );
00538 ASSERT( max_N_of_data > 0 );
00539
00540 if( max_N_of_data < iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] )
00541 {
00542 long indexOfMaxN;
00543 L = -1;
00544 S = -1;
00545
00546
00547 for( i=TheirIndexToOurIndex[maxINDX]+1;
00548 i<iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem]; i++ )
00549 {
00550 L = L_(i);
00551 S = S_(i);
00552
00553 if( L > 4 )
00554 continue;
00555
00556 indexOfMaxN = iso.QuantumNumbers2Index[ipHE_LIKE][nelem][max_N_of_data][L][S];
00557 for(loopindex=0;loopindex<19;loopindex++)
00558 {
00559 iso.DielecRecombVsTemp[ipHE_LIKE][nelem][i][loopindex] =
00560 iso.DielecRecombVsTemp[ipHE_LIKE][nelem][indexOfMaxN][loopindex] *
00561 pow( (double)max_N_of_data/(double)StatesElemNEW[nelem][nelem-ipHE_LIKE][i].n, 3.);
00562 }
00563 }
00564
00565
00566 indexOfMaxN =
00567 iso.QuantumNumbers2Index[ipHE_LIKE][nelem][iso.n_HighestResolved_max[ipHE_LIKE][nelem]][1][1];
00568
00569
00570 for( i=iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem];
00571 i<iso.numLevels_max[ipHE_LIKE][nelem]; i++ )
00572 {
00573 for(loopindex=0;loopindex<19;loopindex++)
00574 {
00575 iso.DielecRecombVsTemp[ipHE_LIKE][nelem][i][loopindex] =
00576 iso.DielecRecombVsTemp[ipHE_LIKE][nelem][indexOfMaxN][loopindex] *
00577 pow( (double)iso.n_HighestResolved_max[ipHE_LIKE][nelem]/
00578 (double)StatesElemNEW[nelem][nelem-ipHE_LIKE][i].n, 3.);
00579 }
00580 }
00581 }
00582 }
00583 }
00584
00585 for( i=0; i<NBLOCK; ++i )
00586 {
00587
00588 data_begin_line[i] = INT_MIN;
00589 }
00590
00591 chFilename = "badnell_dr.dat";
00592 ioDATA = open_data( chFilename, "r" );
00593
00594 count = 0;
00595 number = 0;
00596
00597
00598
00599 while( read_whole_line(chLine, (int)sizeof(chLine), ioDATA) != NULL )
00600 {
00601 count++;
00602
00603 if( chLine[2]=='Z' )
00604 {
00605
00606
00607 data_begin_line[number] = count;
00608 ASSERT( number < NBLOCK );
00609 number++;
00610 }
00611 }
00612
00613
00614 if( lgMustMallocRec )
00615 {
00616 nsumrec = (long *)MALLOC(LIMELM*sizeof(long) );
00617 nDRFitPar = (int**)MALLOC( LIMELM*sizeof( int*) );
00618 lgDRBadnellDefined = (bool **)MALLOC( LIMELM*sizeof(bool*) );
00619 lgDRBadnellDefinedPart2 = (bool **)MALLOC( LIMELM*sizeof(bool*) );
00620 lgRRBadnellDefined = (bool **)MALLOC( LIMELM*sizeof(bool*) );
00621
00622 DRFitParPart1 = (double ***)MALLOC( LIMELM*sizeof(double**) );
00623 DRFitParPart2 = (double ***)MALLOC( LIMELM*sizeof(double**) );
00624 RRFitPar = (double ***)MALLOC( LIMELM*sizeof(double**) );
00625 }
00626
00627 for( nelem=0; nelem<LIMELM; nelem++ )
00628 {
00629 if( lgMustMallocRec )
00630 {
00631 nDRFitPar[nelem] = (int*)MALLOC( (nelem+1)*sizeof( int) );
00632 lgDRBadnellDefined[nelem] = (bool *)MALLOC( (nelem+1)*sizeof(bool) );
00633 lgDRBadnellDefinedPart2[nelem] = (bool *)MALLOC( (nelem+1)*sizeof(bool) );
00634 lgRRBadnellDefined[nelem] = (bool *)MALLOC( (nelem+1)*sizeof(bool) );
00635
00636 DRFitParPart1[nelem] = (double **)MALLOC( (nelem+1)*sizeof(double*) );
00637 DRFitParPart2[nelem] = (double **)MALLOC( (nelem+1)*sizeof(double*) );
00638 RRFitPar[nelem] = (double **)MALLOC( (nelem+1)*sizeof(double*) );
00639 }
00640 for( long ion=0; ion<nelem+1; ++ion )
00641 {
00642 if( lgMustMallocRec )
00643 {
00644 DRFitParPart1[nelem][ion] = (double *)MALLOC( MAX_FIT_PAR_DR*sizeof(double) );
00645 DRFitParPart2[nelem][ion] = (double *)MALLOC( MAX_FIT_PAR_DR*sizeof(double) );
00646 RRFitPar[nelem][ion] = (double *)MALLOC( MAX_FIT_PAR_RR*sizeof(double) );
00647 }
00648 lgDRBadnellDefined[nelem][ion] = false;
00649 lgDRBadnellDefinedPart2[nelem][ion] = false;
00650 lgRRBadnellDefined[nelem][ion] = false;
00651
00652
00653 for( k=0; k<MAX_FIT_PAR_DR; k++ )
00654 {
00655 DRFitParPart1[nelem][ion][k] = 0;
00656 DRFitParPart2[nelem][ion][k] = 0;
00657 }
00658 for( k=0; k<MAX_FIT_PAR_RR; k++ )
00659 {
00660 RRFitPar[nelem][ion][k] = 0;
00661 }
00662 }
00663 }
00664 lgMustMallocRec = false;
00665
00666 count = 0;
00667
00668 fseek(ioDATA, 0, SEEK_SET);
00669
00670 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00671 {
00672 fprintf( ioQQQ, " DISASTER PROBLEM Badnell_rec_init could not read first line of badnell_dr.dat.\n");
00673 cdEXIT(EXIT_FAILURE);
00674 }
00675 count++;
00676
00677
00678 if( (chs = strchr_s(chLine, ')'))==NULL )
00679 {
00680
00681 fprintf( ioQQQ, " DISASTER PROBLEM Badnell_rec_init data file incorrect format.\n");
00682 cdEXIT(EXIT_FAILURE);
00683 }
00684
00685 ++chs;
00686 sscanf(chs, "%4i%2i%2i",&yr, &mo, &dy);
00687
00688 int dr_yr = 2007, dr_mo = 10, dr_dy = 27;
00689 if((yr != dr_yr) || (mo != dr_mo) || (dy != dr_dy))
00690 {
00691 fprintf(ioQQQ,
00692 "DISASTER PROBLEM Badnell_rec_init The version of %s I found (%i %i %i) is not the current version (%i %i %i).\n",
00693 chFilename, yr, mo, dy, dr_yr, dr_mo, dr_dy);
00694 fprintf(ioQQQ," The first line of the file is the following\n %s\n", chLine );
00695 cdEXIT(EXIT_FAILURE);
00696 }
00697
00698 while( read_whole_line(chLine, (int)sizeof(chLine), ioDATA) != NULL )
00699 {
00700 count++;
00701 length_of_line = (int)strlen(chLine);
00702
00703
00704 if( count > data_begin_line[0] && count < data_begin_line[1] && length_of_line >3 )
00705 {
00706
00707 for( i=0; i<MAX_FIT_PAR_DR; i++ )
00708 {
00709 par_C[i] = 0;
00710 }
00711 sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf%lf%lf%lf",
00712 &NuclearCharge, &NumberElectrons, &M_state, &W_state, &par_C[0], &par_C[1], &par_C[2],
00713 &par_C[3], &par_C[4], &par_C[5], &par_C[6], &par_C[7], &par_C[8]);
00714
00715
00716 long int NuclearChargeM1 = NuclearCharge-1;
00717
00718 if(M_state == 1 && NuclearChargeM1 < LIMELM )
00719 {
00720
00721 ASSERT( NumberElectrons < LIMELM );
00722 ASSERT( NuclearChargeM1 < LIMELM );
00723 lgDRBadnellDefined[NuclearChargeM1][NumberElectrons] = true;
00724
00725
00726 nDRFitPar[NuclearChargeM1][NumberElectrons] = 9;
00727 for( i=8; i>=0; i-- )
00728 {
00729 if( par_C[i] == 0 )
00730 --nDRFitPar[NuclearChargeM1][NumberElectrons];
00731 else
00732 break;
00733 }
00734
00735
00736 for( i=0; i<9; i++ )
00737 DRFitParPart1[NuclearChargeM1][NumberElectrons][i] = par_C[i];
00738 }
00739 }
00740 }
00741
00742
00743 fseek(ioDATA, 0, SEEK_SET);
00744 count = 0;
00745 while( read_whole_line(chLine, (int)sizeof(chLine), ioDATA) != NULL )
00746 {
00747 count++;
00748 length_of_line = (int)strlen(chLine);
00749 if( count > data_begin_line[1] && length_of_line > 3 )
00750 {
00751
00752
00753 for( i=0; i<MAX_FIT_PAR_DR; i++ )
00754 {
00755 par_E[i] = 0;
00756 }
00757 sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf%lf%lf%lf",
00758 &NuclearCharge, &NumberElectrons, &M_state, &W_state, &par_E[0], &par_E[1], &par_E[2],
00759 &par_E[3], &par_E[4], &par_E[5], &par_E[6], &par_E[7], &par_E[8]);
00760
00761 long int NuclearChargeM1 = NuclearCharge-1;
00762
00763 if(M_state == 1 && NuclearChargeM1<LIMELM)
00764 {
00765 ASSERT( NumberElectrons < LIMELM );
00766 ASSERT( NuclearChargeM1 < LIMELM );
00767 lgDRBadnellDefinedPart2[NuclearChargeM1][NumberElectrons] = true;
00768
00769
00770 nDRFitPar[NuclearChargeM1][NumberElectrons] = 9;
00771 for( i=8; i>=0; i-- )
00772 {
00773 if( par_E[i] == 0 )
00774 --nDRFitPar[NuclearChargeM1][NumberElectrons];
00775 else
00776 break;
00777 }
00778
00779
00780 for( i=0; i<nDRFitPar[NuclearChargeM1][NumberElectrons]; i++ )
00781 DRFitParPart2[NuclearChargeM1][NumberElectrons][i] = par_E[i];
00782 }
00783 }
00784 }
00785
00786 fclose( ioDATA );
00787
00788
00789 # ifdef PRINT_DR
00790 for( nelem=0; nelem<LIMELM; nelem++ )
00791 {
00792 for( int ion=0; ion<nelem+1;++ion )
00793 {
00794 if( lgDRBadnellDefined[nelem][ion] )
00795 {
00796 fprintf(ofp, "%i %i %e %e %e %e %e %e %e %e %e\n",
00797 nelem, ion, DRFitParPart1[nelem][ion][0],
00798 DRFitParPart1[nelem][ion][1], DRFitParPart1[nelem][ion][2],
00799 DRFitParPart1[nelem][ion][3], DRFitParPart1[nelem][ion][4],
00800 DRFitParPart1[nelem][ion][5], DRFitParPart1[nelem][ion][6],
00801 DRFitParPart1[nelem][ion][7], DRFitParPart1[nelem][ion][8]);
00802 }
00803 }
00804 }
00805 for( nelem=0; nelem<LIMELM; nelem++ )
00806 {
00807 for( int ion=0; ion<nelem+1; ion++ )
00808 {
00809 if( lgDRBadnellDefinedPart2[nelem][ion] )
00810 {
00811 fprintf(ofp, "%i %i %e %e %e %e %e %e %e %e %e\n",
00812 nelem, ion, DRFitParPart2[nelem][ion][0],
00813 DRFitParPart2[nelem][ion][1], DRFitParPart2[nelem][ion][2],
00814 DRFitParPart2[nelem][ion][3], DRFitParPart2[nelem][ion][4],
00815 DRFitParPart2[nelem][ion][5], DRFitParPart2[nelem][ion][6],
00816 DRFitParPart2[nelem][ion][7], DRFitParPart2[nelem][ion][8]);
00817 }
00818 }
00819 }
00820 fclose(ofp);
00821 # endif
00822
00823
00824
00825 bool lgDRBadnellBothDefined = true;
00826 for( nelem=0; nelem<LIMELM; nelem++ )
00827 {
00828 for( int ion=0; ion<nelem+1; ion++ )
00829 {
00830
00831
00832 if( lgDRBadnellDefined[nelem][ion] != lgDRBadnellDefinedPart2[nelem][ion] )
00833 {
00834 fprintf( ioQQQ, "DR %i, RR %i: %c %c\n", nelem, ion,
00835 TorF(lgDRBadnellDefined[nelem][ion]),
00836 TorF(lgDRBadnellDefinedPart2[nelem][ion]));
00837 fprintf( ioQQQ, "PROBLEM ion_recomb_Badnell first and second half of Badnell DR not consistent.\n");
00838 lgDRBadnellBothDefined = false;
00839 }
00840 }
00841 }
00842
00843 if( !lgDRBadnellBothDefined )
00844 {
00845
00846 fprintf(ioQQQ,
00847 "DISASTER PROBLEM The DR data files are corrupted - part 1 and 2 do not agree.\n");
00848 fprintf(ioQQQ," Start again with a fresh copy of the data directory\n" );
00849 cdEXIT(EXIT_FAILURE);
00850 }
00851
00852
00853 chFilename = "badnell_rr.dat";
00854 ioDATA = open_data( chFilename, "r" );
00855
00856
00857 {
00858 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00859 {
00860 fprintf( ioQQQ, " DISASTER PROBLEM Badnell_rec_init could not read first line of badnell_rr.dat.\n");
00861 cdEXIT(EXIT_FAILURE);
00862 }
00863
00864 if( (chs = strchr_s(chLine, ')'))==NULL )
00865 {
00866
00867 fprintf( ioQQQ, " DISASTER PROBLEM Badnell_rec_init data file incorrect format.\n");
00868 cdEXIT(EXIT_FAILURE);
00869 }
00870 ++chs;
00871 sscanf(chs, "%4i%2i%2i", &yr, &mo, &dy);
00872 int rr_yr = 2007, rr_mo = 10, rr_dy = 27;
00873 if((yr != rr_yr)||(mo != rr_mo)||(dy != rr_dy))
00874 {
00875 fprintf(ioQQQ,"DISASTER PROBLEM The version of %s I found (%i %i %i) is not the current version (%i %i %i).\n",
00876 chFilename, yr, mo, dy, rr_yr, rr_mo, rr_dy);
00877 fprintf(ioQQQ," The line was as follows:\n %s\n", chLine );
00878 cdEXIT(EXIT_FAILURE);
00879 }
00880 }
00881
00882 while( read_whole_line(chLine, (int)sizeof(chLine), ioDATA) != NULL )
00883 {
00884
00885 for( i=0; i<MAX_FIT_PAR_RR; i++ )
00886 {
00887 temp_par[i] = 0;
00888 }
00889 if(chLine[0] != '#')
00890 {
00891 sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf",
00892 &NuclearCharge, &NumberElectrons, &M_state, &W_state, &temp_par[0], &temp_par[1],
00893 &temp_par[2], &temp_par[3], &temp_par[4], &temp_par[5]);
00894 long NuclearChargeM1 = NuclearCharge-1;
00895
00896 if(M_state == 1 && NuclearChargeM1<LIMELM)
00897 {
00898 ASSERT( NuclearChargeM1 < LIMELM );
00899 ASSERT( NumberElectrons <= LIMELM );
00900
00901 lgRRBadnellDefined[NuclearChargeM1][NumberElectrons] = true;
00902
00903 for( i=0; i<MAX_FIT_PAR_RR; i++ )
00904 RRFitPar[NuclearChargeM1][NumberElectrons][i] = temp_par[i];
00905 }
00906 }
00907 }
00908
00909
00910 # ifdef PRINT_RR
00911 count = 0;
00912 for( nelem=0; nelem<LIMELM; nelem++ )
00913 {
00914 for( ion=0; ion<nelem+1; ion++ )
00915 {
00916 if( lgRRBadnellDefined[nelem][ion] )
00917 {
00918 fprintf(ofp, "%i %i %e %e %e %e %e %e\n",
00919 nelem, ion, RRFitPar[nelem][ion][0],
00920 RRFitPar[nelem][ion][1], RRFitPar[nelem][ion][2],
00921 RRFitPar[nelem][ion][3],
00922 RRFitPar[nelem][ion][4], RRFitPar[nelem][ion][5]);
00923 count++;
00924 }
00925 }
00926 }
00927 fprintf(ofp, "total lines are %i ", count);
00928
00929 fclose(ofp);
00930 # endif
00931
00932 fclose(ioDATA);
00933
00934 {
00935 enum {DEBUG_LOC=false};
00936 if( DEBUG_LOC )
00937 {
00938 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00939 {
00940 int ion;
00941 fprintf(ioQQQ,"\nDEBUG rr rec\t%i",nelem);
00942 for( ion=0; ion<=nelem; ++ion )
00943 {
00944 fprintf(ioQQQ,"\t%.2e", Badnell_RR_rate_eval(nelem+1 , nelem-ion ) );
00945 }
00946 fprintf(ioQQQ,"\n");
00947 fprintf(ioQQQ,"DEBUG dr rec\t%i",nelem);
00948 for( ion=0; ion<=nelem; ++ion )
00949 {
00950 fprintf(ioQQQ,"\t%.2e", Badnell_DR_rate_eval(nelem+1 , nelem-ion ) );
00951 }
00952 fprintf(ioQQQ,"\n");
00953 }
00954 cdEXIT(EXIT_FAILURE);
00955 }
00956 }
00957 return;
00958 }
00959
00960
00961 void ion_recom_calculate( void )
00962 {
00963 static double TeUsed = -1;
00964 long int ion ,
00965 nelem ,
00966 i;
00967
00968 DEBUG_ENTRY( "ion_recom_calculate()" );
00969
00970
00971 if( fp_equal(phycon.te,TeUsed) )
00972 return;
00973
00974
00975 for( ion=0; ion < LIMELM; ++ion )
00976 nsumrec[ion] = 0;
00977
00978 TeUsed = phycon.te;
00979
00980
00981
00982 multi_arr<double,2> DR_Badnell_rate_stack( LIMELM, LIMELM );
00983 for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
00984 {
00985 for( ion=0; ion < nelem+1; ++ion )
00986 {
00987 long int n_bnd_elec_before_recom ,
00988 n_bnd_elec_after_recom;
00989
00990 n_bnd_elec_before_recom = nelem-ion;
00991 n_bnd_elec_after_recom = nelem-ion+1;
00992
00993 if( (ionbal.DR_Badnell_rate_coef[nelem][ion] =
00994 Badnell_DR_rate_eval(
00995
00996 nelem,
00997
00998
00999 n_bnd_elec_before_recom )) >= 0. )
01000 {
01001 ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = true;
01002
01003 DR_Badnell_rate_stack[ion][nsumrec[ion]++] =
01004 ionbal.DR_Badnell_rate_coef[nelem][ion];
01005 }
01006 else
01007 {
01008
01009 ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = false;
01010 ionbal.DR_Badnell_rate_coef[nelem][ion] = 0.;
01011 }
01012
01013
01014 if( (ionbal.RR_Badnell_rate_coef[nelem][ion] = Badnell_RR_rate_eval(
01015
01016 nelem,
01017
01018 n_bnd_elec_before_recom )) >= 0. )
01019 {
01020 ionbal.lgRR_Badnell_rate_coef_exist[nelem][ion] = true;
01021 }
01022 else
01023 {
01024
01025 ionbal.lgRR_Badnell_rate_coef_exist[nelem][ion] = false;
01026 ionbal.RR_Badnell_rate_coef[nelem][ion] = 0.;
01027 }
01028
01029
01030
01031 ionbal.RR_Verner_rate_coef[nelem][ion] =
01032 t_ADfA::Inst().rad_rec(
01033
01034 nelem+1 ,
01035
01036 n_bnd_elec_after_recom ,
01037 phycon.te );
01038 }
01039 }
01040
01041
01042 double Fe_Gu_c[9][6] = {
01043 { 2.50507e-11, 5.60226e-11, 1.85001e-10, 3.57495e-9, 1.66321e-7, 0. },
01044 { 9.19610e-11, 2.92460e-10, 1.02120e-9, 1.14852e-8, 3.25418e-7, 0. },
01045 { 9.02625e-11, 6.22962e-10, 5.77545e-9, 1.78847e-8, 3.40610e-7, 0. },
01046 { 9.04286e-12, 9.68148e-10, 4.83636e-9, 2.48159e-8, 3.96815e-7, 0. },
01047 { 6.77873e-10, 1.47252e-9, 5.31202e-9, 2.54793e-8, 3.47407e-7, 0. },
01048 { 1.29742e-9, 4.10172e-9, 1.23605e-8, 2.33615e-8, 2.97261e-7, 0. },
01049 { 8.78027e-10, 2.31680e-9, 3.49333e-9, 1.16927e-8, 8.18537e-8, 1.54740e-7 },
01050 { 2.23178e-10, 1.87313e-9, 2.86171e-9, 1.38575e-8, 1.17803e-7, 1.06251e-7 },
01051 { 2.17263e-10, 7.35929e-10, 2.81276e-9, 1.32411e-8, 1.15761e-7, 4.80389e-8 }
01052 },
01053
01054 Fe_Gu_E[9][6] = {
01055 { 8.30501e-2, 8.52897e-1, 3.40225e0, 2.23053e1, 6.80367e1, 0. },
01056 { 1.44392e-1, 9.23999e-1, 5.45498e0, 2.04301e1, 7.06112e1, 0. },
01057 { 5.79132e-2, 1.27852e0, 3.22439e0, 1.79602e1, 6.96277e1, 0. },
01058 { 1.02421e-1, 1.79393e0, 4.83226e0, 1.91117e1, 6.80858e1, 0. },
01059 { 1.24630e-1, 6.86045e-1, 3.09611e0, 1.44023e1, 6.42820e1, 0. },
01060 { 1.34459e-1, 6.63028e-1, 2.61753e0, 1.30392e1, 6.10222e1, 0. },
01061 { 7.79748e-2, 5.35522e-1, 1.88407e0, 8.38459e0, 3.38613e1, 7.89706e1 },
01062 { 8.83019e-2, 6.12756e-1, 2.36035e0, 9.61736e0, 3.64467e1, 8.72406e1 },
01063 { 1.51322e-1, 5.63155e-1, 2.57013e0, 9.08166e0, 3.69528e1, 1.08067e2 }
01064 };
01065
01066 double s_c[5][2] = {
01067 {0.1565e-3, 0.1617e-2},
01068 {0.3026e-3, 0.2076e-1},
01069 {0.3177e-1, 0.6309e-3},
01070 {0.2464e-1, 0.5553e-3},
01071 {0.1924e-1, 0.5111e-3}
01072 },
01073
01074 s_E[5][2] = {
01075 {0.1157e6, 0.1868e6},
01076 {0.9662e5, 0.1998e6},
01077 {0.1928e6, 0.6126e5},
01078 {0.1806e6, 0.3142e5},
01079 {0.1519e6, 0.4978e5}
01080 };
01081
01082
01083 nelem = ipIRON;
01084
01085
01086 double fitSum = 0;
01087
01088
01089 ion = 13;
01090
01091
01092
01093
01094 double fe14_c[6]={7.07E-04,7.18E-03,2.67E-02,3.15E-02,1.62E-01,5.37E-04};
01095
01096
01097
01098 double fe14_E[6]={4.12E-01,2.06E+00,1.03E+01,2.20E+01,4.22E+01,3.41E+03};
01099
01100
01101 for( i=0; i<6; i++ )
01102 {
01103 fitSum += fe14_c[i] * sexp( fe14_E[i]/phycon.te_eV );
01104 }
01105
01106 ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = true;
01107 ionbal.DR_Badnell_rate_coef[nelem][ion] = fitSum / phycon.te32;
01108 DR_Badnell_rate_stack[ion][nsumrec[ion]++] =
01109 ionbal.DR_Badnell_rate_coef[nelem][ion];
01110
01111
01112
01113
01114
01115 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 },
01116 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 };
01117
01118 ion = 12;
01119 if( !ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] )
01120 {
01121 for( i=0; i<10; i++ )
01122 {
01123 fitSum += fe13_c[i] * sexp( fe13_E[i]/phycon.te_eV );
01124 }
01125
01126
01127 ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = true;
01128 ionbal.DR_Badnell_rate_coef[nelem][ion] = fitSum / phycon.te32;
01129 DR_Badnell_rate_stack[ion][nsumrec[ion]++] =
01130 ionbal.DR_Badnell_rate_coef[nelem][ion];
01131 }
01132
01133
01134 double te_eV32 = pow( phycon.te_eV, 1.5 );
01135
01136
01137
01138
01139
01140 for( ion=0; ion<9; ion++ )
01141 {
01142
01143 if( !ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion+5] )
01144 {
01145 fitSum = 0;
01146 for( i=0; i<6; i++ )
01147 {
01148 fitSum += Fe_Gu_c[ion][i] * sexp( Fe_Gu_E[ion][i]/phycon.te_eV );
01149 }
01150 ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion+5] = true;
01151 ionbal.DR_Badnell_rate_coef[nelem][ion+5] = fitSum / te_eV32;
01152 DR_Badnell_rate_stack[ion+5][nsumrec[ion+5]++] =
01153 ionbal.DR_Badnell_rate_coef[nelem][ion+5];
01154 }
01155 }
01156
01157
01158
01159 for( ion=0; ion < LIMELM; ++ion )
01160 {
01161 ionbal.DR_Badnell_rate_coef_mean_ion[ion] = 0.;
01162 if( nsumrec[ion] > 0 )
01163 {
01164
01165
01166
01167
01168
01169
01170
01171
01172 sort( DR_Badnell_rate_stack.ptr( ion, 0 ), DR_Badnell_rate_stack.ptr( ion, nsumrec[ion] ) );
01173
01174 double maxval = DR_Badnell_rate_stack[ion][nsumrec[ion]-1];
01175 if( maxval > 0. )
01176 {
01177
01178
01179
01180
01181 double sum = 0.;
01182 for( int j = nsumrec[ion]/4; j < nsumrec[ion]; ++j )
01183 {
01184 double rate_one = log(max(DR_Badnell_rate_stack[ion][j],1e-5*maxval));
01185 sum += rate_one;
01186 }
01187 sum /= (double)(nsumrec[ion] - nsumrec[ion]/4);
01188 ionbal.DR_Badnell_rate_coef_mean_ion[ion] = exp(sum);
01189 }
01190 }
01191
01192 }
01193
01194
01195
01196 nelem = ipSULPHUR;
01197
01198 for( ion=0; ion<5; ion++ )
01199 {
01200
01201
01202 if( !ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] )
01203 {
01204
01205
01206 fitSum = 0;
01207 for( i=0; i<2; i++ )
01208 {
01209 fitSum += s_c[ion][i] * sexp( s_E[ion][i]/phycon.te );
01210 }
01211 ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = true;
01212 ionbal.DR_Badnell_rate_coef[nelem][ion] = fitSum / phycon.te32;
01213
01214
01215
01216
01217
01218
01219 if( ionbal.nDR_S_guess==0 )
01220 {
01221
01222 ionbal.DR_Badnell_rate_coef[nelem][ion] =
01223 MAX2( ionbal.DR_Badnell_rate_coef[nelem][ion] ,
01224 ionbal.DR_Badnell_rate_coef_mean_ion[ion]);
01225 }
01226 else if( ionbal.nDR_S_guess==1 )
01227 {
01228
01229 ionbal.DR_Badnell_rate_coef[nelem][ion] =
01230 ionbal.DR_Badnell_rate_coef[nelem][ion];
01231 }
01232 else if( ionbal.nDR_S_guess==2 )
01233 {
01234
01235 ionbal.DR_Badnell_rate_coef[nelem][ion] =
01236 ionbal.DR_Badnell_rate_coef[ipOXYGEN][ion]*
01237 ionbal.DR_S_scale[ion];
01238 }
01239 else
01240 TotalInsanity();
01241 }
01242 }
01243
01244
01245 if( ionbal.lgRecom_Badnell_print )
01246 {
01247 fprintf(ioQQQ,"DEBUG Badnell recombination RR, then DR, T=%.3e\n", phycon.te );
01248 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
01249 {
01250 fprintf(ioQQQ,"nelem=%li %s, RR then DR\n",
01251 nelem , elementnames.chElementNameShort[nelem] );
01252 for( ion=0; ion<nelem+1; ++ion )
01253 {
01254 if( ionbal.lgRR_Badnell_rate_coef_exist[nelem][ion] )
01255 {
01256 fprintf(ioQQQ," %.2e", ionbal.RR_Badnell_rate_coef[nelem][ion] );
01257 }
01258 else
01259 {
01260 fprintf(ioQQQ," %.2e", -1. );
01261 }
01262 }
01263 fprintf(ioQQQ,"\n" );
01264 for( ion=0; ion<nelem+1; ++ion )
01265 {
01266 if( ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] )
01267 {
01268 fprintf(ioQQQ," %.2e", ionbal.DR_Badnell_rate_coef[nelem][ion] );
01269 }
01270 else
01271 {
01272 fprintf(ioQQQ," %.2e", -1. );
01273 }
01274 }
01275 fprintf(ioQQQ,"\n\n" );
01276 }
01277
01278 fprintf(ioQQQ,"mean DR recombination ion mean stddev stddev/mean \n" );
01279 for( ion=0; ion<LIMELM; ++ion )
01280 {
01281 double stddev;
01282 stddev = 0.;
01283 for( nelem=ion; nelem<LIMELM; ++nelem )
01284 {
01285 stddev += POW2(
01286 ionbal.DR_Badnell_rate_coef[nelem][ion]-
01287 ionbal.DR_Badnell_rate_coef_mean_ion[ion] );
01288 }
01289 stddev = sqrt( stddev / MAX2( 1 , nsumrec[ion]-1 ) );
01290 fprintf(ioQQQ," %2li %.2e %.2e %.2e\n",
01291 ion ,
01292 ionbal.DR_Badnell_rate_coef_mean_ion[ion] ,
01293 stddev ,
01294 stddev / SDIV(ionbal.DR_Badnell_rate_coef_mean_ion[ion]) );
01295 }
01296 cdEXIT( EXIT_SUCCESS );
01297 }
01298 return;
01299 }