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

Go to the source code of this file.

Functions

void Radiat (double *v, double *rhs)
 
double H_MassFrac (void)
 
double CompEquil (double N, double T, double *v)
 

Function Documentation

double CompEquil ( double  N,
double  T,
double *  v 
)
Parameters
[in]nthe particle number density (not needed, but kept for compatibility)
[in]Tthe temperature (in K) for which equilibrium must be found.
[in,out]van array of primitive variables. On input, only density needs to be defined. On output, fractions will be updated to the equilibrium values.

Compute the equilibrium ionization balance for (rho,T)

Definition at line 205 of file radiat.c.

212 {
213  double cr, ci, st, n_e;
214  st = sqrt(T);
215  cr = 2.6e-11/st; /* recombination / ionization coefficients */
216  ci = 1.08e-8*st*exp(-157890.0/T)/(13.6*13.6);
217  v[X_HI] = cr / (ci + cr); /* compute fraction of neutrals at equilibrium */
218  n_e = (1.0 - v[X_HI] + frac_Z)*N;
219  return n_e;
220 }
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
frac_Z
#define frac_Z
Definition: cooling.h:15
X_HI
Definition: cooling.h:110

Here is the call graph for this function:

double H_MassFrac ( void  )

Compute the mass fraction X of Hydrogen as function of the composition of the gas. 

      f_H A_H

X = -------------— f_k A_k

where

f_k : is the fractional abundance (by number), f_k = N_k/N_tot of atomic species (no matter ionization degree).

A_K : is the atomic weight

Note: In this module, f_H = 1.0

where N_{tot} is the total number density of particles

ARGUMENTS

none

Definition at line 174 of file radiat.c.

200 {
201  return CONST_AH/(CONST_AH + frac_He*CONST_AHe + frac_Z*CONST_AZ);
202 }
CONST_AH
#define CONST_AH
Definition: pluto.h:250
frac_He
#define frac_He
Definition: cooling.h:17
frac_Z
#define frac_Z
Definition: cooling.h:15
CONST_AZ
#define CONST_AZ
Mean atomic weight of heavy elements.
Definition: pluto.h:252
CONST_AHe
#define CONST_AHe
Atomic weight of Helium.
Definition: pluto.h:251
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 4 of file radiat.c.

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

Here is the call graph for this function: