PLUTO
make_tables.c
Go to the documentation of this file.
1 #include "pluto.h"
2 #include "cooling_defs.h"
3 
4 /* ********************************************************************* */
5 void Solve_System (Ion *X, double Ne, double T)
6 /*!
7  * Solve sytem A*x = rhs ( adapted from the GSL library ).
8  * This is used to compute the level populations for the ions.
9  *
10  *********************************************************************** */
11 {
12  int i, j, k, nlev;
13  double scrh,tmpx, d;
14  double **q;
15  double **M, *rhs;
16  int *p;
17 
18  nlev = X->nlev;
19  M = ARRAY_2D(nlev, nlev, double);
20  rhs = ARRAY_1D(nlev, double);
21  p = ARRAY_1D(nlev, int);
22  q = ARRAY_2D (nlev, nlev, double);
23 
24  for (i = 0; i < nlev; i++) {
25  rhs[i] = 0.0;
26  p[i] = 0.0;
27  for (j = 0; j < nlev; j++) {
28  M[i][j] = 0.0;
29  q[i][j] = 0.0;
30  }
31  }
32 
33 /* -------------------------------------------
34  Define qij
35  ------------------------------------------- */
36 
37  for (i = 0; i < nlev; i++) {
38  for (j = i + 1; j < nlev; j++) {
39  if ( (!X->isMAP) && (!X->isCV) && (!X->isH) && (!X->isCHEB) )
40  scrh = lagrange(X->Tom, X->omega[i][j], T, X->nTom, i, j); /* remove i and j after testing */
41  if (X->isCV) scrh = X->omega[i][j][0] * pow( (T / pow(10.,X->omega[i][j][2])), X->omega[i][j][1]);
42  if (X->isMAP) scrh = X->omega[i][j][0] * pow( T/10000., X->omega[i][j][1]);
43  if (X->isH) {
44  if (T<55000) scrh = X->omega[i][j][0] + T*X->omega[i][j][1] + T*T*X->omega[i][j][2] + T*T*T*X->omega[i][j][3];
45  else scrh = X->omega[i][j][4] + T*X->omega[i][j][5] + T*T*X->omega[i][j][6] + T*T*T*X->omega[i][j][7];
46  }
47  if (X->isCHEB) {
48  tmpx = 0.6199646 * log(T) - 6.2803580;
49  scrh = 0.5 * X->omega[i][j][0] + X->omega[i][j][1] * tmpx + X->omega[i][j][2] * (2.*tmpx*tmpx - 1) + X->omega[i][j][3] * (4.*tmpx*tmpx*tmpx - 3*tmpx);
50  scrh = exp(scrh);
51  }
52 
53  q[j][i] = scrh/X->wght[j];
54  q[j][i] *= 8.629e-6/sqrt(T);
55  q[i][j] = (X->wght[j]/X->wght[i])*q[j][i]*exp( - X->dE[i][j]/(kB*T));
56  }}
57  for (i = 0; i < nlev; i++) q[i][i] = 0.0;
58 
59 /* -----------------------------------------
60  compute coefficient matrix
61  ----------------------------------------- */
62 
63  for (i = 0 ; i < nlev ; i++) {
64  for (j = 0 ; j < nlev ; j++) {
65  scrh = 0.0;
66  if (j == i) {
67  for (k = 0 ; k < nlev ; k++) {
68  if (k != i) scrh += Ne*q[i][k];
69  if (k < i) scrh += X->A[i][k];
70  }
71  M[i][j] = -scrh;
72  }
73  if (j < i) M[i][j] = Ne*q[j][i];
74  if (j > i) M[i][j] = Ne*q[j][i] + X->A[j][i];
75  }}
76 
77 /* -------------------------------------
78  Replace 1st eq. with normalization
79  condition and define rhs
80  ------------------------------------- */
81 
82  for (j = 0; j < nlev; j++) {
83  M[nlev-1][j] = 1.0;
84  rhs[j] = 0.0;
85  }
86  rhs[nlev-1] = 1.0;
87 
88 /* ----------------------------------
89  Solve equations by LU decomp
90  ---------------------------------- */
91  LUDecompose (M, nlev, p, &d);
92  LUBackSubst (M, nlev, p, rhs);
93  for (i = 0 ; i < nlev ; i++) X->Ni[i] = rhs[i];
94 }
95 
96 /* ********************************************************************* */
97 void Symmetrize_Coeff (Ion *X)
98 /*
99  *
100  *
101  *
102  *********************************************************************** */
103 {
104  int i,j,n;
105 
106 /* --------------------------------------------------
107  symmetrize dE, A, and Omega
108  in order to prevent swapped index in the
109  ion routines
110  -------------------------------------------------- */
111 
112  for (i = 0; i < X->nlev; i++) {
113  for (j = 0; j < X->nlev; j++) {
114  if (X->dE[i][j] > 0.0) {
115  X->dE[j][i] = X->dE[i][j];
116  X->A[j][i] = X->A[i][j];
117  for (n = 0; n < X->nTom; n++){
118  if ( (X->omega[i][j][n] > 0.0)
119  || (X->isH && X->omega[i][j][n] != 0.0)
120  || (X->isCHEB && X->omega[i][j][n] != 0.0) )
121  X->omega[j][i][n] = X->omega[i][j][n];
122  }
123  }
124  }}
125 
126 /* ----------------------------------------------------
127  2. convert energy to correct units
128  ---------------------------------------------------- */
129 
130  for (i = 0; i < X->nlev; i++){
131  for (j = 0; j < X->nlev; j++) {
132  if (X->dE[i][j] < 0.0) { X->dE[i][j] = 0.0; X->A[i][j] = 0.0; for (n = 0; n < X->nTom; n++) X->omega[i][j][n] = 0.000; }
133  if (X->A[i][j] < 0.0) { X->A[i][j] = 0.0; for (n = 0; n < X->nTom; n++) X->omega[i][j][n] = 0.000; }
134  if (X->dE[i][j] > 0.0) X->dE[i][j] = 1.2399e-4/(1.e-8*X->dE[i][j]) ; /* multply by hc to get energy */
135  if ( (X->omega[i][j][0] < 0.0)&&(!X->isCV)&&(!X->isH)&&(!X->isCHEB)&&(!X->isMAP) ) for (n = 0; n < X->nTom; n++) X->omega[i][j][n] = 0.000;
136  }}
137 
138 
139 /* --------------------------------------------------
140  3. set A(i,j) = 0 when j >= i
141  -------------------------------------------------- */
142 
143  for (i = 0; i < X->nlev; i++) {
144  for (j = i; j < X->nlev; j++) {
145  X->A[i][j] = 0.0;
146  }}
147 
148 /* ---------------------------------------------------
149  4. set diagonal elements of dE and omega to zero
150  --------------------------------------------------- */
151 
152  for (i = 0; i < X->nlev; i++) {
153  for (n = 0; n < X->nTom; n++) X->omega[i][i][n] = 0.000;
154  X->dE[i][i] = 0.0;
155  X->A[i][i] = 0.0;
156  }
157 
158 }
159 
160 /* ********************************************************************* */
161 double lagrange (double *x, double *y, double xp, int n, int ii, int jj)
162 /*
163  * PURPOSE:
164  *
165  * Return lagrange interpolation for a tabulated function.
166  *
167  * x : a vector of n-components with the tabulated values
168  * of the independent variables x;
169  *
170  * y : a vecotr of n-components with the tabulated values
171  * y[n] = y(x[n]);
172  *
173  * xp: the point at which interpolation is desired
174  *
175  * n : the number of points in the table, also
176  * the degree of the intrpolant
177  *
178  *
179  ********************************************************** */
180 {
181  int i, k, j;
182  double scrh, yp;
183 
184  yp = 0.0;
185  if (xp < x[0]) return (y[0]);
186  if (xp > x[n-1]) return (y[n-1]);
187 
188  for (i = 0; i < n; i++){
189  scrh = 1.0;
190  for (k = 0; k < n; k++){
191  if (k == i) continue;
192  scrh *= (xp - x[k])/(x[i] - x[k]);
193  }
194  yp += scrh*y[i];
195  }
196  return(yp);
197 }
198 
199 /* ********************************************************************* */
200 int Create_Ion_Coeff_Tables(double ***tbl)
201 /*!
202  * Compute and save to memory the coefficients depending only
203  * on temperature used at runtime to compute the ionization balance.
204  *
205  *********************************************************************** */
206 #define ZERO_5X 0.0, 0.0, 0.0, 0.0, 0.0
207 {
208  double T, t4, tmprec, lam, f1, f2, x1, x2, P, Q, F[2], Tn[9];
209 
210  long int nv, i, j, ti1, ti2, cindex;
211 
212  static double **coll_ion, **rad_rec, **diel_rec, **chtr_hp, **chtr_h, **chtr_he, **tot_recs;
213 
214 /* Collisional ionization rates , fit parameters - Voronov 1997, ATOMIC DATA AND NUCLEAR DATA TABLES 65, 1-35 */
215 
216 double coll_ion_P[] = { 0., 0., 1.
217  C_EXPAND(0., 1., 1., 1., 1.)
218  N_EXPAND(0., 0., 1., 1., 0.)
219  O_EXPAND(0., 1., 1., 0., 1.)
220  Ne_EXPAND(1., 0., 1., 1., 1.)
221  S_EXPAND(1., 1., 1., 1., 1.)
222  Fe_EXPAND(0., 1., 0.) };
223 double coll_ion_A[] = { 0.291e-7, 0.175e-7, 0.0000
224  C_EXPAND(0.685e-7, 0.185e-7, 0.635e-8, 0.150e-8, 0.0000)
225  N_EXPAND(0.482e-7, 0.298e-7, 0.810e-8, 0.371e-8, 0.0000)
226  O_EXPAND(0.359e-7, 0.139e-7, 0.931e-8, 0.102e-7, 0.0000)
227  Ne_EXPAND(0.150e-7, 0.198e-7, 0.703e-8, 0.424e-8, 0.0000)
228  S_EXPAND(0.549e-7, 0.681e-7, 0.214e-7, 0.166e-7, 0.0000)
229  Fe_EXPAND(0.252e-6, 0.221e-7, 0.0000) }; /* in cm^3/s */
230 double coll_ion_X[] = { 0.232, 0.180, 0.265
231  C_EXPAND( 0.193, 0.286, 0.427, 0.416, 0.0000)
232  N_EXPAND(0.0652, 0.310, 0.350, 0.549, 0.0167)
233  O_EXPAND( 0.073, 0.212, 0.270, 0.614, 0.630)
234  Ne_EXPAND(0.0329, 0.295, 0.0677, 0.0482, 0.0000)
235  S_EXPAND( 0.100, 0.693, 0.353, 1.030, 0.0000)
236  Fe_EXPAND(0.7010, 0.033, 0.0000) };
237 double coll_ion_K[] = { 0.39, 0.35, 0.25
238  C_EXPAND(0.25, 0.24, 0.21, 0.13, 0.02)
239  N_EXPAND(0.42, 0.30, 0.24, 0.18, 0.74)
240  O_EXPAND(0.34, 0.22, 0.27, 0.27, 0.17)
241  Ne_EXPAND(0.43, 0.20, 0.39, 0.58, 0.00)
242  S_EXPAND(0.25, 0.21, 0.24, 0.14, 0.00)
243  Fe_EXPAND(0.25, 0.45, 0.17) };
244 double coll_ion_Tmin[] = { 0., 0., 0.
245  C_EXPAND(0., 0., 0., 0., 0.)
246  N_EXPAND(0., 0., 0., 0., 0.)
247  O_EXPAND(0., 0., 0., 0., 0.)
248  Ne_EXPAND(0., 0., 0., 0., 0.)
249  S_EXPAND(0., 0., 0., 0., 0.)
250  Fe_EXPAND(0., 0., 0.) }; /* in eV */
251 
252 /*
253 double coll_ion_Tmin[28] = { 1., 1., 3., 1., 1., 2., 3., 20.,
254  1., 2., 2., 4., 4., 1., 2., 3., 4., 5.,
255  1., 3., 3., 5., 7., 1., 1., 2., 2., 3. }; */
256 
257 
258 /* Charge transfer with H II ionization rates, fit parameters - Kingdon & Ferland 1996, ApJSS */
259 
260 double chtrH_ion_a[] = { 0.00, 0.00, 0.00
261  C_EXPAND( 0.00, 0.00, 0.00, 0.00, 0.00)
262  N_EXPAND(4.55e-3, 0.00, 0.00, 0.00, 0.00)
263  O_EXPAND( 7.4e-2, 0.00, 0.00, 0.00, 0.00)
264  Ne_EXPAND( 0.00, 0.00, 0.00, 0.00, 0.00)
265  S_EXPAND( 0.00, 0.00, 0.00, 0.00, 0.00)
266  Fe_EXPAND( 5.40, 2.10, 0.00) }; /* in 10^(-9) cm^3/s */
267 double chtrH_ion_b[] = { 0.00, 0.00, 0.00
268  C_EXPAND( 0.00, 0.00, 0.00, 0.00, 0.00)
269  N_EXPAND(-0.29, 0.00, 0.00, 0.00, 0.00)
270  O_EXPAND( 0.47, 0.00, 0.00, 0.00, 0.00)
271  Ne_EXPAND( 0.00, 0.00, 0.00, 0.00, 0.00)
272  S_EXPAND( 0.00, 0.00, 0.00, 0.00, 0.00)
273  Fe_EXPAND( 0.00, 7.72e-2, 0.00) };
274 double chtrH_ion_c[] = { 0.00, 0.00, 0.00
275  C_EXPAND( 0.00, 0.00, 0.00, 0.00, 0.00)
276  N_EXPAND(-0.92, 0.00, 0.00, 0.00, 0.00)
277  O_EXPAND(24.37, 0.00, 0.00, 0.00, 0.00)
278  Ne_EXPAND( 0.00, 0.00, 0.00, 0.00, 0.00)
279  S_EXPAND( 0.00, 0.00, 0.00, 0.00, 0.00)
280  Fe_EXPAND( 0.00, -0.41, 0.00)};
281 double chtrH_ion_d[] = { 0.00, 0.00, 0.00
282  C_EXPAND( 0.00, 0.00, 0.00, 0.00, 0.00)
283  N_EXPAND(-8.38, 0.00, 0.00, 0.00, 0.00)
284  O_EXPAND(-0.74, 0.00, 0.00, 0.00, 0.00)
285  Ne_EXPAND( 0.00, 0.00, 0.00, 0.00, 0.00)
286  S_EXPAND( 0.00, 0.00, 0.00, 0.00, 0.00)
287  Fe_EXPAND( 0.00, -7.31, 0.00) };
288 double chtrH_ion_upT[] = { 0.00, 0.00, 0.00
289  C_EXPAND(1.0, 0.00, 0.00, 0.00, 0.00)
290  N_EXPAND(5.0, 0.00, 0.00, 0.00, 0.00)
291  O_EXPAND(1.0, 0.00, 0.00, 0.00, 0.00)
292  Ne_EXPAND(0.00, 0.00, 0.00, 0.00, 0.00)
293  S_EXPAND(1.0, 0.00, 0.00, 0.00, 0.00)
294  Fe_EXPAND(1.0, 10.0, 0.0) }; /* in 10^4 K */
295 
296 /* Charge transfer with H recombination rates, fit parameters - Kingdon & Ferland 1996, ApJSS */
297 
298 double chtrH_rec_a[] = { 0.00, 7.46e-6, 0.00
299  C_EXPAND(1.76e-9, 1.67e-4, 3.25, 332.46, 0.00)
300  N_EXPAND(1.01e-3, 3.05e-1, 4.54, 3.28, 0.00)
301  O_EXPAND( 1.04, 1.04, 3.98, 0.252, 0.00)
302  Ne_EXPAND( 0.00, 1.00e-5, 14.73, 6.47, 0.00)
303  S_EXPAND(3.82e-7, 1.00e-5, 2.29, 6.44, 0.00)
304  Fe_EXPAND( 0.00, 1.26, 0.00) };
305  /* in 10^(-9) cm^3/s , per Ion - 1 (e.g., the coeffs for He II recomb to He I are written to the He I position */
306 double chtrH_rec_b[] = { 0.00, 2.06, 0.00
307  C_EXPAND( 8.33, 2.79, 0.21, -0.11, 0.00)
308  N_EXPAND( -0.29, 0.60, 0.57, 0.52, 0.00)
309  O_EXPAND(3.15e-2, 0.27, 0.26, 0.63, 0.00)
310  Ne_EXPAND( 0.00, 0.00, 4.52e-2, 0.54, 0.00)
311  S_EXPAND( 11.10, 0.00, 4.02e-2, 0.13, 0.00)
312  Fe_EXPAND( 0.00, 7.72e-2, 0.00) };
313 double chtrH_rec_c[] = { 0.00, 9.93, 0.00
314  C_EXPAND(4278.78, 304.72, 0.19, -0.995, 0.00)
315  N_EXPAND( -0.92, 2.65, -0.65, -0.52, 0.00)
316  O_EXPAND( -0.61, 2.02, 0.56, 2.08, 0.00)
317  Ne_EXPAND( 0.00, 0.00, -0.84, 3.59, 0.00)
318  S_EXPAND(2.57e+4, 0.00, 1.59, 2.69, 0.00)
319  Fe_EXPAND( 0.00, -0.41, 0.00) };
320 double chtrH_rec_d[] = { 0.00, -3.89, 0.00
321  C_EXPAND(-6.41, -4.07, -3.29, -1.58e-3, 0.00)
322  N_EXPAND(-8.38, -0.93, -0.89, -0.19, 0.00)
323  O_EXPAND(-9.73, -5.92, -2.62, -4.16, 0.00)
324  Ne_EXPAND( 0.00, 0.00, -0.31, -5.22, 0.00)
325  S_EXPAND(-8.22, 0.00, -6.06, -5.69, 0.00)
326  Fe_EXPAND( 0.00, -7.31, 0.00) };
327 double chtrH_rec_upT[] = { 0.00, 1.e+5, 1.e+7
328  C_EXPAND(1.e+4, 5.e+4, 1.e+5, 1.e+5, 0.00)
329  N_EXPAND(5.e+4, 1.e+5, 1.e+5, 3.e+5, 0.00)
330  O_EXPAND(1.e+4, 1.e+5, 5.e+4, 3.e+4, 0.00)
331  Ne_EXPAND( 0.00, 5.e+4, 5.e+4, 3.e+4, 0.00)
332  S_EXPAND(1.e+4, 3.e+4, 3.e+4, 3.e+4, 0.00)
333  Fe_EXPAND( 0.00, 1.e+5, 0.00) }; /* in K */
334 
335 /* Charge transfer of O ions with H recombination rates, fit parameters -
336  M. Rakovic', J. G. Wang, D. R. Schultz, and P. C. Stancil (2001)
337  ORNL Charge Transfer Database */
338 
339 double chtrH_rec_o1[8] = {3.348706E+00, 1.281922E+00, -1.870663E+00,-7.195090E-01,-2.473806E-01, 8.577410E-02, 1.721933E-02, 9.213735E-03};
340 double chtrH_rec_o2[8] = {3.266558E+00, 2.264383E+00, -1.019642E+00,-1.090796E+00,-4.403829E-01,-1.369024E-01, 1.606288E-01, 1.925501E-01};
341 double chtrH_rec_o3[8] = {4.860561E+00, 1.680463E+00, -1.110986E+00,-1.346895E+00,-2.697744E-01, 1.329470E-01,-4.511763E-02, 4.614949E-02};
342 double chtrH_rec_o4[8] = {4.680553E+00, 2.456278E+00, -1.026270E+00,-1.587952E+00,-1.318074E-01, 1.451866E-01,-1.526627E-01, 7.596889E-02};
343 double chtrH_rec_o5[8] = {5.723788E+00, 1.771316E+00, -1.254652E+00,-9.009546E-01,-3.558194E-01 -7.138940E-02, 2.846941E-02, 1.313952E-01};
344 double chtrH_rec_Tmin = 1.0e+2, chtrH_rec_Tmax = 1.0e+10;
345 
346 /* Charge transfer with He - ORNL Charge Transfer Database
347  http://www-cfadc.phy.ornl.gov/astro/ps/data/cx/helium/rates/fits.data */
348 
349 /* 1st stage - up to 5000K */
350 double chtrHe_rec_a1[] = { 0.00, 0.00, 0.00
351  C_EXPAND(0.00, 0.00, 1.12, 3.12e-7, 0.00)
352  N_EXPAND(0.00, 4.84e-1, 2.05, 1.26e-2, 0.00)
353  O_EXPAND(0.00, 7.10e-3, 1.12, 0.997, 0.00)
354  Ne_EXPAND(0.00, 1.00e-5, 1.00e-5, 1.77, 0.00)
355  S_EXPAND(0.00, 0.00, 3.58, 7.44e-4, 0.00)
356  Fe_EXPAND(0.00, 0.00, 0.00) };/* in 10^(-9) cm^3/s */
357 double chtrHe_rec_b1[] = { 0.00, 0.00, 0.00
358  C_EXPAND(0.00, 0.00, 0.42, -7.37e-2, 0.00)
359  N_EXPAND(0.00, 0.92, 0.23, 1.55, 0.00)
360  O_EXPAND(0.00, 2.60, 0.42, 0.40, 0.00)
361  Ne_EXPAND(0.00, 0.00, 0.00, 0.14, 0.00)
362  S_EXPAND(0.00, 0.00, 7.77e-3, 0.34, 0.00)
363  Fe_EXPAND(0.00, 0.00, 0.00) };
364 double chtrHe_rec_c1[] = { 0.00, 0.00, 0.00
365  C_EXPAND(0.00, 0.00, -0.69, 3.50e+1, 0.00)
366  N_EXPAND(0.00, 2.37, -0.72, 11.2, 0.00)
367  O_EXPAND(0.00, 8.99, -0.71, -0.46, 0.00)
368  Ne_EXPAND(0.00, 0.00, 0.00, 4.88e-2, 0.00)
369  S_EXPAND(0.00, 0.00, -0.94, 3.74, 0.00)
370  Fe_EXPAND(0.00, 0.00, 0.00) };
371 double chtrHe_rec_d1[] = { 0.00, 0.00, 0.00
372  C_EXPAND(0.00, 0.00, -0.34, 2.40, 0.00)
373  N_EXPAND(0.00, -10.2, -0.19, -7.82, 0.00)
374  O_EXPAND(0.00, -0.78, -1.98e-2, -0.35, 0.00)
375  Ne_EXPAND(0.00, 0.00, 0.00, -3.35, 0.00)
376  S_EXPAND(0.00, 0.00, -0.30, -5.18, 0.00)
377  Fe_EXPAND(0.00, 0.00, 0.00) };
378 /* 2nd stage - 5000K to 10000K */
379 double chtrHe_rec_a2[] = { 0.00, 0.00, 0.00
380  C_EXPAND(0.00, 0.00, 1.12, 3.12e-7, 0.00)
381  N_EXPAND(0.00, 4.84e-1, 2.05, 1.26e-2, 0.00)
382  O_EXPAND(0.00, 7.10e-3, 1.12, 0.997, 0.00)
383  Ne_EXPAND(0.00, 8.48e-3, 1.00e-5, 1.77, 0.00)
384  S_EXPAND(0.00, 0.00, 3.58, 7.44e-4, 0.00)
385  Fe_EXPAND(0.00, 0.00, 0.00) };
386 double chtrHe_rec_b2[] = { 0.00, 0.00, 0.00
387  C_EXPAND(0.00, 0.00, 0.42, -7.37e-2, 0.00)
388  N_EXPAND(0.00, 0.92, 0.23, 1.55, 0.00)
389  O_EXPAND(0.00, 2.60, 0.42, 0.40, 0.00)
390  Ne_EXPAND(0.00, 3.35, 0.00, 0.14, 0.00)
391  S_EXPAND(0.00, 0.00, 7.77e-3, 0.34, 0.00)
392  Fe_EXPAND(0.00, 0.00, 0.00) };
393 double chtrHe_rec_c2[] = { 0.00, 0.00, 0.00
394  C_EXPAND(0.00, 0.00, -0.69, 3.50e+1, 0.00)
395  N_EXPAND(0.00, 2.37, -0.72, 11.2, 0.00)
396  O_EXPAND(0.00, 8.99, -0.71, -0.46, 0.00)
397  Ne_EXPAND(0.00, -1.92, 0.00, 4.88e-2, 0.00)
398  S_EXPAND(0.00, 0.00, -0.94, 3.74, 0.00)
399  Fe_EXPAND(0.00, 0.00, 0.00) };
400 double chtrHe_rec_d2[] = { 0.00, 0.00, 0.00
401  C_EXPAND(0.00, 0.00, -0.34, 2.40, 0.00)
402  N_EXPAND(0.00, -10.2, -0.19, -7.82, 0.00)
403  O_EXPAND(0.00, -0.78, -1.98e-2, -0.35, 0.00)
404  Ne_EXPAND(0.00, -1.50, 0.00, -3.35, 0.00)
405  S_EXPAND(0.00, 0.00, -0.30, -5.18, 0.00)
406  Fe_EXPAND(0.00, 0.00, 0.00) };
407 /* 3rd stage - 10,000K to 50,000K */
408 double chtrHe_rec_a3[] = { 0.00, 0.00, 0.00
409  C_EXPAND(0.00, 0.00, 1.12, 1.49e-5, 0.00)
410  N_EXPAND(0.00, 4.84e-1, 2.05, 1.26e-2, 0.00)
411  O_EXPAND(0.00, 7.10e-3, 1.12, 0.997, 0.00)
412  Ne_EXPAND(0.00, 2.52e-2, 1.34e-4, 1.77, 0.00)
413  S_EXPAND(0.00, 0.00, 3.58, 7.44e-4, 0.00)
414  Fe_EXPAND(0.00, 0.00, 0.00) };
415 double chtrHe_rec_b3[] = { 0.00, 0.00, 0.00
416  C_EXPAND(0.00, 0.00, 0.42, 2.73, 0.00)
417  N_EXPAND(0.00, 0.92, 0.23, 1.55, 0.00)
418  O_EXPAND(0.00, 2.60, 0.42, 0.40, 0.00)
419  Ne_EXPAND(0.00, 0.14, 2.33, 0.14, 0.00)
420  S_EXPAND(0.00, 0.00, 7.77e-3, 0.34, 0.00)
421  Fe_EXPAND(0.00, 0.00, 0.00) };
422 double chtrHe_rec_c3[] = { 0.00, 0.00, 0.00
423  C_EXPAND(0.00, 0.00, -0.69, 5.93, 0.00)
424  N_EXPAND(0.00, 2.37, -0.72, 11.2, 0.00)
425  O_EXPAND(0.00, 8.99, -0.71, -0.46, 0.00)
426  Ne_EXPAND(0.00, -1.99, -2.55, 4.88e-2, 0.00)
427  S_EXPAND(0.00, 0.00, -0.94, 3.74, 0.00)
428  Fe_EXPAND(0.00, 0.00, 0.00) };
429 double chtrHe_rec_d3[] = { 0.00, 0.00, 0.00
430  C_EXPAND(0.00, 0.00, -0.34, -8.74e-2, 0.00)
431  N_EXPAND(0.00, -10.2, -0.19, -7.82, 0.00)
432  O_EXPAND(0.00, -0.78, -1.98e-2, -0.35, 0.00)
433  Ne_EXPAND(0.00, -0.91, -0.37, -3.35, 0.00)
434  S_EXPAND(0.00, 0.00, -0.30, -5.18, 0.00)
435  Fe_EXPAND(0.00, 0.00, 0.00) };
436 /* 4th stage - 50,000K to 100,000K */
437 double chtrHe_rec_a4[] = { 0.00, 0.00, 0.00
438  C_EXPAND(0.00, 0.00, 1.12, 1.49e-5, 0.00)
439  N_EXPAND(0.00, 3.17, 2.05, 1.26e-2, 0.00)
440  O_EXPAND(0.00, 6.21e-1, 1.12, 0.997, 0.00)
441  Ne_EXPAND(0.00, 2.52e-2, 1.34e-4, 2.67e-1, 0.00)
442  S_EXPAND(0.00, 0.00, 0.00, 0.00, 0.00)
443  Fe_EXPAND(0.00, 0.00, 0.00) };
444 double chtrHe_rec_b4[] = { 0.00, 0.00, 0.00
445  C_EXPAND(0.00, 0.00, 0.42, 2.73, 0.00)
446  N_EXPAND(0.00, 0.20, 0.23, 1.55, 0.00)
447  O_EXPAND(0.00, 0.53, 0.42, 0.40, 0.00)
448  Ne_EXPAND(0.00, 0.14, 2.33, 0.54, 0.00)
449  S_EXPAND(0.00, 0.00, 0.00, 0.00, 0.00)
450  Fe_EXPAND(0.00, 0.00, 0.00) };
451 double chtrHe_rec_c4[] = { 0.00, 0.00, 0.00
452  C_EXPAND(0.00, 0.00, -0.69, 5.93, 0.00)
453  N_EXPAND(0.00, -0.72, -0.72, 11.2, 0.00)
454  O_EXPAND(0.00, -0.66, -0.71, -0.46, 0.00)
455  Ne_EXPAND(0.00, -1.99, -2.55, 0.91, 0.00)
456  S_EXPAND(0.00, 0.00, 0.00, 0.00, 0.00)
457  Fe_EXPAND(0.00, 0.00, 0.00) };
458 double chtrHe_rec_d4[] = { 0.00, 0.00, 0.00
459  C_EXPAND(0.00, 0.00, -0.34, -8.74e-2, 0.00)
460  N_EXPAND(0.00, -4.81e-2, -0.19, -7.82, 0.00)
461  O_EXPAND(0.00, -2.22e-2, -1.98e-2, -0.35, 0.00)
462  Ne_EXPAND(0.00, -0.91, -0.37, -1.88e-2, 0.00)
463  S_EXPAND(0.00, 0.00, 0.00, 0.00, 0.00)
464  Fe_EXPAND(0.00, 0.00, 0.00) };
465 /* 5th stage - more than 100,000K */
466 double chtrHe_rec_a5[] = { 0.00, 0.00, 0.00
467  C_EXPAND(0.00, 0.00, 1.12, 1.49e-5, 0.00)
468  N_EXPAND(0.00, 3.17, 2.05, 3.75e-1, 0.00)
469  O_EXPAND(0.00, 6.21e-1, 1.12, 0.997, 0.00)
470  Ne_EXPAND(0.00, 2.52e-2, 0.1, 2.67e-1, 0.00)
471  S_EXPAND(0.00, 0.00, 0.00, 0.00, 0.00)
472  Fe_EXPAND(0.00, 0.00, 0.00) };
473 double chtrHe_rec_b5[] = { 0.00, 0.00, 0.00
474  C_EXPAND(0.00, 0.00, 0.42, 2.73, 0.00)
475  N_EXPAND(0.00, 0.20, 0.23, 0.54, 0.00)
476  O_EXPAND(0.00, 0.53, 0.42, 0.40, 0.00)
477  Ne_EXPAND(0.00, 0.14, 0.24, 0.54, 0.00)
478  S_EXPAND(0.00, 0.00, 0.00, 0.00, 0.00)
479  Fe_EXPAND(0.00, 0.00, 0.00) };
480 double chtrHe_rec_c5[] = { 0.00, 0.00, 0.00
481  C_EXPAND(0.00, 0.00, -0.69, 5.93, 0.00)
482  N_EXPAND(0.00, -0.72, -0.72, -0.82, 0.00)
483  O_EXPAND(0.00, -0.66, -0.71, -0.46, 0.00)
484  Ne_EXPAND(0.00, -1.99, -1.09, 0.91, 0.00)
485  S_EXPAND(0.00, 0.00, 0.00, 0.00, 0.00)
486  Fe_EXPAND(0.00, 0.00, 0.00) };
487 double chtrHe_rec_d5[] = { 0.00, 0.00, 0.00
488  C_EXPAND(0.00, 0.00, -0.34, -8.74e-2, 0.00)
489  N_EXPAND(0.00, -4.81e-2, -0.19, -2.07e-2, 0.00)
490  O_EXPAND(0.00, -2.22e-2, -1.98e-2, -0.35, 0.00)
491  Ne_EXPAND(0.00, -0.91, -2.47e-2, -1.88e-2, 0.00)
492  S_EXPAND(0.00, 0.00, 0.00, 0.00, 0.00)
493  Fe_EXPAND(0.00, 0.00, 0.00) };
494 
495 /* Radiative recombination - Pequignot & al 1991, A&A */
496 double rad_rec_a[] = { 5.596, 0.0000, 0.0000
497  C_EXPAND(5.068, 5.434, 4.742, 4.051, 0.00)
498  N_EXPAND(3.874, 4.974, 4.750, 4.626, 0.00)
499  O_EXPAND(3.201, 4.092, 4.890, 14.665, 0.00)
500  Ne_EXPAND(0.000, 0.000, 0.000, 0.0000, 0.00)
501  S_EXPAND(0.000, 0.000, 0.000, 0.0000, 0.00) }; /* H value NOT used */
502  /* Ne values: 11.800, 5.841, 15.550, 7.538, */
503  /* He I and II values: 8.295, 5.596, */
504 double rad_rec_b[] = { -0.6038, 0.0000, 0.0000
505  C_EXPAND(-0.6192, -0.6116, -0.6167, -0.6270, 0.00)
506  N_EXPAND(-0.6487, -0.6209, -0.5942, -0.9521, 0.00)
507  O_EXPAND(-0.6880, -0.6413, -0.6213, -0.5140, 0.00)
508  Ne_EXPAND( 0.0000, 0.0000, 0.0000, 0.0000, 0.00)
509  S_EXPAND( 0.0000, 0.0000, 0.0000, 0.0000, 0.00)}; /* S, H and Ne to add separately ; prev. H value: -0.6038 */
510  /* Ne values: -0.5404, -0.5921, -0.4825, -0.5540, */
511  /* He I and II values: -0.6193, -0.6038, */
512 double rad_rec_c[] = { 0.3436, 0.0000, 0.0000
513  C_EXPAND(-0.0815, 0.0694, 0.2960, 0.5054, 0.00)
514  N_EXPAND( 0.0000, 0.0000, 0.8452, 0.4729, 0.00)
515  O_EXPAND(-0.0174, 0.0000, 0.0184, 2.7300, 0.00)
516  Ne_EXPAND( 0.0000, 0.0000, 0.0000, 0.0000, 0.00)
517  S_EXPAND( 0.0000, 0.0000, 0.0000, 0.0000, 0.00)};/* prev. H value: 0.3436*/
518  /* Ne values: 3.0300, 0.4498, 3.2740, 1.2960, */
519  /* He I and II values: 0.9164, 0.3436, */
520 double rad_rec_d[] = { 0.4479, 0.0000, 0.0000
521  C_EXPAND(1.2910, 0.7866, 0.6167, 0.6692, 0.00)
522  N_EXPAND(1.0000, 1.0000, 2.8450, -0.4477, 0.00)
523  O_EXPAND(1.7070, 1.0000, 1.5550, 0.2328, 0.00)
524  Ne_EXPAND(0.0000, 0.0000, 0.0000, 0.0000, 0.00)
525  S_EXPAND(0.0000, 0.0000, 0.0000, 0.0000, 0.00) };/* prev. H value: 0.4479*/
526  /* Ne values: 0.2050, 0.6395, 0.3030, 0.3472, */
527  /* He I and II values: 0.2667, 0.4479, */
528 
529 /* Dielectronic recombination - Nussbaumer, Storey, 1983*/
530 double diel_rec_a[] = { 0.0000, 0.0000, 0.0000
531  C_EXPAND(0.0108, 1.8267, 2.3196, 0.0000, 0.00)
532  N_EXPAND(0.0000, 0.0320, -0.8806, 0.4134, 0.00)
533  O_EXPAND(0.0000, -0.0036, 0.0000, 0.0061, 0.00)
534  Ne_EXPAND(0.0000, 0.0000, 0.0000, 0.0000, 0.00)
535  S_EXPAND(0.0000, 0.0000, 0.0000, 0.0000, 0.00) };
536 double diel_rec_b[] = { 0.0000, 0.0000, 0.0000
537  C_EXPAND(-0.1075, 4.1012, 10.7328, 0.0000, 0.00)
538  N_EXPAND( 0.6310, -0.6624, 11.2406, -4.6319, 0.00)
539  O_EXPAND( 0.0238, 0.7519, 21.8790, 0.2269, 0.00)
540  Ne_EXPAND( 0.0000, 0.0000, 0.0000, 0.0000, 0.00)
541  S_EXPAND( 0.0000, 0.0000, 0.0000, 0.0000, 0.00) };
542 double diel_rec_c[] = { 0.0000, 0.0000, 0.0000
543  C_EXPAND(0.2810, 4.8443, 6.8830, 0.0000, 0.00)
544  N_EXPAND(0.1990, 4.3191, 30.7066, 25.9172, 0.00)
545  O_EXPAND(0.0659, 1.5252, 16.2730, 32.1419, 0.00)
546  Ne_EXPAND(0.0000, 0.0000, 0.0000, 0.0000, 0.00)
547  S_EXPAND(0.0000, 0.0000, 0.0000, 0.0000, 0.00) };
548 double diel_rec_d[] = { 0.0000, 0.0000, 0.0000
549  C_EXPAND(-0.0193, 0.2261, -0.1824, 0.0000, 0.00)
550  N_EXPAND(-0.0197, 0.0003, -1.1721, -2.2290, 0.00)
551  O_EXPAND( 0.0349, -0.0838, -0.7020, 1.9939, 0.00)
552  Ne_EXPAND( 0.0000, 0.0000, 0.0000, 0.0000, 0.00)
553  S_EXPAND( 0.0000, 0.0000, 0.0000, 0.0000, 0.00) };
554 double diel_rec_f[] = { 0.0000, 0.0000, 0.0000
555  C_EXPAND(-0.1127, 0.5960, 0.4101, 0.0000, 0.00)
556  N_EXPAND( 0.4398, 0.5946, 0.6127, 0.2360, 0.00)
557  O_EXPAND( 0.5334, 0.2769, 1.1899, -0.0646, 0.00)
558  Ne_EXPAND( 0.0000, 0.0000, 0.0000, 0.0000, 0.00)
559  S_EXPAND( 0.0000, 0.0000, 0.0000, 0.0000, 0.00) };
560 
561 
562 /* Total recombination rate for H, Ne and S - NIFS-DATA-54 - Kato & Asano 1999 */
563 /* recombination coefficient to the X^n ion (from the X^(n+1) is written in array in the position of X^n */
564 double tot_rec_T[5] = { 5.e+3, 1.1e+4, 5.5e+4, 1.1e+5, 2.2e+5 }; /* in K */
565 static int tot_rec_ions[] = { 0, 1, 18, 19, 20, 21, 23, 24, 25, 26};
566 double tot_rec[][5] = { { -12.5, -13.0, -13.2, -13.4, -13.65 }, /* H */
567  { -12.5, -12.6, -12.8, -12.1, -11.75 }, /* He I */
568  /* { -11.4, -11.7, -12.2, -12.45, -12.8 }, */
569  { -12.3, -12.5, -12.6, -11.8, -11.7 }, /* Ne I */
570  { -11.7, -11.9, -11.6, -11.3, -11.3 },
571  { -11.2, -11.4, -11.6, -10.9, -11.1 },
572  { -11.0, -11.1, -11.3, -10.6, -10.7 },
573  { -12.5, -12.4, -10.9, -10.95, -11.4 }, /* S I */
574  { -11.7, -11.8, -10.6, -10.3, -10.6 },
575  { -11.4, -11.6, -10.0, -9.8, -10.1 },
576  { -11.2, -11.3, -9.9, -9.7, -10.0 } }; /* log (Alpha_tot) [cm^3/s]*/
577 
578 double tot_recH_He[2][5] = {{ -12.5, -13.0, -13.2, -13.4, -13.65 }, /* -- H -- */
579  { -12.5, -12.6, -12.8, -12.1, -11.75 }}; /* -- He -- */
580 double tot_recNe[][5] = { { -12.3, -12.5, -12.6, -11.8, -11.7 },
581  { -11.7, -11.9, -11.6, -11.3, -11.3 },
582  { -11.2, -11.4, -11.6, -10.9, -11.1 },
583  { -11.0, -11.1, -11.3, -10.6, -10.7 }};
584 double tot_recS[][5] = { { -12.5, -12.4, -10.9, -10.95, -11.4 },
585  { -11.7, -11.8, -10.6, -10.3, -10.6 },
586  { -11.4, -11.6, -10.0, -9.8, -10.1 },
587  { -11.2, -11.3, -9.9, -9.7, -10.0 }};
588 
589 
590 /* Atomic data for Fe */
591 /* Radiative recombination, Arnaud & Raymond 1992, ApJ 398, 394 */
592 double fe_rr_A [3] = { 1.42e-13, 1.02e-12, 0.00 };
593 double fe_rr_eta [3] = { 0.891, 0.843, 0.00 };
594 /* Dielectronic recombination, Arnaud & Raymond 1992, ApJ 398, 394 */
595 double fe_dr_e1 [3] = { 1.670, 2.860, 0.000};
596 double fe_dr_e2 [3] = { 31.40, 52.1, 0.00};
597 double fe_dr_e3 [3] = { 0.00, 0.00, 0.00};
598 double fe_dr_e4 [3] = { 0.00, 0.00, 0.00};
599 double fe_dr_c1 [3] = { 2.30e-3, 1.50e-2, 0.00};
600 double fe_dr_c2 [3] = { 2.7e-3, 4.7e-3, 0.00};
601 double fe_dr_c3 [3] = { 0.00, 0.00, 0.00};
602 double fe_dr_c4 [3] = { 0.00, 0.00, 0.00};
603 /* Charge transfer with H - above */
604 /* Charge transfer with H+ - above */
605 /* Charge transfer with He - not relevant */
606 
607  print1("> MINeq: creating ionization coefficients tables...\n");
608  print1(" * collisional ionization\n");
609 
610  if (coll_ion == NULL) coll_ion = ARRAY_2D(I_g_stepNumber, NIONS,double);
611 
612  for (i = 0; i < I_g_stepNumber; i++ ) {
613  T = I_TBEG + (float)i*I_TSTEP;
614  for (nv = 0; nv < NIONS; nv++) {
615  if ( coll_ion_Tmin[nv] < to_ev(T) )
616  coll_ion[i][nv] = coll_ion_A[nv] * ( 1. + coll_ion_P[nv] * sqrt(coll_ion_dE[nv]/to_ev(T)) )
617  / ( coll_ion_X[nv] + ( coll_ion_dE[nv]/to_ev(T) ) )
618  * pow( coll_ion_dE[nv]/to_ev(T), coll_ion_K[nv] )
619  * exp ( -coll_ion_dE[nv]/to_ev(T) );
620  else coll_ion[i][nv] = coll_ion_A[nv] * ( 1. + coll_ion_P[nv] * sqrt(coll_ion_dE[nv]/to_ev(T)) )
621  / ( coll_ion_X[nv] + ( coll_ion_dE[nv]/to_ev(T) ) )
622  * pow( coll_ion_dE[nv]/to_ev(T), coll_ion_K[nv] )
623  * exp ( -coll_ion_dE[nv]/to_ev(T) ) * exp( -(coll_ion_Tmin[nv]-to_ev(T))*1.e+2 );
624  }
625  }
626 
627  print1(" * radiative recombination\n");
628  if (rad_rec == NULL) rad_rec = ARRAY_2D(I_g_stepNumber, NIONS,double);
629 
630  for (i = 0; i < I_g_stepNumber; i++ ) {
631  T = I_TBEG + i*I_TSTEP;
632  for (nv = 0; nv < NIONS - Fe_IONS; nv++) {
633  if ( (T*1.e-4/(rad_rec_z[nv]*rad_rec_z[nv])>0.004) && (T*1.e-4/(rad_rec_z[nv]*rad_rec_z[nv])<4.0) )
634  rad_rec[i][nv] = 1.e-13 * rad_rec_z[nv] * rad_rec_a[nv] * pow( 1.e-4*T/pow(rad_rec_z[nv],2.), rad_rec_b[nv] )
635  / ( 1. + rad_rec_c[nv] * pow( 1.e-4*T/pow(rad_rec_z[nv],2.), rad_rec_d[nv] ) ) ;
636  else {
637  if (T*1.e-4/(rad_rec_z[nv]*rad_rec_z[nv])>=4.0 && (T*1.e-4/(rad_rec_z[nv]*rad_rec_z[nv])<5.0)) {
638  rad_rec[i][nv] = 1.e-13 * rad_rec_z[nv] * rad_rec_a[nv] * pow( 1.e-4*T/pow(rad_rec_z[nv],2.), rad_rec_b[nv] )
639  / ( 1. + rad_rec_c[nv] * pow( 1.e-4*T/pow(rad_rec_z[nv],2.), rad_rec_d[nv] ) )
640  * exp((4.0 - T*1.e-4/(rad_rec_z[nv]*rad_rec_z[nv]))*5.0);
641  }
642  else rad_rec[i][nv] = 0.0;
643  }
644  }
645  lam = 157890. / T;
646  rad_rec[i][0] = 5.197e-14 * sqrt(lam) * (0.4288 + 0.5 * log(lam) + 0.469 * pow(lam,-1./3.) );
647 #if Fe_IONS > 0
648  for (nv = NIONS-Fe_IONS; nv < NIONS; nv++) /* data for Fe ions */
649  rad_rec[i][nv] = fe_rr_A[nv - NIONS + 3] * pow( (T*1.e-4), -fe_rr_eta[nv-NIONS+3]);
650 #endif
651  }
652 
653  print1(" * dielectronic recombination\n");
654  if (diel_rec == NULL) diel_rec = ARRAY_2D(I_g_stepNumber,NIONS,double);
655 
656  for (i = 0; i < I_g_stepNumber; i++ ) {
657  T = I_TBEG + i*I_TSTEP;
658  t4 = T/10000.0;
659  for (nv = 0; nv < NIONS - Fe_IONS; nv++) {
660  diel_rec[i][nv] = 1.e-12 * ( diel_rec_a[nv] / t4 + diel_rec_b[nv] + diel_rec_c[nv]*t4 + diel_rec_d[nv] * pow(t4,2.) )
661  * pow( t4, -3./2. ) * exp( -diel_rec_f[nv]/t4 );
662  }
663 #if Fe_IONS > 0
664  for (nv = NIONS-Fe_IONS; nv < NIONS; nv++) {
665  diel_rec[i][nv] = 1.6e-12 * pow(T, -1.5) * ( fe_dr_c1[nv-NIONS+3]*exp(-fe_dr_e1[nv-NIONS+3]/to_ev(T))
666  + fe_dr_c2[nv-NIONS+3]*exp(-fe_dr_e2[nv-NIONS+3]/to_ev(T))
667  + fe_dr_c3[nv-NIONS+3]*exp(-fe_dr_e3[nv-NIONS+3]/to_ev(T))
668  + fe_dr_c4[nv-NIONS+3]*exp(-fe_dr_e4[nv-NIONS+3]/to_ev(T)) );
669  }
670 #endif
671  }
672 
673  print1(" * charge transfer with H+\n");
674  if (chtr_hp == NULL) chtr_hp = ARRAY_2D(I_g_stepNumber,NIONS,double);
675 
676  for (i = 0; i < I_g_stepNumber; i++ ) {
677  T = I_TBEG + i*I_TSTEP;
678  t4 = T/10000.0;
679  for (nv = 0; nv < NIONS; nv++) {
680  if (t4<=chtrH_ion_upT[nv])
681  chtr_hp[i][nv] = 1.e-9 * chtrH_ion_a[nv] * pow( t4, chtrH_ion_b[nv] )
682  * ( 1. + chtrH_ion_c[nv] * exp ( chtrH_ion_d[nv]*t4 ));
683  else chtr_hp[i][nv] = 1.e-9 * chtrH_ion_a[nv] * pow( t4, chtrH_ion_b[nv] )
684  * ( 1. + chtrH_ion_c[nv] * exp ( chtrH_ion_d[nv]*t4 ))
685  * exp( -(t4-chtrH_ion_upT[nv])/(5.e-2*chtrH_ion_upT[nv])); /* smooth to 0 */
686  }
687  }
688 
689  print1(" * charge transfer with H\n");
690  if (chtr_h == NULL) chtr_h = ARRAY_2D(I_g_stepNumber,NIONS,double);
691 
692  for (i = 0; i < I_g_stepNumber; i++ ) {
693  T = I_TBEG + i*I_TSTEP;
694  t4 = T/10000.0;
695  for (nv = 0; nv < NIONS; nv++) {
696  if (T<=chtrH_rec_upT[nv])
697  chtr_h[i][nv] = 1.e-9 * chtrH_rec_a[nv] * pow( t4, chtrH_rec_b[nv] )
698  * ( 1. + chtrH_rec_c[nv] * exp ( chtrH_rec_d[nv]*t4 ));
699  else chtr_h[i][nv] = 1.e-9 * chtrH_rec_a[nv] * pow( t4, chtrH_rec_b[nv] )
700  * ( 1. + chtrH_rec_c[nv] * exp ( chtrH_rec_d[nv]*t4 ))
701  * exp( -(T - chtrH_rec_upT[nv])/(5.e-2*chtrH_rec_upT[nv])); /* smooth to 0 */
702  }
703  }
704 
705  print1(" * charge transfer with He\n");
706  if (chtr_he == NULL) chtr_he = ARRAY_2D(I_g_stepNumber,NIONS,double);
707 
708  for (i = 0; i < I_g_stepNumber; i++ ) {
709  T = I_TBEG + i*I_TSTEP;
710  t4 = T/10000.0;
711  for (nv = 0; nv < NIONS; nv++) {
712 /*
713  if (T<5000.)
714  chtr_he[i][nv] = elem_ab[el_He] * 1.e-9 * chtrHe_rec_a1[nv] * pow( t4, chtrHe_rec_b1[nv] )
715  * ( 1. + chtrHe_rec_c1[nv] * exp ( chtrHe_rec_d1[nv]*t4 ));
716  if ( (T>=5000.) && (T<10000.) )
717  chtr_he[i][nv] = elem_ab[el_He] * 1.e-9 * chtrHe_rec_a2[nv] * pow( t4, chtrHe_rec_b2[nv] )
718  * ( 1. + chtrHe_rec_c2[nv] * exp ( chtrHe_rec_d2[nv]*t4 ));
719  if ( (T>=10000.) && (T<50000.) )
720  chtr_he[i][nv] = elem_ab[el_He] * 1.e-9 * chtrHe_rec_a3[nv] * pow( t4, chtrHe_rec_b3[nv] )
721  * ( 1. + chtrHe_rec_c3[nv] * exp ( chtrHe_rec_d3[nv]*t4 ));
722  if ( (T>=50000.) && (T<100000.) )
723  chtr_he[i][nv] = elem_ab[el_He] * 1.e-9 * chtrHe_rec_a4[nv] * pow( t4, chtrHe_rec_b4[nv] )
724  * ( 1. + chtrHe_rec_c4[nv] * exp ( chtrHe_rec_d4[nv]*t4 ));
725  if (T>=100000.)
726  chtr_he[i][nv] = elem_ab[el_He] * 1.e-9 * chtrHe_rec_a5[nv] * pow( t4, chtrHe_rec_b5[nv] )
727  * ( 1. + chtrHe_rec_c5[nv] * exp ( chtrHe_rec_d5[nv]*t4 ));
728 */
729 /* Replaced temperature intervals with the smooth formula below */
730  chtr_he[i][nv] = elem_ab[el_He] * 1.e-9 * chtrHe_rec_a1[nv] * pow( t4, chtrHe_rec_b1[nv] )
731  * ( 1. + chtrHe_rec_c1[nv] * exp ( chtrHe_rec_d1[nv]*t4 ))
732  * (T>5000.0?exp( -(T-5000.0)/1.e2):1.0)
733  + elem_ab[el_He] * 1.e-9 * chtrHe_rec_a2[nv] * pow( t4, chtrHe_rec_b2[nv] )
734  * ( 1. + chtrHe_rec_c2[nv] * exp ( chtrHe_rec_d2[nv]*t4 ))
735  * (T>10000.0?exp( -(T-10000.0)/2.e2):1.0) * (T<5000.0?exp( -(5000.0-T)/1.e2):1.0)
736  + elem_ab[el_He] * 1.e-9 * chtrHe_rec_a3[nv] * pow( t4, chtrHe_rec_b3[nv] )
737  * ( 1. + chtrHe_rec_c3[nv] * exp ( chtrHe_rec_d3[nv]*t4 ))
738  * (T>50000.0?exp( -(T-50000.0)/1.e3):1.0) * (T<10000.0?exp( -(10000.0-T)/2.e2):1.0)
739  + elem_ab[el_He] * 1.e-9 * chtrHe_rec_a4[nv] * pow( t4, chtrHe_rec_b4[nv] )
740  * ( 1. + chtrHe_rec_c4[nv] * exp ( chtrHe_rec_d4[nv]*t4 ))
741  * (T>100000.0?exp( -(T-100000.0)/2.e3):1.0) * (T<50000.0?exp( -(50000.0-T)/1.e3):1.0)
742  + elem_ab[el_He] * 1.e-9 * chtrHe_rec_a5[nv] * pow( t4, chtrHe_rec_b5[nv] )
743  * ( 1. + chtrHe_rec_c5[nv] * exp ( chtrHe_rec_d5[nv]*t4 ))
744  * (T<100000.0?exp( -(100000.0-T)/2.e3):1.0);
745  }
746  }
747 
748  /* And now, for ions for which we only have the total electron-ion recombination coefficient: */
749 
750  print1(" * total recombination coefficients\n");
751  if (tot_recs == NULL) tot_recs = ARRAY_2D(I_g_stepNumber,NIONS,double);
752  for (i = 0; i < I_g_stepNumber; i++ )
753  for (nv = 0; nv < NIONS; nv++ ) tot_recs[i][nv] = 0.0;
754 
755  for (i = 0; i < I_g_stepNumber; i++ ) {
756  T = (double) (I_TBEG + i*I_TSTEP);
757  if (T <= tot_rec_T[0]) { ti1=-1; ti2=0; } /* Find temperature vector indexes to interpolate - ti1 and ti2 */
758  if (T > tot_rec_T[4]) { ti1=4; ti2=-1; }
759  for (j = 0; j < 4; j++) {
760  if ( (tot_rec_T[j] < T) && (tot_rec_T[j+1]>=T) ) { ti1=j; ti2=j+1; }
761  }
762 
763  /* -- H/He recombination -- */
764 
765  for (nv = 0; nv <= 1; nv++){
766  if (ti1 == -1) tmprec = tot_recH_He[nv][0];
767  else if (ti2 == -1) tmprec = tot_recH_He[nv][4];
768  else tmprec = (T-tot_rec_T[ti1])/(tot_rec_T[ti2]-tot_rec_T[ti1])*tot_recH_He[nv][ti2]
769  + (tot_rec_T[ti2]-T)/(tot_rec_T[ti2]-tot_rec_T[ti1])*tot_recH_He[nv][ti1];
770  tot_recs[i][nv] = pow(10.0,tmprec);
771  }
772 
773  /* -- Ne recombination -- */
774 
775  for (nv = 0; nv < Ne_IONS-1; nv++){
776  if (ti1 == -1) tmprec = tot_recNe[nv][0];
777  else if (ti2 == -1) tmprec = tot_recNe[nv][4];
778  else tmprec = (T-tot_rec_T[ti1])/(tot_rec_T[ti2]-tot_rec_T[ti1])*tot_recNe[nv][ti2]
779  + (tot_rec_T[ti2]-T)/(tot_rec_T[ti2]-tot_rec_T[ti1])*tot_recNe[nv][ti1];
780  tot_recs[i][X_NeI+nv-NFLX] = pow(10.0,tmprec);
781  }
782 
783  /* -- S recombination -- */
784 
785  for (nv = 0; nv < S_IONS-1; nv++){
786  if (ti1 == -1) tmprec = tot_recS[nv][0];
787  else if (ti2 == -1) tmprec = tot_recS[nv][4];
788  else tmprec = (T-tot_rec_T[ti1])/(tot_rec_T[ti2]-tot_rec_T[ti1])*tot_recS[nv][ti2]
789  + (tot_rec_T[ti2]-T)/(tot_rec_T[ti2]-tot_rec_T[ti1])*tot_recS[nv][ti1];
790  tot_recs[i][X_SI+nv-NFLX] = pow(10.0,tmprec);
791  }
792 continue;
793 
794  for (nv = 0; nv < NIONS; nv++) {
795  cindex=-1;
796  for (j = 0; j < 10; j++)
797  if (tot_rec_ions[j] == nv) cindex = j; /* Find the ion index in the total rec. array */
798 
799  if (cindex >= 0) {
800  tmprec=0.; /* make a linear interpolation... */
801  if (ti1 == -1) tmprec = tot_rec[cindex][0];
802  else if (ti2 == -1) tmprec = tot_rec[cindex][4];
803  else tmprec = (T-tot_rec_T[ti1])/(tot_rec_T[ti2]-tot_rec_T[ti1])*tot_rec[cindex][ti2]
804  + (tot_rec_T[ti2]-T)/(tot_rec_T[ti2]-tot_rec_T[ti1])*tot_rec[cindex][ti1];
805  tot_recs[i][nv] = pow(10.0,tmprec); /* and obtain the rate * N(x+1) */
806 /*
807  if (tmprec) tot_recs[i][nv] = pow(10.0,tmprec);
808  else tot_recs[i][nv] = 0.0;
809 */
810  }
811  }
812  }
813 
814  /* -- generate table -- */
815 
816  for (i = 0; i < I_g_stepNumber; i++){
817  for (nv = 0; nv < NIONS; nv++) {
818  tbl[i][0][nv] = coll_ion[i][nv];
819  tbl[i][1][nv] = rad_rec[i][nv];
820  tbl[i][2][nv] = diel_rec[i][nv];
821  tbl[i][3][nv] = chtr_hp[i][nv];
822  tbl[i][4][nv] = chtr_h[i][nv];
823  tbl[i][5][nv] = chtr_he[i][nv];
824  tbl[i][6][nv] = tot_recs[i][nv];
825  }
826 
827  for (j = 0; j < 7; j++) {
828  tbl[i][j][X_CI + C_IONS-1-NFLX] = 0.0;
829  tbl[i][j][X_NI + N_IONS-1-NFLX] = 0.0;
830  tbl[i][j][X_OI + O_IONS-1-NFLX] = 0.0;
831  tbl[i][j][X_NeI + Ne_IONS-1-NFLX] = 0.0;
832  tbl[i][j][X_SI + S_IONS-1-NFLX] = 0.0;
833  }
834  }
835  print1(" * Done.\n");
836  return 0;
837 }
838 #undef ZERO_5X
839 
840 /* ********************************************************************* */
841 int Create_Losses_Tables(double ***losstables, int *nT, int *nN)
842 /*!
843  * Compute and save to disk the radiative
844  * losses tables function of Ne and T
845  *
846  *********************************************************************** */
847 {
848  int ib, ie, first_time, atom_id, i, j;
849  double Ne, T, erg;
850  double tmpN, tmpT;
851  double line_1, line_2, d;
852  char fname[100],fnn[10];
853  Ion atoms[NIONS], *X;
854 
855  erg = 1.602177e-12;
856 
857  Ne = C_NeMIN;
858  *nN = 0;
859  while (Ne < C_NeMAX) {
860  T = C_TMIN;
861  *nT = 0;
862  while (T < C_TMAX) {
863  tmpT = T;
864  T = T*exp(C_TSTEP); /* should be *exp(0.02) */
865  *nT = *nT + 1;
866  }
867  tmpN = Ne;
868  Ne = Ne*exp(C_NeSTEP); /* should be *exp(0.06) */
869  *nN = *nN + 1;
870  }
871 
872  for (i = 0; i < NIONS; i++) atoms[i].dE = NULL;
873 
874  first_time = 1;
875  for (atom_id = 0; atom_id < NIONS; atom_id++) {
876  X = &atoms[atom_id];
877  INIT_ATOM(X,atom_id);
878  Ne = C_NeMIN;
879  *nN = 0;
880  while (Ne < C_NeMAX) {
881  T = C_TMIN;
882  *nT = 0;
883  while (T < C_TMAX) {
884 
885  line_2 = 0.0;
886  Solve_System (X, Ne, T);
887  for (ib = 1; ib < X->nlev; ib++) {
888  for (ie = 0; ie < ib; ie++) {
889  if ( (X->Ni[ib] != X->Ni[ib]) || (X->A[ib][ie] != X->A[ib][ie]) ||
890  (X->dE[ib][ie] != X->dE[ib][ie]) ||
891  (X->Ni[ib] * X->A[ib][ie] * X->dE[ib][ie] > 1.0e+200) ){
892  printf (" Nan found atom = %d ib,ie = %d, %d\n", atom_id,ib,ie);
893  printf (" T, ne = %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n",T, Ne,X->Ni[ib],X->A[ib][ie],X->dE[ib][ie], X->omega[ib][ie][0] );
894  QUIT_PLUTO(1);
895  }
896 
897  if (X->A[ib][ie] > 0.0 && X->Ni[ib] > 0.0)
898  line_2 += X->Ni[ib] * X->A[ib][ie] * X->dE[ib][ie] / Ne ;
899  /* cooling / ion in eV, should be multiplied by ion abundance and Ne to obtain energy/cmc/s*/
900  }}
901  line_2 = line_2*erg;
902 
903  losstables[atom_id][*nN][*nT] = line_2;
904  tmpT = T;
905  T = T*exp(C_TSTEP); /* should be *exp(0.02) */
906  *nT = *nT + 1;
907  }
908  first_time = 0;
909  tmpN = Ne;
910  Ne = Ne*exp(C_NeSTEP); /* should be *exp(0.06) */
911  *nN = *nN + 1;
912  }
913  }
914 
915  print1("> MINEq radiative losses tables generated and saved to memory.\n");
916  return(0);
917 }
Solve_System
void Solve_System(Ion *X, double Ne, double T)
Definition: make_tables.c:5
to_ev
#define to_ev(T)
Definition: cooling_defs.h:65
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
I_TBEG
#define I_TBEG
Definition: cooling_defs.h:171
S_IONS
#define S_IONS
Definition: cooling.h:17
Symmetrize_Coeff
void Symmetrize_Coeff(Ion *X)
Definition: make_tables.c:97
n
static int n
Definition: analysis.c:3
ION::nlev
int nlev
Definition: cooling_defs.h:112
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
ION::isCHEB
int isCHEB
Definition: cooling_defs.h:112
elem_ab
const double elem_ab[8]
Definition: radiat.c:9
configure.scrh
tuple scrh
Definition: configure.py:200
rad_rec_z
const double rad_rec_z[31]
Definition: radiat.c:24
C_TSTEP
#define C_TSTEP
Definition: cooling_defs.h:185
C_IONS
#define C_IONS
Definition: cooling.h:13
Ne_EXPAND
Definition: cooling.h:114
ION::nTom
int nTom
Definition: cooling_defs.h:112
ION::Tom
double Tom[8]
Definition: cooling_defs.h:117
lagrange
double lagrange(double *x, double *y, double xp, int n, int ii, int jj)
Definition: make_tables.c:161
el_He
#define el_He
Definition: cooling_defs.h:72
S_EXPAND
#define S_EXPAND(a, b, c, d, e)
Definition: cooling.h:83
coll_ion_dE
const double coll_ion_dE[31]
Definition: radiat.c:31
N_EXPAND
#define N_EXPAND(a, b, c, d, e)
Definition: cooling.h:41
ION::isH
int isH
Definition: cooling_defs.h:112
O_IONS
#define O_IONS
Definition: cooling.h:15
I_TSTEP
#define I_TSTEP
Definition: cooling_defs.h:170
cooling_defs.h
LUBackSubst
void LUBackSubst(double **a, int n, int *indx, double b[])
Definition: math_lu_decomp.c:84
C_TMAX
#define C_TMAX
Definition: cooling_defs.h:184
X
#define X
Definition: radiat.modified_by_Bhargav.c:3
NFLX
#define NFLX
Definition: mod_defs.h:32
NIONS
#define NIONS
Definition: cooling.h:28
C_EXPAND
#define C_EXPAND(a, b, c, d, e)
Definition: cooling.h:27
C_NeSTEP
#define C_NeSTEP
Definition: cooling_defs.h:181
ION::dE
double ** dE
Definition: cooling_defs.h:116
ION::Ni
double Ni[16]
Definition: cooling_defs.h:115
I_g_stepNumber
#define I_g_stepNumber
Definition: cooling_defs.h:169
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
pluto.h
PLUTO main header file.
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
Ne_IONS
#define Ne_IONS
Definition: cooling.h:16
sod.f1
tuple f1
Definition: sod.py:13
ION::A
double ** A
Definition: cooling_defs.h:116
C_TMIN
#define C_TMIN
Definition: cooling_defs.h:183
i
int i
Definition: analysis.c:2
ION::isMAP
int isMAP
Definition: cooling_defs.h:112
Create_Ion_Coeff_Tables
int Create_Ion_Coeff_Tables(double ***tbl)
Definition: make_tables.c:200
O_EXPAND
#define O_EXPAND(a, b, c, d, e)
Definition: cooling.h:55
INIT_ATOM
void INIT_ATOM(Ion *, int)
Definition: ion_init.c:25
Fe_IONS
#define Fe_IONS
Definition: cooling.h:18
ION::isCV
int isCV
Definition: cooling_defs.h:112
C_NeMAX
#define C_NeMAX
Definition: cooling_defs.h:180
kB
#define kB
Definition: cooling_defs.h:167
ARRAY_2D
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
ION::omega
double *** omega
Definition: cooling_defs.h:117
LUDecompose
int LUDecompose(double **a, int n, int *indx, double *d)
Definition: math_lu_decomp.c:12
ION::wght
double wght[16]
Definition: cooling_defs.h:114
Create_Losses_Tables
int Create_Losses_Tables(double ***losstables, int *nT, int *nN)
Definition: make_tables.c:841
q
static Runtime q
Definition: runtime_setup.c:490
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
Fe_EXPAND
#define Fe_EXPAND(a, b, c)
Definition: cooling.h:91
N_IONS
#define N_IONS
Definition: cooling.h:14
C_NeMIN
#define C_NeMIN
Definition: cooling_defs.h:179