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