PLUTO
radiat.c
Go to the documentation of this file.
1 #include "pluto.h"
2 #include "cooling_defs.h"
3 
4 /* -----------------------------------------------------------------------
5  Global variables definition. They are declared inside cooling_defs.h
6  ----------------------------------------------------------------------- */
7 
8 /* H He C N O Ne S Fe */
9 const double elem_ab[] = { 0.93, 0.074, 3.0e-4, 5.e-5, 4.0e-4, 7.0e-5, 1.5e-5, 2.69e-5 }; /* Number densities */
10 double elem_ab_est[] = { 0.93, 0.074, 3.0e-4, 5.e-5, 4.0e-4, 7.0e-5, 1.5e-5 }; /* Number densities in Orion nebula, Esteban & al 2004 */
11 double elem_ab_lod[] = { 0.92, 0.08, 2.3e-4, 6.2e-5, 4.4e-4, 7.0e-5, 1.25e-5 }; /* Number densities, Solar, Lodders 2003 ApJ */
12 double elem_ab_sol[] = { 0.93, 0.074, 3.0e-4, 9.e-5, 7.e-4, 7.e-5, 1.e-5 }; /* Number fractions, Solar */
13 double elem_ab_uni[] = { 0.93, 0.072, 5.0e-4, 9.e-5, 8.e-4, 8.e-5, 2.e-5 }; /* Number fractions, Universe */
14 const double elem_mass[] = { 1.007, 4.002, 12.01, 14.01, 15.99, 20.18, 32.07, 55.845 }; /* Atomic mass, in a.m.u. */
15 
16 const int elem_part[] = {0, 1, 1
17  C_EXPAND(2, 2, 2, 2, 2)
18  N_EXPAND(3, 3, 3, 3, 3)
19  O_EXPAND(4, 4, 4, 4, 4)
20  Ne_EXPAND(5, 5, 5, 5, 5)
21  S_EXPAND(6, 6, 6, 6, 6)
22  Fe_EXPAND(7, 7, 7) };
23 
24 const double rad_rec_z[] = { 1., 1., 2.
25  C_EXPAND(1., 2., 3., 4., 5.)
26  N_EXPAND(1., 2., 3., 4., 5.)
27  O_EXPAND(1., 2., 3., 4., 5.)
28  Ne_EXPAND(1., 2., 3., 4., 5.)
29  S_EXPAND(1., 2., 3., 4., 5.)
30  Fe_EXPAND(1., 2., 3.) };
31 const double coll_ion_dE[] = { 13.6, 24.6, 54.4
32  C_EXPAND(11.3, 24.4, 47.9, 64.5, 392.1)
33  N_EXPAND(14.5, 29.6, 47.5, 77.5, 97.9)
34  O_EXPAND(13.6, 35.1, 54.9, 77.4, 113.9)
35  Ne_EXPAND(21.6, 41.0, 63.5, 97.1, 126.2)
36  S_EXPAND(10.4, 23.3, 34.8, 47.3, 72.6)
37  Fe_EXPAND(7.87, 16.1879, 30.652) };
38 COOL_COEFF CoolCoeffs;
39 /* ********************************************************************* */
40 void Radiat (double *v, double *rhs)
41 /*!
42  * Cooling for optically thin plasma up to about 200,000 K
43  * Plasma composition: H, HeI-II, CI-V, NI-V, OI-V, NeI-V, SI-V
44  * Assumed abundances in elem_ab
45  * Uses S : Array = Variables vector x line points
46  * rhs : output for the system of ODE
47  * ibeg, iend : begin and end points of the current line
48  *
49  *
50  *********************************************************************** */
51 {
52  int nv, j, k, cooling_nT, cooling_nNe, jj;
53  int ti1, ti2, ne1, ne2, nrt;
54  double mu, sT, scrh, tmpT, tmpNe, tt1, tt2, nn1, nn2, tf1, tf2, nf1, nf2;
55  double N, n_el, rlosst, em, cf1, cf2;
56  double T, *X, *RS;
57  static double ***tab;
58  static double N_rho;
59  static double E_cost, Unit_Time; /* -- for dimensionalization purposes -- */
60  double em2, em3;
61 double prs;
62 
63 /* --------------------------------------------------------------------
64  Load tables from disk at first call.
65  Tables are 2-D with T and Ne being the coordinates.
66  -------------------------------------------------------------------- */
67 
68  if (tab == NULL) {
69 
70  E_cost = UNIT_LENGTH/UNIT_DENSITY/
71  (UNIT_VELOCITY*UNIT_VELOCITY*UNIT_VELOCITY);
72  Unit_Time = UNIT_LENGTH/UNIT_VELOCITY;
73 
74  /* ---------------------------------------
75  Compute cooling function tables
76  --------------------------------------- */
77 
78  for (nv = 0; nv < NIONS; nv++) CoolCoeffs.dLIR_dX[nv] = 0.0;
79 
80  n_el = C_NeMIN;
81  ne1 = 0;
82  while (n_el < C_NeMAX) {
83  T = C_TMIN;
84  ne2 = 0;
85  while (T < C_TMAX) {
86  tmpT = T;
87  T = T*exp(C_TSTEP); /* should be *exp(0.02) */
88  ne2 = ne2 + 1;
89  }
90  tmpNe = n_el;
91  n_el = n_el*exp(C_NeSTEP); /* should be *exp(0.06) */
92  ne1 = ne1 + 1;
93  }
94 
95  tab = ARRAY_3D(NIONS, ne1, ne2, double);
96  Create_Losses_Tables(tab, &cooling_nT, &cooling_nNe);
97  N_rho = find_N_rho();
98  } /* end load tables */
99 
100 /* ---------------------------------------
101  Force species to lie between 0 and 1
102  --------------------------------------- */
103 
104  for (nv = NFLX; nv < (NFLX + NIONS); nv++){
105  v[nv] = MIN(v[nv], 1.0);
106  v[nv] = MAX(v[nv], 0.0);
107  }
108 
109  prs = v[RHOE]*(g_gamma-1.0);
110  if (prs < 0.0) {
111  prs = g_smallPressure;
112  v[RHOE] = prs/(g_gamma - 1.0);
113  }
114 
115  mu = MeanMolecularWeight(v);
116  T = prs/v[RHO]*KELVIN*mu;
117 
118  if (mu < 0.0){
119  print ("! Radiat: negative mu\n");
120  QUIT_PLUTO(1);
121  }
122 
123  /* ---------------------------------
124  offset pointers to work more
125  efficiently.
126  --------------------------------- */
127 
128  X = v + NFLX;
129  RS = rhs + NFLX;
130 
131  sT = sqrt(T);
132  N = v[RHO]*N_rho; /* -- Total number density -- */
133 
134 /* -----------------------------------
135  compute electron number density
136  ----------------------------------- */
137 
138  CoolCoeffs.dnel_dX[0] = -N*elem_ab[el_H];
139  n_el = N*(1.0 - v[X_HI])*elem_ab[el_H]; /* -- contribution from ionized H -- */
140 
141  for (nv = 1; nv < NIONS; nv++) {
142  CoolCoeffs.dnel_dX[nv] = N*(rad_rec_z[nv] - 1.0)*elem_ab[elem_part[nv]];
143  n_el += X[nv]*CoolCoeffs.dnel_dX[nv];
144  }
145 
146  if (n_el/N < 1.e-4) n_el = N*1.e-4; /* OK ????? */
147 
148  CoolCoeffs.Ne = n_el;
149 
150  Find_Rates(T, n_el, N, v); /* find transition rates */
151 
152 /* -----------------------------------------------------------------------------
153  Compute RHS of the system of ODE.
154  Evaluate the transition rates
155  (multiply the rates given by find_rates()
156  by numerical densities)
157  ----------------------------------------------------------------------------- */
158 
159  /* ----------------------------------------------
160  Compute right hand side for all ions
161  ---------------------------------------------- */
162 
163  for (nv = 1; nv <= NIONS - 2; nv++) {
164  RS[nv] = X[nv - 1]*CoolCoeffs.Lrate[nv]
165  - X[nv] *CoolCoeffs.Crate[nv]
166  + X[nv + 1]*CoolCoeffs.Rrate[nv];
167  }
168  RS[0] = (1.0 - X[0])*CoolCoeffs.Rrate[0] - X[0]*CoolCoeffs.Crate[0]; /* separate for H */
169 
170  nv = NIONS - 1;
171  RS[nv] = X[nv - 1]*CoolCoeffs.Lrate[nv] - X[nv]*CoolCoeffs.Crate[nv];
172 
173  for (nv = 0; nv < NIONS; nv++) RS[nv] *= Unit_Time;
174 
175 /* ---------------------------------------------------------------
176  check sums
177  --------------------------------------------------------------- */
178 /*
179 {
180 double sum;
181 
182 sum=fabs(RS[1] + RS[2]);
183 if (sum > 1.e-17){
184  print("> Sum(RHS) != 0 for He!!\n");
185  exit(1);
186 }
187 
188 sum = 0.0;
189 for (nv = 0; nv < C_IONS; nv++) sum += RS[CI-NFLX+nv];
190 if (fabs(sum) > 1.e-10){
191  print("> Sum(RHS) != 0 for C !! (%12.6e)\n", sum);
192  exit(1);
193 }
194 }
195 
196 nv = 8;
197 sum=fabs(RS[nv] + RS[nv+1] + RS[nv+2] + RS[nv+3] + RS[nv+4]);
198 if ( sum > 1.e-10){
199  print("> Sum(RHS) != 0 for N!!\n");
200  exit(1);
201 }
202 nv = 13;
203 sum=fabs(RS[nv] + RS[nv+1] + RS[nv+2] + RS[nv+3] + RS[nv+4]);
204 if ( sum > 1.e-10){
205  print("> Sum(RHS) != 0 for O!!\n");
206  exit(1);
207 }
208 nv = 18;
209 sum = fabs(RS[nv] + RS[nv+1] + RS[nv+2] + RS[nv+3] + RS[nv+4]);
210 if ( sum > 1.e-10){
211  print("> Sum(RHS) != 0 for Ne!!\n");
212  exit(1);
213 }
214 nv = 23;
215 fabs(RS[nv] + RS[nv+1] + RS[nv+2] + RS[nv+3] + RS[nv+4]);
216 if ( sum > 1.e-10){
217  print("> Sum(RHS) != 0 for S!!\n");
218  exit(1);
219 }
220 }
221 */
222 
223 /*
224 print("%12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n", v[HI],
225  RS[0],RS[1], RS[2], RS[3], RS[4]);
226 }
227 exit(1);
228 */
229 
230 /* ********************************************************************************************************* */
231 
232 
233  /* -------------------------------------------------------------
234  Find the indexes in the tables needed for interpolation
235 
236  Constrain both T and n_el to lie in the range provided
237  by the tables.
238  This does not affect the original vector of primitive
239  quantities, but it is more computationally efficient.
240  In other words, when T > Tmax, the cooling function
241  will saturate at Tmax.
242  ------------------------------------------------------------- */
243 
244  tmpT = T; /* keep the double value for after the interpolation */
245  tmpNe = n_el;
246 
247  if (tmpT <= C_TMIN) tmpT = C_TMIN + C_TSTEP*0.5;
248  if (tmpT >= C_TMAX) tmpT = C_TMAX - C_TSTEP*0.5;
249 
250  if (n_el <= C_NeMIN) tmpNe = C_NeMIN + C_NeSTEP*0.5;
251  if (n_el >= C_NeMAX) tmpNe = C_NeMAX - C_NeSTEP*0.5;
252 
253  scrh = log(tmpT/C_TMIN)/C_TSTEP;
254  ti1 = floor(scrh);
255  ti2 = ceil (scrh);
256 
257  scrh = log(tmpNe/C_NeMIN)/C_NeSTEP;
258  ne1 = floor(scrh);
259  ne2 = ceil (scrh);
260 
261  tt1 = C_TMIN*exp(ti1*C_TSTEP);
262  tt2 = C_TMIN*exp(ti2*C_TSTEP);
263 
264  nn1 = C_NeMIN*exp(ne1*C_NeSTEP);
265  nn2 = C_NeMIN*exp(ne2*C_NeSTEP);
266 
267  tf1 = (tt2 - tmpT)/(tt2 - tt1);
268  tf2 = (tmpT - tt1)/(tt2 - tt1);
269  nf1 = (nn2 - tmpNe)/(nn2 - nn1);
270  nf2 = (tmpNe - nn1)/(nn2 - nn1);
271 
272 /* ------------------------------------------------------------
273  get \Delta e[k] coefficients
274  ------------------------------------------------------------ */
275 
276  for (nv = 0; nv < NIONS; nv++) {
277  cf1 = tf1*tab[nv][ne1][ti1] + tf2*tab[nv][ne1][ti2];
278  cf2 = tf1*tab[nv][ne2][ti1] + tf2*tab[nv][ne2][ti2];
279  CoolCoeffs.de[nv] = nf1*cf1 + nf2*cf2;
280  }
281 
282 /* -- Save the slopes of the radiative losses coefficients from the table -- */
283 
284  scrh = 0.5/(nn2 - nn1);
285  for (nv = 0; nv < NIONS; nv++) {
286  CoolCoeffs.de_dne[nv] = ( tab[nv][ne2][ti1] + tab[nv][ne2][ti2]
287  - tab[nv][ne1][ti1] - tab[nv][ne1][ti2])*scrh;
288  }
289 
290 /* ---------------------------------------------------------------
291  For HeII we impose a temperature cut-off above 9.e4 K
292  since HeII --> HeIII
293  --------------------------------------------------------------- */
294 
295  if (T >= 9.e4) {
296  scrh = exp(-(T - 9.e+4)/4.e+4);
297  CoolCoeffs.de[2] *= scrh;
298  CoolCoeffs.de_dne[2] *= scrh;
299  }
300 
301  /* ... and then continue with the radiative losses computation */
302 
303  em = 0.0;
304  for (nv = 0; nv < NIONS; nv++) {
305  em += X[nv]*CoolCoeffs.de[nv]*elem_ab[elem_part[nv]];
306  }
307 
308  if (em < 0.0) print("> 2 negative em => possitive losses???? em = %12.6e\n",em);
309 
310 /* -------------------------------------------------------------------
311  Now finally add the ionization / recombination losses ...
312  ------------------------------------------------------------------- */
313 
314  /* -----------------------------------
315  thermal energy lost by ionization
316  ----------------------------------- */
317 
318  CoolCoeffs.dLIR_dX[0] = 1.08e-8*sT*elem_ab[el_H]/13.6*exp(-157890.0/T)*CONST_eV;
319  em += CoolCoeffs.dLIR_dX[0]*v[X_HI];
320 
321  /*
322  if ( (1.08e-8*sT*v[HI]*elem_ab[el_H]/13.6*exp(-157890.0/T) *CONST_eV) / (v[HI] *CoolCoeffs.Crate[0]*elem_ab[el_H]*coll_ion_dE[0]*Unit_Time*CONST_eV) > 2.0) {
323  print("too different em's: %12.6e - %12.6e\n",
324  (1.08e-8*sT*v[HI]*elem_ab[el_H]/13.6*exp(-157890.0/T) *CONST_eV), (v[HI] *CoolCoeffs.Crate[0]*elem_ab[el_H]*coll_ion_dE[0]*Unit_Time*CONST_eV));
325 
326  }
327  if ( (1.08e-8*sT*v[HI]*elem_ab[el_H]/13.6*exp(-157890.0/T) *CONST_eV) / (v[HI] *CoolCoeffs.Crate[0]*elem_ab[el_H]*coll_ion_dE[0]*Unit_Time*CONST_eV) < 0.5) {
328  print("too different em's: %12.6e - %12.6e\n",
329  (1.08e-8*sT*v[HI]*elem_ab[el_H]/13.6*exp(-157890.0/T) *CONST_eV), (v[HI] *CoolCoeffs.Crate[0]*elem_ab[el_H]*coll_ion_dE[0]*Unit_Time*CONST_eV));
330 
331  }
332 */
333 /*
334  em2 = v[HI]*CoolCoeffs.Crate[0]*elem_ab[el_H]*coll_ion_dE[0]*CONST_eV*N;
335  for (nv=2; nv<NIONS; nv++) {
336  if (CoolCoeffs.Lrate[nv]>0) em2 += X[nv-1]*CoolCoeffs.Lrate[nv]*elem_ab[elem_part[nv]]*coll_ion_dE[nv]*CONST_eV*N;
337  }
338 */
339  /* --------------------------------------
340  thermal energy lost by recombination
341  (2/3 kT per recombination).
342  -------------------------------------- */
343 
344  scrh = 2.6e-11/sT*elem_ab[el_H]*T/11590.0*0.67*CONST_eV;
345  em += scrh*(1.0 - v[X_HI]);
346  CoolCoeffs.dLIR_dX[0] -= scrh;
347 
348 /*
349  if (T>4000) print("ratesHI: old %12.6e new %12.6e at T = %12.6e\n",2.6e-11/sT*(1.0 - v[HI])*elem_ab[el_H]*T/11590.0*0.67*CONST_eV,
350  (1.0 - v[HI])*CoolCoeffs.Rrate[0]*elem_ab[el_H]*0.67*CONST_kB*T,T);
351  */
352 /*
353  em3 = (1.0 - v[HI])*CoolCoeffs.Rrate[0]*elem_ab[el_H]*0.67*CONST_kB*T*N;
354  for (nv=1; nv<NIONS-1; nv++) {
355  if (CoolCoeffs.Rrate[nv]>0) em3 += X[nv+1]*CoolCoeffs.Rrate[nv]*elem_ab[elem_part[nv]]*0.67*CONST_kB*T*N;
356  }
357 */
358 
359  /* --------------------------------------------------------
360  Add cooling contribution from bremsstrahlung
361  -------------------------------------------------------- */
362 
363  scrh = 1.42e-27*sT*elem_ab[el_H];
364  CoolCoeffs.dLIR_dX[2] = 1.42e-27*sT*elem_ab[el_He];
365  em += scrh*(1.0 - v[X_HI]) + CoolCoeffs.dLIR_dX[2]*v[X_HeII];
366 
367  CoolCoeffs.dLIR_dX[0] -= scrh;
368 
369  /* --------------------------------------------------------------------
370  Add cooling contribution from FeII 1.6 & 25 micron
371  and Mg II 2800 and Si II 35 micron
372  emission lines. This is TEMPORARY until the adding
373  of FeI and FeII as ions to NEq !
374  -------------------------------------------------------------------- */
375 
376  em += 1.6e-12*8.63e-6*0.3*0.0495*exp(-575.0/T)/sT * 580.0*sT/(n_el + 580.0*sT) *0.00004;
377  em += 1.6e-12*8.63e-6*0.39*0.775*exp(-8980.0/T)/sT * 1130.0*sT/(n_el + 1130.0*sT) *0.00004;
378 
379  em += 1.6e-12*8.63e-6*8.0*4.43*exp(-5.13e4/T)/sT*1.e10*sT/(n_el + 1.e10*sT)*0.00002*0.7; /* Mg II 2800: Mendoza */
380  em += 1.6e-12*8.63e-6*2.85*0.0354*exp(-410.0/T)/sT*16.8*sT/(n_el + 16.8*sT)*0.000033*0.7; /* Si II 35 micron: Dufton&Kingston, MNRAS 248 */
381 
382  /* ---------------------------------------------------
383  rlosst is the energy loss in units of erg/cm^3/s;
384  it must be multiplied by cost_E in order to match
385  non-dimensional units.
386  Source term for the neutral fraction scales with
387  Unit_Time
388 
389  Pressure source term equals zero when
390  T < g_minCoolingTemp
391  --------------------------------------------------- */
392 
393  rlosst = em*n_el*N;
394 /*
395  rlosst += em3 + em2;
396 */
397 
398 /*
399  if (T > g_minCoolingTemp) rhs[PRS] = -E_cost*rlosst*(g_gamma - 1.0);
400  else rhs[PRS] = 0.0;
401 */
402 
403  rhs[RHOE] = -E_cost*rlosst;
404  rhs[RHOE] *= 1.0/(1.0 + exp( -(T - g_minCoolingTemp)/100.)); /* -- cut cooling function --*/
405 
406  if (rhs[RHOE] > 0.0) {
407  print ("! Error: positive radiative losses! %12.6e %12.6e %12.6e \n", em, n_el, N);
408  QUIT_PLUTO(1);
409  }
410 }
411 
412 /* ********************************************************************* */
413 void Find_Rates(double T, double Ne, double N, double *v)
414 /*
415  *
416  * PURPOSE
417  *
418  * Transition rates computation routine
419  *
420  * Output: CoolCoeffs.Crate, CoolCoeffs.Rrate, CoolCoeffs.Lrate
421  *
422  * dN = ( CoolCoeffs.Lrate * n_(ion-1)
423  * + CoolCoeffs.Rrate * n_(ion+1) - CoolCoeffs.Crate * n_ion ) * dt
424  *
425  *
426  *********************************************************************** */
427 {
428  double dn, lam, t4, tmprec, ft1, ft2, tmpT, scrh;
429  int ti1, ti2, cindex, i, ions, nv, tindex, j, k;
430  static double ***ion_data, **intData;
431 
432  if (ion_data == NULL) { /* -- compute ionization rates tables -- */
433  ion_data = ARRAY_3D(I_g_stepNumber, 7, NIONS, double);
434  intData = ARRAY_2D(7, NIONS, double);
435  Create_Ion_Coeff_Tables(ion_data);
436  }
437 
438  for (nv = 0; nv < NIONS; nv++ ) {
439  CoolCoeffs.Rrate[nv] = 0.0;
440  CoolCoeffs.Lrate[nv] = 0.0;
441  CoolCoeffs.Crate[nv] = 0.0;
442  CoolCoeffs.Ra[nv] = 0.0;
443  CoolCoeffs.La[nv] = 0.0;
444  CoolCoeffs.Ca[nv] = 0.0;
445  CoolCoeffs.Rb[nv] = 0.0;
446  CoolCoeffs.Lb[nv] = 0.0;
447  CoolCoeffs.Cb[nv] = 0.0;
448  CoolCoeffs.Rc[nv] = 0.0;
449  CoolCoeffs.Lc[nv] = 0.0;
450  CoolCoeffs.Cc[nv] = 0.0;
451  }
452 
453  tmpT = T;
454  if (tmpT > I_TEND) tmpT = I_TEND - I_TSTEP*0.5;
455  else if (tmpT < I_TBEG) tmpT = I_TBEG + I_TSTEP*0.5;
456 
457  tindex = floor( (tmpT - I_TBEG)/I_TSTEP);
458 
459  ft1 = ( I_TBEG + I_TSTEP*((double)tindex + 1.0) - tmpT ) / I_TSTEP;
460  ft2 = ( tmpT - I_TBEG - I_TSTEP*((double)tindex) ) / I_TSTEP;
461 
462  for (k = 0; k < 7; k++) {
463  for (ions = 0; ions < NIONS; ions++) {
464  intData[k][ions] = ion_data[tindex][k][ions]*ft1 + ion_data[tindex + 1][k][ions]*ft2;
465  }}
466 
467  /* -------------------------------------------------------------
468  Computing rates for H
469  ------------------------------------------------------------- */
470 
471  lam = 157890. / T;
472  CoolCoeffs.Crate[0] += Ne*intData[0][el_H]; /* collisional ionization */
473  CoolCoeffs.Rrate[0] += Ne*intData[1][el_H]; /* total recombination - NIFS-DATA-54 - Kato & Asano 1999 , from Aldrovandi & Pequignot 1973 */
474 
475  CoolCoeffs.fCH = intData[0][el_H];
476  CoolCoeffs.fRH = intData[1][el_H];
477 
478  /* contribution from charge - exchange (recombination and ionization of all ions) : */
479 /*
480  for (ions=3; ions<28; ions++) {
481 
482  if (T/10000.<=chtrH_ion_upT[ions]) CoolCoeffs.Crate[0] += X[ions+1][ii] * N * elem_ab[elem_part[ions+1]] * 1.e-9 * chtrH_rec_a[ions] * pow( (T/10000.), chtrH_rec_b[ions] ) * ( 1. + chtrH_rec_c[ions] * exp ( chtrH_rec_d[ions]*(T/10000.) ));
483  if (T<=chtrH_rec_upT[ions]) CoolCoeffs.Rrate[0] += 1.e-9 * chtrH_ion_a[ions] * pow( (T/10000.), chtrH_ion_b[ions] ) * ( 1. + chtrH_ion_c[ions] * exp ( chtrH_ion_d[ions]*(T/10000.) ));
484 
485  tmprec = v[ions] * elem_ab[elem_part[ions]] / elem_ab[el_H] / fabs(v[HI]-1.e-13) * N * ion_data[4][ions];
486  CoolCoeffs.Crate[0] += tmprec;
487  CoolCoeffs.Ca[0] += ion_data[4][ions];
488  }
489 */
490 
491 /* -------------------------------------------------------------------
492  Do the other ions now
493  ------------------------------------------------------------------- */
494 
495  for (ions = 2; ions < NIONS - 1; ions++) {
496 
497  /* ------------------------------------------
498  collisional ionization - Voronov 1997
499  ------------------------------------------ */
500 
501  tmprec = Ne*intData[0][ions];
502  CoolCoeffs.Crate[ions] += tmprec;
503  CoolCoeffs.Lrate[ions+1] += tmprec;
504  CoolCoeffs.Ca[ions] += intData[0][ions];
505  CoolCoeffs.La[ions+1] += intData[0][ions];
506 
507  /* -------------------------------------------------------
508  radiative recombination - Pequignot & al, 1991 A&A
509  ------------------------------------------------------- */
510 
511  tmprec = Ne*intData[1][ions-1];
512  CoolCoeffs.Crate[ions] += tmprec;
513  CoolCoeffs.Rrate[ions-1] += tmprec;
514  CoolCoeffs.Ca[ions] += intData[1][ions-1];
515  CoolCoeffs.Ra[ions-1] += intData[1][ions-1];
516 
517  /* -----------------------------------------------------------
518  dielectronic recombination - Nussbaumer & Storey, 1983
519  ----------------------------------------------------------- */
520 
521  tmprec = Ne*intData[2][ions-1];
522  CoolCoeffs.Crate[ions] += tmprec;
523  CoolCoeffs.Rrate[ions-1] += tmprec;
524  CoolCoeffs.Ca[ions] += intData[2][ions-1];
525  CoolCoeffs.Ra[ions-1] += intData[2][ions-1];
526 
527  /* --------------------------------------------------------------------------
528  charge-transfer with H+ - ionization - Kingdon & Ferland 1996, ApJSS
529  -------------------------------------------------------------------------- */
530 
531  scrh = N*elem_ab[el_H]*intData[3][ions];
532  tmprec = (1.0 - v[X_HI])*scrh;
533  CoolCoeffs.Crate[ions] += tmprec;
534  CoolCoeffs.Lrate[ions+1] += tmprec;
535  CoolCoeffs.Cb[ions] -= scrh;
536  CoolCoeffs.Lb[ions+1] -= scrh;
537 
538  /* ---------------------------------------------------------------------------
539  charge-transfer with H - recombination - Kingdon & Ferland 1996, ApJSS
540  --------------------------------------------------------------------------- */
541 
542  scrh = N*elem_ab[el_H]*intData[4][ions-1];
543  tmprec = v[X_HI]*scrh;
544  CoolCoeffs.Crate[ions] += tmprec;
545  CoolCoeffs.Rrate[ions-1] += tmprec;
546  CoolCoeffs.Cb[ions] += scrh;
547  CoolCoeffs.Rb[ions-1] += scrh;
548  }
549 
550  /* -- calculation for last element (ions = NIONS - 1) -- */
551 
552  tmprec = Ne*intData[1][ions-1];
553  CoolCoeffs.Crate[ions] += tmprec;
554  CoolCoeffs.Rrate[ions-1] += tmprec;
555  CoolCoeffs.Ca[ions] += intData[1][ions-1];
556  CoolCoeffs.Ra[ions-1] += intData[1][ions-1];
557 
558  tmprec = Ne*intData[2][ions-1];
559  CoolCoeffs.Crate[ions] += tmprec;
560  CoolCoeffs.Rrate[ions-1] += tmprec;
561  CoolCoeffs.Ca[ions] += intData[2][ions-1];
562  CoolCoeffs.Ra[ions-1] += intData[2][ions-1];
563 
564  scrh = N*elem_ab[el_H]*intData[4][ions-1];
565  tmprec = v[X_HI]*scrh;
566  CoolCoeffs.Crate[ions] += tmprec;
567  CoolCoeffs.Rrate[ions-1] += tmprec;
568  CoolCoeffs.Cb[ions] += scrh;
569  CoolCoeffs.Rb[ions-1] += scrh;
570 
571  /* -- calculation for HeI -- */
572 
573  ions = 1;
574  tmprec = Ne*intData[0][ions];
575  CoolCoeffs.Crate[ions] += tmprec;
576  CoolCoeffs.Lrate[ions+1] += tmprec;
577  CoolCoeffs.Ca[ions] += intData[0][ions];
578  CoolCoeffs.La[ions+1] += intData[0][ions];
579 
580  scrh = N*elem_ab[el_H]*intData[3][ions];
581  tmprec = (1.0 - v[X_HI])*scrh;
582  CoolCoeffs.Crate[ions] += tmprec;
583  CoolCoeffs.Lrate[ions+1] += tmprec;
584  CoolCoeffs.Cb[ions] -= scrh;
585  CoolCoeffs.Lb[ions+1] -= scrh;
586 
587 /* --------------------------------------------------------------------------
588  charge-transfer with He - recombination
589  ORNL Charge Transfer Database
590  http://www-cfadc.phy.ornl.gov/astro/ps/data/cx/helium/rates/fits.data
591  -------------------------------------------------------------------------- */
592 
593  for (ions = 4; ions < NIONS; ions++) { /* 4 ??????????????? */
594  scrh = N*elem_ab[el_He]*intData[5][ions-1];
595  tmprec = v[X_HeI]*scrh;
596  CoolCoeffs.Crate[ions] += tmprec;
597  CoolCoeffs.Rrate[ions-1] += tmprec;
598 
599  CoolCoeffs.Cc[ions] += scrh;
600  CoolCoeffs.Rc[ions-1] += scrh;
601  }
602 
603 /* -------------------------------------------------------
604  And now, for ions for which we only have the total
605  electron-ion recombination coefficient:
606  ( NIFS-DATA-54 - Kato & Asano 1999 )
607  ------------------------------------------------------- */
608 
609  for (ions = 2; ions < NIONS; ions++) {
610  tmprec = Ne*intData[6][ions-1];
611  CoolCoeffs.Crate[ions] += tmprec;
612  CoolCoeffs.Rrate[ions-1] += tmprec;
613 
614  /* -----------------------------------------------
615  Now save the data for the partial derivatives
616  computation
617  ----------------------------------------------- */
618 
619  CoolCoeffs.Ca[ions] += intData[6][ions-1];
620  CoolCoeffs.Ra[ions-1] += intData[6][ions-1];
621  }
622 
623  /* He/He+ contribution from charge-exchange processes: */
624 /*
625  for (ions=3; ions<28; ions++) {
626  if (T<5000.) CoolCoeffs.Rrate[2] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a1[ions] * pow( (T/10000.), chtrHe_rec_b1[ions] ) * ( 1. + chtrHe_rec_c1[ions] * exp ( chtrHe_rec_d1[ions]*(T/10000.) ));
627  if ( (T>=5000.) && (T<10000.) ) CoolCoeffs.Rrate[2] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a2[ions] * pow( (T/10000.), chtrHe_rec_b2[ions] ) * ( 1. + chtrHe_rec_c2[ions] * exp ( chtrHe_rec_d2[ions]*(T/10000.) ));
628  if ( (T>=10000.) && (T<50000.) ) CoolCoeffs.Rrate[2] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a3[ions] * pow( (T/10000.), chtrHe_rec_b3[ions] ) * ( 1. + chtrHe_rec_c3[ions] * exp ( chtrHe_rec_d3[ions]*(T/10000.) ));
629  if ( (T>=50000.) && (T<100000.) ) CoolCoeffs.Rrate[2] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a4[ions] * pow( (T/10000.), chtrHe_rec_b4[ions] ) * ( 1. + chtrHe_rec_c4[ions] * exp ( chtrHe_rec_d4[ions]*(T/10000.) ));
630  if (T>=100000.) CoolCoeffs.Rrate[2] += N * v[ions] *elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a5[ions] * pow( (T/10000.), chtrHe_rec_b5[ions] ) * ( 1. + chtrHe_rec_c5[ions] * exp ( chtrHe_rec_d5[ions]*(T/10000.) ));
631 
632  if (T<5000.) CoolCoeffs.Crate[1] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a1[ions] * pow( (T/10000.), chtrHe_rec_b1[ions] ) * ( 1. + chtrHe_rec_c1[ions] * exp ( chtrHe_rec_d1[ions]*(T/10000.) ));
633  if ( (T>=5000.) && (T<10000.) ) CoolCoeffs.Crate[1] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a2[ions] * pow( (T/10000.), chtrHe_rec_b2[ions] ) * ( 1. + chtrHe_rec_c2[ions] * exp ( chtrHe_rec_d2[ions]*(T/10000.) ));
634  if ( (T>=10000.) && (T<50000.) ) CoolCoeffs.Crate[1] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a3[ions] * pow( (T/10000.), chtrHe_rec_b3[ions] ) * ( 1. + chtrHe_rec_c3[ions] * exp ( chtrHe_rec_d3[ions]*(T/10000.) ));
635  if ( (T>=50000.) && (T<100000.) ) CoolCoeffs.Crate[1] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a4[ions] * pow( (T/10000.), chtrHe_rec_b4[ions] ) * ( 1. + chtrHe_rec_c4[ions] * exp ( chtrHe_rec_d4[ions]*(T/10000.) ));
636  if (T>=100000.) CoolCoeffs.Crate[1] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a5[ions] * pow( (T/10000.), chtrHe_rec_b5[ions] ) * ( 1. + chtrHe_rec_c5[ions] * exp ( chtrHe_rec_d5[ions]*(T/10000.) ));
637  }
638 */
639 
640 /* -------------------------------------------------------------
641  C H E C K
642  ------------------------------------------------------------- */
643 /*
644 for (ions = 0; ions < NIONS; ions++) {
645  CoolCoeffs.Lrate[ions] *= UNIT_LENGTH/UNIT_VELOCITY;
646  CoolCoeffs.Crate[ions] *= UNIT_LENGTH/UNIT_VELOCITY;
647  CoolCoeffs.Rrate[ions] *= UNIT_LENGTH/UNIT_VELOCITY;
648 }
649 
650 print ("H: %12.6e %12.6e %12.6e\n",
651  CoolCoeffs.Lrate[HI-NFLX], CoolCoeffs.Crate[HI-NFLX],CoolCoeffs.Rrate[HI-NFLX]);
652 
653 print ("--------------------------------------------------------------- \n");
654 for (nv = HeI; nv <= HeII; nv++){
655 print ("He(%d): %12.6e %12.6e %12.6e\n",nv-HeI+1,
656  CoolCoeffs.Lrate[nv-NFLX], CoolCoeffs.Crate[nv-NFLX],CoolCoeffs.Rrate[nv-NFLX]);
657 }
658 
659 print ("--------------------------------------------------------------- \n");
660 for (nv = CI; nv <= CV; nv++){
661 print ("C(%d): %12.6e %12.6e %12.6e\n",nv-CI+1,
662  CoolCoeffs.Lrate[nv-NFLX], CoolCoeffs.Crate[nv-NFLX],CoolCoeffs.Rrate[nv-NFLX]);
663 }
664 print ("--------------------------------------------------------------- \n");
665 
666 for (nv = NI; nv <= NV; nv++){
667 print ("N(%d): %12.6e %12.6e %12.6e\n",nv-NI+1,
668  CoolCoeffs.Lrate[nv-NFLX], CoolCoeffs.Crate[nv-NFLX],CoolCoeffs.Rrate[nv-NFLX]);
669 }
670 print ("--------------------------------------------------------------- \n");
671 
672 for (nv = OI; nv <= OV; nv++){
673 print ("O(%d): %12.6e %12.6e %12.6e\n",nv-OI+1,
674  CoolCoeffs.Lrate[nv-NFLX], CoolCoeffs.Crate[nv-NFLX],CoolCoeffs.Rrate[nv-NFLX]);
675 }
676 print ("--------------------------------------------------------------- \n");
677 
678 for (nv = NeI; nv <= NeV; nv++){
679 print ("Ne(%d): %12.6e %12.6e %12.6e\n",nv-NeI+1,
680  CoolCoeffs.Lrate[nv-NFLX], CoolCoeffs.Crate[nv-NFLX],CoolCoeffs.Rrate[nv-NFLX]);
681 }
682 print ("--------------------------------------------------------------- \n");
683 for (nv = SI; nv <= SV; nv++){
684 print ("Se(%d): %12.6e %12.6e %12.6e\n",nv-SI+1,
685  CoolCoeffs.Lrate[nv-NFLX], CoolCoeffs.Crate[nv-NFLX],CoolCoeffs.Rrate[nv-NFLX]);
686 }
687 
688 
689 exit(1);
690 */
691 }
692 
693 /* ****************************************************************************** */
694 double find_N_rho ()
695 /*
696  *
697  * PURPOSE
698  *
699  * Find the number density of atoms / ions (divided by the density so that is
700  * a constant depending on composition only), knowing the ions ratios and density.
701  * We compute a /mu taking into account only the heavy particles (nuclei).
702  *
703  *
704  * LAST MODIFIED: July 18, 2006 by Ovidiu Tesileanu
705  *
706  ******************************************************************************* */
707 {
708  int i, j;
709  double mu, mu1, mu2;
710 
711  mu1 = mu2 = 0.0;
712  for (i = 0; i < 7; i++) {
713  mu1 += elem_ab[i]*elem_mass[i]; /* Numerator part of mu */
714  mu2 += elem_ab[i]; /* Denominator part of mu */
715  }
716  mu = mu1 / mu2;
717  return (UNIT_DENSITY/mu*CONST_NA); /* This is N/rho, with N the total number density of atoms and ions */
718 }
719 
720 /* ********************************************************************* */
721 double H_MassFrac (void)
722 /*!
723  * Compute the mass fraction X of Hydrogen as function of the
724  * composition of the gas.
725  *
726  * f_H A_H
727  * X = ----------------
728  * \sum_k f_k A_k
729  *
730  * where
731  *
732  * f_k : is the fractional abundance (by number),
733  * f_k = N_k/N_tot
734  * of atomic species (no matter ionization degree).
735  *
736  * A_K : is the atomic weight
737  *
738  *
739  * where N_{tot} is the total number density of particles
740  *
741  * ARGUMENTS
742  *
743  * none
744  *
745  *********************************************************************** */
746 {
747  double mmw2;
748  int i, j;
749 
750  mmw2 = 0.0;
751  for (i = 0; i < (Fe_IONS==0?7:8); i++) {
752  mmw2 += elem_mass[i]*elem_ab[i]; /* Denominator part of X */
753  }
754  return ( (elem_mass[0]*elem_ab[0]) / mmw2);
755 }
756 
757 
758 #ifdef CH_SPACEDIM
759 /* ********************************************************************* */
760 void NormalizeIons (double *u)
761 /*!
762  * This function re-normalize the set of conservative variables u in
763  * such a way that the sum of all ion fractions is equal to 1.
764  * Since u is a vector of conserved quantities, the normalization
765  * condition is:
766  *
767  * \sum_j (\rho*X_j) = \rho
768  *
769  * where, for instance, X_j = {OI, OII, OIII, OIV, OV}.
770  * In order to fulfill the previous equation, each conservative
771  * variable is redefined according to
772  *
773  * (\rho*X_j)
774  * (\rho*X_j) --> \rho * -------------------
775  * \sum_j (\rho*X_j)
776  *
777  * This function is necessary only when used with Chombo AMR since
778  * the coarse-to-fine interpolation fails to conserve the correct
779  * normalization. Re-fluxing and fine-to-coarse do not need this
780  * treatment since the normalization is preserved.
781  *
782  *********************************************************************** */
783 {
784  int nv;
785  double phi;
786 
787  phi = 0.0;
788  for (nv = X_HeI; nv <= X_HeII; nv++) phi += u[nv];
789  for (nv = X_HeI; nv <= X_HeII; nv++) u[nv] *= u[RHO]/phi;
790 
791  phi = 0.0;
792  for (nv = X_SI; nv < X_SI+S_IONS; nv++) phi += u[nv];
793  for (nv = X_SI; nv < X_SI+S_IONS; nv++) u[nv] *= u[RHO]/phi;
794 
795  phi = 0.0;
796  for (nv = X_OI; nv < X_OI+O_IONS; nv++) phi += u[nv];
797  for (nv = X_OI; nv < X_OI+O_IONS; nv++) u[nv] *= u[RHO]/phi;
798 
799  phi = 0.0;
800  for (nv = X_NI; nv < X_NI+N_IONS; nv++) phi += u[nv];
801  for (nv = X_NI; nv < X_NI+N_IONS; nv++) u[nv] *= u[RHO]/phi;
802 
803  phi = 0.0;
804  for (nv = X_CI; nv < X_CI+C_IONS; nv++) phi += u[nv];
805  for (nv = X_CI; nv < X_CI+C_IONS; nv++) u[nv] *= u[RHO]/phi;
806 
807  phi = 0.0;
808  for (nv = X_NeI; nv < X_NeI+Ne_IONS; nv++) phi += u[nv];
809  for (nv = X_NeI; nv < X_NeI+Ne_IONS; nv++) u[nv] *= u[RHO]/phi;
810 
811 }
812 #endif
813 
COOLING_COEFFICIENTS::dnel_dX
double dnel_dX[NIONS]
Definition: cooling_defs.h:51
MAX
#define MAX(a, b)
Definition: macros.h:101
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
COOLING_COEFFICIENTS::dLIR_dX
double dLIR_dX[NIONS]
Definition: cooling_defs.h:55
I_TBEG
#define I_TBEG
Definition: cooling_defs.h:171
S_IONS
#define S_IONS
Definition: cooling.h:17
COOLING_COEFFICIENTS::Crate
double Crate[NIONS]
Definition: cooling_defs.h:31
g_gamma
double g_gamma
Definition: globals.h:112
COOLING_COEFFICIENTS::de_dne
double de_dne[NIONS]
Definition: cooling_defs.h:47
UNIT_DENSITY
#define UNIT_DENSITY
Unit density in gr/cm^3.
Definition: pluto.h:369
COOLING_COEFFICIENTS::Cc
double Cc[NIONS]
Definition: cooling_defs.h:41
COOLING_COEFFICIENTS::Rb
double Rb[NIONS]
Definition: cooling_defs.h:43
CONST_eV
#define CONST_eV
Electron Volt in erg.
Definition: pluto.h:256
RHO
#define RHO
Definition: mod_defs.h:19
configure.scrh
tuple scrh
Definition: configure.py:200
elem_ab_uni
double elem_ab_uni[]
Definition: radiat.c:13
COOLING_COEFFICIENTS::de
double de[NIONS]
Definition: cooling_defs.h:46
COOLING_COEFFICIENTS::Lrate
double Lrate[NIONS]
Definition: cooling_defs.h:30
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
NormalizeIons
void NormalizeIons(double *)
Radiat
void Radiat(double *v, double *rhs)
Definition: radiat.c:94
g_smallPressure
double g_smallPressure
Small value for pressure fix.
Definition: globals.h:110
elem_mass
const double elem_mass[]
Definition: radiat.c:14
ARRAY_3D
#define ARRAY_3D(nx, ny, nz, type)
Definition: prototypes.h:172
KELVIN
#define KELVIN
Definition: pluto.h:401
COOLING_COEFFICIENTS::Lb
double Lb[NIONS]
Definition: cooling_defs.h:37
el_He
#define el_He
Definition: cooling_defs.h:72
S_EXPAND
#define S_EXPAND(a, b, c, d, e)
Definition: cooling.h:83
COOLING_COEFFICIENTS::fRH
double fRH
Definition: cooling_defs.h:52
N_EXPAND
#define N_EXPAND(a, b, c, d, e)
Definition: cooling.h:41
elem_ab
const double elem_ab[]
Definition: radiat.c:9
H_MassFrac
double H_MassFrac(void)
Definition: radiat.c:721
O_IONS
#define O_IONS
Definition: cooling.h:15
COOLING_COEFFICIENTS::La
double La[NIONS]
Definition: cooling_defs.h:36
MIN
#define MIN(a, b)
Definition: macros.h:104
I_TSTEP
#define I_TSTEP
Definition: cooling_defs.h:170
cooling_defs.h
rad_rec_z
const double rad_rec_z[]
Definition: radiat.c:24
UNIT_VELOCITY
#define UNIT_VELOCITY
Unit velocity in cm/sec.
Definition: pluto.h:377
elem_part
const int elem_part[]
Definition: radiat.c:16
CONST_NA
#define CONST_NA
Avogadro Contant.
Definition: pluto.h:267
C_TMAX
#define C_TMAX
Definition: cooling_defs.h:184
find_N_rho
double find_N_rho()
Definition: radiat.c:694
X
#define X
Definition: radiat.modified_by_Bhargav.c:3
X_HeI
Definition: cooling.h:110
NFLX
#define NFLX
Definition: mod_defs.h:32
I_TEND
#define I_TEND
Definition: cooling_defs.h:172
UNIT_LENGTH
#define UNIT_LENGTH
Unit Length in cm.
Definition: pluto.h:373
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
I_g_stepNumber
#define I_g_stepNumber
Definition: cooling_defs.h:169
j
int j
Definition: analysis.c:2
X_HI
Definition: cooling.h:110
k
int k
Definition: analysis.c:2
COOLING_COEFFICIENTS::fCH
double fCH
Definition: cooling_defs.h:52
COOLING_COEFFICIENTS
Definition: cooling_defs.h:28
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
coll_ion_dE
const double coll_ion_dE[]
Definition: radiat.c:31
COOLING_COEFFICIENTS::Rrate
double Rrate[NIONS]
Definition: cooling_defs.h:32
COOLING_COEFFICIENTS::Rc
double Rc[NIONS]
Definition: cooling_defs.h:44
g_minCoolingTemp
double g_minCoolingTemp
The minimum temperature (in K) below which cooling is suppressed.
Definition: globals.h:106
pluto.h
PLUTO main header file.
Ne_IONS
#define Ne_IONS
Definition: cooling.h:16
C_TMIN
#define C_TMIN
Definition: cooling_defs.h:183
i
int i
Definition: analysis.c:2
MeanMolecularWeight
double MeanMolecularWeight(double *V)
Definition: radiat.modified_by_Bhargav.c:178
COOLING_COEFFICIENTS::Ra
double Ra[NIONS]
Definition: cooling_defs.h:42
Find_Rates
void Find_Rates(double T, double Ne, double N, double *v)
Definition: radiat.c:413
O_EXPAND
#define O_EXPAND(a, b, c, d, e)
Definition: cooling.h:55
Fe_IONS
#define Fe_IONS
Definition: cooling.h:18
COOLING_COEFFICIENTS::Lc
double Lc[NIONS]
Definition: cooling_defs.h:38
C_NeMAX
#define C_NeMAX
Definition: cooling_defs.h:180
elem_ab_est
double elem_ab_est[]
Definition: radiat.c:10
elem_ab_sol
double elem_ab_sol[]
Definition: radiat.c:12
COOLING_COEFFICIENTS::Ca
double Ca[NIONS]
Definition: cooling_defs.h:39
elem_ab_lod
double elem_ab_lod[]
Definition: radiat.c:11
ARRAY_2D
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
Create_Losses_Tables
int Create_Losses_Tables(double ***, int *, int *)
Definition: make_tables.c:841
Create_Ion_Coeff_Tables
int Create_Ion_Coeff_Tables(double ***)
Definition: make_tables.c:200
COOLING_COEFFICIENTS::Cb
double Cb[NIONS]
Definition: cooling_defs.h:40
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
CoolCoeffs
COOL_COEFF CoolCoeffs
Definition: radiat.c:38
Fe_EXPAND
#define Fe_EXPAND(a, b, c)
Definition: cooling.h:91
el_H
#define el_H
Definition: cooling_defs.h:71
COOLING_COEFFICIENTS::Ne
double Ne
Definition: cooling_defs.h:34
N_IONS
#define N_IONS
Definition: cooling.h:14
C_NeMIN
#define C_NeMIN
Definition: cooling_defs.h:179