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