/home66/gary/public_html/cloudy/c08_branch/source/ion_recomb_Badnell.cpp

Go to the documentation of this file.
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 }
Generated on Mon Feb 16 12:01:17 2009 for cloudy by  doxygen 1.4.7