PLUTO
Macros | Functions
radiat.modified_by_Bhargav.c File Reference
#include "pluto.h"
Include dependency graph for radiat.modified_by_Bhargav.c:

Go to the source code of this file.

Macros

#define X   0.711
 
#define Y   0.2741
 
#define A_Z   30.0 /* mean atomic weight of heavy elements */
 
#define A_He   4.004 /* atomic weight of Helium */
 
#define A_H   1.008 /* atomic weight of Hydrogen */
 

Functions

void Radiat (double *v, double *rhs)
 
double MeanMolecularWeight (double *V)
 

Macro Definition Documentation

#define A_H   1.008 /* atomic weight of Hydrogen */

Definition at line 8 of file radiat.modified_by_Bhargav.c.

#define A_He   4.004 /* atomic weight of Helium */

Definition at line 7 of file radiat.modified_by_Bhargav.c.

#define A_Z   30.0 /* mean atomic weight of heavy elements */

Definition at line 6 of file radiat.modified_by_Bhargav.c.

#define X   0.711

Definition at line 3 of file radiat.modified_by_Bhargav.c.

#define Y   0.2741

Definition at line 4 of file radiat.modified_by_Bhargav.c.

Function Documentation

double MeanMolecularWeight ( double *  V)

Definition at line 178 of file radiat.modified_by_Bhargav.c.

222 {
223  int nv;
224 
225  for (nv = NFLX;nv <(NFLX+NIONS);nv++){
226  if (V[nv] < 0.0) V[nv] = 0.0;
227  if (V[nv] > 1.0) V[nv] = 1.0;
228  }
229 
230  double N_H = V[RHO]*(UNIT_DENSITY/CONST_amu)*(X/A_H); /* n(H) */
231  double N_He = V[RHO]*(UNIT_DENSITY/CONST_amu)*(Y/A_He); /* n(He) */
232  double N_Z = V[RHO]*(UNIT_DENSITY/CONST_amu)*((1.0-X-Y)/A_Z); /*/ n(metals) */
233 
234  double fracHe = N_He/N_H;
235  double fracZ = N_Z/N_H;
236 
237  double fn = V[X_HI];
238 
239  double munum = 1.0 + A_He*(fracHe) + A_Z*(fracZ);
240  double muden = 2 - fn + fracHe + fracZ + 0.5*A_Z*(fracZ);
241 
242  return munum/muden;
243 
244 
245  /*
246  return ( (A_H + frac_He*A_He + frac_Z*A_Z) /
247  (2.0 + frac_He + 2.0*frac_Z - V[X_HI]));
248  */
249 }
UNIT_DENSITY
#define UNIT_DENSITY
Unit density in gr/cm^3.
Definition: pluto.h:369
RHO
#define RHO
Definition: mod_defs.h:19
CONST_amu
#define CONST_amu
Atomic mass unit.
Definition: pluto.h:253
Y
#define Y
Definition: radiat.modified_by_Bhargav.c:4
A_H
#define A_H
Definition: radiat.modified_by_Bhargav.c:8
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
A_He
#define A_He
Definition: radiat.modified_by_Bhargav.c:7
X_HI
Definition: cooling.h:110
A_Z
#define A_Z
Definition: radiat.modified_by_Bhargav.c:6

Here is the caller graph for this function:

void Radiat ( double *  v,
double *  rhs 
)

Cooling for neutral or singly ionized gas: good up to about 35,000 K in equilibrium or shocks in neutral gas up to about 80 km/s. Assumed abundances in ab Uses t : Kelvin dene : electron density cm*-3 fneut : hydrogen neutral fraction (adimensionale) ci,cr : H ionization and recombination rate coefficients

