cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
highen.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*highen do high energy radiation field - gas interaction, Compton scattering, etc */
4 #include "cddefines.h"
5 #include "highen.h"
6 #include "trace.h"
7 #include "heavy.h"
8 #include "radius.h"
9 #include "magnetic.h"
10 #include "hextra.h"
11 #include "thermal.h"
12 #include "ionbal.h"
13 #include "opacity.h"
14 #include "pressure.h"
15 #include "gammas.h"
16 #include "rfield.h"
17 #include "doppvel.h"
18 #include "dense.h"
19 
20 void highen( void )
21 {
22  long int i,
23  ion,
24  nelem;
25 
26  double CosRayDen,
27  crsphi,
28  heatin,
29  sqthot;
30 
31  double hin;
32 
33  DEBUG_ENTRY( "highen()" );
34 
35 
36  /**********************************************************************
37  *
38  * COMPTON RECOIL IONIZATION
39  *
40  * bound electron scattering of >2.3 kev photons if neutral
41  * lgComptonOn usually t, f if "NO COMPTON EFFECT" command given
42  * lgCompRecoil usually t, f if "NO RECOIL IONIZATION" command given
43  *
44  **********************************************************************/
45 
46  /* lgComptonOn turned false with no compton command,
47  * lgCompRecoil - no recoil ionization */
49  {
50  for( nelem=0; nelem<LIMELM; ++nelem )
51  {
52  for( ion=0; ion<nelem+1; ++ion )
53  {
54  double CompRecoilIonRate = 0.,
55  CompRecoilHeatRate = 0.;
56  if( dense.xIonDense[nelem][ion] > 0. )
57  {
58  /* recoil ionization starts at 194 Ryd = 2.6 keV */
59  /* this is the ionization potential of the valence shell */
60  /* >>chng 02 may 27, lower limit is now 1 beyond actual threshold
61  * since recoil energy at threshold was very small, sometimes negative */
62  long limit = MIN2( rfield.nflux , rfield.nPositive );
63  double valence_IP_Ryd = Heavy.Valence_IP_Ryd[nelem][ion];
64  for( i=ionbal.ipCompRecoil[nelem][ion]; i < limit; ++i)
65  {
66  crsphi = opac.OpacStack[i-1+opac.iopcom] * rfield.SummedCon[i];
67 
68  /* direct hydrogen ionization due to compton scattering,
69  * does not yet include secondaries,
70  * last term accounts for number of valence electrons that contribute */
71  CompRecoilIonRate += crsphi;
72 
73  /* recoil energy in Rydbergs
74  * heating modified for suprathermal secondaries below; ANU2=ANU**2 */
75  /* >>chng 02 mar 27, use general formula for recoil energy */
76  /*energy = 2.66e-5*rfield.anu2(i) - 1.;*/
77  double recoil_energy = rfield.anu2(i) *
78  ( EN1RYD /
79  ( ELECTRON_MASS * SPEEDLIGHT * SPEEDLIGHT) ) - valence_IP_Ryd;
80 
81  /* heating is in rydbergs because SecIon2PrimaryErg, SecExcitLya2PrimaryErg, HeatEfficPrimary in ryd */
82  CompRecoilHeatRate += crsphi*recoil_energy;
83  }
84  /* net heating rate, per atom, convert ryd/sec/cm3 to ergs/sec/atom */
85  CompRecoilHeatRate *= EN1RYD;
86 
87  /* this is number of electrons in valence shell of this ion */
88  long nElec = ionbal.nCompRecoilElec[nelem-ion];
89 
90  ionbal.CompRecoilHeatRate[nelem][ion] = nElec*CompRecoilHeatRate;
91  ionbal.CompRecoilIonRate[nelem][ion] = nElec*CompRecoilIonRate;
92 
93  ASSERT( ionbal.CompRecoilHeatRate[nelem][ion] >= 0.);
94  ASSERT( ionbal.CompRecoilIonRate[nelem][ion] >= 0.);
95  }
96  }
97  }
98  }
99  else
100  {
101  for( nelem=0; nelem<LIMELM; ++nelem )
102  {
103  for( ion=0; ion<nelem+1; ++ion )
104  {
105  ionbal.CompRecoilIonRate[nelem][ion] = 0.;
106  ionbal.CompRecoilHeatRate[nelem][ion] = 0.;
107  }
108  }
109  }
110 
111  /**********************************************************************
112  *
113  * COSMIC RAYS
114  *
115  * heating and ionization by cosmic rays, other relativistic particles
116  * CRYDEN=density (1/CM**3), neutral rate assumes 15ev total
117  * energy loss, 13.6 into ionization, 1.4 into heating
118  * units erg/sec/cm**3
119  * iff not specified, CRTEMP is 2.6E9K
120  *
121  **********************************************************************/
122 
123  if( hextra.cryden > 0. )
124  {
125  ASSERT( hextra.crtemp > 0. );
126  /* this is current cosmic ray density, as determined from original density times
127  * possible dependence on radius */
129  {
130  /* >>chng 06 jun 02, add this option
131  * this is case where cr are in equipartition with magnetic field -
132  * set with COSMIC RAY EQUIPARTITION command */
133  CosRayDen = hextra.background_density *
134  /* ratio of energy density in current B to typical galactic
135  * galactic background energy density of 1.8 eV cm-3 is from
136  *>>refer cr background Webber, W.R. 1998, ApJ, 506, 329 */
138  (CR_EDEN_GAL_BACK_EV_CMM3/*1.8eV cm-3*/ * EN1EV/*erg eV-1*/ );
139  }
140  else
141  {
142  /* this is usual case, CR density may depend on radius, usually does not */
143  CosRayDen = hextra.cryden*pow(radius.Radius/radius.rinner,(double)hextra.crpowr);
144  }
145 
146  /* cosmic ray energy density rescaled by ratio to background ion rate and B field */
148  (CR_EDEN_GAL_BACK_EV_CMM3/*1.8eV cm-3*/ * EN1EV/*erg eV-1*/ );
149 
150  /* related to current temperature, when thermal */
151  sqthot = sqrt(hextra.crtemp);
152 
153  /* rate hot electrons heat cold ones, Balbus and McKee 1982
154  * units erg sec^-1 cm^-3,
155  * in sumheat we will multipy this rate by sum of neturals, but for this
156  * term we actually want eden, so mult by eden over sum of neut */
157  ionbal.CosRayHeatThermalElectrons = 5.5e-14/sqthot*CosRayDen;
158 
159  /* ionization rate; Balbus and McKee */
160  ionbal.CosRayIonRate = 1.22e-4/sqthot*
161  log(2.72*pow(hextra.crtemp/1e8,0.097))*CosRayDen;
162 
163  /* option for thermal CRs, first is the usual (and default) relativistic case */
164  if( hextra.crtemp > 2e9 )
165  {
166  /* usual circumstance; relativistic cosmic rays,
167  * cosmic ray ionization rate s-1 cm-3; ext rel limit */
168  ionbal.CosRayIonRate *= 3.;
169 
170  }
171  else
172  {
173  /* option for thermal cosmic rays */
174  ionbal.CosRayIonRate *= 10.;
175  }
176  /* >>chng 04 jan 27, from 0.093 to 2.574 as per following */
177  /* cr heating from Table 1 of
178  *>>refer cr heating Wolfire et al.1995, ApJ, 443, 152
179  * For every ionization due to cosmic rays, ~35eV of heat is added
180  * to the system. This manifests itself in the ionbal.CosRayHeatNeutralParticles term
181  * by the 2.574*EN1RYD term, which is just the energy in ergs in 35 eV.
182  * Change made by Nick Abel and Gargi Shaw, 04 Jan 27. In heatsum
183  * we multiply by the number of secondaries that occur */
185 
186  if( trace.lgTrace )
187  {
188  fprintf( ioQQQ, " highen: cosmic ray density;%10.2e CRion rate;%10.2e CR heat rate=;%10.2e CRtemp;%10.2e\n",
190  }
191  }
192  else
193  {
194  ionbal.CosRayIonRate = 0.;
196  }
197  /* >>chng 06 may 23, Penning ionization
198  ionbal.CosRayIonRate += 1e-9 *
199  iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe2s3S].Pop; */
200 
201  /*fprintf(ioQQQ,"DEBUG cr %.2f %.3e %.3e %.3e\n",
202  fnzone,
203  hextra.cryden ,
204  ionbal.CosRayIonRate ,
205  ionbal.CosRayHeatNeutralParticles );*/
206 
207  /**********************************************************************
208  *
209  * add on extra heating due to turbulence, goes into [1] of [x][0][11][0]
210  *
211  **********************************************************************/
212 
213  /* TurbHeat added with hextra command, DispScale added with turbulence dissipative */
214  if( (hextra.TurbHeat+DoppVel.DispScale) != 0. )
215  {
216  /* turbulent heating only goes into the low-energy heat of this element */
217  /* >>>>chng 00 apr 28, functional form of radius dependence had bee turrad/depth
218  * and so went to infinity at the illuminated face. Changed to current form as
219  * per Ivan Hubeny comment */
220  if( hextra.lgHextraDepth )
221  {
222  /* if turrad is >0 then vary heating with depth */
225 
226  /* >>chng 00 aug 16, added option for heating from back side */
227  if( hextra.turback != 0. )
228  {
231  }
232  }
233  else if( hextra.lgHextraDensity )
234  {
235  /* depends on density */
238  }
239  else if( hextra.lgHextraSS )
240  {
241  /* with SS disk model */
244  powi((double)hextra.HextraSSradius,-3));
245  }
246  /* this is turbulence dissipate command */
247  else if( DoppVel.DispScale > 0. )
248  {
249  double turb = DoppVel.TurbVel * sexp( radius.depth / DoppVel.DispScale );
250  /* if turrad is >0 then vary heating with depth */
251  /* >>chng 01 may 10, add extra factor of length over 1e13 cm */
252  ionbal.ExtraHeatRate = 3.45e-28 / 2.82824 * turb * turb * turb *
253  ( dense.gas_phase[ipHYDROGEN] / 1e10 ) / (DoppVel.DispScale/1e13);
254  }
255  else
256  {
257  /* constant extra heating */
259  }
260  }
261 
262  else
263  {
264  ionbal.ExtraHeatRate = 0.;
265  }
266 
267  /**********************************************************************
268  *
269  * option to add on fast neutron heating, goes into [0] & [2] of [x][0][11][0]
270  *
271  **********************************************************************/
272  if( hextra.lgNeutrnHeatOn )
273  {
274  /* hextra.totneu is energy flux erg cm-2 s-1
275  * CrsSecNeutron is 4E-26 cm^-2, cross sec for stopping rel neutrons
276  * this is heating erg s-1 due to fast neutrons, assumed to secondary ionize */
277  /* neutrons assumed to only secondary ionize */
279  }
280  else
281  {
283  }
284 
285 
286  /**********************************************************************
287  *
288  * pair production in elec field of nucleus
289  *
290  **********************************************************************/
291  t_phoHeat photoHeat;
293  ionbal.PairProducPhotoRate[1] = photoHeat.HeatLowEnr;
294  ionbal.PairProducPhotoRate[2] = photoHeat.HeatHiEnr;
295 
296  /**********************************************************************
297  *
298  * Compton energy exchange
299  *
300  **********************************************************************/
301  rfield.cmcool = 0.;
302  rfield.cmheat = 0.;
303  heatin = 0.;
304  /* lgComptonOn usually t, turns off Compton */
305  if( rfield.lgComptonOn )
306  {
307  for( i=0; i < rfield.nflux; i++ )
308  {
309 
310  /* Compton cooling
311  * CSIGC is Tarter expression times ANU(I)*3.858E-25
312  * 6.338E-6 is k in inf mass Rydbergs, still needs factor of TE */
313  rfield.comup[i] = (double)(rfield.flux[0][i]+rfield.ConInterOut[i]+
314  rfield.outlin[0][i]+ rfield.outlin_noplot[i])*rfield.csigc[i]*(dense.eden*4.e0/
315  TE1RYD*1e-15);
316 
317  rfield.cmcool += rfield.comup[i];
318 
319  /* Compton heating
320  * CSIGH is Tarter expression times ANU(I)**2 * 3.858E-25
321  * CMHEAT is just spontaneous, HEATIN is just induced */
322  rfield.comdn[i] = (double)(rfield.flux[0][i]+rfield.ConInterOut[i]+
323  rfield.outlin[0][i]+ rfield.outlin_noplot[i])*rfield.csigh[i]*dense.eden*1e-15;
324 
325  /* induced Compton heating */
326  hin = (double)(rfield.flux[0][i]+rfield.ConInterOut[i]+rfield.outlin[0][i]+rfield.outlin_noplot[i])*
328  rfield.comdn[i] += hin;
329  heatin += hin;
330 
331  /* following is total compton heating */
332  rfield.cmheat += rfield.comdn[i];
333  }
334 
335  /* remember how important induced compton heating is */
336  if( rfield.cmheat > 0. )
337  {
339  }
340 
341  if( trace.lgTrace && trace.lgComBug )
342  {
343  fprintf( ioQQQ, " HIGHEN: COOL num=%8.2e HEAT num=%8.2e\n",
345  }
346  }
347 
348  if( trace.lgTrace && trace.lgComBug )
349  {
350  fprintf( ioQQQ,
351  " HIGHEN finds heating fracs= frac(compt)=%10.2e "
352  " f(pair)%10.2e totHeat=%10.2e\n",
354  thermal.heating(0,21)/thermal.htot,
355  thermal.htot );
356  }
357  return;
358 }
double Radius
Definition: radius.h:31
realnum * csigh
Definition: rfield.h:271
double HeatHiEnr
Definition: thermal.h:208
double depth
Definition: radius.h:31
double htot
Definition: thermal.h:169
t_thermal thermal
Definition: thermal.cpp:6
long int nCompRecoilElec[LIMELM]
Definition: ionbal.h:181
double * OpacStack
Definition: opacity.h:163
double GammaK(long int ipLoEnr, long int ipHiEnr, long int ipOpac, double yield1, t_phoHeat *photoHeat)
double ** CompRecoilIonRate
Definition: ionbal.h:159
t_opac opac
Definition: opacity.cpp:5
double ** CompRecoilHeatRate
Definition: ionbal.h:165
realnum ** flux
Definition: rfield.h:70
t_Heavy Heavy
Definition: heavy.cpp:5
realnum turrad
Definition: hextra.h:51
double CosRayHeatNeutralParticles
Definition: ionbal.h:128
realnum turback
Definition: hextra.h:51
realnum * outlin_noplot
Definition: rfield.h:191
t_magnetic magnetic
Definition: magnetic.cpp:17
bool lgComBug
Definition: trace.h:37
realnum crtemp
Definition: hextra.h:24
realnum HextraSSradius
Definition: hextra.h:73
t_hextra hextra
Definition: hextra.cpp:5
bool lgNeutrnHeatOn
Definition: hextra.h:81
double * SummedCon
Definition: rfield.h:163
bool lg_CR_B_equipartition
Definition: hextra.h:29
realnum TurbVel
Definition: doppvel.h:21
sys_float sexp(sys_float x)
Definition: service.cpp:1095
double CosRayIonRate
Definition: ionbal.h:124
realnum ** outlin
Definition: rfield.h:191
FILE * ioQQQ
Definition: cddefines.cpp:7
t_DoppVel DoppVel
Definition: doppvel.cpp:5
#define MIN2(a, b)
Definition: cddefines.h:807
t_dense dense
Definition: global.cpp:15
double CrsSecNeutron
Definition: hextra.h:87
double * comup
Definition: rfield.h:238
double ExtraHeatRate
Definition: ionbal.h:135
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
double energydensity
Definition: magnetic.h:44
t_trace trace
Definition: trace.cpp:5
t_ionbal ionbal
Definition: ionbal.cpp:8
double rinner
Definition: radius.h:31
double CosRayHeatThermalElectrons
Definition: ionbal.h:132
double xNeutronHeatRate
Definition: ionbal.h:139
bool lgTrace
Definition: trace.h:12
realnum crpowr
Definition: hextra.h:24
realnum TurbHeat
Definition: hextra.h:42
double cmcool
Definition: rfield.h:277
double cr_energydensity
Definition: hextra.h:32
t_pressure pressure
Definition: pressure.cpp:9
t_rfield rfield
Definition: rfield.cpp:9
realnum DispScale
Definition: doppvel.h:51
realnum * ConInterOut
Definition: rfield.h:156
#define CR_EDEN_GAL_BACK_EV_CMM3
Definition: hextra.h:13
long int iopcom
Definition: opacity.h:222
realnum HextraScaleDensity
Definition: hextra.h:60
double cmheat
Definition: rfield.h:277
double HeatLowEnr
Definition: thermal.h:206
double powi(double, long int)
Definition: service.cpp:690
double anu2(size_t i) const
Definition: mesh.h:115
realnum * OccNumbIncidCont
Definition: rfield.h:119
double PresGasCurr
Definition: pressure.h:46
realnum cryden
Definition: hextra.h:24
t_radius radius
Definition: radius.cpp:5
double heating(long nelem, long ion)
Definition: thermal.h:186
long int ioppr
Definition: opacity.h:222
realnum totneu
Definition: hextra.h:79
bool lgHextraSS
Definition: hextra.h:64
realnum gas_phase[LIMELM]
Definition: dense.h:76
double HextraSS_M
Definition: hextra.h:70
#define ASSERT(exp)
Definition: cddefines.h:617
long int ippr
Definition: opacity.h:222
const int LIMELM
Definition: cddefines.h:307
bool lgHextraDensity
Definition: hextra.h:57
double Valence_IP_Ryd[LIMELM][LIMELM]
Definition: heavy.h:24
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
double * comdn
Definition: rfield.h:238
realnum * csigc
Definition: rfield.h:271
double eden
Definition: dense.h:201
#define MAX2(a, b)
Definition: cddefines.h:828
bool lgCompRecoil
Definition: ionbal.h:150
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
double pow(double x, int i)
Definition: cddefines.h:782
realnum HextraSSalpha
Definition: hextra.h:67
double PairProducPhotoRate[3]
Definition: ionbal.h:142
void highen(void)
Definition: highen.cpp:20
const int ipHYDROGEN
Definition: cddefines.h:348
long int nflux
Definition: rfield.h:48
realnum background_density
Definition: hextra.h:38
long int nPositive
Definition: rfield.h:55
bool lgComptonOn
Definition: rfield.h:281
long int ** ipCompRecoil
Definition: ionbal.h:156
double cinrat
Definition: rfield.h:277
bool lgHextraDepth
Definition: hextra.h:49