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