em(1) = TOTAL EMISSIVITY : (ergs cm**3 s**-1) em(2) = Ly alpha + two photon continuum: Aggarwal MNRAS 202, 10**4.3 K em(3) = H alpha: Aggarwal, Case B em(4) = He I 584 + two photon + 623 (all n=2 excitations): Berrington &Kingston,JPB 20 em(5) = C I 9850 + 9823: Mendoza, IAU 103, 5000 K em(6) = C II, 156 micron: Mendoza, 10,000 K em(7) = C II] 2325 A: Mendoza, 15,000 K em(8) = N I 5200 A: Mendoza, 7500 K em(9) = N II 6584 + 6548 A: Mendoza em(10) = O I 63 micron: Mendoza,2500 K em(11) = O I 6300 A + 6363 A: Mendoza, 7500 K em(12) = O II 3727: Mendoza em(13) = Mg II 2800: Mendoza em(14) = Si II 35 micron: Dufton&Kingston, MNRAS 248 em(15) = S II 6717+6727: Mendoza em(16) = Fe II 25 micron: Nussbaumer&Storey em(17) = Fe II 1.6 micron em(18) = thermal energy lost by ionization em(19) = thermal energy lost by recombination (2/3 kT per recombination. The ionization energy lost is not included here.

Definition at line 11 of file radiat.modified_by_Bhargav.c.

47 {
48 
49  int ii, k;
50  double T, mu, st, rho, pr, fn;
51  double N_H, n_el, cr, ci, rlosst, em[20], src_pr;
52 
53  static double t00[18] = {0.0 , 0.0 , 1.18e5, 1.40e5, 2.46e5, 1.46e4,
54  92.1 , 6.18e4, 2.76e4, 2.19e4, 228.0 , 2.28e4,
55  3.86e4, 5.13e4, 410.0 , 2.13e4, 575.0 , 8980.0};
56 
57  static double ep[18] = {0.0 , 0.0 , 10.2 , 1.89, 21.2 , 1.26,
58  0.0079, 5.33, 2.38 , 1.89, 0.0197, 1.96,
59  3.33 , 4.43, 0.0354, 1.85, 0.0495, 0.775};
60 
61  static double critn[18] = {0.0 , 0.0 , 1.e10, 1.e10, 1.e10, 312.0,
62  0.849, 1.93e7, 124.0, 865.0, 1090.0, 3950.0,
63  177.0, 1.e10 , 16.8, 96.0, 580.0, 1130.0};
64 
65  static double om[18] = {0.0 , 0.0 , 0.90, 0.35, 0.15 , 0.067,
66  0.63, 0.52, 0.90, 0.30, 0.0055, 0.19,
67  0.33, 8.0 , 2.85, 1.75, 0.3 ,0.39};
68 
69  static double ab[18] = {0.0 , 0.0 , 1.0 , 1.0 , 0.1 , 0.0003,
70  0.0003, 0.0003 , 0.0001 , 0.0001 , 0.0006 , 0.0006,
71  0.0006, 0.00002, 0.00004, 0.00004, 0.00004, 0.00004};
72 
73  static double fn1[20] = {0.0, 0.0, 0.0, 0.0, 0.0,
74  0.1, 1.0, 1.0, 0.0, 1.0,
75  0.0, 0.0, 1.0, 1.0, 1.0,
76  1.0, 1.0, 1.0, 0.0, 1.0};
77 
78  static double fn2[20] = {0.0, 0.0, 1.0, 1.0, 1.0,
79  0.0, 0.0, 0.0, 1.0,-1.0,
80  1.0, 1.0,-1.0, 0.0, 0.0,
81  0.0, 0.0, 0.0, 1.0,-1.0};
82  static int first_call = 1;
83  static double E_cost, Unit_Time, N_H_rho, frac_Z, frac_He;
84 
85  if (first_call) {
86 
87  E_cost = UNIT_LENGTH/UNIT_DENSITY/pow(UNIT_VELOCITY, 3.0);
88  Unit_Time = UNIT_LENGTH/UNIT_VELOCITY;
89 
90  /* ------------------------------------------------------
91  conversion factor from total density to hydrogen
92  number density, i.e. nH = N_H_rho * rho
93  ------------------------------------------------------ */
94 
95  N_H_rho = (UNIT_DENSITY/CONST_amu)*(X/A_H); /* n(H) */
96  frac_He = (Y/A_He)*(A_H/X);
97  frac_Z = ((1 - X - Y)/A_Z)*(A_H/X);
98 
99  first_call = 0;
100  }
101 
102  rho = v[RHO];
103  pr = v[PRS];
104 
105 /* -----------------------------------
106  Force fneut to stay between [0,1]
107  ----------------------------------- */
108 
109  v[X_HI] = MAX(v[X_HI], 0.0);
110  v[X_HI] = MIN(v[X_HI], 1.0);
111 
112  if (v[PRS] < 0.0) v[PRS] = g_smallPressure;
113 
114  fn = v[X_HI];
115  mu = MeanMolecularWeight(v);
116  T = pr/rho*KELVIN*mu;
117 
118  /* --------------------------------------------------------
119  Set source terms equal to zero when the temperature
120  falls below g_minCoolingTemp
121  -------------------------------------------------------- */
122 
123 /* if (T < g_minCoolingTemp){
124  a = b = src_pr = 0.0;
125  continue;
126  }
127 */
128  if (mu < 0.0){
129  printf ("Negative mu in radiat \n");
130  exit(1);
131  }
132 
133  st = sqrt(T);
134  N_H = N_H_rho*rho; /* -- number density of hydrogen N_H = N(HI) + N(HII) -- */
135 
136  /* ---- coeff di ionizz. e ricomb. della particella fluida i,j ---- */
137 
138  cr = 2.6e-11/st;
139  ci = 1.08e-8*st*exp(-157890.0/T)/(13.6*13.6);
140 
141  n_el = N_H*(1.0 - fn + frac_Z); /* -- electron number density, in cm^{-3} -- */
142  rhs[X_HI] = Unit_Time*n_el*(-(ci + cr)*fn + cr);
143 
144  em[1] = 0.0;
145  for (k = 2; k <= 17; k++){
146  em[k] = 1.6e-12*8.63e-6*om[k]*ep[k]*exp(-t00[k]/T)/st;
147  em[k] = em[k]*critn[k]*st/(n_el + critn[k]*st);
148  em[k] = em[k]*ab[k]*(fn1[k] + fn*fn2[k]);
149  em[1] = em[1] + em[k];
150  }
151 
152  em[18] = ci*13.6*1.6e-12*fn;
153  em[19] = cr*0.67*1.6e-12*(1.0 - fn)*T/11590.0;
154  em[1] = em[1] + em[18] + em[19];
155 
156  /* ---------------------------------------------------
157  rlosst is the energy loss in units of erg/cm^3/s;
158  it must be multiplied by cost_E in order to match
159  non-dimensional units.
160  Source term for the neutral fraction scales with
161  UNIT_TIME
162  --------------------------------------------------- */
163  double g_gamma;
164 #if EOS == IDEAL
165  g_gamma = 1.4;
166 #elif EOS == PVTE_LAW
167  g_gamma = Gamma1(v);
168 #endif
169 
170  rlosst = em[1]*n_el*N_H;
171  rhs[PRS] = -E_cost*rlosst*(g_gamma - 1.0);
172  rhs[PRS] *= 1.0/(1.0 + exp( -(T - g_minCoolingTemp)/100.0)); /* -- cutoff -- */
173 
174 }
MAX
#define MAX(a, b)
Definition: macros.h:101
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
g_gamma
double g_gamma
Definition: globals.h:112
UNIT_DENSITY
#define UNIT_DENSITY
Unit density in gr/cm^3.
Definition: pluto.h:369
RHO
#define RHO
Definition: mod_defs.h:19
CONST_amu
#define CONST_amu
Atomic mass unit.
Definition: pluto.h:253
frac_He
#define frac_He
Definition: cooling.h:17
Y
#define Y
Definition: radiat.modified_by_Bhargav.c:4
frac_Z
#define frac_Z
Definition: cooling.h:15
g_smallPressure
double g_smallPressure
Small value for pressure fix.
Definition: globals.h:110
A_H
#define A_H
Definition: radiat.modified_by_Bhargav.c:8
KELVIN
#define KELVIN
Definition: pluto.h:401
MIN
#define MIN(a, b)
Definition: macros.h:104
UNIT_VELOCITY
#define UNIT_VELOCITY
Unit velocity in cm/sec.
Definition: pluto.h:377
X
#define X
Definition: radiat.modified_by_Bhargav.c:3
UNIT_LENGTH
#define UNIT_LENGTH
Unit Length in cm.
Definition: pluto.h:373
A_He
#define A_He
Definition: radiat.modified_by_Bhargav.c:7
X_HI
Definition: cooling.h:110
k
int k
Definition: analysis.c:2
Gamma1
double Gamma1(double *v)
Definition: pvte_law.c:231
g_minCoolingTemp
double g_minCoolingTemp
The minimum temperature (in K) below which cooling is suppressed.
Definition: globals.h:106
MeanMolecularWeight
double MeanMolecularWeight(double *V)
Definition: radiat.modified_by_Bhargav.c:178
A_Z
#define A_Z
Definition: radiat.modified_by_Bhargav.c:6

Here is the call graph for this function: