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