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