cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ion_recomb_Badnell.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*ion_recom_calculate calculate radiative and dielectronic recombination rate coefficients */
4 /*Badnell_rec_init This code is written by Terry Yun, 2005 *
5  * It reads rate coefficient fits into 3D arrays and output array.out for testing *
6  * The testing can be commented out */
7 /*Badnell_DR_rate_eval This code is written by Terry Yun, 2005 *
8  * It interpolates the rate coefficients in a given temperature.*
9  It receives ATOMIC_NUM_BIG, NELECTRONS values, temperature and returns the rate coefficient*
10  It returns
11  '-2': initial <= final
12  init < 0 or init >302 or final < 0 or final > 302
13  '-1': the transition is not defined
14  '99': unknown invalid entries */
15 #include "cddefines.h"
16 #include "phycon.h"
17 #include "elementnames.h"
18 #include "atmdat.h"
19 #include "iso.h"
20 #include "ionbal.h"
21 #include "dense.h"
22 
23 static const int MAX_FIT_PAR_DR = 9;
24 static double ***DRFitParPart1;
25 static double ***DRFitParPart2;
26 static int **nDRFitPar;
27 
28 static const int MAX_FIT_PAR_RR = 6;
29 static double ***RRFitPar;
30 static long int *nsumrec;
31 
32 /* flags to recall that we have read the fits from the main data files */
33 static bool **lgDRBadnellDefined ,
36 static bool lgMustMallocRec=true;
37 
38 /* these enable certain debugging print statements */
39 /* #define PRINT_DR */
40 /* #define PRINT_RR */
41 
42 #if defined(PRINT_DR) || defined(PRINT_RR)
43 static const char FILE_NAME_OUT[] = "array.out";
44 #endif
45 
60  /* atomic number on C scale - He - 1 */
61  int nAtomicNumberCScale,
62  /* number of core electrons before capture of free electron */
63  int n_core_e_before_recomb )
64 {
65 
66  double RateCoefficient, sum;
67  int i;
68 
69  DEBUG_ENTRY( "Badnell_DR_rate_eval()" );
70  ASSERT( nAtomicNumberCScale>=0 && nAtomicNumberCScale<LIMELM );
71 
72  if( nAtomicNumberCScale==ipIRON && n_core_e_before_recomb>=12 &&
73  n_core_e_before_recomb<=18 )
74  {
75  /* these data are from
76  *>>refer Fe DR Badnell, N., 2006, ApJ, 651, L73
77  * Fe 8+ to Fe 12+, but also include Fe13+ and Fe 14+,
78  * so these are 26-8=18 to 26-14=12
79  * increasing number of bound electrons, 0 is 14 elec, 1 is 15 elec
80  * Fe 3p^q, q=2-6
81  * this is DR fit coefficients given in table 1 of Badnell 06 */
82  double cFe_q[7][8] =
83  {
84  {5.636e-4, 7.390e-3, 3.635e-2, 1.693e-1, 3.315e-2, 2.288e-1, 7.316e-2, 0.},
85  {1.090e-3, 7.801e-3, 1.132e-2, 4.740e-2, 1.990e-1, 3.379e-2, 1.140e-1, 1.250e-1},
86  {3.266e-3, 7.637e-3, 1.005e-2, 2.527e-2, 6.389e-2, 1.564e-1, 0., 0.},
87  {1.074e-3, 6.080e-3, 1.887e-2, 2.540e-2, 7.580e-2, 2.773e-1, 0., 0.},
88  {9.073e-4, 3.777e-3, 1.027e-2, 3.321e-2, 8.529e-2, 2.778e-1, 0., 0.},
89  {5.335e-4, 1.827e-3, 4.851e-3, 2.710e-2, 8.226e-2, 3.147e-1, 0., 0.},
90  {7.421e-4, 2.526e-3, 4.605e-3, 1.489e-2, 5.891e-2, 2.318e-1, 0., 0.}
91  };
92 
93  /* Table 2 of Badnell 06 */
94  double EFe_q[7][8] =
95  {
96  {3.628e3, 2.432e4, 1.226e5, 4.351e5, 1.411e6, 6.589e6, 1.030e7, 0},
97  {1.246e3, 1.063e4, 4.719e4, 1.952e5, 5.637e5, 2.248e6, 7.202e6, 3.999e9},
98  {1.242e3, 1.001e4, 4.466e4, 1.497e5, 3.919e5, 6.853e5, 0. , 0.},
99  {1.387e3, 1.048e4, 3.955e4, 1.461e5, 4.010e5, 7.208e5, 0. , 0.},
100  {1.525e3, 1.071e4, 4.033e4, 1.564e5, 4.196e5, 7.580e5, 0. , 0.},
101  {2.032e3, 1.018e4, 4.638e4, 1.698e5, 4.499e5, 7.880e5, 0. , 0.},
102  {3.468e3, 1.353e4, 3.690e4, 1.957e5, 4.630e5, 8.202e5, 0. , 0.}
103  };
104  /* nion is for the above block of numbers */
105  long int nion = n_core_e_before_recomb - 12;
106  ASSERT( nion>=0 && nion <=6 );
107 
108  sum = 0;
109  i = 0;
110  /* loop over all non-zero terms */
111  for(i=0; i<8; ++i )
112  {
113  sum += (cFe_q[nion][i] * sexp( EFe_q[nion][i]/phycon.te));
114  }
115 
116  /*RateCoefficient = pow(phycon.te, -1.5) * sum;*/
117  RateCoefficient = sum / phycon.te32;
118 
119  return RateCoefficient;
120  }
121 
122  /*Invalid entries returns '-2':more electrons than protons */
123  else if( nAtomicNumberCScale < n_core_e_before_recomb )
124  {
125  RateCoefficient = -2;
126  }
127  /*Invalid entries returns '-2' if nAtomicNumberCScale and n_core_e_before_recomb are out of the range*/
128  else if( nAtomicNumberCScale >= LIMELM )
129  {
130  RateCoefficient = -2;
131  }
132  /*undefined z and n returns '-1'*/
133  else if( !lgDRBadnellDefined[nAtomicNumberCScale][n_core_e_before_recomb] )
134  {
135  RateCoefficient = -1;
136  }
137  else if( lgDRBadnellDefined[nAtomicNumberCScale][n_core_e_before_recomb] )
138  {
139  /* this branch, recombination coefficient has been defined */
140  sum = 0;
141  i = 0;
142  /* loop over all non-zero terms */
143  for(i=0; i<nDRFitPar[nAtomicNumberCScale][n_core_e_before_recomb]; ++i )
144  {
145  sum += (DRFitParPart1[nAtomicNumberCScale][n_core_e_before_recomb][i] *
146  sexp( DRFitParPart2[nAtomicNumberCScale][n_core_e_before_recomb][i]/phycon.te));
147  }
148 
149  /*RateCoefficient = pow(phycon.te, -1.5) * sum;*/
150  RateCoefficient = sum / phycon.te32;
151  }
152  /*unknown invalid entries returns '-99'*/
153  else
154  {
155  RateCoefficient = -99;
156  }
157 
158  return RateCoefficient;
159 }
160 
166  /* atomic number on C scale - He - 1 */
167  int nAtomicNumberCScale,
168  /* number of core electrons before capture of free electron */
169  int n_core_e_before_recomb )
170 {
171  double RateCoefficient;
172  double B, D, F;
173 
174  DEBUG_ENTRY( "Badnell_RR_rate_eval()" );
175 
176  ASSERT( nAtomicNumberCScale>=0 && nAtomicNumberCScale<LIMELM );
177 
178  if( nAtomicNumberCScale==ipIRON &&
179  n_core_e_before_recomb>=12 && n_core_e_before_recomb<=18 )
180  {
181  /* RR rate coefficients from Table 3 of
182  *>>refer Fe RR Badnell, N. 2006, ApJ, 651, L73
183  * Fe 8+ to Fe 12+, but also include Fe13+ and Fe 14+,
184  * so these are 26-8=18 to 26-14=12
185  * increasing number of bound electrons, 0 is 14 elec, 1 is 15 elec
186  * Fe 3p^q, q=2-6
187  * this is DR fit coefficients given in table 1 of Badnell 06 */
188  double parFeq[7][6] ={
189  {1.179e-9 , 0.7096, 4.508e2, 3.393e7, 0.0154, 3.977e6},
190  {1.050e-9 , 0.6939, 4.568e2, 3.987e7, 0.0066, 5.451e5},
191  {9.832e-10, 0.7146, 3.597e2, 3.808e7, 0.0045, 3.952e5},
192  {8.303e-10, 0.7156, 3.531e2, 3.554e7, 0.0132, 2.951e5},
193  {1.052e-9 , 0.7370, 1.639e2, 2.924e7, 0.0224, 4.291e5},
194  {1.338e-9 , 0.7495, 7.242e1, 2.453e7, 0.0404, 4.199e5},
195  {1.263e-9 , 0.7532, 5.209e1, 2.169e7, 0.0421, 2.917e5}
196  };
197 
198  double temp;
199  /* nion is for the above block of numbers */
200  long int nion = n_core_e_before_recomb - 12;
201  ASSERT( nion>=0 && nion <=6 );
202 
203  temp = -parFeq[nion][5]/phycon.te; /* temp = (-T2/T) */
204  B = parFeq[nion][1] + parFeq[nion][4]*exp(temp);
205  D = sqrt(phycon.te/parFeq[nion][2]); /* D = (T/T0)^1/2 */
206  F = sqrt(phycon.te/parFeq[nion][3]); /* F = (T/T1)^1/2 */
207  RateCoefficient = parFeq[nion][0]/(D*pow((1.+D),(1.-B))*pow((1.+F),(1.+B)));
208 
209  return RateCoefficient;
210  }
211 
212  /*Invalid entries returns '-2':if the z_values are smaller than equal to the n_values */
213  else if( nAtomicNumberCScale < n_core_e_before_recomb )
214  {
215  RateCoefficient = -2;
216  }
217  /*Invalid entries returns '-2' if nAtomicNumberCScale and n_core_e_before_recomb are out of the range*/
218  else if( nAtomicNumberCScale >= LIMELM )
219  {
220  RateCoefficient = -2;
221  }
222  /*undefined z and n returns '-1'*/
223  else if( !lgRRBadnellDefined[nAtomicNumberCScale][n_core_e_before_recomb] )
224  {
225  RateCoefficient = -1;
226  }
227  /* coefficients:A=RRFitPar[0], B=RRFitPar[1], T0=RRFitPar[2], T1=RRFitPar[3], DRFitParPart1=RRFitPar[4], T2=RRFitPar[5] */
228  else if( lgRRBadnellDefined[nAtomicNumberCScale][n_core_e_before_recomb] )
229  {
230 
231  /* RateCoefficient=A*[(T/T0)^1/2*(1+(T/T0)^1/2)^1-B*(1+(T/T1)^1/2)^1+B]^-1
232  where B = B + DRFitParPart1*exp(-T2/T) */
233  double temp;
234 
235  temp = -RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][5]/phycon.te; /* temp = (-T2/T) */
236  B = RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][1] +
237  RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][4]*exp(temp);
238  D = sqrt(phycon.te/RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][2]); /* D = (T/T0)^1/2 */
239  F = sqrt(phycon.te/RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][3]); /* F = (T/T1)^1/2 */
240  RateCoefficient = RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][0]/(D*pow((1.+D),(1.-B))*pow((1.+F),(1.+B)));
241  }
242 
243  /*unknown invalid entries returns '-99'*/
244  else
245  RateCoefficient = -99;
246 
247  return RateCoefficient;
248 }
249 
250 
251 /*Badnell_rec_init This code is written by Terry Yun, 2005 *
252  * It reads rate coefficient fits into 3D arrays and output array.out for testing *
253  * The testing can be commented out */
254 void Badnell_rec_init( void )
255 {
256 
257  double par_C[MAX_FIT_PAR_DR];
258  double par_E[MAX_FIT_PAR_DR];
259  char chLine[INPUT_LINE_LENGTH];
260  int NuclearCharge=-1, NumberElectrons=-1;
261  int i, k;
262  int count, number;
263  double temp_par[MAX_FIT_PAR_RR];
264  int M_state, W_state,
265  nelem;
266 
267  const int NBLOCK = 2;
268  int data_begin_line[NBLOCK];/*it tells you where the data set begins(begins with 'Z')*/
269  int length_of_line; /*this variable checks for a blank line*/
270  FILE *ioDATA;
271  const char* chFilename;
272  int yr, mo, dy;
273  char *chs;
274 
275 #define BIGGEST_INDEX_TO_USE 103
276 
277  /* Declaration of data file name array - done by Kausalya */
278  long TheirIndexToOurIndex[BIGGEST_INDEX_TO_USE];
279  char string[120];
280  double value;
281  bool lgEOL;
282  long int i1, ipISO = ipHE_LIKE;
283  long INDX=0,INDP=0,N=0,S=0,L=0,J=0,maxINDX=0,loopindex=0,max_N_of_data=-1;
284  bool lgFlag = true;
285 
286  static int nCalled = 0;
287 
288  const char* cdDATAFILE[] =
289  {
290  /* the list of filenames for Badnell DR, one to two electron */
291  "",
292  "nrb00#h_he1ic12.dat",
293  "nrb00#h_li2ic12.dat",
294  "nrb00#h_be3ic12.dat",
295  "nrb00#h_b4ic12.dat",
296  "nrb00#h_c5ic12.dat",
297  "nrb00#h_n6ic12.dat",
298  "nrb00#h_o7ic12.dat",
299  "nrb00#h_f8ic12.dat",
300  "nrb00#h_ne9ic12.dat",
301  "nrb00#h_na10ic12.dat",
302  "nrb00#h_mg11ic12.dat",
303  "nrb00#h_al12ic12.dat",
304  "nrb00#h_si13ic12.dat",
305  "nrb00#h_p14ic12.dat",
306  "nrb00#h_s15ic12.dat",
307  "nrb00#h_cl16ic12.dat",
308  "nrb00#h_ar17ic12.dat",
309  "nrb00#h_k18ic12.dat",
310  "nrb00#h_ca19ic12.dat",
311  "nrb00#h_sc20ic12.dat",
312  "nrb00#h_ti21ic12.dat",
313  "nrb00#h_v22ic12.dat",
314  "nrb00#h_cr23ic12.dat",
315  "nrb00#h_mn24ic12.dat",
316  "nrb00#h_fe25ic12.dat",
317  "nrb00#h_co26ic12.dat",
318  "nrb00#h_ni27ic12.dat",
319  "nrb00#h_cu28ic12.dat",
320  "nrb00#h_zn29ic12.dat"
321  };
322  //End of modification
323 
324  DEBUG_ENTRY( "Badnell_rec_init()" );
325 
326  /* must only do this once */
327  if( nCalled > 0 )
328  {
329  return;
330  }
331  ++nCalled;
332 
333 # if defined(PRINT_DR) || defined(PRINT_RR)
334  FILE *ofp = open_data( FILE_NAME_OUT, "w", AS_LOCAL_ONLY );
335 # endif
336 
337  /* Modification done by Kausalya
338  * Start - Try to open all the 29 data files.*/
339  for(nelem=ipHELIUM;nelem<LIMELM;nelem++)
340  {
341  if( nelem < 2 || dense.lgElmtOn[nelem] )
342  {
343  ioDATA= open_data( cdDATAFILE[nelem], "r" );
344 
345  lgFlag = true;
346  ASSERT(ioDATA);
347 
348  for( i=0; i<BIGGEST_INDEX_TO_USE; i++ )
349  TheirIndexToOurIndex[i] = -1;
350 
351  /* Reading lines */
352  while(lgFlag)
353  {
354  if(read_whole_line(string,sizeof(string),ioDATA)!=NULL)
355  {
356  if( nMatch("INDX INDP ",string) )
357  {
358  /* ignore next line of data */
359  if( read_whole_line( string , (int)sizeof(string) , ioDATA ) == NULL )
360  {
361  fprintf( ioQQQ, " Badnell data file appears to be corrupted.\n");
362  cdEXIT(EXIT_FAILURE);
363  }
364 
365  /* This one should be real data */
366  while( read_whole_line(string, (int)sizeof(string), ioDATA) != NULL )
367  {
368  if( strcmp(string,"\n")==0 )
369  {
370  lgFlag = false;
371  break;
372  }
373 
374  i1=3;
375  INDX=(long)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
376  if( INDX >= BIGGEST_INDEX_TO_USE )
377  {
378  INDX--;
379  lgFlag = false;
380  break;
381  }
382 
383  ASSERT( INDX < BIGGEST_INDEX_TO_USE );
384 
385  INDP=(long)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
386  ASSERT( INDP >= 1 );
387 
388  if(INDP==1)
389  {
390  if( (i1=nMatch("1S1 ",string)) > 0 )
391  {
392  i1 += 4;
393  N=(long)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
394  ASSERT( N>=1 );
395  }
396  else
397  {
398  TotalInsanity();
399  }
400 
401  if( (i1=nMatch(" (",string)) > 0 )
402  {
403  i1 += 6;
404  S=(long)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
405  /* S in file is 3 or 1, we need 1 or 0 */
406  ASSERT( S==1 || S==3 );
407  }
408  else
409  {
410  TotalInsanity();
411  }
412 
413  /* move i1 one further to get L */
414  i1++;
415  L=(long)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
416  ASSERT( L >= 0 && L < N );
417 
418  /* move i1 two further to get J */
419  i1 += 2;
420  J=(long)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
421  ASSERT( J <= ( L + (int)((S+1)/2) ) &&
422  J >= ( L - (int)((S+1)/2) ) && J >= 0 );
423 
424  /* if line in data file is higher N than highest considered, stop reading. */
426  TheirIndexToOurIndex[INDX] = iso.QuantumNumbers2Index[ipHE_LIKE][nelem][N][L][S];
427  else
428  {
429  /* Current line is not being used,
430  * decrement INDX so maxINDX is set correctly below. */
431  INDX--;
432  lgFlag = false;
433  break;
434  }
435 
436  /* Must adjust index if in 2^3Pj term */
437  if( N==2 && L==1 && S==3 )
438  {
439  if( J==0 )
440  TheirIndexToOurIndex[INDX] = 3;
441  else if( J==1 )
442  TheirIndexToOurIndex[INDX] = 4;
443  else
444  {
445  ASSERT( J==2 );
446  ASSERT( TheirIndexToOurIndex[INDX] == 5 );
447  }
448  }
449  max_N_of_data = MAX2( max_N_of_data, N );
450  }
451  else
452  {
453  // Stop parsing the tuple since INDP!=1
454  lgFlag = false;
455  }
456  }
457  }
458  }
459  else
460  {
461  // End of file is reached.
462  lgFlag = false;
463  }
464  }
465 
466  maxINDX =INDX;
467  ASSERT( maxINDX > 0 );
468  ASSERT( maxINDX < BIGGEST_INDEX_TO_USE );
469  /* reset INDX */
470  INDX = 0;
471  lgFlag = true;
472  while(lgFlag)
473  {
474  if(read_whole_line(string,sizeof(string),ioDATA)!=NULL)
475  {
476  /* to access the first table whose columns are INDX ,INDP */
477  if( nMatch("INDX TE= ",string) )
478  {
479  lgFlag = false;
480  /* we found the beginning of the data array */
481  /* ignore next line of data */
482  if( read_whole_line( string , (int)sizeof(string) , ioDATA ) == NULL )
483  {
484  fprintf( ioQQQ, " Badnell data file appears to be corrupted.\n");
485  cdEXIT(EXIT_FAILURE);
486  }
487 
488  /* This one should be real data */
489  while( read_whole_line(string, (int)sizeof(string), ioDATA) != NULL )
490  {
491  /* If we find this string, we have reached the end of the table. */
492  if( nMatch("PRTF",string) || INDX >= maxINDX || INDX<0 )
493  break;
494 
495  i1=3;
496  INDX=(long)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
497  if( INDX>maxINDX )
498  break;
499 
500  for(loopindex=0;loopindex<10;loopindex++)
501  {
502  value=(double)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
503  if( TheirIndexToOurIndex[INDX] < iso.numLevels_max[ipHE_LIKE][nelem] &&
504  TheirIndexToOurIndex[INDX] > 0 )
505  {
506  iso.DielecRecombVsTemp[ipHE_LIKE][nelem][TheirIndexToOurIndex[INDX]][loopindex] += value;
507  }
508  }
509 
510  /* data are broken into two lines, read second line here */
511  if( read_whole_line( string , (int)sizeof(string) , ioDATA ) == NULL )
512  {
513  fprintf( ioQQQ, " Badnell data file appears to be corrupted.\n");
514  cdEXIT(EXIT_FAILURE);
515  }
516 
517  /* start of data for second line */
518  i1 = 13;
519  for(loopindex=10;loopindex<19;loopindex++)
520  {
521  value=(double)FFmtRead(string,&i1,INPUT_LINE_LENGTH,&lgEOL);
522  if( TheirIndexToOurIndex[INDX] < iso.numLevels_max[ipHE_LIKE][nelem] &&
523  TheirIndexToOurIndex[INDX] > 0 )
524  {
525  iso.DielecRecombVsTemp[ipHE_LIKE][nelem][TheirIndexToOurIndex[INDX]][loopindex] += value;
526  }
527  }
528  }
529  }
530  }
531  else
532  lgFlag = false;
533  }
534  fclose(ioDATA);
535  ASSERT( maxINDX > 0 );
536  ASSERT( maxINDX < BIGGEST_INDEX_TO_USE );
537  ASSERT( max_N_of_data > 0 );
538 
539  if( max_N_of_data < iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] )
540  {
541  long indexOfMaxN;
542  L = -1;
543  S = -1;
544 
545  /* This loop extrapolates nLS data to nLS states */
546  for( i=TheirIndexToOurIndex[maxINDX]+1;
547  i<iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem]; i++ )
548  {
549  L = L_(i);
550  S = S_(i);
551 
552  if( L > 4 )
553  continue;
554 
555  indexOfMaxN = iso.QuantumNumbers2Index[ipHE_LIKE][nelem][max_N_of_data][L][S];
556  for(loopindex=0;loopindex<19;loopindex++)
557  {
558  iso.DielecRecombVsTemp[ipHE_LIKE][nelem][i][loopindex] =
559  iso.DielecRecombVsTemp[ipHE_LIKE][nelem][indexOfMaxN][loopindex] *
560  pow( (double)max_N_of_data/(double)StatesElem[ipHE_LIKE][nelem][i].n, 3.);
561  }
562  }
563 
564  /* Get the N of the highest resolved singlet P (in the model, not the data) */
565  indexOfMaxN =
567 
568  /* This loop extrapolates nLS data to collapsed n levels, just use highest singlet P data */
569  for( i=iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem];
570  i<iso.numLevels_max[ipHE_LIKE][nelem]; i++ )
571  {
572  for(loopindex=0;loopindex<19;loopindex++)
573  {
574  iso.DielecRecombVsTemp[ipHE_LIKE][nelem][i][loopindex] =
575  iso.DielecRecombVsTemp[ipHE_LIKE][nelem][indexOfMaxN][loopindex] *
576  pow( (double)iso.n_HighestResolved_max[ipHE_LIKE][nelem]/
577  (double)StatesElem[ipHE_LIKE][nelem][i].n, 3.);
578  }
579  }
580  }
581  }
582  }
583 
584  for( i=0; i<NBLOCK; ++i )
585  {
586  /* set to really large negative number - crash if used before being redefined */
587  data_begin_line[i] = INT_MIN;
588  }
589 
590  chFilename = "badnell_dr.dat";
591  ioDATA = open_data( chFilename, "r" );
592 
593  count = 0;
594  number = 0;
595 
596  /*Find out the number line where the data starts
597  * there are two main blocks of data and each starts with a Z in column 2 */
598  while( read_whole_line(chLine, (int)sizeof(chLine), ioDATA) != NULL )
599  {
600  count++;
601 
602  if( chLine[2]=='Z' )
603  {
604  /* number has to be 0 or 1, and indicates the first or second block of data
605  * count is the line number for the start of that block */
606  data_begin_line[number] = count;
607  ASSERT( number < NBLOCK );
608  number++;
609  }
610  }
611 
612  /*set a flag for a undefined data*/
613  if( lgMustMallocRec )
614  {
615  nsumrec = (long *)MALLOC(LIMELM*sizeof(long) );
616  nDRFitPar = (int**)MALLOC( LIMELM*sizeof( int*) );
617  lgDRBadnellDefined = (bool **)MALLOC( LIMELM*sizeof(bool*) );
618  lgDRBadnellDefinedPart2 = (bool **)MALLOC( LIMELM*sizeof(bool*) );
619  lgRRBadnellDefined = (bool **)MALLOC( LIMELM*sizeof(bool*) );
620 
621  DRFitParPart1 = (double ***)MALLOC( LIMELM*sizeof(double**) );
622  DRFitParPart2 = (double ***)MALLOC( LIMELM*sizeof(double**) );
623  RRFitPar = (double ***)MALLOC( LIMELM*sizeof(double**) );
624  }
625 
626  for( nelem=0; nelem<LIMELM; nelem++ )
627  {
628  if( lgMustMallocRec )
629  {
630  nDRFitPar[nelem] = (int*)MALLOC( (nelem+1)*sizeof( int) );
631  lgDRBadnellDefined[nelem] = (bool *)MALLOC( (nelem+1)*sizeof(bool) );
632  lgDRBadnellDefinedPart2[nelem] = (bool *)MALLOC( (nelem+1)*sizeof(bool) );
633  lgRRBadnellDefined[nelem] = (bool *)MALLOC( (nelem+1)*sizeof(bool) );
634 
635  DRFitParPart1[nelem] = (double **)MALLOC( (nelem+1)*sizeof(double*) );
636  DRFitParPart2[nelem] = (double **)MALLOC( (nelem+1)*sizeof(double*) );
637  RRFitPar[nelem] = (double **)MALLOC( (nelem+1)*sizeof(double*) );
638  }
639  for( long ion=0; ion<nelem+1; ++ion )
640  {
641  if( lgMustMallocRec )
642  {
643  DRFitParPart1[nelem][ion] = (double *)MALLOC( MAX_FIT_PAR_DR*sizeof(double) );
644  DRFitParPart2[nelem][ion] = (double *)MALLOC( MAX_FIT_PAR_DR*sizeof(double) );
645  RRFitPar[nelem][ion] = (double *)MALLOC( MAX_FIT_PAR_RR*sizeof(double) );
646  }
647  lgDRBadnellDefined[nelem][ion] = false;
648  lgDRBadnellDefinedPart2[nelem][ion] = false;
649  lgRRBadnellDefined[nelem][ion] = false;
650 
651  /*set fitting coefficients to zero initially*/
652  for( k=0; k<MAX_FIT_PAR_DR; k++ )
653  {
654  DRFitParPart1[nelem][ion][k] = 0;
655  DRFitParPart2[nelem][ion][k] = 0;
656  }
657  for( k=0; k<MAX_FIT_PAR_RR; k++ )
658  {
659  RRFitPar[nelem][ion][k] = 0;
660  }
661  }
662  }
663  lgMustMallocRec = false;
664 
665  count = 0;
666  /*Start from beginning to read in again*/
667  fseek(ioDATA, 0, SEEK_SET);
668  /* read magic number for DR data */
669  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
670  {
671  fprintf( ioQQQ, " DISASTER PROBLEM Badnell_rec_init could not read first line of badnell_dr.dat.\n");
672  cdEXIT(EXIT_FAILURE);
673  }
674  count++;
675 
676  /* look for ')' on the line, magic number comes after it */
677  if( (chs = strchr(chLine, ')'))==NULL )
678  {
679  /* format is incorrect */
680  fprintf( ioQQQ, " DISASTER PROBLEM Badnell_rec_init data file incorrect format.\n");
681  cdEXIT(EXIT_FAILURE);
682  }
683 
684  ++chs;
685  sscanf(chs, "%4i%2i%2i",&yr, &mo, &dy);
686  /* check magic number - the date on the line */
687  int dr_yr = 2007, dr_mo = 10, dr_dy = 27;
688  if((yr != dr_yr) || (mo != dr_mo) || (dy != dr_dy))
689  {
690  fprintf(ioQQQ,
691  "DISASTER PROBLEM Badnell_rec_init The version of %s I found (%i %i %i) is not the current version (%i %i %i).\n",
692  chFilename, yr, mo, dy, dr_yr, dr_mo, dr_dy);
693  fprintf(ioQQQ," The first line of the file is the following\n %s\n", chLine );
694  cdEXIT(EXIT_FAILURE);
695  }
696 
697  while( read_whole_line(chLine, (int)sizeof(chLine), ioDATA) != NULL )
698  {
699  count++;
700  length_of_line = (int)strlen(chLine);
701 
702  /*reading in coefficients DRFitParPart1 */
703  if( count > data_begin_line[0] && count < data_begin_line[1] && length_of_line >3 )
704  {
705  /*set array par_C to zero */
706  for( i=0; i<MAX_FIT_PAR_DR; i++ )
707  {
708  par_C[i] = 0;
709  }
710  sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf%lf%lf%lf",
711  &NuclearCharge, &NumberElectrons, &M_state, &W_state, &par_C[0], &par_C[1], &par_C[2],
712  &par_C[3], &par_C[4], &par_C[5], &par_C[6], &par_C[7], &par_C[8]);
713  /* data files have atomic number on physics scale, convert to C scale
714  * for following code */
715  long int NuclearChargeM1 = NuclearCharge-1;
716 
717  if(M_state == 1 && NuclearChargeM1 < LIMELM )
718  {
719  /*Set a flag to '1' when the indices are defined */
720  ASSERT( NumberElectrons < LIMELM );
721  ASSERT( NuclearChargeM1 < LIMELM );
722  lgDRBadnellDefined[NuclearChargeM1][NumberElectrons] = true;
723 
724  /*counting the number of coefficients */
725  nDRFitPar[NuclearChargeM1][NumberElectrons] = 9;
726  for( i=8; i>=0; i-- )
727  {
728  if( par_C[i] == 0 )
729  --nDRFitPar[NuclearChargeM1][NumberElectrons];
730  else
731  break;
732  }
733 
734  /*assign the values into array */
735  for( i=0; i<9; i++ )
736  DRFitParPart1[NuclearChargeM1][NumberElectrons][i] = par_C[i];
737  }
738  }
739  }
740 
741  /*starting again to read in E values */
742  fseek(ioDATA, 0, SEEK_SET);
743  count = 0;
744  while( read_whole_line(chLine, (int)sizeof(chLine), ioDATA) != NULL )
745  {
746  count++;
747  length_of_line = (int)strlen(chLine);
748  if( count > data_begin_line[1] && length_of_line > 3 )
749  {
750 
751  /*set array par_E to zero*/
752  for( i=0; i<MAX_FIT_PAR_DR; i++ )
753  {
754  par_E[i] = 0;
755  }
756  sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf%lf%lf%lf",
757  &NuclearCharge, &NumberElectrons, &M_state, &W_state, &par_E[0], &par_E[1], &par_E[2],
758  &par_E[3], &par_E[4], &par_E[5], &par_E[6], &par_E[7], &par_E[8]);
759  /* data file is on physics scale but we use C scale */
760  long int NuclearChargeM1 = NuclearCharge-1;
761 
762  if(M_state == 1 && NuclearChargeM1<LIMELM)
763  {
764  ASSERT( NumberElectrons < LIMELM );
765  ASSERT( NuclearChargeM1 < LIMELM );
766  lgDRBadnellDefinedPart2[NuclearChargeM1][NumberElectrons] = true;
767 
768  /*counting the number of coefficients */
769  nDRFitPar[NuclearChargeM1][NumberElectrons] = 9;
770  for( i=8; i>=0; i-- )
771  {
772  if( par_E[i] == 0 )
773  --nDRFitPar[NuclearChargeM1][NumberElectrons];
774  else
775  break;
776  }
777 
778  /*assign the values into array*/
779  for( i=0; i<nDRFitPar[NuclearChargeM1][NumberElectrons]; i++ )
780  DRFitParPart2[NuclearChargeM1][NumberElectrons][i] = par_E[i];
781  }
782  }
783  }
784 
785  fclose( ioDATA );
786 
787  /*output coefficients for defined values for testing */
788 # ifdef PRINT_DR
789  for( nelem=0; nelem<LIMELM; nelem++ )
790  {
791  for( int ion=0; ion<nelem+1;++ion )
792  {
793  if( lgDRBadnellDefined[nelem][ion] )
794  {
795  fprintf(ofp, "%i %i %e %e %e %e %e %e %e %e %e\n",
796  nelem, ion, DRFitParPart1[nelem][ion][0],
797  DRFitParPart1[nelem][ion][1], DRFitParPart1[nelem][ion][2],
798  DRFitParPart1[nelem][ion][3], DRFitParPart1[nelem][ion][4],
799  DRFitParPart1[nelem][ion][5], DRFitParPart1[nelem][ion][6],
800  DRFitParPart1[nelem][ion][7], DRFitParPart1[nelem][ion][8]);
801  }
802  }
803  }
804  for( nelem=0; nelem<LIMELM; nelem++ )
805  {
806  for( int ion=0; ion<nelem+1; ion++ )
807  {
808  if( lgDRBadnellDefinedPart2[nelem][ion] )
809  {
810  fprintf(ofp, "%i %i %e %e %e %e %e %e %e %e %e\n",
811  nelem, ion, DRFitParPart2[nelem][ion][0],
812  DRFitParPart2[nelem][ion][1], DRFitParPart2[nelem][ion][2],
813  DRFitParPart2[nelem][ion][3], DRFitParPart2[nelem][ion][4],
814  DRFitParPart2[nelem][ion][5], DRFitParPart2[nelem][ion][6],
815  DRFitParPart2[nelem][ion][7], DRFitParPart2[nelem][ion][8]);
816  }
817  }
818  }
819  fclose(ofp);
820 # endif
821 
822  /*checking for the match of lgDRBadnellDefined and lgDRBadnellDefinedPart2 -
823  * Both have to be defined*/
824  bool lgDRBadnellBothDefined = true;
825  for( nelem=0; nelem<LIMELM; nelem++ )
826  {
827  for( int ion=0; ion<nelem+1; ion++ )
828  {
829  /* check that first and second half of DR fitting coefficients
830  * are both defined */
831  if( lgDRBadnellDefined[nelem][ion] != lgDRBadnellDefinedPart2[nelem][ion] )
832  {
833  fprintf( ioQQQ, "DR %i, RR %i: %c %c\n", nelem, ion,
834  TorF(lgDRBadnellDefined[nelem][ion]),
835  TorF(lgDRBadnellDefinedPart2[nelem][ion]));
836  fprintf( ioQQQ, "PROBLEM ion_recomb_Badnell first and second half of Badnell DR not consistent.\n");
837  lgDRBadnellBothDefined = false;
838  }
839  }
840  }
841 
842  if( !lgDRBadnellBothDefined )
843  {
844  /* disaster - DR files are not consistent */
845  fprintf(ioQQQ,
846  "DISASTER PROBLEM The DR data files are corrupted - part 1 and 2 do not agree.\n");
847  fprintf(ioQQQ," Start again with a fresh copy of the data directory\n" );
848  cdEXIT(EXIT_FAILURE);
849  }
850 
851  /* now do radiative recombination */
852  chFilename = "badnell_rr.dat";
853  ioDATA = open_data( chFilename, "r" );
854 
855  /* read magic number for RR data */
856  {
857  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
858  {
859  fprintf( ioQQQ, " DISASTER PROBLEM Badnell_rec_init could not read first line of badnell_rr.dat.\n");
860  cdEXIT(EXIT_FAILURE);
861  }
862  /* this is just before date, which we use for magic number */
863  if( (chs = strchr(chLine, ')'))==NULL )
864  {
865  /* format is incorrect */
866  fprintf( ioQQQ, " DISASTER PROBLEM Badnell_rec_init data file incorrect format.\n");
867  cdEXIT(EXIT_FAILURE);
868  }
869  ++chs;
870  sscanf(chs, "%4i%2i%2i", &yr, &mo, &dy);
871  int rr_yr = 2007, rr_mo = 10, rr_dy = 27;
872  if((yr != rr_yr)||(mo != rr_mo)||(dy != rr_dy))
873  {
874  fprintf(ioQQQ,"DISASTER PROBLEM The version of %s I found (%i %i %i) is not the current version (%i %i %i).\n",
875  chFilename, yr, mo, dy, rr_yr, rr_mo, rr_dy);
876  fprintf(ioQQQ," The line was as follows:\n %s\n", chLine );
877  cdEXIT(EXIT_FAILURE);
878  }
879  }
880 
881  while( read_whole_line(chLine, (int)sizeof(chLine), ioDATA) != NULL )
882  {
883  /*read in coefficients - first set array par to zero */
884  for( i=0; i<MAX_FIT_PAR_RR; i++ )
885  {
886  temp_par[i] = 0;
887  }
888  if(chLine[0] != '#')
889  {
890  sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf",
891  &NuclearCharge, &NumberElectrons, &M_state, &W_state, &temp_par[0], &temp_par[1],
892  &temp_par[2], &temp_par[3], &temp_par[4], &temp_par[5]);
893  long NuclearChargeM1 = NuclearCharge-1;
894 
895  if(M_state == 1 && NuclearChargeM1<LIMELM)
896  {
897  ASSERT( NuclearChargeM1 < LIMELM );
898  ASSERT( NumberElectrons <= LIMELM );
899  /*Set a flag to '1' when the indices are defined */
900  lgRRBadnellDefined[NuclearChargeM1][NumberElectrons] = true;
901  /*assign the values into array */
902  for( i=0; i<MAX_FIT_PAR_RR; i++ )
903  RRFitPar[NuclearChargeM1][NumberElectrons][i] = temp_par[i];
904  }
905  }
906  }
907 
908  /*output coefficients for defined values for testing */
909 # ifdef PRINT_RR
910  count = 0;
911  for( nelem=0; nelem<LIMELM; nelem++ )
912  {
913  for( ion=0; ion<nelem+1; ion++ )
914  {
915  if( lgRRBadnellDefined[nelem][ion] )
916  {
917  fprintf(ofp, "%i %i %e %e %e %e %e %e\n",
918  nelem, ion, RRFitPar[nelem][ion][0],
919  RRFitPar[nelem][ion][1], RRFitPar[nelem][ion][2],
920  RRFitPar[nelem][ion][3],
921  RRFitPar[nelem][ion][4], RRFitPar[nelem][ion][5]);
922  count++;
923  }
924  }
925  }
926  fprintf(ofp, "total lines are %i ", count);
927 
928  fclose(ofp);
929 # endif
930 
931  fclose(ioDATA);
932 
933  {
934  enum {DEBUG_LOC=false};
935  if( DEBUG_LOC )
936  {
937  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
938  {
939  int ion;
940  fprintf(ioQQQ,"\nDEBUG rr rec\t%i",nelem);
941  for( ion=0; ion<=nelem; ++ion )
942  {
943  fprintf(ioQQQ,"\t%.2e", Badnell_RR_rate_eval(nelem+1 , nelem-ion ) );
944  }
945  fprintf(ioQQQ,"\n");
946  fprintf(ioQQQ,"DEBUG dr rec\t%i",nelem);
947  for( ion=0; ion<=nelem; ++ion )
948  {
949  fprintf(ioQQQ,"\t%.2e", Badnell_DR_rate_eval(nelem+1 , nelem-ion ) );
950  }
951  fprintf(ioQQQ,"\n");
952  }
953  cdEXIT(EXIT_FAILURE);
954  }
955  }
956  return;
957 }
958 
959 /*ion_recom_calculate calculate radiative and dielectronic recombination rate coefficients */
961 {
962  static double TeUsed = -1;
963  long int ion ,
964  nelem ,
965  i;
966 
967  DEBUG_ENTRY( "ion_recom_calculate()" );
968 
969  /* do not reevaluate if change in temperature is small */
970  if( fabs(phycon.te/TeUsed - 1. ) < 0.001 )
971  {
972  return;
973  }
974 
975  /* save mean rec coefficients - used to derive rates for those ions with none */
976  for( ion=0; ion<LIMELM; ++ion )
977  {
979  nsumrec[ion] = 0;
980  }
981  TeUsed = phycon.te;
982  /* NB - for following loop important to go over all elements and all
983  * ions, not just active ones, since mean DR is used as the guess for
984  * DR rates for unknown species. */
985  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
986  {
987  for( ion=0; ion<nelem+1; ++ion )
988  {
989  long int n_bnd_elec_before_recom ,
990  n_bnd_elec_after_recom;
991 
992  n_bnd_elec_before_recom = nelem-ion;
993  n_bnd_elec_after_recom = nelem-ion+1;
994 # define DR2SMALL 1e-15
995  /* Badnell dielectronic recombination rate coefficients */
996  if( (ionbal.DR_Badnell_rate_coef[nelem][ion] =
998  /* atomic number on C scale */
999  nelem,
1000  /* number of core electrons before capture of free electron,
1001  * for bare ion this is zero */
1002  n_bnd_elec_before_recom )) >= 0. )
1003  {
1004  ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = true;
1005  if( ionbal.DR_Badnell_rate_coef[nelem][ion] > DR2SMALL)
1006  {
1007  /* keep track of mean DR rate for this ion */
1008  ++nsumrec[ion];
1010  log10(ionbal.DR_Badnell_rate_coef[nelem][ion]);
1011  }
1012  }
1013  else
1014  {
1015  /* real rate does not exist, will use mean later */
1016  ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = false;
1017  ionbal.DR_Badnell_rate_coef[nelem][ion] = 0.;
1018  }
1019 
1020  /* Badnell radiative recombination rate coefficients */
1022  /* atomic number on C scale */
1023  nelem,
1024  /* number of core electrons before capture of free electron */
1025  n_bnd_elec_before_recom )) >= 0. )
1026  {
1027  ionbal.lgRR_Badnell_rate_coef_exist[nelem][ion] = true;
1028  }
1029  else
1030  {
1031  /* real rate does not exist, will use mean later */
1032  ionbal.lgRR_Badnell_rate_coef_exist[nelem][ion] = false;
1033  ionbal.RR_Badnell_rate_coef[nelem][ion] = 0.;
1034  }
1035 
1036  /* save D. Verner's radiative recombination rate coefficient
1037  * needed for rec cooling, cm3 s-1 */
1038  ionbal.RR_Verner_rate_coef[nelem][ion] =
1039  t_ADfA::Inst().rad_rec(
1040  /* number number of physics scale */
1041  nelem+1 ,
1042  /* number of protons on physics scale */
1043  n_bnd_elec_after_recom ,
1044  phycon.te );
1045  }
1046  }
1047 
1048  /* this branch starts idiosyncratic single ions */
1049  double Fe_Gu_c[9][6] = {
1050  { 2.50507e-11, 5.60226e-11, 1.85001e-10, 3.57495e-9, 1.66321e-7, 0. },/*fit params for Fe+6*/
1051  { 9.19610e-11, 2.92460e-10, 1.02120e-9, 1.14852e-8, 3.25418e-7, 0. }, /* fitting params for Fe+7 */
1052  { 9.02625e-11, 6.22962e-10, 5.77545e-9, 1.78847e-8, 3.40610e-7, 0. }, /* fitting params for Fe+8 */
1053  { 9.04286e-12, 9.68148e-10, 4.83636e-9, 2.48159e-8, 3.96815e-7, 0. }, /* fitting params for Fe+9 */
1054  { 6.77873e-10, 1.47252e-9, 5.31202e-9, 2.54793e-8, 3.47407e-7, 0. }, /* fitting params for Fe+10 */
1055  { 1.29742e-9, 4.10172e-9, 1.23605e-8, 2.33615e-8, 2.97261e-7, 0. }, /* fitting params for Fe+11 */
1056  { 8.78027e-10, 2.31680e-9, 3.49333e-9, 1.16927e-8, 8.18537e-8, 1.54740e-7 },/*fit params for Fe+12*/
1057  { 2.23178e-10, 1.87313e-9, 2.86171e-9, 1.38575e-8, 1.17803e-7, 1.06251e-7 },/*fit params for Fe+13*/
1058  { 2.17263e-10, 7.35929e-10, 2.81276e-9, 1.32411e-8, 1.15761e-7, 4.80389e-8 }/*fit params for Fe+14*/
1059  },
1060 
1061  Fe_Gu_E[9][6] = {
1062  { 8.30501e-2, 8.52897e-1, 3.40225e0, 2.23053e1, 6.80367e1, 0. }, /* fitting params for Fe+6 */
1063  { 1.44392e-1, 9.23999e-1, 5.45498e0, 2.04301e1, 7.06112e1, 0. }, /* fitting params for Fe+7 */
1064  { 5.79132e-2, 1.27852e0, 3.22439e0, 1.79602e1, 6.96277e1, 0. }, /* fitting params for Fe+8 */
1065  { 1.02421e-1, 1.79393e0, 4.83226e0, 1.91117e1, 6.80858e1, 0. }, /* fitting params for Fe+9 */
1066  { 1.24630e-1, 6.86045e-1, 3.09611e0, 1.44023e1, 6.42820e1, 0. }, /* fitting params for Fe+10 */
1067  { 1.34459e-1, 6.63028e-1, 2.61753e0, 1.30392e1, 6.10222e1, 0. }, /* fitting params for Fe+11 */
1068  { 7.79748e-2, 5.35522e-1, 1.88407e0, 8.38459e0, 3.38613e1, 7.89706e1 }, /*fitting params for Fe+12*/
1069  { 8.83019e-2, 6.12756e-1, 2.36035e0, 9.61736e0, 3.64467e1, 8.72406e1 }, /*fitting params for Fe+13*/
1070  { 1.51322e-1, 5.63155e-1, 2.57013e0, 9.08166e0, 3.69528e1, 1.08067e2 } /* fitting params for Fe+14*/
1071  };
1072 
1073  double s_c[5][2] = {
1074  {0.1565e-3, 0.1617e-2}, /* fitting parameters for S+1 */
1075  {0.3026e-3, 0.2076e-1}, /* fitting parameters for S+2 */
1076  {0.3177e-1, 0.6309e-3}, /* fitting parameters for S+3 */
1077  {0.2464e-1, 0.5553e-3}, /* fitting parameters for S+4 */
1078  {0.1924e-1, 0.5111e-3} /* fitting parameters for s+5 */
1079  },
1080 
1081  s_E[5][2] = {
1082  {0.1157e6, 0.1868e6}, /* fitting parameters for S+1 */
1083  {0.9662e5, 0.1998e6}, /* fitting parameters for S+2 */
1084  {0.1928e6, 0.6126e5}, /* fitting parameters for S+3 */
1085  {0.1806e6, 0.3142e5}, /* fitting parameters for S+4 */
1086  {0.1519e6, 0.4978e5} /* fitting parameters for s+5 */
1087  };
1088 
1089  /* do a series of special cases for Fe DR */
1090  nelem = ipIRON;
1091 
1092  /*the sum of the fitting parameter calculations*/
1093  double fitSum = 0;
1094 
1095  /* Fe+14 - Fe+13 - ion = 0 is A^+ -> A^0 so off by one*/
1096  ion = 13;
1097  /* >>chng Fe+14 -> Fe+13, experimental DR from
1098  * >>refer Fe+14 DR Lukic, D.V. et al 2007, astroph 0704-0905
1099  * following is the MCBP entry from Table 3,
1100  * units of Fe14_c are cm^3 s^-1 K^1.5 */
1101  double fe14_c[6]={7.07E-04,7.18E-03,2.67E-02,3.15E-02,1.62E-01,5.37E-04};
1102  /* units of fe14_E are eV THIS IS DIFFERENT FROM OTHER PAPERS BY
1103  * THESE AUTHORS - THEY CHANGE TEMPERATURE UNITS WITHIN THE SAME
1104  * PARAGRAPH!!! */
1105  double fe14_E[6]={4.12E-01,2.06E+00,1.03E+01,2.20E+01,4.22E+01,3.41E+03};
1106  /* update DR rate - always do this since this reference is more
1107  * recent that previous paper */
1108  for( i=0; i<6; i++ )
1109  {
1110  fitSum += fe14_c[i] * sexp( fe14_E[i]/phycon.te_eV );
1111  }
1112 
1113  ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = true;
1114  ionbal.DR_Badnell_rate_coef[nelem][ion] = fitSum / phycon.te32;
1115  if( ionbal.DR_Badnell_rate_coef[nelem][ion]>DR2SMALL )
1116  {
1117  ++nsumrec[ion];
1119  log10(ionbal.DR_Badnell_rate_coef[nelem][ion]);
1120  }
1121 
1122  /* >>chng 06 jun 21 by Mitchell Martin, added new DR rate coefficient
1123  * for Fe XIV (Fe+13) to Fe XIII (Fe+12) recombination calculation
1124  *>>refer Fe12 DR Schmidt E.W. et al. 2006, ApJ, 641, 157L */
1125  /*fitting parameters used to calculate the DR rate for Fe+13 -> Fe+12*/
1126  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 },
1127  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 };
1128  /* do Fe+13 -> Fe+12 from Savin et al. if not already done */
1129  ion = 12;
1130  if( !ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] )
1131  {
1132  for( i=0; i<10; i++ )
1133  {
1134  fitSum += fe13_c[i] * sexp( fe13_E[i]/phycon.te_eV );
1135  }
1136 
1137  /* update DR rate is not already done */
1138  ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = true;
1139  ionbal.DR_Badnell_rate_coef[nelem][ion] = fitSum / phycon.te32;
1140  if( ionbal.DR_Badnell_rate_coef[nelem][ion] > DR2SMALL )
1141  {
1142  ++nsumrec[ion];
1144  log10(ionbal.DR_Badnell_rate_coef[nelem][ion]);
1145  }
1146  }
1147 
1148  /* this is the temperature in eV evaluated to the 3/2 power */
1149  double te_eV32 = pow( phycon.te_eV, 1.5);
1150 
1151  /* >>chng 06 jul 07 by Mitchell Martin, added DR rate coefficient
1152  * calculations for Fe+6->Fe+5 through Fe+14->Fe+13
1153  * this is still for nelem = ipIRON from the previous calculation
1154  * starts with Fe+6 -> Fe+5 and does the next ion with each iteration */
1155  for( ion=0; ion<9; ion++ )
1156  {
1157  /* only do this rate if not already done by a previous approximation */
1158  if( !ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion+5] )
1159  {
1160  fitSum = 0; /* resets the fitting parameter calculation */
1161  for( i=0; i<6; i++ )
1162  {
1163  fitSum += Fe_Gu_c[ion][i] * sexp( Fe_Gu_E[ion][i]/phycon.te_eV );
1164  }
1165  ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion+5] = true;
1166  ionbal.DR_Badnell_rate_coef[nelem][ion+5] = fitSum / te_eV32;
1167  if( ionbal.DR_Badnell_rate_coef[nelem][ion+5] > DR2SMALL )
1168  {
1169  ++nsumrec[ion+5];
1171  log10(ionbal.DR_Badnell_rate_coef[nelem][ion+5]);
1172  }
1173  }
1174  }
1175  /* this is end of Fe DR rates */
1176 
1177  /* now get mean of the DR rates - may be used for ions with no DR rates */
1178  for( ion=0; ion<LIMELM; ++ion )
1179  {
1180  if( nsumrec[ion] > 0 )
1182  pow(10., ionbal.DR_Badnell_rate_coef_mean_ion[ion]/nsumrec[ion]);
1183  /*fprintf(ioQQQ,"DEBUG %li %.2e\n", ion , ionbal.DR_Badnell_rate_coef_mean_ion[ion] );*/
1184  }
1185 
1186  /* >>chng 06 jun 28 by Mitchell Martin, added DR rate coefficient
1187  * calculations for SII, SIII, SIV, SV, and SVI*/
1188  nelem = ipSULPHUR;
1189  /* starts with S+1 -> S0 and does the next ion with each iteration*/
1190  for( ion=0; ion<5; ion++ )
1191  {
1192 
1193  /* only do this rate if not already done by a previous approximation
1194  * so following used for ion <= 3 */
1195  if( !ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] )
1196  {
1197  /* >>chng 06 jun 28 by Mitchell Martin, added DR rate
1198  * coefficient for SII, SIII, SIV, SV, and SVI*/
1199  fitSum = 0;
1200  for( i=0; i<2; i++ )
1201  {
1202  fitSum += s_c[ion][i] * sexp( s_E[ion][i]/phycon.te );
1203  }
1204  ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = true;
1205  ionbal.DR_Badnell_rate_coef[nelem][ion] = fitSum / phycon.te32;
1206  /*>>chng 07 apr 28, use max of this or mean */
1207  /* change DR data set for S
1208  * three cases for S DR - ionbal.nDR_S_guess
1209  * 0, default larger of guess and Badnell
1210  * 1, pure Badnell
1211  * 3, scaled oxygen */
1212  if( ionbal.nDR_S_guess==0 )
1213  {
1214  /* default larger of guess or Badnell */
1215  ionbal.DR_Badnell_rate_coef[nelem][ion] =
1216  MAX2( ionbal.DR_Badnell_rate_coef[nelem][ion] ,
1218  }
1219  else if( ionbal.nDR_S_guess==1 )
1220  {
1221  /* pure Badnell - compiler will remove this non op */
1222  ionbal.DR_Badnell_rate_coef[nelem][ion] =
1223  ionbal.DR_Badnell_rate_coef[nelem][ion];
1224  }
1225  else if( ionbal.nDR_S_guess==2 )
1226  {
1227  /* scaled oxygen */
1228  ionbal.DR_Badnell_rate_coef[nelem][ion] =
1230  ionbal.DR_S_scale[ion];
1231  }
1232  else
1233  TotalInsanity();
1234  }
1235  }
1236 
1237  /* this set true with PRINT on any of the Badnell set recombination commands */
1239  {
1240  fprintf(ioQQQ,"DEBUG Badnell recombination RR, then DR, T=%.3e\n", phycon.te );
1241  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
1242  {
1243  fprintf(ioQQQ,"nelem=%li %s, RR then DR\n",
1244  nelem , elementnames.chElementNameShort[nelem] );
1245  for( ion=0; ion<nelem+1; ++ion )
1246  {
1247  if( ionbal.lgRR_Badnell_rate_coef_exist[nelem][ion] )
1248  {
1249  fprintf(ioQQQ," %.2e", ionbal.RR_Badnell_rate_coef[nelem][ion] );
1250  }
1251  else
1252  {
1253  fprintf(ioQQQ," %.2e", -1. );
1254  }
1255  }
1256  fprintf(ioQQQ,"\n" );
1257  for( ion=0; ion<nelem+1; ++ion )
1258  {
1259  if( ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] )
1260  {
1261  fprintf(ioQQQ," %.2e", ionbal.DR_Badnell_rate_coef[nelem][ion] );
1262  }
1263  else
1264  {
1265  fprintf(ioQQQ," %.2e", -1. );
1266  }
1267  }
1268  fprintf(ioQQQ,"\n\n" );
1269  }
1270  /* now print mean recombination and standard deviation */
1271  fprintf(ioQQQ,"mean DR recombination ion mean stddev stddev/mean \n" );
1272  for( ion=0; ion<LIMELM; ++ion )
1273  {
1274  double stddev;
1275  stddev = 0.;
1276  for( nelem=ion; nelem<LIMELM; ++nelem )
1277  {
1278  stddev += POW2(
1279  ionbal.DR_Badnell_rate_coef[nelem][ion]-
1281  }
1282  stddev = sqrt( stddev / MAX2( 1 , nsumrec[ion]-1 ) );
1283  fprintf(ioQQQ," %2li %.2e %.2e %.2e\n",
1284  ion ,
1286  stddev ,
1287  stddev / SDIV(ionbal.DR_Badnell_rate_coef_mean_ion[ion]) );
1288  }
1289  cdEXIT( EXIT_SUCCESS );
1290  }
1291 
1292  return;
1293 }

Generated for cloudy by doxygen 1.8.3.1