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