cloudy  trunk
 All Data Structures Namespaces 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-2013 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 #include "taulines.h"
23 
24 static const int MAX_FIT_PAR_DR = 9;
25 static double ***DRFitParPart1;
26 static double ***DRFitParPart2;
27 static int **nDRFitPar;
28 
29 static const int MAX_FIT_PAR_RR = 6;
30 static double ***RRFitPar;
31 
32 /* flags to recall that we have read the fits from the main data files */
33 static bool **lgDRBadnellDefined ,
37 static bool lgMustMallocRec=true;
38 static double RecNoise[LIMELM],
40 
41 static char chDRDataSource[LIMELM][LIMELM][10];
42 static char chRRDataSource[LIMELM][LIMELM][10];
43 
44 /* these enable certain debugging print statements */
45 /* #define PRINT_DR */
46 /* #define PRINT_RR */
47 
48 #if defined(PRINT_DR) || defined(PRINT_RR)
49 static const char FILE_NAME_OUT[] = "array.out";
50 #endif
51 
52 /* This function computes the standard electron density-dependent
53  * suppression factor of the collisional DR rate coefficient of H-like
54  * to Cl-like ions, based on Hugh P. Summers' 1979 report AL-R-5 report.
55  * It is then scalable for other choices of ionic charge and temperature.
56  */
58  /* atomic_number on physics scale = nuclear charge - 6 for C */
59  long int atomic_number,
60  /* ionic_charge = charge before recombination, physics scale, 3 for C+3 -> C+2 */
61  long int ionic_charge,
62  /* eden = electron density */
63  double eden,
64  /*T = temperature (K)*/
65  double T )
66 {
67 
68  /* fitting constants to compute nominal suppression factor function */
69  const double mu = 0.000; /* pseudo Voigt Lorentzian mixture */
70  const double w = 5.64586; /* suppression decay rate */
71  const double x_a0 = 10.1821; /* log10 of the electron density fitting parameter for H-like ions */
72 
73  /* a fitting constant to compute the suppression factor corrected for an
74  * estimate of surviving DR based on the lowest dipole allowed core
75  * excitation energy
76  */
77  const double c = 10.0; /* smaller c means larger fraction will survive, and vice versa */
78 
79  double s, snew, x_a, E_c;
80  double T_0, q_0, A_N; /* seed temperature, charge, and sequence selector*/
81  long int iso_sequence, N_1, N_2;
82 
83  eden = log10(eden);
84 
85  /* the isoelectronic sequence number, iso_sequence=3 for Li-like, etc */
86  iso_sequence = atomic_number - ionic_charge;
87  ASSERT( iso_sequence >= 0 );
88  if( iso_sequence==0 )
89  {
90  snew = 1.;
91  return snew;
92  }
93 
94  /* Temporarily save ionic_charge/10 needed later for excitation energy fits*/
95  realnum ionchar = ionic_charge / 10. ;
96 
97  N_1 = -1;
98  N_2 = -1;
99  /* initiate sequence-wise charge-dependent seed charge */
100  if( (iso_sequence >= 1) && (iso_sequence <= 2) ) /* 1st row sequences */
101  {
102  N_1 = 1;
103  N_2 = 2;
104  }
105  else if( (iso_sequence >= 3) && (iso_sequence <= 10) ) /* 2nd row sequences */
106  {
107  N_1 = 3;
108  N_2 = 10;
109  }
110  else if( (iso_sequence >= 11) && (iso_sequence <= 18) ) /* 3rd row sequences */
111  {
112  N_1 = 11;
113  N_2 = 18;
114  }
115  else if( (iso_sequence >= 19) && (iso_sequence <= 36) ) /* 4th row sequences */
116  {
117  N_1 = 19;
118  N_2 = 36;
119  }
120  else if( (iso_sequence >= 37) && (iso_sequence <= 54) ) /* 5th row sequences */
121  {
122  N_1 = 37;
123  N_2 = 54;
124  }
125  else if( (iso_sequence >= 55) && (iso_sequence <= 86) ) /* 6th row sequences */
126  {
127  N_1 = 55;
128  N_2 = 86;
129  }
130  else if( iso_sequence >= 87 ) /* 7th row sequences */
131  {
132  N_1 = 87;
133  N_2 = 118;
134  }
135  //fprintf(ioQQQ, "DEBUGGG %li %li %.2e %.2e %li %li %li\n",
136  // atomic_number, ionic_charge, eden, T ,
137  // N_1 , N_2 , iso_sequence);
138  ASSERT( N_1>0 && N_2>0 );
139 
140  /* initiate zig-zag approximation for A(N), which must be enveloped by two asymptotes:
141  * Amax = 12 + 10*N and Amin = 12 + 2*N
142  */
143  A_N = 12.0 + 10.0 * N_1 + (10.0 * N_1 - 2.0 * N_2) * (iso_sequence - N_1) / (N_1 - N_2);
144  ASSERT( A_N >= 16.0 );
145 
146  /* Now loop through specific sequences to assign computational estimates
147  * of lowest dipole allowed core excitation energies (these are fits to
148  * NIST statistical weighted energies) and readjust A(N) values for N<=5 sequences.
149  */
150 
151  ASSERT( iso_sequence>0 );
152  if( iso_sequence == 1 ) /* H-like ions */
153  {
154  E_c = 0.0;
155  A_N = 16.0;
156  }
157  else if( iso_sequence == 2 ) /* He-like ions */
158  {
159  E_c = 0.0;
160  A_N = 18.0;
161  }
162  else if( iso_sequence == 3 ) /* Li-like ions */
163  {
164  E_c = 1.96274 + ionchar*(20.30014 + ionchar*(-0.97103 + ionchar*( 0.85453 + ionchar*( 0.13547 + 0.02401*ionchar))));
165  A_N = 66.0;
166  }
167  else if( iso_sequence == 4 ) /* Be-like ions */
168  {
169  E_c = 5.78908 + ionchar*(34.08270 + ionchar*( 1.51729 + ionchar*(-1.21227 + ionchar*( 0.77559 - 0.00410*ionchar))));
170  A_N = 66.0;
171  }
172  else if( iso_sequence == 5 ) /* B-like ions */
173  {
174  E_c = 0.0;
175  A_N = 52.0;
176  }
177  else if( iso_sequence == 7 ) /* N-like ions */
178  {
179  E_c = 11.37092 + ionchar*(36.22053 + ionchar*( 7.08448 + ionchar*(-5.16840 + ionchar*( 2.45056 - 0.16961*ionchar))));
180  }
181  else if( iso_sequence == 11 ) /* Na-like ions */
182  {
183  E_c = 2.24809 + ionchar*(22.27768 + ionchar*(-1.12285 + ionchar*( 0.90267 + ionchar*(-0.03860 + 0.01468*ionchar))));
184  }
185  else if( iso_sequence == 12 ) /* Mg-like ions */
186  {
187  E_c = 2.74508 + ionchar*(19.18623 + ionchar*(-0.54317 + ionchar*( 0.78685 + ionchar*(-0.04249 + 0.01357*ionchar))));
188  }
189  else if( iso_sequence == 15 ) /* P-like ions */
190  {
191  E_c = 1.42762 + ionchar*( 3.90778 + ionchar*( 0.73119 + ionchar*(-1.91404 + ionchar*( 1.05059 - 0.08992*ionchar))));
192  }
193  else
194  {
195  E_c = 0.0; /* forces suppression factor to s for all T */
196  }
197 
198  /* check for low temperatures in sequences below Carbon-like*/
199  if( (iso_sequence <= 5) && (T <= 25000.0 * ionic_charge*ionic_charge) ) /* H-like to B-like ions */
200  A_N *= 2.0;
201 
202  /* initiate charge- and sequence-dependent seed charge qo
203  * qo = (1-sqrt(2/3q))*A(N)/sqrt(q)
204  */
205  q_0 = 1.0 / sqrt((double)ionic_charge);
206  q_0 = A_N * q_0 * (1.0 - 0.816497 * q_0);
207  ASSERT( q_0 > 0.0 );
208 
209  /* initiate charge-dependent seed temperature in K */
210  T_0 = 50000.0 * pow( q_0, 2. );
211 
212  /* scale log activation density to current charge and temperature */
213  x_a = x_a0 + log10( pow( ((double)ionic_charge/q_0), 7. ) * sqrt( T/T_0 ) );
214 
215  /* Now we're going to modify this standard suppression factor curve to
216  * allow for the survival of some fraction of the total DR rate at
217  * generally lower temperatures T, when appropriate.
218  */
219 
220  /* here we compute the standard suppression factor function, s( n_e, T, ionic_charge ) */
221  if( eden >= x_a )
222  {
223  s = ( mu/( 1. + pow((eden-x_a)/w, 2.) ) +
224  (1. - mu) * exp( -LN_TWO * pow((eden-x_a)/w, 2.) ) );
225  }
226  else
227  {
228  s = 1.;
229  }
230  /* converting the standard curve to the revised one allowing for
231  * survival at lower energies
232  */
233  snew = 1. + (s-1.)*exp(-(E_c*EVDEGK)/(c*T));
234 
235  ASSERT( snew >= 0. && snew <= 1. );
236  return snew;
237 }
238 
253  /* atomic number on C scale - He - 1 */
254  int nAtomicNumberCScale,
255  /* number of core electrons before capture of free electron */
256  int n_core_e_before_recomb )
257 {
258 
259  double RateCoefficient, sum;
260  int i;
261 
262  DEBUG_ENTRY( "Badnell_DR_rate_eval()" );
263  ASSERT( nAtomicNumberCScale>=0 && nAtomicNumberCScale<LIMELM );
264 
265  if( nAtomicNumberCScale==ipIRON && n_core_e_before_recomb>=12 &&
266  n_core_e_before_recomb<=18 )
267  {
268  /* these data are from table 1 of
269  *>>refer Fe DR Badnell, N., 2006, ApJ, 651, L73
270  * Fe 8+ to Fe 12+, but also include Fe13+ and Fe 14+,
271  * so these are 26-8=18 to 26-14=12
272  * increasing number of bound electrons, 0 is 14 elec, 1 is 15 elec
273  * Fe 3p^q, q=2-6
274  * these are not in badnell large dat file as of 2011 apr 24 */
275  double cFe_q[7][8] =
276  {
277  {5.636e-4, 7.390e-3, 3.635e-2, 1.693e-1, 3.315e-2, 2.288e-1, 7.316e-2, 0.},
278  {1.090e-3, 7.801e-3, 1.132e-2, 4.740e-2, 1.990e-1, 3.379e-2, 1.140e-1, 1.250e-1},
279  {3.266e-3, 7.637e-3, 1.005e-2, 2.527e-2, 6.389e-2, 1.564e-1, 0., 0.},
280  {1.074e-3, 6.080e-3, 1.887e-2, 2.540e-2, 7.580e-2, 2.773e-1, 0., 0.},
281  {9.073e-4, 3.777e-3, 1.027e-2, 3.321e-2, 8.529e-2, 2.778e-1, 0., 0.},
282  {5.335e-4, 1.827e-3, 4.851e-3, 2.710e-2, 8.226e-2, 3.147e-1, 0., 0.},
283  {7.421e-4, 2.526e-3, 4.605e-3, 1.489e-2, 5.891e-2, 2.318e-1, 0., 0.}
284  };
285 
286  /* Table 2 of Badnell 06 */
287  double EFe_q[7][8] =
288  {
289  {3.628e3, 2.432e4, 1.226e5, 4.351e5, 1.411e6, 6.589e6, 1.030e7, 0},
290  {1.246e3, 1.063e4, 4.719e4, 1.952e5, 5.637e5, 2.248e6, 7.202e6, 3.999e9},
291  {1.242e3, 1.001e4, 4.466e4, 1.497e5, 3.919e5, 6.853e5, 0. , 0.},
292  {1.387e3, 1.048e4, 3.955e4, 1.461e5, 4.010e5, 7.208e5, 0. , 0.},
293  {1.525e3, 1.071e4, 4.033e4, 1.564e5, 4.196e5, 7.580e5, 0. , 0.},
294  {2.032e3, 1.018e4, 4.638e4, 1.698e5, 4.499e5, 7.880e5, 0. , 0.},
295  {3.468e3, 1.353e4, 3.690e4, 1.957e5, 4.630e5, 8.202e5, 0. , 0.}
296  };
297  /* nion is for the above block of numbers */
298  long int nion = n_core_e_before_recomb - 12;
299  ASSERT( nion>=0 && nion <=6 );
300 
301  sum = 0;
302  i = 0;
303  /* loop over all non-zero terms */
304  for(i=0; i<8; ++i )
305  {
306  sum += (cFe_q[nion][i] * sexp( EFe_q[nion][i]/phycon.te));
307  }
308 
309  /*RateCoefficient = pow(phycon.te, -1.5) * sum;*/
310  RateCoefficient = sum / phycon.te32;
311  strcpy(chDRDataSource[nAtomicNumberCScale][nAtomicNumberCScale-n_core_e_before_recomb] ,
312  "Bad06D");
313 
314  return RateCoefficient;
315  }
316 
317  /*Invalid entries returns '-2':more electrons than protons */
318  else if( nAtomicNumberCScale < n_core_e_before_recomb )
319  {
320  RateCoefficient = -2;
321  }
322  /*Invalid entries returns '-2' if nAtomicNumberCScale and n_core_e_before_recomb are out of the range*/
323  else if( nAtomicNumberCScale >= LIMELM )
324  {
325  RateCoefficient = -2;
326  }
327  /*undefined z and n returns '-1'*/
328  else if( !lgDRBadnellDefined[nAtomicNumberCScale][n_core_e_before_recomb] )
329  {
330  RateCoefficient = -1;
331  }
332  else if( lgDRBadnellDefined[nAtomicNumberCScale][n_core_e_before_recomb] )
333  {
334  /* this branch, recombination coefficient has been defined */
335  sum = 0;
336  i = 0;
337  /* loop over all non-zero terms */
338  for(i=0; i<nDRFitPar[nAtomicNumberCScale][n_core_e_before_recomb]; ++i )
339  {
340  sum += (DRFitParPart1[nAtomicNumberCScale][n_core_e_before_recomb][i] *
341  sexp( DRFitParPart2[nAtomicNumberCScale][n_core_e_before_recomb][i]/phycon.te));
342  }
343 
344  strcpy(chDRDataSource[nAtomicNumberCScale][nAtomicNumberCScale-n_core_e_before_recomb] ,
345  "BadWeb");
346 
347  /*RateCoefficient = pow(phycon.te, -1.5) * sum;*/
348  RateCoefficient = sum / phycon.te32;
349  }
350  /*unknown invalid entries returns '-99'*/
351  else
352  {
353  RateCoefficient = -99;
354  }
355 
356  ASSERT( RateCoefficient < 1e-6 );
357 
358  return RateCoefficient;
359 }
360 
366  /* atomic number on C scale - He - 1 */
367  int nAtomicNumberCScale,
368  /* number of core electrons before capture of free electron */
369  int n_core_e_before_recomb )
370 {
371  double RateCoefficient;
372  double B, D, F;
373 
374  DEBUG_ENTRY( "Badnell_RR_rate_eval()" );
375 
376  ASSERT( nAtomicNumberCScale>=0 && nAtomicNumberCScale<LIMELM );
377 
378  if( nAtomicNumberCScale==ipIRON &&
379  n_core_e_before_recomb>=12 && n_core_e_before_recomb<=18 )
380  {
381  /* RR rate coefficients from Table 3 of
382  *>>refer Fe RR Badnell, N. 2006, ApJ, 651, L73
383  * Fe 8+ to Fe 12+, but also include Fe13+ and Fe 14+,
384  * so these are 26-8=18 to 26-14=12
385  * increasing number of bound electrons, 0 is 14 elec, 1 is 15 elec
386  * Fe 3p^q, q=2-6
387  * this is DR fit coefficients given in table 1 of Badnell 06 */
388  double parFeq[7][6] ={
389  {1.179e-9 , 0.7096, 4.508e2, 3.393e7, 0.0154, 3.977e6},
390  {1.050e-9 , 0.6939, 4.568e2, 3.987e7, 0.0066, 5.451e5},
391  {9.832e-10, 0.7146, 3.597e2, 3.808e7, 0.0045, 3.952e5},
392  {8.303e-10, 0.7156, 3.531e2, 3.554e7, 0.0132, 2.951e5},
393  {1.052e-9 , 0.7370, 1.639e2, 2.924e7, 0.0224, 4.291e5},
394  {1.338e-9 , 0.7495, 7.242e1, 2.453e7, 0.0404, 4.199e5},
395  {1.263e-9 , 0.7532, 5.209e1, 2.169e7, 0.0421, 2.917e5}
396  };
397 
398  double temp;
399  /* nion is for the above block of numbers */
400  long int nion = n_core_e_before_recomb - 12;
401  ASSERT( nion>=0 && nion <=6 );
402 
403  temp = -parFeq[nion][5]/phycon.te; /* temp = (-T2/T) */
404  B = parFeq[nion][1] + parFeq[nion][4]*exp(temp);
405  D = sqrt(phycon.te/parFeq[nion][2]); /* D = (T/T0)^1/2 */
406  F = sqrt(phycon.te/parFeq[nion][3]); /* F = (T/T1)^1/2 */
407  RateCoefficient = parFeq[nion][0]/(D*pow((1.+D),(1.-B))*pow((1.+F),(1.+B)));
408  strcpy(chRRDataSource[nAtomicNumberCScale][nAtomicNumberCScale-n_core_e_before_recomb] ,"Bad06");
409 
410  return RateCoefficient;
411  }
412 
413  /*Invalid entries returns '-2':if the z_values are smaller than equal to the n_values */
414  else if( nAtomicNumberCScale < n_core_e_before_recomb )
415  {
416  RateCoefficient = -2;
417  }
418  /*Invalid entries returns '-2' if nAtomicNumberCScale and n_core_e_before_recomb are out of the range*/
419  else if( nAtomicNumberCScale >= LIMELM )
420  {
421  RateCoefficient = -2;
422  }
423  /*undefined z and n returns '-1'*/
424  else if( !lgRRBadnellDefined[nAtomicNumberCScale][n_core_e_before_recomb] )
425  {
426  RateCoefficient = -1;
427  }
428  /* coefficients:A=RRFitPar[0], B=RRFitPar[1], T0=RRFitPar[2], T1=RRFitPar[3], DRFitParPart1=RRFitPar[4], T2=RRFitPar[5] */
429  else if( lgRRBadnellDefined[nAtomicNumberCScale][n_core_e_before_recomb] )
430  {
431 
432  /* RateCoefficient=A*[(T/T0)^1/2*(1+(T/T0)^1/2)^1-B*(1+(T/T1)^1/2)^1+B]^-1
433  where B = B + DRFitParPart1*exp(-T2/T) */
434  double temp;
435 
436  temp = -RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][5]/phycon.te; /* temp = (-T2/T) */
437  B = RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][1] +
438  RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][4]*exp(temp);
439  D = sqrt(phycon.te/RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][2]); /* D = (T/T0)^1/2 */
440  F = sqrt(phycon.te/RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][3]); /* F = (T/T1)^1/2 */
441  RateCoefficient = RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][0]/(D*pow((1.+D),(1.-B))*pow((1.+F),(1.+B)));
442  strcpy(chRRDataSource[nAtomicNumberCScale][nAtomicNumberCScale-n_core_e_before_recomb] ,"Bad06");
443  }
444 
445  /*unknown invalid entries returns '-99'*/
446  else
447  RateCoefficient = -99;
448 
449  return RateCoefficient;
450 }
451 
452 
453 /*Badnell_rec_init This code is written by Terry Yun, 2005 *
454  * It reads rate coefficient fits into 3D arrays and output array.out for testing *
455  * The testing can be commented out */
456 void Badnell_rec_init( void )
457 {
458 
459  double par_C[MAX_FIT_PAR_DR];
460  double par_E[MAX_FIT_PAR_DR];
461  char chLine[INPUT_LINE_LENGTH];
462  int NuclearCharge=-1, NumberElectrons=-1;
463  int count, number;
464  double temp_par[MAX_FIT_PAR_RR];
465  int M_state, W_state;
466 
467  const int NBLOCK = 2;
468  int data_begin_line[NBLOCK];/*it tells you where the data set begins(begins with 'Z')*/
469  int length_of_line; /*this variable checks for a blank line*/
470  FILE *ioDATA;
471  const char* chFilename;
472  int yr, mo, dy;
473  char *chs;
474 
475  const int BIGGEST_INDEX_TO_USE = 103;
476 
477  /* Declaration of data file name array - done by Kausalya */
478  long TheirIndexToOurIndex[BIGGEST_INDEX_TO_USE];
479  char string[120];
480  double value;
481  bool lgEOL;
482  long int i1;
483  long INDX=0,INDP=0,N=0,S=0,L=0,J=0,maxINDX=0,loopindex=0,max_N_of_data=-1;
484  bool lgFlag = true;
485 
486  static int nCalled = 0;
487 
488  const char* cdDATAFILE[] =
489  {
490  /* the list of filenames for Badnell DR, one to two electron */
491  "",
492  "UTA/nrb00_h_he1ic12.dat",
493  "UTA/nrb00_h_li2ic12.dat",
494  "UTA/nrb00_h_be3ic12.dat",
495  "UTA/nrb00_h_b4ic12.dat",
496  "UTA/nrb00_h_c5ic12.dat",
497  "UTA/nrb00_h_n6ic12.dat",
498  "UTA/nrb00_h_o7ic12.dat",
499  "UTA/nrb00_h_f8ic12.dat",
500  "UTA/nrb00_h_ne9ic12.dat",
501  "UTA/nrb00_h_na10ic12.dat",
502  "UTA/nrb00_h_mg11ic12.dat",
503  "UTA/nrb00_h_al12ic12.dat",
504  "UTA/nrb00_h_si13ic12.dat",
505  "UTA/nrb00_h_p14ic12.dat",
506  "UTA/nrb00_h_s15ic12.dat",
507  "UTA/nrb00_h_cl16ic12.dat",
508  "UTA/nrb00_h_ar17ic12.dat",
509  "UTA/nrb00_h_k18ic12.dat",
510  "UTA/nrb00_h_ca19ic12.dat",
511  "UTA/nrb00_h_sc20ic12.dat",
512  "UTA/nrb00_h_ti21ic12.dat",
513  "UTA/nrb00_h_v22ic12.dat",
514  "UTA/nrb00_h_cr23ic12.dat",
515  "UTA/nrb00_h_mn24ic12.dat",
516  "UTA/nrb00_h_fe25ic12.dat",
517  "UTA/nrb00_h_co26ic12.dat",
518  "UTA/nrb00_h_ni27ic12.dat",
519  "UTA/nrb00_h_cu28ic12.dat",
520  "UTA/nrb00_h_zn29ic12.dat"
521  };
522  //End of modification
523 
524  DEBUG_ENTRY( "Badnell_rec_init()" );
525 
526  /* must only do this once */
527  if( nCalled > 0 )
528  {
529  return;
530  }
531  ++nCalled;
532 
533 # if defined(PRINT_DR) || defined(PRINT_RR)
534  FILE *ofp = open_data( FILE_NAME_OUT, "w", AS_LOCAL_ONLY );
535 # endif
536 
537  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
538  {
539  for( long nelem=ipISO; nelem < LIMELM; nelem++ )
540  {
541  if( nelem < 2 || dense.lgElmtOn[nelem] )
542  {
543  for( long ipHi=0; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
544  {
545  for( long k=0; k<NUM_DR_TEMPS; ++k )
546  iso_sp[ipISO][nelem].fb[ipHi].DielecRecombVsTemp[k] = 0.;
547  }
548  }
549  }
550  }
551 
552  /* Modification done by Kausalya
553  * Start - Try to open all the 29 data files.*/
554  for( long nelem=ipHELIUM; nelem<LIMELM; nelem++)
555  {
556  if( nelem < 2 || dense.lgElmtOn[nelem] )
557  {
558  ioDATA= open_data( cdDATAFILE[nelem], "r" );
559 
560  lgFlag = true;
561  ASSERT(ioDATA);
562 
563  for( long i=0; i<BIGGEST_INDEX_TO_USE; i++ )
564  TheirIndexToOurIndex[i] = -1;
565 
566  /* Reading lines */
567  while(lgFlag)
568  {
569  if(read_whole_line(string,sizeof(string),ioDATA)!=NULL)
570  {
571  if( nMatch("INDX INDP ",string) )
572  {
573  /* ignore next line of data */
574  if( read_whole_line( string , (int)sizeof(string) , ioDATA ) == NULL )
575  {
576  fprintf( ioQQQ, " Badnell data file appears to be corrupted.\n");
578  }
579 
580  /* This one should be real data */
581  while( read_whole_line(string, (int)sizeof(string), ioDATA) != NULL )
582  {
583  if( strcmp(string,"\n")==0 )
584  {
585  lgFlag = false;
586  break;
587  }
588 
589  i1=3;
590  INDX=(long)FFmtRead(string,&i1,sizeof(string),&lgEOL);
591  if( INDX >= BIGGEST_INDEX_TO_USE )
592  {
593  INDX--;
594  lgFlag = false;
595  break;
596  }
597 
598  ASSERT( INDX < BIGGEST_INDEX_TO_USE );
599 
600  INDP=(long)FFmtRead(string,&i1,sizeof(string),&lgEOL);
601  ASSERT( INDP >= 1 );
602 
603  if(INDP==1)
604  {
605  if( (i1=nMatch("1S1 ",string)) > 0 )
606  {
607  i1 += 4;
608  N=(long)FFmtRead(string,&i1,sizeof(string),&lgEOL);
609  ASSERT( N>=1 );
610  }
611  else
612  {
613  TotalInsanity();
614  }
615 
616  if( (i1=nMatch(" (",string)) > 0 )
617  {
618  i1 += 6;
619  S=(long)FFmtRead(string,&i1,sizeof(string),&lgEOL);
620  /* S in file is 3 or 1, we need 1 or 0 */
621  ASSERT( S==1 || S==3 );
622  }
623  else
624  {
625  TotalInsanity();
626  }
627 
628  /* move i1 one further to get L */
629  i1++;
630  L=(long)FFmtRead(string,&i1,sizeof(string),&lgEOL);
631  ASSERT( L >= 0 && L < N );
632 
633  /* move i1 two further to get J */
634  i1 += 2;
635  J=(long)FFmtRead(string,&i1,sizeof(string),&lgEOL);
636  ASSERT( J <= ( L + (int)((S+1)/2) ) &&
637  J >= ( L - (int)((S+1)/2) ) && J >= 0 );
638 
639  /* if line in data file is higher N than highest considered, stop reading. */
640  if( N<= iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max + iso_sp[ipHE_LIKE][nelem].nCollapsed_max )
641  TheirIndexToOurIndex[INDX] = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[N][L][S];
642  else
643  {
644  /* Current line is not being used,
645  * decrement INDX so maxINDX is set correctly below. */
646  INDX--;
647  lgFlag = false;
648  break;
649  }
650 
651  /* Must adjust index if in 2^3Pj term */
652  if( N==2 && L==1 && S==3 )
653  {
654  if( J==0 )
655  TheirIndexToOurIndex[INDX] = 3;
656  else if( J==1 )
657  TheirIndexToOurIndex[INDX] = 4;
658  else
659  {
660  ASSERT( J==2 );
661  ASSERT( TheirIndexToOurIndex[INDX] == 5 );
662  }
663  }
664  max_N_of_data = MAX2( max_N_of_data, N );
665  }
666  else
667  {
668  // Stop parsing the tuple since INDP!=1
669  lgFlag = false;
670  }
671  }
672  }
673  }
674  else
675  {
676  // End of file is reached.
677  lgFlag = false;
678  }
679  }
680 
681  maxINDX =INDX;
682  ASSERT( maxINDX > 0 );
683  ASSERT( maxINDX < BIGGEST_INDEX_TO_USE );
684  /* reset INDX */
685  INDX = 0;
686  lgFlag = true;
687  while(lgFlag)
688  {
689  if(read_whole_line(string,sizeof(string),ioDATA)!=NULL)
690  {
691  /* to access the first table whose columns are INDX ,INDP */
692  if( nMatch("INDX TE= ",string) )
693  {
694  lgFlag = false;
695  /* we found the beginning of the data array */
696  /* ignore next line of data */
697  if( read_whole_line( string , (int)sizeof(string) , ioDATA ) == NULL )
698  {
699  fprintf( ioQQQ, " Badnell data file appears to be corrupted.\n");
701  }
702 
703  /* This one should be real data */
704  while( read_whole_line(string, (int)sizeof(string), ioDATA) != NULL )
705  {
706  /* If we find this string, we have reached the end of the table. */
707  if( nMatch("PRTF",string) || INDX >= maxINDX || INDX<0 )
708  break;
709 
710  i1=3;
711  INDX=(long)FFmtRead(string,&i1,sizeof(string),&lgEOL);
712  if( INDX>maxINDX )
713  break;
714 
715  freeBound *fb;
716 
717  if( TheirIndexToOurIndex[INDX] < iso_sp[ipHE_LIKE][nelem].numLevels_max &&
718  TheirIndexToOurIndex[INDX] > 0 )
719  fb = &iso_sp[ipHE_LIKE][nelem].fb[TheirIndexToOurIndex[INDX]];
720  else
721  continue;
722 
723  for(loopindex=0;loopindex<10;loopindex++)
724  {
725  value=(double)FFmtRead(string,&i1,sizeof(string),&lgEOL);
726  fb->DielecRecombVsTemp[loopindex] += value;
727  }
728 
729  /* data are broken into two lines, read second line here */
730  if( read_whole_line( string , (int)sizeof(string) , ioDATA ) == NULL )
731  {
732  fprintf( ioQQQ, " Badnell data file appears to be corrupted.\n");
734  }
735 
736  /* start of data for second line */
737  i1 = 13;
738  for(loopindex=10;loopindex<19;loopindex++)
739  {
740  value=(double)FFmtRead(string,&i1,sizeof(string),&lgEOL);
741  fb->DielecRecombVsTemp[loopindex] += value;
742  }
743  }
744  }
745  }
746  else
747  lgFlag = false;
748  }
749  fclose(ioDATA);
750  ASSERT( maxINDX > 0 );
751  ASSERT( maxINDX < BIGGEST_INDEX_TO_USE );
752  ASSERT( max_N_of_data > 0 );
753 
754  if( max_N_of_data < iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max + iso_sp[ipHE_LIKE][nelem].nCollapsed_max )
755  {
756  long indexOfMaxN;
757  L = -1;
758  S = -1;
759 
760  /* This loop extrapolates nLS data to nLS states */
761  for( long i=TheirIndexToOurIndex[maxINDX]+1;
763  {
764  long ipISO = ipHE_LIKE;
765  L = L_(i);
766  S = S_(i);
767 
768  if( L > 4 )
769  continue;
770 
771  indexOfMaxN = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[max_N_of_data][L][S];
772  for(loopindex=0;loopindex<19;loopindex++)
773  {
774  iso_sp[ipHE_LIKE][nelem].fb[i].DielecRecombVsTemp[loopindex] =
775  iso_sp[ipHE_LIKE][nelem].fb[indexOfMaxN].DielecRecombVsTemp[loopindex] *
776  pow3( (double)max_N_of_data/(double)iso_sp[ipHE_LIKE][nelem].st[i].n());
777  }
778  }
779 
780  /* Get the N of the highest resolved singlet P (in the model, not the data) */
781  indexOfMaxN =
783 
784  /* This loop extrapolates nLS data to collapsed n levels, just use highest singlet P data */
785  for( long i=iso_sp[ipHE_LIKE][nelem].numLevels_max-iso_sp[ipHE_LIKE][nelem].nCollapsed_max;
786  i<iso_sp[ipHE_LIKE][nelem].numLevels_max; i++ )
787  {
788  for(loopindex=0;loopindex<19;loopindex++)
789  {
790  iso_sp[ipHE_LIKE][nelem].fb[i].DielecRecombVsTemp[loopindex] =
791  iso_sp[ipHE_LIKE][nelem].fb[indexOfMaxN].DielecRecombVsTemp[loopindex] *
792  pow3( (double)iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max/
793  (double)iso_sp[ipHE_LIKE][nelem].st[i].n());
794  }
795  }
796  }
797  }
798  }
799 
800  for( long i=0; i<NBLOCK; ++i )
801  {
802  /* set to really large negative number - crash if used before being redefined */
803  data_begin_line[i] = INT_MIN;
804  }
805 
806  chFilename = "badnell_dr.dat";
807  ioDATA = open_data( chFilename, "r" );
808 
809  count = 0;
810  number = 0;
811 
812  /*Find out the number line where the data starts
813  * there are two main blocks of data and each starts with a Z in column 2 */
814  while( read_whole_line(chLine, (int)sizeof(chLine), ioDATA) != NULL )
815  {
816  count++;
817 
818  if( chLine[2]=='Z' )
819  {
820  /* number has to be 0 or 1, and indicates the first or second block of data
821  * count is the line number for the start of that block */
822  data_begin_line[number] = count;
823  ASSERT( number < NBLOCK );
824  number++;
825  }
826  }
827 
828  /*set a flag for a undefined data*/
829  if( lgMustMallocRec )
830  {
831  nDRFitPar = (int**)MALLOC( LIMELM*sizeof( int*) );
832  lgDRBadnellDefined = (bool **)MALLOC( LIMELM*sizeof(bool*) );
833  lgDR_BadWeb_exist = (bool **)MALLOC( LIMELM*sizeof(bool*) );
834  lgDRBadnellDefinedPart2 = (bool **)MALLOC( LIMELM*sizeof(bool*) );
835  lgRRBadnellDefined = (bool **)MALLOC( LIMELM*sizeof(bool*) );
836 
837  DRFitParPart1 = (double ***)MALLOC( LIMELM*sizeof(double**) );
838  DRFitParPart2 = (double ***)MALLOC( LIMELM*sizeof(double**) );
839  RRFitPar = (double ***)MALLOC( LIMELM*sizeof(double**) );
840  }
841 
842  for( long nelem=0; nelem<LIMELM; nelem++ )
843  {
844  if( lgMustMallocRec )
845  {
846  nDRFitPar[nelem] = (int*)MALLOC( (nelem+1)*sizeof( int) );
847  lgDR_BadWeb_exist[nelem] = (bool *)MALLOC( (nelem+1)*sizeof(bool) );
848  lgDRBadnellDefined[nelem] = (bool *)MALLOC( (nelem+1)*sizeof(bool) );
849  lgDRBadnellDefinedPart2[nelem] = (bool *)MALLOC( (nelem+1)*sizeof(bool) );
850  lgRRBadnellDefined[nelem] = (bool *)MALLOC( (nelem+1)*sizeof(bool) );
851 
852  DRFitParPart1[nelem] = (double **)MALLOC( (nelem+1)*sizeof(double*) );
853  DRFitParPart2[nelem] = (double **)MALLOC( (nelem+1)*sizeof(double*) );
854  RRFitPar[nelem] = (double **)MALLOC( (nelem+1)*sizeof(double*) );
855  }
856  for( long ion=0; ion<nelem+1; ++ion )
857  {
858  if( lgMustMallocRec )
859  {
860  DRFitParPart1[nelem][ion] = (double *)MALLOC( MAX_FIT_PAR_DR*sizeof(double) );
861  DRFitParPart2[nelem][ion] = (double *)MALLOC( MAX_FIT_PAR_DR*sizeof(double) );
862  RRFitPar[nelem][ion] = (double *)MALLOC( MAX_FIT_PAR_RR*sizeof(double) );
863  }
864  lgDRBadnellDefined[nelem][ion] = false;
865  lgDRBadnellDefinedPart2[nelem][ion] = false;
866  lgRRBadnellDefined[nelem][ion] = false;
867 
868  /*set fitting coefficients to zero initially*/
869  for( long k=0; k<MAX_FIT_PAR_DR; k++ )
870  {
871  DRFitParPart1[nelem][ion][k] = 0;
872  DRFitParPart2[nelem][ion][k] = 0;
873  }
874  for( long k=0; k<MAX_FIT_PAR_RR; k++ )
875  {
876  RRFitPar[nelem][ion][k] = 0;
877  }
878  }
879  }
880  lgMustMallocRec = false;
881 
882  count = 0;
883  /*Start from beginning to read in again*/
884  fseek(ioDATA, 0, SEEK_SET);
885  /* read magic number for DR data */
886  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
887  {
888  fprintf( ioQQQ, " DISASTER PROBLEM Badnell_rec_init could not read first line of badnell_dr.dat.\n");
890  }
891  count++;
892 
893  /* look for ')' on the line, magic number comes after it */
894  if( (chs = strchr_s(chLine, ')'))==NULL )
895  {
896  /* format is incorrect */
897  fprintf( ioQQQ, " DISASTER PROBLEM Badnell_rec_init data file incorrect format.\n");
899  }
900 
901  ++chs;
902  sscanf(chs, "%4i%2i%2i",&yr, &mo, &dy);
903  /* check magic number - the date on the line */
904  int dr_yr = 2012, dr_mo = 6, dr_dy = 28;
905  if((yr != dr_yr) || (mo != dr_mo) || (dy != dr_dy))
906  {
907  fprintf(ioQQQ,
908  "DISASTER PROBLEM Badnell_rec_init The version of %s I found (%i %i %i) is not the current version (%i %i %i).\n",
909  chFilename, yr, mo, dy, dr_yr, dr_mo, dr_dy);
910  fprintf(ioQQQ," The first line of the file is the following\n %s\n", chLine );
912  }
913 
914  while( read_whole_line(chLine, (int)sizeof(chLine), ioDATA) != NULL )
915  {
916  count++;
917  length_of_line = (int)strlen(chLine);
918 
919  /*reading in coefficients DRFitParPart1 */
920  if( count > data_begin_line[0] && count < data_begin_line[1] && length_of_line >3 )
921  {
922  /*set array par_C to zero */
923  for( long i=0; i<MAX_FIT_PAR_DR; i++ )
924  {
925  par_C[i] = 0;
926  }
927  sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf%lf%lf%lf",
928  &NuclearCharge, &NumberElectrons, &M_state, &W_state, &par_C[0], &par_C[1], &par_C[2],
929  &par_C[3], &par_C[4], &par_C[5], &par_C[6], &par_C[7], &par_C[8]);
930  /* data files have atomic number on physics scale, convert to C scale
931  * for following code */
932  long int NuclearChargeM1 = NuclearCharge-1;
933 
934  if(M_state == 1 && NuclearChargeM1 < LIMELM )
935  {
936  /*Set a flag to '1' when the indices are defined */
937  ASSERT( NumberElectrons < LIMELM );
938  ASSERT( NuclearChargeM1 < LIMELM );
939  lgDRBadnellDefined[NuclearChargeM1][NumberElectrons] = true;
940 
941  /*counting the number of coefficients */
942  nDRFitPar[NuclearChargeM1][NumberElectrons] = 9;
943  for( long i=8; i>=0; i-- )
944  {
945  if( par_C[i] == 0 )
946  --nDRFitPar[NuclearChargeM1][NumberElectrons];
947  else
948  break;
949  }
950 
951  /*assign the values into array */
952  for( long i=0; i<9; i++ )
953  DRFitParPart1[NuclearChargeM1][NumberElectrons][i] = par_C[i];
954  }
955  }
956  }
957 
958  /*starting again to read in E values */
959  fseek(ioDATA, 0, SEEK_SET);
960  count = 0;
961  while( read_whole_line(chLine, (int)sizeof(chLine), ioDATA) != NULL )
962  {
963  count++;
964  length_of_line = (int)strlen(chLine);
965  if( count > data_begin_line[1] && length_of_line > 3 )
966  {
967 
968  /*set array par_E to zero*/
969  for( long i=0; i<MAX_FIT_PAR_DR; i++ )
970  {
971  par_E[i] = 0;
972  }
973  sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf%lf%lf%lf",
974  &NuclearCharge, &NumberElectrons, &M_state, &W_state, &par_E[0], &par_E[1], &par_E[2],
975  &par_E[3], &par_E[4], &par_E[5], &par_E[6], &par_E[7], &par_E[8]);
976  /* data file is on physics scale but we use C scale */
977  long int NuclearChargeM1 = NuclearCharge-1;
978 
979  if(M_state == 1 && NuclearChargeM1<LIMELM)
980  {
981  ASSERT( NumberElectrons < LIMELM );
982  ASSERT( NuclearChargeM1 < LIMELM );
983  lgDRBadnellDefinedPart2[NuclearChargeM1][NumberElectrons] = true;
984 
985  /*counting the number of coefficients */
986  nDRFitPar[NuclearChargeM1][NumberElectrons] = 9;
987  for( long i=8; i>=0; i-- )
988  {
989  if( par_E[i] == 0 )
990  --nDRFitPar[NuclearChargeM1][NumberElectrons];
991  else
992  break;
993  }
994 
995  /*assign the values into array*/
996  for( long i=0; i<nDRFitPar[NuclearChargeM1][NumberElectrons]; i++ )
997  DRFitParPart2[NuclearChargeM1][NumberElectrons][i] = par_E[i];
998  }
999  }
1000  }
1001 
1002  fclose( ioDATA );
1003 
1004  /*output coefficients for defined values for testing */
1005 # ifdef PRINT_DR
1006  for( long nelem=0; nelem<LIMELM; nelem++ )
1007  {
1008  for( int ion=0; ion<nelem+1;++ion )
1009  {
1010  if( lgDRBadnellDefined[nelem][ion] )
1011  {
1012  fprintf(ofp, "%i %i %e %e %e %e %e %e %e %e %e\n",
1013  nelem, ion, DRFitParPart1[nelem][ion][0],
1014  DRFitParPart1[nelem][ion][1], DRFitParPart1[nelem][ion][2],
1015  DRFitParPart1[nelem][ion][3], DRFitParPart1[nelem][ion][4],
1016  DRFitParPart1[nelem][ion][5], DRFitParPart1[nelem][ion][6],
1017  DRFitParPart1[nelem][ion][7], DRFitParPart1[nelem][ion][8]);
1018  }
1019  }
1020  }
1021  for( long nelem=0; nelem<LIMELM; nelem++ )
1022  {
1023  for( int ion=0; ion<nelem+1; ion++ )
1024  {
1025  if( lgDRBadnellDefinedPart2[nelem][ion] )
1026  {
1027  fprintf(ofp, "%i %i %e %e %e %e %e %e %e %e %e\n",
1028  nelem, ion, DRFitParPart2[nelem][ion][0],
1029  DRFitParPart2[nelem][ion][1], DRFitParPart2[nelem][ion][2],
1030  DRFitParPart2[nelem][ion][3], DRFitParPart2[nelem][ion][4],
1031  DRFitParPart2[nelem][ion][5], DRFitParPart2[nelem][ion][6],
1032  DRFitParPart2[nelem][ion][7], DRFitParPart2[nelem][ion][8]);
1033  }
1034  }
1035  }
1036  fclose(ofp);
1037 # endif
1038 
1039  /*checking for the match of lgDRBadnellDefined and lgDRBadnellDefinedPart2 -
1040  * Both have to be defined*/
1041  bool lgDRBadnellBothDefined = true;
1042  for( int nelem=0; nelem<LIMELM; nelem++ )
1043  {
1044  for( int ion=0; ion<nelem+1; ion++ )
1045  {
1046  /* check that first and second half of DR fitting coefficients
1047  * are both defined */
1048  if( lgDRBadnellDefined[nelem][ion] != lgDRBadnellDefinedPart2[nelem][ion] )
1049  {
1050  fprintf( ioQQQ, "DR %i, RR %i: %c %c\n", nelem, ion,
1051  TorF(lgDRBadnellDefined[nelem][ion]),
1052  TorF(lgDRBadnellDefinedPart2[nelem][ion]));
1053  fprintf( ioQQQ, "PROBLEM ion_recomb_Badnell first and second half of Badnell DR not consistent.\n");
1054  lgDRBadnellBothDefined = false;
1055  }
1056  }
1057  }
1058 
1059  if( !lgDRBadnellBothDefined )
1060  {
1061  /* disaster - DR files are not consistent */
1062  fprintf(ioQQQ,
1063  "DISASTER PROBLEM The DR data files are corrupted - part 1 and 2 do not agree.\n");
1064  fprintf(ioQQQ," Start again with a fresh copy of the data directory\n" );
1066  }
1067 
1068  /* now do radiative recombination */
1069  chFilename = "badnell_rr.dat";
1070  ioDATA = open_data( chFilename, "r" );
1071 
1072  /* read magic number for RR data */
1073  {
1074  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1075  {
1076  fprintf( ioQQQ, " DISASTER PROBLEM Badnell_rec_init could not read first line of badnell_rr.dat.\n");
1078  }
1079  /* this is just before date, which we use for magic number */
1080  if( (chs = strchr_s(chLine, ')'))==NULL )
1081  {
1082  /* format is incorrect */
1083  fprintf( ioQQQ, " DISASTER PROBLEM Badnell_rec_init data file incorrect format.\n");
1085  }
1086  ++chs;
1087  sscanf(chs, "%4i%2i%2i", &yr, &mo, &dy);
1088  int rr_yr = 2011, rr_mo = 4, rr_dy = 12;
1089  if((yr != rr_yr)||(mo != rr_mo)||(dy != rr_dy))
1090  {
1091  fprintf(ioQQQ,"DISASTER PROBLEM The version of %s I found (%i %i %i) is not the current version (%i %i %i).\n",
1092  chFilename, yr, mo, dy, rr_yr, rr_mo, rr_dy);
1093  fprintf(ioQQQ," The line was as follows:\n %s\n", chLine );
1095  }
1096  }
1097 
1098  while( read_whole_line(chLine, (int)sizeof(chLine), ioDATA) != NULL )
1099  {
1100  /*read in coefficients - first set array par to zero */
1101  for( long i=0; i<MAX_FIT_PAR_RR; i++ )
1102  {
1103  temp_par[i] = 0;
1104  }
1105  if(chLine[0] != '#')
1106  {
1107  sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf",
1108  &NuclearCharge, &NumberElectrons, &M_state, &W_state, &temp_par[0], &temp_par[1],
1109  &temp_par[2], &temp_par[3], &temp_par[4], &temp_par[5]);
1110  long NuclearChargeM1 = NuclearCharge-1;
1111 
1112  if(M_state == 1 && NuclearChargeM1<LIMELM)
1113  {
1114  ASSERT( NuclearChargeM1 < LIMELM );
1115  ASSERT( NumberElectrons <= LIMELM );
1116  /*Set a flag to '1' when the indices are defined */
1117  lgRRBadnellDefined[NuclearChargeM1][NumberElectrons] = true;
1118  /*assign the values into array */
1119  for( long i=0; i<MAX_FIT_PAR_RR; i++ )
1120  RRFitPar[NuclearChargeM1][NumberElectrons][i] = temp_par[i];
1121  }
1122  }
1123  }
1124 
1125  /*output coefficients for defined values for testing */
1126 # ifdef PRINT_RR
1127  count = 0;
1128  for( long nelem=0; nelem<LIMELM; nelem++ )
1129  {
1130  for( long ion=0; ion<nelem+1; ion++ )
1131  {
1132  if( lgRRBadnellDefined[nelem][ion] )
1133  {
1134  fprintf(ofp, "%li %li %e %e %e %e %e %e\n",
1135  nelem, ion, RRFitPar[nelem][ion][0],
1136  RRFitPar[nelem][ion][1], RRFitPar[nelem][ion][2],
1137  RRFitPar[nelem][ion][3],
1138  RRFitPar[nelem][ion][4], RRFitPar[nelem][ion][5]);
1139  count++;
1140  }
1141  }
1142  }
1143  fprintf(ofp, "total lines are %i ", count);
1144 
1145  fclose(ofp);
1146 # endif
1147 
1148  fclose(ioDATA);
1149 
1150  {
1151  enum {DEBUG_LOC=false};
1152  if( DEBUG_LOC )
1153  {
1154  for( int nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
1155  {
1156  fprintf(ioQQQ,"\nDEBUG rr rec\t%i",nelem);
1157  for( int ion=0; ion<=nelem; ++ion )
1158  {
1159  fprintf(ioQQQ,"\t%.2e", Badnell_RR_rate_eval(nelem+1 , nelem-ion ) );
1160  }
1161  fprintf(ioQQQ,"\n");
1162  fprintf(ioQQQ,"DEBUG dr rec\t%i",nelem);
1163  for( int ion=0; ion<=nelem; ++ion )
1164  {
1165  fprintf(ioQQQ,"\t%.2e", Badnell_DR_rate_eval(nelem+1 , nelem-ion ) );
1166  }
1167  fprintf(ioQQQ,"\n");
1168  }
1170  }
1171  }
1172 
1173  // gaussian noise for dielectronic recombination coefficients guesses
1174  // set with SET DIELECTRONIC RECOMBINATION NOISE command
1175  if( ionbal.guess_noise !=0. )
1176  {
1177  for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
1178  /* log normal noise with dispersion entered on command line */
1179  /* NB the seed for rand was set when the command was parsed */
1180  RecNoise[nelem] = pow(10., RandGauss( 0. , ionbal.guess_noise ) );
1181  }
1182  else
1183  {
1184  for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
1185  RecNoise[nelem] = 1.;
1186  }
1187 
1188  // initialize some products
1189  for( long nelem=0; nelem<LIMELM; ++nelem )
1190  DR_Badnell_rate_coef_mean_ion[nelem] = 0.;
1191 
1192  return;
1193 }
1194 
1195 /*ion_recom_calculate calculate radiative and dielectronic recombination rate coefficients */
1197 {
1198  static double TeUsed = -1 , EdenUsed = -1.;
1199 
1200  DEBUG_ENTRY( "ion_recom_calculate()" );
1201 
1202  /* do not reevaluate if change in temperature is small */
1203  if( fp_equal(phycon.te,TeUsed) && fp_equal( dense.eden , EdenUsed ))
1204  return;
1205 
1206  // collisional suppression factors
1207  //CollisSuppres();
1208 
1209  TeUsed = phycon.te;
1210  EdenUsed = dense.eden;
1211 
1212  for( long nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
1213  {
1214 
1215  for( long ion=0; ion < nelem+1; ++ion )
1216  {
1217  long int n_bnd_elec_before_recom ,
1218  n_bnd_elec_after_recom;
1219 
1220  n_bnd_elec_before_recom = nelem-ion;
1221  n_bnd_elec_after_recom = nelem-ion+1;
1222 
1223  // will insure these are >=0 at end
1224  ionbal.DR_Badnell_rate_coef[nelem][ion] = -1.;
1225  ionbal.RR_rate_coef_used[nelem][ion] = 0.;
1226  strcpy( chDRDataSource[nelem][ion] , "none" );
1227  strcpy( chRRDataSource[nelem][ion] , "none" );
1228 
1229  /* Badnell dielectronic recombination rate coefficients */
1230  if( (ionbal.DR_Badnell_rate_coef[nelem][ion] =
1232  /* atomic number on C scale */
1233  nelem,
1234  /* number of core electrons before capture of free electron,
1235  * for bare ion this is zero */
1236  n_bnd_elec_before_recom )) >= 0. )
1237  {
1238  lgDR_BadWeb_exist[nelem][ion] = true;
1239  }
1240  else
1241  {
1242  /* real rate does not exist, will use mean later */
1243  lgDR_BadWeb_exist[nelem][ion] = false;
1244  }
1245 
1246  /* save D. Verner's radiative recombination rate coefficient
1247  * needed for rec cooling, cm3 s-1 */
1248  ionbal.RR_Verner_rate_coef[nelem][ion] =
1250  /* number number of physics scale */
1251  nelem+1 ,
1252  /* number of protons on physics scale */
1253  n_bnd_elec_after_recom ,
1254  phycon.te );
1255 
1256  /* Badnell radiative recombination rate coefficients */
1258  /* atomic number on C scale */
1259  nelem,
1260  /* number of core electrons before capture of free electron */
1261  n_bnd_elec_before_recom )) >= 0. )
1262  {
1263  ionbal.RR_rate_coef_used[nelem][ion] = ionbal.RR_Badnell_rate_coef[nelem][ion];
1264  }
1265  else
1266  {
1267  strcpy( chRRDataSource[nelem][ion] , "Verner" );
1268  ionbal.RR_rate_coef_used[nelem][ion] = ionbal.RR_Verner_rate_coef[nelem][ion];
1269  }
1270  }
1271  // recombination to bare nuclei has no DR
1272  ionbal.DR_Badnell_rate_coef[nelem][nelem] = 0.;
1273  strcpy(chDRDataSource[nelem][nelem] , "NA");
1274  }
1275 
1276  /* this branch starts idiosyncratic single ions */
1277  double Fe_Gu_c[9][6] = {
1278  { 2.50507e-11, 5.60226e-11, 1.85001e-10, 3.57495e-9, 1.66321e-7, 0. },/*fit params for Fe+6*/
1279  { 9.19610e-11, 2.92460e-10, 1.02120e-9, 1.14852e-8, 3.25418e-7, 0. }, /* fitting params for Fe+7 */
1280  { 9.02625e-11, 6.22962e-10, 5.77545e-9, 1.78847e-8, 3.40610e-7, 0. }, /* fitting params for Fe+8 */
1281  { 9.04286e-12, 9.68148e-10, 4.83636e-9, 2.48159e-8, 3.96815e-7, 0. }, /* fitting params for Fe+9 */
1282  { 6.77873e-10, 1.47252e-9, 5.31202e-9, 2.54793e-8, 3.47407e-7, 0. }, /* fitting params for Fe+10 */
1283  { 1.29742e-9, 4.10172e-9, 1.23605e-8, 2.33615e-8, 2.97261e-7, 0. }, /* fitting params for Fe+11 */
1284  { 8.78027e-10, 2.31680e-9, 3.49333e-9, 1.16927e-8, 8.18537e-8, 1.54740e-7 },/*fit params for Fe+12*/
1285  { 2.23178e-10, 1.87313e-9, 2.86171e-9, 1.38575e-8, 1.17803e-7, 1.06251e-7 },/*fit params for Fe+13*/
1286  { 2.17263e-10, 7.35929e-10, 2.81276e-9, 1.32411e-8, 1.15761e-7, 4.80389e-8 }/*fit params for Fe+14*/
1287  },
1288 
1289  Fe_Gu_E[9][6] = {
1290  { 8.30501e-2, 8.52897e-1, 3.40225e0, 2.23053e1, 6.80367e1, 0. }, /* fitting params for Fe+6 */
1291  { 1.44392e-1, 9.23999e-1, 5.45498e0, 2.04301e1, 7.06112e1, 0. }, /* fitting params for Fe+7 */
1292  { 5.79132e-2, 1.27852e0, 3.22439e0, 1.79602e1, 6.96277e1, 0. }, /* fitting params for Fe+8 */
1293  { 1.02421e-1, 1.79393e0, 4.83226e0, 1.91117e1, 6.80858e1, 0. }, /* fitting params for Fe+9 */
1294  { 1.24630e-1, 6.86045e-1, 3.09611e0, 1.44023e1, 6.42820e1, 0. }, /* fitting params for Fe+10 */
1295  { 1.34459e-1, 6.63028e-1, 2.61753e0, 1.30392e1, 6.10222e1, 0. }, /* fitting params for Fe+11 */
1296  { 7.79748e-2, 5.35522e-1, 1.88407e0, 8.38459e0, 3.38613e1, 7.89706e1 }, /*fitting params for Fe+12*/
1297  { 8.83019e-2, 6.12756e-1, 2.36035e0, 9.61736e0, 3.64467e1, 8.72406e1 }, /*fitting params for Fe+13*/
1298  { 1.51322e-1, 5.63155e-1, 2.57013e0, 9.08166e0, 3.69528e1, 1.08067e2 } /* fitting params for Fe+14*/
1299  };
1300 
1301  /* do a series of special cases for Fe DR */
1302  double te_eV32 = sqrt(pow3(phycon.te_eV));
1303 
1304  /* >>chng 06 jul 07 by Mitchell Martin, added DR rate coefficient
1305  * calculations for Fe+6->Fe+5 through Fe+14->Fe+13
1306  * this is still for nelem = ipIRON from the previous calculation
1307  * starts with Fe+6 -> Fe+5 and does the next ion with each iteration */
1308  for( long ion=0; ion<9; ion++ )
1309  {
1310  /* only do this rate if not already done by a previous approximation */
1311  if( ionbal.DR_Badnell_rate_coef[ipIRON][ion+5]<0. )
1312  {
1313  double fitSum = 0; /* resets the fitting parameter calculation */
1314  for( long i=0; i<6; i++ )
1315  {
1316  fitSum += Fe_Gu_c[ion][i] * sexp( Fe_Gu_E[ion][i]/phycon.te_eV );
1317  }
1318  strcpy(chDRDataSource[ipIRON][ion+5] , "GuPC");
1319  lgDR_BadWeb_exist[ipIRON][ion+5] = true;
1320  ionbal.DR_Badnell_rate_coef[ipIRON][ion+5] = fitSum / te_eV32;
1321  }
1322  }
1323  /* this is end of Fe DR rates */
1324 
1325  // use C08 mean for stability
1326  double BadnelDR_RateSave[LIMELM] =
1327  {
1328  3.78e-13, 1.70e-12, 8.14e-12, 1.60e-11, 2.38e-11,
1329  6.42e-11, 5.97e-11, 1.47e-10, 1.11e-10, 3.26e-10,
1330  1.88e-10, 2.06e-10, 4.14e-10, 3.97e-10, 2.07e-10,
1331  2.46e-10, 3.38e-10, 3.15e-10, 9.70e-11, 6.49e-11,
1332  6.93e-10, 3.70e-10, 3.29e-11, 4.96e-11, 5.03e-11,
1333  2.91e-12, 4.62e-14, 0.00e+00, 0.00e+00, 0.00e+00
1334  };
1335  for( long nelem=0; nelem < LIMELM; ++nelem )
1336  {
1338  BadnelDR_RateSave[nelem] * RecNoise[nelem] *
1339  // default of unity, set with SET DIELECTRONIC RECOMBINATION KLUDGE SCALE command
1340  ionbal.DR_mean_scale[nelem];
1341  }
1342 
1343  // iron is special case with Arnaud & Raymond 1992
1344  // use mean which is low T dr and AR which is high temp
1345  for( long ion=0; ion < ipIRON+1; ++ion )
1346  {
1347  if( ionbal.DR_Badnell_rate_coef[ipIRON][ion] < 0. )
1348  {
1351  + atmdat_dielrec_fe(ion+1, phycon.te );
1352  strcpy(chDRDataSource[ipIRON][ion] , "mean+");
1353  }
1354  }
1355  // this routine will return something for all ions - even if just a guess
1356  for( long nelem=0; nelem < LIMELM; ++nelem )
1357  {
1358  for( long ion=0; ion < nelem+1; ++ion )
1359  if( ionbal.DR_Badnell_rate_coef[nelem][ion] < 0. )
1360  {
1361  strcpy(chDRDataSource[nelem][ion] , "mean");
1363  }
1364  }
1365 
1366  // collisional suppression of DR
1367  for( long nelem=ipLITHIUM; nelem < LIMELM; ++nelem )
1368  {
1369  for( long ion=0; ion < nelem-1; ++ion )
1370  {
1371  // ASSERT(DielSupprsFactor[ion]>=0 && DielSupprsFactor[ion]<=1. );
1372  // old very simple expression
1373  //ionbal.DR_Badnell_rate_coef[nelem][ion] *= DielSupprsFactor[ion];
1374 
1375  // DR collisional suppression based on Badnell rates
1376  ionbal.DR_Badnell_rate_coef[nelem][ion] *= CollisSuppres(
1377  /* This routine takes the following arguments:
1378  * atomic_number = nuclear charge */
1379  nelem+1,
1380  /*ionic_charge = ionic charge*/
1381  ion+1,
1382  /*eden = electron density */
1383  dense.eden,
1384  /*T = temperature (K)*/
1385  phycon.te );
1386 
1387  ASSERT(ionbal.DR_Badnell_rate_coef[nelem][ion] >= 0);
1388  ASSERT(ionbal.RR_rate_coef_used[nelem][ion] >= 0);
1389  }
1390  }
1391 
1392  /* this set true with PRINT RECOMBINATION recombination commands */
1394  {
1395 
1396  fprintf(ioQQQ,"\n\n RR recombination data sources \n" );
1397 
1398  for( long loop=0;loop<30;loop+=10)
1399  {
1400  fprintf(ioQQQ,"\n\n ");
1401  for(long ion=loop; ion<loop+10; ++ion )
1402  {
1403  fprintf(ioQQQ,"&%7li",ion);
1404  }
1405  fprintf(ioQQQ,"\\\\\n" );
1406  for( long nelem=loop; nelem<LIMELM; ++nelem )
1407  {
1408  fprintf(ioQQQ,"%2li %5s ",nelem+1 , elementnames.chElementNameShort[nelem] );
1409  long limit = MIN2(nelem+1,loop+10);
1410  for( long ion=loop; ion<limit; ++ion )
1411  {
1412  fprintf(ioQQQ,"&%7s",chRRDataSource[nelem][ion] );
1413  }
1414  for( long ion=limit; ion<loop+10; ++ion )
1415  {
1416  fprintf(ioQQQ,"&%7s",chRRDataSource[nelem][ion] );
1417  }
1418  fprintf(ioQQQ,"\\\\\n" );
1419  }
1420  }
1421  fprintf(ioQQQ,"\nData sources\n");
1422  fprintf(ioQQQ,"Bad06: Badnell, N., 2006, ApJ, 167, 334B\n");
1423  fprintf(ioQQQ,"Verner: Verner & Ferland, 1996, ApJS, 103, 467\n");
1424 
1425  fprintf(ioQQQ,"\n\n DR recombination data sources \n" );
1426 
1427  for( long loop=0;loop<30;loop+=10)
1428  {
1429  fprintf(ioQQQ,"\n\n ");
1430  for(long ion=loop; ion<loop+10; ++ion )
1431  {
1432  fprintf(ioQQQ,"&%7li",ion);
1433  }
1434  fprintf(ioQQQ,"\\\\\n" );
1435  for( long nelem=loop; nelem<LIMELM; ++nelem )
1436  {
1437  fprintf(ioQQQ,"%2li %5s ",
1438  nelem+1 , elementnames.chElementNameShort[nelem] );
1439  long limit = MIN2(nelem+1,loop+10);
1440  for( long ion=loop; ion<limit; ++ion )
1441  {
1442  fprintf(ioQQQ,"&%7s",chDRDataSource[nelem][ion] );
1443  }
1444  for( long ion=limit; ion<loop+10; ++ion )
1445  {
1446  fprintf(ioQQQ,"&%7s",chDRDataSource[nelem][ion] );
1447  }
1448  fprintf(ioQQQ,"\\\\\n" );
1449  }
1450  }
1451  fprintf(ioQQQ,"\nData sources\nBadWeb: Badnell web site http://amdpp.phys.strath.ac.uk/tamoc/DR/\n");
1452  fprintf(ioQQQ,"Bad06D: Badnell, N., 2006, ApJ, 651, L73\n");
1453  fprintf(ioQQQ,"GuPC: Gu, M. private communication\n");
1454 
1455  fprintf(ioQQQ,"\n\nDEBUG Badnell recombination RR, then DR, T=%.3e\n", phycon.te );
1456  for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
1457  {
1458  fprintf(ioQQQ,"nelem=%li %s, RR then DR\n",
1459  nelem , elementnames.chElementNameShort[nelem] );
1460  for( long ion=0; ion<nelem+1; ++ion )
1461  {
1462  fprintf(ioQQQ," %.2e", ionbal.RR_rate_coef_used[nelem][ion] );
1463  }
1464  fprintf(ioQQQ,"\n" );
1465  for( long ion=0; ion<nelem+1; ++ion )
1466  {
1467  fprintf(ioQQQ," %.2e", ionbal.DR_Badnell_rate_coef[nelem][ion] );
1468  }
1469  fprintf(ioQQQ,"\n\n" );
1470  }
1471  /* now print mean recombination and standard deviation */
1472  fprintf(ioQQQ,"mean DR recombination ion mean \n" );
1473  for( long ion=0; ion<LIMELM; ++ion )
1474  {
1475  fprintf(ioQQQ," %2li %.2e \n",
1476  ion ,
1478  }
1479 
1480  fprintf( ioQQQ, "\n\nCollisSuppres finds following dielectronic"
1481  " recom suppression factors, eden=%10.3e\n", dense.eden );
1482  fprintf( ioQQQ, "nelem ion fac \n" );
1483  for( long nelem=0; nelem<LIMELM; ++nelem )
1484  {
1485  for( long ion=0; ion < nelem+1; ion++ )
1486  {
1487  fprintf( ioQQQ, "%3ld %4ld %10.3e\n", nelem+1 , ion+1,
1488  CollisSuppres(
1489  /* This routine takes the following arguments:
1490  * atomic_number = nuclear charge */
1491  nelem+1,
1492  /*ionic_charge = ionic charge*/
1493  ion+1,
1494  /*eden = electron density */
1495  dense.eden,
1496  /*T = temperature (K) */
1497  phycon.te )
1498  );
1499 
1500  }
1501  fprintf( ioQQQ, "\n");
1502  }
1503 
1504  cdEXIT( EXIT_SUCCESS );
1505  }
1506  return;
1507 }
#define MIN2
Definition: cddefines.h:765
double ** DR_Badnell_rate_coef
Definition: ionbal.h:204
double ** RR_Badnell_rate_coef
Definition: ionbal.h:204
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:616
STATIC double Badnell_DR_rate_eval(int nAtomicNumberCScale, int n_core_e_before_recomb)
double DielecRecombVsTemp[NUM_DR_TEMPS]
Definition: freebound.h:35
const int ipHE_LIKE
Definition: iso.h:63
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
double ** RR_Verner_rate_coef
Definition: ionbal.h:215
const int NISO
Definition: cddefines.h:265
static double DR_Badnell_rate_coef_mean_ion[LIMELM]
char TorF(bool l)
Definition: cddefines.h:714
static const int N
#define MAX2
Definition: cddefines.h:786
long int nCollapsed_max
Definition: iso.h:487
long nMatch(const char *chKey, const char *chCard)
Definition: service.cpp:451
static bool ** lgDRBadnellDefinedPart2
static const int MAX_FIT_PAR_DR
t_phycon phycon
Definition: phycon.cpp:6
static bool ** lgDRBadnellDefined
t_dense dense
Definition: dense.cpp:24
sys_float sexp(sys_float x)
Definition: service.cpp:914
T pow3(T a)
Definition: cddefines.h:942
const double EVDEGK
Definition: physconst.h:186
FILE * ioQQQ
Definition: cddefines.cpp:7
bool lgRecom_Badnell_print
Definition: ionbal.h:209
vector< freeBound > fb
Definition: iso.h:452
double te_eV
Definition: phycon.h:14
static t_ADfA & Inst()
Definition: cddefines.h:179
t_elementnames elementnames
Definition: elementnames.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
static double RecNoise[LIMELM]
STATIC double CollisSuppres(long int atomic_number, long int ionic_charge, double eden, double T)
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:816
#define MALLOC(exp)
Definition: cddefines.h:505
void Badnell_rec_init(void)
t_ionbal ionbal
Definition: ionbal.cpp:5
const int ipIRON
Definition: cddefines.h:334
long int n_HighestResolved_max
Definition: iso.h:505
#define L_(A_)
Definition: iso.h:21
static char chRRDataSource[LIMELM][LIMELM][10]
static double *** DRFitParPart1
#define STATIC
Definition: cddefines.h:101
realnum guess_noise
Definition: ionbal.h:231
double atmdat_dielrec_fe(long int ion, double t)
static int ** nDRFitPar
static char chDRDataSource[LIMELM][LIMELM][10]
double DR_mean_scale[LIMELM]
Definition: ionbal.h:220
static bool lgMustMallocRec
float realnum
Definition: cddefines.h:107
#define EXIT_FAILURE
Definition: cddefines.h:144
const int INPUT_LINE_LENGTH
Definition: cddefines.h:258
void ion_recom_calculate(void)
bool lgElmtOn[LIMELM]
Definition: dense.h:146
#define cdEXIT(FAIL)
Definition: cddefines.h:438
#define S_(A_)
Definition: iso.h:22
static double *** DRFitParPart2
const char * strchr_s(const char *s, int c)
Definition: cddefines.h:1443
static double *** RRFitPar
const double LN_TWO
Definition: physconst.h:50
char chElementNameShort[LIMELM][CHARS_ELEMENT_NAME_SHORT]
Definition: elementnames.h:21
double RandGauss(double xMean, double s)
Definition: service.cpp:1643
#define NUM_DR_TEMPS
Definition: freebound.h:7
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:461
int strcmp(const char *s2)
Definition: parser.h:177
static bool ** lgRRBadnellDefined
#define ASSERT(exp)
Definition: cddefines.h:582
const int ipH_LIKE
Definition: iso.h:62
const int LIMELM
Definition: cddefines.h:262
double rad_rec(long int iz, long int in, double t)
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:688
const int ipHELIUM
Definition: cddefines.h:310
double eden
Definition: dense.h:190
bool lgEOL(void) const
Definition: parser.h:98
STATIC double Badnell_RR_rate_eval(int nAtomicNumberCScale, int n_core_e_before_recomb)
#define S(I_, J_)
long int numLevels_max
Definition: iso.h:493
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
double te
Definition: phycon.h:11
static bool ** lgDR_BadWeb_exist
const int ipHYDROGEN
Definition: cddefines.h:309
const int ipLITHIUM
Definition: cddefines.h:311
double te32
Definition: phycon.h:48
double ** RR_rate_coef_used
Definition: ionbal.h:212
static const int MAX_FIT_PAR_RR
#define EXIT_SUCCESS
Definition: cddefines.h:142
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:381